\section{Introduction}\label{introduction}

In biological and medical research, there is often a need to compare the accuracy of different measurements during diagnosis or the effectiveness of different therapeutic schedules during treatment. To quantitatively compare the differences between two or more methods, the non-inferiority test is used to determine whether the experimental index generated by one method differs from the same index generated by other methods beyond a specified level of magnitude. If the magnitude of the difference is chosen to be zero, the problem becomes a problem of checking whether the same experimental index generated by two or more methods is equivalent, which is called the equivalence test.\cite{liu_ma_wu_tai_2006,guo_luh_2020,sandie_molinari_etal_2022} used a non-inferiority test to compare the accuracy of continuous data from two measurements, and \cite{paul_tiwari_etal_2021,ghosh_guo_ghosh_2022} used non-inferiority test to compare the effectiveness of continuous data from three treatments.

Every sample has a two-dimensional attribute, and each dimension has two possible observations, such data can be described by the bivariate binary distribution, which yields the ordinary $2 \times 2$ contingent table \cite{book_Bishop_2007, fagerland_lydersen_laake_2013_same_unite}.
For example, a patient receives two types of treatments successively, and the results of each treatment may be effective or ineffective. As the observations are from the same sample, there may be a correlation between the two treatment effects \cite{fay_lumbard_2020}. Another example is that the patients suffering from two different types of diseases are treated with the same drug, as the observed data come from two different groups, there may be a correlation between the group attribute of the samples and the drug efficacy\cite{wu_2018}. A third example is that patients with the same disease may have a positive or negative measurement result before receiving treatment, and the result of the same measurement method may also be positive or negative after treatment. Two times of measurement apply to the same sample, and there may be a correlation between the two measurement results \cite{fagerland_lydersen_laak_2013_sct}.
It is important to accurately analyze the bivariate binary data to compare the accuracy of measurement schemes, evaluate the effect of therapeutic schedules, distinguish the application groups of new drugs, and evaluate the adverse effects of drugs. Some researchers assume that there is no correlation between the two dimensions \cite{liu_ma_wu_tai_2006,guo_luh_2020,sandie_molinari_etal_2022,ghosh_guo_ghosh_2022,rohmel_2005}, and this simplification may influence the decision of non-inferiority test or equivalence test \cite{wu_2018,kang_chung_ahn_2005,nam_kwon_2009, C_T_Le_1984}. Therefore, this paper focuses on the equivalence test of $2 \times 2$ contingent table when there is a correlated relationship between the two dimensions.

There are many results on the probability structure and parameter estimation method for the bivariate binary distribution. \cite{C_T_Le_1984} studied a symmetric distribution derived from the bivariate binary distribution, determined the probability structure, and gave a parameter estimation algorithm under the assumption that the two marginal positive probabilities are equivalent. \cite{albert_ingram_1985} discussed a series of distributions derived from the bivariate binary distribution, and pointed out that the $2 \times 2$  contingent table is equal to a four-dimensional multinomial distribution with four joint probabilities.
\cite{manuel_weib_silva_2014,wu_rai_yan_2019} gave the upper bound and lower bound of the correlated coefficient using the marginal positive probabilities, and \cite{wu_rai_yan_2019} proposed a maximum likelihood estimation of the marginal probabilities and the correlated coefficient. Under the hypothesis of equivalence, this paper constructs a two-variable model for the parameter space using the correlated coefficient and the marginal positive probability, discusses the symmetric properties of the marginal positive probability and the discordant probability, establishes the monotonicity of the cumulative marginal distribution function of discordance observation, and determines the relationship between the minimum point of the cumulative marginal distribution function with the correlated coefficient.

Based on the exact or conditional distributions of the bivariate binary distribution, a series of unconditional or conditional p-value methods apply to the non-inferiority/equivalence test.
\cite{rohmel_2005} considered the constrained maximum likelihood estimation and the monotonicity of the proposed estimator.
\cite{kang_chung_ahn_2005} constructed the conditional sample space and proposed the exact/conditional p-value methods by removing the redundant parameters.
\cite{nam_kwon_2009} compared the performances of the adjusted score method, a Wald-type test, and a modification of Obuchowski's method for clustered matched-pair data. And it is observed that sufficient sample size improves the accuracy of the equivalence test.
Using the monotonicity of Barnard convex set, \cite{sidik_2003} proposed two p-value tests and found that the probability of the type I error may exceed the specified nominal significance level for the tests derived from the asymptotic distribution.
\cite{lloyd_moldovan_2011} considered the implementation order of three transforms and found that a proper order of combination improves the power of the p-value test.
\cite{li_liu_goldberg_2011} considered three kinds of p-value tests and calculated their exact size under the conditional monotonicity condition.
\cite{liu_hsueh_hsieh_chen_2002} compared the size and power of two asymptotic tests, a Wald-type test, and a restricted maximum likelihood estimation.
\cite{fagerland_lydersen_laake_2014} found the best tests were the asymptotical McNemar test and the McNemar mid-p test among 24 tests and confidence intervals.
Based on the marginal distribution of the two discordance observations and the estimated confidence region, this paper presents an equivalence test called the margin test. And the influence of sample size, correlation coefficient, and nominal significance level on the size and power of the proposed test is discussed in detail.

The McNemar test is the most frequently used method in the equivalence test for bivariate binary distribution. \cite{wu_2018} studied an adjusted McNemar test and calculated its power with simulated data. \cite{fagerland_lydersen_laak_2013_sct} compared the performances of a mid-p test, an asymptotic McNemar test, an exact conditional McNemar test, and an exact unconditional test through their size and power.
\cite{mohammadi_atashin_hofman_tan_2018} compared the accuracy of two ontology alignment systems over a matching task by the McNemar test. \cite{murtza_khan_akhtar_2019} compared the accuracy of several methods in target detection problems by the McNemar test. The wide application of the McNemar test in biomedical and other fields motivates us to discuss the size and power with different sample sizes, correlation coefficients, and significance levels in detail.
We also compare the performances of the McNemar test and the proposed margin test on two real-world data sets, where the observed data is on the boundary of the rejection regions.

The rest of this paper is organized as follows. Section 2 introduces the bivariate binary model and its notations.
The upper and lower bounds of the correlated coefficient are given in section 3.
Section 4 depicts the parameter space of the null hypothesis with the positive probability and the correlated coefficient. And the symmetric property of the positive probability and the discordant probability is also discussed with the proposed parameter space.
The monotonicity of the cumulative distribution function of the discordant observation is presented in section 4. And the relationship between the minimum point of the cumulative distribution function and the correlated coefficient is also discussed in section 4.
Section 5 gives a confidence region of the two discordance observations by their joint distribution.
And an equivalence test is deduced from the proposed confidence region with the marginal distribution of the discordant observation.
Section 6 briefly introduces the McNemar test and its rejection region.
The performance of the McNemar test and the margin test is compared through the range of acceptance region, size, and power in section 7.
Section 8 compares the performance of the McNemar test and the margin test with two real-world examples.
The summary and discussion are given in section 9.

\section{Bivariate Binary Model}\label{bivariate_binary_model}

\begin{table}[h]
	\caption{The Bivariate Binary Distribution $X=(X_1, X_2)$.}
	\label{bivariate_binary_table_model_1}
	\centering
	\begin{tabular}{l|cc|c}
		\hline
		\diagbox{$X_1$}{$X_2$} & 0 & 1 & \\
		\hline
		0 & $p_{00}$ & $p_{01}$ & $p_{0+}$  \\
		1 & $p_{10}$ & $p_{11}$ & $p_{1+}$  \\
		\hline
		& $p_{+0}$ & $p_{+1}$ &   \\
		\hline
	\end{tabular}
\end{table}

The population $X=(X_1,X_2)$ follows the bivariate binary distribution as Table \ref{bivariate_binary_table_model_1} \cite{book_Bishop_2007}.
Denote the joint distribution of the two-dimensional variable $(X_1, X_2)$ as
\begin{equation}\label{prob:binary_joint}
	p_{ij}=Prob\{X_1=i,X_2=j\}, \quad i,j=0,1.
\end{equation}
From the joint distribution (\ref{prob:binary_joint}), the marginal probabilities of $X_1$ and $X_2$ are
\begin{eqnarray}
	&p_{0+} = p_{00}+p_{01},   \quad & p_{1+} = p_{10}+p_{11}, \\
	&p_{+0} = p_{00}+p_{10},   \quad &p_{+1}=p_{01}+p_{11},  \\
	&p_{0+} = 1 - p_{1+},      \quad &p_{+0} = 1-p_{+1}.
\end{eqnarray}
The marginal distributions of $X_1$ and $X_2$ are determined by the two positive probabilities $p_{1+}$ and $p_{+1}$ respectively. Denote the variance and standard deviation of $X_1$ and $X_2$ as
\begin{eqnarray}
	&\sigma^2_1 = D(X_1) =  p_{1+}(1-p_{1+}), \quad &\sigma_1 = \sqrt{\sigma_1^2}, \\
	&\sigma^2_2 = D(X_2) =  p_{+1}(1-p_{+1}), \quad &\sigma_2 = \sqrt{\sigma_2^2}.
\end{eqnarray}
Denote the covariance of $X_1$ and $X_2$ as
\begin{equation}
	\theta = Cov(X_1,X_2) = E\left[ X_1X_2 \right] - E\left[X_1\right]E\left[X_2\right],
\end{equation}
and denote the correlated coefficient of $X_1$ and $X_2$ as
\begin{equation}\label{rho_exp}
	\rho = \frac{Cov(X_1,X_2)}{\sqrt{D(X_1)\cdot D(X_2)}}=\frac{p_{11}-p_{1+} p_{+1}}{\sqrt{p_{1+}(1-p_{1+}) p_{+1}(1-p_{+1}) }},
\end{equation}
it follows that
\begin{equation}
	\theta = \rho \sigma_1 \sigma_2.
\end{equation}

There are four parameters $ p_{00}$, $p_{01}$, $p_{10}$, $p_{11}$ in the bivariate binary distribution, with the constraint $p_{00} + p_{01} + p_{10} + p_{11} =1$, we need just three of them to describe the joint distribution.
In the following sections, another group of three parameters describes the bivariate binary distribution, including the two positive probabilities $p_{1+}$ and $p_{+1}$, and the correlated coefficient $\rho$.

\section{Bounds of Correlated Coefficient}

From the relationship between the correlated coefficient and the marginal positive probabilities (\ref{rho_exp}), we have
\begin{equation}
	p_{11} = p_{1+}p_{+1} + \rho  \sigma_1 \sigma_2,
\end{equation}
rewrite the two discordant probabilities as
\begin{eqnarray}
	p_{10} = p_{1+} - p_{11} = p_{1+} - p_{1+}p_{+1} - \rho \sigma_1 \sigma_2,  \\
	p_{01} = p_{+1} - p_{11} =  p_{+1}  - p_{1+}p_{+1} - \rho \sigma_1 \sigma_2.
\end{eqnarray}
Notice that $p_{10} > 0 $ and $p_{01} > 0$, we have the following two inequalities
\begin{eqnarray}
	&& \rho \sigma_1 \sigma_2 < p_{1+} - p_{1+}p_{+1},  \label{lower_1}\\
	&& \rho \sigma_1 \sigma_2 < p_{+1} - p_{1+}p_{+1},  \label{lower_2}
\end{eqnarray}
and the two inequalities (\ref{lower_1},\ref{lower_2})  should hold at the same time. If $p_{1+} > p_{+1}$ holds, the upper bound of $\rho$ satisfies
\begin{equation}
	\rho < \frac{p_{+1} - p_{1+}p_{+1}}{\sigma_1 \sigma_2} =  \sqrt{\frac{p_{+1}(1 - p_{1+})}{p_{1+}(1-p_{+1})}},
\end{equation}
and if $p_{1+} < p_{+1}$ holds, the upper bound of $\rho$ satisfies
\begin{equation}
	\rho < \frac{p_{1+} - p_{1+}p_{+1}}{\sigma_1 \sigma_2} = \sqrt{\frac{p_{1+}(1 - p_{+1})}{p_{+1}(1-p_{1+})}}.
\end{equation}

On the other hand, notice that $p_{10} < p_{1+}$ and $p_{01} < p_{0+} =  1 - p_{1+}$ should hold for the two discordant probabilities, therefore we have
\begin{eqnarray}
	&& p_{1+} - p_{1+}p_{+1} -\rho \sigma_1 \sigma_2  < p_{1+},\label{upper_1}\\
	&& p_{+1} - p_{1+}p_{+1} - \rho \sigma_1 \sigma_2 < 1 - p_{1+}, \label{upper_2}
\end{eqnarray}
where the two inequalities (\ref{upper_1},\ref{upper_2}) should hold at the same time. If $p_{1+} + p_{+1} - 1 > 0$ holds, the lower bound of $\rho$ satisfies
\begin{eqnarray}
	\rho \sigma_1 \sigma_2 & > &  p_{1+} + p_{+1} - 1 - p_{1+}p_{+1},   \\
	\rho    & > &  - \frac{(1 - p_{1+})(1 - p_{+1})}{\sqrt{p_{1+}p_{+1}(1-p_{1+})(1-p_{+1})}} \nonumber \\
	& = & - \sqrt{\frac{(1 - p_{1+})( 1 - p_{+1})}
		{p_{1+}p_{+1}}
	}.
\end{eqnarray}
On the other hand, if $p_{1+} + p_{+1} - 1 < 0$ holds, the lower bound of $\rho$ satisfies
\begin{eqnarray}
	\rho \sigma_1 \sigma_2 & > &   - p_{1+} p_{+1},   \\
	\rho & > & - \frac{p_{1+} p_{+1}} {\sqrt{p_{1+} p_{+1} (1 - p_{1+})( 1 - p_{+1})} } \nonumber \\
	& = &  - \sqrt{ \frac{p_{1+} p_{+1}}{(1 - p_{1+})(1 - p_{+1})} }.
\end{eqnarray}

