# Analytical investigation of time-dependent two-dimensional non-Newtonian
boundary layer equations

Imre F. Barna1, G. Bognár2, L. Mátyás3 and K. Hriczó2 1 Hungarian Research
Network, Wigner Research Centre for Physics,
Konkoly-Thege Miklós út 29 - 33, 1121 Budapest, Hungary
2 University of Miskolc, Miskolc-Egyetemváros 3515, Hungary,
3Department of Bioengineering, Faculty of Economics, Socio-Human Sciences and
Engineering, Sapientia Hungarian University of Transylvania, Libertătii sq. 1,
530104 Miercurea Ciuc, Romania

###### Abstract

In this study, five different time-dependent incompressible non-Newtonian
boundary layer models in two dimensions are investigated, including external
magnetic field effects. The power law, the Casson fluid, the Oldroyd-B model,
the Walter fluid B model and the Williamson fluid are analyzed. For the first
two models, analytical results are given for the velocity and pressure
distributions, which can be expressed by different types of hypergeometric
functions. Depending on the parameters involved in the analytical solutions of
the nonlinear ordinary differential equation obtained by the similarity
transformation, a very wide range of solution types is presented.

non-Newtonian fluid, self-similar method, boundary layer, MHD flow, time-
dependent solution

###### pacs:

47.10.-g, 47.10.ab, 47.10.ad, 67.57.De

##  I Introduction

The scientific field of classical fluid mechanics is vast and cannot be
summarised in a few finite arbitrary books. If we restrict ourselves to non-
Newtonian fluids or boundary layer flows, the relevant literature is still
remarkably large. The fundamental physics of such fluid flows can be found in
basic textbooks such as arista , non-newt , patel  schli . In our previous
publication on time-dependent self-similar solutions of compressible and
incompressible heated boundary layer equations barna-bound , we collected the
most relevant literature in the field, which we omit here. Some analytic
results are also available for stationary non-Newtonian boundary layers, for
example, saer , bog , bog2 , mak .

Today, there is a growing need to study the motion of non-Newtonian fluids,
which are often found in industrial applications and in modelling many
manufacturing processes (heat exchangers, pharmaceuticals, food, paper). All
fluids for which the shear stress-shear rate relation cannot be described by
Newton’s law are called non-Newtonian fluids. They can be characterized by
formulas of an extremely diverse nature. Without completeness, we will
investigate five different non-Newtonian fluid models in our study these are
the following: The first is the non-Newtonian power-law model describes the
viscosity of time-independent flow behaviour for both shear thinning (or
pseudoplastic) and shear thickening (or dilatant) cases (see arista , schli ).
This model is used in vast engineering applications, for instance, lubricants,
polymers, and slurries yori , this model is well suited for numerical and
analytical studies of fluid flows bog , bog2 .

The second can be directly derived from the usual Newtonian fluids. This is
the widely used viscosity model for viscoplastic non-Newtonian behaviour and
named the Casson fluid model. Fluids that behave like a solid when the shear
stress is less than the yield stress applied to a liquid, and start to deform
when the shear stress is greater than the yield stress casson .

Our third model is the non-Newtonian Oldroyd-B fluid is used in the literature
for the rheological characterization of a type of viscoelastic fluid. Oldroyd
applied the principle that stresses in a continuous medium can only arise from
deformations and cannot change if the material is only rotated old .

The fourth is the Walters-B model which was proposed by Walters for
viscoelastic fluids walt . It is a generalisation of the Oldroyd-B model and
is used to describe the behaviour of fluids that are mainly used in food and
food processing technologies.

The final model is the Williamson fluid model which can also be used to
describe the viscoelastic shear thinning properties of non-Newtonian fluids
wil . In Williamson’s fluid model, the effective viscosity should decrease
infinitely as the shear rate increases, which is effectively infinite
viscosity at rest (no fluid motion) and zero viscosity as the shear rate
approaches infinity.

The actual existing literature on non-Newtonian boundary layers is enormous,
therefore we cannot give a general overview just mention some relevant recent
publications in the right place in this study.

To analyse the variation in flow characteristics of a fluid in a magnetic
field, an additional magnetic term is added to the momentum equation of the
boundary layer equation system, so that the boundary layer equation system can
describe the magnetic effects.

The aim of this paper is to provide analytic results for the stationary
equations of the non-Newtonian fluid flows in a magnetic field. The above
given five types of non-Newtonian behaviour are examined. The governing
equations are transformed into a system of ordinary differential equations by
a similarity transformation. We analyze the cases in which the fluid flow rate
functions can be analytically given by the method used. In these two cases, we
give the velocity distributions in two space dimensions and point out the
effect of each parameter.

##  II Mathematical models

###  II.1 Power-law viscosity

For power-law viscosity, our partial differential equation system systems
reads as follows:

$\displaystyle\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}$
$\displaystyle=$ $\displaystyle 0,$ (1) $\displaystyle\rho\left(\frac{\partial
u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial
y}\right)$ $\displaystyle=$ $\displaystyle-\frac{\partial p}{\partial
x}+m\frac{\partial}{\partial y}\left(\left|\frac{\partial u}{\partial
y}\right|^{n-1}\cdot\frac{\partial u}{\partial y}\right)+\sigma B_{0}^{2}u,$
(2) $\displaystyle\frac{\partial p}{\partial y}$ $\displaystyle=$
$\displaystyle 0,$ (3)

where the dynamical variables are the velocities components
$u(x,y,t),v(x,y,t)$, the fluid pressure $p(x,y,t)$, $B_{0}=const$ is the
magnetic the magnetic induction, and $\sigma$ is the electrical conductivity.
Additional physical parameters are $\rho,m,n$ the fluid density, consistency
parameter, and the power-law index, respectively. To fix the geometrical
relations in our system consider Fig.(1).

The temporal decay of the strong solutions of the power-law type of non-
Newtonian fluids for variable power-law index was mathematically proven by Ko
power1 . Flow reversal effects in an expanding channel for this type of
viscous fluid were analyzed by power2 . A non-iterative transformation method
to an extended Blasius problem describing a 2D laminar boundary-layer with
power-law viscosity for non-Newtonian fluids was developed by Fazio power3 .
The stationary power-law equations were solved by Pater et al. patel  by the
one-parameter deductive group theory technique. (We have to mention at this
point that the stationary boundary layer equations can be reduced to the
third-order non-linear differential equation which is called the Blasius
equation having an extensive literature which we skip now.)

To obtain the solutions to the system (1)-(3), we apply the following self-
similar Ansatz for $u,v$ and $p$:

$\displaystyle u(x,y,t)=t^{-\alpha}f(\eta),\hskip
28.45274ptv(x,y,t)=t^{-\delta}g(\eta),\hskip
28.45274ptp(x,y,t)=t^{-\gamma}h(\eta),$ (4)

with the argument $\eta=\frac{x+y}{t^{\beta}}$ of the shape functions. All the
exponents $\alpha,\beta,\gamma,\delta$ are real numbers. Solutions with
integer exponents are called self-similar solutions of the first kind, non-
integer exponents generate self-similar solutions of the second kind sedov .

It is necessary to remark that, the heat conduction mechanism to (1) - (3) the
boundary layer equations can exclude the self-similarity using transformation
(4) one gets contradiction among self-similar exponents.

The shape functions $f,g$ and $h$ could be any arbitrary continuous functions
with existing first and second continuous derivatives and will be derived
later on. The physical and geometrical interpretation of the Ansatz were
exhaustively analyzed in former publications imre1 ; imre3 , therefore we skip
it here.

The main points are, that $\alpha,\delta,\gamma$ are responsible for the rate
of decay and $\beta$ is for the rate of spreading of the corresponding
dynamical variable for positive exponents. Negative exponents mean physically
irrelevant cases, blowing up and contracting solutions. The numerical values
of the exponents are considered as follows:

$\alpha=\delta=n/2,\hskip 28.45274pt\beta=1-\alpha=1-n/2,\hskip
28.45274pt\gamma=n.$ (5)

Exponents with numerical values of one-half mean the regular Fourier heat
conduction (or Fick’s diffusion) process. One-half values for the exponent of
the velocity components and unit value exponent for the pressure decay are
usual for the incompressible Navier-Stokes equation imre4 . Note, that the
value of $\gamma$ is responsible for the decay of the pressure field. °

Figure 1: Defining the directions and the velocity components of the
investigated system.

Applying (5), the derived ordinary differential equation (ODE) system reads as

$f^{\prime}+g^{\prime}=0,$ (6)
$\rho\left(-\frac{n}{2}f-\left[1-\frac{n}{2}\right]\eta
f^{\prime}+ff^{\prime}+gf^{\prime}\right)=-h^{\prime}+mnf^{\prime\prime}|f^{\prime}|^{n-1}+\sigma
B_{0}^{2}f,$ (7) $h^{\prime}=0,$ (8)

where the prime means derivation concerning variable $\eta$. Equations (6) and
(8) can be integrated, to get $f+g=c_{1}$ and $h=c_{2}$. However, the dynamic
variable under consideration is the velocity component $u$, which is $f$. From
(7) the derived ODE is the following:

$\frac{mn}{\rho}f^{\prime\prime}|f^{\prime}|^{n-1}+\left(\left[1-\frac{n}{2}\right]\eta-
c_{1}\right)f^{\prime}+Kf=0,$ (9)

where $K=(\rho n/2+\sigma B_{0}^{2}).$ The constant $c_{1}$ can be determined
from possible boundary conditions applied to the fluid flow problem.

Here, we give analytic and numeric solutions to the flow (6)-(8).

i.) First we consider the case $B_{0}=0$. Equation (9) has analytic solutions
only for the cases of $n=1,0$ and $-1$.

For the Newtonian fluid case, when $n=1$, and other parameters
($c_{1},m,\rho$) are general, the solution reads as

$f(\eta)=C_{1}M\left(\frac{\rho}{2},\frac{1}{2},-\frac{\rho(\eta-2c_{1})^{2}}{4m}\right)+C_{2}U\left(\frac{\rho}{2},\frac{1}{2},-\frac{\rho(\eta-2c_{1})^{2}}{4m}\right),$
(10)

where $M(,,)$ and $U(,,)$ are the Kummer M and U functions NIST  and $C_{1}$
and $C_{2}$ stand for the usual integration constants. The parameter $c_{1}$
just shifts the maxima of the shape function and $m$ the consistency parameter
(which is positive) just changes the full width at half the maximum of the
function. The larger the numerical value of $m$ the broader the shape
function. The third parameter, the density of the fluid (which must be also
positive) is the most relevant parameter of the flow which of course, meets
our physical intuition. Figure (2) shows four different velocity shape
functions for different fluid densities. It is important to mention that for
density $0<\rho$ the shape function has a global maxima in the origin and a
decay to zero, (larger densities mean quicker decay.) however for $1<\rho$ the
shape function shows additional oscillations as well. The final velocity field
has the asymptotic form of
$u(x,y=0,t)=t^{-\frac{1}{2}}f\left(\frac{x}{t^{1/2}}\right)$. For
completeness, we give the entire formula for $u$ as well:

