\section{Introduction}
Recent advances in quantum technology have attracted widespread attention in various disciplines. Electronic structure theory is one of the fields expected to potentially benefit from the development of quantum computers because strong correlation can be handled by mapping the wave function directly onto entangled qubits without an exponential increase in computational cost. However, current noisy intermediate-scale quantum (NISQ) devices cannot implement error correction owing to the limited quantum resources; hence, they suffer from errors triggered by noise and short coherent times,\cite{Preskill18, Arute19, McArdle20} which limit the application of general quantum algorithms such as quantum phase estimation. Therefore, several hybrid quantum-classical algorithms have been proposed to address the technical challenges emerging in NISQ hardware. The variational quantum eigensolver (VQE) is one such algorithm,\cite{Peruzzo14, McClean16} in which a quantum state is propagated by a parametrized short quantum circuit that can be optimized in classical post processing, to variationally lower the energy expectation value and obtain the Hamiltonian ground state.

In quantum chemistry, the most conventional ansatz for VQE circuit is the unitary coupled-cluster (UCC),\cite{Kutzelnigg82, Kutzelnigg84, Bartlett89, Romero19, Sokolov20} where fermionic excitation operators with respective to the Fermi vacuum (e.g. Hartree-Fock) are anti-hermitized and exponentiated to form unitary gates. Given the remarkable success of standard CC in a wide range of classical applications, UCC with single and double substitutions (UCCSD) is expected to be a feasible starting point for chemical Hamiltonians. Nevertheless, UCCSD is manifestly incapable of describing strong correlations,\cite{Lee19} which is the primary target of quantum computing. Extending the excitation manifold to include general orbitals, generalized UCCSD captures higher excitation effects in a compact manner, and is thus known to often provide accurate results for strongly correlated systems.\cite{Wecker15, Lee19} However, these UCC-based methods are inhibited by a limited representability because of the fixed structure of unitary gates; in addition, they also require several parameters and deep circuits, which makes them difficult for NISQ devices to handle.

Ansatz-free algorithms are more flexible and can control the accuracy and circuit depth in a highly black-box manner.\cite{Grimsley19, Ryabinkin20, Tang21, Gomes21, Yordanov21} A prominent algorithm in this context is ADAPT-VQE proposed by Grimsley et al.,\cite{Grimsley19} which adaptively adds a set of short unitary gates to the quantum circuit successively. Each unitary is chosen from an operator pool comprising generalized-UCCSD excitation operators, based on the VQE energy gradient. Later, Tang et al. proposed qubit-ADAPT\cite{Tang21}, which employs Pauli rotations instead of a UCC-like operator to drastically reduce the gate depth. Recently, Yordanov et al. demonstrated that ``qubit-excitations,'' each of which is a linear combination of Pauli strings in qubit-ADAPT, can be implemented with only a handful of CNOT gates.\cite{Yordanov20} Based on this result, they employed qubit-excitations to form an operator pool in what they called qubit-excitation based (QEB) ADAPT.\cite{Yordanov21} However, these ADAPT-based methods break important symmetries of chemical Hamiltonians, of which the most significant is spin-symmetry.

It is well known in the quantum chemistry community that breaking spin-symmetry triggers several undesired consequences, which include inaccurate energies and properties, a slow convergence to the exact state, and a difficulty in interpreting obtained results.\cite{Handy85, Andrews91, Tsuchimochi10B} ADAPT-based methods result in spin-symmetry breakages because of the use of fermionic or qubit excitations, for which spin-adaptation is difficult when exponentiated. This is more the rule than the exception. In fact, in strongly correlated systems, spin-symmetry is completely broken because ADAPT typically requires large rotation parameters to represent the multi-configuration character. As will be demonstrated, this causes the slow convergence of ADAPT algorithms and inaccurate results with truncated operators.

To address this issue, this study employs a spin-projection circuit,\cite{Tsuchimochi20} which restores the total spin quantum number to restrict ADAPT and explore the operator candidate in the correct symmetry space. We demonstrate that our algorithm can generally reduce the circuit depth, although it increases the measurement cost. In particular, we show that the spin-projected ADAPT outperforms UCCSD as a better tradeoff between accuracy and circuit depth, which is not always guaranteed in ADAPT. To investigate the importance of spin-symmetry in molecular property calculations, we further derive the first-order energy derivative in the presence of the spin-projection operator.

The remainder of this paper is organized as follows. In Sections \ref{sec:FermionicADAPT} and \ref{sec:qubitADAPT}, we review the original algorithms for ADAPT-VQE, and describe how they break symmetries in Section \ref{sec:spin-contamination}. Section~\ref{sec:SP-ADAPT} presents detailed discussions on the application of spin-projection to each ADAPT scheme, including an explanation on how the simulated qubits can be tapered off in the presence of spin-projection. Subsequently, we derive the energy derivative for ADAPT in Section~\ref{sec:properties}, especially focusing on the importance of orbital response. In Section~\ref{sec:symmetry-breaking}, we demonstrate that the convergence in ADAPT becomes increasingly slower for strongly correlated systems, and in Section~\ref{sec:C}, we present our benchmark results on the spin-projected ADAPT, comparing the accuracy and gate efficiency between different methods. Section~\ref{sec:pools} demonstrates how the number of operators in a pool can be reduced in QEB ADAPT. The calculations for dipole moment and geometry are explained in Section~\ref{sec:property_results}. Finally, we conclude this work in Section~\ref{sec:conclusions}.

\section{Theory}\label{sec:Theory}
\subsection{Fermionic ADAPT-VQE}\label{sec:FermionicADAPT}

First, we review the ADAPT-VQE algorithm originally proposed by Grimsley et al.\cite{Grimsley19} This algorithm is based on the observation that the exact full configuration interaction (FCI) state can be reached from the Hartree-Fock (HF) state using the exponential form,
\begin{align}
	|\Psi^{\rm FCI}\rangle = \prod_k e^{t_k\hat \tau(k)} |\Phi^{\rm HF}\rangle \label{eq:FCI}
\end{align}
where $\hat\tau(k)$ is the $k$th instance of arbitrarily ordered operators comprising the anti-symmetrized single and double excitations,
\begin{subequations}
\begin{align}
&	\hat \tau_q^p = a^\dag_pa_q - a^\dag_q a_p\\
&	\hat \tau_{rs}^{pq} = a^\dag_p a^\dag_q a_r a_s - a^\dag_s a^\dag_r a_q a_p
\end{align}
\label{eq:tau}
\end{subequations}
Here, we have adopted $p>q,r>s$ to denote general spin-orbitals. Note that the same anti-symmetrized excitation operators can appear repeatedly in Eq. \ref{eq:FCI}, but with different variational parameters $t_k$. Although Eq.~(\ref{eq:FCI}) includes only singles and doubles, the effect of higher excitations such as triples and quadruples can be conveniently considered by their exponential products.

Regardless of its exactness, Eq.~\ref{eq:FCI} is extremely cumbersome and long because it requires a large number of $k$; hence, its realization is prohibitively difficult on quantum computers. On the one hand, from a circuit complexity perspective, such as the number of CNOT gates, one can employ only a limited number of unitaries $e^{t_k \hat \tau(k)}$ on NISQ computers. On the other hand, it is expected that such a truncated ansatz with a reasonable number $n$ of unitaries can still outperform fixed ans\"atze such as UCCSD, especially for strongly correlated systems, with  appropriately chosen $\hat \tau(k) \; (k=1,\cdots,n)$ and variationally optimized $t_k$. Therefore, the objective of ADAPT-VQE is to dynamically create an ansatz that approaches FCI, using a maximally compact sequence of $n$ unitary operators, which are successively determined based on the energy gradient. To do so, an operator pool ${\mathscr P}$ is first defined, which comprises $M$ operators that are used to construct the trial state:
\begin{align}
			{\mathscr P} = \{\hat{\tau}^{p_\alpha}_{q_\alpha}+\hat{ \tau}^{p_\beta}_{q_\beta}, \hat{\tau}^{p_\alpha q_\alpha}_{r_\alpha s_\alpha}+\hat{\tau}^{{p}_\beta {q}_\beta}_{{r}_\beta s_\beta} ,\hat{\tau}^{p_\alpha q_\beta}_{r_\alpha s_\beta}+\hat{\tau}^{p_\beta q_\alpha}_{r_\beta s_\alpha}\} \label{eq:Pool}
\end{align}
where we explicitly consider the spins $\alpha$ and $\beta$ and each linear combination of operators is labeled by $m=1,\cdots, M$. Note that each combination of doubles does not comprise the spin-complement form, which we would address later in this paper.

A trial state of ADAPT-VQE at the $n$th cycle is given by
\begin{align}
	\Ket{\psi_n}=e^{\theta_n \hat{\cal A}_n}e^{\theta_{n-1} \hat{\cal A}_{n-1}} \cdots e^{\theta_1 \hat{\cal A}_1}\Ket{\psi^{\rm HF}} \label{eq:ADAPT}
\end{align}
where $\hat {\cal A}_k \in {\mathscr P}$, and the parameters $\{\theta_k\}$ are variationally optimized to minimize the energy expectation value. To select the most contributing operator from the ${\mathscr P}$ for $n+1$th cycle, we create trial states with each $\hat A_m \in {\mathscr P} \; (m=1,\cdots, M)$ as
	\begin{align}
	|\psi_{n+1, m}\rangle = e^{\theta_{n+1} \hat A_m} |\psi_{n}\rangle
	\end{align}
and calculate the energy gradient for each $|\psi_{n+1,m}\rangle$ around $\theta_{n+1} = 0$,
\begin{align}
	R_{m}^{(n)} = \left. \frac{\partial \langle \psi_n | e^{-\theta \hat A_m} \hat H e^{\theta \hat A_m}| \psi_n \rangle}{\partial \theta}\right|_{\theta = 0} = \langle \psi_n | \left[\hat H, \hat A_m\right] | \psi_n \rangle\label{eq:Rm}
\end{align}
which can be evaluated with up to three-body reduced density matrix (3RDM) because of the commutator property.
The operator with the largest gradient $R_m^{(n)}$ is selected as $\hat {\cal A}_{n+1}$, with which the parameters $\{\theta_k: k=1,\cdots,n+1\}$ in $|\psi_{n+1}\rangle$ are all optimized by the standard VQE.

The ADAPT algorithm presented above is converged when the norm of Eq. (\ref{eq:Rm}) is smaller than the threshold $\epsilon$. It is worth noting that such a condition is closely related to the $k$-particle generalized Brillouin theorem.\cite{Kutzelnigg79} Kutzelnigg demonstrated that an exact FCI wave function satisfies $\langle \Psi_{\rm FCI}| \left[\hat H, \hat \tau_\mu\right] |\Psi_{\rm FCI}\rangle = 0$ for all excitation ranks. The convergence in ADAPT-VQE with $\epsilon=0$ only corresponds to the 1- and 2-particle generalized Brillouin conditions, and therefore, ADAPT-VQE is not guaranteed to converge to FCI. However, in practice, it approaches the FCI accuracy as $||{\bf R}||$ decreases.\cite{Claudino20}