The other constraints of the two discordant probabilities are $p_{10} < p_{+0} = 1 - p_{+1}$ and $p_{01} < p_{+1}$, which derive the same upper and lower bounds of $\rho$. Therefore we determine the upper and lower bounds of the correlated coefficient $\rho$ by the two positive probabilities $p_{1+}, p_{+1}$, which is the same as \cite{manuel_weib_silva_2014,wu_rai_yan_2019}.

\section{Parameter Space of the Null Hypothesis}

The purpose of the equivalence test is to determine whether the two positive probabilities $p_{1+}$ and $p_{+1}$ are equal to each other, which yields the following hypothesis test problem
\begin{equation}
	H_{0}: p_{1+} = p_{+1}; \qquad H_{1}:   p_{1+} \neq p_{+1}.
\end{equation}
Denote $p_{1+} = p_{+1} = \pi$, we have
\begin{equation}
	\sigma_1^2 = p_{1+}(1-p_{1+}) = \pi(1-\pi), \quad \sigma_2^2 = p_{+1}(1-p_{+1}) = \pi(1-\pi) \Rightarrow \sigma_1^2 =\sigma_2^2.
\end{equation}
When the null hypothesis holds, the two parameters $\pi$ and $\rho$ completely determine the parameter space.

\subsection{Bounds of Positive Probability}

When the null hypothesis $H_{0}$ holds, which means $p_{1+} = p_{+1} = \pi$, rewrite the upper bound of the correlated coefficient as
\begin{equation}
	\rho < \sqrt{\frac{p_{+1}(1-p_{1+})}{p_{1+}(1-p_{+1})}} = \sqrt{\frac{p_{1+}(1-p_{+1})}{p_{+1}(1-p_{1+})}} = \sqrt{\frac{\pi(1-\pi)}{\pi(1-\pi)}} = 1,
\end{equation}
which imposes no additional constraint on the upper bound of $\rho$.

On the other hand, when $p_{1+}+p_{+1}-1>0$ holds, which means $\pi > 1/2$, rewrite the lower bound of
$\rho$ as
\begin{equation}
	\rho > \frac{p_{1+} + p_{+1} -1 -p_{1+}p_{+1}}{\sigma_1 \sigma_2} = -\frac{1-\pi}{\pi},
\end{equation}
which yields the upper bound of positive probability $\pi$ in terms of the correlated coefficient $\rho$ as
\begin{equation}
	\pi < \frac{1}{1-\rho},
\end{equation}
where $\pi > 1/2$ holds.

When $p_{1+}+p_{+1}-1<0$ holds, which means $\pi < 1/2$, rewrite the lower bound of the correlated coefficient $\rho$ as
\begin{equation}
	\rho > -\frac{p_{1+}p_{+1}}{\sigma_1 \sigma_2} = -\frac{\pi}{1-\pi},
\end{equation}
which yields the lower bound of positive probability $\pi$ in terms of $\rho$ as
\begin{equation}
	\pi > - \frac{\rho}{1-\rho},
\end{equation}
where $\pi<1/2$ holds.

On the other hand, when the correlated coefficient $\rho > 0 $, we always have
\begin{equation}
	- \frac{\rho}{1-\rho} <  0  <  1  < \frac{1}{1-\rho}.
\end{equation}
When $\rho > 0$ holds, the range of positive probability $\pi$ should be
\begin{equation}
	0 < \pi < 1.
\end{equation}

On the other hand, the following inequality (\ref{ineq_rho_boundary}) always holds for a negatively correlated coefficient,
\begin{equation}\label{ineq_rho_boundary}
	0 < -\frac{\rho}{1-\rho} < \frac{1}{1-\rho} < 1,
\end{equation}
which yields the range of positive probability $\pi$ as
\begin{equation}
	-\frac{\rho}{1-\rho} <  \pi < \frac{1}{1-\rho},
\end{equation}
when $\rho<0$ holds. Therefore, we have the following proposition.
\begin{proposition}
	When the null hypothesis $H_0$ holds, the parameter space $\Theta_0(\rho, \pi)$ is determined by $\rho$ and $\pi$, which satisfies
	\begin{equation}\label{parameter_space_H_0}
		\Theta_0(\rho, \pi)
		=
		\left\{
		0< \pi < 1,   0 \leq \rho < 1
		\right\}
		\cup
		\left\{
		-\frac{\rho}{1-\rho} <  \pi < \frac{1}{1-\rho},  -1 < \rho < 0
		\right\}.
	\end{equation}
\end{proposition}

\subsection{Symmetry of Two Probabilities}

If the correlated coefficient $\rho \in (-1,0)$ holds, the summation of the upper and lower bounds of positive probability $\pi$ satisfies
\begin{equation}
	\frac{-\rho}{1 - \rho}  + \frac{1}{1 - \rho } = 1,
\end{equation}
which indicates that the center of the positive probability $\pi$ is at $\pi^* = 1/2$. On the other hand, the positive probability lies in $[0,1]$ when the correlated coefficient is positive, which indicates that the center of $\pi$ is also at $\pi^* = 1/2 $. Therefore we have the following proposition.
\begin{proposition}
	When the null hypothesis $H_{0}$ holds, the center of the positive probability $\pi$ is $\pi^* = \frac{1}{2}$.
\end{proposition}

When the null hypothesis $H_{0}$ holds, the two positive probabilities satisfy $p_{1+} = p_{+1} = \pi $, such that the consistent probability $p_{11}$ satisfies
\begin{equation}
	p_{11} = p_{1+} p_{+1} + \rho \sigma_1 \sigma_2 = \pi^2 + \rho \pi(1 - \pi).
\end{equation}
And the two discordant probabilities are equal to each other when the null hypothesis holds,
which derives $p_{10} = p_{01} = \pi - p_{11}$, thus we have
\begin{equation}
	p_{10} = p_{01} = \pi - p_{11} = \pi -  \left( \pi^2 + \rho \pi (1 - \pi) \right) = \pi (1 - \pi ) (1 - \rho).
\end{equation}
If $\rho < 0$ holds, $\pi$ reaches its minimum at the lower bound $\pi_{min} = -\frac{\rho}{1 - \rho}$, and reaches its maximum at the upper bound $\pi_{max} = \frac{1}{1 - \rho}$. Notice that the discordant probability satisfies
\begin{eqnarray}
	p_{10} (\pi_{min})= p_{01}(\pi_{min}) & = & -\frac{\rho}{1 - \rho} \left(1 + \frac{\rho}{1 - \rho}\right) (1 - \rho) \\
	&  =  &   -\frac{\rho}{1 - \rho},
\end{eqnarray}
and
\begin{eqnarray}
	p_{10} (\pi_{max})= p_{01}(\pi_{max})
	&  = & \frac{1}{1 - \rho} \left(1 - \frac{1}{1 - \rho}\right) (1 - \rho) \\
	&  =  &   -\frac{\rho}{1 - \rho} ,
\end{eqnarray}
which derives the following equations at $\pi_{min}$ and $\pi_{max}$
\begin{eqnarray}
	p_{10} (\pi_{max}) & = &  p_{10}(\pi_{min}), \\
	p_{01} (\pi_{max}) & = &  p_{01}(\pi_{min}).
\end{eqnarray}

Notice that the function $f(\pi) = \pi (1 - \pi)(1 - \rho)$ is symmetric about $\pi^* = \frac{1}{2}$ for any specified $\rho$,  therefore the discordant probability is symmetric about the positive probability $\pi$ for all $ \rho \in (0,1)$. And the discordant probability $p_{10}(\pi)$ reaches its maximum value at $\pi^* = \frac{1}{2}$, we have
\begin{equation}
	p_{10,max} = p_{01,max}= p_{10}\left(\pi^*=\frac{1}{2}\right) = \frac{1}{2} \left( 1 - \frac{1}{2} \right) \left(1 - \rho \right) = \frac{1 - \rho}{4}.
\end{equation}
Notice that the function $f(\rho) = \frac{1 - \rho}{4}$ decreases with increasing $\rho$, we have the following proposition.
\begin{proposition}
	When the null hypothesis $H_{0}$ holds, for any specified $\rho$, the discordant probability $p_{10}$ and $p_{01}$ are symmetric about $\pi^* = \frac{1}{2}$, and we have
	\begin{eqnarray}
		p_{10} & = &  p_{01} = \pi (1 - \pi)(1 - \rho), \\
		p_{10,max} & = & p_{01,max} = \frac{1 - \rho}{4}.
	\end{eqnarray}
	The range of discordant probability satisfies
	\begin{equation}\label{max_p_10}
		0 < p_{10} = p_{01}  \leq \frac{1}{2},
	\end{equation}
	where $p_{10}$ reaches its maximum value when the correlated coefficient $\rho = -1$.
\end{proposition}

\subsection{Monotonicity of Distribution Function}

The bivariate binary distribution can be viewed as the four terms multinomial distribution with observed frequency number $(N_{10},N_{01},N_{00},N_{11})$, where the discordant observations $N_{10}, N_{01}$ follow the binomial distribution with the successful probability $p_{10}$ and $p_{01}$ respectively. When $N_{10} = x_{10}$, the cumulative distribution function of $N_{10}$ is
\begin{equation}
	F(x_{10}; n, p_{10}) = Prob\left\{ N_{10} \leq x_{10} \right\} = \sum_{k=0}^{x_{10}} \frac{n!}{k!(n-k)!}p_{10}^{k}(1 - p_{10})^{n-k},
\end{equation}
where $n$ is the number of all four kinds of observations. For any fixed $x_{10}$ and $n$, the distribution function $F(x_{10};n, p_{10})$ is determined by $p_{10}$. Denote $x_{10} = x, p_{10}=p, C_n^k = \frac{n!}{k!(n-k)!}$ and $F(x_{10};n, p_{10}) = F(p)$, the first order derivation about $p$ is
\begin{eqnarray}\label{comulated_distribution_function}
	\frac{dF(p)}{dp} & = &  \sum_{k=0}^{x} C_n^k k p^{k-1} (1-p)^{n-k} - \sum_{k=0}^{x} C_{n}^k p^k (n-k)(1-p)^{n-k-1} \nonumber \\
	& = & \frac{1}{p} \sum_{k=0}^{x} k C_n^k p^{k} (1-p)^{n-k} -
	\frac{1}{1-p}  \sum_{k=0}^{x} (n-k) C_n^k p^{k} (1-p)^{n-k}.