$u(x,y,t)=\frac{1}{t^{1/2}}\left(C_{1}M\left[\frac{\rho}{2},\frac{1}{2},-\frac{\rho\left\\{\frac{(x+y)}{t^{1/2}}-2c_{1}\right\\}^{2}}{4m}\right]+C_{2}U\left[\frac{\rho}{2},\frac{1}{2},-\frac{\rho\left\\{\frac{(x+y)}{t^{1/2}}-2c_{1}\right\\}^{2}}{4m}\right]\right).$
(11)

We remark that this case was investigated in barna-bound . Figure (3) shows
the velocity field of $v(x,y=0,t)$ for a parameter set. The slight oscillation
and the quick decay are clear to see. It is interesting to mention here, that
during our investigations using the self-similar Ansatz, we regularly get
solutions which contain Kummer’s function for hydrodynamic processes like the
Bénard-Rayleigh convection imre1  or even for regular diffusion equation
laci_roman .

Figure 2: Four different velocity shape functions for Eq. (10). The black,
blue, red and green curves are for $\rho=0.5,1,2$ and $5$, respectively. All
additional parameters are the same for all four curves
($C_{1}=1,C_{2}=0,c_{1}=0$ and $m=2$ ).

Figure 3: The velocity distribution function projection for $y=0$ for the
density of $\rho=2$. All other parameters are the same as below.

For the case $n=K=0$ the solution is a trivial constant $f=c_{2}$.

There are analytic solutions available for $n=-1$ unfortunately, not for the
general case where all the other three parameters $(c_{1},\rho,m)$ are
completely free. The constant $c_{1}$ should be equal to zero and the ratio of
$\rho/m$ should have some special values. The existing possibilities are not
so many.

In the following we present some solutions which came from our numerical
experimental experiences:
For $\rho=1$ and for arbitrary integer or rational m we get the next implicit
relation for the solutions:

$ln(\eta)-ln(f[\eta])-ln\left(C_{1}U\left[-1,\frac{3}{2},\frac{1}{4m}f\\{\eta\\}^{2}\right]+M\left[-1,\frac{3}{2},\frac{1}{4m}f\\{\eta\\}^{2}\right]\right)-C_{2}=0.$
(12)

If the first parameters of the Kummer’s functions are non-positive integers
the series becomes finite

$\displaystyle
M\left[-1,\frac{3}{2},\frac{1}{4m}f\\{\eta\\}^{2}\right]=1-\frac{2}{3m}f(\eta)^{2},$
$\displaystyle
U\left[-1,\frac{3}{2},\frac{1}{4m}f\\{\eta\\}^{2}\right]=-\frac{3}{2}+\frac{1}{m}f(\eta)^{2},$
(13)

considering that $C_{2}>0$ and $\tilde{C_{2}}=ln({C_{2}})$ after some trivial
algebraic steps we get:

$\left(\frac{\tilde{C}_{2}C_{1}}{m}-\frac{2\tilde{C}_{2}}{m}\right)f(\eta)^{3}+\left(\tilde{C}_{2}-\frac{3\tilde{C}_{2}C_{1}}{2}\right)f(\eta)-\eta=0.$
(14)

The velocity shape functions can be seen in figure 4.

Figure 4: The implicit velocity shape functions of Eq. (14) for three
different parameter sets $(C_{1},\tilde{C}_{2},m)$. The black, red and blue
curves are for $(1,4,0.4)$, $(2.4,0.3,1.3)$ and $(2.4,0.3,10)$ .

One can see on Eq. (14), that if $\eta=0$, than $f(\eta=0)=0$ is a solution of
the equation. Further possibilities gives the following rearrangement

$\left[\left(\frac{\tilde{C}_{2}C_{1}}{m}-\frac{2\tilde{C}_{2}}{m}\right)f(\eta)^{2}+\left(\tilde{C}_{2}-\frac{3\tilde{C}_{2}C_{1}}{2}\right)\right]f(\eta)=0.$
(15)

This yields

$f(\eta=0)=\pm\sqrt{\left(\frac{3\tilde{C}_{2}C_{1}}{2}-\tilde{C}_{2}\right)/\left(\frac{\tilde{C}_{2}C_{1}}{m}-\frac{2\tilde{C}_{2}}{m}\right)}.$
(16)

It shows that for $\eta=0$ there is one or there are three corresponding
values to function $f(\eta)$. In the following we try to find the constraints
on the parameters which separate the two cases. If we consider the inverse
function $f(\eta)=y$, we have

$\eta=\left(\frac{\tilde{C}_{2}C_{1}}{m}-\frac{2\tilde{C}_{2}}{m}\right)\cdot
y^{3}+\left(\tilde{C}_{2}-\frac{3\tilde{C}_{2}C_{1}}{2}\right)\cdot y.$ (17)

The derivative of this function is

$\frac{d\eta}{dy}=3\left(\frac{\tilde{C}_{2}C_{1}}{m}-\frac{2\tilde{C}_{2}}{m}\right)\cdot
y^{2}+\left(\tilde{C}_{2}-\frac{3\tilde{C}_{2}C_{1}}{2}\right).$ (18)

If this expression has always the same sign, the function is monotonous, and
then there is just one single root $\eta(y=0)=0$ (i.e. $f(\eta=0)=0$). If this
derivative may change the sign, then multiple roots are possible. The inverse
function is not monotonous and may have three roots if there are real roots of
the equation of derivative $d\eta/dy=0$. This means the following condition

$\frac{\left(\frac{3\tilde{C}_{2}C_{1}}{2}-\tilde{C}_{2}\right)}{\left(\frac{\tilde{C}_{2}C_{1}}{m}-\frac{2\tilde{C}_{2}}{m}\right)}>0.$
(19)

Figure (4) shows two implicit solutions of Eq. (14) for three different
parameter sets. Note, that Eq.(14) is a third-order equation we may get multi-
valued solutions. If the parameter $m$ is much smaller than the two
integration constants then single-valued solutions emerge as well. Whether the
multi-valued solutions have any physical significance is not yet clear, but it
is thought that this property may indicate the existence of finite
oscillations or eddies. Figure (5) shows the projected velocity fields $(y=0)$
for single and double-valued case. (For a better visibility we present the
solution function from an unusual perspective.) The decrease of the velocity
field with time is clearly visible. For the sake of completeness, we give the
final form of the velocity field, which reads:

$\left(\frac{\tilde{C}_{2}C_{1}}{m}-\frac{2\tilde{C}_{2}}{m}\right)t^{3/2}f\left(\frac{(x+y)}{t^{3/2}}\right)^{3}+\left(\tilde{C}_{2}-\frac{3\tilde{C}_{2}C_{1}}{2}\right)t^{1/2}f\left(\frac{(x+y)}{t^{3/2}}\right)-\frac{(x+y)}{t^{3/2}}=0.$
(20)

a) b)

Figure 5: The implicit velocity distribution $u(x,t)$ evaluated from Eq. (21)
for two different parameter sets a) figure is for $(1,4,0.4)$ and b) figure is
for $(0,2.5,0.7)$, respectively.

There are analytic solutions available for $\rho=3$ and for arbitrary real $m$
and $c_{1}$ in the implicit form of:

$\frac{\eta}{f(\eta)}+\frac{C_{1}e^{\frac{9f(\eta)^{2}}{4m}}}{3f(\eta)}-\frac{C_{1}\cdot
erf\left(\frac{3}{2}\sqrt{-\frac{1}{m}}f[\eta]\right)}{2m\sqrt{-\frac{1}{m\pi}}}-\frac{2c_{1}}{3f(\eta)}-C_{2}=0.$
(21)

Equation (21) clearly shows, that for $C_{1}=0$, we get the trivial linear
function as a solution. Figure (6) presents the solutions of Eq. (21) for two
parameter sets. Because of the quadratic appearance of $f(\eta)$, we obtain
dual-valued curves which are not strictly speaking functions. The result
functions increase strongly (or decrease in the case of a negative branch)
very close to the origin and for larger arguments on a nearly horizontal
plateau. The general structure of the solution is quite stable, even a tenfold
change in the parameters does not change the overall feature of the derived
curve. To have a better overview of the properties of the solution, Figure (7)
presents the usual velocity projection of $u(x,y=0,t)$ for a given parameter
set. The physically relevant rapid decay over time is again clearly visible.

Figure 6: The implicit velocity shape functions of Eq.(21) for three different
parameter sets $(C_{1},\tilde{C}_{2},c_{1},m)$ with $\rho=3$. The black, red
and blue curves are for $(1,1,1,0,1)$, $(3,5,0,0,0.2)$ and $(13,11,4,10)$,
respectively

.

Figure 7: The implicit velocity distribution field when the parameters are
$\rho=3$, $(C_{1},\tilde{C}_{2},c_{1},m)$ and $(3,5,0,0.2)$ with $\rho=3$.

ii.) In the second part of our analysis, we investigate the solutions when the
magnetic induction is not equal to zero ($B_{0}\neq 0$). Even now four
different cases exist $n=2,1,\ 0$ and $-1$. We have to examine them one by
one.

First, we consider the dilatant case $n=2$. Equation (9) is reduced to

$\frac{2m}{\rho}f^{\prime\prime}f^{\prime}+Kf=0.$ (22)

multiplying by $f^{\prime}$ we get a total derivative which can be integrated
giving us the ODE of:

$\frac{2m}{3\rho}f^{\prime 3}+\frac{1}{2}(\rho+\sigma B_{0}^{2})f"+c_{0}=0.$
(23)

First considering the $c_{0}=0$ simpler case, the ODE can be directly
integrated and gives multiple solutions: there is a trivial one of $f=0$ there
are two complex conjugated solutions which we skip and a relevant real one.
After some algebraic manipulation, it reads as:

$f=-\frac{6\rho(\rho+\sigma B_{0}^{2})(\eta-c_{1})^{3}}{216\cdot m},$ (24)

which is a third-order parabola in the variable of $\eta$. The final velocity
distribution is

$u(x,y,t)=-\frac{6\rho(\rho+\sigma B_{0}^{2})(x+y-c_{1})^{3}}{216\cdot m\cdot
t},$ (25)

which is divergent at $t=0$ for large spacial coordinates and has quick
$t^{-1}$ decay in time in general. Figure (8) shows the graph of Eq. (25) for
a given parameter set.

Figure 8: The projection of Eq. (25) $u(x,y=0,t)$ for the parameter set of
$(\rho,m,\sigma,B_{0},c_{1})$ equal to $(1,0.1,1,0)$, respectively.

For $c_{0}\neq 0$, there is a formal implicit solution containing an integral:

$\eta-\int^{f(\eta)}\frac{2m}{[-6\rho m^{2}(a^{2}{\rho+a^{2}\sigma
B_{0}^{2}+2c_{0}})]^{1/3}}da-c_{1}=0.$ (26)