\subsection{Qubit-based ADAPT-VQE}\label{sec:qubitADAPT}
The aforementioned Fermionic ADAPT is based on the Jordan-Wigner transformation for mapping fermionic excitation operators Eqs.~(\ref{eq:tau}) to Pauli operators. This introduces the chain of $Z$ operations because of the anti-symmetric character of fermions, e.g.,
\begin{align}
	\hat \tau_{rs}^{pq}&= \frac{i}{8} \big(Y_p X_q X_r X_s + X_p Y_q X_r X_s - X_p X_q Y_r X_s \nonumber\\
	& - X_p X_q X_r Y_s - X_p Y_q Y_r Y_s - Y_p X_q Y_r Y_s \nonumber\\
	&+ Y_p Y_q X_r Y_s + Y_p Y_q Y_r X_s\big) \bigotimes_{t=q+1}^{p-1} Z_{t} \bigotimes_{u=s+1}^{r-1} Z_{u}\label{eq:JW}
\end{align}
where we assume $p>q>r>s$, but similar relations can be obtained for other cases.
The exponential of $\hat \tau_{rs}^{pq}$ is usually implemented by decomposing the eight Pauli strings presented in Eq.~\ref{eq:JW}, each for which requires the ladders of CNOT gates to take into account the parities for $t = q+1, \cdots, p-1$ and $u = s+1, \cdots, r-1$; refer to Fig.~\ref{fig:circ} as an example. Therefore, such a na\"ive implementation of $e^{\theta \hat \tau_{rs}^{pq}}$ necessitates $8(p-q+r-s+2)\ge 48$ CNOTs. However, several studies \cite{Gomes20, Xia20, Zhang21} have demonstrated that even if all the $Z$ operators from the Pauli string in Eq.~(\ref{eq:JW}) are discarded to define the qubit excitation
\begin{align}
\hat{\tilde \tau}_{rs}^{pq} &= \frac{i}{8} \big(Y_p X_q X_r X_s + X_p Y_q X_r X_s - X_p X_q Y_r X_s \nonumber\\
	& - X_p X_q X_r Y_s - X_p Y_q Y_r Y_s - Y_p X_q Y_r Y_s \nonumber\\
	&+ Y_p Y_q X_r Y_s + Y_p Y_q Y_r X_s\big),\label{eq:tau_tilde}
\end{align}
the accuracy obtained by VQE will be similar to that with the full $\hat \tau_{rs}^{pq}$, although it eliminates a considerable number of CNOT gates.
Tang et al. attempted to further decompose $\hat{\tilde \tau}_{rs}^{pq}$ and construct a reduced operator pool for what they called qubit-ADAPT-VQE,\label{Tang21}, noting that each Pauli string has a similar action. Using circuits depicted in Fig.~\ref{fig:circ}(a), each operator $e^{\theta_m\hat A_m}$ requires only six CNOT gates (similarly, two CNOT gates for single excitation variants, $Y_p X_q$). However, as will be discussed, this treatment significantly breaks several symmetries in the Hamiltonian and slows down the algorithm’s convergence for chemical applications. Moreover, because it violates the number symmetry, it cannot be applied to ionized states or half-filled Hubbard models.

Recently, Yordanov et al. demonstrated that the role of $e^{\theta \hat {\tilde \tau}_{rs}^{pq}}$ is to rotate the bit strings between $|0_p0_q1_r1_s\rangle$ and $|1_p1_q0_r0_s\rangle$ by $\theta$, while everything else is unaffected, and proposed the application of a controlled-Ry gate, as illustrated in Fig.~\ref{fig:circ}(b), up to the phase.\cite{Yordanov20} This implementation requires only 13 CNOT gates, although the connectivity between $pqrs$ is assumed. They have also established that a standard double excitation $e^{\theta \tau_{rs}^{pq}}$ can be similarly implemented (Fig.~\ref{fig:circ}(c)), where the number of CNOT gates is reduced to $2(p-q+r-s)+9$. Using Eq.~(\ref{eq:tau_tilde}) and Fig.~\ref{fig:circ}(b), they proposed qubit-excitation-based (QEB) ADAPT, which was shown to perform significantly better than qubit-ADAPT of Ref.~[\onlinecite{Tang21}] in terms of CNOT gate counts and number of parameters.\cite{Yordanov21}

\begin{figure*}
	\includegraphics[width = 55 em]{circ.pdf}
	\caption{Quantum circuits for (a) $e^{-\frac{i\theta}{8} X_pX_qY_rX_s}$, where the CNOT with a dashed line is a shorthand for the CNOT ladder to convey the parities (refer to the inset), (b) $e^{\theta \hat {\tilde \tau}_{rs}^{pq}}$, and (c) $e^{\theta \hat {\tau}_{rs}^{pq}}$ using the controlled-Ry gate.}\label{fig:circ}
\end{figure*}
\subsection{Symmetry-breaking in ADAPT-VQE} \label{sec:spin-contamination}
In the original proposal by Grimsley et al., it appears to be convenient to choose $\hat{\tau}^{p_\alpha q_\alpha}_{r_\alpha s_\alpha}+\hat{\tau}^{{p}_\beta {q}_\beta}_{{r}_\beta s_\beta}$ and $\hat{\tau}^{p_\alpha q_\beta}_{r_\alpha s_\beta}+\hat{\tau}^{p_\beta q_\alpha}_
{r_\beta s_\alpha}$ as the elements of the operator pool for fermionic ADAPT, because operators in each linear combination are symmetric in terms of spin. However, they are neither generally commutative ($\hat \tau^{p_\alpha q_\beta}_{r_\alpha s_\beta}$ and $\hat{\tau}^{p_\beta q_\alpha}_
{r_\beta s_\alpha}$ do not commute if $p=q$ or $r=s$) nor spin-free. This results in two limitations. First, because the operator selection is based on Eq.~(\ref{eq:Rm}) while the VQE part (implicitly) introduces the Trotter approximation, the gradients between these two steps can be inconsistent, which may cause an occasional numerical instability . In our experience, the same operator can be chosen between consecutive ADAPT cycles, and the convergence to FCI can never be achieved in some cases. Second, the application of a spin-dependent operator basis sacrifices the conservation of spin $s$. Although a more recent study has employed the spin-free doubles generator, that is, the sum of these two combinations,\cite{Tang21} it necessitates an undesired Trotter decomposition, as demonstrated in our previous study.\cite{Tsuchimochi20} Using such spin-incomplete operators provokes a serious issue widely known in electronic structure as spin-contamination. This is especially problematic in strongly correlated systems where the initial HF state significantly deviates from the FCI state and therefore large amplitudes $\theta_k$ are essential. Certainly, a lost spin-symmetry would be eventually recovered by the ADAPT algorithm, provided the target spin state is the ground state and the aforementioned instability is circumvented. However, large spin-symmetry breaking is inevitable in the early steps of the ADAPT procedure, which often leads to unphysical gradients and a poor convergence, as will be demonstrated below. The qubit-ADAPT algorithm\cite{Tang21} is also expected to suffer from the same problem. Furthermore, the candidate operators in qubit-ADAPT break other two important symmetries, $\hat S_z$ and $\hat N$, and the problem becomes even more severe. This fact can be easily confirmed by transforming each of the strings in Eq.\ref{eq:tau_tilde} back to their corresponding fermionic representations; for example, for the case where $p=6, q=3, r=1, s=0$, we have
\begin{align}
\frac{i}{8} Y_6 X_3 X_1 X_0=& \frac{1}{8} \big( a_6^\dag a_3^\dag a_1 a_0+ a_6^\dag a_3^\dag a_1^\dag a_0^\dag \nonumber\\&-2a_6^\dag a_4^\dag a_1^\dag a_0^\dag a_4 a_3+\cdots- h.c.\big).
\end{align}

Hence, it is important to maintain a correct symmetry throughout the ADAPT algorithm, as much as possible, to prevent any unnecessary diversion toward FCI. Accordingly, this study investigates the importance of symmetry-conserving sequence of operators, with a special focus on spin-symmetry.
 One may employ the constraint approach where $\lambda \left(\hat O - o\right)^2$ is added to the energy functional as a penalty to enforce the target state to be an (approximate) eigenfunction of the symmetry operator $\hat O$ with the eigenvalue $o$. Here, we take a different approach and restore the broken symmetry in each step of ADAPT by adopting appropriate post treatments, that is, spin-projection.

\subsection{Spin-projected ADAPT}\label{sec:SP-ADAPT}
In electronic structure theory, spin-projection is a vital tool,\cite{Lowdin55A, Scuseria11, Jimenez12, Tsuchimochi16A, Tsuchimochi18} and has recently exhibited a significant potential in quantum computing in the context of VQE.\cite{Tsuchimochi20} It restores the spin of the broken-symmetry quantum state $|\psi\rangle$, where $\hat S^2 |\psi\rangle \ne s(s+1)|\psi\rangle$, by applying a spin-projection operator $\hat P$, such that
\begin{align}
	\hat S^2 \hat P|\psi\rangle &= s(s+1)\hat P|\psi\rangle,\\
	\hat S_z \hat P|\psi\rangle &= m_s\hat P|\psi\rangle,
\end{align}
where $m_s = (N_{\alpha} - N_\beta)/2$.
Following our previous work, we express $\hat P$ as a linear combination of unitary operators,
\begin{align}
	\hat P \approx \sum_g^{N_g} w_g \hat U_g,\label{eq:P}
\end{align}
where $\hat U_g = e^{-i \alpha_g \hat S_z} e^{-i \beta_g \hat S_y} e^{-i \gamma_g \hat S_z}$ and $w_g$ are weights defined by the Wigner D-matrix. In VQE, we are interested in the expectation value of energy,
\begin{align}
	E = \frac{\langle \psi|\hat P^\dag \hat H \hat P|\psi\rangle}{\langle \psi|\hat P^\dag\hat P|\psi\rangle} \approx \frac{\sum_g w_g\langle \psi|\hat H \hat U_g|\psi\rangle}{\sum_g w_g\langle \psi|\hat U_g|\psi\rangle}\label{eq:E}
\end{align}
where we have used $\hat P^2 = \hat P = \hat P^\dag$ and $\left[\hat H, \hat P\right] = 0$. Note that, using transfer operators $|s;m\rangle \langle s;k|$, which generate a spin $s$ state with $m_s = m$ from $m_s = k$, a more general projector $\sum_{m_s'} c_{m_s'} |s;m_s\rangle \langle s;m_s'|$ introduces more variational freedom, provided $|\psi\rangle$ is not an eigenstate of $\hat S_z$, by diagonalizing the Hamiltonian to determine the coefficients $c_k$\cite{Jimenez12}; however, for simplicity, we only consider $|s;m_s\rangle \langle s;m_s|$, which results in Eqs.~\ref{eq:P} and \ref{eq:E}.
When $|\psi\rangle$ is an eigenstate of $\hat S_z$, Eq.~\ref{eq:E} can be further simplified using the commutativity between $\hat H$ and $\hat S_z$, and one only needs to evaluate the Hamiltonian coupling $\langle \psi|\hat H e^{-i\beta_g \hat S_y}|\psi\rangle$ and the overlap $\langle \psi|e^{-i\beta_g \hat S_y}|\psi\rangle$. These expectation values are estimated, for instance, by the Hadamard test with one ancilla qubit, with a constant number of CNOTs ($4n_{\rm qubit}$). With the Gauss-Legendre quadrature, we take $N_g=2$ in the applications below; hence, the overall measurements will be doubled, while the circuit length is almost unchanged with or without $\hat P$.