\end{eqnarray}
Denote the random index function of set $A=\{ k: k \leq x \}$ as
\begin{equation}
	I_{A}(N=k) = \left\{
	\begin{array}{ll}
		1, & k \in A= \{k: k \leq x \};\\
		0, & k \notin A.
	\end{array}
	\right.
\end{equation}

On the other hand, when $N \sim B(n,p)$ holds, the conditional expectation of $N$ under the condition $\{N \leq x\}$ is
\begin{eqnarray}
	E[N | N \leq x ]     & = & \frac{1}{P\{A\}} E[N I_{A}(N)] = \frac{1}{P\{A\}} \sum_{k=0}^{n} k I_{A}(k) C_{n}^k p^k (1-p)^{n-k},   \\
	E[(n-N) | N \leq x ] & = & \frac{1}{P\{A\}} E[(n-N) I_{A}(N)] \nonumber \\
	& = & \frac{1}{P\{A\}} \sum_{k=0}^{n} (n-k) I_{A}(k) C_{n}^k p^k (1-p)^{n-k},
\end{eqnarray}
therefore (\ref{comulated_distribution_function}) can be expressed with conditional expectation as
\begin{eqnarray}
	\frac{dF(p)}{dp} &  = & \frac{P\{ A \}}{p} E[N | N \leq x] - \frac{P\{ A \}}{ 1 - p} E[(n-N) | N \leq x] \\
	& = &
	\frac{ P\{A\} }{p} \left\{  E[N | N \leq x]  - \frac{p}{1-p} E[(n-N) | N \leq x] \right\} \\
	& = &
	\frac{P\{A\}}{p} \left\{  E[N | N \leq x]  - \frac{p}{1-p} n + \frac{p}{1-p} E[N | N \leq x] \right\} \\
	& = &
	\frac{P\{A\}}{p(1-p)} \left\{  E[N | N \leq x]  - n p  \right\} \\
	& = &
	\frac{P\{A\}}{p(1-p)} \left\{  E[N | N \leq x]  - E[N]  \right\} \\
	& = &
	\frac{ P\{A\} }{ p(1-p) }
	\left\{  E[N | N \leq x]  - P\{ N \leq x \}   E[N | N \leq x] \right.\nonumber \\
	& & \left. -   (1 - P\{N \leq x \} ) E[N | N > x ]  \right\} \\
	& = &
	\frac{ P\{ A\} ( 1 - P\{N \leq x \} ) }{ p(1-p) } \left\{ E[N | N \leq x] - E[N | N > x]  \right\}. 	   \label{two_condition_probs}
\end{eqnarray}
Notice that
\begin{eqnarray}
	E[N | N \leq x] & \leq & x, \\
	E[N | N > x ]   & \geq    & x+1,
\end{eqnarray}
subtract to (\ref{two_condition_probs}), we have
\begin{equation}
	\frac{dF(p)}{dp} \leq   \frac{ P\{ A\} ( 1 - P\{ A \} ) }{ p(1-p) } \left[ x - (x+1) \right] \leq   - \frac{ P\{ A\} ( 1 - P\{A\} ) }{ p(1-p) } \leq 0,
\end{equation}
and the following proposition holds.
\begin{proposition}
	For any fixed number of trails $n$ and the observed frequency $x$, the cumulative distribution function of binomial distribution $F(x;n,p)$ decreases with increase of probability $p$, e.g. if $ p_{1} < p_{2} $, we have
	\begin{equation}
		F(x;n, p_{1}) \geq F(x; n, p_{2}).
	\end{equation}
\end{proposition}

Notice (\ref{max_p_10}), the distribution function $F(x_{10};n,p_{10})$ of $N_{10}$ satisfies the following proposition.
\begin{proposition}
	When the null hypothesis $H_{0}$ holds, for any $-1< \rho <1$ and $0< p_{10} < 1$, it follows
	\begin{equation}
		F(x_{10};n,p_{10}) \geq F(x_{10};n,\frac{1-\rho}{4}),
	\end{equation}
	which means that for any fixed $x_{10}$ and $n$, we have
	\begin{equation}
		\min_{p_{10}}F(x_{10};n,p_{10}) = F(x_{10};n,\frac{1}{2}).
	\end{equation}
\end{proposition}
The previous results state that if we consider the influence of the correlated coefficient $\rho$ between the two variables $(X_{1},X_{2})$ on the distribution function of discordant observation $N_{10}$, the distribution function $F(x_{10};n,p_{10})$ will reach its minimum value at the point $p_{10} = \frac{1}{2}$, which actually assumes that $\rho = -1$ holds.

\section{Confidence Region}

For any fixed sample size $n$ and discordant probabilities $p_{10}, p_{01}$, denote the observed discordant numbers $(N_{10}, N_{01})$ as $(x_{10}, x_{01})$, the sample space of the two discordant observations is
\begin{equation}
	\Omega = \left\{ (x_{10}, x_{01}) : 0 \leq x_{10} + x_{01} \leq n, 0 \leq x_{10}, x_{01} \leq n \right\}.
\end{equation}
Denote the significiant level as $1 - \alpha $, and the corresponding confidence region $A_{\alpha}$ as
\begin{equation} \label{trust_interval}
	A_{\alpha} = \left\{ (x_{10} , x_{01}):
	L_{10} \leq x_{10} \leq U_{10},
	L_{01}(x_{10}) \leq x_{01} \leq  U_{01}(x_{10})  \right\},
\end{equation}
where $L_{10},U_{10},L_{01}(x_{10}),U_{01}(x_{10})$ denote the corresponding lower and upper bounds respectively, thus the probability of $\left\{ \left(N_{10},N_{01}\right) \in A_{\alpha} \right\}$ is
\begin{eqnarray}
	& & Prob\left\{ \left(N_{10}, N_{01}\right) \in A_{\alpha} \right\}
	=
	\sum_{n_{10} = L_{10}}^{U_{10}}
	\sum_{n_{01} = L_{01}(n_{10})}^{U_{01}(n_{10})} Prob\left\{N_{10}=n_{10}, N_{01} = n_{01}\right\} \\
	& & =
	\sum_{n_{10} = L_{10}}^{U_{10}}
	\sum_{n_{01} = L_{01}(n_{10})}^{U_{01}(n_{10})} Prob\left\{N_{10}=n_{10}\right\}
	Prob\left\{N_{01} = n_{01} | N_{10}=n_{10} \right\}   \\
	& & =
	\sum_{n_{10} = L_{10}}^{U_{10}}  Prob\left\{N_{10}=n_{10}\right\}
	\sum_{n_{01} = L_{01}(n_{10})}^{U_{01}(n_{10})}
	Prob\left\{N_{01} = n_{01} | N_{10}=n_{10} \right\} \label{two_comb_prob}.
\end{eqnarray}
If the inner summation of the conditional probabilities is a constant  $\omega$, i.e.
\begin{equation}
	\sum_{n_{01} = L_{01}(n_{10})}^{U_{01}(n_{10})}
	Prob\left\{N_{01} = n_{01} | N_{10}=n_{10} \right\}  = \omega,
\end{equation}
and the outer summation of the unconditional probabilities is the same constant $\omega$, i.e.
\begin{equation}
	\sum_{n_{10} = L_{10}}^{U_{10}}  Prob\left\{N_{10}=n_{10}\right\}  = \omega,
\end{equation}
which makes (\ref{two_comb_prob}) can be rewritten as
\begin{equation}
	Prob\left\{ \left( N_{10}, N_{01} \right) \in A_{\alpha} \right\}
	= \omega^2 = 1 - \alpha,
\end{equation}
where  $\alpha$ is the significant level of the confidence region. Therefore we have
\begin{equation}
	\sum_{n_{10} = L_{10}}^{U_{10}}  Prob\left\{N_{10}=n_{10}\right\}  = \omega =\sqrt{1 - \alpha},
\end{equation}
and
\begin{equation}
	\sum_{n_{01} = L_{01}(n_{10})}^{U_{01}(n_{10})}
	Prob\left\{N_{01} = n_{01} | N_{10}=n_{10} \right\}  = \omega =  \sqrt{1 - \alpha}.
\end{equation}

To determine the acceptance region $A_{\alpha}$, we just need to determine the lower bound $L_{10}$ and upper bound  $U_{10}$ for the first discordant variable $N_{10}$, and to determine the corresponding conditional upper and lower bound  $U_{01}(n_{10})$ and $L_{01}(n_{10})$ of the second discordant variable $N_{01}$, which satisfy the following constraints
\begin{eqnarray}
	L_{10} &=& \arg\max_{L} Prob\left\{ N_{10} < L \right\}  \leq   \frac{1- \sqrt{1 - \alpha}}{2}, \label{margin_test_1}\\
	U_{10} &=& \arg\min_{U} Prob\left\{ N_{10} < U \right\}  \geq  1 - \frac{1- \sqrt{1 - \alpha}}{2}, \label{margin_test_2}\\
	L_{01}(n_{10}) &=& \arg\max_{L}	Prob\left\{ N_{01} < L  | N_{10} = n_{10}\right\}  \leq  \frac{1- \sqrt{1 - \alpha}}{2}, \label{margin_test_3}\\
	U_{01}(n_{10}) &=& \arg\min_{U}	Prob\left\{ N_{01} < U | N_{10} = n_{10}\right\}  \geq   1 - \frac{1- \sqrt{1 - \alpha}}{2}. \label{margin_test_4}
\end{eqnarray}
We construct the confidence region $A_{\alpha}$ of the two-dimensional variable $(N_{10},N_{01})$
through their marginal distribution and conditional distribution.

When the null hypothesis $H_0$ holds,  i.e. $p_{10} = p_{01} = p$, the maximum likelihood estimation of the discordant probability $p_{10}$ and $p_{01}$ can choose to be
\begin{equation}
	\hat{p}_{10} = \hat{p}_{01}= \frac{x_{10} + x_{01}}{2n}.
\end{equation}

The frequency $N_{10}$ and $N_{01}$ both follow the binomial distribution $B(n,p)$ respectively. For a fixed sample size $n$ and significant level $\alpha$, denote the left and right endpoints
as $L_{\alpha}$ and $U_{\alpha}$ respectively, which satisfy the following two conditions
\begin{eqnarray}
	& L_{\alpha}   &= \arg\max_{L}Prob\{ N < L \} \leq \frac{1 - \sqrt{1 - \alpha}}{2}, \\
	& U_{\alpha}   &= \arg\min_{U}Prob\{ N < U \} \geq  1 - \frac{1 - \sqrt{1 - \alpha}}{2},
\end{eqnarray}
where $N \sim B(n,\hat{p}_{10})$. Therefore the rejection region for the null hypothesis $H_{0}$ can be
\begin{equation}\label{rejection_margin_test}
	R_{\alpha} = \left\{ \min\left\{x_{10}, x_{01} \right\} < L_{\alpha} \right\} \cup  \left\{\max\left\{x_{10}, x_{01} \right\} > U_{\alpha} \right\}.
\end{equation}
Because the marginal distribution of the two discordant observations $(N_{10}, N_{01})$ is used to construct the rejection region, the proposed method is called the margin test with significant level $\alpha$ thereinafter.

\section{McNemar test}

The most frequently used equivalence test method is the McNemar test \cite{McNemar_test_1947}, and the statistic is
\begin{equation}\label{McNemar}
	M = \frac{\left(x_{10} - x_{01}\right)^2}{x_{10} + x_{01}} \sim \chi^2(1).
\end{equation}
The corresponding rejection region of McNemar test is
\begin{eqnarray}\label{McNemar_rejection}
	R_{\alpha} & = &  \left\{ M > \chi^2_{1-\alpha}(1)\right\}, \\
	& = & \left\{\left(x_{10},x_{01}\right) \in \Omega:  \frac{\left(x_{10} - x_{01}\right)^2}{x_{10} + x_{01}} > \chi^2_{1 - \alpha}(1)   \right\},
\end{eqnarray}
where $\chi^2_{1 - \alpha}(1)$ is the $1 - \alpha$ left quantile of the $\chi^2$-distribution with $1$  degree of freedom.
Although there are many improved versions of the McNemar test, the traditional version (\ref{McNemar}) is considered in this paper, as the focus of this paper is to investigate the influence of the correlated coefficient on the result of the equivalence test.

\section{Evaluation of Tests}

The size and power are the technical indexes to evaluate the performance of the equivalence tests in this paper. When the null hypothesis $H_0$ holds, the two positive probabilities equal to each other, i.e. $p_{1+} = p_{+1} = \pi$, and the possible range of parameter value lies in the space $\Theta_0$, which is determined in (\ref{parameter_space_H_0}). Denote the parameter space of $H_{1}$ as
\begin{equation}\label{parameter_space_H_1}
	\Theta_{1} = \left\{ (p_{10}, p_{01}) :  0 < p_{10} + p_{01} < 1, 0 < p_{10}, p_{01} < 1, p_{10} \neq p_{01} \right\}.
\end{equation}

Denote the rejection region as $R_{\alpha}$, and the corresponding acceptance region is $A_{\alpha}$. The size of test is the probability of making the type I error. For a pair of parameter $(\rho, \pi)\in \Theta_{0}$, the corresponding size is
\begin{equation}\label{size_definition}
	size(\rho,\pi) = Prob\{ (N_{10}, N_{01}) \in R_{\alpha} | (\rho, \pi) \in \Theta_0 \}.
\end{equation}
For any specified rejection region $R_{\alpha}$, (\ref{size_definition}) determines the value of size for each point $(\rho,\pi)$ in the parameter space $\Theta_0$, which yields the surface of size on the parameter space $\Theta_0$.  Taking the advantage of the two-variable model for the parameter space $\Theta_0$, it is possible to intuitively evaluate the performance of the test method through the explicit surface of size in the following sections. The surface of size directly shows the influence of the correlated coefficient on the result of the equivalence test.

On the other hand, the power of test method is the probability of correctly rejecting the null hypothesis for the parameters in the parameter space $\Theta_{1}$, i.e.
\begin{equation}\label{power_definition}
	power(p_{10},p_{01}) = Prob\{ (N_{10}, N_{01}) \in R_{\alpha} | (p_{10}, p_{01}) \in \Theta_{1}\}.
\end{equation}
For a specified rejection region $R_{\alpha}$, calculate the probability following (\ref{power_definition}) for each point $(p_{10},p_{01})$ in the parameter space of the alternative hypothesis $H_1$, which yields the surface of power on the set $\Theta_{1}$.

\subsection{Acceptance Region and Rejection Region}

For a fixed sample size $n$, the McNemar test and the margin test can determine whether a point $\left(x_{10},x_{01}\right)$ in the sample space $\Omega$ belongs to the acceptance region $A_{\alpha }$, or the region of rejection $R_{\alpha}$ respectively. The sample size determines the range of sample space $\Omega$, and the level of significance determines the shape of the acceptance region.

\begin{figure}[h!]
	\begin{center}
		\includegraphics[width = 1.00 \textwidth]{acceptance.eps}
		\caption{Boundaries of Acceptance Regions for $\alpha = 0.05, n=1600, 800, 400, 200, 100, 50$.}
		\label{fig:A_alpha}
	\end{center}
\end{figure}

Under the significant level $\alpha = 0.05$, Fig.\ref{fig:A_alpha} shows the boundaries of acceptance regions of the McNemar test and the margin test for sample sizes $n = 1600,800,400,200,100,50$. The dotted red lines are the boundaries of the McNemar test, and the solid blue lines are the boundaries of the margin test. The shape of the two acceptance regions is similar, and the acceptance region of the margin test always wraps the corresponding domain of the McNemar test. Therefore the probability that a sample follows in the acceptance region of the margin test is slightly larger than the corresponding probability of the McNemar test.

\subsection{Size of Tests}

Denote the acceptance region of $H_{0}$ as
\begin{equation}
	A_{\alpha} = \left\{ (x_{10}, x_{01}): L_{10} \leq x_{10} \leq U_{10},
	L_{01}(x_{10}) \leq x_{01} \leq U_{01}(x_{10}) \right\},
\end{equation}
for any specified value of parameter $(\rho, \pi) \in \Theta_{0}$, the probability that a sample $(N_{10},N_{01})$ lies in the rejection region $R_{\alpha}$ is
\begin{eqnarray}
	size(\rho,\pi) & = &  Prob\{ (N_{10}, N_{01}) \in R_{\alpha} | (\rho,\pi) \in \Theta_{0}\} \nonumber \\
	& = & 1 - Prob\{ (N_{10}, N_{01}) \in A_{\alpha} | (\rho,\pi) \in \Theta_{0}\} \nonumber \\
	& = & 1 - \sum_{n_{10}=L_{10}}^{U_{10}} Prob\left\{ N_{10} = n_{10} \right\}
	\sum_{n_{01}=L_{01}(n_{10})}^{U_{01}(n_{10})} Prob\left\{ N_{01} = n_{01} | N_{10} = n_{10} \right\},
\end{eqnarray}
where $N_{10} \sim B(n,p_{10})$,  $N_{01} | N_{10} = n_{10} \sim B(n - n_{10}, q_{01})$, and the conditional probability is
\begin{equation}
	q_{01} = \frac{p_{01}} {1 - p_{10}}.
\end{equation}

If the resulting size is below the specified value of significant level $\alpha$, the target of controlling the type I error is achieved. Otherwise, the value of size exceeds the significant level, which means the test does not satisfy the desired target of risk control.

\begin{figure}[h!]
	\begin{center}
		\includegraphics[width = 1.00 \textwidth]{size_large_n.eps}
		\caption{Contours of the Size Surfaces for $\alpha = 0.05, n=1600,800,400$.}
		\label{fig:size_large_n}
	\end{center}
\end{figure}

A two-variable model describes the parameter space of $H_0$ in this paper, which makes it possible to see the influence of the correlated coefficient on the value of size. Fig. \ref{fig:size_large_n} and Fig. \ref{fig:size_small_n} show the contours of size surfaces, where the horizontal coordinate indicates the value of correlated coefficient $\rho$, and the vertical coordinate indicates the value of positive probability $\pi$.
When the correlated coefficient changes from $-0.99$ to $0.99$, the type I error shows a downward tendency, and the size is symmetric about $\pi^* = 1/2$. Note that the significant level $\alpha$ is $0.05$, those points covered by the contours below $0.05$ are the correctly accepted parameters.
The acceptance region of the margin test wraps up the acceptance domain of the McNemar test, such that the range of correctly identifying parameters by the margin test is larger than the corresponding region of the McNemar test.

\begin{figure}[h!]
	\begin{center}
		\includegraphics[width = 1.00 \textwidth]{size_small_n.eps}
		\caption{Contours of the Size Surfaces for $\alpha = 0.05, n=200, 100, 50$.}
		\label{fig:size_small_n}
	\end{center}
\end{figure}

For the regions of positively correlated coefficient, the margin test can correctly identify almost all of the parameters in $\Theta_0$, while the McNemar test can only correctly identify a part of the parameters when the sample size $n$ is $100$ or $50$. For the negative correlation coefficients, the two methods can't correctly identify all of the parameters in $\Theta_0$, and the correctly identified region of the margin test is obviously larger than the corresponding domain of the McNemar test. From the Fig.\ref{fig:size_large_n} and Fig. \ref{fig:size_small_n}, we can also find that the size of both tests gradually increases with the increment of sample size, and the risk of the type I error for the McNemar test increases significantly than the increment of the margin test.

We can also find that the surface of size given by the McNemar test fluctuates up and down with the increase of correlation coefficient in three experiments with sample sizes $n=1600,800,400$. The contour line of $0.05$ repeatedly appears three times on the right side of the graph, which indicates that the risk of type I error for the McNemar test fluctuates around $\alpha$ with relatively larger $\rho$.
We also find that the fluctuation of the size surface with the increase of correlation coefficient also appears for a relatively small sample size with a large value of $\alpha$.
However, the surface of size given by the margin test has no such obvious fluctuation, which means that the probability of making the type I error decreases monotonously with the increase of the correlation coefficient. It can conclude that the margin test is more stable for different sample sizes and significance levels.

\subsection{Power of Tests}

For any specified rejection region $R_{\alpha}$ and acceptance region $A_{\alpha} = \Omega - R_{\alpha}$, denote the power of each parameter point $(p_{10},p_{01}) \in \Theta_{1}$ as
\begin{eqnarray}
	& & power(p_{10},p_{01})  =  Prob\{ (N_{10}, N_{01}) \in R_{\alpha} | (p_{10},p_{01}) \in \Theta_{1} \} \nonumber \\
	& & =  1 -  Prob\{ (N_{10}, N_{01}) \in A_{\alpha} | (p_{10},p_{01}) \in \Theta_{1} \}  \nonumber \\
	&  & =   1 - \sum_{n_{10} = L_{10}}^{U_{10}} Prob\{N_{10} = n_{10}\}
	\sum_{n_{01} = L_{01}(n_{10})}^{U_{01}(n_{10})}  Prob\{ N_{01} = n_{01} | N_{10} = n_{10}\},
\end{eqnarray}
where the conditional distribution of $N_{01}$ under the condition $N_{10} = n_{10}$ is a binomial distribution with
\begin{equation}
	N_{01} | N_{10} = n_{10} \sim B(n - n_{10}, q_{01}),
\end{equation}
where the number of trails is $n-n_{10}$, and the successful probability $q_{01}$ satisfies
\begin{equation}
	q_{01} = \frac{p_{01}}{1 - p_{10}}.
\end{equation}

\begin{figure}[h!]
	\begin{center}
		\includegraphics[width = 1.00 \textwidth]{power_large_n.eps}
		\caption{Contours of Power Surfaces with $\alpha = 0.05, n=1600,800,400$.}
		\label{fig:power_large_n}
	\end{center}
\end{figure}

The contours of power surfaces of the McNemar test and the margin test are given in Fig.\ref{fig:power_large_n} and Fig.\ref{fig:power_small_n}, where the level of significant is $0.05$, and the samples sizes $n = 1600,800,400,200,100,50$. Another figure of the contours of two power surfaces is also given in Fig.\ref{fig:power_n_50}  for $\alpha = 0.35, n=50$. The horizontal coordinate indicates the value of $p_{10}$, and the vertical coordinate indicates the value of $p_{01}$ for the parameter $(p_{10},p_{01}) \in \Theta_{1}$.

\begin{figure}[h!]
	\begin{center}
		\includegraphics[width = 1.00 \textwidth]{power_small_n.eps}
		\caption{Contours of Power Surfaces with $\alpha = 0.05, n=200,100,50$.}
		\label{fig:power_small_n}
	\end{center}
\end{figure}

From Fig. \ref{fig:power_large_n} and Fig. \ref{fig:power_small_n}, we can find that the region where the McNemar test and the margin test correctly reject the null hypothesis for the parameters in $\Theta_1$ gradually increases with the increase of sample size, and the region where the two tests accept the null hypothesis always concentrates around the diagonal line direction.

\begin{figure}[h!]
	\begin{center}
		\includegraphics[width = 1.00 \textwidth]{power_n_50_a_35.eps}
		\caption{Contours of Two Power Surfaces with $\alpha = 0.35, n = 50$.}
		\label{fig:power_n_50}
	\end{center}
\end{figure}

On the other hand, the probability of correctly accepting $H_0$ by the McNemar test tends to decrease with the reduction of $p_{10}$ and $p_{01}$. To observe the tendency of reducing, we give the Fig. \ref{fig:power_n_50}, where the significant level is relatively large, $\alpha = 0.35$, and the sample size is relatively small, $n = 50$. We can find that the region enclosed by the contour of $power = 0.35$ concentrates around the diagonal line. And the area given by the McNemar test separates from the original point, while the region given by the margin test closely connects to the original point. The observed connection shows that the margin test tends to correctly accept the null hypothesis for small $p_{10}$ and $p_{01}$, while the McNemar test tends to incorrectly reject the null hypothesis for small discordant probabilities.

\section{Real-World Examples}

To compare the performance of the McNemar test and the margin test with real-world data, especially for small sample sizes and those observed data on the boundary of the acceptance region, we discuss the following two real-world examples.

\subsection{Stem Cell Transplantation and Airway Hypertension}

In the first example, we discuss a small sample case. Twenty-one children underwent airway hypertension testing before and after stem cell transplantation (SCT) \cite{fagerland_lydersen_laak_2013_sct,bentur_lapidot_livnat_2009_SCT_example}, the observed data are shown in Table \ref{tab:21_SCT},
where $n = 21, n_{10} = 7, n_{01} = 1$, AHR means positive response, and NO AHR means negative response.

\begin{table}
	\begin{center}
		\caption{Airway Hypertension Test for 21 Children before and after Stem Cell Transplantation.}
		\label{tab:21_SCT}
		\begin{tabular}{ccccc}
			\toprule
			&  &\multicolumn{2}{c}{After SCT} &  \\
			\cline{3-4}
			&  & AHR    & NO AHR & SUM \\
			\hline
			& AHR    &1  & 1     &2 \\
			Before SCT	& NO AHR            &7  & 12     & 19 \\
			& SUM              &8 & 13   & 21 \\
			\bottomrule
		\end{tabular}
	\end{center}
\end{table}

\begin{table}[h]
	\begin{center}
		\caption{Results of McNemar Test and Margin Test with Sample Disturbance for the Airway Hypertension Data, $\alpha = 0.05$.}
		\label{tab:21_SCT_sample_disturbance}
		\begin{tabular}{llccccc}
			\toprule
			&		&	$n$	&	$n_{10}$	&	$n_{01}$	&	p-value	&	Decision	\\
			\hline
			Original Data	&	McNemar	&	21	&	7	&	1	&	0.0339	&	$H_{1}$	\\
			&	Margin	&	21	&	7	&	1	&	0.1033	&	$H_{0}$	\\
			Add One	&	McNemar	&	22	&	7	&	2	&	0.0596	&	$H_{0}$	\\
			&	Margin	&	22	&	7	&	2	&	0.2057	&	$H_{0}$	\\
			Reduce One	&	McNemar	&	20	&	6	&	1	&	0.0588	&	$H_{0}$	\\
			&	Margin	&	20	&	6	&	1	&	0.1588	&	$H_{0}$	\\
			Adjust One	&	McNemar	&	21	&	6	&	2	&	0.1573	&	$H_{0}$	\\
			&	Margin	&	21	&	6	&	2	&	0.2957	&	$H_{0}$	\\
			\bottomrule
		\end{tabular}
	\end{center}
\end{table}

The results of the equivalence test are presented in Table \ref{tab:21_SCT_sample_disturbance}, where the McNemar test rejects the null hypothesis under the significant level of $\alpha = 0.05$, which means that there is a significant difference between the probabilities of positive response before and after SCT. However, the margin test accepts the null hypothesis with the same level of significance. We propose the following sample disturbance method to  decide whether  to accept the null hypothesis or not, assuming that $n_{01} < n_{10}$ and there is a unit sample disturbance in favor of $H_{0}$:
\begin{enumerate}
	\item[(1)] Add one sample, which increases the sample size by one unit, adding a new observation such that $n = n+1, n_{10} = n_{10}, n_{01} = n_{01 } + 1$;
	\item[(2)] Reduce one sample, which removes an existing observation that is unfavorable to $H_0$ such that $n = n - 1, n_{10} = n_{10} -1, n_{01} = n_{01 }$;
	\item[(3)] Adjust one sample, which maintains the sample size and reduces the difference between the two discordant observations $n = n, n_{10} = n_{10} - 1, n_{01} = n_{01} + 1$.
\end{enumerate}

Since reducing the difference between $n_{10}$ and $n_{01}$ is a change in favor of $H_{0}$, the above three adjustment schemes make the new points move towards the acceptance region. Table \ref{tab:21_SCT_sample_disturbance} gives the results of the two tests with the four kinds of data, where the McNemar test rejects the null hypothesis and the margin test accepts $H_0$ with the original data. However, both tests accept the null hypothesis with the three adjusted observations, which means that the original data is on the boundary of the rejection region of the McNemar test, and the most closed three points to the original data fall into the acceptance region of the McNemar test. On the other hand, the original and deduced data are all in the acceptance region of the margin test. The difference between two discordant observations is not so significant, such that one unit of sample disturbance can lead to a different decision of the McNemar test. Based on the original data, we thus recommend accepting the null hypothesis $H_{0}$ as the difference between the two discordant observations is not sufficiently significant. To make a more accurate judgment, it is also necessary to increase the sample size to keep the observed data far from the boundary of the rejection region.

\begin{table}
	\begin{center}
		\caption{Frequency of Successfully Opening Box B by Male and Female Chimpanzees.}
		\label{tab:open_box}
		\begin{tabular}{ccccc}
			\toprule
			&  &\multicolumn{2}{c}{Gender} &  \\
			\cline{3-4}
			&  & Male   & Fmale & SUM \\
			\hline
			& NO    & 12  & 9     &21 \\
			Open Box B	& YES               &27  & 17     & 44 \\
			& SUM              &39 & 26   & 65 \\
			\bottomrule
		\end{tabular}
	\end{center}
\end{table}

\subsection{Box Opening and Gender}

\cite{smith_ruxton_2020} presented a behavioral ecology example, where $39$ male and $26$ female chimpanzees participated in the experiment to see if they could successfully open box B to get food among two puzzle boxes.
Table \ref{tab:open_box} presents the number of male and female chimpanzees who open box B.

\begin{table}[h]
	\begin{center}
		\caption{Results of McNemar test and Margin test with Sample Disturbance of Chimpanzee Data, $\alpha = 0.05$.}
		\label{tab:monkeys_sample_disturbance}
		\begin{tabular}{llccccc}
			\toprule
			&		&	$n$	&	$n_{10}$	&	$n_{01}$	&	p-value	& Decision	\\
			\hline
			Original Data	&	McNemar	&	65	&	27	&	9	&	0.0027	&	$H_{1}$	\\
			&	Margin	&	65	&	27	&	9	&	0.0120	&	$H_{1}$	\\
			Add One	&	McNemar	&	66	&	27	&	10	&	0.0052	&	$H_{1}$	\\
			&	Margin	&	66	&	27	&	10	&	0.0193	&	$H_{1}$	\\
			Reduce One	&	McNemar	&	64	&	26	&	9	&	0.0041	&	$H_{1}$	\\
			&	Margin	&	64	&	26	&	9	&	0.0165	&	$H_{1}$	\\
			Adjust One	&	McNemar	&	65	&	26	&	10	&	0.0077	&	$H_{1}$	\\
			&	Margin	&	65	&	26	&	10	&	0.0262	&	$H_{0}$	\\
			\bottomrule
		\end{tabular}
	\end{center}
\end{table}

Table \ref{tab:monkeys_sample_disturbance} presents the results of the McNemar test and the margin test with four kinds of data, where the results are the same for the original data, the added data, and the reduced data. And a different result of the two tests is from the adjusted data. It shows that the adjusted sample point $(26,10)$ falls in the acceptance region of the margin test, while the other three points $(27,9)$, $(27,10)$, and $(26,9)$ are all in the rejection regions of the two tests. Therefore the original data $(27,9)$ is actually in the rejection region of both tests, and we can reject the null hypothesis and accept that there is a significant correlation between the ability to open box B and the gender of chimpanzees.

\section{Discussion and Conclusion}

The bivariate binary distribution is a statistical model frequently used in biomedical statistics. Four joint probabilities $(p_{00}, p_{01}, p_{10}, p_{11})$ determine the structure of the model. And the correlation coefficient $\rho$ and two positive probabilities $p_{1+}, p_{+1}$ can describe the same model from the viewpoint of the marginal distribution. When the equivalence assumption holds, the two marginal probabilities are equivalent, such that the correlation coefficient $\rho$ and the positive probability $\pi$ can also determine the structure of the model.

According to the relationship between the joint and the marginal distribution, we determine the upper and lower bounds of the correlation coefficient $\rho$ and find the range of the marginal probability $\pi$ for a fixed correlation coefficient.
Based on the two-variable parameter space of the null hypothesis, we discuss the symmetry of the marginal probability $\pi$ and the discordance probability $p_{10}$.
Using the conditional expectation method, we give the monotonicity of the cumulative distribution function of the discordant observation. And the minimum value of the distribution function is also determined by the correlation coefficient.
Using the joint distribution of two discordant observations, we propose an estimation algorithm for the confidence region of the two discordance observations. And an alternative equivalence test is derived from the resulting confidence region, called the margin test.

The McNear test is frequently used to test the equivalence hypothesis of the bivariate binary model, which assumes that the two dimensions are independent.
To investigate the influence of the correlation between the two dimensions in the bivariate binary model on the result of the equivalence test,  we compare the acceptance region, size, and power of the McNear test and the margin test with different sample sizes and significance levels.
The shape of the two acceptance regions is very similar, and the acceptance region of the margin test is slightly larger than the region of the McNemar test for the same sample size and significant level.
The range that the two tests correctly accept the null hypothesis gradually narrows with the increase in sample size. And the region that the margin test correctly identifies the parameters belonging to the parameter space of the null hypothesis is significantly larger than the corresponding region of the McNemar test for the same sample size, which means that the risk of type I error of the margin test is significantly smaller than the risk of the McNemar test.
And it is found that the risk of type I error of the McNemar test is not monotonically decreasing with the increase of the correlation coefficient $\rho$ for relatively large sample size, while the type I error of the margin test does monotonically decrease with the increase of the correlation coefficient for all of the sample sizes.
From the contours of the power surface, we can find that the margin test tends to accept the null hypothesis when the discordant probability is relatively small, while the McNemar test tends to reject the null hypothesis with the same parameters.

There are two real-world examples in the previous section, one of the observed points lies on the boundary of the rejection region of the McNemar test, and the other point lies on the edge of the rejection region of the margin test.
A unit sample disturbance method determines whether to accept the null hypothesis when the sample point lies between the boundaries of the two rejection regions.

We use the joint distribution of two discordant observations to construct the confidence region. The correlation coefficient is not directly used in the test process, so the influence of the correlation coefficient on the result of the equivalence test can only be observed when the null hypothesis holds. It is necessary to develop new statistical techniques to handle the observations with relatively small discordant probability or negatively correlated coefficients.

\bibliographystyle{plain}


\section{Introduction}

Enzymes are nanometer-size biomolecular complexes that catalyze biochemical reactions and play a crucial role in various life-sustaining processes.
For example, they facilitate molecular transport into the cell and assist chemical reactions that are essential for cellular metabolism and homeostasis~\cite{albertsbook}.
In recent years, nonequilibrium transport phenomena induced by enzymes have attracted a considerable attention due to novel applications such as drug delivery and the design of synthetic nanomachines~\cite{gompper2020}.

In aqueous solutions of enzymes, enhanced diffusion was experimentally reported for certain enzymes~\cite{riedel2015,illien2017}, as well as for passive probe particles~\cite{zhao2017} in the presence of \textit{substrate} molecules.
Here substrates are chemical species that react with enzymes and are converted into product molecules.
When a spacial gradient in the substrate molecule concentration exists, enzymes exhibit a collective motion in the direction of higher or lower concentrations.
This phenomenon is known as chemotaxis~\cite{sengupta2013,sengupta2014,yu2009,dey2014} or antichemotaxis, respectively~\cite{jee2018}.
Although a molecular diffusion enhancement was observed also for organic chemical reactions at much smaller length scales (subnanometer)~\cite{dey2016,wang2020}, the observed diffusion enhancement is still a matter of debate~\cite{macdonald2019,rezaei2022} and more experiments are needed.

To better understand enzyme mobility, various theoretical studies using coarse-grained models have been developed.
In previous works, an enzyme was modeled as built of two sub-units representing its overall structure.
The two units are connected by an elastic spring reflecting inherent enzyme relaxation dynamics~\cite{togashi2007}.
By using equilibrium approaches, it was demonstrated that the internal degrees of freedom leads to enhanced diffusion of the enzyme in solution~\cite{illien2017_epl,hosaka2020}.
On the other hand, in nonequilibrium situations, enzyme conformational dynamics
collectively induces hydrodynamic flows~\cite{hosaka2020_2,hosaka2017,mikhailov2015}.
Furthermore, employing a force dipole model for the enzyme, it was found that such nonequilibrium effects lead to an increase in particle diffusion in solutions as well as for biological membranes~\cite{hosaka2020_2,hosaka2017,mikhailov2015}.
More recently, hydrodynamic interactions between force dipoles were taken into account, and clustering mechanisms of force dipoles were investigated both in flat~\cite{manikantan2020} or curved~\cite{bagaria2022} membrane geometries.

These enzymatic nonequilibrium effects were studied only in passive fluids, whereas at physiological conditions, other sources of nonequilibrium effects can turn the surrounding fluids into active fluids.
For example, membrane proteins such as ATPase can autonomously rotate in the presence of ATP (adenosine triphosphate) or proton gradients, and their induced hydrodynamic flows drive the surrounding membrane into an out-of-equilibrium state~\cite{oppenheimer2019,manneville1999,manneville2001}.
More specifically, ATP consumption and autonomous rotation of membrane proteins can lead to breaking of time-reversal and parity symmetries in membranes~\cite{banerjee2017}.
These ATP consumption and autonomous rotation endow membranes with active and chiral features, respectively.
At length scales larger than the mean distance between rotary proteins, membranes with these proteins can be viewed as two-dimensional (2D) active chiral fluids.
In such active chiral 2D systems, it is known that a dissipationless transport coefficient called the \textit{odd viscosity} emerges~\cite{hosaka2022,fruchart2022}.
To reveal the hydrodynamic effects of odd viscosity, its consequence has been studied for the motion of passive objects~\cite{ganeshan2017,souslov2020,hosaka2021,hosaka2021_2,lier2022}, many-body sedimentation~\cite{khain2022}, and density waves~\cite{banerjee2017,markovich2021}.
One of the peculiar features of active chiral fluids is their nonreciprocal interaction~\cite{hosaka2021,hosaka2021_2,khain2022,lier2022} that is prohibited for passive fluids.
However, despite these intriguing findings, the dynamics of active enzymes in an active chiral environment was not studied, and the role of odd viscosity in biomembranes remains unexplored.

In this paper, we discuss the lateral dynamics of an enzyme pair that interacts via hydrodynamic interactions as it is embedded in a thin layer of an active chiral fluid.
The enzymes are modeled as active particles that induce force dipoles~\cite{mikhailov2015,manikantan2020}.
To investigate collective behavior of \textit{active} particles in an \textit{active} environment, we consider force dipoles in a chiral fluid layer that is modeled as a 2D compressible fluid with odd viscosity.
Since such fluids have a nonreciprocal nature, active particles are expected to exhibit chiral trajectories although the particles themselves are apolar and do no have any preferred direction.

Extending our previous work~\cite{hosaka2021}, we derive the odd viscosity-dependent mobility tensor that includes nonreciprocal hydrodynamic interactions.
As a minimum model to explore hydrodynamic interactions between active particles, we consider pairs of two active particles and discuss the nonlinear two-body dynamics.
In contrast to the passive fluid case without odd viscosity, we find that active particle pairs show spiral trajectories where one particle follows the other.
In the limit of infinitely large odd viscosity, particles show two types of periodic oscillations, including one which is determined by the ratio between the magnitude of the odd viscosity and force dipole.
Our findings reveal that the nonreciprocal response due to odd viscosity leads to complex dynamics of active particles.

The outline of the manuscript is as follows.
In Sec.~\ref{sec:system}, we introduce the hydrodynamic equations for a 2D active chiral fluid and derive its mobility tensor~\cite{hosaka2021}.
In Sec.~\ref{sec:twobody}, we obtain the nonlinear equations for the distance between the two active particles and their relative polar angles~\cite{manikantan2020}.
Using these nonlinear equations, we analyze in Sec.~\ref{sec:result} the dynamics of particle pairs for various values of odd viscosity and classify their behavior into several characteristic states.
Section~\ref{sec:discussion} includes further discussions, and a summary and conclusions are presented in Sec.~\ref{sec:conc}.

\begin{figure}[htb]
\centering
\includegraphics[scale=.3]{./Fig1.eps}
\caption{
Schematic sketch of an active and chiral 2D fluid layer.
The infinitely large, flat, and thin 2D fluid layer (light blue) is located at $z=h$ and has 2D shear, dilatational, and odd viscosities, $\eta_{\rm s}$, $\eta_{\rm d}$, and $\eta_{\rm o}$, respectively.
This fluid layer is in contact with air $(z>h)$ and a 3D fluid (dark blue) underneath $(0<z<h)$ characterized by a 3D shear viscosity $\eta$.
The 3D fluid is bounded from below by an impermeable flat solid surface (brown) located at $z=0$, and the 3D velocity is assumed to vanish at $z=0$.
The active particle that represents an enzyme (orange disk) has radius $a$ and can move laterally within the 2D fluid layer.
The dipolar flow that is induced by particle ``1'' (``2'') is characterized by its direction $\widehat{\mathbf{d}}_1$ $(\widehat{\mathbf{d}}_2)$ and the relative angle $\theta_1$ $(\theta_2)$ as defined in Eq.~(\ref{eq:cos}).
}
\label{fig:system}
\end{figure}

\section{Active chiral 2D layer}
\label{sec:system}

\subsection{Hydrodynamic equations}

We consider an active, chiral, and compressible 2D layer, which is flat, very thin, and infinitely large, as schematically depicted in Fig.~\ref{fig:system}.
The layer is in contact with a 3D incompressible fluid (e.g., water), which is bound from below by a rigid substrate.
This leads to a momentum leakage from the 2D layer to the 3D fluid.
We assume that the 3D fluid thickness $h$ is much smaller than any in-plane characteristic length scales so that
the lubrication approximation will hold for the 3D fluid~\cite{barentin1999,elfring2016}.
The 2D compressible fluid layer can be realized experimentally, for example, by a dilute Gibbs monolayer composed of soluble amphiphiles that can dissolve into the underlying 3D fluid~\cite{barentin1999}.
Such a 2D/3D material transport makes the 2D fluid compressible.
We assume that the adsorption and desorption processes of soluble amphiphiles are instantaneous, and can be regarded as a limiting case of a finite relaxation time~\cite{lier2022}.

The active chiral nature of the 2D fluid is taken into account by introducing the concept of odd viscosity.
Although the regular viscosity (such as shear and dilatational viscosities) is always positive, odd viscosity can be either positive or negative, depending on the chirality direction.
At microscopic scales, the 2D odd viscosity sign is related to the rotational direction of the active constituents such as rotary proteins in biological membranes~\cite{markovich2021,banerjee2017}.
To see clearly the effect of odd viscosity on the active particle dynamics, we consider a compressible 2D fluid~\cite{hosaka2021}.
The reason being that the odd viscosity does not play a role in the velocity field of an incompressible fluid~\cite{ganeshan2017}.
Note that in the incompressible layer limit, the fluid flow becomes independent of the odd viscosity, as will be shown below.

We denote the 2D velocity field by $\mathbf{v}(\mathbf{r})$ with $\mathbf{r}=(x,y)$ being a positional vector in 2D, and $p$ is the hydrostatic pressure of the 3D fluid.
At low Reynolds numbers, the momentum balance equation for the 2D fluid can be written as~\cite{hosaka2021}
\begin{align}
\eta_{\rm s}\nabla^2\mathbf{v} +\eta_{\rm d}\nabla (\nabla\cdot\mathbf{v}) + \eta_{\rm o}\nabla^2 \mathbf{v}^\ast
 -\frac{h}{2}\nabla p
+\mathbf{f}^{\rm 3D}
+\mathbf{F}=0.
\label{eq:stokes}
\end{align}
Here, $\eta_{\rm s}$, $\eta_{\rm d}$, and $\eta_{\rm o}$ are the 2D shear, dilatational, and odd viscosities, respectively, $\nabla=(\partial_x, \partial_y)$ stands for the 2D gradient operator, $\mathbf{f}^{\rm 3D}$ is the vectorial force density exerted on the 2D fluid layer by the underlying 3D fluid, $\mathbf{F}$ is any other force density acting on the 2D fluid, $v_i^\ast=\epsilon_{ij}v_j$ is the velocity vector rotated clockwise by $\pi/2$, and $\epsilon_{ij}$ is the 2D Levi-Civita tensor with $\epsilon_{xx}=\epsilon_{yy}=0$ and $\epsilon_{xy}=-\epsilon_{yx}=1$.
The divergence of the in-plane velocity is given by~\cite{barentin1999,elfring2016}
\begin{align}
\nabla\cdot\mathbf{v} = \frac{h^2}{6\eta}\nabla^2 p,
\label{eq:compress}
\end{align}
where $\eta$ is the shear viscosity of the underlying 3D fluid.

In general, the force density $\mathbf{f}^{\rm 3D}$ caused by the 3D fluid can be expressed as $f_i^{\rm 3D}=-(\zeta_\|\delta_{ij}+\zeta_\perp\epsilon_{ij})v_j$~\cite{soni2019}, where $\zeta_\|$ is the friction coefficient representing the momentum dissipation parallel to $\mathbf{v}$, while $\zeta_\perp$ acts perpendicular to $\mathbf{v}$, and $\delta_{ij}$ is the Kronecker delta.
For a 2D fluid layer supported by a rigid substrate, the parallel friction coefficient is given by $\zeta_\|\simeq
\eta/h$~\cite{evans1988}, in the limit where $h$ is small enough.
In biomembranes, rotary motor proteins protrude out of the membrane plane~\cite{oppenheimer2019}, implying that the 2D odd viscosity $\eta_{\rm o}$ would contribute to the transverse component of the friction coefficient $\zeta_\perp$.
In order to obtain the mobility tensor, we assume an identical hydrodynamic screening length both for the longitudinal $\mathbf{v}$ and transverse $\mathbf{v}^\ast$ velocities.
Hence, we take the transverse friction to be proportional to the ratio between the odd and shear viscosities, $\zeta_\perp=(\eta_{\rm o}/\eta_{\rm s})\zeta_\|$.
This leads to a force density given by
\begin{align}
f_i^{\rm 3D}=-\frac{\eta}{h}\left(\delta_{ij}+\frac{\eta_{\rm o}}{\eta_{\rm s}}\epsilon_{ij}\right)v_j.
\label{eq:f3d}
\end{align}

\subsection{Mobility tensor of an active chiral layer}

The force density $\mathbf{F}$ acting on the fluid layer at position $\mathbf{r}^\prime$ is connected via a second-rank mobility tensor $\mathbf{G}(\mathbf{r})$ with the induced fluid velocity at position $\mathbf{r}$:
\begin{align}
v_{i}(\mathbf{r})=\int {\rm d}^{2} r^{\prime}\, G_{i j}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) F_{j}\left(\mathbf{r}^{\prime}\right).
\end{align}
By solving the coupled hydrodynamic equations~(\ref{eq:stokes}) and (\ref{eq:compress}) in Fourier space, we show in Appendix~\ref{app:GFT} that the mobility tensor $G_{ij}[\mathbf{k}]$ has the following
expression:
\begin{align}
&G_{ij}[\mathbf{k}] =\nonumber\\
&\frac{(\eta_{\rm s}\widehat{k}_i\widehat{k}_j-\eta_{\rm o}\epsilon_{i j})(k^2+\kappa^2)
+(\eta_{\rm s}+\eta_{\rm d})(k^2+\lambda^2)\overline{k}_i\overline{k}_j
}
{\eta_{\rm s}(\eta_{\rm s}+\eta_{\rm d})\left(k^{2}+\kappa^{2}\right)\left(k^{2}+\lambda^{2}\right)+\eta_{\rm o}^2\left(k^{2}+\kappa^{2}\right)^2},
\label{eq:gmobilityk}
\end{align}
where $\mathbf{k} = (k_x ,k_y)$ is the 2D wavevector, $k=|\mathbf{k}|$, $\widehat{k}_i=k_i/k$, and $\overline{k}_i=-\epsilon_{ij}\widehat{k}_j$.
Taking the ratio of the 2D viscosities $\eta_{\rm s}$ and $\eta_{\rm d}$ with the 3D shear viscosity $\eta$ defines two hydrodynamic screening lengths~\cite{barentin1999}
\begin{align}
\kappa^{-1}=\left(\frac{\eta_{\rm s}h}{\eta}\right)^{1/2},~~~~~
\lambda^{-1}=\left[\frac{(\eta_{\rm s}+\eta_{\rm d})h}{4\eta} \right]^{1/2}.
\end{align}
Due to the transverse momentum decay $\zeta_\perp$, the screening length $\kappa^{-1}$ also appears in the term $\eta_{\rm o}(k^2+\kappa^2)\epsilon_{ij}$ of Eq.~(\ref{eq:gmobilityk}).
This term was not considered in our previous work~\cite{hosaka2021}.