Unfortunately, for general arbitrary parameters
$(\rho,m,\sigma,B_{0}^{2},c_{0})$ it cannot be evaluated.

The second case is $n=1,$ the Newtonian fluid case:

$\displaystyle f(\eta)=C_{1}\cdot M\left(\frac{\rho}{2}+\sigma
B_{0}^{2},\frac{1}{2},-\frac{\rho[\eta-2c_{1}]^{2}}{4m}\right)+C_{2}\cdot
U\left(\frac{\rho}{2}+\sigma
B_{0}^{2},\frac{1}{2},-\frac{\rho[\eta-2c_{1}]^{2}}{4m}\right).$ (27)

Note, the interesting feature that both the density and the magnetic induction
give contributions into the first parameter of Kummer’s functions. The density
$\rho$ and the square of the magnetic induction $B_{0}^{2}$ are always
positive the electric conductivity $\sigma$ is also positive for regular
materials. (With the need for completeness we have to note that there are so-
called metamaterials which can have negative electrical conductivity but we
skip that case now. More on materials can be found in meta . For such media
the first parameter of Kummer’s function would be negative which would mean
divergent velocity fields, which contradicts energy conservation.) For regular
materials, the external magnetic field enhances the first positive parameter
of Kummer’s functions which means more oscillations and quicker decay. The
shape function or the final velocity distribution $u(x,y,t)$ is very similar
which was presented in Fig.(2) and Fig.(3). If the first parameter of Kummer’s
M function $(\rho/2+\sigma B_{0}^{2})$ is larger than one then the larger the
magnetic induction the larger the number of oscillations. This is true for the
density and for the electric conductivity as well. Such shape functions are
presented in fig (9).

Figure 9: Five shape functions of Eq. (27) for different parameter sets. The
effect of the magnetic induction is investigated. The black, red, blue and
green lines are for $(C_{1},C_{2},c_{1},\rho,m,\sigma,B_{0}^{2})$ with
numerical values of $(1,0,0,1,1,1,0.1)$, $(1,0,0,1,1,1,0.5)$,
$(1,0,0,1,1,1,1)$, $(1,0,0,1,1,1,2)$, and $(1,0,0,1,1,1,4)$, respectively

.

The third case is $n=0$ which is the simplest one. The original ODE of Eq. (9)
is reduced to

$\eta f^{\prime}+\sigma B_{0}^{2}f=0,$ (28)

with the trivial solution of:

$f(\eta)=C_{1}\eta^{-\sigma B_{0}^{2}}.$ (29)

The velocity distribution has the form

$u(x,y,t)=C_{1}\left(\frac{(x+y)}{t}\right)^{-\sigma B_{0}^{2}}.$ (30)

The pressure field is a pure constant as well. These are simple power laws,
therefore, we skip to present additional figures.

The last case, $n=-1$, is again the most complicated one. There is no analytic
solution available when all parameters are arbitrary. There is however a
special case when the last term in Eq. (9) is zero which defines the
constraint of $\rho=2\sigma B_{0}^{2}$. If we additionally fix $c_{1}=0$ we
get the following solution:

$f(\eta)=\pm\frac{{1}}{3}\sqrt{\frac{6m}{3\rho}}\cdot
arctan\left(\frac{\sqrt{6m\rho}\eta}{\sqrt{-6m\rho\eta^{2}+2mC_{1}}}\right)+C_{2}.$
(31)

It can be easily shown that the solution has a compact support, and the
function is only defined in the region of
$-\sqrt{\frac{C_{1}}{3\rho}}\leq\eta\leq\sqrt{\frac{C_{1}}{3\rho}}$. Note,
that the smaller the fluid density the larger the acceptable velocity range if
$c_{1}$ remains the same. The larger the integral constant $c_{1}$ parameter
the larger the available velocity range. The second integral constant $C_{2}$
just shifts the solution parallel to the $y$ axis. Figure (10) shows the
solution for three different parameter sets. Figure (11) presents the final
typical velocity distribution $v(x,y=0,t)\sim
t^{\frac{1}{2}}f(\frac{x}{t^{3/2}})$ the quick temporal decay is clear to see.
Note, that the prefactor $t^{\frac{1}{2}}$ makes the time backpropagation
impossible for negative values.

Figure 10: The velocity shape functions of Eq. (31) for three different
parameter sets $(C_{1},C_{2},m,\rho)$ with $c_{1}=0$. The black, red and blue
curves are for $(1,0,1,1),(15,0.2,2,1)$ and $(0.4,2,1,0.4)$, respectively

.

Figure 11: The velocity distribution which corresponds to (31) for the
parameter set $(C_{1},C_{2},\rho,m)$ is $(0,0,5,3)$

.

For the sake of completeness, we give the solutions for the pressure as well.
The ODE of the shape function is trivial with the solution of:

$h^{\prime}=0,\hskip 28.45274pth=c_{2}.$ (32)

Therefore, the final pressure distribution reads:

$p(x,y,t)=t^{-\gamma}\cdot h(x,y,t)=\frac{c_{2}}{t^{n}},$ (33)

which means that the pressure is constant in the entire space at a given time,
and the time decay can be different to the velocity field.

###  II.2 Casson fluid

To perform an even more comprehensive analysis to study how the transients
happen in non-Newtonian boundary layers we investigated additional fluid
models.

The first in our line is the simplest one, the so-called Casson fluid model,
now the boundary layer equation has the form of

$\rho\left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}+v\frac{\partial u}{\partial y}\right)=-\frac{\partial p}{\partial
x}+\left(1+\frac{1}{\lambda}\right)\mu\frac{\partial^{2}u}{\partial
x^{2}}+\sigma B_{0}^{2}u,$ (34)

where $\lambda$ is the Casson parameter. Note that for
$\lambda\rightarrow\infty$ the term goes over the regular Newtonian viscous
term. (All five physical parameters $\rho,\lambda,\mu,B_{0}$ and $\sigma$
should should have positive real values and the $\lambda\neq 0$ is also
evident.) The variational approach for the flow of Casson fluids in pipes was
analyzed by Sochi sochi . The blood flow through a stenotic tube was modelled
by the Casson fluid by Tandon et al. tandon .

We still consider the self-similar Ansatz of (4) for the three dynamical
variables. After the usual algebraic manipulations, we arrive at the non-
linear second-order ODE for the shape function of horizontal velocity
component $f(\eta)$ in the next form:

$\rho\left(-\frac{f}{2}-\frac{\eta
f^{\prime}}{2}\right)=\left(1+\frac{1}{\lambda}\right)\mu
f^{\prime\prime}+\sigma B_{0}^{2}f,$ (35)

in this model, all four self-similar exponents have fixed values,
$\alpha=\beta=\delta=1/2,\gamma=1$, which is usual for the Newtonian viscous
fluid equations imre4 . The derived solutions are mainly different to the
first model, therefore we have to give a detailed analysis in the following.
With our usual mathematical program package Maple 12, we can easily derive the
solution in closed form:

$\displaystyle
f(\eta)=e^{-\frac{\rho\eta^{2}}{4\mu\left(\frac{\lambda+1}{\lambda}\right)}}\cdot\eta\cdot\left(C_{1}M\left[-\frac{2\sigma
B_{0}^{2}-\rho}{2\rho},\frac{3}{2},\frac{\rho\eta^{2}}{4\mu\left\\{\frac{\lambda+1}{\lambda}\right\\}}\right]\right.+$
$\displaystyle\left.C_{2}U\left[-\frac{2\sigma
B_{0}^{2}-\rho}{2\rho},\frac{3}{2},\frac{\rho\eta^{2}}{4\mu\left\\{\frac{\lambda+1}{\lambda}\right\\}}\right]\right).$
(36)

It is important to note at this point, that the derived formula shows some
similarities to the results of the regular diffusion equations laci_roman
which was exhaustively discussed in our former studies. The situation is
however a bit different here. Thanks to the positivity of all parameters, the
sign of exponent is always negative (which is a real Gaussian function) which
dictates a rapid decay for large arguments $\eta$ and for any kind of
additional parameter sets. The crucial parameter which qualitatively
classifies the solutions is the numerical value of the first parameter of
Kummer’s functions. We can distinguish three cases:

  * •

$-\left(\frac{2\sigma B_{0}^{2}-\rho}{2\rho}\right)<0$ which is equivalent to
$2\sigma B_{0}^{2}>\rho$ (with the stipulation of $\rho\neq 0$) the solutions
have oscillations. If this parameter is a negative integer the infinite series
of the Kummer’s function breaks down to a finite order polynomial. The smaller
the parameter the larger the number of oscillations.

  * •

$-\left(\frac{2\sigma B_{0}^{2}-\rho}{2\rho}\right)=0$ which is equivalent to
$2\sigma B_{0}^{2}=\rho$ (with the stipulation of $\rho\neq 0$) now both
Kummer’s M and Kummer’s U functions are unity and the solution is reduced to
the Gaussian function times $\eta$ function which has odd symmetry. This is
the limiting solution between the oscillating and the non-oscillating
solutions.

  * •

$-\left(\frac{2\sigma B_{0}^{2}-\rho}{2\rho}\right)>0$ which is equivalent to
$2\sigma B_{0}^{2}<\rho$ (with the stipulation of $\rho\neq 0$) the larger the
parameter the smaller peak value of the solution and the quicker the decay.

Figure (12) presents five shape functions of Eq. (36) for five different fluid
densities, for a clear comparison all the other parameters remain the same.
The larger the density $\rho$ the quicker the decay and the smaller the global
maximum of the shape functions.

Figure 12: Five shape functions of Eq. (36) for different parameter sets. The
role of the fluid density $\rho$ is investigated. The black, red, blue green
and grey lines are for $(C_{1},C_{2},\rho,\lambda,\mu,\sigma,B_{0}^{2})$ with
numerical values of
$(1,0,8,1,1,1,1),(1,0,4,1,1,1,1),(1,0,2,1,1,1,1),(1,0,1,1,1,1,1)$ and
$(1,0,1/2,1,1,1,1)$, respectively

.

Figure 13: Five shape functions of Eq. (36) for different parameter sets. The
role of the magnetic induction $B_{0}^{2}$is investigated. The black, red,
blue and green lines are for $(C_{1},C_{2},\rho,\lambda,\mu,\sigma,B_{0}^{2})$
with numerical values of
$(1,0,1,1,1,1,0.1),(1,0,1,1,1,1,0.5),(1,0,1,1,1,1,1),(1,0,1,1,1,1,2)$ and
$(1,0,1,1,1,1,4)$ respectively

.

In general, the ${\rho}/[{4\mu(\lambda+1)/\lambda]}$ factor is responsible for
the extent (or the full with at half maximum (FWHM) value of the solution) the
larger this parameter the quicker the decay. For changing the Casson parameter
in the physically relevant positive range $\lambda>0$ we found no remarkable
change in the solutions, this is because the $\frac{\lambda+1}{\lambda}$ (is
almost unit) factor appears in the Gaussian together with two other
parameters.