The basic idea here is to introduce $\hat P$ to ADAPT-VQE and constrain the operator search in the correct spin space for a better convergence. Accordingly, we simply define the spin-projected (SP) ADAPT by applying $\hat P$ to Eq.~\ref{eq:ADAPT}.
\begin{align}
	|\psi^{\rm SP}_n\rangle = \hat P \; e^{\theta_n \hat{\cal A}_n}e^{\theta_{n-1} \hat{\cal A}_{n-1}} \cdots e^{\theta_1 \hat{\cal A}_1}\Ket{\psi^{\rm HF}}. \label{eq:SP-fADAPT}
\end{align}
The energy gradient evaluated to find the $n+1$th operator is then similar to Eq.~(\ref{eq:Rm}) as follows:
\begin{align}
R_m^{(n)} &= \left.\frac{\partial }{\partial \theta} \frac{\langle \psi_n| e^{-\theta \hat {\cal A}_m} \hat H \hat P \; e^{\theta \hat {\cal A}_m }|\psi_n\rangle}{\langle \psi_n| e^{-\theta \hat {\cal A}_m} \hat P \; e^{\theta \hat {\cal A}_m }|\psi_n\rangle} \right|_{\theta = 0}\nonumber\\
&= \frac{\langle \psi_n| \left[(\hat H - E) \hat P, \hat A_m\right] |\psi_n\rangle}{\langle \psi_n| \hat P |\psi_n\rangle} \label{eq:RmSP}
\end{align}
The commutator involves higher rank operators owing to the presence of $\hat P$, and $R_m^{(n)}$ requires four-body transition RDM between $|\psi_n\rangle$ and $\hat U_g|\psi_n\rangle$. However, in practice, they can be handled by the (Fermionic-)shift rule without explicitly computing RDMs\cite{Schuld19, Kottmann21, Izmaylov21} consistently, both in the operator selection and VQE steps.

The operator candidates used in the original fermionic ADAPT algorithm are pairwise in terms of spin, because if only either of the spin-dependent excitations is employed in ADAPT-VQE, the other would become almost indispensable to restore the spin symmetry. However, when $\hat P$ is introduced, such a pairwise treatment is no longer necessary. For example, consider the single excitation operators $\hat \tau_{q_\alpha}^{p_\alpha}$ and $\hat \tau_{q_\beta}^{p_\beta}$. Although each operator excites (and de-excites) an electron from the $q$-th orbital to the $p$-th orbital, the role of $\hat \tau_{q_\beta}^{p_\beta}$ is simply to complement the lost spin in $\hat \tau_{q_\alpha}^{p_\alpha}$ and {\it vice versa}, which can be handled by $\hat P$. Based on this consideration, it is advisable to employ completely spin-dependent excitation operators as a pool for Fermionic ADAPT ansatz,
\begin{align}
				{\mathscr P}^{\rm spin} = \{\hat{\tau}^{p_\alpha}_{q_\alpha}, \hat{ \tau}^{p_\beta}_{q_\beta}, \hat{\tau}^{p_\alpha q_\alpha}_{r_\alpha s_\alpha}, \hat{\tau}^{{p}_\beta {q}_\beta}_{{r}_\beta s_\beta} ,\hat{\tau}^{p_\alpha q_\beta}_{r_\alpha s_\beta},\hat{\tau}^{p_\beta q_\alpha}_{r_\beta s_\alpha}\} \label{eq:SPPool}
\end{align}
In the presence of $\hat P$, the role of $\hat{\tau}^{p_\alpha q_\beta}_{r_\alpha s_\beta}$ becomes less significant if its spin-flipped operator $\hat{\tau}^{p_\beta q_\alpha}_{r_\beta s_\alpha}$ is already treated in the ansatz. Therefore, some spin-flipped operators that primarily produce effects similar to $\hat P$ would not be chosen in the algorithm, which will potentially minimize the overall circuit depth. Furthermore, the exponential of each operator in the pool ${\mathscr P}^{\rm spin}$ can be treated without the Trotter decomposition; therefore, the energy derivatives estimated in the VQE part and Eq.~(\ref{eq:RmSP}) are equivalent for the last chosen operator. This eliminates the algorithmic problem discussed in Section \ref{sec:spin-contamination}.

Spin-projection is also applicable to qubit-based ADAPT. However, as we have discussed, qubit-ADAPT, which adopts the decomposed Pauli strings as candidate operators, violates the $\hat S_z$ symmetry, and therefore prevents us from using the simplified $e^{-i\beta_g\hat S_y}$ rotation for spin-projection. The full rotation operator $ e^{-i \alpha_g \hat S_z} e^{-i \beta_g \hat S_y} e^{-i \gamma_g \hat S_z}$ allows one to project the correct $m_s$ state but entails both $\alpha$ and $\gamma$ rotation grids, which results in a drastic increase in the number of measurements, by two orders of magnitudes.
In contrast, in QEB-ADAPT, $\hat {\tilde \tau}_{rs}^{pq}$ preserves the $\hat S_z$ and $\hat N$ symmetries. Because the difference between $\hat {\tilde \tau}_{rs}^{pq}$ and $\hat \tau_{rs}^{pq}$ is simply that the former neglects the parities between qubits, we can express it as
\begin{align}
	\hat {\tilde \tau}_{rs}^{pq} = \hat \tau_{rs}^{pq} + \sum_t c_t \hat \tau_{rst}^{pqt} + \sum_{tu} c_{tu} \hat \tau_{rstu}^{pqtu} +\cdots \label{eq:tau_tilde_expand}
\end{align}
where $t,u$ are the qubit indices appearing as a $Z$ string in Eq.~(\ref{eq:tau}), that is, $t,u = [q+1,p-1], [s+1,r-1]$. As a concrete example, it can be verified that
\begin{align}
	\hat {\tilde\tau}_{10}^{63} = \hat {\tau}_{10}^{63} - 2 \hat \tau_{104}^{634} - 2 \hat \tau_{105}^{635} - 4 \hat \tau_{1054}^{6354}.
\end{align}
Hence, spin-projection can be easily combined with QEB-ADAPT, and we call such an algorithm SP-QEB.

Finally, let us briefly consider the tapering-off technique for spin-projection. As it is known, because chemical Hamiltonians possess number and point-group symmetries, one can identify unitaries that transform a Hamiltonian such that it has only Pauli operators that trivially act on certain qubits, which can thus be discarded.\cite{Bravyi17, Setia20} Such unitaries are identified using the $\mathbb{Z}_2$ symmetry; namely, the parities of $\alpha$ and $\beta$ electron numbers for the number symmetry,\cite{Bravyi17} and the sign changes of the underlying wave function by (Abelian) point-group symmetry operations.\cite{Setia20} However, the spin-rotation operator $e^{-i\beta\hat S_y}$ changes the number parity of $\alpha$ and $\beta$ electrons as follows:
\begin{align}
	\hat S_y = \frac{1}{2i}\sum_{p}\left(a_{p_\alpha}^\dag a_{p_\beta} - a_{p_\beta}^\dag a_{p_\alpha}\right)
\end{align}
Therefore, the $\mathbb{Z}_2$-symmetry can be exploited only for the total number operator $\hat N = \hat N_\alpha +\hat N_\beta$, but not for $\hat S_z = \frac{1}{2}(\hat N_\alpha - \hat N_\beta)$  because $\left[\hat S_y, \hat S_z\right] \ne 0$, resulting in the reduction of only one qubit instead of two. This is a necessary cost for spin-projection to exert its advantages; note that the $\hat S^2$ symmetry is not categorized as a $\mathbb{
Z}_2$ symmetry and the tapering-off scheme is not applicable. From a different perspective, spin-projection can explore a ``hidden'' Hilbert space that is not accessible by standard UCC-like ans\"atze by deliberately breaking and restoring the number symmetry of each electron spin. Although one cannot use number parity to discard the other one qubit from the simulation, there is no such restriction for the point-group symmetry. In some of our calculations discussed below, we taper qubits to ease the computational cost. However, we will assume the number of CNOT gates $N_{\rm CNOT}$, which is central to measuring the circuit complexity in this study, remains the same. In addition, we will not discuss how gate operations would change by the transformation in the tapering-off algorithm.

\subsection{First-order molecular properties}\label{sec:properties}
In molecular systems, one is often interested not only in the total energy but also chemical properties. Time-independent molecular properties are defined as the total energy change with respect to a perturbation $x$ introduced to the Hamiltonian. To obtain properties from quantum computing, several studies have focused on evaluating analytical energy derivatives.\cite{OBrien19, Mitarai20} Although higher-order properties such as polarizability require higher-order derivatives, and therefore, are difficult to compute, first-order properties $dE[x]/ dx$ such as force and dipole moment can be easily obtained in VQE. This is because the quantum state is variational with respect to VQE parameters $\{\theta_k\}$ used in the quantum circuit,
\begin{align}
	\frac{\partial E[{\bm\theta}]}{\partial \theta_k} = 0. \label{eq:dEdtheta}
\end{align}
However, in both ADAPT-VQE and SP-ADAPT-VQE, the canonical orbitals of HF commonly employed as a starting point are generally not optimal. The fully parametrized wave function is given by $ \hat P e^{\hat \kappa} |\psi_{\rm ADAPT}[{\bm\theta}]\rangle$ ($\hat P$ is discarded for standard ADAPT-VQE), where
\begin{align}
\hat \kappa = \sum_{p>q} \kappa_{pq} \hat \tau_q^p.
\end{align}
Hence, the total energy
\begin{align}
	E[{\bm\theta}, {\bm\kappa}, x] = \frac{\langle \psi_{\rm ADAPT}[{\bm\theta}] | e^{-\hat \kappa}\hat H[x] \hat P e^{\hat \kappa} |\psi_{\rm ADAPT}[{\bm\theta}]\rangle}{\langle \psi_{\rm ADAPT}[{\bm\theta}] | e^{-\hat \kappa} \hat P e^{\hat \kappa} |\psi_{\rm ADAPT}[{\bm\theta}]\rangle}
\end{align}
is not usually fully stationary with respect to orbital change
\begin{align}
	\left.\frac{\partial E[{\bm\theta},{\bm\kappa},x]}{\partial \kappa_{pq}} \right|_{{\bm\kappa} = {\bf 0}} \ne 0\label{eq:dEdk}
\end{align}
unless it is fully converged to the FCI state. In other words, the 1-particle generalized Brillouin condition is not satisfied; note that the left side of Eq.~(\ref{eq:dEdk}) is equivalent to the gradient $R_{pq}$ defined in Eqs.~(\ref{eq:Rm}) and (\ref{eq:RmSP}) for standard ADAPT and SP-ADAPT ans\"atze, respectively. To address this issue, we can compute the changes in molecular orbitals induced by $x$, by solving the coupled-perturbed HF equation\cite{Gerratt68, Pople79} or the Z-vector equation,\cite{Handy84} similar to previous studies.\cite{OBrien19, Mitarai20} In this study, we adopt an alternative (but mathematically equivalent) formulation commonly used in quantum chemistry,\cite{Helgaker00, Helgaker12} and rederive the necessary equations to obtain first-order molecular properties in SP-ADAPT.