In the special 2D incompressible limit $\eta_{\rm d}\to\infty$, Eq.~(\ref{eq:gmobilityk}) reduces to the mobility tensor for an incompressible supported 2D fluid~\cite{ramachandran2011,oppenheimer2010}
\begin{align}
G_{ij}^0[\mathbf{k}] =
\frac{\delta_{ij}-\widehat{k}_i\widehat{k}_j}
{\eta_{\rm s}\left(k^{2}+\kappa^{2}\right)}.
\end{align}
Note that $\mathbf{G}^0$ in this limit does not depend on the odd viscosity $\eta_{\rm o}$~\cite{ganeshan2017,khain2022,hosaka2021,souslov2020}.

For simplicity sake, we assume that the two screening lengths are equal $\kappa=\lambda$ (or equivalently, $\eta_{\rm d}=3\eta_{\rm s}$).
Then, Eq.~(\ref{eq:gmobilityk}) can be simplified as
\begin{align}
G_{ij}[\mathbf{k}] =\frac{4\delta_{ij}-3\widehat{k}_i\widehat{k}_j- \mu\epsilon_{i j}}{\eta_{\rm s}(4+\mu^2)\left(k^{2}+\kappa^{2}\right)},
\label{eq:mobilityk}
\end{align}
where $\mu=\eta_{\rm o}/\eta_{\rm s}$ and the relation $\overline{k}_i\overline{k}_j=\delta_{ij}-\widehat{k}_i\widehat{k}_j$ has been used.
The dimensionless parameter $\mu$ is a measure of how far the 2D active chiral fluid departs from its passive analog, due to the active constituents that self-spin at microscopic scales.

