\section{Introduction}

\textcolor{black}{In recent years, periodically driven many-body systems have attracted vast experimental and theoretical efforts.  Ergodicity-breaking phenomena in such systems, especially the concept of discrete time crystal (DTC) \cite{DTC1,DTC2,DTC3,DTC4,DTC5,DTC6,DTC7,DTC8,DTC9,DTC10,DTC11,DTC12,DTC13,DTC14,DTC15,DTC16,DTC17,DTC18,DTC19,DTC20,DTC21,DTC22,DTC23,DTC24,DTC25,DTC26,DTC27,DTC29,DTC30,DTC31,DTC32} have opened up new research avenues for
quantum simulation \cite{DTCqs,DTCqs2}, quantum computation \cite{DTCqc},  and studies of other condensed matter phenomena in the time-domain \cite{DTCcm1,DTCcm2,DTCcm3,DTCcm5,DTCcm6}.  Indeed, owing to the complexity of driven many-body systems,}  \textcolor{black}{theoretical and experimental advances} \textcolor{black}{may potentially lead to vast opportunities towards quantum supremacy, where experimental demonstrations or findings can be beyond the capability of any classical computer.}

\textcolor{black}{The nonequilibirum nature of periodically driven systems can also be exploited to access or generate a variety of topological phases of matters, many of which are not available in equilibrium systems. In particular, studies of non-equilibrium Floquet topological matter \cite{FMF1,FMF2,FMF3,DTCqc,FMF5,FMF6,FMF7,FMF8} have revealed that symmetry-protected topological excitation in periodically many-body systems with period $T$ can have $\pi/T$ quasienergy excitations on top of zero quasienergy excitations.} \textcolor{black}{The topologically protected $\pi/T$ quasienergy excitation are of particular interest because they result in all Floquet eigenstates coming in pairs with $\pi/T$ quasienergy separation. As an immediate consequence, a linear superposition of two $\pi/T$ quasienergy separated eigenstates can be formed, which directly breaks the discrete time translational symmetry and leads to subharmonic oscillations in the system dynamics. This simple physical picture emphasizes the importance of investigating the topological and the spectral (hence also the dynamical) properties together to explore interesting many-body physics. }

\textcolor{black}{Ongoing studies of periodically driven many-body systems are mostly based on one-dimensional systems. For instance, the suppression of thermalization via many-body localization (MBL) \cite{MBL1,MBL2,MBL3,MBL4} is known to be important in one-dimensional (1D) DTCs, but how MBL works in higher dimensions is largely unknown.  This perhaps explains why best known examples of DTCs is established in 1D setting.  Developing more theoretical insights  into periodically driven many-body systems in higher-dimensions is thus important from theoretical point of view and also necessary to guide future experimental designs.}

\textcolor{black}{As a first step towards understanding higher-dimensional Floquet many-body systems, this work aims to reveal some important physics of periodically driven spin ladders at moderate system sizes. Special emphasis is placed on the subharmonic magnetization response to the periodic driving.} \textcolor{black}{Because we study clean systems here, such subharmonic response cannot be taken as a genuine signature of DTCs in the strict sense.} \textcolor{black}{Nevertheless, at finite system sizes, this dynamical feature simulates that of a true DTC over a considerable time scale. Studying subharmonic magnetization dynamics} \textcolor{black}{also provides an important angle to digest how spin-spin coupling in the rung direction of the ladder impacts on the spectral properties of the system.  Referring to a terminology from studies of non-equilibrium Floquet topological matter \cite{FMF1,FMF2,FMF3,DTCqc,FMF5,FMF6,FMF7,FMF8}, the periodically driven spin ladders can accommodate the so-called Majorana $\pi$ modes (MPMs), i.e.,  mutually anticommuting Hermitian operators that also anticommute with the system's one-period (Floquet) propagator. Specifically, by rewriting the system Hamiltonian specified below in its equivalent Majorana lattice model,  these MPMs correspond to, \textcolor{black}{at least for special parameter values}, the topologically protected Majorana modes at the two lattice corners, in the spirit of second-order topological phases \cite{HTI-1,HTI0,HTI1,HTI2,HTI3,HTI4,HTI5,HTI6,HTI7,HTI8,HTI9,HTI10,HTI11,HTI12,HTI12b,HTI12c,HTI13,HTI14,HTI15,HTI16,HTI17,HTI18,HTI19,HTI20,HTI21,HTI22,HTI23,HTI24,HTI25,HTI26,FHTI1,FHTI2,FHTI3,FHTI4,FHTI5,FHTI6,FHTI7,FHTI8,FHTI9}. With this understanding, the spectral properties of a spin ladder at finite system sizes can be investigated in terms of the} \textcolor{black}{fate of these corner MPMs at various system sizes and parameters.}

	This paper is organized as follows. In Sec.~II, we present our model and elucidate the role of MPMs in yielding subharmonic dynamics. In Sec.~III, we compare the proposed model with its 1D counterpart in terms of their Floquet eigenstates structure and computationally investigate its implication on the robustness of the $2T$ magnetization dynamics against variation in the system size and choice of initial states. Finally, we conclude this paper, briefly discuss the potential experimental realization of the proposed model, and present avenues for future work in Sec.~IV.

	\section{Periodically driven spin ladder}
	Consider a quasi-one-dimensional rectangular lattice of size $N_x\times N_y$ subjected to the time-periodic Hamiltonian
	\begin{widetext}
		\begin{equation}
		H(t) = \begin{cases}
		-\sum_{\langle i,j \rangle} \left(J_x \sigma_{z,i,j} \sigma_{z,i+1,j} + J_y\sigma_{z,i,j+1} \sigma_{z,i,j}\right) & \text{ for } nT<t\leq (n+1/2)T \\
		\sum_{i,j} h\sigma_{x,i,j} & \text{ for } (n+1/2)T<t\leq (n+1)T
		\end{cases} \;, \label{sys}
		\end{equation}
	\end{widetext}
	where $\left\lbrace\sigma_{s,i,j}:s=x,y,z\right\rbrace$ is the set of three Pauli matrices acting at lattice site $(i,j)$, $J_x$ and $J_y$ are the ZZ interactions in the $x$- and $y$-directions respectively, $h$ is the Zeeman field strength, $T$ is the driving period, and $n\in \mathbb{Z}$. While Eq.~(\ref{sys}) can generally be interpreted as a two-dimensional (2D) system, we choose to term it a spin ladder for the following reasons: (i) we work in the regime $J_x\ll J_y$ and (ii) we fix $N_y=2$ in the most case studied in this paper.

	As a time-periodic Hamiltonian, Eq.~(\ref{sys}) can be characterized in terms of its one-period time evolution operator (to be called \emph{Floquet operator} onwards) \cite{Flo1,Flo2}
	\begin{equation}
	U \equiv U(T;0)=\mathcal{T} \exp\left(-\mathrm{i} \int_0^T H(t) dt \right)\;,
	\end{equation}
	where $\mathcal{T}$ is the time-ordering operator and $\hbar=1$ is taken here and in the remainder of this paper. The system's quasienergy $\varepsilon \in (-\pi/T,\pi/T]$ and its associated Floquet eigenstate $|\varepsilon\rangle$ are defined from the eigenvalues $(\exp\left[-\mathrm{i}\varepsilon T\right])$ and eigenstates of $U$ respectively. A quasienergy is only defined modulo $2\pi/T$, so that $\varepsilon$ and $\varepsilon+2\pi\ell/T$ with $\ell\in\mathbb{Z}$ represent the same physics.

The Floquet operator associated with Eq.~(\ref{sys}) can be easily obtained as
 \begin{eqnarray}
U &=& e^{-\mathrm{i} \sum_{i,j} \frac{hT}{2} \sigma_{x,i,j}}  \nonumber \\
&& \times e^{\left(\sum_{i,j} \mathrm{i} \left[\frac{J_yT}{2}\sigma_{z,i,j+1}\sigma_{z,i,j}+\frac{J_xT}{2} \sigma_{z,i+1,j}\sigma_{z,i,j}\right]\right)} \;, \label{Flospin}
\end{eqnarray}
By employing the generalized JW transformation detailed in Appendix~\ref{app1}, the different spin operators appearing in Eq.~(\ref{Flospin}) can be written in terms of mutually anticommuting Majorana operators $\gamma_{A,i,j}$ and $\gamma_{B,i,j}$ as
\begin{eqnarray}
\sigma_{x,i,j} &=& \mathrm{i} \gamma_{A,i,j} \gamma_{B,i,j} \;, \sigma_{z,i,j+1}\sigma_{z,i,j} = \mathrm{i} \gamma_{B,i,j} \gamma_{A,i,j+1} \;, \nonumber \\
\sigma_{z,i+1,j}\sigma_{z,i,j} &=& \mathrm{i} \gamma_{B,i,j} \gamma_{A,i+1,j} \prod_{j<n\leq N_y} \left(\gamma_{A,i,n} \gamma_{B,i,n}\right) \nonumber \\
&& \times \prod_{m<j} \left(\gamma_{A,i+1,m} \gamma_{B,i+1,m}\right) \;. \label{opmap}
\end{eqnarray}
Figure 1 illustrates some representative spin operators above when $\gamma_{A,i,j}$ and $\gamma_{B,i,j}$ are arranged into a
2D square lattice. The highly nonlocal terms arising from $\sigma_{z,i+1,j}\sigma_{z,i,j}$ are unavoidable \cite{Chapman} and highlight the strongly interacting nature of our system.
\begin{center}
	\begin{figure}
		\includegraphics[scale=0.4]{pics2.pdf}
		\caption{Schematic of various terms in Eqs.~(\ref{sys}) and (\ref{Flospin}) when written in terms of Majorana operators. There, each term is a product of all Majorana operators contained in its respective dotted line.}
		\label{pic1}
	\end{figure}
\end{center}

\subsection{MPMs as topological quasienergy excitations}

Note first that at $N_x=1$, Eq.~(\ref{sys}) reduces to a clean (disorder-free) variation \cite{DTC8} of the  spin-based 1D DTC models \cite{DTC2,DTC5} that is analytically solvable,
with the presence or absence of Majorana zero modes (MZMs) and MPMs yielding the phase diagram in Fig.~\ref{result1}(a) (See Appendix~\ref{app1} for its analytical derivation).  Each regime in the phase diagram is labeled by the four phases identified in Ref.~\cite{DTC5}. Specifically, regardless of the boundary conditions, the paramagnet (PM), $0$-spin glass ($0$-SG), and $\pi$-spin glass ($\pi$-SG) host, respectively, neither MPM nor MZM, a pair of MZMs, and a pair of MPMs. On the other hand, in the $0\pi$-paramagnet regime, the system hosts a pair of MZMs and a pair of MPMs only under open boundary conditions (OBC); otherwise, it does not support any MZM or MPM. Outside $(h,J_y)\in [0,\pi/2]\times [0,\pi/2]$ (in $T/2=1$ units), the four phases above repeat themselves, e.g., $\pi$-SG phase repeats up to $h=\pi$ at $J_y=\pi/4$, followed by another $0$-SG phase as $h$ further increases.
\begin{center}
	\begin{figure*}
		\includegraphics[scale=1.0]{fig1.pdf}
		\caption{(a) Phase diagram of the periodically driven Ising spin chain Eq.~(\ref{sys}) in the 1D limit, i.e., $N_x=1$. (b) Stroboscopic magnetization dynamics at $N_x=1$, $J_y=0.6\pi/T$, and $h=0.8\pi/T$. (c) The subharmonic response of the system's average magnetization at $\Omega=\pi/T$, $J_x=0.05\pi/T$, $J_y=0.6\pi/T$, and various lattice configurations. There, the system is initialized in a product state $|\uparrow\cdots \uparrow \rangle$ and the average magnetization is evolved up to $2000T$. (d) Perturbation-induced degeneracy splitting may break the $\pi/T$ quasienergy spacing structure.}
		\label{result1}
	\end{figure*}
\end{center}

\textcolor{black}{Since the case of $N_x=1$ can be exactly solved, it is straightforward to digest} the role of MPMs in partially establishing two main dynamical features analogous to those of DTCs \cite{DTC2,DTC3,DTC8}: i) the existence of $2T$-periodic state and ii) its robustness against variation in the system parameters. To this end, we first define a \emph{quasienergy excitation} operator $\gamma_\xi$ satisfying $U\gamma_\xi U^\dagger=\exp\left(-\mathrm{i} \xi\right)\gamma_\xi $. As its name suggests, it ``excites" a given quasienergy eigenstate $|\varepsilon\rangle$ to another $|\varepsilon+\xi\rangle=\gamma_\xi |\varepsilon \rangle $ with quasienergy $\varepsilon+\xi$. In particular, MPMs are defined as quasienergy $\pi/T$ excitations, whose existence implies the presence of two quasienergies $\varepsilon$ and $\varepsilon+\pi/T$ in the system. In this case, a superposition state $|\varepsilon\rangle +|\varepsilon+\pi/T\rangle$ returns to itself upon applying $U$ twice (is $2T$ periodic), thus establishing the first feature above. It is also worthwhile noting that, as further discussed below, each $|\varepsilon\rangle$ and $|\varepsilon+\pi/T\rangle$ ideally takes the form of a symmetric/anti-symmetric cat state. \textcolor{black}{Consequently, an equal weight superposition of $|\varepsilon\rangle$ and $|\varepsilon+\pi/T\rangle$, which demonstrates $2T$ periodicity, yields a product state. Such a product state, and thereby its associated subharmonic dynamics, is naturally robust against decoherence mechanism that tends to break entanglements.}