Accordingly, we introduce the following Lagrangian,
\begin{align}
	{\cal L}[{\bm\theta},{\bm\kappa}, {\bf z}, x] = E[{\bm\theta},{\bm\kappa},x] + \sum_{p>q}z_{pq} F_{pq} \label{eq:L}
\end{align}
where $z_{pq}$ and $F_{pq}$ are Lagrange multipliers to be determined and off-diagonal elements of the canonical Fock matrix, respectively. We note that the canonical HF orbitals require $F_{pq} = 0$ for $p\ne q$ at each $x$, and therefore, ${\cal L}$ is stationary with respect to $z_{pq}$. Evidently, for any $x$, ${\cal L}[{\bm\theta},{\bm\kappa}, {\bf z}, x]$ reproduces the same value as $E[{\bm\theta},{\bm\kappa},x]$, and therefore, $d{\cal L}/dx \equiv dE/dx$; however, ${\cal L}$ has an added advantage in that it can be stationary with respect to all variational parameters: ${\bm \theta}$, ${\bm \kappa}$, and {\bf z} (note that $\partial{\cal L}/\partial \theta_k =0$ is automatically satisfied by VQE, see Eq.~(\ref{eq:dEdtheta})). The stationary condition for ${\bm\kappa}$ can be achieved by solving the linear equation to determine ${\bf z}$:
\begin{align}
\frac{	\partial {\cal L}[{\bm\theta},{\bm\kappa}, {\bf z}, x]}{\partial \kappa_{rs}} = R_{rs} + \sum_{p>q} z_{pq} A_{pq,rs} = 0 \label{eq:R+zA}
\end{align}
where
\begin{align}
	A_{pq,rs} = \frac{\partial F_{pq}}{\partial\kappa_{rs}}
\end{align}
is the Fock derivative with respect to orbital change and can be easily computed by a classical computer. We present the explicit form of {\bf A} and the solution to Eq.~(\ref{eq:R+zA}) in Appendix.

Using the chain rule,
\begin{align}
	\frac{dE}{dx} & = \frac{d{\cal L}}{dx} \nonumber\\ &= \frac{\partial {\cal L}}{\partial x} + \sum_{k} \left(\frac{\partial {\cal L}}{\partial \theta_k}\frac{d\theta_k}{dx}\right) + \sum_{p>q} \left(\frac{\partial {\cal L}}{\partial \kappa_{pq}}\frac{d\kappa_{pq}}{dx} + \frac{\partial {\cal L}}{\partial z_{pq}}\frac{dz_{pq}}{dx} \right) \nonumber\\
	&= \frac{\partial {\cal L}}{\partial x}
\end{align}
where the last term  is straightforward to evaluate once ${\bf z}$ is available. It is noteworthy  that this approach is applicable to higher order derivatives in a mathematically clear way, although we will not go into the details.

The Lagrangian ${\cal L}$ can be expressed by the so-called relaxed density matrices ${\bf D}^{\rm relax}$, which incorporate the response correction ascribed to the orbital change by perturbation (Appendix). It can be easily shown that molecular properties as energy derivatives are computed using the relaxed density matrices, instead of the expectation value of the corresponding observable operator, that is, with the regular density matrices $D_{pq} =\langle \psi|a^\dag_p a_q|\psi\rangle$.
The simple expectation value using the latter is generally less accurate because the orbital response is neglected.
To extend the orbital response correction to SP-ADAPT, we can resort to the derivation in Ref.~[\onlinecite{Tsuchimochi17A}].
In Appendix, we summarize the equations and also review the evaluation of nuclear gradients in the case where HF orbitals are a function of nuclear coordinates.

Although this approach can be applied to any VQE method, in several chemistry-inspired ansatzes such as UCCSD, $R_{rs}$ values are zero or negligible because single excitations may be explicitly treated in the variational optimization. However, we note that the frozen-core (or frozen-virtual) approximation is frequently exercised both in classical and quantum computing, where the lower core (higher virtual) orbitals are considered inert and fixed to those of the reference state that is often HF. Because these frozen orbitals are not explicitly optimized within VQE, a response correction is required. This can be easily achieved by simply expanding the range of orbital space $p,q$ in Eq.~(\ref{eq:L}) to include the frozen orbitals. In such cases, matrix elements with frozen orbitals, denoted by $I, J$, are required. Although the expression for $A_{pq,rI}$ remains the same, the gradient $R_{rI}$ needs to be handled in classical computers because $I$ is not mapped to qubits (refer to Appendix for comprehensive derivations).
Certainly, all these problems would be eliminated if orbitals are optimized,\cite{Takeshita20, Mizukami20}; however, orbital-optimized VQE itself requires these quantities and extra VQE simulations.

As mentioned above, for both fermionic and SP-fermionic ADAPTs, we have access to $R_{rs}^{(n)}$ at each $n$th cycle by construction. Therefore, the energy derivative and first-order molecular properties are readily available by processing density matrices in a classical computer. Conversely, for the qubit-based ADAPT, one needs to explicitly compute $R_{rs}$ unless the quantum state is exact. Otherwise, the property would be computed as an expectation value and can suffer from an uncontrollable error.

\section{Illustrative calculations}

\begin{table}
\caption{Summary of ADAPT protocols.}\label{tb:pool}
    \begin{tabular}{ccccccc}
        \hline\hline
        Method & Pool& Symmetry &  \\
         \hline
        Fermionic& Eq.~(\ref{eq:Pool}) & $\hat S_z, \hat N$\\
        Spin-dependent Fermionic& Eq.~(\ref{eq:SPPool}) &  $\hat S_z, \hat N$\\
        Qubit& $\{X_pY_q, X_pX_qX_rY_s\}$ & None\\
        QEB & $\{{\hat {\tilde \tau}_q^p}, {\hat {\tilde \tau}_{rs}^{pq}}\}$ &  $\hat S_z, \hat N$\\
        \hline\hline
        \end{tabular}
\end{table}

\subsection{Computational details}
\color{black}{}
Before discussing our numerical results, we describe the computational details.
All simulations were conducted using Quantum Unified Kernel for Emulation ({\sc Quket}), developed by us, without the effect of noise.\cite{Quket} {\sc Quket} compiles several open-source libraries to generate a Hamiltonian mapped to a qubit basis and perform quantum simulations with parallel computing. Specifically, it employes {\sc PySCF}\cite{pyscf} to generate molecular orbitals and integrals, and {\sc OpenFermion}\cite{openfermion} to perform the Jordan-Wigner transformation. To simulate quantum circuits, we utilize {\sc Qulacs}.\cite{qulacs} The energy minimization in VQE uses the Broyden-Fletcher-Goldfarb-Shannon (BFGS) algorithm implemented in Scipy. The geometry optimization was performed by interfacing {\sc Quket} with {\sc pyberny}.\cite{pyberny}

For Fermionic-ADAPT, we adopt two different pools, Eqs.~(\ref{eq:Pool}) and (\ref{eq:SPPool}), and the controlled-Ry circuit to perform $e^{\theta \hat\tau_{rs}^{pq}}$ and $e^{\theta \hat\tau_q^p}$ and estimate the number of CNOT gates $N_{\rm CNOT}$ in Fig.~\ref{fig:circ}(c). If the same operator is chosen between two consecutive ADAPT steps because of a numerical error or Trotter error, we choose the operator with the next largest gradient. We will not consider the effect of noise in this work, but simply evaluate the expected performance.

We summarize the ADAPT protocols used in this study and the corresponding operator pools in Table~\ref{tb:pool}. For spin-projection, we always use the pool of either the spin-dependent Fermionic excitations or qubit excitations.

\begin{figure}
	\includegraphics[width = 27 em]{org_fidelity_S2.pdf}

	\caption{Fidelity (solid lines) and $\langle \hat S^2\rangle$ (dotted lines) of ADAPT-VQE, plotted against $N_{\rm CNOT}$ for N$_2$ with different bond lengths: (a) 1.098 {\AA}, (b) 1.8 {\AA}, and (c) 2.5 {\AA}. Circles indicate each ADAPT step.}\label{fig:fidelity}
\end{figure}
\subsection{Spontaneous symmetry breaking in ADAPT}\label{sec:symmetry-breaking}
First, we examine the spontaneous symmetry breaking and its consequence in ADAPT for the N$_2$ molecule with the STO-3G basis set. We have considered three different bond lengths, R$_{\rm N-N} = 1.098, 1.8$, and 2.5 {\AA} to represent weakly and strongly correlated cases. The solid lines in Fig.\ref{fig:fidelity} depict the fidelity of the Fermionic ADAPT-VQE state to the exact FCI state $|\langle \Psi_{\rm FCI}|\psi_{\rm ADAPT} \rangle|^2$ as a function of $N_{\rm CNOT}$, using either ${\mathscr P}$ or ${\mathscr P}^{\rm spin}$, labeled as Fermionic and Fermionic (spin), respectively. In all these cases, ADAPT-VQE approaches the exact state as more operators are added, as expected. The figure indicates that there is no significant difference in convergence between ${\mathscr P}$ and ${\mathscr P}^{\rm spin}$, although the former, the original algorithm proposed in Ref.[\onlinecite{Grimsley19}], is more advantageous in that it requires less VQE parameters, as indicated by the circles.
However, the convergence behavior is significantly different depending on the system. In other words, the stronger the correlation of a system, the slower the convergence of ADAPT-VQE to FCI.

This behavior is attributed to the broken spin characterized by $\langle \hat S^2\rangle$, which is also illustrated in Fig. \ref{fig:fidelity} with the dotted lines. In the strongly correlated case of R$_{\rm N-N} = 2.5$ {\AA} (Fig.~\ref{fig:fidelity}(c)), the HF state exhibits a very small overlap with FCI, with a fidelity of 0.2, which quickly increases to 0.8 after the initial few ADAPT steps with approximately 200 CNOTs. However, simultaneously, $\langle \hat S^2\rangle$ (simply the degree of the spin-symmetry breaking) becomes larger than 1, and at this point, the fidelity convergence gets trapped in a plateau.

Such large spin-symmetry breakings occur initially because, in the VQE step of ADAPT, large VQE parameters are essential for describing strong correlations (i.e., highly excited configurations from the HF reference). Consequently, the $\alpha$ and $\beta$ electrons are treated quite differently, especially in the absence of the Trotter decomposition.\cite{Tsuchimochi20} We should mention that the spin-symmetry is safely conserved if ``paired'' double excitations such as $\hat \tau_{q_\alpha q_\beta}^{p_\alpha p_\beta}$ are chosen. In fact, the first few operators selected in ADAPT are of this type, and $\langle\hat S^2\rangle$ remains as zero (see Figs.~\ref{fig:fidelity} and S1). However, paired double excitations do not have a sufficient ability to describe all types of correlation effects, and spin-unpaired excitations are more suitable in several cases.
With $\langle \hat S^2\rangle>1$, the ADAPT state is a mixture of several different spin states, each with a significant weight. Therefore, the derivative approach expressed in Eq.~(\ref{eq:Rm}) to determine the operator candidate is likely inappropriate, because such an operator search is performed in an incorrect symmetry space.  Of course, by further processing ADAPT with more operators and CNOT gates, the lost spin-symmetry starts to be restored ($\sim$700 CNOTs for Fig.~\ref{fig:fidelity}(c)). At this point, the fidelity also quickly increases to one. Similar but less pronounced results can be confirmed for the tests with shorter bond lengths.