The real-space representation of the mobility tensor can be obtained by the inverse Fourier transform of Eq.~(\ref{eq:mobilityk}).
The derivation is shown in Appendix~\ref{app:G} and the final result is
\begin{align}
G_{i j}(\mathbf{r})=C_{1}(r) \delta_{i j}+C_{2}(r) \widehat{r}_{i} \widehat{r}_{j}+C_{3}(r) \epsilon_{i j},
\label{eq:G}
\end{align}
where $\widehat{\mathbf{r}}=\mathbf{r}/r$ is a unit vector $(r=|\mathbf{r}|)$ and the three position-dependent coefficients are
\begin{align}
C_1(r) &= \frac{1}{2\pi\eta_{\rm s}(4+\mu^2)} \left[-\frac{3}{(\kappa r)^2}+4K_0(\kappa r)+\frac{3K_1(\kappa r)}{\kappa r}\right],
\nonumber\\
C_2(r) &= \frac{3}{2\pi\eta_{\rm s}(4+\mu^2)} \left[\frac{2}{(\kappa r)^2}-K_0(\kappa r)-\frac{2K_1(\kappa r)}{\kappa r}\right],
\label{eq:c123}\\
C_3(r) &=  -\frac{\mu K_0(\kappa r)}{2\pi\eta_{\rm s}(4+\mu^2)}
,\nonumber
\end{align}
and $K_n(x)$ is the modified Bessel function of the second kind~\cite{abramowitzhandbook}.
Note that $C_3$ exists only when $\mu\neq0$ (nonzero odd viscosity).

