\section{Introduction}
\label{Sec:Intro}
Constitutive relations play a major role in modeling response of mechanical systems.
These constitutive relations often can be written as a relation between kinematic
and dynamic variables. This relation can be explicit, meaning that
the value of a dynamic variable (such as force or stress) is given in terms of kinematic
variables (such as displacement or strain). In this case, the
governing differential equations of the mechanical system will be
in the form of Ordinary or Partial Differential Equations (ODEs or PDEs).
However, another possible case is the relation between kinematic and dynamic
variables to be implicit, hence, the value of dynamic variable might not be a function
(one-to-one) relation in terms of kinematic variables \citep{2003_Rajagopal_AoM}.
For example, one might have kinematic variables in terms of dynamic variables. An
example of these material models is the Bingham model \citep{1922_Bingham,1966_Peryzna_AAM}.
This case, renders the governing equations of the mechanical system as Differential-Algebraic
Equations (DAEs) or Partial Differential-Algebraic Equations (PDAEs) \citep{2010_Darbha_JFI}.
Obviously, these equations may not be solved exactly and one needs to resort to numerical
approximations.

Historically, numerical methods have been developed for ordinary (partial)
differential equations. Classical examples of time-integration algorithms
are designed to perform well when applied to ODEs. The efforts by
researchers in this case has been tremendously successful and forms the
basis for most simulations in scientific computing. However, numerical
integration of DAEs is still in a state of transition and development.
The reason why numerical integration of DAEs can be more
challenging than ODEs is that the presence of algebraic variables or
relations can invoke instabilities \citep{1995_Petzold_DAE,1991_Hairer_book}.
In fact, a numerical time-integrator can be successful in approximating
the solution to ODEs, but to fail in estimating the solution to a DAE.
Considering the importance of DAEs to mathematical modeling of mechanical
systems such as multi-body dynamics \citep{1997_Schiehlen_MSD,2010_BrulsCardona_JCND},
domain decomposition for PDEs \citep{2013_Bursi_CM,2014_Karimi_JCP,2015_Karimi_CMAME},
as well as implicit constitutive relations in mechanical systems \citep{2010_Rajagopal_MRC,2010_Darbha_JFI,2012_Prazak_AM}
and many others, many research endeavors have been centered around development of reliable
numerical integration of DAEs for various applications. It is worthy to note
that the algebraic constraints or relations in a DAE can be of different types.
The algebraic constraint can be smooth (differentiable) in which case
classical ODE solvers need minor adjustments to give accurate results. However,
as the case in this paper, the algebraic constraint (the implicit constitutive
relation) is continuous but non-smooth (not differentiable). Therefore,
numerical time-integration needs to be designed with more care.

In this paper, time-integration of a lumped system (see figure \ref{Fig:sdof})
governed by the implicit Bingham-Kelvin constitutive relation is studied. The
reason why we have chosen this rather simple system is to examine
the basic performance of the derived time-integration methods while avoiding
other numerical complexities. The system considered in this paper, has
a non-smooth response to excitations. Based on findings in \citep{2010_Darbha_JFI},
and generalized trapezoidal integration \citep{1990_Wood,2012_Hughes},
we derive an integration method that includes two independent parameters.
A complete description of numerical integration under the proposed methods for
cases of linear and nonlinear Bingham models is given. Finally,
the performance of the derived methods is compared to the results
from algorithm given in \citep{2010_Darbha_JFI} using several
numerical examples.

\subsection{Outline of the paper}
The rest of the paper is organized as follows. In section \ref{Sec:Method},
the governing equations and the new numerical scheme is detailed
and is followed by numerical examples in section \ref{Sec:Numerical}.
In section \ref{Sec:Conclusion} the conclusions are drawn and suggestions
for future research are given.

\section{Method for numerical Integration}
\label{Sec:Method}
\subsection{Model problem and governing equations}
Consider the SDOF system shown in figure \ref{Fig:sdof}, where
the block of mass $m$ is connected to a linear elastic spring
of stiffness $k$ and a dashpot with Bingham constitutive relation. The
governing equations for the motion of the block can be written as:
\begin{subequations}
\label{Eqn:SDOF_Equation}
	\begin{align}
		&\dot{v}(t) = \frac{1}{m}\left( f_{\mathrm{ext}}(t)
		- f_{s}(t) - f_{d}(t)\right) \\
		&f_{s}(t) = k u(t) \\
		&v(t) = \varphi(f_{d}(t)) \\
		&\dot{u}(t) = v(t).
	\end{align}