Although fermionic ADAPT follows the $\hat S_z$ and $\hat N$ symmetries, qubit-ADAPT breaks all, and as a result, exhibits a significant slow-down of convergence to FCI. This is illustrated by Fig.~\ref{fig:N2_25}, where we have plotted the fidelity and symmetry expectation values, $\langle \hat S^2\rangle, \langle \hat S_z\rangle,$ and $\langle \hat N\rangle$ using the problematic case of $N_2$ at R$_{\rm N-N} = 2.5$ {\AA}. We note that, in this particular case, both algorithms coincidentally result in similar CNOT gate counts to achieve the FCI state. However, the number of VQE parameters for Fermionic-ADAPT to achieve a fidelity of 0.99 is 30, while that for qubit-ADAPT is 128, implying that more measurements are required for qubit-ADAPT. Interestingly, the symmetry-breaking of qubit-ADAPT in $\hat S_z$ and $\hat N$ is rather moderate; however, the error in $\langle\hat S^2\rangle$ is substantially larger than that of Fermionic-ADAPT. Apparently, such a large error in $\langle\hat S^2\rangle$ severely affects the fidelity, and it is evident that the plateaus that emerge in the fidelity and $\langle\hat S^2\rangle$ plots correspond to each other.

\begin{figure}
	\includegraphics[width = 25em]{N2_25A_qubit.pdf}
	\caption{Comparison between fermionic-ADAPT and qubit-ADAPT for N$_2$ at R$_{\rm N-N}=$ 2.5 {\AA}. Circles indicate each ADAPT step.} \label{fig:N2_25}
\end{figure}

\begin{figure*}
\includegraphics[width = 50em]{Error_vs_CNOT.pdf}
\caption{Energy error from FCI (in Hartree) at each CNOT gate count. (a) N$_2$ at equilibrium geometry (R$_{\rm N-N} = 1.098$ {\AA}), (b) stretched N$_2$ (R$_{\rm N-N} = 2.5$ {\AA}), (c) stretched H$_6$ (R$_{\rm H-H}=2$ {\AA}), (d) half-filled 6-site Hubbard chain ($U=8$). The green and red plots indicate standard UCCSD and broken-symmetry UCCSD, respectively.}\label{fig:Error_vs_CNOT}
\end{figure*}

\begin{figure*}
	\includegraphics[width = 50em]{Error_vs_param.pdf}
	\caption{Same as Fig.~\ref{fig:Error_vs_CNOT} but as a function of the number of parameters.}\label{fig:Error_vs_param}
\end{figure*}
\begin{figure}
	\includegraphics[width = 25em]{H6_Nparam_vs_CNOT.pdf}
	\caption{Increase in the number of CNOT gates with the number of parameters for H$_6$.
	}\label{fig:Nparam_vs_CNOT}
\end{figure}

\subsection{Accuracy and gate efficiency in ADAPT and SP-ADAPT}\label{sec:C}
As demonstrated above, symmetry-breaking is associated with strong correlation, and occurs spontaneously in the ADAPT algorithms, thereby slowing down the convergence. Therefore, one would naturally expect that maintaining the correct symmetry can mitigate the problem. In this section, we study the effect of spin-projection using several systems, which include the N$_2$ molecule (with R$_{N-N} = 1.098, 2.5$ {\AA}), linear chain of H$_6$ equally separated by 2 {\AA}, and half-filled one-dimensional periodic Fermi-Hubbard model with six sites and $U=8$. For the molecular systems, we have used the STO-3G basis set and six electrons in six orbitals. Therefore, all the models considered here are composed of 12 qubits.

We choose the initial state to be the HF determinant for the molecular systems; however, for the Hubbard model, we intentionally set the charge localized state where the first three sites are doubly occupied and the remaining three sites are empty to break the spatial symmetry.

Fig.~\ref{fig:Error_vs_CNOT} presents the energy error from FCI and the number of CNOT gates in the quantum circuit for several methods, where the shaded area in orange corresponds to chemical accuracy ($< 1$ mHartree error). We first consider the performance of qubit-ADAPT.\cite{Tang21} For the molecular systems (Figs.~\ref{fig:Error_vs_CNOT}(a-c)), no significant advantage is found over Fermionic ADAPT in terms of CNOT gate counts. This is because the current implementation of Fermionic ADAPT relies on the controlled-Ry circuit of Ref.~[\onlinecite{Yordanov20}], which can eliminate a large number of CNOT gates. In contrast, we find that qubit-ADAPT is disadvantageous in two aspects. First, the energy lowering in each VQE step is considerably small, as plotted in Fig.~\ref{fig:Error_vs_param}, and it requires several VQE parameters (measurements), which also makes it difficult to converge each VQE step because of a highly non-linear parameter space. Second, it breaks the electron number symmetry; the application of qubit-ADAPT to ionic cases is not straightforward, as illustrated by the example of the half-filled Hubbard model, where qubit-ADAPT falls into the {\it true} ground state of $\langle \hat N\rangle =4$ and $\langle \hat S^2\rangle = 2$.

Now let us focus on the results obtained by spin-projection. Evidently, by eliminating the incorrect spin-symmetry components, SP-ADAPT methods (both for Fermionic and QEB) require much less CNOT gates to achieve the same accuracy as those without spin-projection. SP-Fermionic ADAPT is already effective for the weakly-correlated N$_2$ at equilibrium geometry, as illustrated in Fig.~\ref{fig:Error_vs_CNOT}(a). The CNOT gate reduction in this system with spin-projection (by a factor of $0.6\sim0.7$) is almost solely attributed to the fact that spin-complement excitations are handled by the superposition of spin-rotated states; in contrast, with Fermionic ADAPT, one needs to explicitly construct a quantum circuit to perform spin-complement excitations, except for paired doubles. This can be verified by plotting the energy error against the number of VQE parameters (Fig.~\ref{fig:Error_vs_param}). From Fig.~\ref{fig:Error_vs_param}(a), the number of VQE parameters required to achieve the same accuracy is quite similar between Fermionic and SP-Fermionic ADAPT algorithms, which implies the two algorithms yield similar quantum states for each ADAPT iteration.

The efficacy of spin-projection becomes more distinct for strongly correlated systems. For the stretched N$_2$ and H$_6$ (Figs.~\ref{fig:Error_vs_CNOT}(b) and (c)), the number of CNOT gates required to reach the FCI ground state with SP-ADAPT is less than half of that with the corresponding broken-symmetry ADAPT. Having said that, spin-projection may not seem to be advantageous because it employs up to 200 CNOT gates. For instance, QEB-ADAPT exhibits a lower energy than SP-QEB-ADAPT in the stretched N$_2$ case, and the convergence of qubit-ADAPT is initially even faster. However, this is an artifact due to the unphysical spin contamination effect, and the quality of the quantum state is far from satisfactory. A large spin contamination also leads to incorrect properties and state assignment. Furthermore, as seen in Fig.~\ref{fig:N2_25}, again, it often causes the plateau.

The number of required parameters in SP-ADAPT (Fig.~\ref{fig:Error_vs_param}) is less than half of standard ADAPT, as expected. As illustrated in Fig.~\ref{fig:Nparam_vs_CNOT}, the scaling of CNOT gate counts with the number of used parameters is essentially identical between ADAPT and SP-ADAPT when the same operator pool is used. Therefore, a lower number of parameters in the latter simply implies a more gate efficient circuit at each energy accuracy. SP-Fermionic ADAPT converges slightly faster than SP-QEB ADAPT; this may simply be because the former contains more operators in the pool. However, it is interesting to observe that both the SP-Fermionic and SP-QEB algorithms require almost the same number of parameters to represent the exact FCI. This suggests we can further reduce the number of operators to form a pool, which can reduce the measurements for derivative evaluations, while expecting an unchanged fine convergence behavior.

Finally, we also performed UCCSD to evaluate the advantages of ADAPT and SP-ADAPT. In our UCCSD implementation, only the gates with non-zero parameters are constructed to minimize the gate depth. Because the molecular systems we consider have a high spatial symmetry, most UCCSD excitations are symmetry-forbidden, thereby permitting gate-efficient circuits. With such a treatment, UCCSD is already as powerful as ADAPT for N$_2$ at equilibrium and H$_6$ (see the green plots in Figs.~\ref{fig:Error_vs_CNOT}(a) and (c)). For the stretched N$_2$, ADAPT appears to be significantly more efficient and accurate; however, this may well be attributed to the spin-symmetry breaking. For a fair comparison, we also present the broken-symmetry UCCSD\cite{Tsuchimochi20} in red, which provides a fully variational solution at the cost of large spin-contamination. In addition, note that the ``full'' UCCSD, which performs symmetry-forbidden gate operations, that is, even if the amplitudes are zero, can be prepared with 2001 CNOT gates for all cases. For the Hubbard model, UCCSD naturally requires all excitations, and thus, 2001 CNOT gates, as we have initiated the calculation with no spatial-symmetry; UCCSD’s performance is similar to QEB-ADAPT. In conclusion, although ADAPT is imminently more flexible than the fixed ansatz of UCCSD, it ``passes through'' a UCCSD-like state in the process of converging to the FCI state; this is quite logical because UCCSD is usually considered an accurate method that is both chemically and theoretically well established.

\begin{figure}
		\includegraphics[width = 25em]{qeb.pdf}
	\caption{Comparison of the CNOT performance between different operator pools in QEB and SP-QEB.}\label{fig:qeb}
\end{figure}

\begin{figure}
		\includegraphics[width = 25em]{qeb_param.pdf}
	\caption{Same as Fig.~\ref{fig:qeb} but as a function of the number of parameters.}\label{fig:qeb_param}
\end{figure}

\subsection{Comparison between different operator pools}\label{sec:pools}
The previous section suggested that, not only is the number of CNOT gates required to obtain the FCI state constantly smaller for QEB-ADAPT than for Fermionic-ADAPT, but also the number of required variational parameters is very similar with spin-projection. We also performed calculations with spin-dependent Fermionic ADAPT using ${\mathscr P}^{\rm spin}$ defined in Eq.~(\ref{eq:SPPool}) and arrived at the same conclusion for broken-symmetry standard ADAPT algorithms.

The authors of Ref.~[\onlinecite{Tang21}] verified that only $2n-2$ Pauli operators can span the $n$-qubit real Fock space. Although this theorem appears very promising, much more operators are usually required in the ADAPT algorithm because it chooses only one operator at a time; thus, the pool has to be physical in some sense (i.e., it should contain fermionic excitations higher than singles to examine the $k$-particle generalized Brillouin theorem). Nevertheless, the theorem suggests the number of operators in the QEB pool can be further reduced, thereby providing an opportunity to save the computational efforts to choose the next operator. Here, we propose the following three additional qubit-excitation pools and compare their performances:
\begin{itemize}
	\item Scheme 1. $\hat {\tilde \tau}_{q_\alpha}^{p_\alpha}$, $\hat {\tilde \tau}_{r_\beta s_\alpha}^{p_\beta q_\alpha}$ ($\beta$ singles and doubles with the same spins are omitted).
	\item Scheme 2. Same as Scheme 1 but with $p \ge q$.
	\item Scheme 3. Same as Scheme 2 but with $r \ge s$.