The role of the magnetic induction is also important to know. Figure (13)
shows the effect of the magnetic induction $B_{0}^{2}$ where all the other
parameters and unchanged. If the first parameter of Kummer’s M function is
smaller than unity the solutions has a global maxima and a quick decay, if
this parameter is larger than the unity then the larger the magnetic induction
the larger the oscillations of the solutions.

Most of these solutions have odd symmetry. It can be easily seen because the
Gaussian is an even function $\eta$ is an odd function and the series of
Kummer’s M and U functions which have a quadratic argument are again even
infinite or finite series which together given odd symmetry. However, using
the series expansion it can be shown that if the first parameter of the
Kummer’s U (only for U ) function $-\left(\frac{2\sigma
B_{0}^{2}-\rho}{2\rho}\right)$ is a negative half-integer then the solutions
get even symmetry. The properties of such kinds of solutions are exhaustively
discussed in laci .

To complete our analysis we give the solutions for the pressure as well. The
ODE of the shape function is trivial with the solution of:

$h^{\prime}=0,\hskip 28.45274pth=c_{2}.$ (37)

Therefore, the final pressure distribution is:

$p(x,y,t)=t^{-1}\cdot h(x,y,t)=\frac{c_{2}}{t^{-1}},$ (38)

which means that the pressure is constant in the entire space at a given time,
and the temporal decay follows the simple inverse law.

###  II.3 Oldroyd-B model

The next fluid flow model is the Oldroyd-B model old , where the momentum
equation is the following:

$\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}+v\frac{\partial v}{\partial
y}+\frac{\lambda_{1}}{\rho}\left(u^{2}\frac{\partial^{2}u}{\partial
x^{2}}+2uv\frac{\partial^{2}u}{\partial x\partial
y}+v^{2}\frac{\partial^{2}u}{\partial
y^{2}}\right)=-\frac{1}{\rho}\frac{\partial p}{\partial
x}+\frac{\nu}{\rho}\frac{\partial^{2}u}{\partial y^{2}}$
$\displaystyle+\frac{\nu\lambda_{2}}{\rho}\left(u\frac{\partial^{3}u}{\partial
x\partial y^{2}}+v\frac{\partial^{3}u}{\partial y^{3}}-\frac{\partial
u}{\partial y}\frac{\partial^{2}u}{\partial y^{2}}-\frac{\partial u}{\partial
y}\frac{\partial^{2}v}{\partial y^{2}}\right)+\sigma B_{0}^{2}u.$ (39)

Even without the external term with the magnetic induction, the self-similar
Ansatz of (4) leads to contradiction among the four self-similar exponents
indicating that the system has no self-similar symmetry with power-law time
decay. The main cause of the lack of symmetry is the appearance of the third
spatial derivatives together with the second ones. The Oldroyd-B model was
exhaustively investigated in the last decades from many points of views,
optimal time decay rates for the higher order spatial derivatives of solutions
were analysed by Wang wang . The global wellposedness of the model was
investigated by Elgindi and Liu elg . The global existence results of some
Oldroyd-B models were proven by hieber . Finally, we mention the review paper
of rev  which summarizes most of the known mathematical results.

###  II.4 Walters’ Liquid B model

The third possible non-Newtonian candidate is the so-called Walters’ Liquid B
model, walt  where the momentum equation has the form of:

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial
v}{\partial y}=-\frac{1}{\rho}\frac{\partial p}{\partial
x}+\frac{\mu}{\rho}\frac{\partial^{2}u}{\partial
y^{2}}+\frac{k_{0}}{\rho}\left(u\frac{\partial^{3}u}{\partial x\partial
y^{2}}+v\frac{\partial^{3}u}{\partial y^{3}}-\frac{\partial u}{\partial
y}\frac{\partial^{2}u}{\partial x\partial y}+\frac{\partial u}{\partial
x}\frac{\partial^{2}u}{\partial y^{2}}\right)+\sigma B_{0}^{2}u.$ (40)

Practically, we observed the same as in the previous case, the existence of
the terms with the third derivatives destroyed the self-similar symmetry
therefore no solutions can be derived.

Flow and heat transfer of Walter’s liquid model with stretching walls were
applied to haemodynamics by Misra, Shit and Rath mis . Khrisna kris  studied
the Hall and ion slip effects on MHD laminar flow of in a Walter’s‐B fluid.

###  II.5 Williamson fluid

In the last model, we took the Williamson fluid wil  with the equation of
motion of:

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial
v}{\partial y}=-\frac{1}{\rho}\frac{\partial p}{\partial
x}+\frac{\partial}{\partial y}\left(\frac{\mu}{\rho}\frac{\partial u}{\partial
y}+\mu\frac{\Gamma}{2}\left(\frac{\partial u}{\partial
y}\right)^{2}\right)+\sigma B_{0}^{2}u.$ (41)

Even if this model lacks the property of time-dependent self-similar symmetry,
therefore we cannot present analytical results. Malik, et al. malik  however
numerically studied the effects of variable thermal conductivity and heat
generation/absorption on Williamson fluid flow and heat transfer. The swimming
effects of microbes in blood flow of nano-bioconvective Williamson’s fluid was
investigated by Rana et all. rana . We found these recent results which we
think are worth mentioning.

##  III Summary and Outlook

In our study, we investigated five different time-dependent non-Newtonian
boundary layer problems using the self-similar approach. We show that
analytical results exist and give them. For the non-Newtonian power-law fluid,
we found different types of solutions depending on the value of the power-law
exponent, which can be expressed by Kummer functions or other implicit
functions. The effects of density and other parameters in the nonlinear
ordinary differential equation have been investigated. The effects of the
magnetic field can be described by including an additional magnetic term in
the equation. Analytical solutions have also been found for such cases. For
non-Newtonian Casson fluid flow, the effect of the Casson parameter is also
analyzed. For the other three non-Newtonian cases, the Oldroyd-B model,
Walter’s Liquid B model, and Williamson fluid, the self-similar transformation
(4) cannot be performed, for these cases the solution according to the power
law and the Casson fluid cannot be given. All five models might be
investigated with the travelling wave Ansatz in the far future because the
spatial and temporal symmetry shift is always present in these equations.

##  IV Acknowledgments

One of us (I.F. Barna) was supported by the NKFIH, the Hungarian National
Research Development and Innovation Office. This work was supported by project
no. 129257 implemented with the support provided by the National Research,
Development and Innovation Fund of Hungary, financed under the $K\\_18$
funding scheme.
Author Contributions: I.F.B. Conceptualization, analytic calculations,
manuscript writing; G.B. Conceptualization, Literature; K.H. Literature
collection, writing the final manuscript; L.M. Correction of the manuscript
Data Availability Statement: The data that supports the findings of this study
are all available within the article.
Conflicts of Interest: The authors declare no conflict of interest.

## References

  * (1) G. Astarita and G. Marucci. Principles of non-Newtonian fluid mechanics, McGraw-Hill, London, New York, 1974.
  * (2) I. Fridtjov Rheology and Non-Newtonian Fluids, Springer, 2014.
  * (3) M. Patel and M. Timol, Non-Newtonian fluid models and Boundary Layer flow , LAP LAMBERT Academic Publishing, 2020.
  * (4) H. Schlichting and K. Gersten Boundary-layer theory Berlin Heidelberg New York: Springer; 2016.
  * (5) Y. Hori Hydrodynamic lubrication Tokyo: Springer; 2006.
  * (6) I. F. Barna, G. Bognár, L. Mátyás and K. Hriczó, Journal of Thermal Analysis and Calorimetry 147, (2022) 13625.
  * (7) C. Saengow and A.J. Giacomin, Phys Fluids. 30, (2018) 030703.
  * (8) G. Bognár, Int J Nonlinear Sci Numer Simul. 10, (2009) 1555.
  * (9) G. Bognár and K. Hriczó, Acta Polytechnica Hungarica, 8, (2011) 131.
  * (10) T. M. Ajayi, A. J. Omowaye, and I. L. Animasaun, Journal of Applied Mathematics 2017, Article ID 1697135.
  * (11) N. Casson,  Rheology of Disperse Systems (Edited by Mill, D. C.) pp. 84–102. Pergamon Press, Oxford, 1959.
  * (12) J. Oldroyd, Proc. Roy. Soc. London Series A, Math. and Phys. Sciences 200, (1950) 523.
  * (13) K. Walters, Quart. J. Mech. Appl. Math. 6, (1963) 63.
  * (14) R.V. Williamson, Ind. Eng. Chem. 21, (1929) 1108.
  * (15) S. Ko, J. Math. Phys. 63, (2022) 041508.
  * (16) R.S. Herbst, C. Harley, K.R. Rajagopal, International Journal of Non-Linear Mechanics, 154, (2023) 104445.
  * (17) R. Fazio, Calcolo, 59, (2022) 43.
  * (18) M. Patel, H. Surati and M. G. Timol, Math. J. Interdiscip. 9, (2021) 35.
  * (19) L. Sedov, Similarity and Dimensional Methods in Mechanics, CRC Press, 1993.
  * (20) I.F. Barna, M.A. Pocsai, S. Lökös and L. Mátyás, Chaos Solitons and Fractals 103, (2017) 336.
  * (21) I.F. Barna, L. Mátyás and M.A. Pocsai, Fluid. Dyn. Res. 52, (2020) 015515.
  * (22) I.F. Barna, Commun. Theor. Phys. 56, (2011) 745.
  * (23) F. W. J. Olver, D. W. Lozier, R. F. Boisvert and C. W. Clark NIST Handbook of Mathematical Functions Cambridge University Press, 2010.
  * (24) F. Cappolino Theory and Phenomena of Metamaterials CRP Press, 2009.
  * (25) L. Mátyás and I.F. Barna Romanian Journal of Physics 67, (2022) 101.
  * (26) L. Mátyás and I.F. Barna, Universe 9, (2023) 264.
  * (27) T. Sochi, International Journal of Modeling, Simulation, and Scientific Computing, 7, (2016) 1650007.
  * (28) P.N. Tandon, U.V.S. Rana, M. Kawahara, V.K. Katiyar, International Journal of Bio-Medical Computing, 32, (1993) 61.
  * (29) Y. Wang, Z. Angew. Math. Phys. 74, (2023) 3.
  * (30) T.M. Elgindi, J.L. Liu, Differential Equations, 259, (2015) 1958
  * (31) M. Hieber, Y. Naito, Y. Shibata, J. Differential Equations, 252, (2012) 2617.
  * (32) M. Renardy, B. Thomases, J. Non-Newton. Fluid Mech., 293, (2021) 104573.
  * (33) J. C. Misra, G. C. Shit and H. J. Rath, Computers & Fluids, 37, (2008) 1.
  * (34) M. Veera Krishna, Heat Transfer 49, (2020) 2311.
  * (35) M.Y. Malik, M. Bibi, F. Khan and T. Salahuddin, AIP Advances, 6, (2016) 035101.
  * (36) B.M.J. Rana, S.M. Arifuzzaman, Saiful Islam, Sk. Reza-E-Rabbi, Abdullah Al-Mamun, Malati Mazumder, Kanak Chandra Roy, Md. Shakhaoath Khan, Thermal Science and Engineering Progress, 25, (2021) 101018.