Expanding the mobility tensor, Eq.~(\ref{eq:G}), for $\kappa r\ll1$, we obtain
\begin{align}
G_{i j}(\mathbf{r}) &\approx \frac{1}{8\pi\eta_{\rm s}(4+\mu^2)}
\left[
\left(-3-10\gamma+ 10\ln\frac{2}{\kappa r}\right)\delta_{ij}\right.\nonumber\\
&\left.+6 \widehat{r}_{i} \widehat{r}_{j}
+4\mu\left(\gamma - \ln\frac{2}{\kappa r}\right)\epsilon_{ij}
\right].
\label{eq:Gsmall}
\end{align}
In this limit, the mobility tensor depends logarithmically on the distance $r$.
In the opposite limit of $\kappa r\gg1$, we obtain
\begin{align}
&G_{i j}(\mathbf{r})\approx  \frac{1}{2\pi\eta_{\rm s}(4+\mu^2)} \nonumber\\
&\times\left[
\frac{3}{(\kappa r)^2} \left(-\delta_{ij}+2\widehat{r}_{i} \widehat{r}_{j}\right)
-\mu
\sqrt{\frac{\pi}{2\kappa r}}e^{-\kappa r}\epsilon_{ij}
\right],
\label{eq:Glarge}
\end{align}
where $\gamma\approx0.5772$ is Euler's constant.
The term due to the regular viscosities in $\mathbf{G}$ decays algebraically ${\sim}(\kappa r)^{-2}$, while the term due to odd viscosity decreases exponentially ${\sim}e^{-\kappa r}$.

In the following, we concentrate on the regime $\kappa r\ll1$ in order to investigate the effect of odd viscosity on the collective dynamics of active particles.
The $\kappa r\ll1$ limit is justified when the product of the 2D shear viscosity $\eta_{\rm s}$ and the 3D fluid thickness $h$ is much larger than the 3D viscosity $\eta$, as occurring in physiological conditions.
Using typical values such as $\eta_{\rm s}\approx\SI[parse-numbers=false]{10^{-9}}{\pascal\second\metre}$, $h\approx\SI[parse-numbers=false]{10^{-8}}{\metre}$, and $\eta\approx\SI[parse-numbers=false]{10^{-3}}{\pascal\second}$~\cite{saffman1976,ramachandran2010}, we find $\kappa^{-1}\approx\SI[parse-numbers=false]{10^{-7}}{\metre}$.
When active particles move across length scales $\SI{1}{\nano\metre}\leq r\leq\SI{10}{\nano\metre}$, we have $10^{-2}\leq\kappa r\leq10^{-1}$ and the condition $\kappa r\ll1$ is satisfied.

\section{Two-body hydrodynamic interactions}
\label{sec:twobody}

\subsection{Hydrodynamic force dipole}

In the 2D active chiral compressible fluid introduced above, we consider two active particles of radius $a$ that interact hydrodynamically with each other, as depicted in Fig.~\ref{fig:system}.
A single active particle can represent an enzymatic molecule.
Since enzymes generate dipolar flows while changing their conformations~\cite{mikhailov2015}, they can be modeled as active particles that induce hydrodynamic force dipole~\cite{manikantan2020,hosaka2017}.
Without any hydrodynamic interactions between other active particles, a single particle with a force dipole does not show any motility because the force-free condition is imposed on each particle.

\begin{figure}[htb]
\centering
\includegraphics[scale=.7]{Fig2.eps}
\caption{
Streamlines of the 2D velocity $\mathbf{v}(x,y)$ generated by a force dipole.
The force dipole along the $x$-axis [$\widehat{\mathbf{d}}=(1,0)$] is centered at the origin (the black double arrow) for (a) $\mu=0$ (no odd viscosity, $\eta_{\rm o}=0)$ and (b) $\mu=1$ [see Eq.~(\ref{eq:velocity})].
The blue arrows indicate the flow direction.
}
\label{fig:dipole}
\end{figure}

When an active particle with a force dipole resides at the origin and its dipole is directed along a given unit vector $\widehat{\mathbf{d}}$, it induces a velocity field given by~\cite{mikhailov2015}
\begin{align}
v_{i}(\mathbf{r})=-\sigma \widehat{d}_{k} \partial_{k} G_{i j}(\mathbf{r}) \widehat{d}_{j}.
\label{eq:dipole}
\end{align}
Here, $\sigma=fa$ is the magnitude of the force dipole where $f$ is the force magnitude and $a\ll r$ is the distance between the two point forces (the particle size).

\begin{figure*}[htb]
\centering
\includegraphics[scale=.45]{Fig3.eps}
\caption{
Plots of the rescaled distance between the two particles, $r=|\mathbf{r}_{12}|$, as a function of the rescaled time $t$ for $\eta_{\rm o}=0$.
(a) Convergence state (CS, cross in Fig.~\ref{fig:phase}) for $(\theta_1,\theta_2)=(\pi/2,\pi/2)$.
The black circle corresponds to the collision time between the two particles.
(b) Monotonic divergence state (MDS, closed triangle in Fig.~\ref{fig:phase}) for $(\theta_1,\theta_2)=(\pi/4,\pi/4)$.
(c) Nonmonotonic divergence state (NDS, open triangle in Fig.~\ref{fig:phase}) for $(\theta_1,\theta_2)=(\pi/2,\pi/4)$.
}
\label{fig:mu0}
\end{figure*}

Returning to the particle pair, the vector connecting particle ``1'' to particle ``2'' is denoted as $\mathbf{r}_{12}$, and the two relative angles $\theta_1,\theta_2$ are defined as
\begin{align}
\cos\theta_1 = \frac{\widehat{\mathbf{d}}_1\cdot\mathbf{r}_{12}}{|\mathbf{r}_{12}|},~~~~~
\cos\theta_2 = \frac{\widehat{\mathbf{d}}_2\cdot\mathbf{r}_{21}}{|\mathbf{r}_{21}|},
\label{eq:cos}
\end{align}
where $\widehat{\mathbf{d}}_1$ $(\widehat{\mathbf{d}}_2)$ is the dipolar direction of particle $1$ $(2)$.
Substituting Eq.~(\ref{eq:Gsmall}) into Eq.~(\ref{eq:dipole}), we obtain the translational velocities of particle $2$ relative to particle $1$
\begin{align}
&\mathbf{v}_{21}(\mathbf{r})= \frac{\sigma_1}{4\pi\eta_{\rm s}(4+\mu^2)r}
\nonumber\\
&\times\left[
3\cos(2\theta_1)\widehat{\mathbf{r}}_{12}+2\cos\theta_1\widehat{\mathbf{d}}_1-2\mu\cos\theta_1\widehat{\mathbf{d}}_1^\ast
\right],
\label{eq:velocity}
\end{align}
with $\widehat{d}_i^\ast=\epsilon_{ij} \widehat{d}_j$.
The derivation of Eq.~(\ref{eq:velocity}) is shown in Appendix~\ref{app:relative}.

The velocity field that is induced by a force dipole is plotted in Figs.~\ref{fig:dipole}(a) and (b) for $\mu=0$ $(\eta_{\rm o}=0)$ and $\mu=1$, respectively.
Here, the force dipole is located at the origin and is directed along the $x$-direction with $\widehat{\mathbf{d}}_1=(1,0)$.
When $\mu=0$ [as in Fig.~\ref{fig:dipole}(a)], the flow is symmetric with respect to both the $x$- and $y$- axes, and the azimuthal component $\widehat{\mathbf{d}}_1$ causes surrounding fluids to flow away from the origin.
This leads to an outward flow around the force dipole.
On the other hand, when $\mu=1$ [as in Fig.~\ref{fig:dipole}(b)], one sees that the emerging flow lines are tilted along the diagonal line $x=y$.
This symmetry breaking in the $x$- and $y$- axes is due to the perpendicular contribution $\widehat{\mathbf{d}}^\ast_1$ in Eq.~(\ref{eq:velocity}) when $\mu\neq0$.