\end{itemize}

For the 12-qubit systems, the number of operators in each pool is 855, 555, 400, 190, and 155 for spin-dependent Fermionic, QEB, Scheme 1, Scheme 2, and Scheme 3, respectively. We have plotted the energy convergence of the earlier H$_6$ example against CNOT gates in Fig.~\ref{fig:qeb} and parameter numbers in Fig.~\ref{fig:qeb_param}. From these figures, it can be observed that, at convergence (an energy error of 10$^{-6}$), all QEB schemes perform similarly in terms of both the numbers of CNOT gates and VQE parameters. However, each plot behaves slightly differently. For broken-symmetry ADAPT algorithms (indicated by the lines without points), full QEB provides a slightly better description than Scheme 3 before the convergence. This could be because the entire operator pool in QEB contains more choices than Scheme 3. In particular, the latter does not contain spin-complement operators, which are subsequently chosen in several instances because spin-complement excitations result in similar energy gradients. It is also clear from Fig.~\ref{fig:qeb_param} that spin-dependent Fermionic ADAPT generally exhibits a better energy convergence for a given number of VQE parameters, although it is inefficient in terms of CNOT gate counts. For example, with 100 parameters, the energy obtained by Fermionic ADAPT is one order of magnitude more accurate than those by QEB schemes.

With spin-projection, the performances of different QEB schemes are almost equivalent with a certain CNOT gate count and parameter number. In particular, its accuracy is comparable to SP-Fermionic ADAPT when the number of parameters is fixed, because spin-projection can partially replicate the spin-complement effects. However, SP-Fermionic ADAPT is almost always more accurate than SP-QEB-ADAPT, although slightly. In the particular case of H$_6$, their differences are marginal, and may be rather attributed to the fact that the current approach to determining the next operator is not optimal; operators that are not chosen could identify a lower energy than the chosen one when VQE is performed. However, this behavior is general, as is clearly observed in Fig.~\ref{fig:Error_vs_param}, and can seemingly become more significant for larger systems, as will be demonstrated later. Nevertheless, the improved similarity between the SP-Fermionic and SP-QEB algorithms over standard ADAPT is encouraging for practical applications where the actual operator number (or CNOT gate counts) that can be implemented is assumed to be limited.
\begin{figure}
	\includegraphics[width = 27em]{H2O_F.pdf}
	\caption{Error from FCI in calculated dipole moment of H$_2$O by altering the number of operators in ADAPT.
	(a) Spin-dependent Fermionic ADAPT without orbital response,
	(b) Spin-dependent Fermionic ADAPT with orbital response,
	(c) SP-Fermionic ADAPT without orbital response, and
	(d) SP-Fermionic ADAPT with orbital response.}\label{fig:H2O}
	\end{figure}

\begin{figure*}
	\includegraphics[width = 55em]{O3.pdf}
	\caption{(a) Optimized bond length of O$_3$ as a function of the maximum number of CNOTs used in ADAPT. (b) Same as (a) but for the angle. (c) Potential energy curve of O$_3$ in Hartree, computed with $N_{\rm CNOT}=1000$. (d) Energy error from FCI.}\label{fig:O3}
\end{figure*}

\subsection{Molecular properties}\label{sec:property_results}
In this section, we study the molecular properties computed by several ADAPT-based algorithms using the method discussed in Sec.~\ref{sec:properties}. In a practical application of the ADAPT algorithm, we cannot expect to obtain FCI accuracy, but have to make a tradeoff between accuracy and either the circuit depth or the number of measurements. To benchmark their accuracy, we compare the results obtained with the fixed number of VQE parameters, $N_{\rm param}$, or the number of CNOT gates, $N_{\rm CNOT}$.

\subsubsection{Dipole moment of H$_2$O}
Dipole moment is an important first-order property of a molecule, providing a measure of polarity of  the system. To assess the accuracy of ADAPT in computing this quantity, we have employed the H$_2$O molecule with a fixed angle of $104.5^\circ$ and varied the OH bond length symmetrically. Fig.~\ref{fig:H2O} summarizes the dipole moment error from the FCI value calculated at different ADAPT cycles ($N_{\rm param} = 6, 9,$ and $12$). The results of spin-dependent Fermionic ADAPT without and with the orbital response correction are presented in Figs.~\ref{fig:H2O}(a) and (b), respectively. As can be observed from the figures, the orbital response correction is important for describing the dipole moment correctly at short bond lengths ($\sim$1.5 {\AA}). Note that the correction effect is almost absent in UCCSD because it takes care of single excitations explicitly in VQE by playing the role of orbital relaxation. For the larger bond distance, the orbital response correction sometimes worsens the result; however, this means the ADAPT state is not stationary with respect to orbital rotation and the good dipole moments obtained in the unrelaxed results (Fig.~\ref{fig:H2O}(a)) are fortuitous. Having said that, because we are using the HF canonical orbitals in ADAPT, the orbital Hessian matrix {\bf A} in Eq.~(\ref{eq:R+zA}) can be nearly singular at a certain bond length (around 1.9 {\AA}), where the first-order relaxed density matrix of ADAPT is overcorrected. This singular behavior is circumvented by increasing the number of operators such that the residual $R_{rs}$ is close enough to zero.

When spin-projection is applied, we essentially observe a similar trend (see Figs.~\ref{fig:H2O}(c,d)). Although spin-symmetry breaking in ADAPT does not affect the estimated dipole moment because both the singlet and contaminated triplet states give almost zero dipole moments at stretched bond lengths, SP-Fermionic ADAPT shows promising results, yielding better convergence with respect to $N_{\rm param}$, especially in a more strongly correlated region (Fig.~\ref{fig:H2O}(d)).

\subsubsection{Geometry and potential curve of O$_3$}
It is difficult to investigate the accuracy in estimating geometrical parameters of strongly correlated systems as we are restricted to small molecules or toy models because of a limited computational budget, and most small molecules are only weakly correlated at equilibrium. In fact, previous studies have considered such simple molecules,\cite{OBrien19, Mitarai20, Delgado21} which do not allow one to properly evaluate the potential of quantum computers. Another issue is that, for adaptive algorithms like ADAPT, the potential energy surface is not smooth, and determining the minimum can be challenging. However, it is interesting to benchmark how effective the geometry optimization in ADAPT methods can be. We think it is also pedagogical to compare the results of spin-dependent Fermionic and QEB ADAPT algorithms in predicting potential energy surfaces when the same quantum resource is available.

Accordingly, we choose the ozone molecule with the STO-3G basis as our test case, using an active space comprising 9 orbitals and 12 electrons, resulting in an 18-qubit system. To make a direct comparison between the accuracy of each ADAPT method, we have assumed that the number of CNOT gates that one can handle for each calculation is limited (except for UCCSD, where 2534 CNOT gates were required). Namely, all ADAPT results were obtained with less than a certain $N_{\rm CNOT}$ ranging from 300--1000. In all simulations, we tapered qubits to reduce the computational cost but assumed the total number of CNOT gates remain unchanged.
\begin{table}
\tabcolsep=0.7em
\caption{Non-parallelity error in Hartree.}\label{tb:NPE}
    \begin{tabular}{ccccccc}
        \hline\hline
        Fermionic & QEB & SP-Fermionic& SP-QEB & UCCSD \\
         \hline
        0.014 & 0.066 & 0.009 & 0.021 & 0.012\\
        \hline\hline
        \end{tabular}
\end{table}

Figs.~\ref{fig:O3}(a) and (b) depict the optimized bond length ({\AA}) and angle ($^\circ$), where the FCI result is depicted by the red line (R$_{\rm O-O} = 1.400$ {\AA} and $\angle {\rm OOO} = 109.9^\circ$). Because O$_3$ is a two-determinant system, UCCSD is reasonably accurate, as indicated by the dotted gray line. Thus, symmetry-breaking in ADAPT is not so significant for both the Fermionic and QEB states. Our calculations indicate that $\langle \hat S^2\rangle$ in ADAPT is approximately $0.1$ at most, at the equilibrium geometry estimated by FCI. From the figures, the bond length is more accurately predicted by Fermionic ADAPT than by QEB ADAPT, with or without spin-projection, while the results for the angle exhibit the opposite trend. Without spin-projection, the two methods yielded significantly different results for geometry optimization, especially at a larger $N_{\rm CNOT}$. Increasing $N_{\rm CNOT}$ does not always lead to more accurate results, at least up to $N_{\rm CNOT} = 1000$. However, by applying spin-projection, it appears the geometries predicted by the SP-Fermionic and SP-QEB algorithms are improved with $N_{\rm CNOT}$. They are also similar to each other, implying the two methods also result in similar states.

To further investigate these results, we have plotted in Fig.~\ref{fig:O3}(c) the potential energy curve of symmetric bond dissociation of O$_3$ with a fixed angle of 109.9$^\circ$, using $N_{\rm CNOT} = 1000$. In Fig.~\ref{fig:O3}(d), we have also presented the energy error from FCI in Hartree. In all the methods, to a certain degree, the energy is more accurate at a shorter (weakly correlated) distance and becomes worse as the bond is stretched (strongly correlated). While this result is quite reasonable, we found that the change in QEB energy is rather  significant compared to that in the Fermionic methods, especially without spin-projection. When comparing the Fermionic and QEB results, their energy accuracy is inverted at 1.40 {\AA}. We took the non-parallelity error (NPE) of the potential energy curve, which is defined as the difference between the maximum and minimum errors from FCI throughout the potential curve, and tabulated in Table \ref{tb:NPE}. As this table shows, Fermionic ADAPT provides a more parallel curve to that of FCI than QEB ADAPT. This explains why the latter is less accurate than the former in predicting the bond length, as illustrated in Fig.~\ref{fig:O3}(a).

 At this point, we are uncertain whether this failure of QEB can be attributed to its fundamental deficiency or the possibly inappropriate choice of operators in the ADAPT algorithm. Because the deviation between Fermionic and QEB ADAPT results are significantly mitigated by spin-projection (Figs.~\ref{fig:O3}(c,d)), we suspect symmetry breaking to be behind the different behaviors. Therefore, in Fig.~\ref{fig:O3_S2}, we have summarized the $\langle \hat S^2\rangle$ values obtained from the calculations in Fig.~\ref{fig:O3}(c). In fact, we find that, although Fermionic ADAPT retains the same degree of symmetry breaking throughout all bond distances, QEB is more spin-contaminated when the molecule is stretched. As we have seen several times, symmetry breaking can slow down the convergence in ADAPT; therefore, the degree of spin-contamination is directly related to the energy accuracy with a fixed $N_{\rm CNOT}$. Hence, spin-projection plays an important role in equalizing the two schemes to some extent.

\begin{figure}
		\includegraphics[width = 25em]{O3_S2.pdf}
	\caption{Spin-contamination of Fermionic and QEB ADAPT in O$_3$.}\label{fig:O3_S2}