\end{subequations}
where $u(t)$ and $v(t)$ are the displacement and velocity of the
block. The force of dashpot and spring are shown by $f_{d}(t)$ and
$f_{s}(t)$, respectively. The function $\varphi$
is defined as follows:
\begin{align}
	\label{Eqn:dashpotFun}
	\varphi (f) := \left\{ \begin{array}{l l}
		0 & \big| f \big| \leq f_{y} \\
		\gamma \left( \big| f \big| - f_{y}\right)^{N}
		\mathrm{Sign}[f] & \big| f \big| \geq f_{y}
		\end{array}
	\right.,
\end{align}
where $\gamma$, $N$ and $f_{y}$ are constitutive parameters.
Equations \eqref{Eqn:SDOF_Equation} and \eqref{Eqn:dashpotFun}
form a DAE system of equations. In this problem, the velocity is
given as a function of the force of dashpot, hence, it is referred
to as an implicit constitutive relation. Figure \ref{Fig:BinghamNorton}
shows the relation between the dashpot force and the velocity for both
$N = 1$ and $N \neq 1$. From figure \ref{Fig:BinghamNorton} it can be
observed that the function $\varphi(f_{d})$ is a nondecreasing function.
Herein, inspired by \citep{2012_Prazak_AM}, we show that the solution to
the equation \eqref{Eqn:SDOF_Equation} is unique.

\emph{Proposition}. Consider initial conditions $u(t=0)=u_{0}$ and
$\dot{u}(t=0)=v_{0}$ and let $\varphi$ be a nondecreasing function,
then the solution to equation \eqref{Eqn:SDOF_Equation} is unique.

\emph{Proof}. Suppose $u$ and $\tilde{u}$ are both solutions to equation
\eqref{Eqn:SDOF_Equation}, with the initial conditions $u(t=0)=u_{0}$ and
$\dot{u}(t=0)=v_{0}$. Consider the \emph{arbitrary} time-interval $[0,\delta]$.
Since the initial conditions for both solutions $u$ and $\tilde{u}$ are
the same, then $u - \tilde{u} = 0$ for $t = 0$. From equation \eqref{Eqn:SDOF_Equation}
we can derive:
\begin{align}
	\label{Eqn:Proof1}
	\ddot{u} - \ddot{\tilde{u}} + \frac{1}{m}(f_{s} - \tilde{f}_{s})
	+ \frac{1}{m}(f_{d} - \tilde{f}_{d}) = 0,
\end{align}
where $f_{s} = k u$ and $\tilde{f}_{s} = k \tilde{u}$. Respectively,
$f_{d}$ and $\tilde{f}_{d}$ correspond to the dashpot force for
values of $\dot{u}$ and $\dot{\tilde{u}}$. Multiplying equation \eqref{Eqn:Proof1}
by $2(\dot{u} - \dot{\tilde{u}})$ gives
\begin{align}
	\label{Eqn:Proof2}
	2(\dot{u} - \dot{\tilde{u}})(\ddot{u} - \ddot{\tilde{u}}) +
	\frac{2}{m}(\dot{u} - \dot{\tilde{u}})(f_{s} - \tilde{f}_{s})
	+ \frac{2}{m}(\dot{u} - \dot{\tilde{u}})(f_{d} - \tilde{f}_{d}) = 0.
\end{align}
Note that the nondecreasing assumption on $\varphi$ of equation \eqref{Eqn:dashpotFun}
gives
\begin{align}
	\label{Eqn:Proof3}
	(\dot{u} - \dot{\tilde{u}})(f_{d} - \tilde{f}_{d}) \geq 0.
\end{align}
Denoting $w = u - \tilde{u}$ and using relations \eqref{Eqn:Proof2} and
\eqref{Eqn:Proof3} result in
\begin{align}
	\label{Eqn:Proof4}
	\frac{\mathrm{d} \dot{w}^{2}}{\mathrm{d}t} \leq -\frac{k}{2} \frac{\mathrm{d} w^{2}}{\mathrm{d}t}.
\end{align}
Note that based on definition of $w$ we have $w(t=0) = 0$ and $\dot{w}(t=0)=0$.
Integrating \eqref{Eqn:Proof4} over time interval $[0,\delta]$ gives
\begin{align}
	\label{Eqn:Proof5}
	\left(\dot{w}(t = \delta)\right)^{2} + \frac{k}{2} \left( w(t=\delta)\right)^2 \leq 0.
\end{align}
The only way relation \eqref{Eqn:Proof5} can be satisfied is to have
\begin{align}
	\label{Eqn:Proof6}
	\dot{w}(t = \delta) = 0 \quad \text{and} \quad w(t=\delta) = 0.
\end{align}
This result shows that the value of the solution to \eqref{Eqn:SDOF_Equation}
for a give initial condition is unique.$\Box$

Note that the proposition given above applies to both linear and nonlinear Bingham
models as in both cases $\varphi(f_{d})$ is a nondecreasing function.

Here, we will differentiate equation \eqref{Eqn:SDOF_Equation}b
with respect to time in order to reduce the index of this DAE system.
This reduction in DAE index results in a more stable numerical integration.
Doing so, results in the following system of equations:
\begin{subequations}
	\label{Eqn:SDOF_Equation2}
	\begin{align}
		&\dot{v}(t) = \frac{1}{m}\left( f_{\mathrm{ext}}(t)
		- f_{s}(t) - f_{d}(t)\right) \\
		&\dot{f}_{s}(t) = k v(t) \\
		&v(t) = \varphi(f_{d}(t)).
	\end{align}
\end{subequations}
Arguments on existence and uniqueness of a solution to this system
can be found in \citep{2010_Darbha_JFI,2012_Prazak_AM}, as well as
using the uniqueness agument given above. In the following section,
the efforts in reference \citep{2010_Darbha_JFI}
will be extended and applied to solve the system given in
\eqref{Eqn:SDOF_Equation2}.

\subsection{Numerical time integration}
Suppose that the time interval of interest is $(0,T]$ and is
divided into sub-intervals of equal size (time-step) $\Delta t$.
The value of variable $\Box$ at time level $n\Delta t$ will be
shown as $\Box^{(n)}$. We will begin by integrating the linear
momentum balance equation:
\begin{align}
	\label{Eqn:method1}
	\int_{t^{(n)}}^{t^{(n+1)}} \dot{v}(t) \mathrm{d}t =
	\frac{1}{m} \int_{t^{(n)}}^{t^{(n+1)}} \left( f_{\mathrm{ext}}(t)
	 - f_{s}(t) - f_{d}(t)\right) \mathrm{d}t.
\end{align}
Using trapezoidal integration, equation \eqref{Eqn:method1} results
in:
\begin{align}
	\label{Eqn:method2}
	v^{(n+1)} - v^{(n)} = \frac{\Delta t}{m} \left( (1 - \alpha)
	\left( f_{\mathrm{ext}}^{(n)} - f_{s}^{(n)} - f_{d}^{(n)}\right)
	+ \alpha \left( f_{\mathrm{ext}}^{(n+1)} - f_{s}^{(n+1)} -
	f_{d}^{(n+1)}\right)\right),
\end{align}
where $\alpha$ is the integration parameter ($0\leq\alpha\leq1$). The
equation \eqref{Eqn:SDOF_Equation2}b can be integrated in a similar
fashion to give:
\begin{align}
	\label{Eqn:method3}
	f_{s}^{(n+1)} - f_{s}^{(n)} = k\Delta t \left( (1-\beta)v^{(n)}
	+ \beta v^{(n+1)}\right),
\end{align}
where $\beta$ is an integration parameter in $[0,1]$. Equations
\eqref{Eqn:method2} and \eqref{Eqn:method3} can be combined to give:
\begin{align}
	\label{Eqn:method4}
	\left( 1 + \frac{\alpha \beta \Delta t^{2} k}{m}\right) v^{(n+1)}
	= \frac{\alpha \Delta t}{m}\left( \hat{f}^{(n+1)} - f_{d}^{(n+1)}\right),
\end{align}
where $\hat{f}^{(n+1)}$ is defined as:
\begin{align}
	\label{Eqn:method5}
	\hat{f}^{(n+1)} = f_{\mathrm{ext}}^{(n+1)} +
	\left( \frac{1}{\alpha} - 1\right)f_{\mathrm{ext}}^{(n)} -
	\frac{1}{\alpha} f_{s}^{(n)} - \left( \frac{1}{\alpha} - 1\right)f_{d}^{(n)}
	+ \left( \frac{m}{\alpha \Delta t} - k \Delta t (1 - \beta)\right)v^{(n)}.
\end{align}
The constitutive relation for dashpot given in equation \eqref{Eqn:dashpotFun}
shows us that $\mathrm{Sign}[v] = \mathrm{Sign}[f_{d}]$ at all time levels.
Also, from equation \eqref{Eqn:method4} we can derive that $\mathrm{Sign}[v^{(n+1)}]
= \mathrm{Sign}[\hat{f}^{(n+1)}]$, because
\begin{align}
	\label{Eqn:method6}
	\underbrace{\left(\left( 1 + \frac{\alpha \beta \Delta t^{2} k}{m}\right)
	\big|v^{(n+1)}\big| + \frac{\alpha \Delta t}{m} \big| f_{d}^{(n+1)}\big|
	\right)}_{\geq 0}\mathrm{Sign}[v^{(n+1)}]
	= \underbrace{\frac{\alpha \Delta t}{m} \big| \hat{f}^{(n+1)}}_{\geq 0}
	\big| \mathrm{Sign}[\hat{f}^{(n+1)}].
\end{align}
Equation \eqref{Eqn:method6} indicates that the value of $\hat{f}^{(n+1)}$
can be used as a predictor to determine the relation between $v^{(n+1)}$
and $f_{d}^{(n+1)}$. This equation implies that if $\big| f_{d}^{(n+1)} \big| \geq f_{y}$,
then we also have $\big| \hat{f}^{(n+1)} \big| \geq \big| f_{d}^{(n+1)}\big|$.
Therefore, if $\big| \hat{f}^{(n+1)} \big| \leq f_{y}$
then $v^{(n+1)} = 0$ and from equation \eqref{Eqn:method4} we get
$f_{d}^{(n+1)} = \hat{f}^{(n+1)}$. Otherwise, the force of dashpot needs
to be solved for using the following equation:
\begin{align}
	\label{Eqn:method7}
	\left( 1 + \frac{\alpha \beta \Delta t^{2} k}{m}\right)\gamma
	\left( \big| f_{d}^{(n+1)} \big| - f_{y}\right)^{N} \mathrm{Sign}[\hat{f}^{(n+1)}]
	- \frac{\alpha \Delta t}{m}\left( \hat{f}^{(n+1)} - f_{d}^{(n+1)}\right)=0.
\end{align}
Equation \eqref{Eqn:method7} is obtained from equations \eqref{Eqn:method4}
and \eqref{Eqn:dashpotFun}. Obviously, this equation is in general a
nonlinear equation. In the case of $N = 1$, equation
\eqref{Eqn:method7} can be simplified as:
\begin{align}
	\label{Eqn:method8}
	f_{d}^{(n+1)} = \frac{1}{1 + \frac{\alpha}{\gamma m \left( \frac{1}{\Delta t}
	+ \frac{\alpha \beta \Delta t k}{m}\right)}} \left( \frac{\alpha}{\gamma m
	\left( \frac{1}{\Delta t} + \frac{\alpha \beta \Delta t k}{m}\right)}
	\hat{f}^{(n+1)} + \mathrm{Sign}[\hat{f}^{(n+1)}] f_{y}\right).
\end{align}
Hence a nonlinear solver is not required. Once the value of the $f_{d}^{(n+1)}$
is found, the value of $v^{(n+1)}$ can be found to be
\begin{align}
	\label{Eqn:method9}
	v^{(n+1)} = \gamma \left( \big| f_{d}^{(n+1)} \big| - f_{y}\right)^{N}
	\mathrm{Sign}[f_{d}^{(n+1)}].
\end{align}
The value of the spring force can be updated using equation \eqref{Eqn:method3}
and the displacement of the mass can be found from the following relation:
\begin{align}
	\label{Eqn:method10}
	u^{(n+1)} = \frac{1}{k} f_{s}^{(n+1)}.
\end{align}
To complete the description of time-stepping algorithm, we need to clarify
two points:
\begin{enumerate}[(a)]
\item \textsf{Initialization}:~In this paper, we will assume all the
	initial displacements to be purely elastic. In other words, given an
	initial displacement $u_{0}$, the initial spring force will be taken
	to be $f_{s,0} = k u_{0}$ and initial dashpot force $f_{d,0}$ should be
	found from
	\begin{align}
		\label{Eqn:method11}
		\big| v_{0} \big| = \gamma \left( \big| f_{d,0}\big| - f_{y}\right)^{N}
	\end{align}
	if initial velocity $v_{0}$ is nonzero. Note that $\mathrm{Sign}[f_{d,0}] =
	\mathrm{Sign}[v_{0}]$. Consequently, equation \eqref{Eqn:method11} for the
	case of $N = 1$ gives:
	\begin{align}
		\label{Eqn:method12}
		f_{d,0} = \left( f_{y} + \frac{\big| v_{0} \big|}{\gamma} \right) \mathrm{Sign}[v_{0}].
	\end{align}

\item \textsf{Initial guess for nonlinear solver}:~Equation \eqref{Eqn:method7}
	is a nonlinear equation with respect to $f_{d}^{(n+1)}$. Hence, as with
	any nonlinear solver, an initial guess is needed. Based on our numerical
	experiments, the initial guess of $f_{y} \mathrm{Sign}[\hat{f}^{(n+1)}]$
	seems to be a reasonable choice.

\end{enumerate}
This concludes the description of the new time-stepping algorithm. A pseudocode
for this method is given in algorithm \ref{Alg:NumMethod}. In the
following section, numerical examples will be presented to demonstrate
the effect of various time-steps and time-integration parameters $\alpha$
and $\beta$.
\begin{algorithm}
	\caption{\textsf{Numerical time-integration algorithm}:~The pseudocode
	for time-stepping algorithm for elstic-viscoplastic SDOF system.}
	\label{Alg:NumMethod}
	\begin{algorithmic}[1]
		\STATE Initialize $f_{s,0}$, $f_{d,0}$ for given $u_{0}$ and $v_{0}$.
		\STATE Set values for $\Delta t$, $\alpha$ and $\beta$.
		\STATE Set $t = 0$.
		\WHILE{$t \leq T$}
			\STATE Set $t = t + \Delta t$.
			\STATE Calculate the predictor $\hat{f}$ from equation \eqref{Eqn:method5}.
			\IF{$\big| \hat{f} \big| \leq f_{y}$}
				\STATE Set the value of velocity $v = 0$.
				\STATE Set $f_{d} = \hat{f}$.
			\ELSE
				\STATE Find the value of dashpot force $f_{d}$ from equation \eqref{Eqn:method7}.
				\STATE Find the value of velocity from equation \eqref{Eqn:method9}.
			\ENDIF
			\STATE Find the value of the spring force $f_{s}$ from equation \eqref{Eqn:method3}.
			\STATE Update the value of displacement from equation \eqref{Eqn:method10}.
		\ENDWHILE
	\end{algorithmic}
\end{algorithm}

\section{Numerical Results}
\label{Sec:Numerical}
In this section, we will showcase the performance of the derived method in
comparison to the method proposed in \citep{2010_Darbha_JFI}. The amount of
damped (dissipated) energy will be denoted by $E_{d}(t)$ and is defined as
\begin{align}
	\label{Eqn:Edissip}
	E_{d}(t) := \int_{0}^{t} v(t) f_{d}(t) \;\mathrm{d}t.
\end{align}
Note that because $v(t)$ and $f_{d}(t)$ have the same sign, the value
of $E_{d}(t)$ will be nonnegative at all time-levels. The error in
kinematic quantities will be denoted by $e_{p}$ for $p=u,v$, and is
defined as
\begin{align}
	\label{Eqn:error}
	e_{p} := \frac{1}{M}\sqrt{\sum_{i=1}^{M}(p^{(i)} - p_{\mathrm{ref}}^{(i)})^{2}}\;,
\end{align}
where $M$ is the total number of time-steps and $p_{\mathrm{ref}}$ is
the reference (benchmark) values. The numerical results are obtained
from implementation in MATLAB \citep{MATLAB}.
\subsection{The case of $N = 1$}
In this problem, consider the SDOF system shown in figure \ref{Fig:sdof}
subject to external force given as:
\begin{align}
	\label{Eqn:fext}
	f_{\mathrm{ext}}(t) = 2 \sin (2 \pi t) e^{(-0.2 t)}.
\end{align}
The initial conditions are $u_{0} = 0$ and $v_{0} = 0$. The physical parameters
for this problem are given in Table \ref{Tbl:PhysBingN=1}.
To demonstrate the effect of time-integration parameters $\alpha$ and $\beta$,
we use different combinations of values of either 1 or 1/2 for each of them.
Time-integration parameters used in all cases are given in Table \ref{Tbl:BingN=1}.
Note that the benchmark in Table \ref{Tbl:BingN=1} refers to the method
proposed in \citep{2010_Darbha_JFI}, which coincides with $\alpha=1$ and $\beta=1$
under the method given in this paper. All the numerical solutions will be compared
to that of benchmark method.

The numerical results for displacement and velocity, against time
are shown in figures \ref{Fig:Bing_DispVel}. The dashpot
force $f_{d}$ and spring force $f_{s}$ is shown in figure \ref{Fig:Bing_force}.
It can be seen that the numerical results from the benchmark solution and case 1
(i.e., $\alpha=1,\;\beta=1$ and $\alpha=1,\;\beta=1/2$) are very different from the
cases 2 and 3. The numerical solution from case 1 is indeed very close to the solution
from the benchmark method. Repeating the same simulation with cases 2 and 3 of Table
\ref{Tbl:BingN=1} with smaller time-steps show that the damping properties in
values of dashpot force $f_{d}$ do not improve (see figure \ref{Fig:Bingham_case23}).
In fact the numerical values are very close to the ones with time-step of $10^{-4}$.
An important physical quantity is the energy damping in equation \eqref{Eqn:Edissip},
which shown in figure \ref{Fig:BingDissip} for benchmark case and case 1 of Table \ref{Tbl:BingN=1}.
The amount of energy dissipated calculated from each of these cases are very similar,
while cases 3 and 4 suggest zero energy damping. The solution from the benchmark
method and case 1 converge to one another in order $\mathcal{O}(\Delta t)$,
shown in figure \ref{Fig:Bing_error}. However, cases 2 and 3 underestimate
the dashpot force $f_{d}$ resulting in inaccurate values for displacement and
velocity. In fact, the value of $f_{d}$ in cases 2 and 3 seems to dampen excessively.

From this numerical experiment, we find that among the choices for parameters $\alpha$
and $\beta$, the values for case 1 and the benchmark method are the only viable ones.
It seems that choosing $\alpha = 1/2$, introduces excessive damping to the value of
$f_{d}$. Due to the non-smooth nature of the constitutive relation considered in this paper,
this property results in wrong values for other variables.

\begin{table}
	\centering
	\caption{\textsf{Parameters for case of $N=1$}:~Values for physical parameters
	are given.}
	\label{Tbl:PhysBingN=1}
	\begin{tabular}{ c c c c c c } \hline
	Parameter & $m$ & $k$ & $f_{y}$ & $\gamma$ & $N$ \\ \hline
	Value & 1 & 100 & 1 & 1 & 1 \\ \hline
	\end{tabular}
\end{table}
\begin{table}
	\centering
	\caption{\textsf{Parameters for case of $N=1$}:~The time-integration
	parameters for various cases are given.}
	\label{Tbl:BingN=1}
	\begin{tabular}{ c c c c }\hline
		Case & $\Delta t$ & $\alpha$ & $\beta$ \\ \hline
		benchmark & $10^{-6}$ & 1 &  1 \\
		1 & $10^{-4}$ & 1 & 1/2 \\
		2 & $10^{-4}$ & 1/2 & 1 \\
		3 & $10^{-4}$ & 1/2 & 1/2 \\ \hline
	\end{tabular}
\end{table}

\subsection{The case of $N \neq 1$}
In this case, the value of dashpot force $f_{d}$ is the solution of the nonlinear
equation \eqref{Eqn:method7}. For this case, the physical parameters are given in
Table \ref{Tbl:PhysNort}. Time-integration parameters are given in Table \eqref{Tbl:Nort}.
Here, we take the value of $N = 3$ and the benchmark numerical solution is from the
method in \citep{2010_Darbha_JFI}. The external force is the same as in equation
\eqref{Eqn:fext}, and the initial displacement and velocity are both zero.

The numerical results are given in figures \ref{Fig:Nort_DispVel}--\ref{Fig:Nort_force}.
Similar to the case of $N = 1$, the benchmark numerical solution and case 1 are
close and compatible with one another. However, the parameters selected for cases 2 and
3 are very different and predict a value of zero throughout the time-interval of interest
for velocity and displacement. In fact, in cases where $\alpha = 1/2$ the smoothing
results in damping in values of $f_{d}$. It is worth noting that unlike the case of $N=1$,
the displacement seems to converge to a point other than $u_{0}$ as time proceeds.
In the case of $N\neq1$, the amount of energy damped shown in figure \ref{Fig:NortDissip}
suggests a similar conclusion to the case of $N=1$. In this example, the energy dissipation
predicted by case 1 is very similar to the benchmark solution. Other cases of $\alpha$
and $\beta$ predict zero energy dissipation in this system.

This experiment concludes that the only viable value for $\alpha$ under the proposed
method is 1. The numerical solution to this nonlinear and non-smooth system is sensitive
to smoothing effects of trapezoidal integration in linear momentum balance equation.
The numerical results seem to be much less sensitive to the choice of $\beta$ between
1/2 and 1.

\begin{table}
	\centering
	\caption{\textsf{Parameters for case of $N\neq1$}:~Physical
	parameters in the case of $N \neq 1$ are given.}
	\label{Tbl:PhysNort}
	\begin{tabular}{ c c c c c c } \hline
	Parameter & $m$ & $k$ & $f_{y}$ & $\gamma$ & $N$ \\ \hline
	Value & 1 & 10 & 1 & 1 & 3 \\ \hline
	\end{tabular}
\end{table}
\begin{table}
	\centering
	\caption{\textsf{Parameters for case of $N\neq1$}:~Values of the
	time-integration parameters.}
	\label{Tbl:Nort}
	\begin{tabular}{ c c c c }\hline
		Case & $\Delta t$ & $\alpha$ & $\beta$ \\ \hline
		benchmark & $10^{-7}$ & 1 &  1 \\
		1 & $10^{-7}$ & 1 & 1/2 \\
		2 & $10^{-7}$ & 1/2 & 1 \\
		3 & $10^{-7}$ & 1/2 & 1/2 \\ \hline
	\end{tabular}
\end{table}

\subsection{Implicit-explicit time-integration}
Thus far, the numerical experiments suggest that the only viable value for
$\alpha$ is 1. However, for $\alpha=1$, the numerical results for various
values of $\beta$ seem to be very close. In this section, we showcase the
numerical results for the case of $\alpha=1$ and $\beta = 0$. For $\beta=0$,
the equation \eqref{Eqn:method3} will become an explicit Euler integration
scheme while equation \eqref{Eqn:method2} remains implicit Euler integration.
In this case, the constitutive relation of the elastic spring is integrated
using an explicit scheme and the momentum balance is integrated via an
implicit relation.

The problem is solved for the external force in equation \eqref{Eqn:fext} using
time-steps $\Delta t = 10^{-4}, 10^{-5}$ and $10^{-6}$. The error (compared to
the benchmark numerical result) is shown in figure \ref{Fig:ImExError}. It can be
seen that the numerical solutions from the two methods converge as $\mathcal{O}(\Delta t)$,
which implies that the choice of $\alpha=1$ and $\beta=0$ is a viable choice of
time-integration parameters. This method, therefore, does not suffer from excessive
numerical damping for values of dashpot force $f_{d}$.

\section{Conclusion}
\label{Sec:Conclusion}
In this paper the governing differential-algebraic equations for a
SDOF system obeying Bingham constitutive model were reviewed and
a proof for uniqueness of the solution was presented. Using the
generalized trapezoidal integration a method for numerical integration
of the given equations was derived with two independent time-integration
parameters. A pseudocode for this numerical
algorithm along with initialization and initial guess for the
nonlinear equations was provided. This derivation was followed by
numerical examples showcasing its performance compared to a
benchmark method. The findings from these numerical experiments
can be summarized as below:
\begin{enumerate}[(1)]
	\item Because of the non-smooth nature of the physical system
	considered in this paper, accuracy in estimation of dashpot force
	is of utmost importance. Choice of time-integration parameters
	can lead to excessive numerical damping in value of dashpot force
	, which in turn gives completely different numerical values for
	other kinematic and dynamic variables.
	\item Under the presented numerical scheme in this paper, the
	only viable choice of integration parameter $\alpha$ is 1.
	Other values of $\alpha$ result in additional (unphysical)
	smoothing in the numerical solution.
	\item Using the method of numerical time integration in this paper
	the value of integration parameter $\beta$ can be anything in range
	$[0,1]$. It seems that the values of $\beta$ do not introduce any
	more significant damping to the estimated value of dashpot force.
	\item Given $\alpha=1$ and $\beta \in [0,1]$ under the numerical
	scheme here, the numerical solution converges to the solution
	from the benchmark method in order $\mathcal{O}(\Delta t)$.
	\item Under the proposed method in this paper, one can use the
	$\alpha=1$ and $\beta=0$ values for integration parameters. This
	means that one can have a convergent implicit-explicit time-integration
	of the non-smooth system considered here.
\end{enumerate}
Future developments in this area can be extension of such methods to
other implicit and non-smooth material models Maxwell type
\citep{ChristensenBook}, extension to 2 and 3-dimensional problems,
and extension of Newmark-type integration methods \citep{1959_Newmark,1993_Hulbert_JAM} for this class of problems.


\section{Introduction}

Let $n$ be a positive integer,
and let $\mathbb{F}$ be any field.
The \emph{subspace lattice} $\mathcal{P}_n(\mathbb{F})$
(also known as the \emph{finite projective space})
is the poset of linear subspaces of $\mathbb{F}^n$.
The general linear group $\mathrm{GL}_n(\mathbb{F})$
acts on $\mathcal{P}_n(\mathbb{F})$.
The natural grading structure
of $\mathcal{P}_n(\mathbb{F})$
is given by $\mathrm{GL}_n(\mathbb{F})$-action,
and
each fiber (i.e., orbit) is called a \emph{Grassmannian}.
In other words,
the Grassmannian $\mathrm{Gr}(m,n)$ is the set of
$m$-dimensional subspaces in $\mathbb{F}^n$, where $0 \le m \le n$.
Let $\mathcal{B}$ denote the Borel subgroup of $\mathrm{GL}_n(\mathbb{F})$.
Then, the $\mathcal{B}$-action defines
a finer ``hyper-cubic'' grading structure
of $\mathcal{P}_n(\mathbb{F})$,
and
each fiber contained in $\mathrm{Gr}(m,n)$
is called a \emph{Schubert cell of $\mathrm{Gr}(m,n)$}.
See \cite{MR2474907} for details.
The author showed in \cite{W}
that
the algebra defined from the ``hyper-cubic'' grading structure of
$\mathcal{P}_n(\mathbb{F})$ together with its incidence structure
has a close relation to the quantum affine algebra
$U_q(\widehat{\mathfrak{sl}}_2)$, if $\mathbb{F}$ is a finite field of $q^2$ elements.
In this paper,
we study the Schubert cells of a Grassmannian
from the combinatorial point of view of \emph{association schemes}.
More precisely,
we show that the association scheme
defined by the $\mathcal{B}$-action
on each Schubert cell is a
\emph{generalized wreath product} of
one-class association schemes
with the base set $\mathbb{F}$.
The concept of
a generalized wreath product
of association schemes
was introduced by
R.~A.~Bailey~\cite{MR2206477} in 2006.
The (usual) wreath product
of association schemes has
been actively studied (see e.g.,
\cite{MR2712017, MR2610292, MR2002219, MR2964723, MR3612420, MR2747797, MR3047011, MR2802177}),
and we may view the result of this paper as demonstrating
the fundamental importance of Bailey's generalization as well.

Before the main discussion,
we briefly recall the notion of the generalized wreath product of association schemes.
For the definition of association schemes,
see \cite{MR882540, MR2535398, MR2184345},
and
for the theory of posets,
see \cite{MR2868112}.
Let $(X,\le)$ be a nonempty finite poset.
A subset $Y$ in $X$ is called an \emph{anti-chain}
if any two elements in $Y$ is incomparable.
For an anti-chain $Y$ in $X$,
define the \emph{down-set} (also known as the \emph{order ideal}) by
\[
\mathrm{Down}(Y) = \lbrace x \in X \mid \text{$x < y$ for some $y \in Y$}\rbrace.
\]
Note that this definition follows \cite{MR2206477}
and it is different from \cite{MR2868112},
where
$\mathrm{Down}(Y) \cup Y$ is called
the down-set of $Y$.
For each $x \in X$,
let $\mathcal{Q}_x$ denote an $r_x$-class association scheme
on a set $\Omega_x$.
We do not assume either $\mathcal{Q}_x$ is symmetric
or $\Omega_x$ is finite.
Let $R_{x,i}$ denote the $i$-th associate class
for $i \in \lbrace 0,1,\ldots,r_x\rbrace$.
By convention,
we choose the index
so that $R_{x,0} = \lbrace (\omega, \omega) \mid \omega \in \Omega_x\rbrace$.
We set $\Omega = \prod_{x \in X} \Omega_x$.
For each anti-chain $Y$ in $X$
and for each $(i_x)_{x \in Y} \in \prod_{x \in Y} \lbrace 1,2,\ldots,r_x\rbrace$,
let
$R(Y,(i_x)_{x \in Y})$ denote the set of
$((\alpha_x)_{x \in X}, (\beta_x)_{x \in X}) \in \Omega \times \Omega$
satisfying
(i) $\alpha_x = \beta_x$ if $x \in X \setminus (Y \cup \mathrm{Down}(Y))$,
and
(ii) $(\alpha_x, \beta_x) \in R_{x,i_x}$ if $x \in Y$.
Let $\mathcal{R}$ denote the set of
$R(Y,(i_x)_{x \in Y})$ for all anti-chains $Y$ in $X$ and
$(i_x)_{x \in Y} \in \prod_{x \in Y} \lbrace 1,2,\ldots,r_x\rbrace$.

\begin{thm}[cf. {\cite[Theorem 3]{MR2206477}}]\label{Bailey}
The pair $(\Omega, \mathcal{R})$ is an association scheme.
\end{thm}

The association scheme $(\Omega, \mathcal{R})$ in Theorem \ref{Bailey} is called the
\emph{generalized wreath product of $\mathcal{Q}_x$
over the poset $X$}.
We remark that Bailey~\cite[Theorem 3]{MR2206477}
assumes that each base set $\Omega_x$ is finite
and each association scheme $\mathcal{Q}_x$ is symmetric.
However, the theorem is still true if
we drop both of these assumptions.

\section{Subspace lattices}

Throughout this paper, we fix a positive integer $n$ and a field $\mathbb{F}$.
Let $\mathbb{F}^n$ denote the $n$-dimensional column vector space over $\mathbb{F}$.
By the \emph{subspace lattice}, denoted by
$\mathcal{P}_n(\mathbb{F})$,
we mean the poset consisting of all subspaces
in $\mathbb{F}^n$ with partial order given by inclusion.
We fix
the sequence $\lbrace V_i \rbrace_{i = 0}^n$ in $\mathcal{P}_n(\mathbb{F})$ such that
each $V_i$ consists of vectors whose
bottom $n-i$ entries are zero.
We remark that
$\lbrace V_i \rbrace_{i = 0}^n$ is a \emph{(complete) flag} (i.e., a maximal chain) in $\mathcal{P}_n(\mathbb{F})$.

Let $\mathrm{Mat}_n(\mathbb{F})$
denote the set of $n \times n$ matrices with entries in $\mathbb{F}$.
Let $\mathrm{GL}_n(\mathbb{F})$ denote the set of
invertible matrices in $\mathrm{Mat}_n(\mathbb{F})$.
Observe that
$\mathrm{GL}_n(\mathbb{F})$ acts on $\mathbb{F}^n$
by left multiplication,
and hence it acts on $\mathcal{P}_n(\mathbb{F})$.
Let $\mathcal{B}$ denote the (Borel) subgroup in $\mathrm{GL}_n(\mathbb{F})$ stabilizing
$\lbrace V_i \rbrace_{i = 0}^n$.
In other words,
$\mathcal{B}$ consists of all
upper triangular invertible matrices in $\mathrm{Mat}_n(\mathbb{F})$.
In this paper,
we consider the $\mathcal{B}$-action on $\mathcal{P}_n(\mathbb{F})$.

A matrix in $\mathrm{Mat}_n(\mathbb{F})$
is said to be in \emph{reverse column echelon form}
if the following two conditions are met:
\begin{enumerate}
\item[(CE1)] Any zero columns are right of all nonzero columns.
\item[(CE2)] The last nonzero entry of a nonzero column is always
strictly below of the last nonzero entry of its right column.
\end{enumerate}
A matrix in $\mathrm{Mat}_n(\mathbb{F})$
is said to be in \emph{reduced reverse column echelon form}
if it is in reverse column echelon form and
the following third condition is also met:
\begin{enumerate}
\item[(CE3)] Every last nonzero entry of a nonzero column is $1$ and
is the only nonzero entry in its row.
\end{enumerate}

By elementary linear algebra, we have the following.
\begin{prop}\label{prop:bij}
There exists a bijection between the following two sets:
\begin{enumerate}
\item the set of all matrices in $\mathrm{Mat}_n(\mathbb{F})$
in reduced reverse column echelon form,
\item the set $\mathcal{P}_n(\mathbb{F})$ of all subspaces in $\mathbb{F}^n$,
\end{enumerate}
that sends a matrix in $\mathrm{Mat}_n(\mathbb{F})$
to its column space.
\end{prop}
\begin{proof}
See for instance \cite{MR3013937}.
\end{proof}

\begin{ex}[$n=7$]\label{ex1}
Let $e_1, e_2, \ldots, e_7$ denote the standard basis for $\mathbb{F}^7$.
Suppose $U$ is the $4$-dimensional subspace
in $\mathbb{F}^7$ given by
$U = \mathrm{Span}\lbrace
8e_1+6e_3+4e_6+2e_7,
8e_1+9e_3+e_4+e_5,
4e_1+e_2+5e_3+e_4,
3e_1+e_2\rbrace$.
Then the following is the matrix corresponding to $U$
by the bijection in Proposition \ref{prop:bij}:
\[
M =
\begin{pmatrix}
4 & 7 & 1 & 3 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 & 0 \\
3 & 4 & 5 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 & 0 & 0 & 0 \\
2 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 0
\end{pmatrix}.
\]
\end{ex}
Observe that
for
$M, N \in \mathrm{Mat}_n(\mathbb{F})$
in reduced reverse column echelon form
and for $G \in \mathcal{B}$,
the column space of $N$ moves to
that of $M$ by the $G$-action
if and only if
$M$ and $GN$ are column equivalent.
Since $GN$ is in reverse column echelon form
(but not necessarily reduced),
these conditions are equivalent to
$M = GNH^{T}$ for some $H \in \mathcal{B}$.
For notational convenience,
we write $M \sim N$
if there exist $G, H \in \mathcal{B}$ such that $M = GNH^T$.
Observe that $\sim$ is an equivalence relation on $\mathrm{Mat}_n(\mathbb{F})$.

For the rest of this paper,
we will identify $\mathcal{P}_n(\mathbb{F})$ with the set of
all matrices in $\mathrm{Mat}_n(\mathbb{F})$ in reduced reverse column echelon form
by the bijection in Proposition \ref{prop:bij}.

\section{The $\mathcal{B}$-action on $\mathcal{P}_n(\mathbb{F})$}

For a positive integer $m$, we write
$[m] = \lbrace 1,2,\ldots, m\rbrace$.
We define a partial order in
the index set $[n] \times [n]$ of matrices in $\mathrm{Mat}_n(\mathbb{F})$ by
$(i,j) \le (k,l)$ if $i \le k$ and $j \le l$.
This is known as the \emph{direct product order}
in \cite[Section~3.2]{MR2868112}.
For $M \in \mathrm{Mat}_n(\mathbb{F})$,
by the \emph{support} of $M$, denoted by $\mathrm{Supp}(M)$, we mean
the subposet of $[n] \times [n]$
consisting of all indices $(i,j) \in [n] \times [n]$ with $M_{i,j} \neq 0$.
The \emph{pivot-set} of $M$,
denoted by $\mathrm{Piv}(M)$,
is the set of all maximal elements in $\mathrm{Supp}(M)$.
Each element in the pivot-set is called a \emph{pivot}.
Observe that
$(i,j) \in \mathrm{Piv}(M)$
if and only if
$M_{i,j} \neq 0$ and $M_{k,l} = 0$ if $(k,l) > (i,j)$.
We remark that
every entry indexed by a pivot of a matrix in $\mathcal{P}_n(\mathbb{F})$
must be $1$ by the condition (CE3).

\begin{lem}\label{lem:wequiv}
For $M, N \in \mathcal{P}_n(\mathbb{F})$,
the following are equivalent:
\begin{enumerate}
\item $M \sim N$,
\item $\mathrm{Piv}(M) = \mathrm{Piv}(N)$.
\end{enumerate}
\end{lem}
\begin{proof}
(i) $\Rightarrow$ (ii)
Suppose $M \sim N$.
There exist
$G, H \in \mathcal{B}$ such that $M = GNH^T$.
It suffices to show $\mathrm{Piv}(N) \subseteq \mathrm{Piv}(M)$.
For $(i,j) \in \mathrm{Piv}(N)$,
we have
$N_{i,j} = 1$ and $N_{k,l} = 0$ if $(k,l) > (i,j)$.
Since
$G, H$ are upper triangular,
we have
\[
M_{k,l} = \sum_{s = k}^n\sum_{t = l}^n G_{k,s}N_{s,t}H_{l,t}
=
\begin{cases}
G_{i,i}H_{j,j} & \text{if $(k,l) = (i,j)$},\\
0 & \text{if $(k,l) > (i,j)$}.
\end{cases}
\]
Since $G, H$ are invertible,
$G_{i,i}H_{j,j} \neq 0$.
These imply $(i,j) \in \mathrm{Piv}(M)$
and hence $\mathrm{Piv}(N) \subseteq \mathrm{Piv}(M)$.

(ii) $\Rightarrow$ (i)
Suppose $\mathrm{Piv}(M) = \mathrm{Piv}(N)$.
Take $X \in \mathcal{P}_n(\mathbb{F})$ with
$X_{i,j} = 1$ if $(i,j) \in \mathrm{Piv}(M)$ and $X_{i,j} = 0$ otherwise.
Observe that
for each $j \in [n]$,
there exists at most one $k$ such that $(j,k) \in \mathrm{Piv}(M)$
and then we
define $G \in \mathrm{Mat}_n(\mathbb{F})$ by
\[
G_{i,j} =
\begin{cases}
M_{i,k} & \text{if $(j,k) \in \mathrm{Piv}(M)$ for some $k$},\\
\delta_{i,j} & \text{if there is no $k$ such that  $(j,k) \in \mathrm{Piv}(M)$}
\end{cases}
\]
for $i,j \in [n]$.
Then we have $G_{i,i} = 1$ for $i \in [n]$
and $G_{i,j} = 0$ if $j < i$ for $i, j \in [n]$.
Thus
$G \in \mathcal{B}$.
By the direct calculation, we have $M = GX$
and
hence, $M \sim X$.
Similarly we have $N \sim X$ and so $M \sim N$.
\end{proof}

Let $1 \le m \le n-1$ and $M \in \mathcal{P}_n(\mathbb{F})$
with $\rank M = m$.
Note that
we avoid the trivial cases $m = 0$ and $m = n$.
Since the pivots of $M$ lie in the first $m$ columns,
$\mathrm{Piv}(M)$
is an anti-chain in $[n] \times [m]$ of size $m$.
For $1 \le m \le n-1$
and for an anti-chain $\alpha$ in $[n] \times [m]$
 of size $m$,
we set
\begin{equation}
\mathcal{O}_{\alpha} = \lbrace M \in \mathcal{P}_n(\mathbb{F}) \mid
\mathrm{Piv}(M) = \alpha
\rbrace.
\label{OA}
\end{equation}
For each $1 \le m \le n-1$ and
each anti-chain $\alpha$ in $[n] \times [m]$ of size $m$,
consider
$M \in \mathcal{P}_n(\mathbb{F})$ with
$M_{i,j} = 1$ if $(i,j) \in \alpha$ and $M_{i,j} = 0$ otherwise.
Then we have $M \in \mathcal{O}_\alpha$
and in particular $\mathcal{O}_\alpha \neq \emptyset$.

\begin{prop}\label{prop:OA}
The rank of any matrix in \eqref{OA}
is $m = |\alpha|$.
Moreover,
each subset \eqref{OA}
is an orbit of the $\mathcal{B}$-action on $\mathcal{P}_n(\mathbb{F})$.
\end{prop}
\begin{proof}
Immediate from the construction and Lemma \ref{lem:wequiv}.
\end{proof}

Recall the Grassmannian $\mathrm{Gr}(m,n)$
and we identify $\mathrm{Gr}(m,n)$ with a set of matrices
by the bijection in Proposition \ref{prop:bij}.
In other words,
\[
\mathrm{Gr}(m,n) = \lbrace M \in \mathcal{P}_n(\mathbb{F}) \mid
\rank M = m
\rbrace.
\]
By Proposition \ref{prop:OA},
each $\mathcal{O}_\alpha$ in \eqref{OA} is
a $\mathcal{B}$-orbit in $\mathrm{Gr}(m,n)$,
where $m = |\alpha|$.
Thus,
it is called a \emph{Schubert cell of a Grassmannian} \cite{MR2474907}.

\begin{ex}[$n=7$, $m = 4$]\label{ex2}
Take $M \in \mathcal{P}_7(\mathbb{F})$ as in Example \ref{ex1}.
Then we have
$\mathrm{Piv}(M) = \lbrace
(2,4), (4,3), (5,2), (7,1)
\rbrace$.
Moreover,
$\mathcal{O}_{\mathrm{Piv}(M)}$ is
the set of matrices of the form
\begin{equation}
\begin{pmatrix}
\ast & \ast & \ast & \ast & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 & 0 \\
\ast & \ast & \ast & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 & 0 & 0 & 0 \\
\ast & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 0
\end{pmatrix},
\label{ex:Oalpha}
\end{equation}
where the symbol $\ast$ denotes an arbitrary element
in $\mathbb{F}$.
\end{ex}

\begin{lem}\label{lem:diag}
Let $1 \le m \le n-1$
and let $\alpha$ denote an anti-chain in $[n] \times [m]$ of size $m$.
For $M, N, M', N' \in \mathcal{O}_\alpha$,
the following are equivalent:
\begin{enumerate}
\item $(M,N)$ moves to $(M',N')$ by the diagonal $\mathcal{B}$-action,
\item $\mathrm{Piv}(M-N) = \mathrm{Piv}(M'-N')$.
\end{enumerate}
\end{lem}
\begin{proof}
(i) $\Rightarrow$ (ii)
Suppose there exist $G, H, K \in \mathcal{B}$
such that $M' = GMH^T$ and $N' = GNK^T$.
Then we have
\begin{equation}
\mathrm{Piv}(GM) = \mathrm{Piv}(GN) = \alpha,
\label{*}
\end{equation}
\begin{equation}
\mathrm{Piv}(GM-GN) = \mathrm{Piv}(M-N)
\label{**}
\end{equation}
since $G \in \mathcal{B}$ (cf. Lemma \ref{lem:wequiv}).
We write $\alpha = \lbrace(k_r,r) \mid r \in [m] \rbrace$
and observe that $k_1 > k_2 > \cdots > k_m$
and that $\mathrm{Piv}(M-N) \subseteq \mathrm{Down}(\alpha)$.

Take $(i,j) \in \mathrm{Piv}(M-N)$.
Observe that $j \in [m]$ and $k_j > i$
and hence
there exists $m' = \max\lbrace r \in [m] \mid k_r > i\rbrace$.
For $r, l \in [m]$ with $j \le r \le m'$ and $j \le l \le m'$,
since $H, K$ are upper triangular,
we have
\begin{align*}
M'_{k_r,l} &= \sum_{t = l}^n (GM)_{k_r,t}H_{l,t}, \\
N'_{k_r,l} &= \sum_{t = l}^n (GN)_{k_r,t}K_{l,t}.
\end{align*}
In the above equations, we have the following:
For each $r,l$, we have $M'_{k_r,l} = N'_{k_r,l} = \delta_{r,l}$ by (CE3);
For each $r,t$, we have $(GM)_{k_r,t} = (GN)_{k_r,t}$
since $(k_r,t) > (i,j)$ and by \eqref{**};
For each $r, t$ with $r < t$, we have $(GM)_{k_r,r} = (GN)_{k_r,r} = G_{k_r,k_r} \neq 0$
and $(GM)_{k_r,t} = (GN)_{k_r,t} = 0$
since $(k_r,r) \in \alpha$ and by \eqref{*}.
By these comments,
for each $l \in [m]$ with $j \le l \le m'$,
both $(H_{l,l}, H_{l,l+1}, \ldots, H_{l,m'})$
and $(K_{l,l}, K_{l,l+1}, \ldots, K_{l,m'})$
are solutions to the same system of $m'-j+1$ independent linear equations.
Hence, $H_{l,t} = K_{l,t}$ for $j,t \in [m]$ with
$j \le l \le t \le m'$.

For $(k,l) \in [n] \times [n]$ with $(k,l) > (i,j)$,
we have
$(GM)_{k,t} = (GN)_{k,t}$ if $t \ge l$ since $(k,t) > (i,j)$ and by \eqref{**},
and we also have
$(GM)_{k,t} = (GN)_{k,t} = 0$ if $t > m'$ by the definition of $m'$ and by \eqref{*}.
Recall that we have shown $H_{l,t} = K_{l,t}$ if $j \le l \le t \le m'$.
By these comments, we have
\begin{align*}
M'_{k,l} = \sum_{t = l}^{n} (GM)_{k,t}H_{l,t}
= \sum_{t = l}^{n} (GN)_{k,t}K_{l,t} = N'_{k,l}.
\end{align*}
Similarly,
we have
\begin{align*}
M'_{i,j} - N'_{i,j} = \sum_{t = j}^{n} (GM)_{i,t}H_{j,t}
- \sum_{t = j}^{n} (GN)_{i,t}K_{j,t}
= ((GM)_{i,j} - (GN)_{i,j})H_{j,j}.
\end{align*}
In the above equations,
we have $(GM)_{i,j} \neq (GN)_{i,j}$ by \eqref{**},
and we also have $H_{j,j} \neq 0$ since $H \in \mathcal{B}$.
Therefore we obtain $M'_{i,j} \neq N'_{i,j}$.
These imply $(i,j) \in \mathrm{Piv}(M'-N')$
and hence $\mathrm{Piv}(M-N) \subseteq \mathrm{Piv}(M'-N')$.
Since $G, H, K$ are invertible,
we also have
$\mathrm{Piv}(M'-N') \subseteq \mathrm{Piv}(M-N)$.
Consequently, we have
$\mathrm{Piv}(M-N) = \mathrm{Piv}(M'-N')$.

(ii) $\Rightarrow$ (i)
Let
$M, N, M', N' \in \mathcal{O}_\alpha$
with $\mathrm{Piv}(M-N) = \mathrm{Piv}(M'-N')$.
Take $X \in \mathcal{O}_\alpha$ with
$X_{i,j} = 1$ if $(i,j) \in \alpha \cup \mathrm{Piv}(M-N)$ and $X_{i,j} = 0$ otherwise,
and $Y \in \mathcal{O}_\alpha$ with
$Y_{i,j} = 1$ if $(i,j) \in \alpha$ and $Y_{i,j} = 0$ otherwise.
Define $G \in \mathrm{Mat}_n(\mathbb{F})$ by
\[
G_{i,j} = \begin{cases}
N_{i,k} & \text{if $(j,k) \in \alpha$ for some $k$},\\
M_{i,k} - N_{i,k} & \text{if $(j,k) \in \mathrm{Piv}(M-N)$ for some $k$}, \\
\delta_{i,j} & \text{if there is no $k$ such that $(j,k) \in \alpha \cup \mathrm{Piv}(M-N)$}
\end{cases}
\]
for $i,j \in [n]$.
Then we have $G \in \mathcal{B}$
and $M = GX$ and $N = GY$.
Similarly there exists $G' \in \mathcal{B}$
such that
$M' = G'X$ and $N' = G'Y$.
Therefore
$(M,N)$ moves to $(M',N')$ by the diagonal action of
$G'G^{-1}$.
\end{proof}

Let $1 \le m \le n-1$
and let $\alpha$ denote an anti-chain in $[n] \times [m]$ of size $m$.
Let $M, N \in \mathcal{O}_\alpha$.
Observe that
$\mathrm{Piv}(M-N)$
is an anti-chain in
\begin{equation}
\mathcal{D}(\alpha)
=
\lbrace (i,j) \in \mathrm{Down}(\alpha) \mid
\text{there is no $k$ such that $(i,k) \in \alpha$}\rbrace.
\label{Dalpha}
\end{equation}
For $1 \le m \le n-1$
and for an anti-chain $\alpha$ in $[n] \times [m]$ of size $m$
and for an anti-chain $\beta$ in $\mathcal{D}(\alpha)$,
we set
\begin{equation}
\mathcal{R}_{\alpha,\beta} = \lbrace (M,N) \in \mathcal{O}_\alpha \times \mathcal{O}_\alpha \mid
\mathrm{Piv}(M-N) = \beta
\rbrace.
\label{Rab}
\end{equation}
For each $1 \le m \le n-1$ and
each anti-chain $\alpha$ in $[n] \times [m]$ of size $m$
and each anti-chain $\beta$ in $\mathcal{D}(\alpha)$,
consider
$M \in \mathcal{P}_n(\mathbb{F})$ with
$M_{i,j} = 1$ if $(i,j) \in \alpha \cup \beta$ and $M_{i,j} = 0$ otherwise,
and $N \in \mathcal{P}_n(\mathbb{F})$ with
$N_{i,j} = 1$ if $(i,j) \in \alpha$ and $N_{i,j} = 0$ otherwise.
Then we have $(M, N) \in \mathcal{R}_{\alpha,\beta}$
and in particular $\mathcal{R}_{\alpha,\beta} \neq \emptyset$.

\begin{prop}\label{prop:Rab}
Let $1 \le m \le n-1$
and let $\alpha$ denote an anti-chain in $[n] \times [m]$ of size $m$.
Each subset \eqref{Rab}
is an orbital of the $\mathcal{B}$-action on $\mathcal{O}_\alpha$.
\end{prop}
\begin{proof}
Immediate from Lemma \ref{lem:diag}.
\end{proof}

Let $1 \le m \le n-1$.
For an anti-chain $\alpha$ in $[n] \times [m]$ of size $m$,
consider
\begin{equation}
\mathcal{D}_1(\alpha) = \lbrace i \mid \text{there is no $k$ such that $(i,k) \in \alpha$}\rbrace.
\label{D1alpha}
\end{equation}
Then we have $|\mathcal{D}_1(\alpha)| = n-m$ since
$|\alpha| = m$.
For $1 \le i \le n-m$,
we define
$\lambda_i = |\lbrace j \mid (d_i,j) \in \mathcal{D}(\alpha)\rbrace|$,
where $d_i$ denotes the $i$-th smallest element in $\mathcal{D}_1(\alpha)$.
Then $\lambda = (\lambda_1, \lambda_2, \ldots, \lambda_{n-m}) \in \mathbb{N}^{n-m}$ is
an integer partition (i.e., a non-increasing sequence) with largest part at most $m$,
where
\[
\mathbb{N} = \lbrace 0,1,\ldots\rbrace.
\]
Consider the map $\varphi_m$ which sends $\alpha$
to $\lambda$.

For an integer partition $\lambda = (\lambda_1, \lambda_2, \ldots, \lambda_l)$,
the \emph{Ferrers board} of shape $\lambda$
is defined by
\[
\lbrace (i,j) \in \mathbb{N} \times \mathbb{N} \mid 1 \le i \le l, 1 \le j \le \lambda_i \rbrace.
\]
We endow the Ferrers board
with direct product order in $\mathbb{N} \times \mathbb{N}$.

\begin{lem}\label{al}
For $1 \le m \le n-1$,
the map $\varphi_m$ is a bijection between the following two sets:
\begin{enumerate}
\item the anti-chains in $[n] \times [m]$ of size $m$;
\item the integer partitions in $\mathbb{N}^{n-m}$ with largest part at most $m$.
\end{enumerate}
\end{lem}
\begin{proof}
Let $1 \le m \le n-1$.
It is clear that $\varphi_m$ is a map from (i) to (ii).
We define the map $\varphi_m'$ from (ii) to (i) as follows.
For a given integer partition $\lambda$ in $\mathbb{N}^{n-m}$ with largest part at most $m$,
we define
$\alpha$ as the set of maximal elements
in the Ferrers board of shape $\mu = \lambda \cup (m, m-1,\ldots,1)$,
which is the integer partition obtained by rearranging
parts of both $\lambda$ and $(m, m-1,\ldots,1)$ in non-increasing order.
Since $\mu$ is in $\mathbb{N}^{n}$ and its largest part is $m$,
$\alpha$ is an anti-chain in $[n] \times [m]$ of size $m$.
The map $\varphi_m'$ is defined to send $\lambda$
to $\alpha$.
By construction,
$\varphi_m$ and $\varphi_m'$ are inverses and hence bijections.
\end{proof}

\begin{lem}\label{bmu}
Let $1 \le m \le n-1$
and
let $\alpha$ denote an anti-chain in $[n] \times [m]$ of size $m$.
The poset
$\mathcal{D}(\alpha)$ in \eqref{Dalpha} is isomorphic to
the Ferrers board of shape $\varphi_m(\alpha)$.
Moreover,
there is a one-to-one correspondence
between
the anti-chains in $\mathcal{D}(\alpha)$
and the subpartitions of $\varphi_m(\alpha)$.
\end{lem}
\begin{proof}
Recall the set $\mathcal{D}_1(\alpha)$ in \eqref{D1alpha}.
Then define the map $\psi$
from
the Ferrers board of shape $\varphi_m(\alpha)$
to $\mathcal{D}(\alpha)$ by
$\psi(i,j) = (d_i,j)$,
where $d_i$ denotes the $i$-th smallest element in $\mathcal{D}_1(\alpha)$.
It is obvious that $\psi$ is an order-preserving bijection.
The first assertion follows.
To show the second assertion,
we define the map $\rho$ from
the subpartitions of $\varphi_m(\alpha)$
to
the anti-chains in $\mathcal{D}(\alpha)$
by
$\rho(\mu) = \psi(\max(\mu))$,
where $\max(\mu)$ is the set of all maximal elements in
the Ferrers board of shape $\mu$.
From the construction,
the map $\rho$ is also a bijection.
The second assertion follows.
\end{proof}

\begin{ex}[$n=7$, $m=4$]
Take the anti-chain $\alpha = \lbrace
(2,4), (4,3), (5,2), (7,1)
\rbrace$ as in Example \ref{ex2}.
Recall that $\mathcal{O}_\alpha$ is the set of
matrices of the form \eqref{ex:Oalpha}.
Then $\varphi_4(\alpha) = (4,3,1)$.
We remark that
each number in $(4,3,1)$
equals the number of $\ast$'s in each row without a pivot.
\end{ex}

\section{The association scheme on each Schubert cell}

For $1 \le m \le n-1$
and for an anti-chain $\alpha$ in $[n] \times [m]$ of size $m$,
by Propositions \ref{prop:OA} and \ref{prop:Rab},
the pair
\begin{equation}
\mathfrak{X}_\alpha = (\mathcal{O}_\alpha, \lbrace\mathcal{R}_{\alpha, \beta}\rbrace_\beta)
\label{Xalpha}
\end{equation}
becomes an association scheme,
where $\beta$ runs over all anti-chains in $\mathcal{D}(\alpha)$ in \eqref{Dalpha}.
See \cite[Preface]{MR2184345}.
We remark that by Lemma \ref{al},
the family of association schemes $\lbrace \mathfrak{X}_\alpha \rbrace_\alpha$
can be indexed by integer partitions $\lambda \in \mathbb{N}^{n-m}$ with largest part at most $m$.
In this case,
the associate classes of $\mathfrak{X}_\alpha = \mathfrak{X}_\lambda$ are indexed by
the subpartitions of $\lambda$
by Lemma \ref{bmu}.
\begin{thm}
Let $1 \le m \le n-1$
and let $\alpha$ denote an anti-chain in $[n] \times [m]$ of size $m$.
The association scheme $\mathfrak{X}_\alpha$ in \eqref{Xalpha} is symmetric.
\end{thm}
\begin{proof}
Immediate from the definition of $\mathcal{R}_{\alpha, \beta}$.
\end{proof}

\begin{thm}
For $1 \le m \le n-1$
and for an anti-chain $\alpha$ in $[n] \times [m]$ of size $m$,
the association scheme $\mathfrak{X}_\alpha$ in \eqref{Xalpha} is
the generalized wreath product of the one-class association schemes
with the base set $\mathbb{F}$ over the poset $\mathcal{D}(\alpha)$ in \eqref{Dalpha}.
\end{thm}
\begin{proof}
For
$M, N \in \mathcal{O}_\alpha$
and for
an anti-chain $\beta$ in $\mathcal{D}(\alpha)$,
we have
$(M,N) \in \mathcal{R}_{\alpha, \beta}$
if and only if
$M_{i,j} = N_{i,j}$ if $(i,j) \not\in \beta \cup \mathrm{Down}(\beta)$, $M_{i,j} \neq N_{i,j}$ if $(i,j) \in \beta$.
Therefore,
this associate relation
is the same as that of
the generalized wreath product of the one-class association schemes
with the base set $\mathbb{F}$ over the poset $\mathcal{D}(\alpha)$.
So the result follows.
\end{proof}

\section{Concluding remarks}

This paper focuses on the Schubert cells
of a Grassmannian.
It would be an interesting problem to find similar results on Schubert cells for other types of BN-pairs.

The Terwilliger algebra, introduced by P.~Terwilliger \cite{MR1203683},
of the wreath product of one-class association schemes
is discussed in several papers \cite{MR2610292, MR2747797, MR2802177}.
We will consider
the Terwilliger algebra of the generalized wreath product of one-class association schemes
in a future paper.

\section*{Acknowledgments}
The author gratefully acknowledges the many helpful suggestions of his advisor, Hajime Tanaka
during the preparation of the paper.
The author thanks Paul Terwilliger
for giving valuable comments
and Motohiro Ishii for drawing the author's attention to the theory of Schubert cells.

\begin{bibdiv}
\begin{biblist}

\bib{MR2206477}{article}{
   author={Bailey, R. A.},
   title={Generalized wreath products of association schemes},
   journal={European J. Combin.},
   volume={27},
   date={2006},
   number={3},
   pages={428--435},
   issn={0195-6698},
}

\bib{MR882540}{book}{
   author={Bannai, Eiichi},
   author={Ito, Tatsuro},
   title={Algebraic combinatorics. I},
   subtitle={Association schemes},
   publisher={The Benjamin/Cummings Publishing Co., Inc., Menlo Park, CA},
   date={1984},
   pages={xxiv+425},
   isbn={0-8053-0490-8},
}

\bib{MR2712017}{thesis}{
   author={Bhattacharyya, Gargi},
   title={Terwilliger algebras of wreath products of association schemes},
   type={Thesis (Ph.D.)--Iowa State University},
   publisher={ProQuest LLC, Ann Arbor, MI},
   date={2008},
   pages={84},
   isbn={978-0549-68830-3},
}

\bib{MR2610292}{article}{
   author={Bhattacharyya, Gargi},
   author={Song, Sung Y.},
   author={Tanaka, Rie},
   title={Terwilliger algebras of wreath products of one-class association
   schemes},
   journal={J. Algebraic Combin.},
   volume={31},
   date={2010},
   number={3},
   pages={455--466},
   issn={0925-9899},
}

\bib{MR2002219}{article}{
   author={Hanaki, Akihide},
   author={Hirotsuka, Kaoru},
   title={Irreducible representations of wreath products of association
   schemes},
   journal={J. Algebraic Combin.},
   volume={18},
   date={2003},
   number={1},
   pages={47--52},
   issn={0925-9899},
}

\bib{MR3013937}{collection}{
   title={Handbook of linear algebra},
   series={Discrete Mathematics and its Applications (Boca Raton)},
   editor={Hogben, Leslie},
   edition={2},
   publisher={CRC Press, Boca Raton, FL},
   date={2014},
   pages={xxx+1874},
   isbn={978-1-4665-0728-9},
}

\bib{MR2964723}{article}{
   author={Kim, Kijung},
   title={Terwilliger algebras of wreath products by quasi-thin schemes},
   journal={Linear Algebra Appl.},
   volume={437},
   date={2012},
   number={11},
   pages={2773--2780},
   issn={0024-3795},
}

\bib{MR2474907}{book}{
   author={Lakshmibai, V.},
   author={Brown, Justin},
   title={Flag varieties},
   series={Texts and Readings in Mathematics},
   volume={53},
   subtitle={An interplay of geometry, combinatorics, and representation
   theory},
   publisher={Hindustan Book Agency, New Delhi},
   date={2009},
   pages={xiv+272},
   isbn={978-81-85931-92-0},
}

\bib{MR2535398}{article}{
   author={Martin, William J.},
   author={Tanaka, Hajime},
   title={Commutative association schemes},
   journal={European J. Combin.},
   volume={30},
   date={2009},
   number={6},
   pages={1497--1525},
   issn={0195-6698},
}

\bib{MR3612420}{article}{
   author={Song, Sung Y.},
   author={Xu, Bangteng},
   author={Zhou, Shenglin},
   title={Combinatorial extensions of Terwilliger algebras and wreath
   products of association schemes},
   journal={Discrete Math.},
   volume={340},
   date={2017},
   number={5},
   pages={892--905},
   issn={0012-365X},
}

\bib{MR2868112}{book}{
   author={Stanley, Richard P.},
   title={Enumerative combinatorics. Volume 1},
   series={Cambridge Studies in Advanced Mathematics},
   volume={49},
   edition={2},
   publisher={Cambridge University Press, Cambridge},
   date={2012},
   pages={xiv+626},
   isbn={978-1-107-60262-5},
}

\bib{MR2747797}{article}{
   author={Tanaka, Rie},
   title={Classification of commutative association schemes with almost
   commutative Terwilliger algebras},
   journal={J. Algebraic Combin.},
   volume={33},
   date={2011},
   number={1},
   pages={1--10},
   issn={0925-9899},
}

\bib{MR3047011}{article}{
   author={Tanaka, Rie},
   author={Zieschang, Paul-Hermann},
   title={On a class of wreath products of hypergroups and association
   schemes},
   journal={J. Algebraic Combin.},
   volume={37},
   date={2013},
   number={4},
   pages={601--619},
   issn={0925-9899},
}

\bib{MR1203683}{article}{
   author={Terwilliger, Paul},
   title={The subconstituent algebra of an association scheme. I},
   journal={J. Algebraic Combin.},
   volume={1},
   date={1992},
   number={4},
   pages={363--388},
   issn={0925-9899},
}

\bib{W}{article}{
   author={Watanabe, Yuta},
   title={An algebra associated with a flag in a subspace lattice over a finite field and the quantum affine algebra $U_q(\widehat{\mathfrak{sl}}_2)$},
   eprint={arXiv:1709.06329 [math.CO]},
   date={2017},
}

\bib{MR2802177}{article}{
   author={Xu, Bangteng},
   title={Characterizations of wreath products of one-class association
   schemes},
   journal={J. Combin. Theory Ser. A},
   volume={118},
   date={2011},
   number={7},
   pages={1907--1914},
   issn={0097-3165},
}

\bib{MR2184345}{book}{
   author={Zieschang, Paul-Hermann},
   title={Theory of association schemes},
   series={Springer Monographs in Mathematics},
   publisher={Springer-Verlag, Berlin},
   date={2005},
   pages={xvi+283},
   isbn={978-3-540-26136-0},
   isbn={3-540-26136-2},
}
\end{biblist}
\end{bibdiv}

\bigskip

\noindent
Yuta Watanabe \\
Graduate School of Information Sciences \\
Tohoku University \\
Sendai, 980-8579 Japan\\
email: \texttt{watanabe@ims.is.tohoku.ac.jp}

\end{document}


\section{Introduction}
The general matrix element for the decay $t \to b+W^+$ including the leading
order (LO) standard model (SM) contribution is usually written as (see e.g.\
Ref.~\cite{Bernreuther:2008us})
\begin{equation}
\label{genme}
M_{tbW^+}=-\frac{g_W}{\sqrt{2}}\varepsilon^{\mu\,\ast} \bar u_b
\Big[(V^\ast_{tb}+f_L)\gamma_\mu P_L +f_R \gamma_\mu P_R
  +\frac{i\sigma_{\mu\nu}\,q^\nu}{m_W}\Big(g_L\,P_L+g_R\,P_R\Big)\Big],
\end{equation}
where $P_{L,R}=(1\mp\gamma_5)/2$. The SM structure of the $tbW^+$ vertex is
obtained by dropping all terms except for the contribution proportional to
$V^\ast_{tb}\sim1$.

The angular decay distribution for polarized top quark decay
$t(\uparrow)\to b+\ell^+ +\nu_\ell$ in the top quark rest frame is given by
\begin{eqnarray}
\label{angdis}
\frac{d\Gamma}{d\cos\theta d\phi}&=&
A+BP_t\cos\theta_P+CP_t\sin\theta_P\cos\phi+DP_t\sin\theta_P\sin\phi
\end{eqnarray}
which corresponds to the decay distribution introduced in
Refs.~\cite{Korner:1998nc,Groote:2006kq} augmented by the last $T$-odd term.
At LO of the SM one has $A=B$ and $C=D=0$. The second azimuthal term
proportional to $D$ corresponds to a $T$-odd contribution. This can be seen by
rewriting the angular factor as a triple product according to
\begin{equation}
\sin\theta_P\sin\phi=\hat p_\ell \cdot (\hat p_b \times \hat s_t)
\end{equation}
where (see Fig.~\ref{D-azi})
\begin{figure}\begin{center}
\epsfig{figure=D-azi.eps, scale=0.8}
\caption{\label{D-azi}Definitions of polar and azimuthal angles for the
process $t \to b+W^+(\to \ell^++\nu_\ell)$}
\end{center}\end{figure}
\begin{equation}
\hat p_\ell=(0,0,1) \qquad \hat p_b=(\sin\theta_b,0,\cos\theta_b)
\qquad \hat s_t=(\sin\theta_P \cos\phi,\cos\theta_P \sin\phi,\cos\theta_P)
\end{equation}
Let us repeat the arguments presented in Ref.~\cite{Groote:2006kq} that led us
to the conclusion that the equality $A=B$ already implies the vanishing of the
$T$-even azimuthal contribution $C$ confirming the LO result $C=0$. Consider
Eq.~(\ref{angdis}) for $P_t=1$ and $\sin\phi=0$, i.e.\ $\cos\phi=\pm 1$. Then
factor out the unpolarized rate term $A$. Assume first that $C/A$ is positive
and choose $\cos\phi=-1$. Expand the trigonometric functions around
$\theta_P=\pi$ for positive values of $\delta$, i.e.\
$\cos(\pi-\delta)\approx-1+\frac12\delta^2$ and
$\sin(\pi-\delta)\approx \delta$. The differential rate is then proportional to
\begin{equation}
1+\cos\theta_P-\frac{C}{A}\sin\theta_P\approx\frac12\delta^2
-\frac{C}{A}\delta=\frac{1}{2}\delta\left(\delta-\frac{2C}{A}\right).
\end{equation}
The differential rate can be seen to be negative for $\delta$ in the interval
$[0,2C/A]$. The interval can be shrunk to zero by setting $2C/A=0$, i.e.\ by
setting $C=0$. If $C/A$ is assumed to be negative, one chooses $\cos\phi=+1$
and again arrives at the conclusion $2C/A=0$ by expanding the trigonometric
functions around $\theta_P=\pi$.

The same chain of arguments but this time with $\cos\phi=0$ leads to the LO
positivity constraint for the $T$-odd structure, $D=0$.

At next-to-leading order (NLO) of QCD one no longer has $A=B$. However, the
relative difference $(A-B)/A$ is quite small which, as we will see, in turn
implies useful positivity constraints for the $T$-odd structure coefficient
$D$. As concerns the $T$-even azimuthal structure, the NLO corrections to the
LO result $C=0$ are so small that the positivity of the differential rate is
not endangered~\cite{Groote:2006kq}.

We now derive the NLO positivity constraint for the $T$-odd structure
coefficient $D$. We shall work in the approximation $m_b=0$ which implies that
the coupling terms $f_R$ and $g_L$ in Eq.~(\ref{genme}) are zero. The NLO
forms of the integrated $m_b=0$ rates are listed in
Refs.~\cite{Czarnecki:1990pe,Czarnecki:1993gt,Czarnecki:1994pu}. They read
\begin{eqnarray}
\frac{d\Gamma}{d\cos\theta d\phi}&=& (A^{(0)}+A^{(1)})\Big( 1+
\frac{A^{(0)}+B^{(1)} }{A^{(0)}+A^{(1)}}P_t\cos\theta_P \nonumber\\
&&+\frac{C}{A^{(0)}+A^{(1)}}P_t\sin\theta_P \cos\phi
+\frac{D}{A^{(0)}+A^{(1)}}P_t\sin\theta_P \sin\phi \Big),
\end{eqnarray}
where
\begin{eqnarray}
\frac{A^{(1)}}{A^{(0)}}&=&\frac{\alpha_sC_F}{4\pi}\,\frac1{(1-x^2)^2(1+2x^2)}
\bigg( (1-x^2)(5+9x^2-6x^4) \nonumber\\
&&-2(1-x^2)^2(1+2x^2)\Big[\frac{2\pi^2}{3} +4\ln(1-x^2)\ln x
  +4\mathop{\rm Li}\nolimits_2(x^2) \Big]\nonumber\\
&&-8x^2(1+x^2)(1-2x^2)\ln x -2(1-x^2)^2(5+4x^2)\ln(1-x^2)\bigg)\nonumber\\
&=&-0.0846955,
\end{eqnarray}
and
\begin{eqnarray}
\frac{B^{(1)}}{A^{(0)}}&=&\frac{\alpha_sC_F}{4\pi}\,\frac1{(1-x^2)^2(1+2x^2)}
\bigg( -(1-x^2)(15-x^2+2x^4) \nonumber\\
&&+(1-x^2)(1+x^2+4x^4)\frac{2\pi^2}{3} -2(1-x^2)^2(5+4x^2)\ln(1-x^2)\nonumber\\
&&-16x^2(2+x^2-x^4)\ln x
-16(1-x^2)(2+2x^2-x^4)\ln x \ln (1-x^2)\nonumber\\
&&-4(1-x^2)(5+5x^2-4x^4)\mathop{\rm Li}\nolimits_2(x^2)\bigg) \nonumber\\
&=&-0.0863048,
\end{eqnarray}
where $x=m_W/m_t$. Here we have also listed the numerical values for the two
ratios using $\alpha_s(m_t)=0.1062$, $m_t=173.21{\rm\,GeV}$ and
$m_W=80.385{\rm\,GeV}$~\cite{Patrignani:2016xqp}. The ratio expressions
$A^{(1)}/A^{(0)}$ and $B^{(1)}/A^{(0)}$ have been rechecked in
Ref.~\cite{Groote:2006kq}. Reference~\cite{Groote:2006kq} also contains results
on the azimuthal rate coefficient $C$. This coefficient, however, will be of
no concern in the derivation of the positivity bounds for the $T$-odd rate
coefficient $D$. In fact, setting $sin\phi=\pm1$ will eliminate the
contribution of $C$. This will be our choice.

Next we must determine the contribution of the imaginary part of the coupling
factor $g_R$ to the $T$-odd azimuthal rate term $D$. The relevant contribution
arises from the interference of the coupling factor $g_R$ with the Born term
contribution. It is for this reason that there is no $\mathop{\rm Im}\nolimits f_L$ contribution
to the $T$-odd rate coefficient $D$ since the coupling term is
self-interfering. After some algebra one finds
\begin{equation}
\frac{D}{A^{(0)}}=\frac{3\pi(1-x^2)}{4(1+2x^2)}\mathop{\rm Im}\nolimits g_R
\end{equation}
where we have only kept the contribution linear in $\mathop{\rm Im}\nolimits g_R$. Further, we
assume $\mathop{\rm Im}\nolimits g_R $ to be positive, and set $P_t=1$ and $\sin\phi=-1$. We
expand around $\theta_P=\pi$ for small positive values of $\delta$ which gives
$\cos(\pi-\delta)=-1+\frac{1}{2}\delta^2$ and $\sin(\pi-\delta)=\delta$ to
obtain
\begin{equation}\label{WthetaP}
W(\theta_P)\sim 1+(1-\Delta)\cos\theta_P
-\frac{D}{A^{(0)}+A^{(1)}}\sin\theta_P
=\Delta-\frac{D\delta}{A^{(0)}+A^{(1)}}+\frac{1-\Delta}2\delta^2,
\end{equation}
where we have defined the small quantity
\begin{equation}
\Delta=\frac{A^{\rm NLO}-B^{\rm NLO}}{A^{\rm NLO}}
  =\frac{A^{(1)}-B^{(1)}}{A^{(0)}+A^{(1)}}
=\frac{A^{(1)}-B^{(1)}}{A^{(0)}(1+A^{(1)}/A^{(0)})},
\end{equation}
keeping in mind that $A^{(0)}=B^{(0)}$. Numerically one has $\Delta=0.001758$
where the small difference to the numerical results in
Ref.~\cite{Groote:2006kq} results from having used updated values
$m_W=80.385{\rm\,GeV}$ and $m_t=173.21{\rm\,GeV}$~\cite{Patrignani:2016xqp}.

The rate proportional to $W(\theta_P)$ in Eq.~(\ref{WthetaP}) becomes negative
if the contribution proportional to $\mathop{\rm Im}\nolimits g_R$ becomes larger than the
remaining terms. However, this is no longer the case if the quadratic
equation~(\ref{WthetaP}) in $\delta$ has no real-valued zeros. The pertinent
condition for the discriminant reads
\begin{equation}
\frac{3\pi(1-x^2)}{4(1+2x^2)}|\mathop{\rm Im}\nolimits g_R|
  \le\sqrt{2\Delta(1-\Delta)}\left(1+\frac{A^{(1)}}{A^{(0)}}\right).
\end{equation}
Numerically one obtains
\begin{equation}
\mathop{\rm Im}\nolimits g_R \le\,0.0420
\end{equation}
The same chain of arguments, but now for negative values of $\mathop{\rm Im}\nolimits g_R$ and
$\sin\phi=1$ leads to $\mathop{\rm Im}\nolimits g_R \ge\,-0.0420$ such that one has the two-sided
constraint
\begin{equation}
\label{bound0}
\mathop{\rm Im}\nolimits g_R \in [-0.0420,\,0.0420]\,.
\end{equation}
The angle $\delta_0$ for which the quadratic form~(\ref{WthetaP}) becomes
degenerate can be calculated to be
$\delta_0=\sqrt{2\Delta/(1-\Delta)}=\pm 0.0593$. The corrections to the
expansion of $\cos(\pi - \delta)$ and $\sin(\pi - \delta)$ are of the
order $O(\delta^2=0.00352)$ and, therefore, quite small.

In Ref.~\cite{Fischer17} we have calculated the SM absorptive electroweak
contributions to $\mathop{\rm Im}\nolimits g_R$ with the result $\mathop{\rm Im}\nolimits g_R=-2.175\times 10^{-3}$
(see also Refs.~\cite{GonzalezSprinberg:2011kx,Arhrib:2016vts}). This value is
easily accommodated in the positivity bound~(\ref{bound0}).

The ATLAS Collaboration has recently published the bound~\cite{Aaboud:2017aqp}
\begin{equation}
\mathop{\rm Im}\nolimits g_R \in [-0.18,\,0.06]
\end{equation}
based on the analysis of sequential polarized two-body top quark decays
$t(\uparrow) \to b + W^+ (\to \ell^+ + \nu)$.
A somewhat tighter bound has been published in Ref.~\cite{Aaboud:2017yqf}
using also sequential polarized two-body top quark decays. The bound reads
\begin{equation}
\label{bound2}
\mathop{\rm Im}\nolimits\left(\frac{g_R}{f_L}\right) \in [-0.07,\,0.06]
\end{equation}
which we translate into a bound on $\mathop{\rm Im}\nolimits g_R$ by substituting the LO result
$f_L=1$ in Eq.~(\ref{bound2}). Both bounds are not far away from the
positivity bound on $\mathop{\rm Im}\nolimits g_R$ derived in this paper. Using the same chain of
arguments one can establish the corresponding bound for $m_b\ne 0$. Using NLO
$m_b\ne 0$ results from Ref.~\cite{Czarnecki:1994pu} on the unpolarized and
polarized rate functions $A^{(1)}$ and $B^{(1)}$ we find that for $m_b=4.8{\rm\,GeV}$
the bound is marginally strengthened to $\mathop{\rm Im}\nolimits g_R\in[-0.0418,0.0418]$. The
condition for obtaining this bound reads
\begin{equation}
\frac{3\pi(1-x^2+y^2)\sqrt\lambda}{4\left(\lambda+3x^2(1-x^2+y^2)\right)}
  |\mathop{\rm Im}\nolimits g_R|\le\sqrt{2\Delta(1-\Delta)}\left(1+\frac{A^{(1)}}{A^{(0)}}\right),
\end{equation}
where $y=m_b/m_t$ and $\lambda=\lambda(1,x^2,y^2)=1+x^4+y^4-2x^2-2y^2-2x^2y^2$
is the K\"all\'en function. $A^{(0)}$, $A^{(1)}$ and the small quantity
$\Delta$ are evaluated for $m_b\ne 0$.

\subsection*{Acknowledgments}
We would like to thank J.~Mueller for discussions.
This work was supported by the Estonian Science Foundation under Grant
No.~IUT2-27. S.G.\ acknowledges the hospitality of the theory group THEP at
the Institute of Physics at the University of Mainz and the support of the
Cluster of Excellence PRISMA at the University of Mainz.


\section{Introduction} \label{intro}

In this paper, we prove the following:
\begin{thm}\label{thm:main}
Let $p \ge 5$ be a prime. The only integer solutions to the equation
\begin{equation}\label{eq:main}
(x-r)^3 + x^3 + (x+r)^3 = y^p
\qquad \gcd(x,r)=1
\end{equation}
with $0 < r\le 10^6$ are the trivial ones satisfying $xy=0$.
\end{thm}
The restriction $\gcd(x,r)=1$ is natural one, for
otherwise it is easy to construct artificial solutions by
scaling.
Our approach combines Frey curves and modularity with
a result of Mignotte based on linear form in logarithms
and various elementary techniques. This approach can easily
be extended to treat larger ranges of the parameter $r$.

Both our theorem and our methods extend theorems and ideas already found
in the literature. For example,
Cassels \cite{Cassels} and Uchiyama \cite{Uchiyama} independently,
and using only elementary arguments, solved \eqref{eq:main}
with $p=2$ and $r=1$.
Zhongfeng Zhang \cite{Zhang} extended the work of Cassels and Uchiyama and determined all perfect powers that
sums of three consecutive cubes; that is he solved \eqref{eq:main} with $r=1$ and  $p\geq 2$.

Power values of sums of powers of consecutive integers, or consecutive terms in
an arithmetic progression, have been considered by many authors from different
points of view. Techniques from algebraic number
theory and Diophantine approximation have allowed the resolution of such equations with small exponents, as well as proofs of general
theorems.
This can be seen in the work of  B\'{e}rces,  Pink,  Sava\c{s} and  Soydan \cite{BPSSoydan}, Brindza \cite{Br}, Cassels \cite{Cassels}, Gy\H{o}ry, Tijdeman and  Voorhoeve \cite{GTV},
Hajdu \cite{Ha}, Patel \cite{Patelsquares}, Patel and Siksek \cite{PatelSiksek}, Pint\'er \cite{P}, Sch\"affer \cite{Sc}, Schinzel and Tijdeman \cite{ScTi}, Soydan \cite{Soydan}, Urbanowicz \cite{U1},
and Zhang and Bai \cite{ZB} among many others.

Contributions from Bennett, Gy\H{o}ry and Pint\'er \cite{BGP}, Bennett Patel and Siksek \cite{BPS1}, \cite{BPS2} Pint\'er  \cite{P2} and Zhang \cite{Zhang},
resolved many such equations using Frey--Hellegouarch curves and Galois representations.

Recently, the equations
\[
(x-r)^k+x^k+(x+r)^k=y^n
\]
have been studied for $k=2$ by Koutsianas \cite{Koutsianas}, Koutsianas and Patel \cite{KoutsianasPatel} and for $k=4$ by Zhongfeng Zhang \cite{Zhang2}, via the study of Frey-Hellegouarch curves and Galois representations.

\section{Small exponents: $p=2$ and $3$} \label{smallp}

For $p=2$, we have the following subfamily of equation~\eqref{eq:main}:
\[
(x-r)^3 + x^3 + (x+r)^3 = y^2.
\]
For a fixed value of $r$,
since this defines an elliptic curve,
we may find all integer points via an appropriate computer package (for example \texttt{Magma}).

As $r$ varies, we immediately see that there are infinite families of solutions. One such family can be described by the following equations:
\begin{align*}
x &= 24a^2b^2, \\
r &= 9a^4 - 8b^4, \\
y &= 12ab(9a^4 - 8b^4),
\end{align*}
where $a$ and $b$ are integers. We easily see that $(x,r,y)$ are integral and satisfy the equation $(x-r)^3 + x^3 + (x+r)^3 = y^2$.

\medskip

For $p=3$, we consider the following subfamily of equation~\eqref{eq:main}:
\[
(x-r)^3 + x^3 + (x+r)^3 = y^3.
\]
We can define a smooth curve $C$ in $\mathbb{P}^2$:
\[
C \; : \; 3x(x^2+2r^2) = y^3.
\]
As this is a plane cubic, it is a curve of genus $1$. Let
$P=(x,y,r)=(0,0,1) \in C(\mathbb{Q})$. Then $C$ is isomorphic
to an elliptic curve $E$ where the isomorphism $\phi \; : C \rightarrow E$
takes $P$ to the point at infinity on $E$. Using \texttt{Magma},
we find that $E$ has the model
\[
E \; : \; Y^2=X^3-648 \, .
\]
The map $\phi$ is given by $\phi(x:y:r)=(6y/x^3,-36r/x^2)$.
We find that $E({\mathbb Q}) \cong {\mathbb Z}^2$ with Mordell--Weil basis
given by $(18,72)$ and $(9,9)$.

In particular, as $r$ varies, both cases $p=2$ and $p=3$ of equation~\eqref{eq:main} give rise to  infinitely many integer solutions.

\section{Four cases: equations of signature $(p,p,p)$ and $(p,p,2)$} \label{cases}

We may rewrite equation~\eqref{eq:main} as:
\[
3x(x^2+2r^2) = y^p \quad x,r,y,p \in {\mathbb Z}, \quad p\geq 2, \quad \gcd(x,r) = 1.
\]
Since $3 \mid y$ we can write $y=3w$. Thus
\begin{equation} \label{eq:mainfact}
x(x^2+2r^2)=3^{p-1} w^p.
\end{equation}
Note that $\gcd(x,x^2+2r^2)=1$ or $2$ according to
whether $2 \mid x$ or $2 \nmid x$.
Thus, we consider four cases and apply a simple descent argument in each case
to obtain ternary equations.

\begin{itemize}
\item {\bf Case 1:} $2 \mid x$ and $3 \mid x$.
We rewrite equation~\eqref{eq:mainfact} as:
\[
x=2^{p-1} \cdot 3^{p-1} w_1^p, \qquad x^2+2r^2=2 w_2^p.
\]
Substituting and simplifying we obtain:
\begin{equation}\label{eqn:case1}
w_2^p-2^{2p-3} \cdot 3^{2p-2} w_1^{2p}=r^2.
\end{equation}
\item {\bf Case 2:} $2 \mid x$ and $3 \nmid x$.
We rewrite equation~\eqref{eq:mainfact} as:
\[
x=2^{p-1} w_1^p, \qquad x^2+2r^2= 2 \cdot 3^{p-1} w_2^p.
\]
Substituting and simplifying we obtain:
\begin{equation}\label{eqn:case2}
3^{p-1} w_2^p-2^{2p-3} w_1^{2p}=r^2.
\end{equation}
\item {\bf Case 3:} $2 \nmid x$ and $3 \mid x$.
We rewrite equation~\eqref{eq:mainfact} as:
\[
x=3^{p-1} w_1^p, \qquad x^2+2r^2= w_2^p.
\]
Substituting and simplifying we obtain:
\begin{equation}\label{eqn:case3}
w_2^p-3^{2p-2} w_1^{2p}=2r^2.
\end{equation}
\item {\bf Case 4:} $2 \nmid x$ and $3 \nmid x$.
We rewrite equation~\eqref{eq:mainfact} as:
\[
x=w_1^p, \qquad x^2+2r^2= 3^{p-1} w_2^p.
\]
Substituting and simplifying we obtain:
\begin{equation}\label{eqn:case4}
3^{p-1} w_2^p-w_1^{2p}=2r^2.
\end{equation}

\end{itemize}

Cases 1 and 2 give rise to ternary equations of signature $(p,p,2)$ in the variables $(x,r,y,p)$.
Here, we are able to prove Theorem~\ref{thm:main} without any restriction on $r$
using Chabauty's method for $p=5$ and the modular approach for $p \geq 7$.
This is expanded upon in Sections~\ref{sec:p=5} and \ref{sec:case12} respectively.

Cases 3 and 4 also produce equations of signature $(p,p,2)$ in the variables
$(x,r,y,p)$, however, these equations currently don't seem tractable under the modular method.
Therefore, we consider their ternary equation as having signature $(p,p,p)$ in the variables $(x,y,p)$
when fixing values of $r$. Moving away from the modular method and taking a more classical approach yields
results for $r \leq 10^6$ with minimum computational effort. This is expanded upon in
the remaining sections.

\section{The case $p=5$}
\label{sec:p=5}

For $p=5$,
equations \eqref{eqn:case1}, \eqref{eqn:case2}, \eqref{eqn:case3}, \eqref{eqn:case4} can
be rewritten as curves of genus $2$ defined over ${\mathbb Q}$.
For a curve $C$ of genus $2$ defined over ${\mathbb Q}$, the method of Chabauty and Coleman
can be used to determine the rational points provided the rank of the Mordell--Weil
group of the Jacobian is at most $1$. For an introduction to the method
of Chabauty and Coleman we recommend \cite{MP} and \cite{SiksekChab}.
Fortunately the method is implemented in the \texttt{Magma} software
package and we were able to solve equations \eqref{eqn:case1}, \eqref{eqn:case2}
and \eqref{eqn:case3} for $p=5$ using this implementation.

\subsection{Case 1}
Let $p=5$ and consider \eqref{eqn:case1}. Write
$X=w_2/w_1^2$ and $Y=r/w_1^5$. Then $(X,Y)$ is a rational
point on the genus $2$ curve
\[
C \; : \; Y^2=X^5-2^7 \cdot 3^8.
\]
Let $J$ be the Jacobian on $C$. Using \texttt{Magma}, which uses
a $2$-descent algorithm of Stoll \cite{Stoll}, we find that the Mordell--Weil
rank of $J$ is $0$. The \texttt{Magma} implementation of Chabauty gives
 $C({\mathbb Q})=\{\infty\}$. Thus $w_1=0$ and $w_2=1$, which gives the
trivial solution to \eqref{eq:main}.

\subsection{Case 2}
We now consider \eqref{eqn:case2} with $p=5$. With $X$, $Y$ as above we see
that $(X,Y)$ is a rational point on
\[
C \; : \; Y^2=3^4 X^5-2^7.
\]
The Mordell--Weil group of the Jacobian is infinite cyclic,
generated by the point, given in Mumford representation,
\[
\left(X^2 + \frac{20}{81} X + \frac{8}{9}, \; \frac{1124}{81} X + \frac{32}{9} \right).
\]
The \texttt{Magma} Chabauty implementation again gives $C({\mathbb Q})=\{\infty\}$
and this gives the trivial solution to \eqref{eq:main}.

\subsection{Case 3}
Let $p=5$ and consider \eqref{eqn:case3}.
We let $X=w_2/w_1^2$ and $Y=2r/w_1^5$. Then $(X,Y)$ is a rational point
on
\[
C \; : \; Y^2=2(X^5-3^8).
\]
Let $P=(9,324) \in C({\mathbb Q})$. The Mordell--Weil group of the Jacobian
is infinite cyclic, generated by $P-\infty$. The Chabauty implementation
now gives
$\infty$, $(9,324)$, $(9,-324)$ as the only rational points. Thus
$w_1=0$, $w_2=1$ or $w_1=1$, $w_2=9$, $r=\pm 162$. The former gives a
trivial solution to \eqref{eq:main} and the latter can be discarded
as it does not
satisfy $\gcd(x,r)=1$.

\subsection{Case 4}
We mention in passing that for $p=5$ the Mordell--Weil rank of the Jacobian
of the genus $2$ curve corresponding to equation~\eqref{eqn:case4} is $2$, and so
Chabauty's method is inapplicable.

\section{Cases 1 and 2: a modular approach} \label{sec:case12}

In this section, we are able to give a complete resolution to Cases 1 and 2 for any
$r \in {\mathbb Z}$ and $p \ge 7$.
In order to achieve this,  we make use of the work of Bennett and Skinner
\cite{BenS} (see also \cite{Siksek}). That work
builds upon the celebrated modularity theorem of Wiles, Breuil, Conrad,
Diamond and Taylor \cite{Wiles}, \cite{modularity}, on Ribet's level lowering theorem
\cite{Ribet}, and on Mazur's isogeny theorem \cite{Mazur}.

\begin{lem}\label{lem:cases12mod}
Let $p\geq 7$ be a prime. Suppose $2 \mid x$. Then equation~\eqref{eq:main} has no integer solutions.
\end{lem}
The reader will note that the condition $2 \mid x$ corresponds to Cases 1 and 2 above.
\begin{proof}
Suppose first that $2 \mid x$ and $3 \mid x$ and thus we are in Case 1. Thus we have a
non-trivial solution to a ternary equation of signature $(p,p,2)$, namely \eqref{eqn:case1}.
Our condition $\gcd(x,r)=1$ readily implies that the three terms in \eqref{eqn:case1}
are coprime.
Changing the sign of $r$ if necessary, we may suppose that
$r \equiv 1 \pmod{4}$.
Using the recipes of Bennett and Skinner \cite{BenS}
we attach to this the Frey-Hellegouarch curve:
\[
E \; : \; Y^2+XY=X^3+\left(\frac{r-1}{4} \right) X^2-2^{2p-9} \cdot 3^{2p-4} w_1^{2p} X.
\]
By \cite[Lemma 3.3]{BenS},
$E$ arises modulo $p$ from a newform of weight $2$ and level $6$.
As there are no such newforms we have a contradiction.

Hence, we suppose that $2 \mid x$ and $3 \nmid x$ (Case 2). We now consider the
ternary equation \eqref{eqn:case2} where again the three terms are coprime
and again we may suppose $r \equiv 1 \pmod{4}$.
The recipes
of Bennett and Skinner \cite{BenS} attach to this the following Frey-Hellegouarch curve:
\[
E \; : \; Y^2+XY=X^3+\left( \frac{r-1}{4}\right) X^2-2^{2p-9} w_1^{2p} X.
\]
This too arises from a newform of weight $2$ and level $6$, and again, we arrive at a contradiction.
\end{proof}

\section{Cases 3 and 4: Bounding $p$} \label{Boundp}

The aim of this section is to bound $p$ in Cases 3 and 4.
We require the following theorem of Mignotte, which is found in \cite[Chapter~12, p.~423]{Cohen2}.

\begin{thm}[Mignotte]\label{thm:Mignotte}
Assume that the exponential Diophantine inequality
\[
|ax^n - by^n | \leq c, \quad \text{ with } a,b,c \in {\mathbb Z}_{\ge 0} \text{ and }  a\neq b,
\]
has a solution in strictly positive integers $x$ and $y$ with $\max\{x,y\} > 1$. Let $A = \max \{a,b,3\}$. Then
\[
n \leq \max \left\{ 3 \log\left(1.5| c/b| \right),\; \frac{7400\log A}{\log\left(1+\log A/\log\left(\mid a/b\mid\right)\right)}\right\}.
\]
\end{thm}

A simple application of Mignotte's Theorem yields workable bounds for
$p$ for equations~\eqref{eqn:case3} and
\eqref{eqn:case4}, and hence for Cases 3 and 4 respectively.
We first deal with Case 3.
\begin{lem}
Suppose $2 \nmid x$ and $3 \mid x$ and $ |r| \leq 1.55\times 10^{1697}$.
If equation~\eqref{eq:main} has an integral non-trivial solution, then $p<24,000$.
\end{lem}
\begin{proof}
Multiplying \eqref{eqn:case3} by $9$ gives
\begin{equation}\label{eqn:Linlog3}
9w_2^p - (3 w_1)^{2p}=18r^2.
\end{equation}
We apply Theorem~\ref{thm:Mignotte} to equation~\eqref{eqn:Linlog3} with $a=9$, $b=1$, $c=18r^2$ and we obtain
\begin{align*}
p &\leq \max \left\{ 3 \log\left(3^3 \cdot r^2  \right),\; 7400\log 9 / \log 2 \right\}\\
 & \leq \max \left\{ 23457.4384\ldots ,\; 23457.4450\ldots \right\}.
\end{align*}
We deduce that $p < 24,000$, which completes the proof.
\end{proof}

Next we deal with Case 4.
\begin{lem}
Suppose $2 \nmid x$ and $3 \nmid x$. Suppose $ |r| \leq 2.99\times 10^{848}$.
If equation~\eqref{eq:main} has an integral non-trivial solution, then $p\leq 12,000$.
\end{lem}
\begin{proof}
Multiplying \eqref{eqn:case4} by $3$ gives
\begin{equation}\label{eqn:Linlog4}
(3w_2)^p- 3\cdot w_1^{2p}=6r^2.
\end{equation}
We apply Theorem~\ref{thm:Mignotte} to equation~\eqref{eqn:Linlog4} with $a=3$, $b=1$, $c=6r^2$
to obtain
\begin{align*}
p &\leq \max \left\{ 3 \log\left( 3^2\cdot r^2  \right),\; 7400\log 3 / \log 2 \right\}\\
 & \leq \max \left\{ 11727.7611\ldots ,\; 11728.7225\ldots \right\}.
\end{align*}
Hence,if equation~\eqref{eq:main} has an integral non-trivial solution, then $p < 12,000$.
\end{proof}

\section{Criteria for eliminating equations of signature $(p,2p,2)$} \label{sec:criteria}

To prove Theorem~\ref{thm:main} we would like to solve equation \eqref{eq:main}
for $0 < r \le 10^6$. We have seen that it is enough to solve equations \eqref{eqn:case3}
and \eqref{eqn:case4} for $0<r \le 10^6$. In the previous section we saw that
the assumption $0<r\leq 10^6$ forces
$p <24000$ in \eqref{eqn:case3} and $p<12000$ in \eqref{eqn:case4}.
This leaves us with a finite, albeit very large, number of equations
to solve. These all have form
\begin{equation}\label{eqn:tern}
ay_2^{p} - by_1^{2p} = c
\end{equation}
where $p$ is a fixed odd prime and $a$, $b$ and $c$ are positive integers
satisfying $\gcd(a,b,c)=1$.
In this section, we state various criteria for equations of the form \eqref{eqn:tern}
to have no integer solutions.
Proofs of these criteria can be found  in \cite{BPS2} and \cite{Patelthesis}; the latter providing extensive details
and motivation.
In later sections, we subject equations \eqref{eqn:case3} and \eqref{eqn:case4}
to these criteria.

\subsection{Criteria for the non-existence of solutions} \label{SG}

\begin{lem}\label{lem:Sophiecriterion}
Let $p \ge 3$ be a prime. Let $a$, $b$ and $c$ be positive integers such that
$\gcd(a,b,c)=1$. Let $q=2k p+1$ be a prime that does
not divide $a$. Define
\begin{equation}\label{eqn:mu}
\mu(p,q)=\{ \eta^{2p} \; : \; \eta \in {\mathbb F}_q \}
=\{0\} \cup \{ \zeta \in {\mathbb F}_q^* \; : \; \zeta^{k}=1\}
\end{equation}
and
$$
B(p,q)=\left\{ \zeta \in \mu(p,q) \; : \; ((b \zeta+c)/a)^{2k} \in \{0,1\} \right\} \, .
$$
If $B(p,q)=\emptyset$,
then equation~\eqref{eqn:tern}
does not have integral solutions.
\end{lem}

The proof is straightforward and can be found in \cite{BPS2} and \cite{Patelthesis}.

\subsection{Local Solubility}\label{sub:locsol}
In this subsection, we appeal to classical local solubility methods
for Diophantine equations in order
to conclude non-existence of solutions for many tuples $(a,b,c,p)$
in equation~\eqref{eqn:tern}.

In \eqref{eqn:tern}, recall the condition  $\gcd(a,b,c)=1$.
Let $g=\Rad(\gcd(a,c))$ and suppose that $g > 1$. Then $g \mid y_1$,
and we can write $y_1=g y_1^\prime$. Thus
\[
a y_2^p- b g^{2 p} {y_1^\prime}^{2p}=c.
\]
Removing a factor of $g$ from the coefficients,
we obtain
\[
a^\prime {y_2}^p - b^\prime {y_1^\prime}^{2p}=c^\prime,
\]
where $a^\prime=a/c$ and $c^\prime=c/g<c$. Similarly, if $h=\gcd(b,c)>1$, we obtain
\[
a^\prime {y_2^\prime}^p - b^\prime {y_1}^{2p}=c^\prime,
\]
where $c^\prime=c/h<c$. Applying these operations repeatedly,
we arrive at an equation of the form
\begin{equation}\label{eqn:predescent}
A \rho^p-B \sigma^{2p}=C
\end{equation}
where $A$, $B$, $C$ are now pairwise coprime.
A necessary condition for the existence of solutions is that
for any odd prime $q \mid A$, the residue $-BC$ modulo $q$ is a square.
Besides this basic test, we also check for local solubility
at the primes dividing $A$, $B$, $C$, and all primes $q \le 19$.

\subsection{A Descent}\label{sub:furtherdesc}

If local techniques of Section~\ref{sub:locsol} fail to rule out solutions
to equation \eqref{eqn:tern} for particular coefficients and exponent $(a,b,c,p)$
then we may perform a further descent to rule out solutions.
With $A$, $B$, $C$ as in \eqref{eqn:predescent}
we let
\[
B^\prime=\prod_{\text{$\ord_q(B)$ is odd}} q.
\]
Thus $B B^\prime=v^2$. Write $A B^\prime=u$ and $C B^\prime=m n^2$
with $m$ squarefree. Rewrite \eqref{eqn:predescent}
as
\[
(v \sigma^p+n \sqrt{-m})(v \sigma^p-n \sqrt{-m})=u \rho^p.
\]
Let $K={\mathbb Q}(\sqrt{-m})$ and ${\mathcal O}$ be its ring of integers.
Let $\mathfrak{S}$ contain the prime ideals of ${\mathcal O}$ that
divide $u$ or $2n \sqrt{-m}$.
Clearly
$(v\sigma^p+n \sqrt{-m}) {K^*}^p$ belongs to the ``$p$-Selmer group''
\[
K(\mathfrak{S},p)=\{\epsilon \in K^*/{K^*}^p \; : \;
\text{$\ord_\mathcal{P}(\epsilon) \equiv 0 \mod{p}$ for all $\mathcal{P} \notin \mathfrak{S}$}
\}.
\]
This is an ${\mathbb F}_p$-vector space of finite dimension
(\cite[Proof of Proposition VIII.1.6]{Silverman}), and
can be computed by \texttt{Magma} using the command \texttt{pSelmerGroup}.
Let
\[
\mathcal{E}=\{ \epsilon \in K(\mathfrak{S},p) \; : \; \Norm(\epsilon)/u \in {{\mathbb Q}^*}^p \}.
\]
It follows that
\begin{equation}\label{eqn:furtherdescent}
v \sigma^p+n \sqrt{-m}=\epsilon \eta^p,
\end{equation}
where $\eta \in K^*$ and $\epsilon \in \mathcal{E}$.

\begin{lem}\label{lem:valuative}
Let $\mathfrak{q}$ be a prime ideal of $K$. Suppose one of the following
holds:
\begin{enumerate}
\item[(i)] $\ord_\mathfrak{q}(v)$, $\ord_\mathfrak{q}(n\sqrt{-m})$, $\ord_\mathfrak{q}(\epsilon)$
are pairwise distinct modulo $p$;
\item[(ii)] $\ord_\mathfrak{q}(2v)$,
$\ord_\mathfrak{q}(\epsilon)$, $\ord_\mathfrak{q}(\overline{\epsilon})$
are pairwise distinct modulo $p$;
\item[(iii)] $\ord_\mathfrak{q}(2 n \sqrt{-m})$,
$\ord_\mathfrak{q}(\epsilon)$, $\ord_\mathfrak{q}(\overline{\epsilon})$
are pairwise distinct modulo $p$.
\end{enumerate}
Then there is no $\sigma \in {\mathbb Z}$
and $\eta \in K$ satisfying \eqref{eqn:furtherdescent}.
\end{lem}

Proof and further details may be found in \cite{BPS2} and \cite{Patelthesis}.

\begin{lem}\label{lem:furtherdescent}
Let $q=2k p+1$ be a prime.
Suppose $q{\mathcal O}=\mathfrak{q}_1 \mathfrak{q}_2$ where $\mathfrak{q}_1$, $\mathfrak{q}_2$
are distinct, and such that $\ord_{\mathfrak{q}_j}(\epsilon)=0$
for $j=1$, $2$. Let
\[
\chi(p,q)=\{ \eta^p \; : \; \eta \in {\mathbb F}_q \}.
\]
Let
\[
C(p,q)=\{\zeta \in \chi(p,q) \; : \;
((v \zeta+n\sqrt{-m})/\epsilon)^{2k} \equiv \text{$0$ or $1 \mod{\mathfrak{q}_j}$
for $j=1$, $2$}\}.
\]
Suppose $C(p,q)=\emptyset$.
Then there is no $\sigma \in {\mathbb Z}$
and $\eta \in K$ satisfying \eqref{eqn:furtherdescent}.
\end{lem}

The proof of Lemma~\ref{lem:furtherdescent} is a simple modification
of the proof of Lemma~\ref{lem:Sophiecriterion}, the latter can be found in \cite{BPS2} and \cite{Patelthesis}.

\begin{comment}
\section{Applications}\label{sec:App}

In this section, we subject equations~\eqref{eqn:case1},~\eqref{eqn:case2},~\eqref{eqn:case3} and \eqref{eqn:case4} to the tests of Section~\ref{sec:criteria} to hopefully eliminate many tuples $(a,b,c,p)$.
We record the results of the computation for each case separately, and each table records the data as follows:
the second column
gives a breakdown of the number of equations surviving the criterion of Lemma~\ref{lem:Sophiecriterion}, according to the exponent $p$;
the third column gives the number of these equations surviving the
local solubility tests of Section~\ref{sub:locsol}; the fourth
column gives the number that then survive the further descent
of Section~\ref{sub:furtherdesc}.
\section{Application to Cases 1 and 2}\label{sec:Appcase12}

We know that equations \eqref{eqn:case1} and \eqref{eqn:case2} have no solutions with $p \ge 7$
thanks to Section~\ref{sec:case12}. We still have consider that case $p=5$.
In view of our final aim of proving Theorem~\ref{thm:main} we restrict to $r<10^6$. This
gives us two million equations to solve of the form \eqref{eqn:case1} and \eqref{eqn:case2}.
We can rule out some of these straightaway. Indeed,
Equation~\eqref{eqn:case1} is under the assumption that $2,3 \mid x$ and $\gcd(x,r)=1$. Hence $2,3 \nmid r$.  Thus we only
consider the values $r <10^6$ also satisfying $2$, $3 \nmid r$.
We wrote a \texttt{Magma} implementation of the tests of the previous section which we applied
to equation \eqref{eqn:case1} with $p=5$ and the remaining values of $r$. We found that
$63,013$ survived the criterion of Lemma~\ref{lem:Sophiecriterion}. Of these only
$17,569$ survived the local solubility tests from Section~\ref{sub:locsol},
and of these only $8$ survived the further descent (Lemmas~\ref{lem:valuative} and~\ref{lem:furtherdescent}).
The full time of computations took roughly $18$ hours. The remaining $8$ equations will be considered in Section~\ref{Thue}.

We also consider Equation~\eqref{eqn:case2} with $p=5$ and $r<10^6$. Here again we know that $2,3 \nmid r$.
We found that
$17,439$ values of $r$ survived the criterion of Lemma~\ref{lem:Sophiecriterion},
and of these  $7,599$ survived the local solubility tests and of these $28$ survived the further descent.
The $28$ remaining equations will also be dealt with in Section~\ref{Thue}.
\end{comment}

\begin{comment}
\begin{table}[!htbp]
\centering

\begin{tabu}{|c|c|c|c|}
\hline
\multirow{3}{*}{Exponent $p$} & equations~\eqref{eqn:case1} surviving  & number surviving  & number surviving\\
				 &  Lemma~\ref{lem:Sophiecriterion}&   after local & after further\\
				 &  with exponent $p$ & solubility tests  & descent\\
\hline\hline
$5$ & $63,013$ &  $17,569$ & $8$ \\
\hline\hline
\end{tabu}
}
\vskip1ex
\caption{}
\label{table:1}
\end{table}
\section{Application to Case 2}\label{sec:Appcase2}

In this section, we subject equation~\eqref{eqn:case2} with $p = 7$ and $r \leq 10^6$ to Lemma~\ref{lem:Sophiecriterion}, the local solubility tests from Section~\ref{sub:locsol} as well as lemmata~\ref{lem:valuative} and \ref{lem:furtherdescent}. The full time of computations took roughly $6$ hours.

\begin{table}[!htbp]
\centering

\begin{tabu}{|c|c|c|c|}
\hline
\multirow{3}{*}{Exponent $p$} & equations~\eqref{eqn:case2} surviving  & number surviving  & number surviving\\
				 &  Lemma~\ref{lem:Sophiecriterion} &   after local & after further\\
				 &  with exponent $p$ & solubility tests  & descent\\
\hline\hline
$5$ & $17,439$ &  $7,599$ & $28$ \\
\hline\hline
\end{tabu}
}
\vskip1ex
\caption{}
\label{table:2}
\end{table}

\end{comment}

\section{Application to Case 3}\label{sec:Appcase3}

We wrote a \texttt{Magma} program that subjected equation~\eqref{eqn:case3} with $7 \le p<24,000$
and $0 <r\leq 10^6$
to Lemma~\ref{lem:Sophiecriterion}, local solubility from subsection~\ref{sub:locsol} as well as lemmata~\ref{lem:valuative} and \ref{lem:furtherdescent}.
Recall the assumption that $2 \nmid x,\; 3 \mid x$ and $\gcd(x,r)=1$ for Case 3. Thus $3 \nmid r$.  This allowed us to reduce the number of values of $r$ that we subject to the tests of Section~\ref{sec:criteria}.
The full time of computations took roughly $66$ hours.
The results are summarised in Table~\ref{table:3}.

\begin{table}[!htbp]
\centering

\begin{tabu}{|c|c|c|c|}
\hline
\multirow{3}{*}{Exponent $p$} & equations~\eqref{eqn:case3} surviving  & number surviving  & number surviving\\
				 &  Lemma~\ref{lem:Sophiecriterion} &   after local & after further\\
				 &  with exponent $p$ & solubility tests  & descent\\
\hline\hline
$7$ & $67,600$  &  $29,201$ & $4$ \\
\hline
$11$ & $7,233$ &  $5,334$ &  $2$\\
\hline
$13$ & $1,718$  &  $808$ & $2$  \\
\hline
$17$ & $848$  &   $453$ & $1$  \\
\hline
$19$ & $2,032$  & $1,023$ &  $1$  \\
\hline
$23$ & $19$ &  $12$ &  $1$\\
\hline
$29$ & $8$  &  $5$ & $1$  \\
\hline
$31$ & $47$  &   $24$ & $1$  \\
\hline
$37$ & $5$  & $5$ &  $1$  \\
\hline
$41 \le p <24,000$ & $0$ & $0$ & $0$\\
\hline\hline
Total & $79,510$ & $36,865$ & $14$ \\
\hline\hline
\end{tabu}
}
\vskip1ex
\caption{}
\label{table:3}
\end{table}

\section{Application to Case 4}\label{sec:Appcase4}

We subjected equation~\eqref{eqn:case4} with $5 \le p<12,000$ and $0 < r\leq 10^6$ to
Lemma~\ref{lem:Sophiecriterion}, local solubility from Section~\ref{sub:locsol}
as well as lemmata~\ref{lem:valuative} and \ref{lem:furtherdescent}.
Equation~\eqref{eqn:case4} is under the assumption that $2,3 \nmid x$ and $\gcd(x,r)=1$.
Since $x^2 + 2r^2 = 3^{p-1}w_1^p$, we easily conclude that $3 \nmid r$.
This allowed us to significantly cut down the number of values of $r$ that we subject to the tests
of Section~\ref{sec:criteria}.
The full time of computations took roughly $70$ hours.
The results are summarised in Table~\ref{table:4}.

\begin{table}[!htbp]
\centering

\begin{tabu}{|c|c|c|c|}
\hline
\multirow{3}{*}{Exponent $p$} & equations~\eqref{eqn:case4} surviving  & number surviving  & number surviving\\
				 &  Lemma~\ref{lem:Sophiecriterion} &   after local & after further\\
				 &  with exponent $p$ & solubility tests  & descent\\
\hline\hline
$5$ & $219,921$ &  $65,821$ & $777$ \\
\hline
$7$ & $25,308$  &  $12,994$ & $0$ \\
\hline
$11$ & $3,201$ &  $1,834$ &  $0$\\
\hline
$13$ & $948$  &  $439$ & $0$  \\
\hline
$17$ & $385$  &   $201$ & $0$  \\
\hline
$19$ & $1,825$  & $1,108$ &  $0$  \\
\hline
$23$ & $7$ &  $3$ &  $0$\\
\hline
$29$ & $1$  &  $1$ & $0$  \\
\hline
$31$ & $3$  &   $2$ & $0$  \\
\hline
$37$ & $3$  & $2$ &  $0$  \\
\hline
$41$ & $1$  & $1$ &  $0$  \\
\hline
$43 \le p <12,000$ & $0$ & $0$ & $0$\\
\hline\hline
Total & $251,603$ & $82,406$ & $777$ \\
\hline\hline
\end{tabu}
}
\vskip1ex
\caption{}
\label{table:4}
\end{table}

\section{Modularity, Thue equations and completing the proof} \label{Thue}

We notice that in Table~\ref{table:3}, equations remain to be solved with an unusually high exponent $p$, and thus Magma's Thue Solver will not yield results.

On closer inspection, one notices that the tuples $(a,b,c)$ that remain (which describe equation~\eqref{eqn:tern}) are of the form $(a,b,2^k)$ with $k \geq 5$. Hence, these cases may be eliminated by once again, appealing to the modular method.

\begin{lem}\label{lem:modpow2}
Suppose $p\geq 5$ is a prime. Suppose $2 \nmid x$ and $3 \mid x$. Suppose $r = 2^k$ with $k \geq 2$. Then equation~\eqref{eq:main} has no integer solutions.
\end{lem}

\begin{proof}
Since $2 \nmid x$ and $3 \mid x$, recall that we obtain the ternary equation~\eqref{eqn:case3} of signature $(p,p,p)$:
\[
w_2^p-3^{2p-2} w_1^{2p}=2^{2k+1},
\]
with all three terms coprime (this follows readily from our imposed condition $\gcd(x,r)=1$).
Following Kraus \cite{KrausFermat}, we attach to this the Frey-Hellegouarch curve
\[
E \; : \; Y^2 = X(X - 3^{2p-2}w_1^{2p} )(X + 2^{2k+1}).
\]
By \cite[Section 4]{KrausFermat}
(which again builds on the aforementioned deep theorems of Mazur, Ribet,
Wiles, Conrad, Diamond and Taylor)
this arises modulo $p$ from a newform of weight $2$ and level $6$.
As there are no such newforms we arrive at a contradiction.
\end{proof}

Lemma~\ref{lem:modpow2} allows us to eliminate $13$ equations from Table~\ref{table:3}, all with large exponent $p$. Thus,
for Case 3 it remains to solve only $1$ equation with exponent $p=7$.

\subsection{Solving the final equations}
To complete the proof of Theorem~\ref{thm:main} it remains to solve $1+777=778$ equations of the form \eqref{eqn:tern}. Letting
$\sigma = w_2$ and $\tau = w_1^2$,
these can be considered as Thue equations,
\[
a\sigma^p - b \tau^p = c \, .
\]
where the exponent $p=5$ in all but one case where $p=7$.
Finally, we appeal to \texttt{Magma's} Thue solver \cite{magma}; details of the algorithms employed by this
solver can be found in \cite[Chapter VII]{Smart}. Solving these Thue equations we have found no non-trivial
solutions.
This completes the proof of Theorem~\ref{thm:main}.