Second, the presence of MPMs in the above model can be traced back from the existence of \emph{topological edge modes} arising in their equivalent Majorana fermionic description via a Jordan-Wigner (JW) transformation or its generalization (see Eq.~(A1) in Appendix~\ref{app1}). Consequently, such MPMs are not a result of fine tuning and exist for a range of system parameters, thereby preserving the $2T$ periodic state identified above. Moreover, due to the nonlocal nature of the JW-like transformation, the original spin system can always be mapped to its corresponding Majorana fermionic description under OBC, possibly with additional nonlocal terms which may not always lead to the pairwaise annihilation of MPMs. This in turn yields the $\pi$-SG regime scenario in which MPMs exist \emph{regardless} of the boundary conditions applied to the original spin system.

\subsection{Robustness of corner MPMs}

At $N_x>1$ and $hT/2=\pi/2$, a quick inspection of Fig.~\ref{pic1} allows the identification of $\gamma_{A,1,1}=\sigma_{z,1,1}$ and $\gamma_{B,N_x,N_y}=\left(\prod_{(l<N_x),(m\leq N_y)} \sigma_{x,l,m} \right) \left(\prod_{(l<N_y)} \sigma_{x,N_x,l} \right) \sigma_{y,N_x,N_y}$ as the system's two MPMs.  {Note that in the equivalent Majorana description of our system, $\gamma_{A,1,1}$ and $\gamma_{B,N_x,N_y}$ are precisely located at the two corners (see Fig.~\ref{pic1}), suggesting that they may share a similar origin as the topological corner modes in typical Floquet second-order topological superconductors \cite{FHTI8,FHTI9}}.

\textcolor{black}{In general conditions, cases with $N_x\ne 1$} \textcolor{black}{are no longer integrable.} \textcolor{black}{This is precisely why the present study for limited system sizes is still necessary and useful.}  Though the fate of the above identified corner MPMs is no longer analytically tractable at $h\neq \pi/2$, we may still numerically probe their existence by defining and evaluating the two corner-state spectral functions at zero and $\pi/T$ quasienergy, i.e.,
\begin{eqnarray}
s_{0,1} &=& \mathcal{N}_1  \sum_{n\in \mathcal{X}} \int_{-\epsilon}^{\epsilon} S_1(\varepsilon_n, \eta) d\eta  \;, \nonumber \\
s_{0,2} &=& \mathcal{N}_2 \sum_{n\in \mathcal{X}} \int_{-\epsilon}^{\epsilon} S_2(\varepsilon_n, \eta) d\eta \;, \nonumber \\
s_{\pi,1} &=& \mathcal{N}_1 \sum_{n\in \mathcal{X}} \left( \int_{-\pi/T-\epsilon}^{-\pi/T+ \epsilon} S_1(\varepsilon_n, \eta) d\eta + \int_{\pi/T-\epsilon}^{\pi/T+ \epsilon} S_1(\varepsilon_n, \eta)  d\eta \right) \nonumber \;, \\
s_{\pi,2} &=& \mathcal{N}_2 \sum_{n\in \mathcal{X}} \left( \int_{-\pi/T-\epsilon}^{-\pi/T+ \epsilon} S_2(\varepsilon_n, \eta) d\eta + \int_{\pi/T-\epsilon}^{\pi/T+ \epsilon} S_2(\varepsilon_n, \eta)  d\eta \right) \nonumber \;,
\end{eqnarray}
where $\epsilon\ll 1$, $\mathcal{X}$ represents a rather arbitrarily chosen integer between $1$ and the size of the system's Hilbert space, $\mathcal{N}_1$ and $\mathcal{N}_2$ are normalization constants such that $\sum_{n\in \mathcal{X}} \int_{-\pi/T}^{\pi/T} S_1(\varepsilon_n,\eta) d\eta = \sum_{n\in \mathcal{X}} \int_{-\pi/T}^{\pi/T} S_2(\varepsilon_n,\eta) d\eta =1$, and
\begin{eqnarray}
S_1(\varepsilon_n, \eta) &=& \sum_{\varepsilon_m} \delta(\varepsilon_n-\varepsilon_m -\eta) |\langle \varepsilon_n | \gamma_{A,1,1} | \varepsilon_m \rangle |^2 \;, \nonumber \\
S_2(\varepsilon_n, \eta) &=& \sum_{\varepsilon_m} \delta(\varepsilon_n-\varepsilon_m -\eta) |\langle \varepsilon_n | \gamma_{B,N_x,N_y} | \varepsilon_m \rangle |^2 \;.
\end{eqnarray}
In particular, the quantities $s_{0,1}$, $s_{0,2}$, $s_{\pi,1}$, and $s_{\pi,2}$ are adapted from the edge spectral functions defined in Ref.~\cite{pzmandppm} and can be understood as follows.