\end{figure}

\section{Conclusions}\label{sec:conclusions}
One of the unsolved problems in quantum chemistry is a consensus way of treating strong correlations. Quantum computing is expected to provide a solution to this problem because it can, in principle, handle entangled quantum states directly mapped onto qubits. The recently proposed ADAPT algorithm has paved the way to construct an efficient quantum circuit adaptively; however, its convergence is slow for strongly correlated systems because of a significant symmetry breaking.

To address this issue, we introduced spin-projection to ADAPT. We demonstrated that spin-projection can be quite effective when combined with ADAPT, offering shallower circuits with lesser CNOT gates and VQE parameters to achieve the same accuracy at the cost of increased measurements. This is especially the case if a gate-efficient circuit is adopted for fermionic and qubit excitations, whose performances are often comparable. Our calculations also indicated that the pool of QEB ADAPT can be reduced by discarding certain classes of spin-dependent qubit excitations. However, we concluded that it is not worthwhile to apply symmetry-projection to qubit-ADAPT because it would require the restoration of lost $\hat S_z$ (and number) symmetry, which considerably increases the number of measurements in evaluating the energy.

We have also derived the first-order energy derivative in the presence of spin-projection, which enabled the calculation of dipole moment and geometry optimization with SP-ADAPT. It was demonstrated that the orbital response correction is important both for ADAPT and SP-ADAPT in the calculation of dipole moment, because they are far from stationary with respect to orbital changes, unlike UCCSD, which is less sensitive to orbital rotation because of the presence of explicit single excitations. Furthermore, the method was applied to the geometry optimization of O$_3$ to quantify the capabilities of Fermionic and QEB schemes. We found that QEB ADAPT is less stable in achieving a constant accuracy throughout the potential energy surface and thus is less predictable than Fermionic ADAPT in terms of optimized geometry. The deviation between the two schemes was caused by a larger spin-contamination in QEB ADAPT and could be largely mitigated by performing spin-projection. However, the reason behind this unfavorable behavior of QEB ADAPT remains unclear, and further studies will be required to elucidate its cause.

\section*{Acknowledgements}
This work was supported by JST, PRESTO Grant Number JPMJPR2016, Japan and also by JSPS KAKENHI, Grant Numbers JP20K15231 and JP18H03900. We are grateful for the computational resources provided by ECCSE, Kobe University.


\section{Introduction}
\begin{figure*}[t]
    \centering
    \includegraphics[width=\linewidth]{fig1_new3.pdf}
    \caption{\textbf{Experimental setup and illustrations of the cascaded Brillouin laser induced backaction.} \textbf{a,} Pump 1 (counterclockwise, grey) and Pump 2 (clockwise, black) waves are generated from an external cavity diode laser (ECDL) amplified by an Erbium doped fiber amplifier (EDFA). Both the pump 1 and pump 2 frequencies are shifted using acousto-optic modulators (AOM) to create a relative frequency offset. The frequency of Pump 1 is further Pound-Drever-Hall (PDH) locked to the cavity resonance using a phase modulator (PM). Green (blue) arrow refers to Brillouin laser waves SBL1 (SBL2) discussed in panel {\bf b}. Pump 1 has a slightly higher power so that the cascaded SBL is generated in the counter-clockwise direction (red). The beat signal of SBL1 and SBL2 is generated by a photodetector (PD) and its frequency is measured by a frequency counter (FC) and then frequency discriminated using a frequency tracking circuit (F-V circ) for subsequent phase sensitive detection by a lock-in amplifier (LIA), and temperature control using a light emitting diode (LED) and thermal electrical cooler (TEC).
    \textbf{b,} Upper panel (low temperature case, T$_{\rm L}$): the Brillouin shift $\Omega_B ({ \rm T}_{\rm L})$ is smaller than the mode frequency difference ($m\times$FSR, where $m$ the mode number difference and FSR is the free spectral range in angular frequency). Backaction on SBL1 pushes its frequency away from the backaction absorption (purple) maximum, thereby increasing its frequency. This decreases the SBL1-SBL2 beating frequency ($\Delta \omega_{s}/ 2 \pi$), when the cascaded laser power increases. Note: SBL2's frequency is not affected by the mode-pushing effect because the backaction absorption is directional according to the phase-matching condition. Lower panel (high temperature case, { \rm T}$_{\rm H}$): the Brillouin shift $\Omega_B ({\rm T}_{\rm H})$ is larger than the mode frequency difference so that SBL1 has its frequency pushed lower. This increases the SBL1-SBL2 beating frequency ($\Delta \omega_{s}/ 2 \pi$), when the cascaded laser power increases.}
\label{Fig1}
\end{figure*}

\begin{figure*}[t]
    \centering
    \includegraphics[width=\linewidth]{Figure2_2.pdf}
    \caption{\textbf{Measurement of Brillouin backaction.}
    \textbf{a,} The beating frequency of the counter-propagating SBLs is plotted versus pump 1 power at a series of resonator temperatures as indicated. Frequencies at all temperatures closely track one another in the noncascaded regime with a slight power dependence is induced by the Kerr effect. In the cascaded regime, Brillouin backaction resolves the temperature differences with beating frequency showing a distinct dependence upon pump power. The slope of this dependence changes sign at ${\rm T_0}$ (approximately 26.561 $^{\circ}$C).

    \textbf{b,} Measured SBL beating frequency change per unit pump power plotted versus temperature tuning. The temperature at perfect phase matching condition ($T_0$) is indicated.  \textbf{c,} Measured SBL beating frequency change plotted versus temperature. Slope is 41 kHz/K at 6.5 mW pump power.
}

\label{Fig2}
\end{figure*}

The realization of microresonator-based Brillouin lasers \cite{SBS_CaF2,Carmon,lee2012chemically,Loh2016_Brillouin_Rb_Cell,Rakich2017Laser,LiGyro,KYang,SiNx_SBS} has generated interest in their potential application to compact and potentially integrated Brillouin systems\cite{SBS_Review_Integrated_Circuits}. Moreover, high-coherence Brillouin lasers, featuring short-term linewidths below 1 Hz \cite{lee2012chemically,SBS_Ref,Loh2016_Brillouin_Rb_Cell,SiNx_SBS}, have been used for precision measurement and signal generation. This includes microwave synthesis \cite{SBLSynth,li2014electro}, interrogation of atomic clocks \cite{Loh2020Atomic_Clock_SBS}, and rotation measurement \cite{LiGyro,Earth_Gyro_CPSBL,SiNx_SBS}; all of which require excellent temperature stability of the laser mode volume\cite{Loh2019Temp_SBS}. Here, we demonstrate a new method for temperature stabilization of Brillouin lasers based on backaction produced by the Brillouin anti-Stokes process. This process is shown to provide for phase sensitive locking to a temperature set point $T_0$ given by the following condition, which expresses the Brillouin phase matching condition (see Fig. \ref{Fig1}) \cite{lee2012chemically, SBS_Ref}:
\begin{equation}
    \Omega_B(T_0) = m \times {\rm FSR}(T_0)
\label{eqn0}
\end{equation}
where $\Omega_B(T)$ is the Brillouin shift at temperature $T$ (see Supplementary Information) and $m \times{\rm FSR}(T)$ is an integer multiple ($m$) of the resonator free-spectral-range (${\rm FSR}(T)$) at temperature $T$. The actual temperature $T_0$ can be set by micro-fabrication control of $FSR$ at a specified operating wavelength.

The short term frequency stability of a stimulated Brillouin laser (SBL) is set by fundamental noise associated with the thermal occupancy of photons involved in the Brillouin process \cite{SBS_Ref}. In high-Q microresonators, this adds a white noise contribution to frequency noise with an equivalent short time linewidth less than 1 Hz \cite{SBS_Ref,LiGyro,SiNx_SBS}. However, on longer time scales the frequency stability is most often set by temperature variations. Here, thermo-optical and thermo-mechanical effects change the cavity resonant frequency \cite{Matsko2007_WGMR_freq_ref_I, Savchenkov2007_WGMR_freq_ref_II,SBS_gain_BW}, while the temperature dependence of the sound velocity causes drifts in the Brillouin frequency shift \cite{silica_sound}. The latter couples temperature to the lasing frequency through mode pulling \cite{SBS_Ref} and can also induce short term linewidth variation through the Brillouin $\alpha$ parameter \cite{Yuan2020Linewidth_enhancement_SBS}. To compensate temperature drift, measurement of cavity temperature using modes belonging to different polarizations was demonstrated \cite{Strekalov2011Temp_birefringent,Fescenko2012Temp_Dual_Mode} and has been recently employed in fiber optic Brillouin laser systems and silicon-nitride chip-resonator systems to stabilize frequency \cite{Loh2019Temp_SBS, Loh2020Atomic_Clock_SBS,UCSB_thermometry}. These dual-polarization modes feature different frequency tuning rates versus temperature, thereby providing a way to convert change in modal temperature to measurement of a frequency change. In contrast to this method which relies upon an external frequency reference to establish locking, the backaction method described here features an intrinsic reference temperature $T_0$ given by Eq. (\ref{eqn0}). It's sensitivity limit is also determined by the fundamental white frequency noise of the laser as opposed to the integrated laser linewidth.

\begin{figure*}
    \centering
    \includegraphics[width=17cm]{Figure3.pdf}
    \caption{\textbf{Thermal tuning of the backaction with power dithering.}
    \textbf{a,} The open-loop SBL beat frequency (backaction regime of Fig. \ref{Fig2}) is measured versus time at a series of chip temperatures. A weak sawtooth power modulation is applied to Pump 1.
    \textbf{b,} Zoom-ins of the plot in \textbf{a} showing how the polarity and amplitude of the sawtooth beat frequency modulation depends on temperature. \textbf{Inset,} Temperature stabilization of the resonator with servo control activated.}
    \label{Fig3}
\end{figure*}

We consider a Brillouin ring laser geometry shown in Fig. \ref{Fig1}a wherein two pumping waves (dark and light grey arrows) are coupled from a waveguide into clockwise and counter-clockwise directions of the resonator. The frequencies of the pumping waves are close to a resonance, but are not necessarily on resonance. Exact details of this geometry and explanations of the Brillouin process are provided in reference \cite{lai2019observation}. Briefly, each pump wave provides power that is sufficient to excite corresponding Stokes laser waves (SBL1 and SBL2 shown as green and blue arrows) that propagate opposite to their pumping wave and with a lower (Brillouin-shifted) frequency as a result of the phase matching condition. Power transfer from the pumping waves occurs by way of the Stokes scattering process. After laser action occurs, however, a strong anti-Stokes process occurs that is driven by the lasing fields. This creates absorption near the pumps which is proportional to laser power, and this absorptive backaction clamps their circulating powers and hence the laser gain. The exact lasing frequencies can be controlled by tuning of the pumping frequencies, which causes frequency pulling of the laser frequencies such that about 1 MHz of pump frequency tuning induces 10s of kHz of Brillouin laser frequency pulling ($\Delta \omega_s / 2 \pi$). This configuration (non-degenerate, counter-propagating SBLs) is the starting point for implementation of the back-action temperature control method and a schematic of the key spectral features is provided in Fig. \ref{Fig1}b.

Now suppose that the power of SBL1 is increased so that it begins to function as a pumping wave for a new laser wave (the cascaded laser wave shown in red in Fig. \ref{Fig1}b). The power of SBL2 is intentionally kept below threshold so that no cascading occurs. As noted above for the original pumping waves, the onset of laser action, now in the cascaded wave, induces absorptive backaction on its pump, SBL1, that is proportional to the cascaded laser power. It is important to note that on account of the phase matching condition, this backaction acts only SBL1 (i.e., it is directional and does not affect SBL2). The absorptive backaction experienced by SBL1 is shown in purple in Fig. \ref{Fig1}b and compensates SBL1 optical gain (not shown in the figure) provided by the original Pump 1. As a result, SBL1 power is held constant, which is equivalent to the gain of the cascaded laser field (shown in orange in Fig. \ref{Fig1}b) being clamped (even if Pump 1 power is increased). The backaction absorption spectrum has a similar spectral profile to the Brillouin gain and the location of its maximum relative to the frequency of SBL 1 is determined by the phase matching condition for the backaction process. Specifically, maximum nonlinear absorption occurs for the condition of perfect phase matching given by Eq. (\ref{eqn0}).

The use of this backaction for temperature control is now briefly summarized, but a detailed analysis is given in Supplementary Information. Associated with the backaction absorption is a dispersion contribution that pushes the frequency of SBL1 away from the absorption center. At $T=T_0$ the pushing is zero since this temperature corresponds to perfect phase matching for which the SBL 1 frequency is at the absorption maximum. However, cases where $T<T_0$ and $T>T_0$ (upper and lower panels in Fig. \ref{Fig1}b) result in slight phase mismatch. Defining a phase mismatch parameter as $\Delta \omega \equiv m \times{\rm FSR} - \Omega_B(T)$, for $\Delta \omega > 0$ (upper panel) and $\Delta \omega < 0$ (lower panel) the SBL 1 frequency is pushed higher and lower, respectively. Since the frequency of SBL2 is not affected by this process, measurement of the frequency difference of SBL1 and SBL2 gives a direct way to measure $T-T_0$. As an aside, in reaching this conclusion, it is important to note that cavity resonant frequency tuning with respect to temperature is common mode for SBL1 and SBL2 as they share the same cavity longitudinal mode.

To implement this measurement it is convenient to use the power dependence of the frequency pushing. Specifically, a weak modulation of the pump 1 power ($P$) will induce a corresponding frequency modulation of SBL1 through modulation of the backaction dispersion. The resulting pump power ($P$) dependence of the SBL2-SBL1 beat frequency ($\Delta\omega_s$) is given by the following equation (see SI):
\begin{equation}
    \frac{\partial \Delta\omega_s}{\partial P} \approx \frac{4g_0\gamma_{ex}}{\hbar\omega_p\gamma^2\Gamma}\frac{d\Omega_B}{dT} (T-T_0)
\label{eqn1}
\end{equation}
where $\Delta\omega_{s}\equiv\omega_{r}-\omega_{s}$, $\omega_{s}$ ($\omega_{r}$) is the absolute angular frequency of SBL1 (SBL2), $g_0$ is the Brillouin gain, $\omega_p$ is the pump angular frequency, $\gamma_{ex}$ is the external coupling rate, $\gamma$ is the total cavity loss rate, and $\Gamma$ is the Brillouin gain bandwidth. This result shows that frequency discrimination of $\Delta\omega_{s}$ combined with subsequent phase sensitive detection will provide an error signal whose magnitude and sign vary as $T-T_0$.

In the experiment, we used a 36 mm-diameter silica wedge resonator on silicon \cite{lee2012chemically} with 8 $\mu$m thickness and wedge angle of 30 degrees. The resonator is packaged (similar to ref. \cite{Earth_Gyro_CPSBL}) with a thermal electrical cooler (TEC), a light emitting diode, and a thermistor.  The TEC and thermistor are used for coarse temperature control and monitoring. The ultra-high-quality factor of the microcavity and the precisely controlled resonator size enable efficient generation of stimulated Brillouin laser action in the opposite propagation directions (operating wavelengths close to $\lambda\approx1553.3$ nm). The intrinsic Q factor is 300 million and the SBL threshold is 0.9 mW. Details on the optical pumping as well as generation of SBL1, SBL2 and the cascaded laser wave are provided in the Figure \ref{Fig1}a caption.

To setup the phase-sensitive servo temperature control, the beat frequency of SBL1 and SBL2 is dithered by modulating pump 1 (and in turn the cascaded SBL power) using a power modulator (orange dashed line). The frequency of the beat signal is monitored using a frequency counter. A frequency tracking circuit demodulates the dithered signal and phase sensitive measurement is performed using a lock-in amplifier (Stanford Research SR830) to generate the error signal. Temperature control applies the feedback signal to a light emitting diode (fine control) and a thermal electrical cooler (coarse control).

Figure \ref{Fig2} shows the measured SBL1/SBL2 beating frequency ($\Delta \omega_s / 2 \pi$) when sweeping the pump power 1 from below to above cascaded laser operation. Sweeps are performed at a series of temperatures as indicated. As expected, when the pump power is low in the non-cascaded regime, all temperatures provide identical traces. The observed slope on all of these traces is the result of frequency shift provided by the Kerr effect \cite{wang2019petermann}. On the other hand, for higher pump powers in the cascaded regime, the back-action effect is apparent with each temperature showing a distinctly different linear dependence on power. Significantly, the slope of this dependence is observed to change sign as discussed above, corresponding to temperatures above and below $T_0$ in Eq. (\ref{eqn1}). This sign change is essential for implementation of the phase sensitive detection servo control. An experimental plot showing the sign change in slope is provided in Fig. \ref{Fig2}b, where the absolute temperature reference T$_0$=26.561 $^{\circ}$C is also measured. The beating frequency is observed to show a linear power dependence on temperature over the narrow range measured (see Fig. \ref{Fig2}c measured at 6.5 mW pump power). The temperature tuning rate is smaller than the corresponding value provided by the dual-polarization approach \cite{Loh2019Temp_SBS}. Nonetheless, it should be possible to substantially increase this rate in resonators designed for forward-Brillouin scattering, wherein the reduced Brillouin-shift is accompanied by much narrower linewidths (see SI). By fitting these data sets, the experimental back-action strength $\partial^2\Delta\omega_s/(\partial P\partial T)$ is measured to be $2\pi\times27.6$ MHz/(W$\cdot$K), and is consistent with the theory from Eq. (\ref{eqn1}) ($2\pi\times22.5$ MHz/(W$\cdot$K), see SI).

The linear power and temperature dependence of the backaction are further illustrated in Figure \ref{Fig3}a where the beating frequency is measured versus time at a series of temperatures.  To illustrate the change in slope with power at each temperature a weak and slowly-varying saw-tooth power modulation is applied. Zoom-in views of the corresponding saw-tooth modulation in the beat frequency are presented in Fig. \ref{Fig3}b. The change in polarity and amplitude of the backaction-induced modulation are apparent as the temperature is set to values above and below T$_0$.

Finally, long term temperature stabilization of the system is demonstrated by closing the servo control loop. Here, a faster sinusoidal power modulation (200 to 500 Hz) is used to generate a small frequency dither on the SBL beating frequency. It is demodulated by a frequency to voltage conversion circuit and the error signal is generated by phase-sensitive detection with a lock-in amplifier as before. The output of a proportional-integral (PI) servo drives a 1 W white LED (see Figure \ref{Fig1}a) for faster fine-control of temperature. The temperature is also controlled by a TEC that provides slower feedback. The temperature feedback result is shown in the inset of Fig. \ref{Fig3}a. After an initial relaxation oscillation, the long term temperature drift ($>10$ seconds) is stabilized.

With the servo-control loop disconnected, but with the power dither active, the beating frequency exhibits a continuous drift as large as 2.3 kHz in an hour, corresponding to around 0.13 $^\circ$C temperature change per hour. This is apparent in both the measured SBL beat frequency (see Fig. \ref{Fig4} inset) and its Allan Deviation (ADEV) measurement presented in the main panel in Fig. \ref{Fig4}). However, with the servo-control loop connected there is no observable drift in the beating frequency over an hour of measurement (see Fig. \ref{Fig4} inset). Also, over 1500 s averaging time (limited by data size of 1 hour) the Allan Deviation remains around 2 Hz, corresponding to about 0.1 mK temperature variation. In the short term, the drift suppression is believed to be limited by the thermal response of the cavity to the LED. In the future, a faster form of temperature feedback (e.g., an integrated resistive heater placed near the resonator) should further reduce this response time. Feed forward frequency correction could also be employed.

