\section{Introduction}
The Legendre curves play a fundamental role in $3$-dimensional contact geometry. Let $(M;\varphi,\xi,\eta, g)$ be an almost contact metric $3$-manifold. Then an integral curve of the contact distribution $\ker\eta=\{X\in\Gamma(\mathcal{T}M)\,\vert\,\eta(X)=0\}$ is known as Legendre curve; $\Gamma\left(\mathcal{T}M\right)$ being the section of tangent bundle $\mathcal{T}M$ of $M$ \cite{BD,DEB}. The authors in \cite{camci, MB} have studied these curves in $3$-dimensional Sasakian pseudo-metric manifolds. In the present study, we will focus on some specific Legendre curves in quasi-Sasakian pseudo-metric $3$-manifolds. \\
From \cite{SB}, it is known that a regular curve $\upsilon:I\to\mathbb{E}^3$ in a $3$-dimensional Euclidean space $\mathbb{E}^3$ is said to be a spherical curve if and only if
\begin{align}\label{de}
\frac{\tau(s)}{\kappa(s)}+\left(\frac{1}{\kappa(s)}\left(\frac{1}{\tau(s)}\right)'\right)'=0,
\end{align}
where curvature $\kappa(s)$ and torsion $\tau(s)$ of $\upsilon$ are non-zero smooth functions on open interval $I$, and the prime $'$ denotes differentiation with respect to arc length parameter $s$. The integral representation of Eq. \eqref{de} is given by
\begin{align}\label{ie}
\kappa(s)^{-1}=C\cos\left(\int_{s_0}^{s}\tau(t)\,dt\right)+D\sin\left(\int_{s_0}^{s}\tau(t)\,dt\right).
\end{align}
In this study, we present spherical characterization of Legendre curves in quasi-Sasakian pseudo metric $3$-manifolds. Moreover, we obtain that null Legendre curves in this class of manifolds are geodesic.

\section{Preliminaries}\label{pre}
Let $M$ be a $C^\infty$ and paracompact $(2n+1)$-manifold, then $M$ is said to be an almost contact manifold if its structure group ${\rm GL}_{2n+1}\mathbb{R}$ of tangent bundle $\mathcal{T}M$ is reducible to ${\rm U}(n)\times\{1\}$. This is equivalent to the existence of $(\varphi,\xi,\eta)$-structure satisfying
\begin{align}\label{2.1}
\varphi^2 = -\mathcal{I} +\eta\otimes\xi,\quad\eta\left(\xi\right)=1,
\end{align}
where $\varphi$ is an endomorphism, $\xi$ is a Reeb vector field, $\eta$ is a global $1$-form such that $\eta\wedge (d\eta)^n\ne 0$, $d$ being the exterior differential operator and $\mathcal{I}$ is the identity of $\mathcal{T}M$. We call $\eta$ the contact form. Eq. \eqref{2.1} yields $\varphi\xi = 0$, $\eta \circ\varphi =0$ and ${\rm rank}(\varphi)= 2n$. A pseudo-metric $g$ is said to be compatible with $(\varphi,\xi,\eta)$-structure if
\begin{align}\label{2.3}
g\left(\varphi\cdot, \varphi\cdot \right)=g(\cdot, \cdot) -\varepsilon \eta(\cdot)\eta(\cdot),
\end{align}
where $\varepsilon^2=1$ and signature of $g$ is either $(2p, 2n-2p +1)$ or $(2p+1, 2n-2p)$ according as Reeb vector field is  timelike or spacelike respectively. Such $(M;\varphi,\xi,\eta, g)$ is called an almost contact pseudo-metric $(2n+1)$-manifold. Here, $\eta (X)= \varepsilon g(X, \xi)$ and $g(\xi,\xi)=\varepsilon$. Thus $\xi$ is never lightlike. The fundamental $2$-form $\Phi$ is defined by
$\Phi (\cdot, \cdot) = \varepsilon g(\cdot, \varphi \cdot)$.\\

On the direct product manifold $M\times\mathbb{R}$ : $\left(X, \ell\frac{d}{dt}\right)$ is an arbitrary tangent vector, where  $\ell$ is a smooth function on $M\times\mathbb{R}$, $t$ being standard coordinate on $\mathbb{R}$  and $X\in\Gamma(\mathcal{T}M)$. The almost complex structure $J$ on this direct product is defined by
\begin{align*}J\left(X, \ell\frac{d}{dt}\right)=\left(\varphi X - \ell\xi, \eta(X)\frac{d}{dt}\right).\end{align*}
Then $M$ is said to be normal if and only if $J$ is integrable. Equivalently, $M$ is normal if and only if
\begin{align}\label{2.7}
N_{\varphi}(X, Y) + 2 d\eta(X, Y)\xi = 0,
\end{align}
where $N_{\varphi}$ denotes the Nijenhuis torsion of endomorphism $\varphi$ which is given as follows:
$$N_{\varphi}(X, Y) = [\varphi X, \varphi Y] - \varphi[\phi X, Y]+ \varphi^2[X, Y]  - \varphi[X, \varphi Y],$$
for any tangent vectors $X$, $Y$ on $M$ (cf, \cite{DEB, DP2019, GC2010, KS}).\\
In this paper, we confine to dimension three. Following \cite{ZO}, we give certain results relevant to this case. Let $M$ be an almost contact pseudo-metric $3$-manifold, then we receive that
\begin{align}\label{nablax}
(\nabla_{X}\varphi)Y = - \eta(Y)\varphi\nabla_{X}\xi+\varepsilon g(\varphi\nabla_{X}\xi, Y)\xi ,\ X\in\Gamma (\mathcal{T}M).
\end{align}
\begin{prop}\label{equiv}
In an almost contact pseudo-metric $3$-manifold $M$, we have the following mutually equivalent conditions:
\begin{itemize}
\item[(i)] $M$ is normal,
\item[(ii)] $\nabla_{\varphi X}\xi=\varphi\nabla_X \xi$,
\item[(iii)] $\nabla_{X}{\xi} = - \varepsilon\alpha \varphi X + \varepsilon\beta \left(X-\eta\left(X\right)\xi \right)$,\end{itemize}
where $X\in\Gamma (\mathcal{T}M)$, $\alpha$ and $\beta$ being smooth functions on $M$ for which we have that
\begin{align}\label{3.14}
2 \alpha = {\rm trace}\left\lbrace X \rightarrow \varphi\nabla_{X}\xi\right\rbrace, \,\,\, 2\beta = {\rm trace}\left\lbrace X \rightarrow\nabla_{X}\xi\right\rbrace.
\end{align}
\end{prop}
\noindent Using Eq.\eqref{nablax} and above proposition, we have
\begin{align}\label{3.1}
(\nabla_{X}\varphi)Y = \beta(g(\varphi X, Y)\xi - \varepsilon\eta(Y)\varphi X)+\alpha(g(X, Y)\xi - \varepsilon\eta (Y)X).
\end{align}
Let us denote by $\mathbb{N}^3$ a normal almost contact pseudo-metric $3$-manifold. Analogous to \cite{Inoguchi15}, $\mathbb{N}^3$ is called a quasi-Sasakian pseudo-metric $3$-manifold (denoted by $\mathbb{Q}^3$) if $\beta=0$ and $\xi(\alpha)=0$, in particular Sasakian if $\beta=0$ and $\alpha=1$. Below, we give an example of $\mathbb{N}^3$.
\begin{ex}
Let $(x,y,z)$ be the standard Cartesian coordinates on $\mathbb{R}^3$ and consider the $1$-form $\eta=2ydx+dz$. We put $\xi=\partial_3$ and consider the endomorphism $\varphi$ defined by
$\varphi\partial_1=\partial_2,\ \varphi\partial_2=2y\partial_3-\partial_1,\ \varphi\partial_3=0$, where $\partial_1=\frac{\partial}{\partial x}$, $\partial_2=\frac{\partial}{\partial y}$ and $\partial_3=\frac{\partial}{\partial z}$. Then it follows that $\eta(\xi)=1$ and $\varphi^2=-\mathcal{I}+\eta\otimes\xi$. Hence $(\varphi,\xi,\eta)$ is an almost contact structure on $\mathbb{R}^3$. By straightforward computations,we obtain that
\begin{align*}
N_{\varphi}(X, Y) \left(\partial_i ,\partial_j\right)+2d\eta\left(\partial_i ,\partial_j\right)\xi=0,\quad i,j\in\{1,2,3\}
\end{align*}
which implies that the structure is normal.

\noindent Let $\mathcal{N}^3_{\varepsilon}:=\mathbb{R}^3$ and consider a normal almost contact structure $(\varphi,\xi,\eta)$ restricted to $\mathcal{N}^3_{\varepsilon}$. Next, we consider the metric tensor $g=\exp(2z)(dx^2+dy^2)+\varepsilon\eta\otimes\eta$, with matrix representation with respect to the standard basis
 \begin{align*}
 g[(\partial_i, \partial_j)]=\begin{bmatrix}4\varepsilon y^2 +\exp(2z)&0& 2\varepsilon y \\ 0&\exp(2z)&0\\ 2\varepsilon y&0&\varepsilon \end{bmatrix},
 \end{align*}
where $\varepsilon^2=1$. Using the endomorphism $\varphi$ and the metric $g$, we obtain $g\left(\varphi X, \varphi Y\right)=g(X, Y) -\varepsilon \eta(X)\eta(Y)$ and $\eta(X)=\varepsilon g(X,\xi)$, and hence that $(\mathcal{N}^3_{\varepsilon};\varphi,\xi,\eta,g)$ is a normal almost contact pseudo-metric $3$-manifold. For the Levi-Civita connection $\nabla$ with respect to this metric, we have
\begin{align*}
\nabla_{\partial_1}\partial_1=&2y\partial_1 -4\varepsilon y\exp(-2z)\partial_2 -\varepsilon(4\varepsilon y^2+\exp(2z))\partial_3,\\
\nabla_{\partial_1}\partial_2 =& 2\varepsilon y\exp(-2z)\partial_1 +(-4\varepsilon y^2+\exp(2z))\exp(-2z)\partial_3=\nabla_{\partial_2}\partial_1, \\
\nabla_{\partial_1}\partial_3=& \partial_{1}-\varepsilon \exp(-2z)\partial_2 -2y\partial_3=\nabla_{\partial_3}\partial_1,\\
\nabla_{\partial_2}\partial_2=&2y\partial_1 -\varepsilon(4\varepsilon y^2+\exp(2z))\partial_3, \\
 \nabla_{\partial_2}\partial_3=&\varepsilon \exp(-2z)\partial_1+\partial_2 -2\varepsilon y\exp(-2z)\partial_3=\nabla_{\partial_3}\partial_2, \nabla_{\partial_3}\partial_3=0.\nonumber
\end{align*}
Using Eq. \eqref{3.1} and the above expressions, we find $\alpha=\exp(-2z)$ and $\beta=\varepsilon$. \end{ex}

\section{Legendre curves}\label{lc1}
Let $ M$ be an almost contact pseudo-metric $3$-manifold with Levi-Civita connection $\nabla$ and $\upsilon : I\longrightarrow M$ a unit speed curve parametrized by the arc-length $s$ in $M$,  $I$ being an open interval, then the Frenet frame (or Frenet $3$-frame) $\{T\coloneqq\upsilon',\ N,\ B\}$ of $\upsilon$ satisfies the following {\it Frenet-Serret} formulas:
\begin{align}
  \nabla_{T}T = \kappa N, \,\,\, \nabla_{T}N = -\kappa T + \varepsilon \tau B \,\,\, {\rm{and}} \,\,\, \nabla_{T}B = -\tau N,\label{4.1}
\end{align}
where $\kappa=\vert\nabla_{T}{T}\vert$  and $\tau$ are the geodesic curvature and geodesic torsion of $\upsilon$, respectively. The vector fields $T$, $B$ and $N$ are known as the tangent, binormal and principal normal vector fields of $\upsilon$, respectively.
If $\nabla_{\upsilon '}\upsilon ' = 0$ then $\upsilon$ is said to be a geodesic and non-geodesic if $\kappa>0$ everywhere on $I$.  \\
\noindent If $\upsilon$ is an integral curve of the contact distribution $\mathcal{D}={\rm Ker}\, \eta$, equivalently, $\eta(\upsilon ') = 0$ then we say that $\upsilon$ is a {\it Legendre curve} in $M$. These curves have been intensively studied by several authors (see \cite{BD, camci, CM14, CM15, Inoguchi, JW07, MB}).\\

\subsection{Legendre curves in $\mathbb{N}^3$}
Let $\upsilon:I\rightarrow \mathbb{N}^3$ be a Legendre curve in $\mathbb{N}^3$. Then $\{\upsilon',\varphi\upsilon',\xi\}$ is an orthonormal frame along $\upsilon$, where $g(\upsilon',\upsilon')=g(\varphi\upsilon',\varphi\upsilon')=1$ and $g(\xi,\xi)=\varepsilon$. Let us put $\nabla_{T}T=a_1\upsilon'+\vartheta\varphi\upsilon'+c_1\xi,$ where $a_1,\vartheta$ and $c_1$ are some smooth functions on $I$. Then we find that $a_1=g(\nabla_{\upsilon'}\upsilon', \upsilon')=0$, $\vartheta=g(\nabla_{\upsilon'}\upsilon', \varphi\upsilon')$ and $-\varepsilon c_1=g(\upsilon',\nabla_{\upsilon'}\xi)$. Employing Eqs. \eqref{2.3} and \eqref{3.1}, we have $g(\upsilon',\nabla_{\upsilon'}\xi)=\varepsilon\beta$. Thus,
\begin{align}\label{nablaTT}
\nabla_{T}T=-\beta \xi+\vartheta \varphi \upsilon',
\end{align}
which yields $\kappa=\sqrt{\vert\vartheta^2+\varepsilon\beta^2\vert}$. If $\kappa>0$ everywhere on $I$ then differentating $N=\frac{1}{\kappa}\nabla_T T=-\frac{\beta}{\kappa}\xi+\frac{\vartheta}{\kappa}\varphi \upsilon'$ along $\upsilon$ and using Eqs. \eqref{3.1} and \eqref{nablaTT}, we obtain that
 \begin{align}
 \nabla_{T}N=\frac{c}{\kappa}\left(\vartheta \xi+\epsilon \beta \varphi \upsilon'\right)-\kappa\upsilon', \label{nablaN}
 \end{align}
where $c=\alpha+\frac{\beta \vartheta'-\beta' \vartheta}{\epsilon\beta^2+\vartheta^2}$. In light of \eqref{4.1},  we have from Eq. \eqref{nablaN} that $\varepsilon\tau B=\frac{c}{\kappa}(\vartheta \xi+\epsilon \beta \varphi \gamma')$, which provides $\tau=\vert c\vert$. This leads to the following result which generalize Theorem 1 of \cite{JW07} to the case of spacelike or timelike Reeb vector field.
\begin{prop}\label{prop}
Let $\upsilon : I\longrightarrow \mathbb{N}^3$ be a non-geodesic Legendre curve in $\mathbb{N}^3$. Then its curvature $\left(\kappa\right)$ and torsion $\left(\tau\right)$ are given by
\begin{align}\label{4.2}
\kappa=&\sqrt{\vert\vartheta^2+\varepsilon\beta^2\vert},\\
\label{4.3} \tau=&\left\arrowvert\alpha +\left(\frac{\beta\vartheta'-\beta'\vartheta}{\kappa^2}\right)\right\arrowvert,
\end{align}
where $\vartheta$ is a smooth function defined by $\vartheta=g\left(\nabla_{\gamma'}\gamma', \varphi\gamma'\right)$
on $I$, and $\alpha$, $\beta$ being the same as in \eqref{3.14}.
\end{prop}
\begin{rem}
\begin{itemize}
\item[(a)]In $\mathbb{N}^3$, the decomposition of Reeb vector field $\xi$ in the Frenet frame of a non-geodesic Legendre curve is given by
\begin{align}
\xi=\frac{\epsilon}{\sqrt{\vert\vartheta^2+\varepsilon\beta^2\vert}}\left\{-\beta N+\vartheta B\right\}.\label{xi}
\end{align}
\item[(b)] Let $\upsilon$ be a non-geodesic Legendre curve in $\mathbb{N}^3$ with $\xi$ is parallel to binormal $B$. Then $\mathbb{N}^3$ is a quasi-Sasakian pseudo-metric $3$-manifold.
\item[(c)] Let $\upsilon$ be a non-geodesic Legendre curve in $\mathbb{N}^3$ such that $\xi$ is parallel to principal normal $N$ then it has torsion $\vert\alpha\vert$. In particular, the Legendre curves of a $\beta$-Kenmotsu manifold with $\xi$ is spacelike and parallel to $N$ are circles \cite{CMM}.
\end{itemize}
\end{rem}
\begin{ex}
\noindent If we define a curve $\upsilon:I\subset\mathbb{R}\rightarrow\mathcal{N}^3_{\varepsilon}$ by $\upsilon(s)=\left(\upsilon_1(s),\upsilon_2(s),\upsilon_3(s)\right)$. Then $\upsilon$ is a Legendre curve if and only if
\begin{align}\label{con}
\begin{cases}
2\upsilon_1'\upsilon_2+\upsilon_3'=0,\\
\upsilon_1'^2+\upsilon_2'^2=\exp(-2\upsilon_3).\end{cases}
\end{align}
It can be easily seen that
\begin{align*}
\upsilon'=\upsilon_1'\partial_1+\upsilon_2'\partial_2-2\upsilon_1'\upsilon_2\partial_3\ {\rm and}\
\varphi\upsilon'=-\upsilon_2'\partial_1+\upsilon_1'\partial_2+2\upsilon_2\upsilon_2'\partial_3.
\end{align*}
Below, we present certain Legendre curves in $\mathcal{N}^3_{\varepsilon}$.
\begin{itemize}
\item[(i)] Let us consider a curve
\begin{align*}
\upsilon^1(s)=(1,s,0)
\end{align*}
in $\mathcal{N}^3_{\varepsilon}$, where $s\in I\subseteq\mathbb{R}$ is arc length parameter of $\upsilon^1$. Then this curve satisfies Eq. \eqref{con}, so it is a Legendre curve in $\mathcal{N}^3$. For such a curve $\alpha(\upsilon(s))=1$, $\beta(\upsilon(s))=\varepsilon$, $\vartheta=g\left(\nabla_{\upsilon'}\upsilon',\varphi\upsilon'\right)=-2s$, $$\kappa=\sqrt{\vert1+4\varepsilon s^2\vert} \ {\rm and}\ \tau=\left\vert1-\frac{2\varepsilon}{\kappa^2}\right\vert.$$
\item[(ii)] Let us consider a curve
\begin{align*}
\upsilon^2(s)=(-\ln s,1/2,\ln s)
\end{align*}
in $\mathcal{N}^3_{\varepsilon}$, where $s>0$ is arc length parameter of $\upsilon^2$. Then this curve satisfies Eq. \eqref{con}, so it is a Legendre curve in $\mathcal{N}^3$. For such a curve $\alpha(\upsilon(s))=1/{s^2}$, $\beta(\upsilon(s))=\varepsilon$, $\vartheta=0$, $$\kappa=1 \ {\rm and}\ \tau=\frac{1}{s^2}.$$
\end{itemize}
Further, the spacelike and timelike situations of these Legendre curves are interesting. The Euclidean pictures of $\upsilon^1$ and $\upsilon^2$ are rendered in Fig. $1$ and Fig. $2$ respectively.
\begin{figure}[ht]
\begin{minipage}{0.3\textwidth}
\includegraphics[height=5 cm]{1.pdf}\\
\center(a) $\upsilon^1$ in $\mathcal{N}^3_{1}$
\end{minipage}
\hspace*{2 cm}
\begin{minipage}{0.3\textwidth}
\includegraphics[height=5 cm]{2.pdf}\\
\center(b) $\upsilon^1$ in $\mathcal{N}^3_{-1}$
\end{minipage}
\caption{Legendre curve $\upsilon^1$ in $\mathcal{N}^3_{\varepsilon}$}
\end{figure}

\begin{figure}[ht]
\begin{minipage}{0.3\textwidth}
\includegraphics[height=5 cm]{33.pdf}\\
\center(c) $\upsilon^2$ in $\mathcal{N}^3_{1}$
\end{minipage}
\hspace*{2 cm}
\begin{minipage}{0.3\textwidth}
\includegraphics[height=5 cm]{4.pdf}\\
\center(d) $\upsilon^2$ in $\mathcal{N}^3_{-1}$
\end{minipage}
\caption{Legendre curve $\upsilon^2$ in $\mathcal{N}^3_{\varepsilon}$}
\end{figure}
\end{ex}

\noindent Let $\upsilon:I\rightarrow \mathbb{N}^3$ be a unit speed curve such that $\eta(\upsilon')=m$, where $m$ is a smooth function on $I$. Let us consider the vector fields $\upsilon', \varphi\upsilon'$ and $\xi$ for which, we have that $g(\upsilon', \upsilon')=1$, $g(\upsilon', \xi)=\varepsilon m$, $g(\varphi\upsilon', \varphi\upsilon')=1-\varepsilon m^2$, $g(\xi, \xi)=\varepsilon$ and $g(\upsilon',\varphi\upsilon')=g(\xi,\varphi\upsilon')=0$. Then, $\upsilon', \varphi\upsilon', \xi$ are linearly independent. Further, $\{\upsilon', \varphi\upsilon', \xi\}$ forms a basis of $\mathcal{T}_{\upsilon(s)}M^3$ for any $s\in I$ if and only if $1-\varepsilon m^2\ne 0.$ In this case, one can define orthonormal frame $\{V_1,V_2,V_3\}$ along $\upsilon$ as follows:
\begin{align}\label{frame}
V_1=\upsilon',\quad\quad V_2=\frac{\varphi\upsilon'}{\delta},\quad\quad V_3=\frac{\xi-\varepsilon m\upsilon'}{\delta}
\end{align}
where $g(V_1, V_1)=g(V_2, V_2)=1$, $g(V_3, V_3)=\varepsilon$ and $\delta=\sqrt{\vert1-\varepsilon m^2\vert}$.

\noindent It is important to note that $\{\upsilon', \varphi\upsilon', \xi\}$ is linearly dependent if and only if either $\upsilon'=m\xi$ or $\upsilon'=m\xi+\varepsilon\varphi\upsilon'$. In reality, if these vector fields are linearly dependent then $\xi$ is spacelike and $\vert m\vert=1$. In this case, the curve $\upsilon$ is necessarily a geodesic. Therefore, for the non-geodesic curve $\upsilon$, we must have that neither $\varepsilon= 1$ nor $\vert m\vert=1$, that is $\delta\ne 0$. The decomposition of characteristic vector field with respect to the frame $\left\{V_1,V_2,V_3\right\}$ is given by
\begin{align}
\xi=\varepsilon\left(mV_1+\delta V_3\right).\label{dxi}
\end{align}
 Below we present some preparatory lemmas for later use.
\begin{lem}
Let $\upsilon:I\rightarrow \mathbb{N}^3$ be a non-geodesic unit speed curve in $\mathbb{N}^3$ such that $\eta(\upsilon')=m$. Then the covariant derivative of $V_1, V_2, V_3$ along $\upsilon$ are given by:
\begin{align}
\nabla_{\upsilon'}V_1=&\delta\vartheta_1 V_2-\left(\beta\delta-\frac{m'}{\delta}\right)V_3,\label{E1}\\
\nabla_{\upsilon'}V_2=&-\delta\vartheta_1V_1+(\alpha+m\vartheta_1)V_3,\label{E2}\\
\nabla_{\upsilon'}V_3=&\varepsilon\left( \left(\beta\delta-\frac{m'}{\delta}\right) V_1-(\alpha+m\vartheta_1)V_2\right),\label{E3}
\end{align}
where
\begin{align}\label{abcm}
\delta=\sqrt{\vert1-\varepsilon m^2\vert} \ {\rm and}\ \vartheta_1=\frac{1}{\delta^2}g(\nabla_{\upsilon'}\upsilon',\varphi\upsilon').
\end{align}
\end{lem}
\begin{proof}
Let us express $\nabla_{\upsilon'}V_1=a_2V_1+b_2V_2+c_2V_3$, where $a_2,\ b_2, \ c_2$ are some smooth functions on $I$. Then $a_2=g(\nabla_{\upsilon'}V_1, V_1)=0$.

\noindent Using Eq. \eqref{3.1}, we find
\begin{align*}
-b_2&=g(V_1, \nabla_{\upsilon'}V_2)=\frac{1}{\delta}g((V_1,(\nabla_{\upsilon'}\varphi)\upsilon'+\varphi\nabla_{\upsilon'}\upsilon')\\
&=-\frac{1}{\delta}g(\varphi V_1, \nabla_{\upsilon'}\upsilon'))=-\delta\vartheta_1.
\end{align*}

\noindent Employing Proposition \eqref{equiv}, we have
\begin{align*}
-\varepsilon c_2&=g(V_1,\nabla_{\upsilon'}V_3)=\frac{1}{\delta}g(V_1,\nabla_{\upsilon'}\xi-\varepsilon m'\upsilon')\\
&=\frac{1}{\delta}g(V_1,\varepsilon\beta(\upsilon'-m\xi)-\varepsilon m'\upsilon')
=\frac{\varepsilon}{\delta}(\beta \delta^2-m').
\end{align*}
In light of above expressions, we obtain that $a_2=0$, $b_2=\delta\vartheta_1$, $c_2=-\frac{\varepsilon}{\delta}(\beta \delta^2-m')$ and this leads to Eq. \eqref{E1}. Similarly, we compute $\nabla_{\upsilon'}V_2$ and $\nabla_{\upsilon'}V_3$.
\end{proof}

\begin{lem}
Let $\upsilon:I\rightarrow \mathbb{N}^3$ be a non-geodesic unit speed curve in $\mathbb{N}^3$ such that $\eta(\upsilon')=m$. Then its curvature and torsion are given by
 \begin{align}\label{kt1}
  \begin{cases}
 \kappa=&\delta\sqrt{\left\vert\vartheta_1^2+\varepsilon\left(\beta-\frac{m'}{\delta^2}\right)^2\right\vert},\\
\tau=&\left\vert \alpha+m\vartheta_1+\frac{\left(\beta\vartheta_1'-\beta'\vartheta_1\right)-2\left(\frac{m'\vartheta_1'}{\delta^2}\right)+\left(\frac{m'\vartheta_1}{\delta^2}\right)'}{\vartheta_1^2+\varepsilon\left(\beta -\frac{m'}{\delta^2}\right)^2}\right\vert,
 \end{cases}
 \end{align}
 where $\delta$ and $\vartheta_1$ are the same as in Eq. \eqref{abcm}.
\end{lem}
\begin{proof}
\noindent Using Eq. \eqref{E1} and computing the length of $\nabla_{\upsilon'}\upsilon'$, we receive the curvature function $\kappa$. By the virtue of Eqs. \eqref{4.1} and \eqref{E1}, we find
\begin{align}\label{N}
N=\frac{\delta\vartheta_1}{\kappa}V_2-\left(\frac{\beta\delta}{\kappa}-\frac{m'}{\kappa\delta}\right)V_3.
\end{align}

\noindent Let us denote $a_2=\delta\vartheta_1$ and $b_2=\left(\beta\delta-\frac{m'}{\delta}\right)$. Now, we compute
\begin{align}
\nabla_{\upsilon'}N=&\left(\frac{a_2}{\kappa}\right)'V_2+\left(\frac{a_2}{\kappa}\right)\nabla_{\upsilon'}V_2-\left(\frac{b_2}{\kappa}\right)'V_3-\left(\frac{b_2}{\kappa}\right)\nabla_{\upsilon'}V_3\nonumber\\
=&\left(-\frac{{a_2}^2}{\kappa}-\frac{\varepsilon {b_2}^2}{\kappa}\right)V_1+\left(\frac{{a_2}'\kappa-a_2{\kappa}'}{\kappa^2}+\frac{\varepsilon b_2}{\kappa}\left(\alpha+\frac{a_2 m}{\delta}\right)\right)V_2 \nonumber\\
& +\left(-\frac{{b_2}'\kappa-b_2{\kappa}'}{\kappa^2}+\frac{a_2}{\kappa}\left(\alpha+\frac{a_2 m}{\delta}\right)\right)V_3\nonumber\\
=&-\kappa V_1+\frac{\varepsilon b_2}{\kappa}\left(\alpha+\frac{a_2 m}{\delta}+\frac{{a_2}'b_2-a_2{b_2}'}{\kappa^2}\right)V_2+\frac{a_2}{\kappa}\left(\alpha+\frac{a_2 m}{\delta}+\frac{{a_2}'b_2-a_2{b_2}'}{\kappa^2}\right)V_3\nonumber\\
=&-\kappa V_1+\left(\alpha+\frac{a_2 m}{\delta}+\frac{{a_2}'b_2-a_2{b_2}'}{\kappa^2}\right)\left(\frac{\varepsilon b_2V_2+a_2V_3}{\kappa}\right).\label{t}
\end{align}
In light of the Frenet-Serret formulas, we have from Eq. \eqref{t} that $$\tau=\left\vert \alpha+\frac{a_2 m}{\delta}+\frac{{a_2}'b_2-a_2{b_2}'}{\kappa^2}\right\vert.$$  Using the values of $a_2$ and $b_2$ in the above expression, we find the torsion function of $\upsilon$ as mentioned in Eq. \eqref{kt1}. This completes the proof.
\end{proof}
\begin{thm}
Let $\upsilon:I\rightarrow \mathbb{N}^3$ be a non-geodesic unit speed curve in $\mathbb{N}^3$ such that $\eta(\upsilon')=m$. Then, the decomposition of Reeb vector field $\xi$ in the Frenet frame $\{T:=\upsilon',N,B\}$ of $\upsilon$ is given by
\begin{align}\label{decomp}
\xi={\varepsilon}(m T+\eta(N)N)+\eta(B)B,
\end{align}
 where $\eta(N)=(m'-\beta\delta^2)/\kappa$, $\eta(B)=(\varepsilon\mathop{{sgn}}(\tau)\delta^2\vartheta_1)/\kappa$, $\delta$ and $\vartheta_1$ are the same as in Eq. \eqref{abcm}.
\end{thm}
\begin{proof}\ Let us put $\xi=pT+qN+rB$, where $p,q$ and $r$ are smooth functions on $I$. Then $p=g(\xi,T)=\varepsilon m$, $q=g(\xi,N)=\varepsilon\eta(N)$ and $r=\varepsilon g(\xi,B)=\eta(B)$. The values of $\eta(N)$ and $\eta(B)$ directly follows from Eqs. \eqref{4.1}, \eqref{dxi}, \eqref{N} and \eqref{t}. This completes the proof.

\end{proof}
\noindent In view of Eq. \eqref{decomp}, we find that $\eta(B)^2+\varepsilon\eta(N)^2=\delta^2$. This leads to the following remark.
\begin{rem}
For a non-geodesic unit speed curve $\upsilon$ in $\mathbb{N}^3$ such that $\eta(\upsilon')=m$, we have
\begin{align*}
\delta^2-\varepsilon\eta(N)^2\ge 0.
\end{align*}
\end{rem}

\subsection{Legendre curves in $\mathbb{Q}^3$}
Let $\upsilon$ be a non-geodesic Legendre curve in a quasi-Sasakian pseudo-metric $3$-manifold $\mathbb{Q}^3$. Then by the virtue of proposition \ref{prop}, curvature and torsion of $\upsilon$ are given by
\begin{align}\label{quasikt}
\kappa=\vartheta,\ \tau=\vert\alpha\vert,
\end{align}
where $\vartheta=g\left(\nabla_{\upsilon'}\upsilon', \varphi\upsilon'\right)>0$ and $\alpha$ being the same as in \eqref{3.14}.
\begin{thm} \label{thm Br}
For a non-geodesic unit speed curve $\upsilon$ in $\mathbb{Q}^3$ with $\eta(\upsilon')=m$. If $ \tau=\vert\alpha\vert$ and at one point $m=m'=m'' =0$, where $\alpha$ is non-zero at every point of $M$. Then $\upsilon$ is a Legendre curve.
\end{thm}
\begin{proof}
For spacelike vector $\xi$, this is Theorem 3 of \cite{JW07}, so we concentrate on the timelike situation.

\noindent In this class of manifold, we have from Eq. \eqref{kt1} that
\begin{align}\label{tauquasi}
\tau=&\left\vert \alpha+m\vartheta_1+\left(\left(\frac{m'\vartheta_1}{\delta^2}\right)'-2\left(\frac{m'\vartheta_1'}{\delta^2}\right)\right)\left(\vartheta_1^2-\frac{m'^2}{\delta^4}\right)^{-1}\right\vert,
\end{align}
where $\delta=\sqrt{1+ m^2}\ge 1$ and $\vartheta_1\delta^2=g(\nabla_{\upsilon'}\upsilon',\varphi\upsilon')$. Since $m(s_0)=m'(s_0)=m''(s_0)=0$ at certain $s_0\in I$, therefore $\tau(s_0)=\vert\alpha(s_0)\vert\ne 0$. Further, $\vartheta_1(s_0)\ne 0$ since $\upsilon$ is non-geodesic. Taking these into consideration and using $\tau=\vert\alpha\vert$, Eq. \eqref{tauquasi} yields
\begin{align}\label{vartheta1m}
\left(m''\vartheta_1-m'\vartheta_1'\right)-3mm'\delta^{-2}\vartheta_1+m\delta^2\vartheta_1^3=0.
\end{align}
Obviously $m=0$ is a solution and $m'=0$ implies that $m=0$. Hence we now assume that $m$ is not constant. Setting $t=m'/\vartheta_1$, Eq. \eqref{vartheta1m} becomes
$\frac{d}{dm}\left(t^2\right)-\frac{6m}{\delta^2}t^2+2m\delta^2=0$.
Integration gives ${m'}^2=(1/3+C\delta^3)\delta^3\vartheta_1^2$, where $C$ is constant of integration. Using initial conditions, we have $C=-1/3$. Thus ${m'}^2=-(\delta^3-1)\delta^3\vartheta_1^2/3$. Since $\delta\ge 1$ and $\vartheta_1\ne 0$, we have $m=0$, a contradiction. This completes the proof.
\end{proof}

\begin{ex}
Let $\mathcal{Q}^3_{\alpha}:=\mathbb{R}_{+}\times\mathbb{R}^2$ with standard Cartesian coordinates $(x,y,z)$, $\eta=-2xdy+dz$, $\xi=\partial_3$, $\varphi$ defined by: $\varphi\partial_1=\partial_2+2x\partial_{3},\ \varphi\partial_2=-\partial_1,\ \varphi\partial_3=0$ and metric tensor $g=x^{2}(dx^2+dy^2)+\varepsilon\eta\otimes\eta$, where $\partial_1=\frac{\partial}{\partial x}$, $\partial_2=\frac{\partial}{\partial y}$ and $\partial_3=\frac{\partial}{\partial z}$. Then $(\mathcal{Q}^3_{\alpha};\varphi,\xi,\eta,g)$  is a normal almost contact pseudo-metric $3$-manifold.  For the Levi-Civita connection $\nabla$ with respect to this metric, we find
\begin{align}\label{levi}
&\nabla_{\partial_1}\partial_1=\frac{1}{x}\partial_1, \nabla_{\partial_2}\partial_2=-\frac{(1+4\varepsilon)}{x}\partial_1, \nabla_{\partial_3}\partial_3=0, \nabla_{\partial_2}\partial_3=\nabla_{\partial_3}\partial_2=\frac{\varepsilon}{x^{2}}\partial_1, \nonumber
\\
&\nabla_{\partial_1}\partial_2=\nabla_{\partial_2}\partial_1=\frac{(2\varepsilon+1)}{x}\partial_2+(1+4\varepsilon)\partial_3, \nabla_{\partial_1}\partial_3=\nabla_{\partial_3}\partial_1=-\frac{\varepsilon}{x^{2}}\partial_2-\frac{2\varepsilon}{x}\partial_3.
\end{align}
In view of Eqs. \eqref{3.1} and \eqref{levi}, we have $\alpha=1/x^{2}$ and $\beta=0$. Thus $\mathcal{Q}^3_{\alpha}$ is a quasi-Sasakian $3$-manifold.
\end{ex}
\noindent Using definition of Legendre curve, we find that $\upsilon:I\rightarrow\mathcal{Q}^3_{\alpha}$ defined by $\upsilon(s)=\left(\upsilon_1(s),\upsilon_2(s),\upsilon_3(s)\right)$ is a Legendre curve if and only if
\begin{align}\label{k1}
\begin{cases}
\upsilon_3'=2\upsilon_1\upsilon'_2,\\
\upsilon_1'^2+\upsilon_2'^2=\upsilon_1^{-2}.\end{cases}
\end{align}
In this case, we have
\begin{align*}
\upsilon'=\upsilon_1'\partial_1+\upsilon_2'\partial_2+2\upsilon_1\upsilon'_2\partial_3\ {\rm and}\
\varphi\upsilon'=-\upsilon_2'\partial_1+\upsilon_1'\partial_2+2\upsilon_1\upsilon'_1\partial_3.
\end{align*}
From Eq. \eqref{k1}, we have
\begin{align*}
\upsilon_1'=\upsilon_1^{-1}\cos\phi(s), \upsilon_2'=\upsilon_1^{-1}\sin\phi(s), \upsilon_3'=2\sin\phi(s),
\end{align*}
where $\phi\in C^{\infty}(I)$. Thus, we have the general expression of Legendre curves in $\mathcal{Q}^3_{\alpha}$ as follows:
\begin{prop}\label{lq3}
Let $\upsilon:I\rightarrow\mathcal{Q}^3_{\alpha}$ be a non-geodesic Legendre curve, then $\upsilon$ is given by
\begin{align}\label{k2}
\upsilon(s)=\left(\int_{s_0}^{s}{\mu^{-1}}(t)\zeta(t)dt, 2\int_{s_0}^{s}\sin\phi(s) \right),
\end{align}
where $\phi$ being smooth function on $I$, $\zeta(s)=(\cos\phi(s), \sin\phi(s))$ is any parametrization of circle $\mathbb{S}^1$ and $\mu$ is a non-zero smooth function defined by
$\mu^2(s)=2\int_{s_0}^{s}\cos\phi(t)dt$.
\end{prop}
\begin{thm}
The curvature and torsion of non-geodesic Legendre curve \eqref{k2} in $\mathcal{Q}^3_\alpha$ are given by
\begin{align}\label{kt3}
\kappa=\phi'+\frac{\sin\phi}{\mu^2},
\tau=\frac{1}{\mu^2}.
 \end{align}
\end{thm}

\begin{proof}
 Using Eqs. \eqref{levi} and \eqref{k2}, we have
\begin{align}
\nabla_{\upsilon'}\upsilon'=\left(\phi'+\frac{\sin\phi}{\mu^2}\right)\left(-\frac{\sin\phi}{\mu}\partial_1+\frac{\cos\phi}{\mu}\partial_2+2(\cos\phi)\partial_3\right).
\end{align}
Employing Frenet formula, the above equation yields $\kappa=\phi'+\frac{\sin\phi}{\mu^2}$. If $\kappa>0$ everywhere on $I$ then differentiating normal vector field $N$ along $\upsilon$, we get
\begin{align}\label{quasin}
\nabla_{\upsilon'}N&=-\left(\phi'+\frac{\sin\phi}{\mu^2}\right)\left(\frac{\cos\phi}{\mu}\partial_1+\frac{\sin\phi}{\mu}\partial_2\right)+\left(\frac{\cos2\phi}{\mu^2}-2\phi'\sin\phi\right)\partial_3\nonumber\\
&=-\kappa\upsilon'+\frac{1}{\mu^2}\partial_3.
\end{align}
By the virtue of \eqref{4.1}, we have from Eq. \eqref{quasin} that $\varepsilon\tau B=\frac{1}{\mu^2}\partial_3$, which provides $\tau$. This completes the proof.
\end{proof}

\section{Null Legendre curves in $\mathbb{Q}^3$}
If we consider the Reeb vector field to be timelike, i.e. $\varepsilon=-1$ in $\mathbb{Q}^3$ then the signature of the associated metric \eqref{2.3} is $(2,1)$. In this situation we have lightlike (null) vector fields in $\mathbb{Q}^3$. Here we give characterization of null Legendre curves in $\mathbb{Q}^3$.

\noindent An arbitrary curve $\upsilon:I\to \mathbb{Q}^3$  is said to be a null curve if
\begin{align}\label{nc}
g(\upsilon'(s),\upsilon'(s))=0\ {\rm and}\ \upsilon'(s)\ne 0\ \forall s\in I.
\end{align}
The general Frenet frame $\{T:=\upsilon', U, V\}$ along null curve $\upsilon$ is determined by
\begin{align}\label{nf}
\begin{cases}
g(U,U)=1,\ g(T,T)=g(V,V)=0,\\ g(T,V)=-1,\ g(T,U)=g(U,V)=0,
\end{cases}
\end{align}
where $T\times U=-T,\, U\times V=-V,\,V\times T=U$. This frame is positively oriented if $det(T,U,V)=1$. The frame  $\{T:=\upsilon', U, V\}$ satisfying the following general Frenet equations:
\begin{align}\label{nf1}
  \nabla_{T}T =\tilde{h}T+\tilde{ \kappa}_1 U, \,\,\, \nabla_{T}U = -\tilde{\tau}_1  T+\tilde{ \kappa}_1 V \,\,\, {\rm{and}} \,\,\, \nabla_{T}V = -\tilde{h}V-\tilde{\tau}_1 U,
\end{align}
where $\tilde{h},\tilde{ \kappa}_1$ and $\tilde{ \tau}_1$ are smooth functions on $I$. Furthermore, $\tilde{ \kappa}_1$ and $\tilde{ \tau}_1$ are known as curvature and torsion of null curve $\upsilon$, respectively. Since frame $\{T:=\upsilon', U, V\}$ and Frenet formulas \eqref{nf1} depend on the choice of screen distribution and parametrization of $\upsilon$. Therefore, they are not unique in general (for more details see \cite{KLD}).

\begin{prop}
Let $\upsilon:I\to \mathbb{Q}^3$ be a null Legendre curve in $\mathbb{Q}^3$. Then $\upsilon$ is a geodesic.
\end{prop}
\begin{proof}
From Eq. \eqref{nf}, we can express
\begin{align}\label{express}
\nabla_{\upsilon'}\upsilon'= a_3\upsilon'+b_3 U,
\end{align}
where $a_3$ and $b_3$ are certain smooth functions on $I$. Using Eqs. \eqref{3.1} and \eqref{express} in the differentiation of $g(\xi, \upsilon')=0$ along $\upsilon$, we find that $b_{3}g(\xi, U)=0$. For $b_{3}\neq 0$ we obtain $g(\xi, U)=0$. Then Eq. \eqref{nf} clearly indicates that $\xi=a_{4}\upsilon'+b_{4}V$ along $\upsilon$ for some smooth functions $a_4$ and $b_4$ on $I$. Now using above expression of $\xi$ and Eq. \eqref{nf}, we get $g(\xi, \xi)=-2 a_{4} b_{4}=-1$ which is not possible due to the fact that  $g(\xi, \upsilon')=-b_{4}=0$. Hence, we must have $b_{3}=0$ and therefore \eqref{express} gives $\nabla_{\upsilon'}\upsilon'= a_3\upsilon'$, this implies that $\upsilon$ is a geodesic after a reparametrization. This completes the proof.
\end{proof}

\noindent Below we give the spherical characterization of Legendre curve in a quasi-Sasakian pseudo-metric $3$-manifold $\mathbb{Q}^3$.
\section{Spherical Legendre curves in $\mathbb{Q}^3$}
Let $\mathbb{Q}^3$ be a quasi-Sasakian pseudo-metric $3$-manifold for which $\alpha\ne 0$. Then the sphere $\mathbb{S}_{\pm}^2$ in $\mathbb{Q}^3$ is defined by $\mathbb{S}_{\pm}^2:=\mathbb{S}_{\circ}^2\cup\mathbb{H}_{\circ}^2$, where $$\mathbb{S}_{\circ}^2=\{P\in M^3\,\vert\, g(P,P)=r^2\}$$ and $$\mathbb{H}_{\circ}^2=\{P\in M^3\,\vert \,g(P,P)=-r^2\}.$$
\begin{prop}\label{prop1}
Let $\upsilon:I\to \mathbb{Q}^3$ be a non-geodesic Legendre curve in $\mathbb{Q}^3$ such that $\alpha\ne 0$. Then the center of osculating sphere of $\upsilon$ at the point $\upsilon(s)$ is given by
\begin{align}\label{cent}
c(s)=\upsilon(s)+\frac{1}{\vartheta}N-\frac{\vartheta'}{\vartheta^2\vert\alpha\vert}B.
\end{align}
\end{prop}
\begin{proof}
Let us consider the sphere $\mathbb{S}_{\pm}^2$ passing through four neighbouring points of the Legendre curve $\upsilon$. Let $c(s)$ be the center of sphere and $r$ be its radius. Then the function $\gamma:I\to\mathbb{R}$ defined by $$\gamma(s)=g(c(s)-\upsilon(s),c(s)-\upsilon(s))\mp r^2$$
satisfies the following conditions:
\begin{align*}
\gamma(s)=\gamma'(s)=\gamma''(s)=\gamma'''(s)=0.
\end{align*}
In light of Eqs. \eqref{4.1} and \eqref{quasikt}, the above equation leads to the following relations:
\begin{align}\label{os}
\begin{cases}
g(c(s)-\upsilon(s),c(s)-\upsilon(s))=\pm r^2,\\
g(c(s)-\upsilon(s),T)=0,\\
g(c(s)-\upsilon(s),N)=\frac{1}{\vartheta},\\
g(c(s)-\upsilon(s),B)=-\frac{\varepsilon\vartheta'}{\vartheta^2\vert\alpha\vert}.
\end{cases}
\end{align}
Since $c(s)-\upsilon(s)\in {\rm Span}\{T,N,B\}$. Therefore, in view of Eq. \eqref{os}, we obtain \eqref{cent}.
\end{proof}
Let the sphere $\mathbb{S}_{\pm}^2$ in $\mathbb{Q}^3$ be centered at origin and Legendre curve $\upsilon$ lie on it. Then in view of above proposition,  osculating sphere of $\upsilon$ is $\mathbb{S}_{\pm}^2$, $\forall s\in I$.
\begin{prop}\label{prop2}
Let $\upsilon:I\to \mathbb{Q}^3$ be a non-geodesic Legendre curve in $\mathbb{Q}^3$ such that $\alpha\ne 0$. Let $\vartheta(s)$ be non-constant if Reeb vector field $\xi$ is spacelike and $\vartheta(s)\ne \vartheta(s_0)\exp\left(\pm\int_{s_0}^{s}\vert\alpha(t)\vert dt\right)$ if $\xi$ is timelike. Then the radius of $\mathbb{S}_{\pm}^2$ is constant, $\forall s\in I$ if and only if the centers of osculating spheres are the same constants.
\end{prop}
\begin{proof}
Let us assume that the radius of osculating sphere is constant. Then from Eq. \eqref{cent}, we find
\begin{align}\label{os1}
\left(\frac{1}{\vartheta}\right)^2+\varepsilon\left(\frac{\vartheta'}{\vartheta^2\vert\alpha\vert}\right)^2=\pm r^2.
\end{align}
 Differentiating above expression along $\upsilon$, we have
\begin{align*}
\frac{\varepsilon \vartheta'}{\vartheta^{2}\vert\alpha\vert}\left\{\left(\frac{\vartheta'}{\vartheta^{2}\vert\alpha\vert}\right)'-\frac{\varepsilon\vert\alpha\vert}{\vartheta}\right\}=0.
\end{align*}
Consequently, if $\vartheta'\neq 0$ then we get
\begin{align}\label{os2}
\left(\frac{\vartheta'}{\vartheta^{2}\vert\alpha\vert}\right)'-\frac{\varepsilon\vert\alpha\vert}{\vartheta}=0.
\end{align}
Further, differentiation of \eqref{cent} along $\upsilon$ gives
\begin{align}\label{os3}
\nabla_{\upsilon'}c(s)=\left\{\left(\frac{\vartheta'}{\vartheta^{2}\vert\alpha\vert}\right)'-\frac{\varepsilon\vert\alpha\vert}{\vartheta}\right\}B.
\end{align}
In view of Eqs.\eqref{os2} and \eqref{os3}, we find that $c(s)=$ constant for all $s \in I$.\\
Conversely, let $\nabla_{\upsilon'}c(s)=0\  \forall s\in I$ then using \eqref{os3} we obtain
\begin{align}\label{os4}
\left(\frac{\vartheta'}{\vartheta^{2}\vert\alpha\vert}\right)'-\frac{\varepsilon\vert\alpha\vert}{\vartheta}=0.
\end{align}
Now differentiating \eqref{os1} along $\upsilon$, we have
\begin{align*}
\frac{\varepsilon \vartheta'}{\vartheta^{2}\vert\alpha\vert}\left\{\left(\frac{\vartheta'}{\vartheta^{2}\vert\alpha\vert}\right)'-\frac{\varepsilon\vert\alpha\vert}{\vartheta}\right\}=\pm r(s)\left(\nabla_{\upsilon'}r(s)\right).
\end{align*}
Employing \eqref{os4} in the above equation, we get
\begin{align*}
 r(s)\left(\nabla_{\upsilon'}r(s)\right)=0
\end{align*}
which implies that $r(s)=$ constant for all $s\in I$. This completes the proof.
\end{proof}
\noindent By the consequences of Propositions \eqref{prop1} and \eqref{prop2}, we have the following result:
\begin{thm}
Let $\upsilon:I\to \mathbb{Q}^3$ be a non-geodesic Legendre curve in $\mathbb{Q}^3$ such that $\alpha\ne 0$. Suppose $\vartheta(s)$ is non-constant if $\xi$ is spacelike and $\vartheta(s)\ne \vartheta(s_0)\exp\left(\pm\int_{s_0}^{s}\vert\alpha(t)\vert dt\right)$ if $\xi$ is timelike. Then the following statements are mutually equivalent:
\begin{itemize}
\item[(a)] $\upsilon\subset \mathbb{Q}^3$ is a spherical curve.
\item[(b)] the centers of osculating spheres $\forall s\in I$ at the point $\upsilon(s)$ are the same constants.
\item[(c)] $\left(\frac{\vartheta'}{\vartheta^{2}\vert\alpha\vert}\right)'-\frac{\varepsilon\vert\alpha\vert}{\vartheta}=0$.
\end{itemize}
\end{thm}
\subsection{Integration of Eq. \eqref{os2}}
\noindent Eq. \eqref{os2} can be rewritten as
\begin{align}\label{os5}
\left(\left(\frac{1}{\vartheta}\right)'\frac{1}{\vert\alpha\vert}\right)'+\varepsilon\left(\frac{1}{\vartheta}\right)\vert\alpha\vert=0.
\end{align}
If we consider $\frac{1}{\vert\alpha\vert}=x(s)$ and $\frac{1}{\vartheta}=y(s)$, then above equation takes the form
\begin{align}\label{os6}
\left(x(s)y'(s)\right)'+\varepsilon\left(\frac{y(s)}{x(s)}\right)=0.
\end{align}
Changing the variable in above expression by substituting $\int_{s_0}^{s}\frac{1}{x(t)}dt=z(s)$, then we have
 \begin{align}\label{os7}
\frac{1}{x(s)}\frac{d^2y}{dz^2}=\left(x(s)y'(s)\right)'.
\end{align}
From Eqs. \eqref{os6} and \eqref{os7}, we find
 \begin{align}\label{os8}
\frac{d^2y}{dz^2}+\varepsilon y=0.
\end{align}
If $\xi$ is spacelike then general solution of differential equation \eqref{os8} is given by
\begin{align}\label{os9}
y(s)=A_1\cos\left(\int_{s_0}^{s}\vert\alpha(t)\vert dt\right)+A_2\sin\left(\int_{s_0}^{s}\vert\alpha(t)\vert dt\right).
\end{align}
If $\xi$ is timelike then general solution of \eqref{os8} is given by
\begin{align}\label{os10}
y(s)=B_1\cosh\left(\int_{s_0}^{s}\vert\alpha(t)\vert dt\right)+B_2\sinh\left(\int_{s_0}^{s}\vert\alpha(t)\vert dt\right).
\end{align}
Since $\vartheta(s)\ne \vartheta(s_0)\exp\left(\pm\int_{s_0}^{s}\vert\alpha(t)\vert dt\right)$, therefore we find that $B_1\ne\pm B_2$.
\begin{rem}
In light of Eqs. \eqref{os9} and \eqref{os10}, it can be seen that Legendre curve given by \eqref{k2} in $\mathcal{Q}^3_{\alpha}$ is not spherical.
\end{rem}


\section{Introduction}

Cosmological systems are often modeled as perturbations around a
homogeneous, isotropic Friedmann-Lema\^{i}tre-Robertson-Walker (FLRW) spacetime where the background dynamics are described by homogeneous and isotropic stress-energy sources.
Yet, the homogeneity of the stress-energy tensor is manifestly
broken on smaller length scales where discrete objects exist,
and inhomogeneous structures and well-isolated astrophysical systems predominate.
A notion of homogeneity and isotropy can still be recovered through averaging,
in which small-scale structures in the Universe are coarse-grained and an
effective FLRW cosmology emerges.
Observations support the picture that we live in an approximately homogeneous and isotropic FLRW universe\cite{Aghanim:2018eyx,Komatsu_2011};
what remains uncertain is the precise relationship between the actual Universe, with its extreme inhomogeneity on many scales, and the perturbed homogeneous isotropic cosmology we use as a model.

In recent years, numerical studies have begun to explore some of the differences
between these pictures\cite{Giblin:2015vwq,Bentivegna:2015flc}, giving rise to many questions including: how do beams of light
and gravitational waves propagate in a warped cosmological vacuum rather than
a perturbed perfect fluid\cite{0907.4109,1204.2411}?
What microphysics best describes the
manner in which isolated objects contribute to global cosmological expansion?

In this work, we begin to examine the latter of these questions by simulating a lattice of
spinning black holes and examining properties of the spacetime, which is shown to develop
an overall average FLRW-like cosmological expansion.
Black hole lattice models have been employed as toy models for cosmological systems in order to ask
such questions in the past \cite{0907.4109,1204.2411,1204.3568,1306.1389,1306.4055,1307.7673,1309.2876,1404.1435}
(and see \cite{1801.01083} for a recent review), and these models and similar
semianalytic models have been able to provide insights into both the physics of spatial
hypersurfaces in these models and, more recently, observables \cite{Bentivegna:2016fls,Sanghai:2017yyn,Fleury:2017owg}.
Such models have been found to reproduce FLRW-type behavior with varying degrees
of fidelity, with properties that depend on the precise details of the
inhomogeneous structure \cite{Korzynski:2013tea,Korzynski:2014nna,Fleury:2018cro,Durk:2017rky}. This dependence
is interesting in and of itself, as it suggests that the cosmological properties
of inhomogeneous spacetimes do not always provide insight into the more
fundamental small-scale properties of a spacetime. For example, different measures
of the mass contained within these spacetimes have been found to disagree by
orders of magnitude: some definitions of mass appear to coincide
with FLRW expectations, while others do not \cite{1203.6478,1204.3568,1405.3197}.

Here we extend these models to a lattice of spinning (rather than purely static) black holes.
We lay down initial conditions and follow the subsequent evolution of the spacetime
using numerical general-relativistic simulations. For black holes parametrized by a
given mass and spin, we examine the expansion rate within the box, and explore
how different energy components contribute to cosmological
expansion. Although we do not calculate cosmological observables here,
this work lays down the foundation for a series of future work regarding
observational consequences.

We first describe our procedure for setting initial conditions with a
spinning black hole in a periodic spacetime in Section~\ref{sec:ICs}.
This is motivated by previous solutions
found within the conformal-transverse-traceless decomposition \cite{1204.2411},
now extended to obtain a solution similar to the Bowen-York solution
in the vicinity of the black hole. We provide details of the numerical
scheme used to evolve the spacetime in Section~\ref{sec:evolution}, and review
definitions of mass useful for characterizing properties of the spacetime.
In Section~\ref{sec:results}
we describe the different contributions to the Hamiltonian constraint equation,
showing how various terms contribute to cosmological expansion.
We find that our initial conditions contain a substantial anisotropic energy
density exterior to the black hole that quickly decays; the remaining spinning
black hole sources curvature which continues to give rise to expansion within
the lattice.
We then evaluate the behavior of different statistical measures of lattice
properties, concluding that volume-averaged properties appropriately describe the behavior,
and finding that the horizon mass, including both the irreducible (bare) mass and the
angular momentum, is sufficient for describing the observed expansion.
Lastly, we examine the averaged expansion rate, and compare this to the
cosmological expansion rate one might infer based on the mass of the black hole.
We find that the expansion rate initially behaves as a mixture of matter and radiative
content, consistent with residual radiation present in the initial data, with
the radiative content decaying and matter-dominated behavior emerging.

\section{Creating a spinning-black-hole lattice cosmology}
\label{sec:ICs}

We begin by reviewing the 3+1 decomposition of Einstein's equations, and
writing the constraint equations from this formalism in a form suitable
for numerically setting initial conditions with spinning black holes. We
restrict this discussion to vacuum solutions, although this formalism can
be generalized to include stress-energy sources. We will in particular make
use of the conformal transverse-traceless (CTT) decomposition of Einstein's
equations \cite{York:1978}, which extends the standard 3+1 decomposition,
in order to obtain solutions on spatial hypersurfaces.

We begin by writing the line element as
\begin{equation}
  ds^{2}=-\alpha^{2}dt^{2}+\gamma_{ij}\left(dx^{i}+\beta^{i}dt\right)\left(dx^{j}+\beta^{j}dt\right)\,.
\end{equation}
The non-dynamical Einstein's equations, projected onto spatial hypersurfaces
described by this metric, can be written as
\begin{eqnarray}
\label{eq:constraints}
  R + K^2 - K_{ij}K^{ij} &=& 0, \\
  D_j K^{j}_{\;i} - D_i K & = & 0,\nonumber
\end{eqnarray}
respectively known as the Hamiltonian and momentum constraint equations.
The derivatives $D_i$ are covariant with respect to the 3-metric $\gamma_{ij}$
  and $R$ is the associated 3-dimensional Ricci scalar.
The extrinsic curvature, $K_{ij}$, can be further decomposed into its trace, $K$,
and a traceless tensor, $A_{ij}$,
\begin{eqnarray}
\label{eq:2}
K_{ij} = A_{ij} + \frac13 \gamma_{ij} K.
\end{eqnarray}
In an FLRW model, the trace, $K$, parameterizes the Hubble expansion rate with $H_{\rm FLRW} = -K/3$,
while $K$ is instead zero on time-symmetric hypersurfaces, for example
asymptotically flat spacetimes in appropriate coordinates\cite{Bowen:1980yu}.

The 3-metric can also be conformally decomposed,
$\gamma_{ij} = \Psi^4 \tilde{\gamma}_{ij}$, $A^{ij} = \Psi^{-10} \hat{A}^{ij}$,
allowing us to rewrite the constraint equations (\ref{eq:constraints}) in terms
of these new variables,
\begin{eqnarray}
\label{eq:basic_constraints}
  &\tilde{D}_i \tilde{D}^i \Psi - \frac{1}{8} \tilde{R} \Psi + \frac{1}{8} \hat{A}_{ij} \hat{A}^{ij} \Psi^{-7} + 2 \pi \Psi^5 - \frac{1}{12} K^2 \Psi^5 = 0 \\
  &\tilde{D}_j \hat{A}^{ij} - \frac{2}{3} \Psi^6 \tilde{D}^i K = 8 \pi \Psi^{10} S^i\,, \nonumber
\end{eqnarray}
where $\tilde{D}_i$ and $\tilde{R}$ are now associated with the conformal metric
$\tilde{\gamma}_{ij}$.
The CTT decomposition
further breaks $\hat{A}^{ij} $ into longitudinal and transverse pieces,
\begin{eqnarray}
\label{eq:15}
\hat{A}^{ij} = \hat{A}^{ij}_{L} + \hat{A}^{ij}_{TT}.
\end{eqnarray}
Here $\hat{A}^{ij}_{TT}$ is transverse and traceless, satisfying $\tilde{D}_j\hat{A}^{ij}_{TT} = 0$.
The longitudinal piece $\hat{A}^{ij}_{L}$ can be written in terms of a vector $X^i$ as
\begin{eqnarray}
\label{eq:ALij}
\hat{A}^{ij}_{L}  = \tilde{D}^i X^j + \tilde{D}^j X^i - \frac{2}{3} \tilde{D}_k X^k \tilde{\gamma}^{ij} \equiv \left( \tilde{L} X \right)^{ij}.
\end{eqnarray}
The transverse-traceless component $\hat{A}^{ij}_{TT}$
contains information about transverse gravitational radiation, and can be set to zero
in order to minimize the gravitational radiation content of a solution. However,
this will not completely eliminate gravitational radiation, which can be sourced
nonlinearly, especially in a strong-gravity regime such as we are considering here. The
longitudinal component, on the other hand, contains information about the ``vector mode''
content of the spacetime, including frame-dragging and anisotropic effects. Generally,
vector modes are ignored in a cosmological setting, but their presence here will be
important for obtaining spinning-black-hole solutions.

Based on the above arguments, we further simplify the constraints by
setting $\hat{A}_{ij}^{TT} = 0$ , and choosing the metric to be 3-conformally flat,
$\tilde{\gamma}_{ij} = \delta_{ij}$. We then obtain
\begin{eqnarray}
\label{eq:simplifiedconstraints}
  &\nabla^2\Psi + \frac{1}{8} \left( \tilde{L} X \right)_{ij} \left( \tilde{L} X \right)^{ij} \Psi^{-7} - \frac{1}{12} K^2 \Psi^5 = 0 \nonumber \\
  &\nabla^2 X^{i} + \frac13 \partial^i \partial_j X^j- \frac{2}{3} \Psi^6 \partial^i K = 0,
\end{eqnarray}
where $\nabla^2$ is the Cartesian Laplacian and
\begin{eqnarray}
\label{eq:6}
\left( \tilde{L} X \right)^{ij} = \partial^i X^j  + \partial^j X^{i} - \frac23 \delta^{ij} \partial_k X^k.
\end{eqnarray}

When $K = 0$, and in an asymptotically flat spacetime, a solution for $X^i$ known as the Bowen-York solution is given by
\begin{eqnarray}
\label{eq:approximateKerr}
X^i=\tilde{\epsilon}^{ijk}\frac{x_j J_k}{r^3},
\end{eqnarray}
where $x^i$ are the Cartesian coordinates, $r$ is the coordinate distance from the origin, and
$J^i$ is a vector satisfying $\tilde{D}_i J_j = 0$.
Here $\tilde{\epsilon}^{ijk} \equiv \sqrt{\tilde{\gamma}} \epsilon^{ijk}$ is the 3D Levi-Civita tensor associated with the conformal metric
$\tilde{\gamma}_{ij}$, so that $\tilde{D}_i \tilde{\epsilon}^{ijk}=0$.

Substituting this solution (\ref{eq:approximateKerr}) into expression (\ref{eq:ALij}),
we obtain
\begin{eqnarray}
\label{eq:ALijBowenYork}
\hat{A}^{ij}_{L} = \left( \tilde{L} X \right)^{ij}=\frac{6}{r^3}x^{(i}\tilde{\epsilon}^{j)kl}J_k x_l / r^2.
\end{eqnarray}
This solution is commonly considered to contain a spinning black hole with spin $J_k$
\cite{Bowen:1980yu}.
The value of $\hat{A}^{ij}_{L}$ given by (\ref{eq:ALijBowenYork})
agrees with that of a
Kerr black hole at spatial infinity, implying this is true, however near the black hole this solution is not
equivalent to the Kerr metric. The Bowen-York solution has been found to contain some residual gravitational
radiation, and a maximum possible spin of $\|J\|=0.93$ \cite{gr-qc/9710096}.
We can nevertheless use this solution as inspiration for constructing initial conditions
in a cosmological setting, where the spacetime is no longer asymptotically flat.

Due to the discontinuity at the boundary, the Bowen-York solution (\ref{eq:approximateKerr}) is incompatible with periodic boundary conditions.
However, following a procedure similar to \cite{1204.2411, 1204.3568},
we can regularize the solution by multiplying parameters $M$ and $J^{i}$ in the metric fields by a {\it transition function}
\begin{eqnarray}
\label{eq:transnfn}
W(r; \sigma, l) = \left\{
             \begin{array}{lr}
             0 &  0 \leq r < l\\
             ((r-l-\sigma)^6\sigma^{-6}-1)^6 & l\leq r < l+\sigma\\
             1 & l+\sigma \leq r
             \end{array}                                                            \right.,
\end{eqnarray}
such that $W = 0$ at the origin, and transitions to $W = 1$ over some
distance scale $\sigma$ beginning at $r=l$. The vector $X^i$ is then regularized as
\begin{eqnarray}
\label{eq:XiBYinPB}
X^i\approx\tilde{\epsilon}^{ijk}\frac{x_j J_k}{r^3} (1-W(r)).
\end{eqnarray}
We can also regularize the solution to the conformal factor $\Psi$ as
\begin{eqnarray}
\label{eq:PsiBYinPB}
\Psi\approx 1 + \frac{ M }{2r} (1-W(r)).
\end{eqnarray}
Eqs. \ref{eq:XiBYinPB} and \ref{eq:PsiBYinPB} can then be used as an initial guess for
solving the constraint equations.

\begin{figure}[htb]
  \centering
    \begin{subfigure}{0.47\textwidth}
      \includegraphics[width=0.99\textwidth]{u_ps.png}
    \end{subfigure}
    \begin{subfigure}{0.47\textwidth}
      \includegraphics[width=0.99\textwidth]{X__i_ps.png}
    \end{subfigure}
    \caption{2D slices show the corrections to the initial guesses when solving Eq. \ref{eq:simplifiedconstraints} with a relaxation scheme.
      The initial guesses are given by Eqs. \ref{eq:XiBYinPB} and \ref{eq:PsiBYinPB} and the differences between those guesses and the
      exact solution are shown for field $\Psi$ (left) and $X^i$ (right).
      These also correspond to the fields $u$ and
      $X'^i$ as defined in Eq.~\ref{eq:16}. The choice of parameters for this solution correspond to R3 in Table \ref{t1}.
      } \label{init_u_X}
\end{figure}

In order for $K$ to be a negative constant at the boundaries, corresponding
to FLRW-like cosmological expansion, and zero in the center, corresponding to
a black hole solution, we also modulate the extrinsic curvature using $W(r)$,
\begin{eqnarray}
\label{eq:KinPB}
K = K_c W(r),
\end{eqnarray}
where $K_c$ is a constant, similar to \cite{1306.1389}.

By plugging Eq. \ref{eq:KinPB} into the Hamiltonian and momentum
constraint equations and taking the approximate solution from Eqs.
\ref{eq:XiBYinPB} and \ref{eq:PsiBYinPB} to be an initial guess with $J_k = (0, 0, a)$
in Cartesian coordinates, we can proceed to solve for $\Psi$ and $X^i$.
The singularities in the solution are avoided by employing the so-called {\it puncture approach};
further details regarding this can be found in \ref{sec:init}.
In Fig. \ref{init_u_X} we show snapshots of the absolute difference (fields $u$ and $X'^i$
from Eqs. \ref{eq:16}) between the exact solution to Eqs. \ref{eq:simplifiedconstraints} and
the approximate solution of Eqs. Eq. \ref{eq:XiBYinPB} and \ref{eq:PsiBYinPB}. The main change
seen for $\Psi$, which is initially $\mathcal{O}(1)$ at the boundaries and much larger near the black
hole as per Eq.~\ref{eq:PsiBYinPB},
is an overall distortion of the physical volume of the spacetime, with additional radial corrections.
For $X^i$, the predominant correction is a large radial contribution in the transition region.

\section{Lattice Evolution}
\label{sec:evolution}

We solve the initial constraints and evolve the spacetime using
the grid-based numerical relativity code \textsc{CosmoGRaPH} \cite{Mertens:2015ttp}.
We first solve the constraint equations using an integrated elliptical-equation
solver, which employs a standard Full Multigrid (FMD) iteration scheme and an
inexact-Newton-relaxation method \cite{Press:2007:NRE:1403886}. We verify the resulting initial
conditions by checking that the Hamiltonian and momentum constraint equations are satisfied
with increasing precision as resolution is increased.

After setting initial conditions, spatial slices are advanced using the
BSSNOK formulation of numerical relativity \cite{Nakamura:1987zz,Shibata:1995we,Baumgarte:1998te},
with 4th order Runge-Kutta timestepping. All fields are discretized as
cell-centered data, and centered 4th-order finite-difference stencils are
used for all derivatives except for advection terms $\sim \beta^i\partial_i$, where
upwind derivatives are used instead. Note that the BSSNOK scheme evolves
$\tilde{A}_{ij} = \Psi^{-4} A_{ij}$,
rather than the $\hat{A}_{ij}$ defined when setting initial conditions.

The gauge condition used in our simulation is a revised version
of the widely employed ``1+log'' and ``Gamma-driver'' gauge condition:
\begin{eqnarray}
\label{eq:gaugecondition}
  \partial_t \alpha &=& -2 \eta \alpha (K - \langle K \rangle_{\rm edge}) + \beta^i \partial_i \alpha, \nonumber \\
  \partial_t \beta^i  & = & B^i,  \\
  \partial_t B^i & = & \frac34 \partial_t \Gamma^i - B^i. \nonumber
\end{eqnarray}
This differs from the usual ``1+log'' gauge by introducing a reference
expansion rate, $\langle K \rangle_{\rm edge}$, which is
the conformal average of the extrinsic curvature $K$ along all edges of
the computational domain box defined in Eq. \ref{eq:21}.
This gauge choice has been demonstrated to have powerful singularity-avoidance properties\cite{Brown:2007pg};
the modification we make by subtracting $K$ relative
to the average boundary value allows the spatial slice to be driven
towards FLRW-like expansion away from the black hole.

Because of the collapsing nature of the black hole,
we have also integrated an Adaptive Mesh Refinement (AMR) framework
into the time evolution,
provided by the code SAMRAI \cite{DBLP:conf/sc/WissinkHKSE01},
an open-source structured adaptive-mesh-refinement application infrastructure.
By building hierarchies of grid levels with different resolution, dividing and
distributing patches into computational nodes, SAMRAI realizes high-efficiency
adaptive-mesh refinement and parallelization. To synchronize data on different levels,
we use tri-cubic Hermite interpolation\cite{kreyszig2019advanced}, which we find results in a high
degree of numerical stability.

We also require a technique to locate apparent horizons. Because of the aspherical
nature of spinning black holes, there is no symmetry of the apparent horizon
that makes it simple to locate. We therefore use the AHFinderDirect package
\cite{Thornburg:2003sf} to find the apparent horizon on a given spatial hypersurface.
The definition and the method of extracting angular momentum from
an isolated horizon come from
\cite{Dreyer:2002mx}, where the angular momentum of a black hole $J$ is defined as
\begin{eqnarray}
\label{eq:19}
J = \frac{1}{8\pi} \oint (\varphi^a R^b K_{ab}) \mathrm{d}^2V.
\end{eqnarray}
Note here $R^b$ is an outgoing vector normal to the horizon and
$\varphi^a$ is not a Killing vector of the full spacetime but a symmetry
vector defined locally on the horizon that preserves the induced metric $q_{ab}$,
so that
\begin{eqnarray}
\label{eq:induced}
 \mathcal{L}_{\varphi} q_{ab} = 0
\end{eqnarray}
(see \cite{Dreyer:2002mx} for more detail). The eigenvalue closest to unity
associated with the symmetry vector is within a percent of unity, indicating
the spacetime is very close to axisymmetric in the vicinity of the black hole.

We then track the black hole's irreducible mass, spin, and horizon mass, as well as
the expansion history of the spacetime.
The irreducible and horizon masses are defined as
\begin{eqnarray}
\label{eq:24}
M_{\rm H}^2 & \equiv & M_{\rm irr}^2 + \frac{J^2}{4M_{\rm irr}^2} \nonumber \\
M_{\rm irr}^2 & \equiv & \mathcal{A}/16\pi\,.
\end{eqnarray}
Here $\mathcal A$ is the area of the horizon, defined as
$\mathcal{A} \equiv \oint \sqrt{q} \mathrm{d}^2V$, where $q$ is the determinant
of the induced metric on the horizon.

In order to examine how well the spinning-black-hole-lattice universe
corresponds to a FLRW universe with similar expansion properties, or to
check how well the lattice obeys a Friedmann-like equation, we need to
define an effective density of the spacetime, $\rho_{\rm eff}$, and an average spacetime
expansion rate, $\left<K\right>$.
We then define a dimensionless parameter
\begin{equation}
\label{eq:FLRW}
\mathcal{C} \equiv \frac{\rho_{\rm eff} }{\langle K \rangle^2 / 24\pi}\,.
\end{equation}
According to the Friedmann equation, one should have $\mathcal{C}=1$
for an appropriately chosen ${\rho_{\rm eff} }$
and $\langle K \rangle$ as the effective Hubble parameter equals $\langle K \rangle /3$.

We first consider whether, in defining ${\rho_{\rm eff} }$
and $\langle K \rangle$, it is more appropriate to use a volume-averaging operation or to
average over edges of the box.
We can set $\rho_{\rm eff} = \rho_{\rm edge,\,eff}$, given by
\begin{equation}
\label{eq:23}
\rho_{\rm edge,\,eff} \equiv M_{\rm eff} / (D_{\bot}^2 D_{\parallel}).
\end{equation}
$M_{\rm eff}$ can be either $M_{\rm H}$ or $M_{\rm irr}$, while
$D_{\bot}$ and $D_{\parallel}$ are distances along edges of the box in directions that
are perpendicular and  parallel to the spin direction respectively,
$D_i \equiv \int \mathrm{d}x^i \sqrt{\gamma_{ii}}$ for a Cartesian direction $i$
(where no sum over $i$ is implied).

An alternative is to make $\rho_{\rm eff} = \rho_{\rm vol,\,eff}$, with
\begin{equation}
\label{eq:rhovol}
\rho_{\rm vol,\,eff} \equiv M_{\rm eff} / V,
\end{equation}
where
\begin{equation}
\label{eq:vol_defn}
V \equiv \int_{r>r_H} \mathrm{d}^3 x \sqrt{\gamma}
\end{equation}
is the conformal volume exterior to the black-hole horizon.

More generically,
we define an averaged physical quantity $Q$ on the edge or volume
\begin{eqnarray}
\label{eq:21}
  \left< Q\right>_{\rm edge} &\equiv&  \frac{\sum_{\rm all\;edges} \int \mathrm{d}x^i \sqrt{\gamma_{ii}} Q }{\sum_{\rm all\;edges} D_i} \\ \nonumber
  \left< Q \right>_{\rm vol} &\equiv& \frac{ 1 }{V} \int_{r>r_H} \mathrm{d}^3 x \sqrt{\gamma} Q\,.
\end{eqnarray}
We can now define the ratio of the left and right hand sides of the effective Friedmann equations as
\begin{eqnarray}
\label{eq:20}
  \mathcal{C}_{\rm edge,\,eff} = \frac{\rho_{\rm edge,\,eff}}{\langle K^2 \rangle_{\rm edge}  / 24 \pi}, \;\;\;\;\;\; \mathcal{C}_{\rm vol,\,eff} = \frac{\rho_{\rm vol,\,eff}}{\langle K^2 \rangle_{\rm vol}  / 24 \pi}.
\end{eqnarray}

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

In this section, we will present our main result. We will mainly focus on the
the expansion history and effects of spins.
We will also investigate time evolution of dimensionless parameter $\mathcal{C}$ and evaluate
the effect of statistics.
\subsection{Initial condition effects on physical lattice properties}

The free parameters in our setup are the box size $L$, the mass scale $M$, the spin $a$,
the extrinsic curvature at the periodic boundary $K_c$,
and the parameters appearing in the transition function (\ref{eq:transnfn}), $l$, and $\sigma$.
To demonstrate the impact that varying these parameters has on the physical
properties of the spacetime (namely the black hole masses, densities, and $\mathcal{C}$),
we have listed these properties and  the corresponding
parameters in Table \ref{t1} for ten representative simulations.

\begin{table}[htbp]
  \caption{\label{t1} Parameters of initial setups}
  \begin{center}
  \item[]\begin{tabular}{l|ccccc|ccccc|c}
    \mr
    Runs &$L$& $a$& $K_c$ &$l$&$\sigma$& $M_{\rm irr}$ &$M_{\rm H}$       & $M_{\rm H}/V$&$\mathcal{C}_{{\rm edge},\,H}$ & $\mathcal{C}_{{\rm vol},\,H}$ & $c$\\
    \mr
    R1   &10  & 0  & -0.21 & 1 & 3.5    & 1.023    & 1.023  & 0.000414&0.897& 0.829 & 7.7\\
    R2   &10  & 0.6& -0.21 & 1 & 3.5    & 1.107    & 1.140  & 0.000404&0.894& 0.817 & 7.7\\
    R3   &10  & 0.9& -0.21 & 1 & 3.5    & 1.180    & 1.240  & 0.000395&0.891& 0.807 & 7.4\\
    R4   &10  & 0  & -0.21 & 0.1 & 3    & 0.982    & 0.982  & 0.000467&0.996& 0.844 & 8.3\\
    R5   &10  & 0.6& -0.21 & 0.1 & 3    & 1.072    & 1.108  & 0.000456&0.994& 0.828 & 7.4\\
    R6   &10  & 0.9& -0.21 & 0.1 & 3    & 1.147    & 1.212  & 0.000446&0.993& 0.815 & 8.1\\
    R7   &10  & 0.6& -0.15 & 1 & 3.5    & 1.246    & 1.269  & 0.000216&0.900& 0.842 & 7.6\\
    R8   &10  & 0.6& -0.1  & 1 & 3.5    & 1.453    & 1.468  & 0.000101&0.904& 0.865 & 7.6\\
    R9   &11  & 0.6& -0.21 & 1 & 3.5    & 1.058    & 1.096  & 0.000420&0.921& 0.827 & 7.4\\
    R10  &12  & 0.6& -0.21 & 1 & 3.5    & 1.018    & 1.060  & 0.000432&0.941& 0.832 & 7.4
  \end{tabular}
  \end{center}
\end{table}
In this table, only the first 2nd-6th columns are free parameters that were
chosen initially, while the 7th-11th columns are derived parameters that can only be calculated
after initial constraints are fully solved. The $\rho_{\rm eff}$ used
in columns for $\mathcal{C}_{{\rm edge},\,H}$ and $\mathcal{C}_{{\rm vol},\,H}$ are calculated using $M_{\rm H}$.
The convergence rate for each run on initial slice is represented by parameter $c$ whose definition can be found in Eq. \ref{eq:c}.
In all runs, we do not vary $M$, instead choosing to work in units where $M=1$.

Examining these initial configurations, we can observe the following:
\begin{itemize}
\item Although the input parameter $M$ is equal to $M_{\rm H}$ for a Kerr spacetime
  with asymptotic flat boundary, the resulting $M_{\rm H}$ in the table only roughly tracks $M$,
  depending on other parameters as well.
\item Both  $\mathcal{C}_{{\rm edge},\,H}$
  and $\mathcal{C}_{{\rm vol},\,H}$ are somewhat less than
  1 initially, and change very little when
  the spin parameter $a$, box size $L$, or boundary extrinsic curvature $K_c$ are varied. This
  implies that the initial spatial slice is always ``close'' to FLRW.
\item Only by changing the combination of $l$ and $\sigma$ does $\mathcal{C}_{{\rm edge},\,H}$
  change significantly; However, the value of $\mathcal{C}_{{\rm vol},\,H}$ still does
  not change.
\item Increasing the box size $L$ (comparing R9 and R10 to R2)  does not increase the physical size of the box.
\item Changing the boundary extrinsic curvature $K_c$ (comparing R7 and R8 to R2) will
  change the effective density and physical box length significantly, but still keeps the ratio in the last
  two columns unchanged.
\end{itemize}

To summarize, the parameters that predominantly determine physical properties of the system are $K_c$, $M$, and $a$. These
strongly affect the simulation volume and black hole mass and spin. $L$, $l$, and $\sigma$
instead affect the coordinate description of the spacetime, and only weakly affect physical properties.

The initial value of the ratio $\mathcal{C}$ quantifies the deviation from the FLRW universe,
and is found to be relatively independent of our parameter choices.
We now wish to study its time evolution as well as the best way to fit
$\rho_{\rm eff}$ to a FLRW universe. Those topics are our main interests and will be
discussed in the following sections.

\subsection{Expansion properties and energy content}

We can now analyze the different contributions to the Hamiltonian constraint
equation, or the different ``energy'' contributions in the spacetime contributing
to expansion. Because we work in a vacuum spacetime, there is  no actual stress-energy
contribution, and all expansion must be a result of either curvature or kinetic terms
in the constraint equations. We can decompose these terms as in Hamilotnian constraint in Eq.~\ref{eq:basic_constraints} to get
\begin{eqnarray}
\label{modified_constraint}
 R / 8  +  \tilde{A}_{ij} \tilde{A}^{ij} / 8 +  K^2 / 12 = 0,
\end{eqnarray}
and analyze the average behavior of: the curvature $\langle R \rangle / 8$; the anisotropic expansion term
$\langle \tilde{A}_{ij} \tilde{A}^{ij}\rangle / 8$, which contains contributions from
vector and tensor modes and their interactions; and the expansion itself,
$\left< K^2\right> / 12$.

\begin{figure}[htbp]
  \centering
  \begin{subfigure}{0.3\textwidth}
      \includegraphics[width=0.99\textwidth]{a_0_AijAij_t_0_v5_ps.png}
    \end{subfigure}
    \begin{subfigure}{0.3\textwidth}
      \includegraphics[width=0.99\textwidth]{a_06_AijAij_t_0_v5_ps.png}
    \end{subfigure}
    \begin{subfigure}{0.37\textwidth}
      \includegraphics[width=0.99\textwidth]{a_09_AijAij_t_0_v5_ps.png}
    \end{subfigure}
    \begin{subfigure}{0.3\textwidth}
      \includegraphics[width=0.99\textwidth]{a_0_AijAij_t_10_v5_ps.png}
    \end{subfigure}
    \begin{subfigure}{0.3\textwidth}
      \includegraphics[width=0.99\textwidth]{a_06_AijAij_t_10_v5_ps.png}
    \end{subfigure}
    \begin{subfigure}{0.37\textwidth}
      \includegraphics[width=0.99\textwidth]{a_09_AijAij_t_10_v5_ps.png}
    \end{subfigure}

    \caption{The conformally related trace-free part of the extrinsic curvature,
      $\tilde{A}_{ij}\tilde{A}^{ij}$, on two-dimensional slices that intersect
      the black hole through spatial hypersurfaces. The black circles are
      corresponding slices of apparent horizons.
      $\tilde{A}_{ij}\tilde{A}^{ij}$ is shown
      runs R1, R2 and R3 with spins $a = 0, 0.6, 0.9$ (left to right), and at times $t=0, 10$  (top to bottom).
      All quantities are in units where $M=1$.
    }
  \label{AijAij}
\end{figure}

In Fig.~\ref{AijAij}, we examine the contribution of the $\tilde{A}_{ij}\tilde{A}^{ij}$
term, which contains information about gravitational-wave and vector-mode energy content.
It shows that the vector and tensor modes are concentrated
near the black hole horizon, especially as the spacetime evolves and relaxes away from the
na\"ive initial conditions that we set. This is unsurprising in that vector and tensor modes
can be sourced nonlinearly in a strong-gravity regime, in contrast to the linear regime in a
cosmological setting, where they are expected to be negligible.
The absence of this contribution further away from the black hole shows that scalar
curvature is the dominant contribution to cosmological expansion.

\begin{figure}[htbp]
  \centering
      \begin{subfigure}{0.3\textwidth}
      \includegraphics[width=0.99\textwidth]{a_0_Ksq_t_0_v4_ps.png}
    \end{subfigure}
        \begin{subfigure}{0.3\textwidth}
      \includegraphics[width=0.99\textwidth]{a_06_Ksq_t_0_v4_ps.png}
    \end{subfigure}
    \begin{subfigure}{0.37\textwidth}
      \includegraphics[width=0.99\textwidth]{a_09_Ksq_t_0_v4_ps.png}
    \end{subfigure}
    \begin{subfigure}{0.3\textwidth}
      \includegraphics[width=0.99\textwidth]{a_0_Ksq_t_10_v4_ps.png}
    \end{subfigure}
        \begin{subfigure}{0.3\textwidth}
      \includegraphics[width=0.99\textwidth]{a_06_Ksq_t_10_v4_ps.png}
    \end{subfigure}
    \begin{subfigure}{0.37\textwidth}
      \includegraphics[width=0.99\textwidth]{a_09_Ksq_t_10_v4_ps.png}
    \end{subfigure}
    \begin{subfigure}{0.3\textwidth}
      \includegraphics[width=0.99\textwidth]{a_0_Ksq_t_15_v4_ps.png}
    \end{subfigure}
        \begin{subfigure}{0.3\textwidth}
      \includegraphics[width=0.99\textwidth]{a_06_Ksq_t_15_v4_ps.png}
    \end{subfigure}
    \begin{subfigure}{0.37\textwidth}
      \includegraphics[width=0.99\textwidth]{a_09_Ksq_t_15_v4_ps.png}
    \end{subfigure}
    \caption{Time evolution snapshots of $K^2$ in R1, R2 and R3 with
$a = 0, 0.6, 0.9$ (from left to right), and $t=0 (M), 10 (M), 15(M)$
(from top to bottom).}
  \label{Ksq}
\end{figure}

Fig. \ref{Ksq} similarly depicts the time evolution of $K^2$, which is the volume expansion
rate. A behavior very similar to $\tilde{A}_{ij} \tilde{A}^{ij}$ is identified --
deviations from cosmological-type expansion are found near the black hole, that gradually
become smooth far away from the black hole, especially as the simulation progresses.

Residual oscillations can be seen in the expansion rate, both spatially varying, and
as a function of time. The behavior of these oscillations depends on both the initial
conditions and the gauge we choose, and we therefore do not consider these to be indicative of
an expansion rate that is physically oscillatory, ie. that would strongly impact the way
a geodesic observer would view the spacetime. We leave this speculation to future work,
although see also \cite{1801.01083}.

\begin{figure}[htbp]
  \centering
    \includegraphics[width=0.99\textwidth]{p1.pdf}
  \caption{Behavior of the contributions $\left<R\right> / 8$, $\langle \tilde{A}_{ij}\tilde{A}^{ij}\rangle/8$ and $\left<K^2\right>/12$ to the Hamiltonian constraint in Eq. \ref{modified_constraint} in R2. The left panel shows edge-averaged terms, and the right panel shows results with volume-averaged terms. }
  \label{Aij_avg}
\end{figure}

We examine the different contributions to the Hamiltonian constraint equation more
quantitatively for R2 in Fig. \ref{Aij_avg}. This demonstrates that the vector and tensor
contributions $\langle \tilde{A}_{ij}\tilde{A}^{ij}\rangle$
are relatively small near the edge,
but are appreciable when averaged over the volume exterior to the horizon.
It is important to note that this interpretation will be affected by gauge choice:
for example, the first-order gauge-invariant vector mode usually considered in a cosmological
setting, as well as true observables (eg. properties integrated along geodesics according
to observers), will contain a contribution from the shift that is not shown here.

\subsection{Expansion-mass correspondence}

As we evolve the spacetime using the gauge choice of Eq.~\ref{eq:gaugecondition},
lengths of edges as well as volume of
the spatial slice expand in the intuitively expected manner.
The coordinate
size of the black hole apparent horizon initially expands as the solution stabilizes, then
shrinks due to cosmological expansion, while the area stays the same.
We run the code until the black-hole horizon
becomes too small to be resolved accurately. At this moment, lengths in spatial slices
have roughly grown by a factor of $e^2$.

Because of the asymmetric setup, one might expect to see a different expansion rate in
different directions.
However, we find less than a $0.1\%$ difference in lengths along different edges of
the computational box ($D_{\parallel}$ and $D_{\bot}$),
and thus will ignore this discrepancy and focus on the quantities averaging on both
parallel and perpendicular edges (see Eq. \ref{eq:21}).

The masses $M_{\rm irr}$ and $M_{\rm H}$ on the initial slice of each runs are shown in Table \ref{t1}.
During the evolution, their time dependence is nearly negligible: the
relative fluctuation in their values is as small as $0.1\%$ and dominated by
numerical uncertainty (see \ref{sec:convergence} for more detail),
consistent with the area theorem and conservation of angular momentum.

\begin{figure}[htbp]
  \centering
    \includegraphics[width=0.99\textwidth]{p3.pdf}
    \caption{$\mathcal{C} = \rho_{\rm eff} / (\left<K\right>^2/24\pi)$ for runs R3 and R6.
      $\mathcal{C}$ is evaluated using edge and
      volume averages, with both $M_{\rm irr}$ and $M_{\rm H}$ used in $\rho_{\rm eff}$.
    }
  \label{ratios_M_H}
\end{figure}
\begin{figure}[htbp]
  \centering
    \includegraphics[width=0.99\textwidth]{p4.pdf}
    \caption{Behavior of the ratio $\mathcal{C}_{\rm vol}$ with different spins
    and masses. The horizon mass is used in the left panel, while the irreducible mass
    is used in the right.}
  \label{ratios_spins}
\end{figure}

From Fig. \ref{ratios_M_H}, we can see that
$\mathcal{C}_{\rm edge}$ fluctuates with a large amplitude, consistent with
behavior seen in \cite{1306.4055}, while $\mathcal{C}_{\rm vol}$ gently increases to unity.
The difference between choices of irreducible mass and horizon mass results in a constant
shift between curves, which we investigate below.

We can attempt to reduce the amplitude of oscillations in the edge-averaged case
by adjusting $l$ and $\sigma$ to obtain  $\mathcal{C}_{\rm edge}\simeq 1$ on the initial slice.
Comparing the panels in Fig.~\ref{ratios_M_H}, we see that this does not help in eliminating
the fluctuations.  The large amplitude fluctuations we see in $\mathcal{C}_{\rm edge}$ apparently arise from a combination of the way we slice the initial spacetime and the gauge condition, indicating
volume-averaged quantities appear to be a more appropriate representation of the physical behavior
of the system.

Fig. \ref{ratios_spins} provides us with more insight into
the spin dependence of $\mathcal{C}_{\rm vol}$ for different definitions of mass.
$\mathcal{C}_{\rm vol}$ is approximately spin-independent when the horizon mass is used to construct $\rho_{\rm eff}$. It appears to approach the expected matter-dominated FLRW value of unity as the simulation evolves.
Neither of these features are maintained when choosing $\rho_{\rm eff} = M_{\rm irr} / V$ (right panel).
The diminishing value of $\mathcal{C}_{\rm vol}$
as the spin is increased implies that an energy contribution to the Hamiltonian
constraint is not being accounted for.
The horizon mass $M_{\rm H}$ is thus a better
choice than $M_{\rm irr}$ when considering the effective mass in a spinning-black-hole-lattice universe.

\begin{figure}[tbp]
  \centering
    \includegraphics[width=0.99\textwidth]{p2.pdf}
    \caption{Dependence of $M_{\rm H} / M_{\rm H, 0}$
      and $\langle K^2 \rangle_{\rm vol}/\langle K^2 \rangle_{\rm vol, 0}$ as a function of scale factor.
      We also compare the averaged expansion rate to an FLRW model, according to a best fit function with radiation (top 3 lines) and without (bottom 3 lines), showing the actual expansion is well-described by a mix of matter and radiation.
    }
  \label{rho_K_vs_a}
\end{figure}

Lastly, we consider how the averaged values we consider here map to corresponding
FLRW spacetimes in Fig.~\ref{rho_K_vs_a}. As a function of
FLRW scale factor, $a \equiv V^{1/3}$, with V as in Eq.~\ref{eq:vol_defn}, we note that the horizon mass $M_{\rm H}$
is conserved, implying matter-domination-like expansion with $\rho_{\rm eff} \propto a^{-3}_{\rm FLRW}$.
There is no spin dependence to within numerical uncertainty. However, the time-dependence of the effective
Hubble parameter relative to an FLRW model shows a mix of matter and radiation contributions.
The bottom 3 lines in this figure imply that a purely matter dominated expansion $(a/a_0)^{-3}$ is not a good fit,
while the top 3 lines show the behavior is well-described by including radiative content with
a best fit function $0.19(a/a_0)^{-4} + 0.81(a/a_0)^{-3}$.

\section{Discussion and Conclusions}
\label{sec:conclusions}

In this paper, we have built a new kind of black-hole-lattice universe by
solving the constraint equations with a conformal-transverse-traceless decomposition with periodic boundary conditions.
A series of space-like hypersurfaces corresponding to expanding universes with spinning black holes
were identified. We found that the expansion features of those initial slices were very
close to an effective FLRW universe regardless of the initial parameters choices.

We then evolved the initial slices with a singularity-avoiding gauge choice,
finding no significant difference between the expansion in directions parallel and perpendicular to the spin.
The effective density described by the black hole mass evolved similarly to a matter-dominated universe,
ie. the mass of the black hole was conserved,
while the effective Hubble parameter only followed matter-dominated behavior at late times.
When quantifying the deviation of the expansion rate from FLRW-like behavior,
we found that averages taken over edges of the simulation coordinate box displayed large fluctuations, while volume averaged quantities showed much smaller deviations and approached FLRW asymptotically.
By fitting the spinning-black-hole-lattice universe to the FLRW universe, we were able
to identify the effective mass that governs the expansion as the horizon mass of the black hole,
rather than the irreducible mass.

In future work, we can track physical observables through the spacetime,
to better characterize the effects of highly non-linear non-stationary
perturbations on universes that, like are own, appear to be on-average
homogeneous and statistically isotropic on large scales.

It is noteworthy that the spacetimes that we have considered have a preferred direction, determined by the orientation of the spin of the single black hole in the fundamental domain.
We anticipate exploring more general initial conditions with multiple spinning black holes
and $0$ total angular momentum.
Finally, although the spurious gravitational waves
introduced by initial gauge fixing (conformally flat and $\hat{A}^{TT}_{ij} = 0$) have been shown to not critically affect the late-time evolution in both previous work \cite{gr-qc/9710096, Burko:2005fa} and in our observation of $\tilde{A}_{ij} \tilde{A}^{ij}$, it may nevertheless be helpful to explore other
schemes, like conformal thin-sandwich (CTS) decomposition, to set more general initial conditions with reduced spurious gravitational wave content.

\section*{Acknowledgments}

We thank Tim Clifton, Mikołaj Korzyński, and Eloisa
Bentivegna for insightful discussions that helped shape this work.
This work benefited from the Sexten
Center for Astrophysics workshop on GR effects in cosmological
large-scale structure, and made use of the High Performance
Computing Resource in the Core Facility for Advanced Research Computing
at Case Western Reserve University.
JBM acknowledges support as a CITA National Fellow and from Perimeter Institute
for Theoretical Physics. Research at Perimeter Institute is supported by
the Government of Canada through the Department of Innovation, Science
and Economic Development Canada and by the Province of Ontario through
the Ministry of Research, Innovation and Science.  JTG is supported
by the National Science Foundation Grant No. PHY-1719652.
GDS and CT were supported in part by grant DE-SC0009946 from the US DOE.

\section*{\refname}