Since $\gamma_{A,1,1}$ and $\gamma_{B,N_x,N_y}$ represent exact $\pi/T$ quasienergy excitations in the ideal $h=\pi/2$ case, any $|\varepsilon_m\rangle $ is mapped onto some $|\varepsilon_n\rangle $ with $|\varepsilon_n-\varepsilon_m|=\pi/T$. It then follows that $s_{0,1}=s_{0,2}=0$, \textcolor{black}{but $s_{\pi,1}=s_{\pi,2}=1$}. Away from the ideal case, $\gamma_{A,1,1}$ and $\gamma_{B,N_x,N_y}$ no longer represent the system's true corner MPMs. However, if corner MPMs still exist, they are expected to have large support on $\gamma_{A,1,1}$ and $\gamma_{B,N_x,N_y}$, thus yielding $s_{0,1},s_{0,2}\approx0$ and $s_{\pi,1},s_{\pi,2}\approx 1$. \textcolor{black}{To this end, since these corner MPMs, if exist, are not exactly $\pi/T$ quasienergy excitations at finite system sizes, a state $|\varepsilon_n\rangle$ may be mapped onto $|\varepsilon_m\rangle$ with  $|\varepsilon_n-\varepsilon_m|=\Delta \varepsilon_{n,m}\neq \pi/T$. This motivates the integration of $S_1(\varepsilon_n,\eta)$ and $S_2(\varepsilon_n,\eta)$ with respect to $\eta$ over a small quasienergy window centered around $\pm \pi/T$ in the definition of $s_{\pi,1}$ and $s_{\pi,2}$ above. Moreover, since $\Delta \varepsilon_{n,m}$ may not be the same over all Floquet eigenstates, the sampling of $S_1(\varepsilon_n,\eta)$ and $S_2(\varepsilon_n,\eta)$ over several Floquet eigenstates is necessary. This justifies the summation over $\chi$ appearing in $s_{\pi,1}$ and $s_{\pi,2}$ above. Similar consideration is taken into account in defining $s_{0,1}$ and $s_{0,2}$ to properly probe the presence of zero quasienergy excitations.} We numerically evaluate these corner spectral functions in Fig.~\ref{spectral} for a system of size $4\times 2$ under periodic boundary conditions (PBC) in both directions. In particular, not only do we verify the robustness of the corner MPMs over a considerable range of parameter values, we also find evidence that the range of parameter values over which these corner MPMs exist matches the phase diagram of the system's 1D counterpart, i.e., Fig.~\ref{result1}(a). That is, by fixing $J_y=0.6\pi/T$ in panel (a), we observe that $s_{\pi,1}$ and $s_{\pi,2}$ will only become close to zero for $h<0.6\pi/T\ (T=2)$. Again consistent with the phase diagram in Fig.~\ref{result1}(a), for $h<0.4\pi/T\ (T=2)$, the values of $s_{0,1}$ and $s_{0,2}$ shoot up, thus hinting the existence of corner MZMs.  Likewise, by fixing $h=0.8\pi/T$ in panel (b), we observe that $s_{\pi,1}$ and $s_{\pi,2}$ remain finite for $0.2\pi/T <J_y < 0.8\pi/T\ (T=2)$.  \textcolor{black}{These detailed analysis hence clearly indicates that the MPMs (as seen below, MZMs are not of much interest in the present study)  in the original one-dimensional spin chain survives in the spin ladder system, but now in the form of second-order topological corner modes over a considerable range of system parameters.}

\begin{center}
	\begin{figure}
		\includegraphics[scale=0.51]{cornerspectral.pdf}
		\caption{The four corner spectral functions (with $\chi=16$ and $\epsilon=0.01$) plotted as a function of (we set $T=2$) (a,c) $h$ and (b,d) $J\equiv J_y$. The system size is taken as $4\times 2$ and the other systems parameters are $J_x=0.05\pi/2$, (a,c) $J_y=0.3\pi$ and (b,d) $h=0.4\pi$.}
		\label{spectral}
	\end{figure}
\end{center}

It is then interesting to connect the above-identified corner MPMs with the subharmonic response of the system, specifically, the $2T$ magnetization dynamics. This angle of analysis is certainly motivated by how DTCs are examined.  Further, in a recent experiment on an acoustic platform \cite{Zhangexp}, it is observed that a topological corner $\pi$ mode can be used to generate stable period-doubling oscillations after a corner excitation.  Certainly, the subharmonic response in our spin ladder system is of many-body nature; it involves the bulk many-body states in general because the MPMs are highly nonlocal in the spin representation.
To present our computational findings regarding  the 2T magnetization dynamics,  we first define the system's magnetization as
\begin{equation}
\langle S_z (t)\rangle = \sum_{i=1}^{N_x} \sum_{j=1}^{N_y} \langle \sigma_{z,i,j} \rangle \;. \label{magnetization}
\end{equation}

\subsection{Effect of the system size on the sub-harmonic response in 1D}

\textcolor{black}{The subharmonic response associated with genuine DTCs is known to exhibit an infinite life-time in the thermodynamic limit.  With this insight it would be interesting to discuss how the $2T$ magnetization dynamics in our clean spin ladders changes upon increasing the system size.
To that end we start with the 1D limit} of Eq.~(\ref{sys}), i.e., by taking $N_x=1$, and analyze the effect of increasing the system size $N_y$. First, note that MPMs generally have a finite correlation length. This results in some overlap between two MPMs localized at the opposite ends, which prevents them from being true $\pi/T$ quasienergy excitations. This explains the imperfect $2T$ magnetization dynamics with a finite lifetime observed at $h\neq \pi/2$ and small system size (see e.g. Fig.~\ref{result1}(b)). Increasing the system size reduces the overlap between two MPMs, which in turn brings them closer to being $\pi/T$ quasienergy excitations and leads to improvement in the $2T$ magnetization dynamics and its lifetime up to a critical system size $N_{y,\rm crit}$ (see Fig.~\ref{result1}(c)).

Remarkably, increasing the system size beyond $N_{y,\rm crit}$ is found to weaken the $2T$ magnetization dynamics.  This finding may be partially understood as follows. At $hT/2=\pi/2$, the initial state $|\uparrow\cdots \uparrow \rangle $ can be exactly written as a superposition of two Floquet eigenstates with quasienergies $-N_y J_y$ and $-N_y J_y+\pi/T$. In particular, at sufficiently small system sizes, these two quasienergies are generally nondegenerate. However, due to the periodicity of the quasienergy Brillouin zone $\varepsilon \in \left( -\pi/T, \pi/T\right]$, the quasienergies $-N_y J_y$ and $-N_y J_y+\pi/T$ are only defined modulo $2\pi/T$. As the system size ($N_y$) becomes sufficiently large, these values may loop back in the quasienergy Brillouin zone. This, coupled with the fact that the number of quasienergy eigenstates scales up significantly with the system size,  strengthens the chance of degeneracy or near-degeneracy between many pairs of Floquet eigenstates. If such type of unavoidable degeneracy involves states differing by a single spin excitation, it becomes extremely sensitive to perturbation. Indeed, a small deviation from $hT/2=\pi/2$ may result in considerable hybridization between two such Floquet eigenstates, thus splitting their quasienergy by an arbitrary value. This in turn results in the destruction of $\pi/T$ quasienergy pairing and, consequently, $2T$ magnetization dynamics (see Fig.~\ref{result1}(d)).

\textcolor{black}{One may also appreciate the loss of \textcolor{black}{period doubling dynamics} in terms of  the movement of the domain wall along the spin chain.  At $hT/2=\pi/2$, as long as the state has non-vanishing magnetization, e.g. the product state in the $\sigma_z$ basis $|\psi_0\rangle=|\downarrow\uparrow\cdots\uparrow\uparrow\rangle$, then the period doubling dynamics is self-evident. Indeed, in this case the Floquet operator $U$ can rewritten as $U=P_xU_p,$ where $P_x$ represents a perfect spin flip operator $$P_x=e^{-i\pi/2\sum_j\sigma_{x,j}}$$ and the dynamics emanating from one eigenstate of $U_p$ with nonvanishing magnetization apparently displays $2T$ oscillation dynamics.   For the same state but with $h=(\pi/T+\delta h)T/2\neq \pi/2$, we may first consider the zero-th order approximation to the effective Floquet Hamiltonian of $U_p=e^{-iH_{p,F}T}$, namely,

 $$H_{p,F}^{(0)}= -J_y/2\sum_j \sigma_{z,1,j} \sigma_{z,1,j+1} - \delta h/2\sum_j\sigma_{x,j},$$ which is the transverse field Ising Hamiltonian (TFIH). For  $\delta h<J_y$, the TFIH is in the ferromagnetic phase, where we have an almost degenerate ground state separated by a finite energy gap from all other excited states.  With the presence of the perfect spin flip operator $P_x$ in $U$, the originally degenerate states in $H_{p,F}^{(0)}$ take us to other Floquet states of $U$ with a quasienergy separation of $\pi/T$ in $U$.   The $\delta h$ term however induces the movement of the domain wall, initially introduced in the state $|\psi_0\rangle$,  along the spin chain.  Indeed, using again Jordan-Wigner transformation and Bogoliubov transformation, the $\delta h$ term behaves likes a hopping term for a domain wall.
 Thus, tuning away $h$ from the peculiar value $\pi/T$ generates a mobile domain wall, without a direct energy cost \textcolor{black}{to the first order approximation}.  This indicates that $|\psi_0\rangle$ is no longer an eigenstate of $U_p$, and
 $|\psi_0\rangle$ under the action of $U_p$ alone can hence diffuse into a linear combination of different states with the domain wall at different locations. Therefore, the joint action of $P_x$ and $U_p$ cannot bring the state back after two periods, hence equivalently there is loss of $\pi/T$ pairing. \textcolor{black}{In addition, these delocalized domain walls fluctuate and thus destroy the long-range order in the state, leading to a vanishing value of magnetization. Therefore, only the state with spins all up or all down in the 1D model can display the period doubling dynamics for a suffiecintly long time. However, with hybridization between different Floquet states, i.e. Floquet resonance, as discussed previously, domain walls inevitably emerge in the state under Floquet driving.}}

The above-analyzed loss of \textcolor{black}{persistent period doubling dynamics} can generally be avoided in the presence of strong disorder \cite{MBL1,MBL2,MBL3,MBL4} which \textcolor{black}{leads to localization of the domain walls and prevents Floquet heating}, or nonlocal interaction \cite{DTC9}. Together with the robust $\pi/T$ quasienergy pairing due to Majorana $\pi$ modes and their topological protection, one may then obtain genuine DTC phases with disorder.  Though the absence of any of these mechanisms precludes the clean systems under consideration here from being a true DTC, a better understanding of the subharmonic responses in a clean setting is still of great interest, especially to the quantum simulation of clean periodically driven systems.