\begin{figure}[t]
    \centering
    \includegraphics[width=8.5cm]{figure4.pdf}
    \caption{\textbf{Allan Deviation of SBL beating frequency $\Delta\omega_s$.}
    Allan deviations of SBL beating frequency with (blue) and without (red) servo control. The error bar gives the standard deviation.
    \textbf{Inset,} SBL beating frequency versus time with (blue) and without (red) servo control.
}\label{Fig4}
\end{figure}

In summary, we have investigated Brillouin backaction and shown that it provides a way to phase sensitively lock an optical resonator to an absolute temperature defined by the phase matching condition. The back-action has been shown to induce both linear power and temperature dependences in a readily measured optical beat frequency. The polarity of the power dependence depends upon operation above or below the phase matching temperature.  This feature and the high stability of the beat frequency were used to servo control the optical mode temperature to 0.1 mK stability levels. While not yet at the temperature stability level of the cross polarization method when applied to microcavities (0.008 mK) \cite{OEwavesTstability}, cross-polarization stabilization has existed for a decade. And we believe the initial results presented here can be substantially improved. An important aspect of the phase sensitive control is that the noise limit is determined by SBL frequency noise at large offset frequencies set by backaction modulation. Here, the rate is 500 Hz, but could be set even higher if necessary to avoid laser frequency noise. Moreover, characterization of the absolute optical frequency stability that is achievable by this method is a possible area of investigation. The range of material systems on which Brillouin laser action has been demonstrated (including integrable platforms) suggests that this method can find wide use.

\vbox{}
\noindent \textbf{\large Acknowledgments} \\
This project was supported by the Defense Advanced Research Projects Agency (DARPA) through SPAWAR (grant no. N66001-16-1-4046), the Air Force Office of Scientific Research (FA9550-18-1-0353) and the Kavli Nanoscience Institute. Yu-Kun Lu would like to thank the Caltech SURF program for financial support.

\vbox{}
\noindent \textbf{\large Author Contributions} \\
Y.H.L. and K.V. conceived the microresonator cascaded Brillouin thermometer. Y.H.L., Y.K.L., Z.Y., H. W. and K.V. constructed the theoretical model. M.G.S. fabricated the ultra-high-Q silica microresonator. Y.H.L., Z.Y., and M.G.S. performed the experiment and analyzed the data. All authors contributed to writing the manuscript. K.V. supervised the project.

\vbox{}
\noindent \textbf{\large Competing interests} \\
The authors declare no competing interests.

\vbox{}
\noindent \textbf{\large Data Availability} \\ The data that support the plots within this paper and other findings of this study are available from the corresponding authors upon reasonable request.