\subsection{Pair dynamics of active particles}

We proceed by examining the hydrodynamic interactions between a pair of active particles, as depicted in Fig.~\ref{fig:system}.
Each particle moves laterally and rotates in the 2D fluid as a result of the fluid velocity induced by the other particle.
The separation between the two particles evolves as $\dot{r}=dr/dt=(\mathbf{v}_{21}-\mathbf{v}_{12})\cdot\widehat{\mathbf{r}}$~\cite{manikantan2020}, where $\mathbf{v}_{21}$ and $\mathbf{v}_{12}$ are the translational velocities of particle $2$ relative to $1$ [see Eq.~(\ref{eq:velocity})] and vice versa, and $\mathbf{r}=\mathbf{r}_{12}$.
When $\sigma_1=\sigma_2=\sigma$, the distance between the two particles evolves as
\begin{align}
&\dot{r}= \frac{\sigma}{4\pi\eta_{\rm s}(4+\mu^2)r} \nonumber\\
&\times
\biggl[2+4\cos(2\theta_1)+4\cos(2\theta_2)
-\mu\left[\sin(2\theta_1)+\sin(2\theta_2)\right] \biggl].
\label{eq:r}
\end{align}

Each particle rotates at a rate equal to half the vorticity of the flow velocity induced by the other particle: $\Omega_1=(\nabla\times\mathbf{v}_{12})_z/2$ and $\Omega_2=(\nabla\times\mathbf{v}_{21})_z/2$~\cite{manikantan2020}, where $\Omega_1$ and $\Omega_2$ are the $z$-components of vortices for the 2D rotation.
An additional effect of rotating $\mathbf{r}$ emerges when the induced velocity has an azimuthal component~\cite{manikantan2020}.
In total, the two angles evolve as
\begin{align}
\begin{aligned}
\dot{\theta}_1 &=-\frac{\sigma}{4\pi\eta_{\rm s}(4+\mu^2)r^2}[4\sin(2\theta_2)+\mu\cos(2\theta_2)]+\dot{\phi},
\\
\dot{\theta}_2 &=-\frac{\sigma}{4\pi\eta_{\rm s}(4+\mu^2)r^2}[4\sin(2\theta_1)+\mu\cos(2\theta_1)]+\dot{\phi},
\label{eq:theta12}
\end{aligned}
\end{align}
where $\dot{\phi}$ is an angular rotation due to the azimuthal components $\widehat{\mathbf{d}}$ and $\widehat{\mathbf{d}}^\ast$ of the velocity in Eq.~(\ref{eq:velocity}):
\begin{align}
&\dot{\phi}= -\frac{\sigma}{4\pi\eta_{\rm s}(4+\mu^2)r^2}\nonumber\\
&\times
\biggl[ \sin(2\theta_1)+\sin(2\theta_2)+\mu \left[2+\cos(2\theta_1)+\cos(2\theta_2)\right] \biggl].
\label{eq:phi}
\end{align}
More details about the derivation of Eq.~(\ref{eq:theta12}) are given in Appendix~\ref{app:relative}.
From Eqs.~(\ref{eq:theta12}) and (\ref{eq:phi}), one can see that the angular evolution is determined not only by the angle of the other particle, but also by its own angle because there is a coupling between $\theta_1$ and $\theta_2$, Eq.~(\ref{eq:phi}).
The asymptotic expressions of Eqs.~(\ref{eq:r}) and (\ref{eq:theta12}) for the two limits of $\mu\ll1$ and $\mu\gg1$ are also given in Appendix~\ref{app:relative}.

\begin{figure*}[htb]
\centering
\includegraphics[scale=.55]{Fig4.eps}
\caption{
Pair dynamics of two active particles (denoted as red and green) initially located at $(x, y)=(0,0)$ and $(4a,0)$.
The arrows indicate the trajectory direction.
(a) Trajectories of a particle pair whose initial orientation is $(\theta_1,\theta_2)=(3\pi/4, \pi/2)$ for $\mu=3$ (SP1, orange diamond in Fig.~\ref{fig:phase}).
(b) Trajectories of a particle pair whose initial orientation is $(\theta_1,\theta_2)=(\pi/2,\pi/16)$ for $\mu=3$ (SP2, blue inverted triangle in Fig.~\ref{fig:phase}).
(c) Trajectories of a particle pair whose initial orientation is $(\theta_1,\theta_2)=(3\pi/4,\pi/2)$ for $\mu=10$.
}
\label{fig:mu3}
\end{figure*}

\section{Results for numerical simulations}
\label{sec:result}

We investigate the dynamics of two active particles by solving the nonlinear equations~(\ref{eq:r})-(\ref{eq:phi}) for $r$, $\theta_1$, and $\theta_2$.
First, we analyze the active particle trajectories for the three cases ($\eta_{\rm o}=0$, $\eta_{\rm o}\neq0$, and $\eta_{\rm o}\to\infty$), and classify their behaviors into several characteristic states.
Then, we perform numerical simulations with different initial relative angles, and plot the phase diagram for the two angles while keeping the initial inter-particle distance $r$ constant.
To understand the mechanism of the particle collective behavior, we plot the phase space for $\theta_1$ and $\theta_2$ using Eq.~(\ref{eq:theta12}), and perform a linear stability analysis to obtain the characteristic time scale for the pair dynamics.

\subsection{Collective behavior of the active particle pair}
\label{sec:collective}

When $\eta_{\rm o}=0$ (or equivalently, $\mu=0$), we examine the typical temporal evolution of the pair inter-particle distance $r$ for various values of the initial angles, $\theta_1$ and $\theta_2$, while keeping the initial distance fixed at $r=4a$.
When $(\theta_1,\theta_2)=(\pi/2,\pi/2)$, Fig.~\ref{fig:mu0}(a) shows that the two particles approach each other and collide at $r=a$.
For other values of $(\theta_1,\theta_2)$, they move away monotonically and nonmonotonically as shown in Fig.~\ref{fig:mu0}(b) [$(\theta_1,\theta_2)=(\pi/4,\pi/4)$] and Fig.~\ref{fig:mu0}(c) [$(\theta_1,\theta_2)=(\pi/2,\pi/4)$].
We classify these three behaviors as ``convergence state'' (CS), ``monotonic divergence state'' (MDS), and ``nonmonotonic divergence state'' (NDS), shown respectively in Figs.~\ref{fig:mu0}(a), (b), and (c).

Finite values of the odd viscosity $(\mu\neq0)$ result in a perpendicular (nonreciprocal) component $\widehat{\mathbf{d}}^\ast_1$ in the induced velocity, as in Eq.~(\ref{eq:velocity}).
Consequently, it leads to the ``spiral state 1'' (SP1) when $(\theta_1,\theta_2)=(3\pi/4,\pi/2)$, as shown in Fig.~\ref{fig:mu3}(a).
For the angles $(\theta_1,\theta_2)=(\pi/2,\pi/16)$ [Fig.~\ref{fig:mu3}(b)], the active particles show orbiting behavior following attractive or repulsive trajectories before the spiral behavior, which we call the ``spiral state 2'' (SP2).
As the odd viscosity increases, spiral patterns become more evident as shown in Fig.~\ref{fig:mu3}(c) for $\mu=10$.

\begin{figure*}[htb]
\centering
\includegraphics[scale=.55]{Fig5.eps}
\caption{
Pair dynamics of active particles in the limit of $\mu\gg1$.
(a) The particles are initially located at $(x, y)=(0,0)$ and $(4a,0)$ with dipolar orientations $\theta_1=0.64\pi$ and $\theta_2=0.35\pi$, respectively.
Snapshots of the particle dynamics at (b) $t/(|\eta_{\rm o}|a^2/\sigma)=10^3$ and (c) $t/(|\eta_{\rm o}|a^2/\sigma)=2\times10^3$.
The pair dynamics corresponds to the circular oscillation state (COS, green circle) in phase space, Fig.~\ref{fig:phase}(c).
}
\label{fig:mu_infinity}
\end{figure*}

\begin{figure}[htb]
\centering
\includegraphics[scale=.5]{Fig6.eps}
\caption{
Inter-particle distance $r/a$ of a particle pair (black line) and the two dipolar orientations $\theta_1,\theta_2$ (red dashed and green dotted lines), respectively as a function of the dimensionless time $t/(|\eta_{\rm o}|a^2/\sigma)$.
The initial values are $r/a=4$, $\theta_1=2$, and $\theta_2=1.1$.
}
\label{fig:mu_infinity_rel}
\end{figure}

\begin{figure*}[htb]
\centering
\includegraphics[scale=.55]{Fig7.eps}
\caption{
Pair dynamics of active particles in the limit of $\mu\gg1$.
(a) The particles are initially located at $(x, y)=(0,0)$ and $(4a,0)$ with dipolar orientations $\theta_1=3\pi/4$ and $\theta_2=\pi/4$, respectively.
Snapshots of the particle dynamics at (b) $t/(|\eta_{\rm o}|a^2/\sigma)=10^3$ and (c) $t/(|\eta_{\rm o}|a^2/\sigma)=4\times10^3$.
The pair dynamics corresponds to the radial oscillation state (ROS, red star) in the phase space in Fig.~\ref{fig:phase}(c).
}
\label{fig:mu_infinity_star}
\end{figure*}

\begin{figure}[htb]
\centering
\includegraphics[scale=.5]{Fig8.eps}
\caption{
Inter-particle distance $r/a$ of a particle pair (black line), and the two dipolar orientations $\theta_1,\theta_2$ (red dashed and green dotted lines) as a function of the dimensionless time $t/(|\eta_{\rm o}|a^2/\sigma)$.
The initial values are $r/a=4$, $\theta_1=3\pi/4$, and $\theta_2=\pi/4$.
}
\label{fig:mu_infinity_star_rel}
\end{figure}

In the limit of $\mu\gg1$, two types of oscillatory behaviors are seen and they depend on the angles.
Their trajectories and relative dynamics are presented in Figs.~\ref{fig:mu_infinity}, \ref{fig:mu_infinity_rel}, \ref{fig:mu_infinity_star}, and \ref{fig:mu_infinity_star_rel}.
For the dipolar angles $(\theta_1,\theta_2)=(0.64\pi,0.35\pi)$, Fig.~\ref{fig:mu_infinity} shows that one particle follows the other, and they make a circular orbit around a common point located at $(x,y)\approx(2a,4a)$.
Figure~\ref{fig:mu_infinity_rel} shows that the distance and the angles periodically oscillate in time.
This means that the active particles do not remain too close or too far away from each other along the circular trajectory.
For $(\theta_1,\theta_2)=(3\pi/4, \pi/4)$, as shown in Fig.~\ref{fig:mu_infinity_star}, the particles manifest a reciprocal rotation that leads to radial trajectories.
Figure~\ref{fig:mu_infinity_star_rel} exhibits an oscillation period in $r$, $\theta_1$, and $\theta_2$, which is smaller than the period for the circular oscillation shown in Fig.~\ref{fig:mu_infinity_rel}.
We coin the oscillatory behaviors in Figs.~\ref{fig:mu_infinity} and \ref{fig:mu_infinity_star} as the ``circular oscillation state'' (COS) and ``radial oscillation state'' (ROS), respectively.

\subsection{State diagram in the $(\theta_1,\theta_2)$ plane}

To examine the dependence of the active particle dynamics on the polar angles, we numerically integrate Eqs.~(\ref{eq:r}) and (\ref{eq:theta12}) with $13$ different initial values of $\theta_1$ for $\theta_2$.
Figures~\ref{fig:phase}(a), (b), and (c) show the state diagrams in the $(\theta_1,\theta_2)$ plane for $\mu=0$, $3$, and $\mu\gg1$, respectively.
State diagrams are constructed by varying initial angles and categorizing corresponding states of the active particle pair, as shown in Sec.~\ref{sec:collective}.
When $\mu=0$, Fig.~\ref{fig:phase}(a) shows that the particle pair exhibits a ``monotonic divergence state'' (MDS, closed triangle) along the diagonal ($\theta_1=\theta_2$), except for $(\theta_1,\theta_2)=(\pi/2,\pi/2)$.
However, the particle pair along the other diagonal ($\theta_1=-\theta_2$) exhibits a ``convergence state'' (CS, cross).
Note that around the points of $(\theta_1,\theta_2)=(\pi/2,0)$, $(0,\pi/2)$, $(\pi/2,\pi)$, and $(\pi,\pi/2)$, the particles show a MDS behavior (closed triangle), whereas the ``nonmonotonic divergence state'' (NDS, open triangle) is observed for most of the other points in the $(\theta_1,\theta_2)$ plane.

When the odd viscosity is finite, as shown in Fig.~\ref{fig:phase}(b), the ``spiral state 1'' (SP1, orange diamond) and ``spiral state 2'' (SP2, blue inverted triangle) emerge, and the active particle pair shows the two types of spiral behavior, as in Figs.~\ref{fig:mu3}(a) and (b).
In the limit of $\mu\gg1$, as shown in Fig.~\ref{fig:phase}(c), the spiral states disappear altogether and the characteristic circular and radial oscillation states appear in turn.
One can see that around the symmetrical point, $(\theta_1,\theta_2)=(\pi/2,\pi/2)$, the ``circular oscillation state'' (COS, green circle) dominates, while the particles show the ``radial oscillation state'' (ROS, red star) in other points of the $(\theta_1,\theta_2)$ plane.