# Astrophysical Bounds on Mirror Dark Matter, Derived from Binary Pulsars
Timing Data

Itzhak Goldman1,2
1 Department of Physics, Afeka Engineering College
Tel Aviv, Israel
<EMAIL_ADDRESS>
2 Department of Astrophysics, Tel Aviv University
Tel Aviv, Israel

###### Abstract

Mirror Dark Matter (MDM) has been considered as an elegant framework for a
particle theory of Dark Matter (DM). It is supposed that there exists a dark
sector which is mirror of the ordinary matter. Some MDM models allow particle
interactions mirror and ordinary matter, in addition to the gravitational
interaction. The possibility of neutron to mirror neutron transition has
recently been discussed both from theoretical and experimental perspectives.
This paper is based on a previous work in which we obtained stringent upper
limits on the possibility of converting neutrons to mirror neutrons in the
interiors of neutron stars, by using timing data of binary pulsars. Such a
transition would imply mass loss in neutron stars leading to a significant
change of orbital period of neutron star binary systems. The observational
bounds on the period changes of such binaries, therefore put strong limits on
the above transition rate and hence on the neutron – mirror-neutron mixing
parameter $\epsilon^{\prime}$. Our limits are much stronger than the values
required to explain the neutron decay anomaly via $n-n^{\prime}$ mixing.

###### keywords:

dark matter, neutron stars, pulsars

95.35.+d, 97.60.Jd, 97.60.Gb

##  1 Introduction

The idea of mirror dark matter (MDM) has been considered by various authors
[1]. In this class of models each ordinary particle $p$ has a mirror particle
$p^{\prime}$ counterpart. There is, a generic problem in that the extra mirror
neutrinos ($\nu^{\prime}$) and mirror photon ($\gamma^{\prime}$) contribute
too much to the number of degrees of freedom at the Big Bang Nucleosynthesis
(BBN) epoch, destroying the success of the big bang nucleosynthesis
predictions. A possible way out is to assume a breaking of the $Z_{2}$
symmetry in the early universe so as to have asymmetric inflationary reheating
in the two sectors (ordinary matter and MDM) resulting in a lower reheat
temperature ($T^{\prime}$) of the mirror sector compared to ( $T$) of the
visible one [4]. This breaking eventually trickles down to the low energies
leading in general to a splitting the mirror and visible fermion masses [5].
The above mentioned symmetric picture could however remain almost exact if the
asymmetric inflation picture is carefully chosen. Cosmologies of such
scenarios have been discussed in [2] and references therein

.

An interesting new phenomenon is possible in almost exact mirror models, if
there are interactions mixing the neutron with the mirror neutron state (
denoted by $\epsilon_{n-n^{\prime}}\equiv\epsilon^{\prime}$). In such a case,
one can expect $n\to n^{\prime}$ oscillations to take place in the laboratory
[7] and indeed there are ongoing and already completed searches for such
oscillations [3] at various neutron facilities.

We note that $n-n^{\prime}$ oscillation is similar to neutron-anti-neutron
oscillation suggested very early [9] and extensively discussed in the
literature; For recent reviews, see [10]. The rate governing the $n-\bar{n}$
oscillations was limited by laboratory experiments to be

$\displaystyle\epsilon_{n\bar{n}}=1/{\tau_{n\bar{n}}}<10^{-8}~{}Sec^{-1}~{}{\rm
or}~{}10^{-23}~{}{\rm eV}$ (1)

We focus on the possibility that $n\to n^{\prime}$ transitions can lead to
mass loss of a neutron star and if the latter is a member of a binary s
system, then this mass loss affects the binary period. Indeed [12] the mass
loss of any kind of such neutron stars implies an increase of the orbital
period of the binary system $P_{b}$. In the following sections we find that
observational limits on $|\dot{P}_{b}/P_{b}|$ of binary pulsars yield
stringent bounds on the $n-n^{\prime}$ mixing parameter $\epsilon^{\prime}$.

##  2 Transition of a neutron star to a mixed neutron-mirror neutron star
induced by $n-n^{\prime}$ mixing

We note that $n\to n^{\prime}$ transitions, kinematically forbidden in nuclei,
can occur in neutron stars. The neutrons in neutron stars are mainly bound by
gravity and _not_ by nuclear forces. Let us, suppose that an $n\to n^{\prime}$
conversion occurred at some point in the star. Under the pressure a
neighboring neutron then will be pushed to ”hole” generated by the converted
neutron, gaining in the process kinetic energy which is of order of the Fermi
energy $E_{F}$. Additional energy gain is obtained since the produced
$n^{\prime}$ gravitates to the center of the star and a surface neutron
replaces the neutron which went into the above mentioned ”hole”. Therefore the
net (eventually radiated via neutrino and mirror neutrino emission) energy
stemming from an early, single $nn^{\prime}$ transition is:

$m_{n}c^{2}\left(e^{\phi(R)}-e^{\phi(0)}\right)+<E_{F}>$ (2)

where $e^{\phi(r)}=g_{oo}^{1/2}(r)$, $e^{\phi(R)}=(1-2GM/c^{2})^{1/2}$ and
$<E_{F}>$ is the average Fermi energy of the disappeared neutron. Thus the
gravitational mass of the neutron star will decrease. This has led [13] to
propose that $(n\to n^{\prime})$ oscillations generating completely mixed
$nn^{\prime}$ stars, may explain the observed mass distribution of neutron
stars.

##  3 The rate of $n\rightarrow n^{\prime}$ transition

We express the $n\rightarrow n^{\prime}$ transition rate
$\Gamma(n\rightarrow{n^{\prime}})$ by

$\Gamma(n\rightarrow{n^{\prime}})=\Gamma(nn)P_{nn^{\prime}}$ (3)

where $\Gamma_{nn}$ is the rate of $nn$ collisions and $P_{nn^{\prime}}$ is
the probability of having $n^{\prime}$ (rather than $n$) at the time of the
collision.

Eq.(3) is based on the assumption [14] in which one assumes that:

a) The coherent buildup of the $|n^{\prime}>$ component in the initial purely
$|n>$ state of the two component system, proceeds unimpeded by nuclear
interactions during the time of flight between two consecutive collisions.

b) The coherent build-up stops upon collision and the $n^{\prime}$ part is
released as out-going mirror neutron particles.

Denoting by $t_{nn}$ the mean free time for neutron neutron collisions and by
$\epsilon^{\prime}$ the rate for $n$ to $n^{\prime}$ oscillation, the
Hamiltonian in the two dimensional $|n>,|n^{\prime}>$ Hilbert space leads to

$P_{nn^{\prime}}=[\epsilon^{\prime}\cdot t_{nn}]^{2})$

Substituting the above $P_{nn^{\prime}}$ and $\Gamma_{nn}=t_{nn}^{-1}$ in Eq.
(3) one gets

$\displaystyle\Gamma_{n\to n^{\prime}}=t_{nn}{\epsilon^{\prime}}^{2}$ (4)

Taking $t_{nn}\approx{10^{-23}}~{}sec$, the flight time of a neutron a
$O(Fermi)$ distance at a speed $\sim 1/3c$ yields

$\displaystyle\Gamma_{n\to{n^{\prime}}}=6\times
10^{-8}[\epsilon^{\prime}/{10^{-11}eV}]^{2}~{}{yr}^{-1}$ (5)

In what follows we shall use astrophysical data to obtain bounds on
$\Gamma_{n\to{n^{\prime}}}$ and employ Eq. (5) to derive bounds on
$\epsilon^{\prime}$.

##  4 Neutron star models and their descendant fully mixed neutron - mirror
neutron stars

To obtain the resulting mass and radius decrease, we solved numerically the
TOV equations with a commonly used equation of state [21] . We first solve the
TOV structure equations for pure neutrons. Then for the same baryon number we
solve the TOV equations for a totally mixed neutron - mirror neutron star. We
employ the method used by us in [6]. In this specific case we have two fluids
obeying the same nuclear equation of state where the only interaction is
through the gravitational field.

We considered 3 different models each characterized by its total baryon number
(that does not change in the transition). For each baryon mass we calculated
the initial and final mass and radius. We also calculated the initial
$e^{\phi(R)}-e^{\Phi(0)}$. The resulting models are:

  * $M_{b}=1.81M_{\odot}$

$\ \ \ M=1.57M_{\odot},\ \ \ \ R=12.2km,\ \ \ e^{\phi(R)}-e^{\Phi(0)}=0.15$

$\ \ \ M_{nn^{\prime}}=1.43M_{\odot}=0.91M\ \ \ \ \ \
R_{nn^{\prime}}=8.8km=0.72R$

  * $M_{b}=2.17M_{\odot}$

$\ \ \ M=1.82M_{\odot}\ \ \ \ R=12.37km,\ \ \ e^{\phi(R)}-e^{\Phi(0)}=0.17$

$M_{nn^{\prime}}=1.62M_{\odot}=0.89M\ \ \ \ R_{nn^{\prime}}=8.8km=0.71R,\ \ \
$

  * $M_{b}=2.44M_{\odot}$

$\ \ M=1.97M_{\odot},\ \ \ \ R=12.4,\ \ \ e^{\phi(R)}-e^{\Phi(0)}=0.19$

$M_{nn^{\prime}}=1.7M_{\odot}=0.86M\ \ \ \ R_{nn^{\prime}}=8.65km=0.7R,\ \ \ $

##  5 The relation between the mixing parameter $\epsilon^{\prime}$ and the
mass loss rate

Returning to the main goal of the paper, namely limiting $\epsilon^{\prime}$,
we need to relate the rate of neutron to mirror neutron transition to the
stellar $\frac{\dot{M}}{M}$. Since neutrons from the entire volume of the
neutron star can transform to $n^{\prime}$, the resulting total mass loss of
the neutron star due to this is proportional to $\Gamma_{n\to n^{\prime}}$.

A great advantage is that all the different binary pulsars (with variety of
masses, spin down ages and companions ) should conform to a single fundamental
parameter: $\epsilon^{\prime}$. Thus, we can use the youngest pulsars to set
stringent limits on $\epsilon^{\prime}$. In turn, this implies that also the
older pulsars are still in the process of transition. In particular, for small
enough values of $\epsilon^{\prime}$ even these older pulsars may be at the
very initial stages of the pure to mixed star transition.

We next present two estimates of $\frac{\dot{M}}{M}$.

First Estimate