\section{$2T$-magnetization dynamics of a driven spin ladder}
 \textcolor{black}{In this section, we shall show that a periodically driven spin ladder allows for an inherent mechanism for avoiding the above-mentioned dangerous degeneracy up to moderate system sizes. This new insight is interesting, because it suggests that the dynamics of driven spin ladders resembles that of DTCs at larger system sizes than in  the corresponding 1D spin chains. }

\subsection{Comparison with the 1D model}
To establish a simple and fair comparison between the spin ladder dynamics and the associated 1D model, we consider the model of Eq.~(\ref{sys}) at $N_x,N_y=2,2$ and $N_x,N_y=1,4$ respectively.

We start by considering the system configuration of $N_x,N_y=1,4$. At $hT/2=\pi/2$, the system's $16$ quasienergies can be explicitly obtained as $\varepsilon_{s_1,s_2,s_3,+}= -J_y \sum_{j=1} s_j$ and $\varepsilon_{s_1,s_2,s_3,-}= -J_y \sum_{j=1} s_j +\pi/T$, where $s_1,s_2,s_3=\pm 1$. At sufficiently small $J_y$, the quasienergies $\varepsilon_{1,1,1,\pm}$, which are associated with eigenstates $|\varepsilon_{1,1,1,\pm}\rangle \propto |\uparrow \uparrow \uparrow \uparrow \rangle \pm |\downarrow \downarrow \downarrow \downarrow \rangle$, may be nondegenerate. On the other hand, the other Floquet eigenstates are at least twofold degenerate, e.g., $\varepsilon_{+,-,-,+} = \varepsilon_{-,+,-,+}$ (see Fig.~\ref{qspec}(a, b) for the complete quasienergy spectrum). In particular, these degenerate Floquet eigenstates involve states differing by a single spin excitation, e.g., $|\uparrow \downarrow \uparrow \uparrow \rangle$ and $|\uparrow \downarrow \downarrow \uparrow \rangle$, which can be easily split by perturbation (see Fig.~\ref{qspec}(b)). As an immediate consequence numerically confirmed below, while choosing $|\uparrow\cdots \uparrow \rangle$ or $|\downarrow\cdots \downarrow \rangle$ as an initial state yields a robust $2T$ magnetization dynamics with respect to variation in the parameter $h$ (see Fig.~\ref{result1}(b)), other choices of initial states, which have significant support on the above degenerate Floquet eigenstates, no longer display the same robustness.

\begin{center}
	\begin{figure}
		\includegraphics[scale=0.5]{qspec.pdf}
		\caption{Quasienergy levels (and their $\pi/T$-shifted values) associated with system configurations (a,b) $N_x, N_y=1, 4$ and (c,d) $N_x,N_y=2,2$. We take (a,c) a solvable parameter value of $h=\pi/T$ associated with perfect $2T$ magnetization dynamics and (b,d) a perturbed value of $h=0.85 \pi/T$.}
		\label{qspec}
	\end{figure}
\end{center}

We now turn our attention to the system configuration of $N_x,N_y=2,2$. At $hT/2=\pi/2$, such a system is also analytically solvable with quasienergies $\varepsilon_{s_1,s_2,s_3,+}= -J_y (s_1+s_2)- J_y'(1+s_1s_2)s_3$. Similarly to the $N_x, N_y=1,4$ setting, Floquet eigenstates degeneracy exist except for $|\varepsilon_{-1,-1,-1,\pm}\rangle = |\uparrow \downarrow \downarrow \uparrow \rangle \pm |\downarrow \uparrow \uparrow \downarrow \rangle$ and $|\varepsilon_{1,1,1,\pm}\rangle = |\uparrow \uparrow \uparrow \uparrow \rangle \pm |\downarrow \downarrow \downarrow \downarrow \rangle$. However, in the $N_x,N_y=2,2$ configuration, the quasienergy degeneracy only involves states differing by at least two spins excitations. Since the source of perturbation in the present model comes from variation in the $h$ parameter characterizing a single spin rotation in the $x$-direction, such a degeneracy cannot be easily lifted, e.g., compare Fig.~\ref{qspec}(c) with Fig.~\ref{qspec}(d).

In Table~\ref{T1}, we compare the minimum and maximum quasienergy spacing (relative to $\pi/T$) at various system sizes. Consistent with the above analysis, the spin ladder configurations outperform their 1D counterpart in terms of having quasienergy spacing closer to $\pi/T$ at moderate imperfection $h\neq \pi/T$. Increasing the system size generally brings the quasienergy spacing closer to $\pi/T$. Exceptions to this pattern can be attributed to the presence of accidental degeneracy between states differing by a single spin excitation, which becomes more probable at larger system sizes due to the periodicity of the quasienergy Brillouin zone and the larger number of Floquet eigenstates.

\begin{center}
\begin{table}[]
    \centering
    \begin{tabular}{|c|c|c||c|c|c|}
    \hline
    Size & Min & Max & Max & Min & Size \\ \hline
        $2 \times 2$ & $1.20\times 10^{-3}$ & $0.023$ & $0.15$ & $1.81\times 10^{-3}$ & $1\times 4$\\
        $3\times 2$ & $5.42\times 10^{-6}$ & $0.016$ & $0.086$ & $2.08\times 10^{-4}$ & $1\times 6$\\
        $4\times 2$ & $2.06\times 10^{-5}$ & $0.017$ & $0.050$ & $2.49\times 10^{-5}$ & $1\times 8$ \\
        $5\times 2$ & $2.14\times 10^{-6}$ & $0.013$ & $0.032$ & $1.38\times 10^{-4}$ & $1\times 10$ \\
        $6\times 2$ & $4.53\times 10^{-7}$ & $0.015$ & $0.017$ & $1.37\times 10^{-5}$ & $1\times 12$ \\ \hline
    \end{tabular}
    \caption{Maximum and minimum deviation of the quasienergy spacing from $\pi/T$ at various system sizes. The system parameters are chosen as $h=0.85\pi/T$, $J_x=0.05\pi/T$, and $J_y=2/T$.}
    \label{T1}
\end{table}
\end{center}

\subsection{Computational results}
To support the above argument, we plot in Fig.~\ref{result5} the stroboscopic magnetization dynamics (see Eq.~(\ref{magnetization}))and its associated power spectrum, i.e., $\langle \tilde{S}_z\rangle = \frac{1}{M}\sum_{n=1}^M \langle S_z (nT) \rangle e^{-\mathrm{i} \Omega n T/M}$, at configurations $1\times 8$ and $4\times 2$ under the initial state of $|\uparrow \downarrow \uparrow \cdots \uparrow\rangle$. Indeed, while $2T$-periodicity remains robust in the $4\times 2$ setting, it is no longer visible in the purely 1D system.

\begin{center}
	\begin{figure}
		\includegraphics[scale=0.5]{result5.pdf}
		\caption{The stroboscopic magnetization dynamics (a, c) and its associated power spectrum (b, d) under the initial state of $|\uparrow \downarrow \uparrow \cdots \uparrow\rangle$. Panel (a, b) and (c, d) correspond to the configurations $N_x, N_y= 4\times 2$ and $1\times 8$ respectively. All parameter values are taken as $J_x=0.05\pi/T$, $h=0.85\pi/T$, and $J_y=2/T$.}
		\label{result5}
	\end{figure}
\end{center}

In Fig.~\ref{result3}(a, b), we further demonstrate the robustness of the $2T$ magnetization dynamics in the proposed $N_x\times N_y$ system by scanning its $\pi/T$ subharmonic peak over the parameter $h$ at two different initial states. First, we consider an initial state in which all but one spins is in the up ($+\hat{z}$) position, with the subharmonic response of the system's averaged magnetization shown in Fig.~\ref{result3}(a).  In this case, even a considerably large 1D lattice no longer demonstrates robust subharmonic response, as evidenced by the very narrow window of subharmonic peaks with respect to variation in a system parameter observed in Fig.~\ref{result3}(a). By contrast, in the spin ladder configurations, a significantly large subharmonic response is still observed over a much broader range of the parameter $h$.

Second, we consider an initial state in which all spins are aligned in a diagonal direction, e.g., $+\hat{n}_{\pi/8}=\cos(\pi/4)\hat{z}+\sin(\pi/4)\hat{y}$, then stroboscopically evaluate the magnetization in the same direction and present our results in Fig.~\ref{result3}(b). Here, we find that $h=\pi/T$ generically no longer represents the value at which the subharmonic response is maximum.  The relatively large subharmonic response of 1D configuration $1\times 16$ at $h=\pi/T$ is accidental, because the associated magnetization profile actually exhibits periodicity larger than $2T$ (See Appendix~\ref{app2}). Such a $>2T$ periodicity is merely a result of some fine tuning of parameter values and quickly disappears as $h$ is tuned away from $\pi/T$. As all spins are aligned at an angle of $\pi/4$ from the $+\hat{z}$-axis, the genuine maximum $2T$ subharmonic response is instead observed around $h=1.25\pi/T$, which can intuitively be understood as follows. Ignoring the interaction terms, a Floquet operator at $h=5\pi/(4T)$ effectively rotates the system state towards $|\uparrow\cdots \uparrow \rangle$, which is known to have considerable support on two nondegenerate Floquet eigenstates with quasienergy separation of $\pi/T$.  In the vicinity of $h=1.25\pi/T$, the magnetization profile vs time and the associated power spectrum are also shown in Fig.~\ref{result3}(c) and (d).