\begin{figure*}[htb]
\centering
\includegraphics[scale=.42]{Fig9.eps}
\caption{
State diagram in the $(\theta_1,\theta_2)$ plane for (a) $\mu=0$, (b) $\mu=3$, and (c) $\mu\gg1$.
Different symbols represent different pair dynamics: monotonic divergence state (MDS, closed triangle), nonmonotonic divergence state (NDS, open triangle), convergence state (CS, cross), spiral state 1 (SP1, orange diamond), spiral state 2 (SP2, blue inverted triangle), circular oscillation state (COS, green circle), and radial oscillation state (ROS, red star).
The phase behavior of $(\theta_1,\theta_2)$ for a fixed value of $r$ for (d) $\mu=0$ [see Eqs.~(\ref{eq:theta1s}) and (\ref{eq:theta2s})], (e) $\mu=3$ [see Eq.~(\ref{eq:theta12})], and (f) $\mu\gg1$ [see Eqs.~(\ref{eq:theta1inf}) and (\ref{eq:theta2inf})].
Shaded regions indicate that the particle pair approaches each other, i.e., when $\dot{r}<0$.
}
\label{fig:phase}
\end{figure*}

\subsection{Phase space in dynamical systems}

By connecting the state diagrams with the dynamics of the particle angles, we can analyze more quantitatively the behavior of the active particle pair in the $(\theta_1,\theta_2)$ plane.
The phase space of $\theta_1$ and $\theta_2$ are plotted in Figs.~\ref{fig:phase}(d), (e), and (f) for $\mu=0$, $3$, and $\mu\gg1$, respectively.
The shaded regions in those figure parts indicate that the active particles approach each other ($\dot{r}<0$).
The phase plane for $\mu=0$ is shown in Fig.~\ref{fig:phase}(d), where the vector field ($\dot{\theta}_1,\dot{\theta}_2$) is calculated from Eqs.~(\ref{eq:theta1s}) and (\ref{eq:theta2s}) in Appendix~\ref{app:relative}, while keeping $r$ const.
One can see that there are closed orbits whose centers are located at $(\theta_1,\theta_2)=(\pi/2,0)$, $(0,\pi/2)$, $(\pi/2,\pi)$, and $(\pi,\pi/2)$~\cite{strogatz}.
The point $(\pi/2,\pi/2)$ is unstable along the diagonal line ($\theta_1=\theta_2$), leading to trajectories that move away from that point.
On the other hand, the trajectory starting on the opposite diagonal line ($\theta_1=-\theta_2$) moves towards the $(\pi/2,\pi/2)$ point~\cite{strogatz}.
This point is called a \textit{hyperbolic saddle-point}~\cite{strogatz,manikantan2020}.
Any perturbation of these diagonal points kicks a system into a closed orbit around one of four centers at $(\pi/2,0)$, $(0,\pi/2)$, $(\pi/2,\pi)$, or $(\pi,\pi/2)$ (fixed points)~\cite{strogatz,manikantan2020}.

Using linear stability analysis, we can calculate the time scale of the nonlinear oscillations.
We assume that the angles $\theta_1$ and $\theta_2$ evolve exponentially as $e^{\alpha t}$ when the active particles at a separation $d$ are perturbed from one of the fixed points.
The obtained eigenvalues $\alpha=\pm \sqrt{6}\sigma i/(4\pi\eta_{\rm s}d^2)$ are purely imaginary, leading to an oscillation period $T_0=2\pi/|\alpha|=8\pi^2\eta_{\rm s}d^2/(\sqrt{6}\sigma)$.
Note that most of the initial conditions $\theta_1$, $\theta_2$ in Fig.~\ref{fig:phase}(a) show the divergence state (MDS or NDS) for which $d$ is not constant in time.
Hence, the particles do not show closed orbits characterized by $T_0$.

When $\mu$ is finite, the eigenvalue $\alpha$ has a negative real part, ${\rm Re}(\alpha)<0$, leading to exponentially decaying behavior.
As seen in Fig.~\ref{fig:phase}(e), the fixed points at $(\pi/2,\pi)$ and $(\pi,\pi/2)$ in Fig.~\ref{fig:phase}(d) change to stable spirals that are located approximately at $(\pi/2,3\pi/4)$ and $(3\pi/4,\pi/2)$.
Note that the hyperbolic saddle point remains at ($\pi/2,\pi/2$).
Due to the emerging spirals, most initial values converge to one of the fixed points, leading to constant angles $\theta_1,\theta_2$ in the long-time limit.
Since the spiral centers are not inside the shaded region, particles show the divergence spiral state rather than the convergence spiral state.
Apart from the spirals, some of the initial values do not directly converge to the spiral centers located approximately at $(3\pi/4,\pi/2)$ and $(\pi/2,3\pi/4)$.
Thus, the particles show reciprocating trajectories at an early stage, while they converge to the spirals in the long-time limit with trajectories as in Fig.~\ref{fig:mu3}(b).

In the limit of infinitely large odd viscosity $(\mu\to\infty)$, ${\rm Re}(\alpha)\to0$, and the fixed points with closed orbits again emerge at $(\theta_1,\theta_2)\approx(0.37\pi,0.64\pi)$ and $(0.64\pi,0.37\pi)$, as can be seen in Fig.~\ref{fig:phase}(f).
Active particles show the ``monotonic divergence state'' (MDS) along the diagonal line $\theta_1=\theta_2$, and any perturbation from this line kicks the particles away into either a ``circular oscillation state'' (COS, green circle) or ``radial oscillation state'' (ROS, red star).
The centers are characterized by the period $T_\infty=4\sqrt{3}\pi^2|\eta_{\rm o}|d^2/(\sqrt{5}\sigma)\sim|\eta_{\rm o}|$ and the oscillations are governed by the magnitude of odd viscosity.
This oscillation period is now determined by the ratio between the odd-viscous resistance $|\eta_{\rm o}|d$ and the characteristic force due to hydrodynamic interactions $\sigma/d$.

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

When the odd viscosity is zero, the active particle pair mostly shows ``monotonic divergence state'' [MDS, Fig.~\ref{fig:mu0}(b)] and ``nonmonotonic divergence state'' [NDS, Fig.~\ref{fig:mu0}(c)].
Generally speaking, one can see that the azimuthal component $\widehat{\mathbf{d}}$ in Eq.~(\ref{eq:velocity}) due to the induced flow causes the particle pair to unbind and run away from each other.
This leads to the above diverging states.
When the odd viscosity is finite, the particle pair displays spiral states (Fig.~\ref{fig:mu3}) that are due to both the azimuthal $\widehat{\mathbf{d}}_1$ and nonreciprocal components $\widehat{\mathbf{d}}^\ast_1$.
However, when the odd viscosity dominates, the $\widehat{\mathbf{d}}$-term vanishes and the nonreciprocal term $\widehat{\mathbf{d}}^\ast$ causes the pair to exhibit an oscillatory behavior (Figs.~\ref{fig:mu_infinity} and \ref{fig:mu_infinity_star}).
These results demonstrate that the nonreciprocal interaction due to odd viscosity gives rise not only to chiral spiral behavior, but also to oscillatory pair dynamics at sufficiently large values of the odd viscosity.

To see the effect of odd viscosity on the dynamics of a pair of active particles, we compare our results with those for a passive incompressible fluid without odd viscosity~\cite{manikantan2020}.
At length scales smaller than the hydrodynamic screening length, Manikantan showed~\cite{manikantan2020} that pairs of particles display oscillatory dynamics in a 1D coordinate system where one particle follows the other one.
In contrast, our results revealed that the particles show an oscillatory behavior in 2D (as in Figs.~\ref{fig:mu_infinity} and \ref{fig:mu_infinity_star}) due to the nonreciprocal interaction term $\epsilon_{ij}$ in Eq.~(\ref{eq:G}).
On the other hand, at length scales larger than the hydrodynamic screening length, the odd-viscosity term becomes exponentially smaller than the regular viscosity terms, as seen in Eq.~(\ref{eq:Glarge}).
This leads to vanishing chiral and oscillatory dynamics.
Furthermore, clustering behavior due to momentum leakage can emerge, as was observed before~\cite{manikantan2020}.

The inter-particle distance $r$ shown in Fig.~\ref{fig:mu_infinity_rel} represents the oscillatory behavior whose period is determined from the eigenvalue, $T_\infty/(|\eta_{\rm o}|a^2/\sigma)\approx490$, where $d=4a$ is used.
Note that the ``radial oscillation state'' (ROS, shown in Fig.~\ref{fig:mu_infinity_star}) is not determined by $T_\infty$ because the oscillatory behavior is not due to the closed orbits around the centers located approximately at $(0.37\pi, 0.64\pi)$ or $(0.64\pi,0.37\pi)$.
This can be understood by the fact that the area with green circles in Fig.~\ref{fig:phase}(c) almost overlaps with the area around these centers in Fig.~\ref{fig:phase}(f).
Note that such an oscillatory behavior is not observed in passive systems, where the oscillations are governed solely by closed orbits around fixed points in dynamical systems~\cite{manikantan2020}.

It is important to give numerical estimates of the physical quantities in our model.
Using typical values for enzymatic molecules, such as a molecular size $a\approx\SI[parse-numbers=false]{10^{-9}}{\metre}$ and typical exerted force $f\approx\SI[parse-numbers=false]{10^{-11}}{\newton}$, we obtain $\sigma\approx\SI[parse-numbers=false]{10^{-20}}{\newton\metre}$.
When the extent of the active chirality in the 2D fluid is high, we can set $\eta_{\rm o}\approx\SI[parse-numbers=false]{10^{-8}}{\pascal\second\metre}$ and $\eta_{\rm s}\approx\SI[parse-numbers=false]{10^{-9}}{\pascal\second\metre}$ for the 2D viscosities, leading to the limiting condition $\mu\gg1$.
Then, the period for the circular oscillations can be estimated as $T_\infty\approx3\times\SI[parse-numbers=false]{10^{-5}}{\second}$ when $d=a$.
Since this time scale is comparable to the period of biomolecular chemical reactions catalyzed by fast enzymes such as catalase or urease~\cite{hosaka2020_2}, odd viscosity would affect the chemical reactions through the enzymatic collective behavior.
Experimentally, such circular pair dynamics can be investigated by observing the response of the surrounding passive inclusions~\cite{svetlizky2021}.
Using optical tweezers, Svetlizky \textit{et al.}\ examined the response of a 2D colloidal suspension to a localized circular motion.
Their result suggests that the circular motion of a pair particle in the presence of odd viscosity can be evaluated by the symmetry breaking in the correlation function of the surrounding passive particles~\cite{svetlizky2021}.

\section{Conclusion and outlook}
\label{sec:conc}

In this paper, we have presented a theory of the dynamics of a pair of active force dipoles in a 2D active chiral fluid characterized by an odd viscosity $\eta_{\rm o}$.
The 2D active chiral fluid is described by a mobility tensor in Eq.~(\ref{eq:G}) with an asymmetric (nonreciprocal) part that is a direct consequence of a finite odd viscosity $(\eta_{\rm o}\neq0)$.
Without the odd viscosity, the particle pair shows both convergence and divergence states, as shown in Fig.~\ref{fig:mu0}.
However, with finite values of the odd viscosity, the particles start to exhibit various chiral pair dynamics, such as spiral (Fig.~\ref{fig:mu3}), circular (Fig.~\ref{fig:mu_infinity}), and radial trajectories (Fig.~\ref{fig:mu_infinity_star}), and this chiral dynamics is due to the antisymmetric (nonreciprocal) component of the mobility tensor.
The circular trajectory can be understood as arising from closed orbits at sufficiently large odd viscosity.
Furthermore, our results demonstrate that the nonreciprocal interaction due to the odd viscosity leads to a wealth of pair dynamics of active particles.

In this work, we have considered a transverse friction coefficient $\zeta_\perp$ between a 2D fluid in contact with an underlying 3D fluid.
In other studies, such transverse momentum leakage from the 2D layer to the 3D fluid beneath it was attributed to an anisotropic substrate~\cite{soni2019} or a Coriolis force in geophysical systems~\cite{tauber2019,tauber2020}.
Here, we show that transverse friction can originate from the odd viscosity of the 3D fluid.
When such an odd viscosity is present in the underlying 3D fluid, an additional term $\eta_{\rm o}^{\rm 3D}\partial_z \nabla^{\rm 3D}\times\mathbf{u}$ will appear in the 3D Stokes equation~\cite{markovich2021,khain2022}, where $\nabla^{\rm 3D}$ is the 3D gradient operator and $\mathbf{u}$ is the 3D fluid velocity.
By employing the lubrication approximation for the 3D fluid, we can obtain the force arising from the 3D odd viscosity on the 2D fluid as $\eta_{\rm o}^{\rm 3D}\mathbf{v}^\ast/h$.
This means that the 3D odd viscosity gives rise to a transverse friction coefficient, $\zeta_\perp=\eta_{\rm o}^{\rm 3D}/h$.
A more detailed discussion of such transverse flows due to odd viscosity will be given elsewhere.

The 2D odd viscosity can vary in the 2D plane, although such an effect was not considered in the present work.
In actual biomembranes, active proteins are often inhomogeneously distributed in the membrane and form active protein-rich domains that are called \textit{lipid rafts}.
The odd viscosity can then be different inside and outside the raft domain~\cite{hosaka2021_2}, and the spatial variation of the odd viscosity needs to be considered.
Some micro-organisms that are known to exhibit directional motions along viscosity gradients are called viscotaxis~\cite{liebchen2018}.
Hence, it will be of value to explore how the nonreciprocal flow field due to odd viscosity gradients couples with the viscotaxis convection.
This interesting question is left for future investigations.

\begin{acknowledgments}

We thank A.\ Kakugo and K.\ Yasuda for fruitful discussions and helpful suggestions.
D.A.\ acknowledges support from the Israel Science Foundation (ISF) under Grant No.\ 213/19 and support by the National Natural Science Foundation of China (NSFC) - ISF joint program under Grant No.\ 3396/19.
S.K.\ acknowledges the support by the National Natural Science Foundation of China  (NSFC)  (Nos.\ 12274098 and
12250710127) and the startup grant of Wenzhou Institute, University of Chinese Academy of Sciences
(No.\ WIUCASQD2021041).
\end{acknowledgments}

\section*{Data availability}

The data that support the findings of this study are available from the corresponding author upon reasonable request.

\begin{widetext}