The first estimate uses the average value
$\frac{\dot{M}}{M}\approx\Gamma_{nn^{\prime}}\frac{\Delta M}{M}$ during the
complete transition to a mixed star where, $\Delta M$ is the total mass
reduction. The $\frac{\Delta M}{M}=0.09\div 0.14$ obtained in the previous
section then yields a representative value

$\left|\frac{\dot{M}}{M}\right|\approx 0.12\Gamma_{nn^{\prime}}$ (6)

Second Estimate

Here we use the results of the numerical solutions of the neutron stars
presented in section 4. We apply Eq. (2) to the very beginning of the
transition process and find

$\left|\frac{\dot{M}}{M}\right|\geq\left(e^{a}nd{\phi(R)}-e^{\phi(0)}\right)\Gamma_{nn^{\prime}}=(0.15\div
0.19)\Gamma_{nn^{\prime}}$ (7)

It reassuring that the estimates by the two methods agree up to a factor of
1.5. It also makes sense that at the very beginning the process is faster than
that derived from the first method which represents a time average over the
entire transition.

Thus we adopt

$\left|\frac{\dot{M}}{M}\right|\geq 0.14\Gamma_{nn^{\prime}}$ (8)

as a representing value during the entire transition episode.

##  6 Limits on $\epsilon^{\prime}$ derived from timing observations of binary
pulsars

Jeans (1924) [11], pointed out that the mass of the star that emits
electromagnetic radiation decreases with time and therefore the orbital
elements of a binary system should evolve with time. Assuming that in the
local frame of each star the radiation emission is spherically symmetric, he
obtained

$Ma=\ {\rm constant}$ (9)

where $a$ is the semi-major axis and $M=m_{1}+m_{2}$ is the total mass of the
system. This and the expression for the binary period

$P_{b}=2\pi\sqrt{\frac{a^{3}}{GM}}$ (10)

imply

$\frac{\dot{P}_{b}}{P_{b}}=-2\frac{\dot{M}}{M}$ (11)

Since $\dot{M}<0$, $\dot{P}_{b}>0$ so that the orbital period keeps
increasing.

Except for the nature of emission process this is the same as the situation
addressed here. We consider next four binary pulsars and obtain the the limits
on the present rate of the mass change for each system.

In general, only neutron stars which are still in the process of transiting to
a mixed star will exhibit mass loss. A priory, for extremely ”high”
$\Gamma_{nn^{\prime}}$ some of the older observed pulsars in binaries may be
”too old”.

PSR 1916+13

We use the data of [15]. The spin down age of the pulsar is $1.1\times
10^{8}yr$. It is commonly assumed that the pulsar companion is a neutron star,
so that there is no mass transfer between the binary members. Also one
neglects the mass accretion from the interstellar medium. After accounting for
the expected gravitational radiation, Galactic acceleration, and further
dynamical corrections there is still some limited room for a positive change
of binary period that could have followed from the mass decrease.

$\frac{\dot{P}_{b}}{P_{b}}<1.36\times 10^{-11}yr^{-1}$ (12)

implying

$\left|\frac{\dot{M}}{M}\right|<6.8\times 10^{-12}yr^{-1}$

for the two neutron stars losing mass. Using Eqs. (5) and (14) the bound

$\epsilon^{\prime}<2.8\times 10^{-13}eV$

follows.

PSR J1141-6545

This is a young pulsar of age about $2\times 10^{6}yr$. The companion is a
massive white dwarf of mass $0.98M_{\odot}$ and the neutron star has a mass
$1.3M_{\odot}$ [19, 20]