\begin{center}
	\begin{figure}
		\includegraphics[scale=0.5]{result3.pdf}
		\caption{(a,b) The subharmonic response of the system's average magnetization at $\Omega=\pi/T$ and various lattice configurations under two different choices of initial states (a) $|\uparrow \downarrow \uparrow \cdots \uparrow \rangle $ and (b) $|\nearrow \cdots \nearrow \rangle$. (c,d) The stroboscopic evolution of the average magnetization in the $\hat{n}_{\pi/8}$ direction and its associated power spectrum for the initial state $|\nearrow \cdots \nearrow \rangle$. There, green/black, magenta/cyan, and red/blue respectively correspond to system sizes $1\times 16$, $4\times 4$, and $8\times 2$. In all panels, the other systems are chosen as $J_x=0.05\pi/T$, $J_y=0.6\pi/T$, and (c,d) $h=1.225\pi/T$.}
		\label{result3}
	\end{figure}
\end{center}

\section{Concluding remarks}

\textcolor{black}{We have investigated the topological and dynamical features of  periodically driven spin ladder configurations up to moderate system sizes.  Aligned with how the dynamics of genuine DTCs is analyzed,  we examine the robustness and time scale of subharmonic response to the periodic driving,  in connection with the presence of MPMs at two corners in the equivalent Majorana lattice models. While such MPMs may become extremely nonlocal in the original spin representation, their role as quasienergy $\pi/T$ excitations ensures the system's full Floquet eigenstates to form $\pi/T$ quasienergy spacing, which then leads to $2T$ magnetization dynamics over a finite but considerable time scale.  The topological nature of such corner MPMs leads to the robustness of the obtained $2T$ periodic feature over a range of parameter values. The subharmonic response of spin ladders can survive at larger system sizes than its 1D counterpart, with less sensitivity to the choice of initial states.  We have also highlighted the role of spin-spin interaction in the rung direction in removing the degeneracy between states differing by a single spin excitation.  }

Due to the simplicity of our proposed model and its similarity with existing 1D DTC models, its experimental realization is straightforward. In particular, there exists various platforms on which periodically driven spin systems are already experimentally realized, which include trapped ions \cite{DTCexp1,DTCexp7}, Nitrogen-Vacancy centers in diamonds \cite{DTCexp2}, and nuclear magnetic moments \cite{DTCexp3,DTCexp4,DTCexp5}. Among these available platforms, trapped ions of Ref.~\cite{DTCexp1,DTCexp7} are the most suitable platform due to the flexibility in arranging the $^{171}Yb^+$ ions into various different lattice configurations, whereas their experimental execution may simply follow directly the method presented in Ref.~\cite{DTCexp1,DTCexp7}.  Alternatively, the emerging Sycamore processor \cite{DTCqs2,Sycamore} represents another promising platform to simulate various properties of DTCs or their related models. In particular, the ZZ interaction required to implement Eq.~(\ref{sys}) can be implemented using the native gates within the Sycamore processor \cite{DTCqs2}.

This work opens up interesting avenues for future work. First, following the intricate interplay between MPMs and subharmonic response unveiled in this work, it is of interest to investigate subharmonic response in a class of spin systems associated with a type of periodically driven topological superconductor, whose number of MPMs is tunable by some system parameters \cite{FMF6,FMF7,FMF8}. Second, exploring the effect of increasing the number of connectivities for each spin is another direction worth pursuing. In particular, with each spin possessing either three or four connectivities (depending on the boundary conditions) in the model considered in this paper, accidental degeneracy involving states differing by two spin excitations may still remain. It is expected that a modification to Eq.~(\ref{sys}) by adding more spin connectivities will further remove this less dangerous degeneracy and hence may have significant impact on the dynamical features.  Finally, our main construction of disorder-free system above can be naturally extended to simulating other types of subharmonic response characterized by $nT$-periodic observables \cite{repDTC}, where $n\in \mathbb{Z}$.

\vspace{0.3cm}

	\begin{acknowledgements}
		{\bf Acknowledgement}: R.W.B. is supported by the Australian Research Council Centre of Excellence for Engineered Quantum Systems (EQUS, CE170100009).  R.W.B. thanks Isaac Kim and Adrian Chapman for useful discussions.  J.G. acknowledges fund support by the Singapore NRF Grant No. NRFNRFI2017-04 (WBS No. R-144-000-378-281).
	\end{acknowledgements}


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

The setting we consider is the two-arm non-sequential randomized study with a continuous endpoint (e.g. a pill-placebo double-blind clinical trial assessing blood glucose improvement) that seeks inference for an additive treatment effect. There are $2n$ individuals which are sample-size balanced between the two groups. Each subject $i$ is placed into the treatment or control group and this information is encoded in $w_i \in \braces{-1, +1}$ where -1 denotes assignment to control and +1 denotes assignment to treatment. The vector $\bv{w} := \bracks{w_1~\ldots~w_{2n}}$ is called an \textit{assignment}, \textit{allocation}, or a \textit{randomization}. The collection of the legal $\bv{w}$ vectors and the probability distribution on the legal vectors in the experimental setting is known as an experimental \textit{design}, \textit{strategy}, \textit{algorithm} or \textit{procedure}.

There are $p$ covariates with the $j$th's value for subject $i$ denoted $x_{ij}$ that are assumed known in advance of the randomization (the setting of known measurements is sometimes called \textit{offline} to distinguish it from the \textit{sequential} setting, the latter not being the setting considered in this paper). After any randomization, the values of the subjects' covariates in each group are approximately the same; and it is this fact that is the main reason randomization is employed when seeking causal inference \citep{Cornfield1959}. Thus, \textit{imbalance} (of which there are many principled metrics) in the covariate values between the treatment subjects and the control subjects should be small. Upon completion, the sample responses $y_i$, the result of an unknown process we call the \emph{response model}, are assessed. These sample responses are used to compute a causal estimate of the population average treatment effect $\beta_T$.

Our contribution is a new experimental design that provides (1) lower error in the measurement of the average treatment effect and (2) is robust when the response model is linear and/or non-linear. If the response model is linear, the best experimental design is one that optimizes the treatment assignment for minimal covariate imbalance which we call \emph{imbalance-optimizing designs}. If the response model is nonlinear, a good strategy is to optimize the treatment assignment via creating binary matches (i.e. nonbipartite pairings of subjects) with small intramatch covariate distance. If the response model is a combination of both linear and nonlinear, then both criteria (minimize covariate imbalance and matching) become important. Our hybrid design does both: first by creating optimal nonbipartite matches on covariate distance and then optimizing intramatch assignments to provide small covariate imbalance across the whole sample. The latter can be accomplished by many methods; we investigate rerandomization and a pairwise greedy-switching heuristic (the subject of our previous work) where this heuristic now operates on the differences of the covariate values within each match in contrast to the raw subject covariate values.

To illustrate the advantages of our hybrid experimental design that both matches and imbalance-optimizes, consider linear and nonlinear response models, respectively $y_i = x_i + \frac{1}{2}\beta_T w_i + \mathcal{E}_i$ and $y = x^2 + \frac{1}{2} \beta_T w_i + \mathcal{E}$ where $\beta_T$ we set to 1, $x_i$'s are drawn iid from $U(0, 3)$ and the $\mathcal{E}_i$'s are drawn iid from $\normnot{0}{0.1^2}$. We use a small sample size of $2n = 40$ and consider three designs: (G) pure imbalance optimization via the greedy pair switching design which will be detailed in Section~\ref{sec:methodology}, (M) optimal non-bipartite matching and (MG) optimal non-bipartite matching followed by imbalance optimization via the greedy pair switching design. We simulate the above system 1,000 times and measure covariate imbalance using the metric log$_{10}(|\bar{\bv{x}}_T - \bar{\bv{x}}_C|)$ and squared error between the estimate $\ybar_T - \ybar_C$ and $\beta_T$. We average over the simulations to estimate both means and display the results in Table~\ref{tab:basic_demo}.

\begin{table}[h]
\centering
\begin{tabular}{cccc|ccc}
&&Average & Average & \multicolumn{3}{c}{p value for the average} \\
Response && Log$_{10}$  & Squared & \multicolumn{3}{c}{squared error comparison} \\
 Model & Design & Imbalance  &  Error & G & M & MG \\ \hline
Linear & G  & -3.90 & 0.00099 & - & $\approx 0$ & 0.74  \\
Linear & M  & -1.87 & 0.00135 & - & - & $\approx 0$  \\
Linear & MG  &  -4.55 & 0.00105 & - & - & - \\ \hline
Nonlinear & G &  - &0.04350 & - & $\approx 0$  & $\approx 0$ \\
Nonlinear & M &  - &0.00614 & - & - & 0.08\\
Nonlinear & MG &  - & 0.00273 & - & - & - \\
\end{tabular}
\caption{Results of the illustrative simulation explained in the text.}
\label{tab:basic_demo}
\end{table}

In the linear model, the pure imbalance-optimizing design (G) provides better MSE performance than matching (M) and employing the hybrid design combining matching with imbalance-optimization (MG) does not penalize performance (i.e. G and MG have MSE's that are statistically equal). When the model is non-linear, the covariate imbalance is less relevant, and  thus (M) provides much better MSE performance than (G) and the hybrid design (MG) outperforms (M) although not statistically significantly here (we will see in the results section that this advantage is likely only for very small sample sizes). Thus the hybrid (MG) design is robust and adaptive as it has the best performance in both linear and nonlinear settings.

When examining the imbalance, M achieves small covariate imbalance but not when compared to G / MG which is expected as that is not the former's primary objective in contrast to the latters' primary objective. MG has smaller covariate imbalance than G as it can combine the reduction in imbalance from M with the largely independent reduction in imbalance from G, the pure imbalance-optimizing procedure. The nearly order of magnitude difference in imbalance between MG and G does not translate to a significantly lower MSE as there is diminishing returns for imbalance optimization; the theoretical reasons for this will become clear in the following section.

\section{Methodology}\label{sec:methodology}

\subsection{Setup and Assumptions}\label{subsec:setup}

We denote the responses $\bv{y} = \bracks{y_1, \ldots, y_{2n}}^\top$ where the number of subjects $n$ is assumed even. The \emph{assignment} or \emph{allocation} vector is $\bv{w} = \bracks{w_1, \ldots, w_n}^\top$ whose entries are either +1 (the subject received $T$) or -1 (the subject received a $C$) and $\bv{w} \in \braces{-1,+1}^{2n}$. We define a \emph{design} $D$ as $W_D$, a discrete uniform random variable with support $\mathbb{W}_D \subseteq \braces{-1,+1}^{2n}$ and number of vectors $\abss{\mathbb{W}_D} \leq 2^{2n}$. We restrict the designs we consider to those that have the mirror property \citet[Assumption 2.2]{Kapelner2020}, \citet[Assumption 1b]{Kallus2018}, \citet{Nordin2020} i.e. for any $\bv{w}$, if treatment assignments and control assignments were switched then that resulting assignment would also be in the design: $\bv{w} \in \mathbb{W} \Rightarrow \bv{1} - \bv{w} \in \mathbb{W}$. We further restrict the designs we consider to be \emph{forced balance procedures} where all allocations have the same number of treated and control subjects \citep[Chapter 3.3]{Rosenberger2016}. This is a minor restriction is a standard design known as the balanced complete randomization design (BCRD) or completely random design, denoted $\mathbb{W}_{BCRD} = \braces{\bv{w}\,:\, \sum w_i = 0}$ which has $\binom{2n}{n}$ assignments. If a design further restricted this set, i.e. $\mathbb{W}_D \subset \mathbb{W}_{BCRD}$, we call $D$ a \emph{restricted design}.

We assume the only source of the randomness in the response is thus the treatment assignments $\bv{w}$. This assumption on the source of randomness is termed the \emph{randomization model} \citep[Chapter 6.3]{Rosenberger2016}, the \textit{Fisher model} or the \textit{Neyman model} whereby \qu{the $2n$ subjects are the population of interest; they are not assumed to be randomly drawn from a superpopulation} \citep[page 297]{Lin2013}.

As in \citet{Kallus2018}, we employ the simple \emph{differences-in-means} estimator,

\vspace{-0.25cm}\begin{eqnarray}\label{eq:estimator}
\hat{\beta} := \frac{\bv{w}^\top \bv{y}}{2n} = \frac{1}{2}(\Ybar_T - \Ybar_C)
\end{eqnarray}

which is population average treatment estimator regardless of the response model and the effect of the treatment \citep[Example 4.7]{Lehmann1998}.

\subsection{Previous Literature}

Our design idea stems from a combination of two long-standing streams of research: first, the pure-imbalance optimizing design literature and second, the binary-matching bias-reducing literature.

For the pure-imbalance optimizing design, we first reference the rerandomization design R. Here, one begins with vectors from BCRD and discards those whose covariate imbalance is beyond a desired threshold. This idea that goes back to Fisher, the father of randomized experiments who was aware of imbalanced allocations and ironically warned against pure unrestricted randomization, the procedure that he is famous for.

Here, $\mathbb{W}_{R(a)} := \braces{\bv{w}\,:\, d_M(\bv{w}) \leq a, \bv{w} \in \mathbb{W}_{BCRD}}$ where

\vspace{-0.25cm}\begin{eqnarray}\label{eq:mahal}
d_M(\bv{w}) := n (\bar{\bv{x}}_T - \bar{\bv{x}}_C)^\top \hat{\Sigma}_X^{-1} (\bar{\bv{x}}_T - \bar{\bv{x}}_C),
\end{eqnarray}

a popular metric of covariate imbalance among the treatment and control groups known as the Mahalanobis distance, $a$ is an upper-bound acceptability threshold of the covariate imbalance, $\bar{\bv{x}}_T$ and $\bar{\bv{x}}_C$ are the averages across the $p$ covariates in the treatment and control groups respectively and $\hat{\Sigma}_X$ is the sample variance-covariance matrix of all $2n$ subjects. Since $a < \infty$ implies that $\mathbb{W}_{R(a)} \subset \mathbb{W}_{BCRD}$, rerandomization is a restricted design by construction. Assuming multivariate-normally distributed covariates and an additive treatment effect but no assumption on the nature of the response model, the multiplicative reduction in MSE of the difference-in-averages treatment estimator when using the R design is $\eta(p, a) R^2$ where

\vspace{-0.25cm}\begin{eqnarray}\label{eq:rerand_reduction}
\eta(p, a) := 1 - \frac{2}{p}  \frac{\gamma(p/2 + 1, a/2)}{\gamma(p/2, a/2)},
\end{eqnarray}

$\gamma(\cdot, \cdot)$ denotes the lower incomplete gamma function and $R^2$ is variance explained in an OLS linear regression of the response on the covariates \citep{Morgan2012}.

This reduction in treatment estimator error comes from two sources (1) a reduction of the imbalance in the $p$ covariates which is at their control and (2) the degree of linearity of the covariates in the response as measured by $R^2$.

As for their first source of estimator error reduction, lowering the imbalance in rerandomization is only slightly more impressive than BCRD. There are many ideas for achieving better reduction in covariate imbalance. First, \citet[Section 3.3]{Kallus2018} conjectured using a heuristic argument that the optimal reduction (i.e. for the best possible vector) is exponential $2^{-\Omega(n)}$. Numerical optimization such as in \citet{Bertsimas2015} employ heuristics to approximate the optimal vector. Some heuristics come with theoretical guarantees, e.g. the greedy pair-switching design of \citet{Krieger2019} (henceforth denoted KAK19) uses a greedy heuristic that begins with a draw from BCRD, considers switching every T/C pair and retains the one with the greatest reduction in the imbalance metric and stops switching when it is no longer fruitful. The resulting vectors have a provably very low imbalance $O_p(n^{-(1 + 2/p)})$ which can then be further enhanced by generating many $\bv{w}$'s and retaining only the best of those akin to rerandomization. Moreover, since there is a provably low number of switches, the degree of randomness is nearly that of BCRD. We denote this restricted design as G.

Better rates than offered by G are unlikely to matter as the imbalance reduction multiple (Equation~\ref{eq:rerand_reduction}) has diminishing returns in $a$. After a procedure such as KAK19, the performance of the estimator is much more sensitive to the second source of error, i.e. the strength of the covariates' affect on the response (as gauged by $R^2$). If the $R^2$ is low (in the case of a non-linear response model and/or a high noise setting), covariate imbalance reduction will not be too fruitful. Moreover, \citet[Section 2.3.3]{Kallus2018} proves that minimizing $a$ in rerandomization is fruitful for minimax estimator error reduction only in the case of a linear response model. \citet{Kallus2018} explains the reason why this advantage is only limited to linear response models: even if $d_M(\bv{w}) = 0$ (i.e. the first moments are matched perfectly), the covariate distributions in the two arms could be very different. In a nonlinear response model, the estimator performance can depend on more subtle differences in the covariates' distributions (e.g. their tail behavior).

The second stream of research, binary matching, implicitly attempts to equalize the covariate distributions in both arms. Here, the $2n$ subjects' indices are first organized into a \textit{matched pair structure} $\mathcal{M}_d$, a set whose elements are \textit{pairs} (a set of two subject indices). The structure $\mathcal{M}_d$ is created by first assuming a distance function $d(\bv{x}_r, \bv{x}_s) : \mathbb{R}^p \times \mathbb{R}^p \rightarrow \mathbb{R}_{\geq 0}$ whose inputs are subject covariates. A $2n \times 2n$ nonnegative symmetric matrix is then computed consisting of distances for all pairs of subjects. The optimal match structure $\mathcal{M}_d^\star$ would be one of the $\binom{2n}{n}$ sets that minimizes the sum of the resulting intramatch distances i.e. $\sum_{\braces{r,s} \in \mathcal{M}} d(\bv{x}_r, \bv{x}_s)$. The algorithm that produces the optimal match structure is known as \emph{optimal nonbipartite matching}. This algorithm has been shown to be reducible to a polynomial-time algorithm (see \citealp{Lu2011} for details and history). To create an assignment vector $\bv{w}$, the individual assignments within each of the $n$ pairs are randomly allocated to T/C or C/T via $n$ iid Bernoulli draws. We denote the optimal binary matching design with the Mahalanobis distance metric (Equation~\ref{eq:mahal}) as M and its allocation space by $\mathbb{W}_{M} = \{\bv{w}\,:\,w_r = -w_s, \braces{r,s} \in \mathcal{M}_{d_M}^\star\} \subset \mathbb{W}_{BCRD}$.

Can binary matching play a role in experimental design? This M design was investigated by \citet{Greevy2004} who reported better imbalance and higher power using randomization tests. However, matching in order to have robust estimation under a nonlinear response model has a long literature in observational studies and we recommend \citet{Stuart2010} for a broad overview. (Without the ability to manipulate $\bv{w}$, the matching procedure is quite different as it is limited to fixed $w_i$'s and is called \emph{bipartite matching}). \citet{Rubin1979} was the first to show this robustness and recommended that $d$ be specified as the Mahalanobis distance of Equation~\ref{eq:mahal}.

Further, \citet[Section 2.3.2]{Kallus2018} proves that binary matching via $d$ is the minimax variance minimizing experimental design for response functions that are Lipschitz continuous with respect to metric $d$. It is reasonable to conclude in a real-world experimental setting (for instance in a clinical trial in medicine) that the response would be continuous with derivatives not changing too quickly (i.e. satisfying the Lipschitz conditions). Thus, to reduce the variance of the treatment estimator in nonlinear response models, we elect to use a binary matching design.

Our contribution is relatively simple: we first compute $\mathcal{M}_{d_M}^\star$ as part of an M design and this binary pairing structure will provide robustness to nonlinear response models. Then instead of assigning the specific pairs' subjects via Bernoulli draws as explained above, we further restrict the pair assignments to insist on small imbalance by using R (resulting in the matching-then-rerandomization design we term MR) or by using G (resulting the matching-then-greedy-pair-switching design we term MG).  We turn to the details of MR and MG now.

\subsection{Our Algorithms}

First, $\mathcal{M}_{d_M}^\star$ is computed by running the optimal non-bipartite matching algorithm implemented by \citet{Beck2016} in the \texttt{R} package \texttt{nbpMatching}.

For R after M, $\mathbb{W}_{MR(a)} := \braces{\bv{w}\,:\, d_M(\bv{w}) \leq a, \bv{w} \in \mathbb{W}_M}$ i.e. we draw many different assignments from $\mathbb{W}_M$ and retain those whose Mahalanobis distances meet the threshold $a$ thereby finding small imbalances subject to the match structure.

For G after M, we cannot simply input $\bv{w} \in \mathbb{W}_M$ directly into KAK19 as this would violate the $\mathcal{M}_{d_M}^\star$ structure, breaking the binary pairings. Instead, we temporarily view the matched pairs as individual subjects (of which there are $n$) as the input into KAK19. This means we temporarily view the intramatch covariate vector differences as the subjects' covariates which are inputted into KAK19. Put another way, the T arm now refers to a matched pair assigned T/C where the covariate difference vector is computed by taking the first subject's covariates minus the second subject's covariates. The C arm would be the opposite. In one iteration of KAK19, all \textit{pair of pairs switches} are considered: from (T/C)/(C/T) to (C/T)/(T/C) or vice versa. One of these will result in the largest reduction in imbalance as measured among the pair of pairs which is algebraically the same imbalance as measured among the entire $\bv{w}$ vector for the $2n$ subjects. This best pair of pairs switch is retained and the algorithm stops when an iteration cannot reduce the imbalance any longer. This scheme thereby preserves $\mathcal{M}_{d_M}^\star$.

Under our designs, to evaluate hypotheses concerning the value of $\beta_T$ we can use the randomization test and to construct approximate confidence intervals for $\beta_T$ we can invert the randomization test. For details, see \citet[Section 2.2]{Morgan2012} and KAK19, Section 6; both provide a careful step-by-step explanation and a literature review.

One can consider an alternative implementation of the G in MG as follows. After M, each iteration of G can perform one switch of one pair (instead of a switch of a pair of pairs). Here, we now have fewer switch possibilities, i.e. $n$ and not $\binom{n}{2}$, and we conjecture this procedure would not perform as well as the one proposed above.

One can also consider a GM (or RM) procedure i.e. to first imbalance-optimize via G (or R) and then employ binary matching. In that scenario, there will be a different $\mathcal{M}$ for each initial imbalance optimized vector based on a bipartite matching. The gain here would be that the imbalance-optimization would be slightly better but it comes at a cost: bipartite matches are much more constrained than nonbipartite matches and thus we conjecture that these designs will have much lower robustness to nonlinear models.

\subsection{Theoretical Investigation of MG and MR}\label{sec:theory}

\subsubsection{Imbalance}

The key idea of our designs is to retain the match structure while optimizing imbalance. Instead of employing G or R using $\mathbb{W}_{BCRD}$ as a base, we now employ G or R using $\mathbb{W}_{M}$ as a base. We are not aware of theoretical results for the imbalance after optimal nonbipartite matching when $p > 1$ and thus we focus on the case when $p = 1$.

In this univariate case, consider the order statistics $X_{(1)}, \ldots, X_{(2n)}$. Optimal matching creates $\mathcal{M} = \braces{\braces{2n, 2n-1}, \ldots, \braces{2,1}}$, i.e. pairing the first largest covariate with the second largest covariate value, then pairing the third largest with the fourth largest, etc. We will consider the case where $X_1,\ldots,X_{2n}$ are iid $\uniform{0}{1}$. Define $\Delta_i = X_{(2i)} - X_{(2i-1)}$ for $i = 1, 2, \ldots, n$. By \citet[Theorem 6.6c]{Dasgupta2011}, $\bracks{\Delta_1,\ldots,\Delta_n}$ has the same distribution as $\oneover{S}\bracks{E_1,\ldots,E_n}$, where $E_1,\ldots,E_n,\ldots,E_{2n+1}$ are iid $\exponential{1}$ and $S := \sum_{i=1}^{2n+1} E_i = O_p(n)$. Since the denominator $S$ is common to all entries of the vector, applying the greedy pair-switching of KAK19 to $\bracks{\Delta_1,\ldots,\Delta_n}$ is the same as (in terms of distribution) applying the algorithm to $\bracks{E_1,\ldots,E_n}$ and then dividing by $S$.

The imbalance of assignments in G is $O_p(n^{-3})$ by KAK19 (Theorem 1), a distribution-free result. When dividing by $S$, which is $O_p(n)$, the resultant rate of our MG design is $O_p(n^{-4})$.

The imbalance of assignments in R is less than $a$ by construction. However, if $a$ is too small, then practical computational constraints dictate that all BCRD assignments one generates will have imbalance greater than $a$, resulting in the inability to locate any assignments in $\mathbb{W}_{R(a)}$. KAK19 reframes the discussion of R's imbalance by making this computational constraint explicit. Let $n_R$ denote the number of BCRD assignments one searches through for each resultant assignment in R. For every collection of $n_R$, the assignment that provides the minimum imbalance in that collection is retained. These retained vectors then have imbalance which is order $n_R^{-1} O_p(n^{-1/2})$ where $O_p(n^{-1/2})$ is the imbalance of BCRD. It then follows by the same argument above that when dividing by $S$, the resultant imbalance of MR is then $n_R^{-1} O_p(n^{-3/2})$.

Extending the theoretical results for imbalance in MG beyond the assumption of uniformly distributed covariates is difficult. The main result needed to prove that the order of the imbalance of the assignments in G is $O_p(n^{-3})$ is to show that there exists a pair of covariate values to switch $X_i, X_j$ where $w_i = 1$ and $w_j = -1$ that result in a necessary reduction in imbalance. The analogous result here would require the locating of $\Delta_i, \Delta_j$ where $w_i = 1$, $w_{i+1} = -1$, $w_j = -1$ and $w_{j+1} = 1$ that results in the same reduction in imbalance after switching both pairs. Showing the necessary reduction in imbalance assumed that the $X_i$'s are continous and iid. We were able to use this result in this paper when the $X_i$'s were uniformly distributed, but if the $X_i$'s are normally distributed, it is much more complicated as the $\Delta_i$'s corresponding to covariate values in a tail are stochastically larger than the $\Delta_i$'s corresponding to covariate values in the center. The imbalance of MG in the normal case is shown by simulation to be lower than G (see Figure~\ref{fig:mse_density_plots_normal} in the Supporting Information).

\subsubsection{Degree of Randomness}

We now turn to degree of randomness in MG and MR. We reiterate that in experimental design, it is important for the assignments to be highly random as this randomness provides insurance against a large variance in the treatment effect estimator due to unobserved covariates \citep{Kapelner2020}. Both G and R are less random than BCRD since they are forms of restricted randomization. In the case of G, KAK19 (Theorem 2 and Proposition 1) shows that G's assignments are asymptotically as random as BCRD according to three degree-of-randomness metrics: (a) the pairwise entropy metric of assignments, (b) the standard error metric of probability of pairwise assignments and (c) the accidental bias metric of \citet[Section 5]{Efron1971} defined as the maximum eigenvalue of the variance covariance matrix of $W$, the random variable producing assignments in any strategy (see Section~\ref{subsec:setup}). The reason for G's high degree of randomness is that the number of pairwise switches made in G is small, being only $O_p(\sqrt{n})$. Thus, since G starts with BCRD, and its assignments are not changed significantly, the resulting $\mathbb{W}_G$ is very large. In the case of R, we directly limit the number of designs based on the cutoff $a$. As, $a$ decreases, $\mathbb{W}_R$ shrinks relative to $\mathbb{W}_{BCRD}$, making R less random. Since $a$ can never be too small due to computational limitations, it follows that randomness is not substantially deteriorated in R. Thus, MG is as random as M asymptotically and MR-in-practice is likely as random as M asymptotically as well. And M is highly random according to all three degree-of-randomness metrics, only slightly less random than BCRD.

\section{Simulation Results}\label{sec:simulations}

To illustrate our main results we turn to simulations. All simulations herein were performed with \texttt{GreedyExperimentalDesign}, an \texttt{R} package available on \texttt{CRAN} whose core is implemented in \texttt{Java} for speed. We measure mean squared error (MSE) of the difference-in-averages estimator for an additive average treatment effect when the response model depends on underlying covariate(s).

We conjecture that (a) the MSE performance of MR/MG will be similar to imbalance-optimizing designs G/R when the response model is purely linear, (b) the MSE performance of MR/MG will be similar to the pure nonbipartite matching design M when the response model is purely nonlinear and (c) the MSE performance of MR/MG will be better than either exclusively imbalance-optimizing G/R or exclusively nonbipartite matching M when the response model is a hybrid of a strong linear component and strong nonlinear component.

To demonstrate the conjecture we consider $n=100$, two covariates and five response models that take the form $\bv{Y} = \beta_T\bv{w} + f(x_1, x_2) + \bv{\errorrv}$ where $x_1, x_2$ are independent draws from $U(-\sqrt{3},+\sqrt{3})$ so that their variance is one and the $n$ components of $\bv{\errorrv}$ are iid draws from $\normnot{0}{0.5^2}$. The five models are (Z) a null zero model $f = 0$ for calibration and ensuring the simulation is working (L) a purely linear model $f = 3x_1 + 3x_2$, (LsNL) a mostly linear but slightly nonlinear model $f = 3x_1 + 3x_2 + x_1^2$, (LNL) a linear and nonlinear model $f = 3x_1 + 3x_2 + x_1^2 + x_2^2 + x_1 x_2$ and (NL) a purely nonlinear model $f = x_1^2 + x_2^2 + x_1 x_2$.

We then consider the six designs of the previous section all of which feature an equiprobable draw from a set of mirrored allocations: (BCRD) balanced and completely randomized, (R) rerandomization of \citet{Morgan2012} where only the top 1\% of BCRD vectors are retained, (G) the greedy pair switching of KAK19, (M) optimal nonbipartite matching, and our two new hybrid approaches of (MG) optimal nonbipartite matching followed by greedy pair switching on the within-pair covariate differences and (MR) optimal nonbipartite matching followed by rerandomization on the within-pair covariate differences where the top 1\% of pairings are retained. Thus, we have $5 \times 6$ = 30 sample size-model-design simulation settings.

For each sample size-model-design simulation setting, we simulate 50 different covariate and noise value settings, $x_1, x_2, \mathcal{E}_1, \ldots, \mathcal{E}_n$. Within each specific covariate-noise realization setting, we draw 500 unique allocation vectors from each of the six design strategies. We define MSE as $\expesubnostr{X_1, X_2, \bv{\errorrv}}{\cexpesubnostr{W_D}{(\hat{\beta}_T - \beta_T)^2}{X_1, X_2}}$ for each sample size-model-design setting. This MSE is estimated by first taking an arithmetic average of the 500 squared errors within the covariate-noise setting and then averaging these over the 50 covariate-noise settings. Estimates for matching designs (M and MR) usually are computed by the average paired differences. This estimate is algebraically equivalent to the naive difference-in-averages estimate and thus estimates are computed uniformly for all six designs.

\begin{figure}[h]
\centering
\includegraphics[width=6.5in]{fig1.pdf}
\caption{Density plots for MSE (on a log scale) over the covariate distributions, noise and assignments when $n=100$ for the five models (the vertically stacked panes) and six designs (depicted as different colors) explained in the text. The vertical lines indicate the MSE estimate (explained in the text) by design. Lines directly adjacent to one another are the result of a slight jitter. Statistically significant differences are tabulated in Table~\ref{tab:all_pvals}.}
\label{fig:mse_density_plots}
\end{figure}

Figure~\ref{fig:mse_density_plots} provides density plots for MSE when $n=100$ for all five models and six designs and Table~\ref{tab:all_pvals} tabulates the statistically significance of each design comparison. The null response model (Z) serves as a check on the integrity of our simulation: since the covariates are disjointed from the response, all designs perform equally as seen in the top pane of Figure~\ref{fig:mse_density_plots} and for this reason it was omitted from Table~\ref{tab:all_pvals}. For all other response models where the covariates inform the response, the figure and tables illustrates the main thrust of this paper's contribution: (1) Our designs, MR and MG, are always the best performers (or tied for the best performance) indicating robustness to linear, nonlinear and mixtures of linear and nonlinear response models. Other observations support this main thrust. (2) BCRD is the worst performer (or tied for worst performance) which is expected given that it does not allocate subjects with the benefit of knowing the covariate values. (3) In the exclusively linear model (L), the pure imbalance optimizing designs perform the best i.e. G, R, MR, MG all perform equally with M lagging. (4) In the exclusively nonlinear model (NL), the matching designs perform the best i.e. M, MR, MG all perform equally. (5) In the linear-nonlinear models (LsNL, LNL), our matching-plus-imbalance-optimizing designs MR and MG are equally the best, the matching M is second best and the pure imbalance optimizing designs G and R lag behind M. (6) R lags behind G because its imbalance is a much smaller rate (see KAK19 for details why). MR also lags behind MG for the same reason but the difference is not appreciable in these simulations.

\begin{table}[h]
\centering

\begin{subtable}{.5\linewidth}{
\begin{tabular}{rlllll}
 & R & G & M & MR & MG \\
  \hline
BCRD & *** & *** & *** & *** & *** \\
  R &  &  &  &  &  \\
  G &  &  & ** &  &  \\
  M &  &  &  & ** & ** \\
  MR &  &  &  &  &  \\
\end{tabular}
}
\caption{Model: L}
\label{tab:pvals1}
\end{subtable}~
\begin{subtable}{.5\linewidth}{
\begin{tabular}{rlllll}
 & R & G & M & MR & MG \\
  \hline
BCRD & *** & *** & *** & *** & *** \\
  R &  &  & *** & *** & *** \\
  G &  &  & ** & *** & *** \\
  M &  &  &  & * & * \\
  MR &  &  &  &  &  \\
\end{tabular}
}
\caption{Model: LsNL}
\label{tab:pvals2}
\end{subtable}\\
~\\~\\

\begin{subtable}{.5\linewidth}{
\begin{tabular}{rlllll}
 & R & G & M & MR & MG \\
  \hline
BCRD & *** & *** & *** & *** & *** \\
  R &  &  & *** & *** & *** \\
  G &  &  & *** & *** & *** \\
  M &  &  &  & * & * \\
  MR &  &  &  &  &  \\
\end{tabular}
}
\caption{Model: LNL}
\label{tab:pvals3}
\end{subtable}~
\begin{subtable}{.5\linewidth}{
\begin{tabular}{rlllll}
 & R & G & M & MR & MG \\
  \hline
BCRD &  & *** & *** & *** & *** \\
  R &  & *** & *** & *** & *** \\
  G &  &  & *** & *** & *** \\
  M &  &  &  &  &  \\
  MR &  &  &  &  &  \\
\end{tabular}
}
\caption{Model: NL}
\label{tab:pvals4}
\end{subtable}

\caption{Tukey-Kramer comparisons of the MSE for each design pair by model. * indicates $p < 0.05$, ** indicates $p < 0.01$ and *** indicates $p < 0.001$. The Z model is not displayed as none of its design pairs' MSEs tested statistically significant.}
\label{tab:all_pvals}
\end{table}

We also repeated the simulation for three other sample sizes ($n=32$, $n=132$ and $n=200$) and they all have qualitatively similar results. The tabulation of MSE estimates for all sample sizes, designs and response models can be found in Tables~\ref{fig:ols_n_32}-\ref{fig:ols_n_200} in the Supporting Information. The entire simulation was also repeated for the case of normally distributed covariates and the results were not substantially different (see Figure~\ref{fig:mse_density_plots_normal} in the Supporting Information).

\section{Concluding Remarks}\label{sec:conclusion}

We proposed a randomized experimental design that is a hybrid between two well studied strategies: those that optimize covariate imbalance and those that use binary matching. The former is known to be minimax optimal when the response model is linear in the covariates and the latter is known to be minimax optimal when the response model is continuous with limits on its derivative with respect to all $p$ covariates. Additionally, we have theory that shows that the fusing of both strategies fortuitously enhances the covariate imbalance all while retaining a degree of randomization close to that of the classic completely random designs.

Because we fuse both types of design together without sacrificing the performance of either, our designs provides very low MSE performance when estimating a population average treatment effect in the purely linear case (because we optimize covariate imbalance), the purely nonlinear case (as we our assignments are optimal nonbipartite matched pairs) and response models that feature both linear and nonlinear components. Thus we expect our design to be very powerful when employed in real-world settings such as clinical trials and Internet-based experimentation.

There are many extensions to this work. First, there are different types of matching such as ratio matching \citep[Section 3.1.2]{Stuart2010}, matching in more than two arms (ibid, Section 6.1.4), matching with discrete covariates (the Mahalanobis distance function is not the most appropriate here as discussed in \citealt{Gu1993}) and matching with some covariates being more important than others. As some matches may not be acceptable, one can caliper match; this would result in a subset of the subjects being matched and a subset being unmatched. In this case, the G algorithm would have to be redesigned but the R procedure would be fairly straightforward. Further, one may wish to use the OLS estimator instead of the classical difference-in-averages estimator (Equation~\ref{eq:estimator}). One can argue that restricted randomization can leave you vulnerable to high MSE estimation due to imbalance in unobserved covariates \citep[Sections 2.2.1, 2.2.3, 2.3.1 and 2.3.3]{Kapelner2020}. Therein, they proved that restricted designs with a high degree of randomness, such as our hybrid designs in this work (see Section~\ref{sec:theory}), are highly random and thus are not susceptible. Further, \citet[Equation 14]{Kapelner2020} showed that the OLS estimator improves asymptotic MSE by an order of magnitude in sample size. However, imbalance-optimizing designs (such as R and G) still provide improvement in the linear case. Covariate adjustment cannot aid with the nonlinear component of the response model and thus our hybrid designs will perform very well in conjunction with the OLS estimator on nonlinear models. And the OLS estimator is actually preferred in matching designs \citep{Rubin1979}. Future work can also adapt this design to the sequential setting, where individuals in the experiment arrive one by one and must be assigned quickly to an arm. We can apply the on-the-fly matching procedure in \citet{Kapelner2014} with the rerandomization approach in \citet{Zhou2018} in the same vein as our hybrid procedure.

\subsection*{Replication}

All figures and tables can be reproduced by running the \texttt{R} code found at \url{https://github.com/kapelner/GreedyExperimentalDesign/blob/master/hybrid_paper}.

\subsection*{Acknowledgements}

This research was supported by Grant No 2018112 from the United States-Israel Binational Science Foundation (BSF).

\bibliographystyle{apa}