This system is a superb laboratory for testing GR. Using the data from the
above observational papers one finds that the residual (subtracting from the
measured value the gravitational radiation terms as well as the galactic
acceleration and the kinematic effect and allowing for the uncertainties of
all the above) positive possible value of the orbital period rateThe mirror
neutrinos ($\nu^{\prime}$) and mirror photon ($\gamma^{\prime}$) contribute
too much to the number of degrees of freedom at the Big Bang Nucleosynthesis
(BBN of change is very small

$\frac{\dot{P}_{b}}{P_{b}}<1.84\times 10^{-12}yr^{-1}$

implying for the given masses and taking into account that only the neutron
star is undergoing the mass loss

$\left|\frac{\dot{M}}{M}\right|=\frac{M+M_{c}}{2M}\frac{\dot{P}_{b}}{P_{b}}<1.6\times
10^{-12}yr^{-1}$

where $M_{c}$ is the mass of the white dwarf companion.

This implies a limit

$\epsilon^{\prime}<1.4\times 10^{-13}eV$

This result is very important, as for the other three pulsars pulsars, one may
have argued that the the transition $n\to n^{\prime}$ has already finished.
This young pulsar closes the door on this argument.

PSR J0437-4715

This is a neutron star- white dwarf binary with masses of $1.76M_{\odot}$ and
$0.25M_{\odot}$, respectively. The spin-down age is $1.6\times 10^{9}\ yr$.
Using the data from the above observational papers one finds that the residual
(subtracting from the measured value the gravitational radiation terms as well
as the galactic acceleration and the kinematic effect and allowing for the
uncertainties of all the above) positive possible value of the orbital period
rate of change is [20]

$\frac{\dot{P}_{b}}{P_{b}}<2.8\times 10^{-11}yr^{-1}$

implying

$\left|\frac{\dot{M}}{M}\right|<1.6\times 10^{-11}yr^{-1}$

and therefore

$\epsilon^{\prime}<4.4\times 10^{-13}eV$

PSR J1952+2630

This pulsar is in a binary orbit with a $(0.93\div 1.4)M_{\odot}$ white dwarf
companion [18] . Its spin down age is $7.7\times 10^{7}$ yr, the orbital
period is $0.39$ days and during 800 days of follow-up the error on the period
is $7\times 10^{-12}$ days. This leads to $\frac{\dot{P}_{b}}{P_{b}}<8.2\times
10^{-12}yr^{-1}$.

Taking into account that only the pulsar losses mass we get

$\left|\frac{\dot{M}}{M}\right|=\left(\frac{M+M_{c}}{2M}\right)\frac{\dot{P}_{b}}{P_{b}}<7\times
10^{-12}yr^{-1}$

where $M_{c}$ is the mass of the white dwarf companion.

Thus one finds

$\epsilon^{\prime}<2.9\times 10^{-13}eV$

The expected period change due to gravitational radiation in this pulsar is 2
orders of magnitude smaller, and thus does not interfere with the derived
limit.

##  7 Discussion

In this paper, used astrophysical data, notably, precision pulsar timing
measurements to strongly constrain a putative $n\to n^{\prime}$ tThe mirror
neutrinos ($\nu^{\prime}$) and mirror photon ($\gamma^{\prime}$) contribute
too much to the number of degrees of freedom at the Big Bang Nucleosynthesis
(BBNransition in neutron stars.

We solved numerically the geneandral relativistic structure equations for
neutron stars with three different baryon mass. Using a realistic nuclear
matter equation of stat, we first solved for a pure ordinary neutron star and
then (for the same baryon mass) we obtained the solution for a fully mixed
neutron- mirror neutron star. In this way we found the average mass reduction
rate over the the time span of the transition. From the solution of the
ordinary neutron star we found the mass reduction at the very beginning of the
process. We found that two methods yield similar relations between the rate of
mass change of the star and the rate of the microscopic process.

This allows us to restrict the mixing parameter between. We find that the key
parameter $\epsilon^{\prime}$ responsible for $n\to n^{\prime}$ transition is
restricted to be below $4.4\times 10^{-13}$ eV.

All four binary pulsars considered here yield quite similar limits on
$\epsilon^{\prime}$ . Of particular importance is the the limits are quite the
same even though systems ages vary between $2\times 10^{6}\ yr$ and $1.6\times
10^{9}yr$. The presence of the youngest binary pulsar refutes the possibility
that for the older pulsars the process has been already finished.

Our limits on $\epsilon^{\prime}$ also exclude the possibility that
$nn^{\prime}$ oscillations can explain the neutron decay anomaly [22]

## Acknowledgments

This work has been supported by the Afeka College Research Fund.

## References

  * [1] Lee, T. D. and Yang,C. N. Phys. Rev. 104, 254 (1956) communication; Kobzarev, y. ,Okun , L., and Pomeranchuk, I. Ya, Yad. Fiz. 3, 1154 (1966); Pavsic, M. Int. J. T. P. 9, 229 (1974); Kolb, E. W. , Seckel, D. and Turner, M. Nature, 514, 415 (1985); Foot, R. , Lew, H. and Volkas, R. Phys. Lett. B 272, 67 (1991);
  * [2] Foot, R. Mirror dark matter: Cosmology, galaxy structure and direct detection Int. J. Mod. Phys. A 29, 1430013 (2014)
  * [3] Broussard , L. and Y. Kamyshkov, Talk at BLV 2017; Berezhiani, Z., Frost, M. ,Kamyshkov, Y. ,Rybolt, B. and Varriano,L., Neutron Disappearance and Regeneration from Mirror State Phys. Rev. D 96, 035039 (2017).
  * [4] Berezhiani, Z., Dolgov A. and R. N. Mohapatra, R. N., Asymmetric inflationary reheating and the nature of mirror universe, Phys. Lett. B 375, 26 (1996).
  * [5] Mohapatra, R. N.and Nussinov, S., Constraints on Mirror Models of Dark Matter from Observable Neutron-Mirror Neutron Oscillation, Phys. Lett. B 776, 22 (2018)
  * [6] I. Goldman, R. N. Mohapatra, S. Nussinov, D. Rosenbaum and V. Teplitz, Possible Implications of Asymmetric Fermionic Dark Matter for Neutron Stars Phys. Lett. B 725, 200 (2013)
  * [7] Berezhiani, Z. and Bento, L., Neutron - mirror neutron oscillations: How fast might they be?, Phys. Rev. Lett. 96, 081801 (2006); Berezhiani, Z. More about neutron - mirror neutron oscillation Eur. Phys. J. C 64, 421 (2009).
  * [8] Goldman I., Mohapatra R. N., Nussinov S., bounds on neutron mirror neutron mixing from pulsar timing, PhRvD, 100, 123021,(2019)
  * [9] Kuzmin, V. A. Cp violation and baryon asymmetry of the universe, Pisma Zh. Eksp. Teor. Fiz. 12, 335 (1970); Glashow, S,L. The Future of Elementary Particle Physics, NATO Sci. Ser. B 61, 687 (1980); Mohapatra, R. N. and R. E. Marshak,R. E. Local B-L Symmetry of Electroweak Interactions, Majorana Neutrinos and Neutron Oscillations, Phys. Rev. Lett. 44, 1316 (1980) Erratum: [Phys. Rev. Lett. 44, 1643 (1980)].
  * [10] Mohapatra,R. N., Neutron-Anti-Neutron Oscillation: Theory and Phenomenology, J. Phys. G 36, 104006 (2009) ; Phillips, ,D. G. et al. ., Neutron-Antineutron Oscillations: Theoretical Status and Experimental Prospects, Phys. Rept. 612, 1 (2016)
  * [11] Jeans, J. H. MNRAS 85, 2 (1924).
  * [12] Goldman, I. and Nussinov, S. Limit on Continuous Neutrino Emission from Neutron Stars, JHEP 1008, 091 (2010)
  * [13] Mannarelli, M., Berezhiani, Z. , Biondi, R. and Tonnelli, F., Talk at the Nordita ESS workshop, December (2018).
  * [14] Cowsik, R. and Nussinov, S. some constraints on $\Delta B=2$ n anti-n oscillations Phys. Lett. 101B, 237 (1981).
  * [15] Weisberg, J. M., Nice, D. J., Taylor, J. H. . Timing Measurements of the Relativistic Binary Pulsar PSR B1913+16. The Astrophysical Journal 722, 1030 (2010)
  * [16] Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., and Hessels, J. W. T., A two-solar-mass neutron star measured using Shapiro delay, nature, 467, 1081 (2010)
  * [17] Ng, C., Kruckow, M. U., Tauris, T. M., et al., a young radio pulsar with a massive, compact companion, MNRAS, 476, 4315 (2018)
  * [18] Lazarus P., et al., Timing of a young mildly recycled pulsar with a massive white dwarf companion, MNRAS, 437, 1485 (2014).
  * [19] Bhat, N. D. R., Bailes, M. and Verbiest,J. P. W., Gravitational-radiation losses from the pulsar-white-dwarf binary PSR J1141-6545, Phys. Rev. D 77, 124017 (2008).
  * [20] Verbiest, J. P. W., et al. , Precision Timing of PSR J0437-4715: An Accurate Pulsar Distance, a High Pulsar Mass, and a Limit on the Variation of Newton’s Gravitational Constant, Astrophys. J. 679, 675 (2008)
  * [21] Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2010, The Equation of State from Observed Masses and Radii of Neutron Stars, astrophysical journal, 722, 33(2010)
  * [22] Berezhiani, Z. Neutron lifetime puzzle and neutron-mirror neutron oscillation. European Physical Journal C 72, 4849 (2019)


# Uniqueness of an association scheme related to the Witt design on 11 points

Alexander L. Gavrilyuk  Sho Suda

###### Abstract

It follows from Delsarte theory that the Witt $4$-$(11,5,1)$ design gives rise
to a $Q$-polynomial association scheme $\mathcal{W}$ defined on the set of its
blocks. In this note we show that $\mathcal{W}$ is unique, i.e., defined up to
isomorphism by its parameters.

##  1 Introduction

A $t$-$(v,k,\lambda)$ design is a set of $k$-subsets (called blocks) of $v$
points such that every $t$-subset is contained in exactly $\lambda$ blocks. A
$t$-design with $\lambda=1$ is called a Steiner system, and the most
celebrated ones have parameters $4$-$(11,5,1)$, $5$-$(12,6,1)$,
$3$-$(22,6,1)$, $4$-$(23,7,1)$, and $5$-$(24,8,1)$, which are often referred
to as the Mathieu designs or the Witt systems. In particular, the
$4$-$(11,5,1)$ design $W_{11}$ arises from a 4-transitive action of the
Mathieu group $M_{11}$ on 11 points, and its existence and uniqueness was
first shown by Witt [3]. The design $W_{11}$ has 66 blocks, and every two
distinct blocks $B$ and $B^{\prime}$ have $1$, $2$, or $3$ points in common.
Let us define a binary symmetric relation $R_{i}$ on $W_{11}$ by

$(B,B^{\prime})\in R_{i}\Leftrightarrow|B\cap B^{\prime}|=4-i,$

for $i=1,2,3$ and let $R_{0}=\\{(B,B)\mid B\in W_{11}\\}$. Then the pair
$\mathcal{W}=(W_{11},\\{R_{0},R_{1},R_{2},R_{3}\\})$ is an association scheme
(of 3 classes) by [7, Theorem 5.25]. Moreover, this scheme is $Q$-polynomial
(see [2] for the definitions and more results about $P,Q$-polynomial
association schemes). Williford [8, 9] compiled the tables of feasible
parameters of primitive 3-class $Q$-polynomial association schemes (on up to
2800 vertices), where the uniqueness of the scheme $\mathcal{W}$ was left
blank (see also [6, Appendix B]). In this note, we show that $\mathcal{W}$ is
unique, i.e., it is determined up to isomorphism by its parameters. The proof
is computer-assisted by Mathematica and relies on a spherical representation
of the scheme [1].

##  2 The parameters of $\mathcal{W}$

Let us recall some standard facts from the theory of association schemes (see
[2]). Let $A_{i}$ denote the logical matrix of the relation $R_{i}$, for
$i=0,1,2,3$. Then $A_{0}$ is the identity matrix of size $66$ and:

  1. (1)

$\sum_{i=0}^{3}A_{i}=J$, the square all-one matrix of size $66$,

  2. (2)

$A_{i}^{\top}=A_{i}$ ($0\leq i\leq 3$),

  3. (3)

$A_{i}A_{j}=\sum_{k=0}^{3}p_{i,j}^{k}A_{k}$, where $p_{i,j}^{k}$ are
nonnegative integers ($0\leq i,j,k\leq 3$), called the intersection numbers of
the scheme, which we refer to as the parameters of the scheme.

The intersection numbers of $\mathcal{W}$, written in the form of matrices
$\left(L_{i}\right)_{kj}=(p_{i,j}^{k})$, are found in the tables by Williford
[9]:

$\displaystyle L_{1}$ $\displaystyle=\left[\begin{array}[]{cccc}0&30&0&0\\\
1&15&10&4\\\ 0&15&6&9\\\ 0&8&12&10\\\ \end{array}\right],\quad
L_{2}=\left[\begin{array}[]{cccc}0&0&20&0\\\ 0&10&4&6\\\ 1&6&10&3\\\
0&12&4&4\\\ \end{array}\right],\quad
L_{3}=\left[\begin{array}[]{cccc}0&0&0&15\\\ 0&4&6&5\\\ 0&9&3&3\\\ 1&10&4&0\\\
\end{array}\right],$

and they determine ([2, Theorem 4.1]) the first and second eigenmatrices $P$
and $Q=\frac{1}{66}P^{-1}$ of $\mathcal{W}$:

$\displaystyle P$ $\displaystyle=\left[\begin{array}[]{cccc}1&30&20&15\\\
1&8&-2&-7\\\ 1&-1&-2&2\\\ 1&-6&8&-3\\\ \end{array}\right],\quad
Q=\left[\begin{array}[]{cccc}1&10&44&11\\\
1&\frac{8}{3}&-\frac{22}{15}&-\frac{11}{5}\\\
1&-1&-\frac{22}{5}&\frac{22}{5}\\\
1&-\frac{14}{3}&\frac{88}{15}&-\frac{11}{5}\\\ \end{array}\right],$

where the $(P)_{ij}$-entry ($0\leq i,j\leq 3$) is the eigenvalue of $A_{j}$ on
the $i$-th maximal common eigenspace of the matrices $A_{0},\ldots,A_{3}$
(which commute and hence can be simultaneously diagonalized).

######  Theorem 2.1.

An association scheme with the above parameters is isomorphic to
$\mathcal{W}$.

In what follows, we assume that
$\mathcal{X}=(X,\\{R_{0},R_{1},R_{2},R_{3}\\})$ is an association scheme with
the same parameters as $\mathcal{W}$. Let $E_{1}$ denote the orthogonal
projection matrix onto the 1st maximal common eigenspace (of multiplicity
$m_{1}=Q_{01}=10$) of the matrices $A_{i}$’s of $\mathcal{X}$. Since $E_{1}$
is positive semidefinite, we may regard $\frac{|X|}{m_{1}}E_{1}$ as the Gram
matrix of vectors in the unit sphere $S^{m_{1}-1}$ in $\mathbb{R}^{m_{1}}$,
and write $\frac{|X|}{m_{1}}E_{1}=F_{1}F_{1}^{\top}$ where $F_{1}$ is a
$|X|\times m_{1}$ matrix. We now identify a vertex $x\in X$ with the $x$-th
row $\mathbf{x}\in\mathbb{R}^{m_{1}}$ of $F_{1}$, and such a map is said to be
a spherical representation of the scheme $\mathcal{X}$ into its first
eigenspace. Define the angle set
$A(X):=\left\\{\langle\mathbf{x},\mathbf{y}\rangle\mid{x},{y}\in
X,{x}\neq{y}\right\\}$ and observe that
$A(X)=\left\\{\frac{Q_{j1}}{Q_{01}}\mid 1\leq j\leq 3\right\\}$. Note that the
map is injective, since $Q_{j1}\neq Q_{01}$ provided that $j\neq 0$; thus, it
defines all relations of the scheme. Therefore, to prove Theorem 2.1 it
suffices to show that the Gram matrix of the vertex set $X$ embedded into the
unit sphere $S^{9}$ in $\mathbb{R}^{10}$ is unique, up to orthogonal
transformation.

##  3 Proof of Theorem 2.1

Observe that the graph $(X,R_{2})$ is strongly regular with parameters
$(66,20,10,4)$. Such a strongly regular graph is isomorphic to the triangular
graph $T(12)$ by [5]; thus, we may identify the point set $X$ with the vertex
set of $T(12)$, which is ${[12]\choose 2}$, where $[12]:=\\{1,\ldots,12\\}$,
and $R_{2}$ with the edge set of $T(12)$, which is
$\left\\{\\{v,v^{\prime}\\}\mid|v\cap
v^{\prime}|=1,v,v^{\prime}\in{[12]\choose 2}\right\\}$. The next lemma can
easily be seen from this description of $T(12)$.

######  Lemma 3.1.

Each vertex of the triangular graph $T(12)$ is contained in two maximum
cliques of order $11$. Each vertex outside such a clique has exactly two
neighbors in it.

Note that a clique of order $11$ in $T(12)$ is a Delsarte clique, as it
attains the Delsarte bound [4, Proposition 1.3.2]. Fix a vertex $x=\\{1,2\\}$
of $T(12)$. The two Delsarte cliques containing $x$ are

$\displaystyle C_{1}=\\{\\{1,j\\}\mid 2\leq j\leq 12\\}\quad\text{ and }\quad
C_{2}=\\{\\{2,j\\}\mid j=1,3\leq j\leq 12\\}.$

Consider the image of $C_{1}$ in the spherical representation in
$\mathbb{R}^{10}$: since $C_{1}$ is a clique, the angle between any two
vectors is the same; thus, these 11 vectors form a regular simplex. Let
$V_{1}$ be an $11\times 10$ matrix whose row vectors correspond to the
vertices of $C_{1}$. Up to orthogonal transformation, we may assume that the
matrix $V_{1}$ has the following form:

$\displaystyle V_{1}=\left[\begin{smallmatrix}1&0&0&0&0&0&0&0&0&0\\\
-\frac{1}{10}&\frac{3\sqrt{11}}{10}&0&0&0&0&0&0&0&0\\\
-\frac{1}{10}&-\frac{\sqrt{11}}{30}&\frac{2\sqrt{55}}{15}&0&0&0&0&0&0&0\\\
-\frac{1}{10}&-\frac{\sqrt{11}}{30}&-\frac{\sqrt{55}}{60}&\frac{\sqrt{385}}{20}&0&0&0&0&0&0\\\
-\frac{1}{10}&-\frac{\sqrt{11}}{30}&-\frac{\sqrt{55}}{60}&-\frac{\sqrt{385}}{140}&\frac{\sqrt{1155}}{35}&0&0&0&0&0\\\
-\frac{1}{10}&-\frac{\sqrt{11}}{30}&-\frac{\sqrt{55}}{60}&-\frac{\sqrt{385}}{140}&-\frac{\sqrt{1155}}{210}&\frac{\sqrt{33}}{6}&0&0&0&0\\\
-\frac{1}{10}&-\frac{\sqrt{11}}{30}&-\frac{\sqrt{55}}{60}&-\frac{\sqrt{385}}{140}&-\frac{\sqrt{1155}}{210}&-\frac{\sqrt{33}}{30}&\frac{\sqrt{22}}{5}&0&0&0\\\
-\frac{1}{10}&-\frac{\sqrt{11}}{30}&-\frac{\sqrt{55}}{60}&-\frac{\sqrt{385}}{140}&-\frac{\sqrt{1155}}{210}&-\frac{\sqrt{33}}{30}&-\frac{\sqrt{22}}{20}&\frac{\sqrt{330}}{20}&0&0\\\
-\frac{1}{10}&-\frac{\sqrt{11}}{30}&-\frac{\sqrt{55}}{60}&-\frac{\sqrt{385}}{140}&-\frac{\sqrt{1155}}{210}&-\frac{\sqrt{33}}{30}&-\frac{\sqrt{22}}{20}&-\frac{\sqrt{330}}{60}&\frac{\sqrt{165}}{15}&0\\\
-\frac{1}{10}&-\frac{\sqrt{11}}{30}&-\frac{\sqrt{55}}{60}&-\frac{\sqrt{385}}{140}&-\frac{\sqrt{1155}}{210}&-\frac{\sqrt{33}}{30}&-\frac{\sqrt{22}}{20}&-\frac{\sqrt{330}}{60}&-\frac{\sqrt{165}}{30}&\frac{\sqrt{55}}{10}\\\
-\frac{1}{10}&-\frac{\sqrt{11}}{30}&-\frac{\sqrt{55}}{60}&-\frac{\sqrt{385}}{140}&-\frac{\sqrt{1155}}{210}&-\frac{\sqrt{33}}{30}&-\frac{\sqrt{22}}{20}&-\frac{\sqrt{330}}{60}&-\frac{\sqrt{165}}{30}&-\frac{\sqrt{55}}{10}\\\
\end{smallmatrix}\right],$

where the first row corresponds to $x$.

Now the problem is to add the remaining $55$ vectors such that together with
the vectors from $C_{1}$ their Gram matrix is permutation equivalent to
$E_{1}$. Since $A(X)=\\{\alpha_{1},\alpha_{2},\alpha_{3}\\}$, where
$\alpha_{1}=\frac{4}{15},\alpha_{2}=-\frac{1}{10},\alpha_{3}=-\frac{7}{15}$,
every such a vector $\mathbf{u}\in S^{9}$ in question satisfies
$\displaystyle{V_{1}\mathbf{u}^{\top}\in\\{\alpha_{1},\alpha_{2},\alpha_{3}\\}^{10}}$,
and there are at most $3^{10}$ candidates for $\mathbf{u}$.

We first determine the coordinates of the vectors in the image of $C_{2}$. Let
$z\in C_{1}\setminus\\{x\\}$ and $y\in C_{2}$. Then $(x,z)\in R_{2}$,
$(x,y)\in R_{2}$ and $(y,z)\in R_{1}\cup R_{2}\cup R_{3}$. Since
$p_{2,1}^{2}=6,p_{2,2}^{2}=10,p_{2,3}^{2}=3$ and, by Lemma 3.1, for each
vertex $z\in C_{1}\setminus\\{x\\}$, there is exactly one vertex $y\in
C_{2}\setminus\\{x\\}$ such that
$\langle\mathbf{y},\mathbf{z}\rangle=\alpha_{2}$, it follows that the vertices
in $C_{2}$ are taken from the following set:

$\displaystyle Y_{1}$ $\displaystyle=\left\\{\mathbf{y}\in S^{9}\mid
V_{1}\mathbf{y}^{\top}=(\alpha_{2},\mathbf{v})^{\top},\mathbf{v}=\left(\\{[\alpha_{1}]^{6},[\alpha_{2}]^{1},[\alpha_{3}]^{3}\\}\right)\right\\},$

where $\left(\\{[\alpha_{1}]^{6},[\alpha_{2}]^{1},[\alpha_{3}]^{3}\\}\right)$
denotes a vector of length 10 having $6$ entries equal to $\alpha_{1}$, $1$
entry $\alpha_{2}$, and $3$ entries $\alpha_{3}$. Note that $|Y_{1}|=840$.
Consider a graph with vertex set $Y_{1}$ and edge set $E_{1}$ defined by
$\\{\mathbf{y},\mathbf{y}^{\prime}\\}\in E_{1}$ if and only if
$\langle\mathbf{y},\mathbf{y}^{\prime}\rangle=\alpha_{2}$. Then every clique
of order $10$ in the graph $(Y_{1},E_{1})$ is a candidate for $C_{2}$. Let us
define a $10\times 10$ matrix $V_{2}=C\cdot V_{1}$, where

$C=\frac{1}{3}\setcounter{MaxMatrixCols}{11}\begin{bmatrix}[r]0&1&1&1&1&1&1&0&-1&-1&-1\\\
0&1&1&1&1&-1&-1&-1&1&1&0\\\ 0&1&1&0&-1&1&-1&1&1&-1&1\\\
0&1&1&-1&0&-1&1&1&-1&1&1\\\ 0&1&-1&1&-1&1&1&-1&0&1&1\\\
0&1&-1&-1&1&1&0&1&1&1&-1\\\ 0&0&-1&1&1&-1&1&1&1&-1&1\\\
0&-1&1&1&-1&0&1&1&1&1&-1\\\ 0&-1&1&-1&1&1&1&-1&1&0&1\\\
0&-1&0&1&1&1&-1&1&-1&1&1\\\ \end{bmatrix},$

and $V=\left[\begin{matrix}V_{1}\\\ V_{2}\end{matrix}\right]$. We can prove
the following lemma with the aid of a computer.

######  Lemma 3.2.

There exist exactly $30240$ cliques of order $10$ in the graph
$(Y_{1},E_{1})$. For every such clique $C_{2}$, the Gram matrix of $C_{1}\cup
C_{2}$ is permutationally equivalent to that of the rows of $V$.

By Lemma 3.2, we may fix the set $Y_{2}$ of the row vectors of $V_{2}$
representing $C_{2}$ and extend the above argument to other Delsarte cliques,
which yields that $X\setminus(C_{1}\cup C_{2})$ can be represented by a subset
of the following set:

$\displaystyle Y=\left\\{\mathbf{u}\in S^{9}\mid
V_{1}\mathbf{u}^{\top}=\mathbf{v}^{\top},\mathbf{v}=\left(\\{[\alpha_{1}]^{6},[\alpha_{2}]^{2},[\alpha_{3}]^{3}\\}\right)\right\\}\setminus
Y_{2}.$

Note that $|Y|=4610$, which is somewhat smaller than $3^{10}$. We determine
the set $Y$ and select only those of its vectors whose inner products with the
vectors from $Y_{2}$ belong to $A(X)$. Let $Z$ denote the set of such vectors.
With the aid of a computer, we obtain $|Z|=90$. We proceed by finding a
maximal subset of $Z$ such that the inner product of every two of its distinct
vectors is in $A(X)$.

######  Lemma 3.3.

The set $Z$ splits into $Z_{1}\cup Z_{2}$ in such a way that
$|Z_{1}|=|Z_{2}|=45$, $Z_{1}\cap Z_{2}=\emptyset$,
$A(Z_{1})=A(Z_{2})=\\{\alpha_{1},\alpha_{2},\alpha_{3}\\}$, and
$\langle\mathbf{z}_{1},\mathbf{z}_{2}\rangle\not\in\\{\alpha_{1},\alpha_{2},\alpha_{3}\\}$
for all $\mathbf{z}_{1}\in Z_{1}$ and all $\mathbf{z}_{2}\in Z_{2}$. Moreover,
the Gram matrices $C_{1}\cup C_{2}\cup Z_{1}$ and $C_{1}\cup C_{2}\cup Z_{2}$
are permutationally equivalent.

Thus, we may assume that the $66$ vertices of $X$ are represented by the row
vectors from $V_{1}$, $V_{2}$, and $Z_{1}$, and we can directly verify that,
together with the binary relations defined by their inner products, they form
an association scheme of $3$ classes with the same parameters as
$\mathcal{W}$. This completes the proof of Theorem 2.1.

### Acknowledgments

The authors would like to thank the reviewers for valuable comments. The
research of Alexander Gavrilyuk is supported by JSPS KAKENHI Grant Number
22K03403. The research of Sho Suda is supported by JSPS KAKENHI Grant Number
22K03410.

## References

  * [1] E. Bannai, E. Bannai, H. Bannai. Uniqueness of certain association schemes, European J. Combin., 25:261–267, 2004.
  * [2] E. Bannai and T. Ito. Algebraic Combinatorics I: Association Schemes. Benjamin/Cummings, Menlo Park, CA, 1984.
  * [3] T. Beth, D. Jungnickel. Mathieu groups, With designs, and Golay codes. Lecture Notes in Math. 893:157–169, 1981.
  * [4] A.E. Brouwer, A.E. Cohen, A. Neumaier. Distance-regular graphs. Springer-Verlag, Berlin, 1989. xviii+495.
  * [5] L.C. Chang. The uniqueness and nonuniqueness of the triangular association scheme, Sci. Record, 3:604–613, 1959.
  * [6] E.R. van Dam. Three-Class Association Schemes, J. of Alg. Combin., 10:69–107, 1999.
  * [7] P. Delsarte. An algebraic approach to the association schemes of coding theory. Philips Res. Rep. 10 (Suppl.), 1973.
  * [8] A.L. Gavrilyuk, J. Vidali, J.S. Williford. On few-class $Q$-polynomial association schemes: feasible parameters and nonexistence results, Ars Math. Contemp., 20(1):103–127, 2021.
  * [9] J.S. Williford. Tables of feasible parameter sets for primitive 3-class $Q$-polynomial association schemes.