\section{Introduction}

An intriguing discrepancy between the estimated values of the cosmological parameters from Planck CMB measurements and some local observations has
persisted over the last couple of years. While measurements of the Hubble constant from the local Universe is converging at $H_0\sim73~{\rm km/s/Mpc}$ with
progressively narrower errorbars~\cite{2001ApJ...553...47F,2016ApJ...826...56R,0004-637X-855-2-136,2018ApJ...861..126R}, Planck CMB measurements
consistently estimate $H_0=67.8 \pm 0.9~{\rm km/s/Mpc}$
for the case of the concordance $\Lambda$CDM model~\cite{Ade:2015xua,Aghanim:2018eyx}. On the other hand we are also confronted with a discrepancy between growth
of structure measurements from weak lensing surveys such as KiDS~\cite{2017MNRAS.465.1454H}, CFHT~\cite{2017MNRAS.465.2033J}, DES~\cite{2017arXiv170801530D}, HSC~\cite{Hikage:2018qbn} and
that deducted from Planck CMB measurements as reflected in the estimated values of $S_8=\sigma_8 \sqrt{\frac{\Omega_{0m}}{0.3}}$. While most of these constraints on the cosmological
parameters from different surveys are in fact model dependent estimations and, in particular, assume the concordance $\Lambda$CDM model of the Universe, the persistence
of the discrepancy has made it imperative to investigate and reveal the nature of these disagreements. In fact, if these
tensions are real and become more statistically significant with future observations (and not attributable to different systematics), we may be forced to seek
physics beyond the concordance model of cosmology that
effectively require exciting extensions or modifications beyond the current model~\cite{Sahni:2014ooa, Zhao:2017cud, Shafieloo:2018gin}.
In this work we seek a modified form of the primordial spectrum deviating from its conventional minimal power-law form within
the context of the concordance model that satisfies different astrophysical and cosmological observations and also alleviates the tensions in the
different estimations of key cosmological parameters. We show that it is indeed possible to project the discrepancies onto the form of the primordial power spectrum (PPS)
and the reconstructed form of the PPS might potentially hint toward some specific physics in the early Universe. Nevertheless, the shape of the reconstructed form of the
PPS might prove to be a guide toward the identification of some systematics in the pipelines of the data reduction and analysis of Planck CMB data as well. First we describe
the formalism used for our analysis and reconstruction, then we present our results and state our conclusion at the end.

\section{Formalism}~\label{sec:formalism}
The primary goal behind this work is to reconstruct a form of the primordial power spectrum that results in the same CMB observables of the best fit
standard $\Lambda$CDM model while also being consistent with substantially different combinations of the background cosmological parameters
to accommodate some discrepancies observed.
To achieve that we invoke strong priors on the key background parameters such as Hubble constant $H_0$, matter density $\Omega_{\rm m}$ and $\sigma_8$
consistent with low redshift observations and then reconstruct the form of the PPS that results in the same CMB temperature observables corresponding to
the best fit $\Lambda$CDM model. For the purpose of reconstruction, we use the modified Richardson-Lucy (MRL) formalism that has been
introduced in~\cite{HSSWMAP9} and has been subsequently used in~\cite{HSSFF,HSSPlanck} for parameter estimation and  consistency check
besides aiming at finding the features in the primordial power spectrum. For early works on the Richardson-Lucy algorithm and its
variants, see~\cite{Richardson:72,Lucy:1974yx,Shafieloo:2003gf,Shafieloo:2006hs,Shafieloo:2007tk}.

Planck baseline is defined using standard $\Lambda$CDM model with power law form of primordial power spectrum. In addition the
sum of neutrino mass is fixed to be $0.06$eV and the effective number of relativistic species is fixed to be 3.046. We use Planck
baseline best-fit power spectrum to Planck 2015 TT and lowT data and we reconstruct the form of primordial power spectra that matches
the Planck best fit power spectrum while accommodating different and specifically chosen set of background cosmological parameters. Therefore,
we use the smooth best fit angular power spectrum ($\Lambda$CDM model to Planck TT data) as an input data in our algorithm. As the
input power spectrum is always positive, we do not use the binned reconstruction part even when the signal-to-noise ratio falls~\cite{HSSWMAP9}.

In our analysis we assume that the signatures of non-standard physics that may give rise to
the discordance in estimated parameters between Planck CMB and other observations can be substantially captured
by the features in the reconstructed PPS. This assumption is motivated from some of our previous analyses.
In~\cite{Shafieloo:2007tk,HSSFF} we demonstrated that allowing a free form of the PPS increases the degeneracy
in the background cosmological parameters substantially. In~\cite{HSSPlanck}, we showed, when effects of CMB lensing is not
included in the analysis, MRL projects its partial imprints in the PPS as some specific form of features.

Since we are using theoretical best fit and not the data, we do not have the disadvantage of fitting noise. In particular, for
Planck 2013 data combined from all frequency channels, we had demonstrated in~\cite{HSSPlanck} that with higher
iterations, CMB lensing effect can be significantly mimicked by the reconstructed features.

To start with and to calculate the radiative transport kernel we fix $H_0=73.48 {~\rm km/s/Mpc}$~\cite{0004-637X-855-2-136}.
To be consistent with results from weak lensing surveys we consider the matter density to be lower than the derived
values for the case of the standard $\Lambda$CDM model fitting the Planck data. We find that using the Planck TT + lowT best fit
for $\omega_{\rm b}=\Omega_{\rm b}h^2$ and $\omega_{\rm CDM}=\Omega_{\rm CDM}h^2$ provides a value of $\Omega_{\rm m}=0.259$ if
$H_0=73.48 {~\rm km/s/Mpc}$ is used which is suitable for our purpose. While heights the CMB acoustic
peaks and dips are sensitive to $\omega_{\rm b}$ and $\omega_{\rm CDM}$~\cite{Hu:2000ti, Aghamousa:2013jvj, Aghamousa:2014nca},
using the same values for our kernel results in minimal features in the form of the primordial spectrum at large scales.

Once we have the reconstructed form of the PPS from our MRL algorithm, we include scope for some fine tuning (in order to maximize the likelihood)
by varying $A_{\rm scale}$ the overall amplitude of the reconstructed PPS, allowing a lateral shift of the feature positions using
$\Delta\ln k$ and also performing a Gaussian smoothing of the features with a smoothing width of $\Delta_{\rm smooth}$.  Next we
use the reconstructed PPS for cosmological background parameter estimation. While in the case of the standard $\Lambda$CDM model
one assumes a power-law form of the PPS to perform parameter estimation, here we fix the form of the PPS to be what we have reconstructed
and then establish how much cosmological parameters can vary while satisfying Planck temperature observations. In particular we are
interested to see how much $\Omega_{\rm m}$ and $H_0$ can vary around their fixed assumed values used to reconstruct the primordial power spectrum.

For background parameters we use the conventional parameters $\Omega_{\rm b}h^2$, $\Omega_{\rm CDM}h^2$, $\theta$ and $\tau$ that denotes baryon and cold dark matter densities, the ratio of the sound horizon to the angular diameter distance at decoupling and the Thomson scattering optical depth respectively. We use {\tt CAMB}~\cite{Lewis:1999bs} and {\tt CosmoMC}~\cite{Lewis:2002ah} for parameter
estimation. We calculate the angular power spectra at all multipoles without interpolation. We use publicly available Planck 2015 data and likelihood codes~\cite{Aghanim:2015xee,Ade:2015xua}. In particular we
use Planck TT + lowT data for comparison. We do not use the polarization data at this stage to follow a conservative approach.

One important issue here to note is the CMB lensing. In our reconstruction procedure and the likelihood analysis we considered the
lensing effect for different forms of the PPS and the choice of the background parameters.

In MRL formalism the angular power spectrum data is used in unlensed format (${\cal C}_{\ell}^{\rm unlensed}$) since the effect of lensing needs to be subtracted
in order to connect the primordial spectrum to angular power spectrum via our kernel which is a function of the background parameters.
\begin{equation}
{\cal C}_{\ell}^{\rm lensed}={\cal C}_{\ell}^{\rm unlensed}+{\cal C}_{\ell}^{\rm lensing\hspace{1mm}template}
\end{equation}

In our work we aim to have the same CMB temperature angular power spectrum observables of the best fit $\Lambda$CDM model for the case of our reconstruction. Hence to consider the lensing effect in our analysis, we should have the equality relation,

\begin{widetext}
\begin{equation}
{\cal C}_{\ell}^{\rm unlensed}[\rm Planck~best~fit]+{\cal C}_{\ell}^{\rm lensing\hspace{1mm}template}[\rm Planck~best~fit] ={\cal C}_{\ell}^{\rm unlensed}[\rm assumed\hspace{1mm}model]+{\cal C}_{\ell}^{\rm lensing\hspace{1mm}template}[\rm assumed\hspace{1mm}model]\end{equation}
\end{widetext}

Note that in our analysis the ${\cal C}_{\ell}^{\rm unlensed}[\rm assumed\hspace{1mm}model]$ is used for the reconstruction and our assumed model described
described earlier has the combination of $\Omega_{\rm m}=0.259$ and $H_0=73.48 {~\rm km/s/Mpc}$ where $\omega_{\rm b}$, $\omega_{\rm CDM}$ and $\tau$ are the
same values of the best fit $\Lambda$CDM model.

\section{Results}\label{sec:results}
\begin{figure}
\resizebox{250pt}{200pt}{\includegraphics{ps_HST-eps-converted-to.pdf}}
\caption{\label{fig:psk}
The shape of the reconstructed primordial power spectra assuming $H_0=73.48 {~\rm km/s/Mpc}$ and $\Omega_{\rm m}=0.259$.
While for this combination $\omega_{\rm b}$, $\omega_{\rm CDM}$ have the same values as the case of the best fit $\Lambda$CDM model,
these parameters can satisfy local $H_0$ measurement as well as the constraints from weak lensing surveys. Different power spectra
are drawn from the samples that are consistent with the data after lateral shifts and smoothing to the original reconstructed spectrum.
}

\end{figure}
\begin{figure}
\resizebox{240pt}{240pt}{\includegraphics{compare_omm_H0.pdf}}
\caption{\label{fig:omm-h0} Marginalized 68\% and 95\% confidence contours in the $\Omega_{\rm m}$ and $H_0$ plane. Planck 2015 constraints using power law as PPS (in red) and using reconstructed PPS (in blue) are shown. Note that the reconstruction prefers a higher value of $H_0$ that matches with the local Hubble parameter measurement.}
\end{figure}
\begin{figure}
\resizebox{240pt}{240pt}{\includegraphics{compare_omm_s8.pdf}}
\caption{\label{fig:omm-s8} Marginalized 68\% and 95\% confidence contours of
$\Omega_{\rm m}$ and $\sigma_8$ normalization. Planck 2015 constraints using power law PPS (in red) and using reconstructed PPS
(in blue) are shown while KiDS-450 weak lensing constraints are shown in grey. Note that the discrepancy is completely removed in the case of our reconstruction where there is an overlap between the 1$\sigma$ regions.
}
\end{figure}
\begin{figure*}
\resizebox{160pt}{120pt}{\includegraphics{H0-eps-converted-to.pdf}}
\resizebox{160pt}{120pt}{\includegraphics{omegam-eps-converted-to.pdf}}
\resizebox{160pt}{120pt}{\includegraphics{sigma8-eps-converted-to.pdf}}
\caption{\label{fig:1Dandcorrelation} Marginalized probability distribution of $H_0$ [left], $\Omega_{\rm m}$ [centre], $\sigma_8$ [right] are plotted. The Hubble
parameter constraints are shifted to higher values that agrees with the local measurement of $H_0$~\cite{0004-637X-855-2-136}. The mean value of matter density is substantially lower in the case of our reconstruction compared to the baseline case. Note that in both the cases the constraints are tighter in the case of reconstruction as its form has been fixed while there is a flexibility in the tilt of the power-law PPS. The $\sigma_8$ normalization constraints are comparable in both cases as we allow an overall amplitude shift of the PPS which is directly connected to the amplitude of the $\sigma_8$.}
\end{figure*}

Fig.~\ref{fig:psk} shows the reconstructed forms of the PPS assuming $H_0=73.48 {~\rm km/s/Mpc}$ and $\Omega_{\rm m}=0.259$ where we also allow some
lateral shifts or smoothing of the features after the reconstruction to maximize the likelihood. Results obtained overlap with
each other to generate a very narrow band. The reconstructed form of the PPS show some prominent features at different scales. While
we see a suppression of power at large scales, we notice sharp fluctuations in the form of the PPS at wavenumbers larger than 0.02 $\rm Mpc^{-1}$.

We present a comparison of the parameter constraints between the standard $\Lambda$CDM model with power law PPS and the case of our study in Figs.~\ref{fig:omm-h0}
and~\ref{fig:omm-s8}. The comparison of the marginalized probabilities of these three parameters are provided in Fig.~\ref{fig:1Dandcorrelation}.
The power law case is marked with Planck2015. Using the reconstructed PPS leads to high $H_0$ and low matter
density around the assumed initial values (where we fixed the kernel for reconstruction) as expected. Using the reconstructed form of the PPS reduces the size of the confidence ball considerably because of the fine
tuned features and no flexibility in the tilt of the reconstructed PPS.

Constraints on $\sigma_8$ however remains similar, and even slightly larger than the case of the standard model with power-law form of the PPS as we allow the overall
amplitude of the power spectrum to vary and $\sigma_8$ mainly depends on this  variation. We also plot the marginalized constraints on $\Omega_{\rm m}$ and $\sigma_8$
in Fig.~\ref{fig:omm-s8} from the KiDS observation~\cite{Joudaki:2016kym}
in the background. It is evident that the 3$\sigma$ tension can be removed by our specific form of PPS from reconstruction. In fact the marginalized estimated value of
$S_8=0.811\pm0.03$ from our reconstruction fitting Planck CMB data is clearly more consistent with $S_8=0.745\pm0.038$ from KiDS,  $S_8=0.783^{+0.025}_{-0.021}$ from DES and $S_8=0.737\pm0.038$
from CFHTLenS (in comparison with the case of the power law PPS). It is intriguing that our reconstructed form of the PPS combined with the choice of the parameters we assumed, can indeed satisfy Planck CMB data as well
as the local Hubble measurement and also resolves the tension between the Planck measured values of matter density and matter power normalization $\sigma_8$ with that inferred using the lensing data from KiDS-450, DES and other weak lensing surveys.

\section{Discussion}\label{sec:discussion}
In this paper we project the inconsistencies in the estimated values of the key
cosmological parameters assuming the concordance cosmological model from
different surveys, to the form of the primordial spectrum. In other words,
we reconstruct a form of the primordial spectrum that for a
particular set of cosmological parameters (that are consistent with
different low redshift observations), can result identically as the case
of the concordance $\Lambda$CDM model with power-law form of the PPS
fitting Planck CMB temperature data. In the result section we show that
this is possible. The particular form of the reconstructed PPS derived by enforcing
the higher values of $H_0$ from local measurements also satisfies
the lower value of $S_8$ estimated by various weak lensing surveys. While
the main purpose of this paper has been to demonstrate this possibility,
there are few important issues to note:

\begin{enumerate}
 \item The features at high $k$ values are very similar to the features we
 reconstructed previously in~\cite{HSSPlanck} when we did not consider
 CMB lensing (trying to project the effect on the form of the PPS). While
 in our analysis we have considered the lensing effect for each point in
 the parameter space, our results appear to suggest that the CMB lensing
 templates prefer a background cosmology very close to the concordance
 model with its best fit parameters. It is not so straightforward to
 interpret this observation however, it might be worthy of further
 investigation. Our results indirectly suggest that the lensing
 templates which we use to analyze CMB data can have substantial
 effect on the constraints on the form of the PPS and background
 cosmological parameters. This might not be evident when we use a
 power-law form of the PPS, but the issue brought to light more clearly when
 features in the form of the primordial spectrum are allowed.

 \item At face value it appears to be unnatural to generate the complex
 form of the reconstructed PPS within an
 inflationary scenario without extreme fine
tuning. However, we do not provide any conclusive reason to close the
possibility of a physical early Universe explanation. In fact the
reconstructed form of the PPS at large scales might provide a hint
for a specific phenomenological model while at the small scales we
might have to consider other issues such as
the lensing effects indicated in the earlier point.

 \item Using polarization data it should be possible to validate further the
 possibility of the reconstructed form of the PPS. Likewise, using polarization
 data we might be able to look for a more optimized form of the PPS to remove
 tensions from different observations. This is deferred to future study.

 \item A wider exploration of the underlying parameter space of the cosmological
 model would be essential to reveal potential routes to ameliorate the disagreements
 in cosmological parameters inferred. In this work we fixed two of the key cosmological parameters,
 $H_0$ and $\Omega_{\rm m}$ for reconstruction. An interesting possible
 extension of this work where the parameters are held fixed for the MRL
 reconstruction, is then to allow them to vary within a particular range of interest. This is planned for future works.

 \item We considered a combination of key cosmological parameters in order to
 satisfy the local $H_0$ measurement as well as estimations of the $\sigma_8$ and
 $S_8$ from weak lensing surveys. However, these constraints from weak lensing surveys
 are generally based on assumption of the concordance model itself where power-law
 form of the primordial spectrum is assumed. While we have a different form of the
 PPS from our reconstruction, observational constraints of $\sigma_8$ and $S_8$ might become
 slightly different from what they have been reported! In fact a comprehensive analysis requires an iterative approach
 considering all the theoretical effects as well as the effects on observational
 constraints. This is also beyond the scope of this work which requires a more
 generalized pipeline to estimate cosmological parameters.

 \item Though we did not discuss in the paper, our result may also alleviate the tension
 between the cosmological constant dark energy and the Ly-$\alpha$ forest BAO measurement
 at $z=2.33$~\cite{Bautista:2017zgn}. The tension has been reported as a serious problem
 for the cosmological constant and a hint for an evolving dark energy~\cite{Sahni:2014ooa,Zhao:2017cud}.
 By allowing a lower matter density (with respect to the standard model cosmology fitting Planck data),
 the expansion history $h(z)$ in the case of $\Lambda$ dark energy can have substantially lower values
 at redshifts higher than two that may help to make this model more consistent with the Ly-$\alpha$ BAO observations.

\end{enumerate}

\section*{Acknowledgments}
The authors would like to thank Shahab Joudaki for discussions regarding the KiDS results.
The authors would like to acknowledge the use of APC cluster (\href{https://www.apc.univ-
paris7.fr/FACeWiki/pmwiki.php?n=Apc-cluster.Apc-cluster}{https://www.apc.univ-
paris7.fr/FACeWiki/pmwiki.php?n=Apc-cluster.Apc-cluster}). A.S. would like to acknowledge the support
of the National Research Foundation of Korea (NRF-2016R1C1B2016478).


\section{Introduction}

Spinor Bose gases - in particular, ensembles of ultracold Bose atoms with internal spin degrees of freedom - offer a remarkable platform for the study of quantum fluids and phenomena such as quantum phase transitions and superfluidity \cite{Kawaguchi12,StamperKurn13}. One reason for this is the rich variety of collision-induced, spin-mixing dynamics that are possible in spinor Bose-Einstein condensates (BECs), together with the exotic quantum states that may result from these dynamics. Of particular interest in this context, driven in part by potential application to quantum-enhanced metrology, have been highly entangled (e.g., spin squeezed) states \cite{Pu00,Duan02,Mustecaplioglu02,Sau10,Zhang13}, and remarkable experimental progress has been made in this field of research \cite{Lucke11,Gross11,Hamley12,Hoang13,Peise15,Hoang16NatComm,Linnemann16,Kruse16,Hoang16PNAS,Luo17,Zou18}.

One system that has been of interest is a small, tightly confined BEC of spin-1 atoms in the presence of a magnetic field. In the single mode approximation, the interplay between collisions and the quadratic Zeeman shift is modeled by the Hamiltonian
\begin{equation}
\hat{H} = \frac{\Lambda}{N} \mathbf{\hat{S}}^2 + q \hat{N}_0 , \label{hamiltonianeqn}
\end{equation}
where $\mathbf{\hat{S}} = (\hat{S}_x,\hat{S}_y,\hat{S}_z)$ and $\hat{S}_i$ are collective angular momentum operators for the ensemble of $N$ spin-1 atoms, $\hat{N}_0$ is the population operator for the $m=0$ state, and $\Lambda$ and $q$ characterize the interaction strength and quadratic Zeeman shift, respectively.

Despite its apparent simplicity, the Hamiltonian (\ref{hamiltonianeqn}) admits a variety of ground states dependent on the signs and relative magnitudes of $\Lambda$ and $q$. One of those ground states is the macroscopic spin singlet. This is a state of $N$ atoms where the \emph{collective} angular momentum is zero. This state features strong entanglement and is given by the superposition
\begin{align}
\ket{S=0} &= \sum\limits_{j=0}^{N/2} c_{j} \ket{j,N-2j,j},
\end{align}
using the notation $\ket{n_{-1},n_0,n_{+1}}$ for $n_i$ atoms in the magnetic state $m=i$ and with coefficients given by
\begin{align}
c_0 = \frac{1}{\sqrt{N+1}}, &~~~~ c_j = - \sqrt{\frac{N-2j+2}{N-2j+1}} c_{j-1}.
\end{align}
This state features genuine multipartite entanglement of the entire ensemble \cite{Vitagliano11,Vitagliano14}. This entanglement is useful for quantum metrology \cite{UrizarLanz13} and could be of use in fields including quantum memory \cite{Lidar98} and quantum information processing \cite{Cabello02,Bartlett03}.

Preparing the spin singlet of a BEC as the ground state of Hamiltonian (\ref{hamiltonianeqn}) is, in practice, a very challenging prospect, given the extremely small energy scales involved \cite{StamperKurn13}. However, other methods have been proposed to produce the spin-singlet state. One such method involves a sequence of quantum non-demolition measurements using pulses of light to probabilistically prepare a highly entangled macroscopic spin singlet \cite{Behbood13,Toth10,Hauke13}.

In this work, we propose a method to produce the macroscopic spin singlet via interactions mediated by cavity-assisted Raman transitions. The imbalanced Dicke model for spin-1 atoms \cite{Zhiqiang17,Masson17}, or more simply a Tavis-Cummings model \cite{Masson18PhotonCounting}, can probabilistically produce the spin singlet heralded by the {\it absence} of measured photons in the cavity output. This protocol works with an efficiency equal to the overlap between the initial state and the spin singlet, which, however, is $1/(N+1)$. We thus propose the use of quasi-adiabatic sweep techniques - an established method for these systems in spinor BECs - to enhance that initial overlap. We then show that a protocol of a quasi-adiabatic sweep followed by the probabilistic distillation of the state using a Tavis-Cummings model offers a method to produce a singlet state with very high fidelity and a reasonably high efficiency.

\section{Setup and model}

We consider an ensemble of $^{87}$Rb atoms held tightly within an optical cavity. We assume that the ensemble is sufficiently dilute that we can ignore any direct atom-atom interactions, and instead engineer effective interactions via cavity-assisted Raman transitions \cite{Dimer07,Morrison08PRL,Morrison08PRA,Zhiqiang17,Masson17,Zhiqiang18,Masson18PhotonCounting}, as illustrated in Fig.~\ref{leveldiagram}. Here, as demonstrated recently in \cite{Zhiqiang17}, we consider transitions within the complete $F=1$ hyperfine ground state. In the limit that the detunings of the Raman transitions are much larger than the energy separations of the excited state hyperfine levels, this gives an effective open Dicke model for an ensemble of spin-1 particles. In that case, the evolution of the density operator can be described by the master equation \cite{Zhiqiang17,Masson17,Masson18PhotonCounting}
\begin{equation}
\dot{\rho} = -i[\hat{H},\rho] + \kappa \mathcal{D}[\hat{a}] \rho \label{dickemodelmasterequation},
\end{equation}
where $\hat{a}$ is the cavity mode annihilation operator, $\mathcal{D}[\hat{a}]\rho = 2\hat{a}\rho\hat{a}^\dagger - \rho \hat{a}^\dagger \hat{a} - \hat{a}^\dagger\hat{a}\rho$, and
\begin{eqnarray}\label{dickemodelH}
\hat{H} &&= \omega \hat{a}^\dagger \hat{a} + \omega_0 \hat{S}_z \nonumber
\\
&& ~+ \frac{\lambda_-}{\sqrt{2N}} (\hat{a}\hat{S}_+ + \hat{a}^\dagger \hat{S}_-) + \frac{\lambda_+}{\sqrt{2N}} (\hat{a}\hat{S}_- + \hat{a}^\dagger\hat{S}_+).
\end{eqnarray}
The coefficients of the various terms are determined by light shifts and Raman transition rates; in particular,
\begin{align}
\omega &= \omega_c - \frac{\omega_- + \omega_+}{2} + \frac{Ng^2}{3\Delta} , \\
\omega_0 &= \omega_z - \frac{\omega_- - \omega_+}{2} + \frac{\Omega_+^2 - \Omega_-^2}{24\Delta} ,\\
\lambda_\pm &= \frac{\sqrt{N}g\Omega_\pm}{12\Delta}.
\end{align}
Here, $\omega_c$ is the frequency of the cavity mode, $\omega_\pm (\Omega_\pm)$ are the bare frequencies (Rabi frequencies) of the $\sigma_\pm$-polarized laser fields, $\omega_z$ is the linear Zeeman splitting of the $F=1$ levels, $g$ is the single-atom-cavity coupling strength (for the $^{87}$Rb $D_2$ line cycling transition), $\kappa$ is the cavity field decay rate, and $\Delta$ is the detuning of the fields from the atomic resonance. Note that we assume that the atoms all couple to the cavity mode with the same strength $g$, which can in practice be achieved by confining the atoms tightly at cavity mode antinodes in an optical lattice potential (which may be created via another, far-detuned cavity mode at twice the wavelength of the other mode \cite{Zhiqiang17,Zhiqiang18}).

\begin{figure}[t]
\includegraphics[width=0.35\textwidth]{dickemodelspinorlevelschemerb87.png}
\caption{Implementation of an effective spin-1 Dicke model using the $F=1$ ground state of $^{87}$Rb. Interactions are mediated by detuned Raman transitions on the $D_1$ line composed of a cavity mode (red) and $\sigma_-$- (blue) and $\sigma_+$- (green) polarized lasers.\label{leveldiagram}}
\end{figure}

By turning off the $\sigma_+$-polarized laser and choosing $\omega =\omega_0$, the above model reduces to the resonant Tavis-Cummings model for spin-1 atoms. We will consider this model in Section IV in the context of a dissipative, probabilistic scheme for preparing the spin-singlet state.

Alternatively, we can consider the above model in a dispersive limit, such that $\omega\gg\omega_0, \lambda_\pm$, in which case the cavity mode can be adiabatically eliminated. Considering again the case where $\Omega_+=0$, the model can then be reduced to \cite{Masson17,Davis18,Marino18}
\begin{equation}
\dot{\rho} = -i[\hat{H},\rho] + \frac{\Gamma}{N} \mathcal{D}[\hat{S}_-]\rho,\label{equationdickemodelspinor}
\end{equation}
with
\begin{equation}\label{eq:H}
\hat{H} = \omega_0' \hat{S}_z + \frac{\Lambda}{N} (\hat{S}_x^2 + \hat{S}_y^2) ,
\end{equation}
and in terms of the Dicke model parameters given above, we define new parameters
\begin{equation}
\omega_0' = \omega_0 + \frac{\Lambda}{N}, ~~~ \Lambda = - \frac{\omega\lambda^2}{2(\omega^2 + \kappa^2)}, ~~~ \Gamma = -\frac{\kappa}{\omega}\Lambda.
\end{equation}
By manipulating the microscopic parameters, it is possible to (at least approximately) set $\omega_0' = 0$, but given that we consider a system initiated with all atoms in the $m=0$ level (i.e., with $\braket{\hat{S}_z} = 0$), and that $\hat{S}_z$ is conserved by the Hamiltonian evolution, this term does not impact the evolution in the absence of photon detections, which we show is the heralding condition for the production of the singlet. Hence, we see that the cavity-mediated, coherent spin interactions described by (\ref{eq:H}) can emulate the collisional interactions of BECs in the single mode approximation.

Further to this, an artificial quadratic Zeeman shift can, for example, be produced by a weak, auxiliary $\pi$-polarized laser field acting near the $F'=1$ excited manifold. This shift is considered to be time dependent, and could be adjusted either by moving it closer or further from resonance, or by adjusting the power of the weak field. This leaves us with a Hamiltonian of the form
\begin{equation}
\hat{H} = \frac{\Lambda}{N} (\hat{S}_x^2 + \hat{S}_y^2) - q(t) \hat{N}_0.\label{hamiltoniantimedep}
\end{equation}

\section{Entanglement criteria}

Due to the additional degrees of freedom for particles with spin $>1/2$, entanglement in spinor particles can be quantified by a range of different inequalities \cite{Vitagliano11,Vitagliano14}. From \cite{Vitagliano14}, we have that an ensemble of spin-1 particles with no entanglement satisfies
\begin{equation}
(\Delta \hat{S}_x)^2 + (\Delta \hat{S}_y)^2 + (\Delta \hat{S}_z)^2 \geq N.
\end{equation}
Therefore, if an ensemble breaks the less strict bound
\begin{equation}
\braket{\hat{S}_x^2} +\braket{ \hat{S}_y^2} + \braket{\hat{S}_z^2} \geq N \rightarrow \braket{\hat{\mathbf{S}}^2} \geq N,\label{vitaglianocriterion}
\end{equation}
then that ensemble is entangled. This measure thus has two purposes: it shows us how the spin length decreases as well as acting as an entanglement witness. Furthermore, $\braket{\hat{\mathbf{S}}^2}$ gives a bound on the maximum number of atoms that are not entangled \cite{Toth10}.

\section{Dissipative evolution\label{dissipation}}

We consider the master equation (\ref{dickemodelmasterequation}) with $\omega = \omega_0 = 0$. An initial state with all the atoms in the $m=0$ state can be decomposed into a superposition of Dicke states $\ket{S,0}$ with different spin lengths as (for even $N$ \footnote{Our scheme also works for the $F=2$ hyperfine ground state, i.e., an ensemble of spin-2 atoms, where the initial state $\ket{m=0}^{\otimes N}$ has a finite overlap with the spin-singlet state regardless of whether the number of atoms is even or odd. Spin-1 atoms are considered here for simplicity.})
\begin{align}\label{eq:psi0}
\ket{n=0} \otimes \ket{0,N,0} &= \ket{n=0} \otimes \sum\limits_{k=0}^{N/2} d_k \ket{2k,0} ,
\end{align}
where the distribution $\{ d_k\}$ is strongly peaked around $k\simeq\sqrt{N}$ \cite{Masson18PhotonCounting}.
Excitations are produced in the cavity in conjunction with a spin ladder operator, i.e., via terms in the Hamiltonian of the form $\hat{a}^\dagger \hat{S}_\pm$, and that ladder operator does not operate on the spin singlet (i.e., $\hat{S}_\pm \ket{0,0}=0$). This means that entanglement is generated between the cavity mode and the atoms: an empty cavity with the spin singlet, and non-zero photon numbers with the states of non-zero spin length occurring in (\ref{eq:psi0}). Any photon emitted from the cavity must therefore collapse the state into a superposition of spin states not containing the spin singlet, while a null measurement will project the atomic state into the spin singlet. Monitoring the cavity output thus gives some probability of projecting the state into the spin singlet; the probability is simply the overlap of the initial state with the spin singlet, which is $1/(N+1)$ for the state (\ref{eq:psi0}).

If we turn off one of the lasers; in particular, if we set $\lambda_+ = 0$, then an initial atomic state $\ket{2k,0}$ will evolve, subject to Eqs.~(\ref{dickemodelmasterequation}) and (\ref{dickemodelH}) (i.e., the damped Tavis-Cummings model, when $\lambda_+ = 0$), to a steady state $\ket{2k,- 2k}$, with the emission of $2k$ photons from the cavity \cite{Masson18PhotonCounting}. An ideal way in which to study this behavior is to use the method of Monte Carlo wave function simulations (or quantum trajectories) \cite{Carmichael1993}, and in Fig.~\ref{dissipationplot} we present results of this approach applied to the model of (\ref{dickemodelmasterequation}) and (\ref{dickemodelH}), where each ``jump'' in the simulations corresponds to the emission of a photon from the cavity. Fig.~\ref{dissipationplot}(a) shows that, after a few cavity lifetimes, the state is always cleanly projected to a definite spin length. After a time $\kappa t = 10$ with $N=10$ atoms, the overlap with the spin singlet is either (essentially) one or zero, depending on the output photon record (Fig.~\ref{dissipationplot}(b)).

The projection, for the parameters of Fig.~\ref{dissipationplot}, occurs on a timescale on the order of the photon lifetime in the cavity, $1/\kappa$. This is because the production of photons in the cavity due to coupling to the atoms happens on a much faster timescale than cavity loss, and the rate of photon emission follows the rate of the slower process. For this to hold for arbitrary atom number we require that photon production happens for \emph{all} states at a rate faster than $\kappa$. The slowest state to produce photons is $\ket{S=2,0}$, and so we require $\lambda_- \sqrt{3/N} \gg \kappa$. If that is not the case, the rate will instead be governed by the intracavity photon production rate, and so the singlet would be heralded by the absence of photon detections over some timescale more than a few cavity lifetimes.

\begin{figure}[t]
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{taviscummings}
\caption{}
\end{subfigure}

\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{taviscummings_hist}
\caption{}
\end{subfigure}
\caption{Spin length, found by calculating $\braket{\hat{\mathbf{S}}^2} = S(S+1)$ and solving for $S$, (a) over time and (b) at $\kappa t = 10$, in 1000  quantum Monte-Carlo trajectories for $N=10$ atoms under master equation (\ref{dickemodelmasterequation}) with $\lambda_- = 6\kappa$, and $\omega = \omega_0 = \lambda_+ = 0$.\label{dissipationplot}}
\end{figure}

With a more realistic photon detection scheme that has an efficiency $\eta$, the absence of emitted photons from the Tavis-Cummings system would project to a mixed state (a normalizing factor is omitted),
\begin{align}
\rho = &|d_0|^2 \ket{0,0}\bra{0,0} \notag \\&+ \sum\limits_{k=1}^{N/2} (1-\eta)^{2k} |d_k|^2 \ket{2k,-2k} \bra{2k,-2k}.
\end{align}
For reasonably high efficiencies (in particular high enough such that $(1-\eta)^{2}|d_2|^2 \ll |d_0|^2$), this should still be dominated by the spin singlet. Alternatively, the Dicke model ($\lambda_+=\lambda_-$), which produces a continuous stream of photons, could be used. This would allow for very high fidelity singlet production even without perfect detection efficiency. It should be noted that very high effective atom-cavity coupling, $\lambda_\pm$, or very long measurement times would be necessary to project out the lower spin state contributions.

The probability of this procedure working is simply the initial overlap with the singlet state, $1/(N+1)$. For small numbers of atoms this process offers an efficient method to prepare a highly entangled state. However, for 1000 atoms, the maximum efficiency would be 0.1$\%$. Other experimental considerations, such as photon detection efficiency, would further reduce the overall efficiency.

This may be mitigated somewhat by multiple runs with a single atomic ensemble. Using a feedback system conditioned on a photon detection (or some threshold depending on dark count rates), switching back on the repumping system to reinitialize the state into the state $\ket{0,N,0}$ would reintroduce an overlap with the spin singlet. Since the time required for each run is relatively short, this should allow for multiple runs over the course of the lifetime of the ensemble.

\section{Quasi-adiabatic methods in spinor BECs}

If we want to improve the efficiency of our procedure, we need to improve the overlap between the initial state and the spin singlet. To do so, we can make use of the various ground states admitted by the spin-collisional Hamiltonian (\ref{hamiltonianeqn}). In the limit where the quadratic Zeeman shift dominates over the collisional interaction, i.e., $|q| \gg |\Lambda|$, and $q$ is negative, the ground state of the system has all atoms in the $m=0$ state, i.e., $\ket{0,N,0}$, which is of course readily prepared via suitable optical pumping. If, instead, $q$ is positive, then the ground state is degenerate between all states for which there are no atoms in the $m=0$ state. Depending on the spread through the $\hat{S}_z$ states that satisfy this, the degeneracy includes entangled states, such as the twin Fock state $\ket{N/2,0,N/2}$, as well as completely classical states such as $\ket{N,0,0}$ or coherent combinations.

 In the other limit, where instead the collisional interaction dominates over the quadratic Zeeman shift, i.e., $|\Lambda| \gg |q|$, there are again two ground states. With ferromagnetic interactions ($\Lambda < 0$, as for $^{87}$Rb) the ground state is degenerate for all states with a maximum spin length $N$. As above, this can range from a highly entangled Dicke state to a completely classical state, depending on how $\hat{S}_z$ is constrained. For anti-ferromagnetic interactions ($\Lambda > 0$, as for $^{23}$Na), there is only one ground state -- the spin singlet.

The existence of these different ground states allows, in theory, for preparation of entangled states by adiabatic passage. For example, an ensemble can be prepared in the (unentangled) state $\ket{0,N,0}$, which is the ground state for large, negative $q$. An adiabatic sweep from there to $q=0$ produces, depending on the sign of $\Lambda$, either the spin singlet \cite{Zhang13,Hoang16PNAS} or, due to the conservation of $\hat{S}_z$, the Dicke state $\ket{S = N, S_z = 0}$ \cite{Zou18}. Alternatively, an adiabatic sweep through to large, positive $q$ would prepare the twin Fock state $\ket{N/2,0,N/2}$ \cite{Zhang13,Luo17}.

However, a key issue with performing this experimentally is maintaining adiabaticity throughout the sweep. At the phase transition points, the energy gap between the ground  and first excited states becomes extremely small (zero in the limit $N\rightarrow\infty$). Since an adiabatic sweep requires parameters to change on a time scale very slow compared to the inverse of the energy gap, the sweep has to be extremely slow through the transition. Hence, a true adiabatic sweep typically faces severe challenges associated with achievable experimental run times.

Hoang et al. \cite{Hoang16PNAS} used an ensemble of 40000 atoms and a sweep time of 35~s (the minimum time for true adiabaticity). After such a long time, less than 25\% of the original BEC remained trapped. Due to the noise that atom loss induces in the magnetization, no measurable entanglement was left in the ensemble. The alternate approach is to ramp faster than adiabatic, introducing a small amount of energy to the system. This moves the state out of the true ground state, but into states near the ground state that still feature high entanglement. Luo et al. \cite{Luo17} instead ramped in 3~s to high negative $q$ attempting to produce the twin Fock state, finding just $4\%$ of the atoms remained in the $m=0$ state. Such a quasi-adiabatic sweep has also been shown to produce metrologically useful entanglement in a BEC \cite{Zou18}.

\section{Hamiltonian evolution\label{hamiltonian}}

\begin{figure}[b!]
\begin{subfigure}{0.495\textwidth}
\includegraphics[width=\textwidth]{expo_length_q7}
\caption{\label{exposweepa}}
\end{subfigure}

\begin{subfigure}{0.495\textwidth}
\includegraphics[width=\textwidth]{expo_overlap_q7}
\caption{\label{exposweepb}}
\end{subfigure}
\caption{Simulations for $N=40$ atoms. (a) Spin length squared and (b) overlap with the spin-singlet state for a sweep $q(t)/\Lambda = 7.0\mathrm{e}^{-0.08\Lambda t}$ with and without losses. Results with losses, i.e. non-zero $\Gamma$, are plotted as an ensemble average of 1000 trajectories. Inset zooms in on the overlap at later times conditioned on zero jumps.\label{exposweep}}
\end{figure}

We now consider evolution with the engineered master equation (\ref{equationdickemodelspinor}) with the Hamiltonian (\ref{hamiltoniantimedep}). Here, the intention is to emulate results from BECs by quasi-adiabatically sweeping the quadratic Zeeman shift, characterized by $q$, from some large value to essentially zero. To truly maximize the overlap, we should optimize the sweep such that it is as close to adiabatic as possible in the given timeframe \cite{Hoang16PNAS,Luo17,Zou18}. However, in the interest of simplicity, we have instead considered just sweeps that are straightforward in form, but that still greatly enhance the overlap with the spin singlet. We find that an exponential decay, $q(t)/\Lambda = q_0 \mathrm{e}^{-\xi \Lambda t}$, produces a much higher fidelity spin singlet compared to linear-in-time or reciprocal-in-time decays.

In Fig.~\ref{exposweep}, we use an average over an ensemble of quantum trajectories to approximate the master equation solution and show that the sweep greatly increases the overlap with the spin singlet, and, concomitantly, greatly reduces the spin length. We also see that the non-adiabaticity adds energy to the system, resulting in large oscillations of the spin length during the sweep. These oscillations, and, in particular, their phase as they are ``frozen out'' by the quadratic Zeeman shift settling at zero, add considerable noise to the resultant spin length and overlap.

\begin{figure}[b!]
\begin{subfigure}{0.46\textwidth}
\includegraphics[width=\textwidth]{expo_overlap_qvary}
\caption{\label{exposweepc}}
\end{subfigure}

\vspace{4mm}

\begin{subfigure}{0.46\textwidth}
\includegraphics[width=\textwidth]{expo_overlap_xi}
\caption{\label{exposweepe}}
\end{subfigure}
\caption{Simulations for $N=40$ atoms. (a,b) Final overlap with the spin singlet with zero emitted photons and (insets) post-selection efficiency (i.e., no photons emitted) given (a) $q(t)/\Lambda = q_0 \mathrm{e}^{-0.08 \Lambda t}$, or (b) $q(t)/\Lambda = q_0 \mathrm{e}^{-\xi \Lambda t}$ with $\Gamma / \Lambda = 0.001$. Post-selection efficiencies are the percentage of 1000 trajectories that would produce zero photons. \label{exposweep2}}
\end{figure}

Using quantum trajectories, we find that the results can be split into two categories: with jumps and without jumps. With even a single jump, corresponding to the emission of a cavity photon, the overlap with the spin singlet becomes zero. This is because the spin singlet is a dark state to the jump operator, $\hat{S}_-$. Trajectories without a jump exhibit high overlap with the singlet state; in fact, for some sweep parameters, the overlap is higher than without losses. This is because the null measurement back-action increases the overlap with the dark state of the system: the spin singlet.

In line with the work in Section \ref{dissipation}, there is the possibility of post-selecting runs with a null photon output measurement, i.e., no jumps. As the spin length increases, so does the number of photons such a state produces. In particular, due to the spread of spin length, certain states in the initial superposition can release huge numbers of photons, greatly increasing $\braket{\hat{\mathbf{S}}^2}$ for that trajectory. For example, with $\Gamma/\Lambda = 0.001$ in Fig.~\ref{exposweep}, we find that 829 of the 1000 trajectories have no jumps, and so have very high overlap with the spin singlet. Of the other 171 trajectories, there are 139 with $\braket{\hat{\mathbf{S}}^2} \approx 6$, i.e., $S\approx2$, commensurate with the production of one or two photons. All bar two trajectories have $\braket{\hat{\mathbf{S}}^2} < 40$, with the remaining two having a final value $\braket{\hat{\mathbf{S}}^2} \sim 1600$. These states therefore each contribute a significant portion of the ensemble-averaged result, but, even with a realistic single photon detector, could easily be discarded by post-selection, since they emit $\sim 40$ photons to attain such a spin length.

Considering Fig.~\ref{exposweep}, we obtain, with $\Gamma = 0$, a state with all 40 atoms entangled, while with $\Gamma>0$ (i.e., finite cavity decay), we have more than 30 atoms entangled according to the criterion (\ref{vitaglianocriterion})  introduced earlier in the paper. Notably, the entanglement bound is strongly violated even for the master equation result, i.e., with no post-selection, and even with the highest considered cavity-mediated, damping $\Gamma / \Lambda = 0.005$.

Choosing the parameters for the sweep involves striking a balance between a slow ramp, finishing with the quadratic Zeeman shift close to zero, and limiting the cavity-mediated losses. However, clear dependences are not obvious due to the oscillations introduced by the non-adiabatic sweep. Instead, wide regions of sweep parameters offer reasonable results. We can see in Fig.~\ref{exposweepc} that it is difficult to pinpoint an ideal value of the initial quadratic Zeeman shift, $q_0$, with a general trend that $q_0 \gtrsim 4\Lambda$ produces higher overlap. However, higher $q_0$ reduces the probability of a trajectory without photon loss. Fig.~\ref{exposweepe} shows that the rate of the exponential decay, $\xi$, is also difficult to choose so as to optimize the overlap, though, broadly speaking, higher $\xi$ gives a higher post-selection efficiency.

We note that, of course, optimising the quasi-adiabatic sweep in more complex ways should enhance the overlap of post-selected trajectories. However, such optimized sweeps might not necessarily be ideal for the probability of zero photons being produced during the sweep and so that would need to be taken into consideration for the optimisation process.

\section{Combination}

We now show that for large atom numbers, the dissipative scheme can offer reasonable efficiency if it is preceded by the sort of quasi-adiabatic singlet preparation discussed above. The total efficiency of the scheme is now the post-selection efficiency of the sweep (i.e., the probability of no photons being emitted during the sweep) multiplied by the efficiency of the dissipative scheme, which is the overlap at the end of the sweep,
\begin{equation}
p = p_s |\braket{\psi | S=0}|^2 .
\end{equation}

\begin{figure}[t!]
\begin{subfigure}{0.235\textwidth}
\includegraphics[width=\textwidth]{bignumberefficiency}
\caption{}
\end{subfigure}
\begin{subfigure}{0.235\textwidth}
\includegraphics[width=\textwidth]{bignumberoverlap}
\caption{}
\end{subfigure}

\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{bignumberoverall}
\caption{}
\end{subfigure}
\caption{Simulations for $N=1000$ atoms, with a sweep $q(t)/\Lambda = q_0 \mathrm{e}^{-0.05\Lambda t}$, to a final time $\Lambda t = 200$. (a) Overlap with the spin singlet, (b) probability of no photons being emitted during the sweep, and (c) the product of those two, giving the efficiency of a combination of sweep and dissipative distillation.\label{largeatoms}}
\end{figure}

Since we are interested in only the final overlap with the singlet state, we can reduce our basis to those states accessible before a jump has occurred; that is, we consider only states with an exact number of pairs in the $m=\pm1$ states, $\ket{k,N-2k,k}$. Such a reduction in basis allows us to greatly increase the number of atoms we simulate. We model the back-action of the null measurement by adding a non-Hermitian term to the Hamiltonian,
\begin{equation}
\hat{H} = \frac{\Lambda-i\Gamma}{N}(\hat{S}_x^2 + \hat{S}_y^2) + q(t)\hat{N}_0 ,
\end{equation}
and calculate the probability of no photons being produced during the sweep from the jump operator expectation value,
\begin{equation}
p_s = \prod\limits_{i=0}^{t_{\mathrm{max}}/\mathrm{d}t} \left\{ 1-2\Gamma \braket{\hat{S}_x^2 + \hat{S}_y^2}_{t_i} \mathrm{d}t \right\} ,
\end{equation}
where $\braket{\hat{\mathcal{O}}}_{t_i}$ is the expectation of the operator at time $t = i\,\mathrm{d}t$. This means that not only do we vastly reduce the basis size, but we also have only to integrate the Schr\"{o}dinger equation once, rather than running a large number of trajectories.

Fig.~\ref{largeatoms} shows the overlap and post-selection efficiency for $N=1000$ atoms. We can see a smooth relationship between the sweep parameters and the post-selection efficiency. Simply, the sweep should be fast and from a small, non-zero $q_0$. Optimisation of the overlap is more complicated, though we see a wide range of parameters for which that overlap is significantly enhanced. The product of these thus gives a wide region of sweep parameters for which the total efficiency to produce a spin singlet of near perfect fidelity for 1000 atoms would be $\sim 10-20\%$. Even allowing for experimental reductions to this efficiency, such a method would produce the state frequently enough to allow for study and, potentially, use of the many-body entangled state.

For potential parameters, we consider parameters for the Tavis-Cummings model of $\set{\lambda_-,\kappa}/(2\pi) \simeq \set{300,10}$ kHz. Large $\omega/2\pi \simeq 10$ MHz would then produce the spinor collisional model necessary for the quasi-adiabatic sweep. Switching to the dissipative method can then be performed by rapidly changing $\omega \rightarrow 0$. From a microscopic perspective, this involves shifting the frequency of the $\sigma_-$-polarized laser closer to the cavity frequency. In this limit, we have an open Tavis-Cummings model, which would perfectly project the singlet state on the order of 1~ms or less. With these parameters, the scenario of Fig.~\ref{largeatoms} corresponds to a sweep that lasts 3.2~ms, and hence to a time for the total procedure of approximately 4~ms.

\section{Conclusion}

We have proposed a method to prepare the spin-singlet state of an ensemble of integer-spin atoms. Despite great interest in such a state, the production of the singlet is an open problem. Its preparation in spinor BEC experiments faces significant challenges; one potential method is that of adiabatic transformation, but this is made very difficult by the tiny energy gap at the phase transition, which greatly inflates the required timescales, as well as by the need to thoroughly minimize any residual magnetic field.

Our methods bypass these issues by using an alternative approach based upon engineered dynamics and projective measurement in cavity QED.
Using a scheme of cavity-assisted Raman transitions to produce an effective spin-1 Tavis-Cummings model, we make use of the fact that the spin singlet is a unique dark state of the system. This means that the production of the singlet state is heralded by the \emph{absence} of photons detected in the cavity output channel.

Using a slightly modified scheme of cavity-assisted Raman transitions enables us to emulate spinor BEC dynamics and thereby implement a quasi-adiabatic transformation to enhance the overlap of the atomic state with the singlet state prior to implementing the effective Tavis-Cummings model, thus enhancing the probabilistic distillation of the singlet state through monitoring of the dissipative output channel. We believe that this approach offers a realistic possibility of reliably producing the spin-singlet state experimentally, and that it highlights the potential of augmenting spinor interaction models with dissipative channels such as cavity modes.

\acknowledgments
The authors acknowledge the contribution of NeSI high-performance computing facilities to the results of this research. New Zealand's national facilities are provided by the New Zealand eScience Infrastructure and funded jointly by NeSI's collaborator institutions and through the Ministry of Business, Innovation and Employment's Research Infrastructure program.

\bibliographystyle{apsrev4-1}


\section{Introduction}

  Recently, methods based on attention mechanisms have been used in various fields of natural language processing \citep{attention, attention-summarization}. In tasks which deal with \textbf{binary relations} that take two words as arguments, e.g., dependency parsing, a good number of these models have achieved high performance \citep{mlp-attention, bilinear-attention, deepbiaffine}.

  \textbf{Biaffine transformation} is a method to incorporate an attention mechanism into binary relations proposed by \cite{deepbiaffine} (following them, we call this method \textbf{biaffine classifier}). They achieved the state-of-the-art performance among graph-based dependency parsers for the English Penn Treebank. In addition, the state-of-the-art transition-based parser on English Penn Treebank uses a biaffine classifier to evaluate the probability distribution of a word coming into the stack at a step point~\citep{stackpointer}.

  While biaffine transformation has rich expressiveness in modeling binary relations, its number of parameters in the weight matrix (bilinear term) is $O(n^2)$ (where $n$ is the number of dimensions). This redundant parameters can result in high degree of freedom of the model, thus causing overfitting especially when a large number of training samples are not available~\citep{hole}.

  In this paper, we attempt to reduce the redundancy by introducing the assumption of either symmetry or circularity in the weight matrix at a biaffine classifier. With either assumption, we can vectorize the matrix and reduce the space complexity to $O(n)$. Additionally, the time complexity becomes $O(n)$ in the case of symmetry and $O(n\log n)$ in the case of circularity with the fast Fourier transform.
Furthermore, while the expressiveness of the model based on symmetry is restricted\footnote{Denote the score of the bilinear term by $f(a,b)$ when a word pair $a$, $b$ has a dependency relation. In this case, if the bilinear term has symmetry, $f(a,b)=f(b,a)$ . Therefore, this is not appropriate for expression of a directed edge. },  one based on circularity is able to express asymmetry relations

In our experiments, we imposed constraints on the biaffine classifiers of the deep biaffine parser and examined the effect on the accuracy of the model. For our experiments, we used the dependency parsing dataset from the CoNLL 2017 shared task~\citep{conll2017}. We chose four languages which have relatively rich training examples and four languages which have fewer. From our experiments, we found that the method with the circularity assumption outperformed the baseline in most of the languages.

\section{Deep Biaffine Parser}

Models introduced in this paper are based on the \textbf{Deep Biaffine Parser} proposed by \citep{deepbiaffine}. They achieved the state-of-the-art accuracy on the CoNLL 2017 shared task for Universal Dependencies~\citep{conll2017}.

This model receives a sequence of words and POS tags, and calculates the probability of an arc between each pair of words as well as a syntactic function label for each arc. For evaluation of scores, it uses Long Short-Term Memories (LSTMs), Multi-Layer Perceptrons (MLPs) and biaffine classifiers.

In the following sections, first, we explain the biaffine transformation which is the essential part of a biaffine classifier while skipping explanation about LSTM and MLP for the sake of simplicity. Then, we describe the overview of the model. It is worth noting that the structure of the model is different from that in \citep{deepbiaffine} in that it utilizes character level information. This is because we used an updated version of the model that was made for the shared task \citep{dozat2017stanford}.

\subsection{Biaffine Transformation}

For the dependency parsing score functions, we use the biaffine transformation shown below to model binary relations. Here, $\oplus$ stands for concatenation of vectors. The first term on the right side represents relatedness score, and the second term the score of $v_i$ and $v_j$ appearing independently.  $b$ is bias term.

\begin{equation}\label{eq:biaffine}
\begin{split}
g(\mathbf{v}_i, \mathbf{v}_j)
&=
\mathbf{v}_i^{\mathrm{T}}\mathbf{Av}_j + (\mathbf{v}_i \oplus \mathbf{v}_j) ^{\mathrm{T}} \mathbf{b} + b
\end{split}
\end{equation}

\subsection{Structure of the Model}

We show the structure of this model below.

\begin{enumerate}

\item
First, this model takes two sequences: words and POS tags. It uses a unidirectional LSTM to encode each words' character-level information into a vector. It then sums this vector with a separate token-level word embeddings.

 \item
It then concatenates the vectors obtained above with POS embeddings and encodes them via a three layer bidirectional LSTM. $y_i$ denotes a vector made by concatenating the hidden states (not including cell states) from the LSTMs of both directions which corresponds to the $i$ th word $w_i$.

\item
These outputs are then transformed with MLPs. Here, we use distinct MLPs for dependents and heads in the prediction of arcs and labels.
\begin{equation}
\begin{split}
\mathbf{v}_{i}^{arc-head} &= \mathrm{MLP}^{arc-head}(\mathbf{y}_i),\\
\mathbf{v}_j^{arc-dep} &= \mathrm{MLP}^{arc-dep}(\mathbf{y}_j), \\
\mathbf{v}_{i}^{label-head} &= \mathrm{MLP}^{label-head}(\mathbf{y}_i),\\
\mathbf{v}_j^{label-dep} &= \mathrm{MLP}^{label-dep}(\mathbf{y}_j).
\end{split}
\end{equation}

where $\mathbf{v}_{i}^{arc-head}, \mathbf{v}_{j}^{arc-dep} \in \mathbb{R}^{n}$ \\
and $\mathbf{v}_{i}^{label-head}, \mathbf{v}_{j}^{label-dep} \in \mathbb{R}^{m}$ \\

\item
The scores of arcs between each word pair are calculated using a biaffine transformation.
\begin{equation} \label{eq:arc}
s_{i,j}^{(arc)} = {\mathbf{v}_i^{arc-head}}^{\mathrm T} \mathbf W \mathbf{v}_j^{arc-dep} + {\mathbf{v}_i^{arc-head}}^{\mathrm T} \mathbf b^{(arc)}.
\end{equation}

\item
Running Chu-Liu/Edmonds algorithm \citep{mst, edmonds} on the scores calculated in (3), we obtain a tree structure that maximizes the total score.

\item
It evaluates a score $s_{i,j}^{(l)}$  of assigning a label $l$ $(l \in \{1, 2, ..., L\};$ $L$: the number of labels) on the arc between the dependent word $w_i$ and its predicted head word $w_j$. The equation is defined below.

\begin{equation} \label{eq:label}
\begin{split}
s_{i,j}^{(l)}
&= {\mathbf{v}_{i}^{label-head}}^{\mathrm{T}} \mathbf{U}^{[l]} {\mathbf v_j^{label-dep}} \\
&+ {({\mathbf{v}_{i}^{label-head}} \oplus \mathbf v_j^{label-dep})}^{\mathrm{T}} \mathbf{b}^{[l]} \\
&+ b^{[l]}.
\end{split}
\end{equation}

Here, a distinct weight matrix $\mathbf{U}^{[l]}$, weight vector $\mathbf{b}^{[l]}$ and bias $b^{[l]}$ are used for each label. The first term on the right side of the equation (4) represents the score of assigning the label $l$ to the arc between dependent $w_i$ and head $w_j$. The second term expresses the score of the label when the dependent and head are given independently.

\end{enumerate}

In our experiment, $\mat{W}\in\mathbb{R}^{n\times n}$ and $\mat{U}^{[l]}\in\mathbb{R}^{m\times m}$ account for about 17\% of the parameters in the deep biaffine parser. By reducing these parameters, we can expect not only improved memory efficiency but also less overfitting.

\section{Proposed Methods}
\label{sec:proposed}
In this section, we introduce our two proposed methods to reduce the number of parameters of the model. We impose either a symmetry or circularity constraint on the weight matrix $\mathbf{W}$ of (\ref{eq:arc}) and $\mathbf{U}^{[l]}~(\forall l\in\{1,2,\dots,L\})$ of (\ref{eq:label}).

\subsection{Symmetric Matrix Constraint}
This method assumes that the weight matrices of the bilinear terms are \textbf{symmetric matrices} and thus are diagonalizable. As a result, we can transform the bilinear term of the score functions into a ``triple inner product'' of two input vectors and a weight vector.
\subsubsection{Diagonalization of the weight matrix in the bilinear term}
When a matrix $\mathbf{W} \in \mathbb{R}^{n \times n}$ is symmetric, it can be diagonalized by an orthogonal matrix $\mathbf{O} \in \mathbb{R}^{n \times n}$ as below:
\[
\mathbf{W} = \mathbf{O} {\mathop{\text{diag}}(\mathbf{w})} \mathbf{O}^{\mathrm{T}}
\]
where $\mathbf{w} \in \mathbb{R}^{n}$ consists of the eigenvalues of $\mathbf{W}$ and $\mathop{\text{diag}}(\mathbf{w}) \in \mathbb{R}^{n \times n}$ represents the diagonal matrix whose diagonal elements are $\mathbf{w}$. With this property, we can rewrite the bilinear term as follows:
\begin{equation}
\begin{split}
\label{eq:diag}
\mathbf{v}_i^{\mathrm{T}} \mathbf{W} \mathbf{v}_j
&=
\mathbf{v}_i^{\mathrm{T}} \mathbf{O} \mathop{\text{diag}}(\mathbf{w}) \mathbf{O} ^{\mathrm{T}} \mathbf{v}_j  \\ &=
{\mathbf{v}^{\prime}_{i}}^{\mathrm{T}} \mathop{\text{diag}}(\mathbf{w}) \mathbf{v}^{\prime}_{j} \\
&=
\langle\mathbf{v}^{\prime}_{i}, \mathbf{w}, \mathbf{v}^{\prime}_{j}\rangle.
\end{split}
\end{equation}
where, $\mathbf{v}^{\prime}_{i}=\mathbf{O}^{\mathrm{T}} \mathbf{v}_i$ and $\mathbf{v}^{\prime}_{j} = \mathbf{O}^{\mathrm{T}} \mathbf{v}_j$, assuming $\mat{O}$ is learned implicitly.
$\langle\mathbf{v}^{\prime}_{i}, \mathbf{w}, \mathbf{v}^{\prime}_{j}\rangle$ is a ``triple inner product'' of $\mathbf{v}^{\prime}_{i}, \mathbf{w}$ and $\mathbf{v}^{\prime}_{j}$ defined by $\langle \mathbf{a}, \mathbf{b}, \mathbf{c} \rangle = \sum_{k=1}^n a_k b_k c_k$. Consequently, the symmetry constraint on the matrix can reduce the number of weight parameters from $n^2$ to $n$.

\subsubsection{Simultaneous Diagonalization}
\label{sim-diag}
When a set of symmetric matrices forms a commuting family, they can be diagonalized by the same orthogonal matrix~\citep{liu2017analogy}.
So we assume the weight matrices $\mathbf{U}^{[1]}, \mathbf{U}^{[2]}, \dots, \mathbf{U}^{[L]}$ of the scoring functions (\ref{eq:label}) form a commuting family. Namely, we assume:
\[
\mathbf{U}^{[p]}\mathbf{U}^{[q]} = \mathbf{U}^{[q]}\mathbf{U}^{[p]},~\forall p, q \in \{1,2, \dots, L\}.
\]
With this assumption, all $L$ weight matrices can be diagonalized simultaneously. Therefore, the vectors $\mathbf{v}_{i}^{label-head}$ and $\mathbf{v}_{j}^{label-dep}$ can be mapped by the same orthogonal matrix for all score functions.
\subsubsection{Score Functions}
Based on the above, under the assumption that all weight matrices are symmetric, we substitute the bilinear term in the biaffine transformations with a triple inner product. First, the score function for arc is defined as follows:
\begin{equation}
\label{eq:1arc}
\begin{split}
s_{i,j}^{(arc)} &= \langle {\mathbf{v}_i^{arc-head}}, \mathbf{w}, \mathbf{v}_j^{arc-dep} \rangle \\
&+ ({\mathbf v_{i}^{arc-head} \oplus \mathbf v_j^{arc-dep}})^{\mathrm T} \mathbf{b}
.
\end{split}
\end{equation}
where, $\mathbf{w} \in \mathbb{R}^n,~\mathbf{b} \in \mathbb{R}^{2n}$. Unlike in (\ref{eq:arc}), the second term of this function contains the arc-dep vector because we confirm that it improves the performance.

Scoring functions for labels are defined as follows:
\begin{equation} \label{eq:1label}
\begin{split}
s_{i,j}^{(l)} &= \langle {\mathbf v_{i}^{label-head}}, \mathbf{u}^{[l]}, {\mathbf v_j^{label-dep}} \rangle  \\
&+ ({\mathbf v_{i}^{label-head} \oplus \mathbf v_j^{label-dep}})^{\mathrm T} \mathbf {b}^{[l]}
.
\end{split}
\end{equation}
where, $\mathbf{u}^{[l]} \in \mathbb{R}^{m}~,\mathbf{b}^{[l]} \in \mathbb{R}^{2m}$. We eliminate the bias term $b^{[l]}$ in (\ref{eq:label}) because we confirm that it does not affect the performance.

\subsection{Circulant Matrix Constraint}
\cite{hole} used a circulant matrix for the bilinear transformation in knowledge graph completion model~\citep{rescal} to reduce the number of parameters and improve the computational efficiency. Following this method, we assume the weight matrices of the bilinear term in the biaffine transformations are circulant and propose new scoring functions based on that.

\subsubsection{Bilinear Transformation Using a Circulant Matrix}
We define the circulant matrix $C(\mathbf{w}) \in \mathbb{R}^{n \times n}$ for a vector $\mathbf{w} \in \mathbb{R}^{n} $ as follows:
\begin{equation}
C(\mathbf{w})=
\begin{bmatrix}
    w_1       & w_n & \dots & w_3 & w_2 \\
    w_2       & w_1 & w_n &  & w_3 \\
    \vdots    & w_2 & w_1 & \ddots & \vdots \\
    w_{n-1}   &     & \ddots & \ddots & w_{n} \\
    w_{n}     &  w_{n-1} & \dots & w_{2} & w_{1}
\end{bmatrix}
.
\end{equation}
where, $\mathbf{w} ^{\mathrm{T}} = (w_1, \dots, w_n)$. . Then, we replace the bilinear term with one where the weight matrix is a circulant matrix $C(w)$ with $n$ parameters:

\begin{equation}
\label{eq:circ-bilinear}
\mathbf{v}_{i} ^{\mathrm{T}} C(\mathbf{w})  \mathbf{v}_{j}.
\end{equation}

\subsubsection{Score Functions}
We propose score functions that employ (\ref{eq:circ-bilinear}) as the bilinear term in the biaffine transformation. The score function for an arc is then defined as follows:

\begin{equation}
\label{eq:2arc}
\begin{split}
s_{i,j}^{(arc)} &= {\mathbf{v}_i^{arc-head}} C(\mathbf{w}) \mathbf{v}_j^{arc-dep}  \\
&+ ({\mathbf v_{i}^{arc-head} \oplus \mathbf v_j^{arc-dep}})^{\mathrm T} \mathbf{b}
.
\end{split}
\end{equation}
where, $\mathbf{w} \in \mathbb{R}^n,~\mathbf{b} \in \mathbb{R}^{2n}$.

The score functions for labels are defined as follows:
\begin{equation} \label{eq:2label}
\begin{split}
s_{i,j}^{(l)} &= {\mathbf v_{i}^{label-head}}C(\mathbf{u}^{[l]}){\mathbf v_j^{label-dep}} \\
&+ ({\mathbf v_{i}^{label-head} \oplus \mathbf v_j^{label-dep}})^{\mathrm{T}} \mathbf{b}^{[l]}
.
\end{split}
\end{equation}
where, $\mathbf{u}^{[l]} \in \mathbb{R}^{m}~,\mathbf{b}^{[l]} \in \mathbb{R}^{2m}$.

\subsubsection{Efficient Computation Using Fast Fourier Transformation}
In this section, we explain how to compute (\ref{eq:circ-bilinear}) efficiently using a fast Fourier transformation (FFT). We denote an $n$-point discrete Fourier transformation (DFT) matrix as $\mathop{\mathfrak{F}} _n\in \mathbb{C}^{n \times n}$. Then, any circulant matrix $C(\mathbf{w}) \in \mathbb{R}^{n \times n}$ can be diagonalized as follows~\citep{gray2006toeplitz}:

\[
C(\mathbf{w}) = {\mathop{\mathfrak{F}}}_n^{-1}\mathop{\text{diag}}({\mathop{\mathfrak{F}}}_n \mathbf{w}){\mathop{\mathfrak{F}}}_n.
\]

With this property, we can rewrite (\ref{eq:circ-bilinear}) as a triple hermitian inner product ~\citep{liu2017analogy}:
\begin{align}
\label{eq:prop2}
\mathbf{v}_{i} ^{\mathrm{T}} C(\mathbf{w}) \mathbf{v}_{j}
&= \mathbf{v}_{i} {\mathop{\mathfrak{F}}}_n^{-1}\mathop{\text{diag}}({\mathop{\mathfrak{F}}}_n \mathbf{w}){\mathop{\mathfrak{F}}}_n \mathbf{v}_{j} \nonumber\\
&= \frac{1}{n} {\overline{{\mathop{\mathfrak{F}}}_n \mathbf{v}_{i}}}^{\mathrm{T}} \mathop{\text{diag}}({\mathop{\mathfrak{F}}}_n \mathbf{w}){\mathop{\mathfrak{F}}}_n \mathbf{v}_{j} \nonumber\\
&= \langle \mathbf{v}^{\prime}_{i}, \mathbf{w}^{\prime}, \mathbf{v}^{\prime}_{j} \rangle \nonumber\\
&= \Re (\langle \mathbf{v}^{\prime}_{i}, \mathbf{w}^{\prime}, \mathbf{v}^{\prime}_{j} \rangle).
\end{align}
Here, $\mathbf{v}^{\prime}_{i}=\overline{{\mathop{\mathfrak{F}}}_n\mathbf{v}_i}$, $\mathbf{v}^{\prime}_j={{\mathop{\mathfrak{F}}}_n\mathbf{v}_j}$,
${\mathbf{w}^{\prime}}=\frac{1}{n}\mathop{\text{diag}}({\mathop{\mathfrak{F}}}_n\mathbf{w})$, and all of them are $n$ dimensional complex vectors. $\Re(\cdot)$ is the operator which takes the real parts of its argument. With this transformation, we can compute the bilinear transformation with a circulant matrix in $O(n \log n)$ using a FFT.

The DFT of an $n$-dimensional vector $\mathbf{x}$,  ${\mathop{\mathfrak{F}}}_n\mathbf{x}$, is conjugate symmetric if and only if $\mathbf{x}$ is a real vector~\citep{hayashi:17-2}.
In our experiment, we initialize $\mathbf{w}'$ with the DFT of a real vector and update it in complex space.  We update ``frequency'' domain (complex space) vectors using only the operations which have correspondence to ``time domain'' vectors. Thus, as described in \cite{hayashi:17-2}, the conjugate symmetry of vectors are kept while learning because their initial values satisfy it.

\subsubsection{Expressiveness of the Bilinear Transformation Using Circulant Matrices}
In this section, we explain about the expressiveness of circulant matrices in relation to an arbitrary matrix $\mat{W}\in\mathbb{R}^{n\times n}$. \cite{complex-journal} show that for any $\mat{W}$, there exists a normal matrix $\mat{W}'\in\mathbb{C}^{n\times n}$ such that $\mat{W}=\Re(\mat{W}')$. Further, as with a symmetric matrix, a normal matrix can be diagonalized as follows:
\begin{equation*}
\mat{W}=\Re(\mat{W}')=\Re(\mat{O}\mathop{\text{diag}}(\mat{w}')\mat{O}^{*}).
\end{equation*}
Here, $\mat{O}\in\mathbb{C}^{n\times n}$ is a unitary matrix, $\mat{O}^{*}$ is the conjugate transpose of $\mat{O}$ and $\mat{w}'\in\mathbb{C}^{n}$ is the complex vector which consists of the eigenvalues of $\mat{W}'$. The bilinear transformation whose weight matrix $\mat{W}$ is replaced with this, can be transformed into (\ref{eq:prop2}). The unitary matrix $\mat{O}$ is a bijective function, so the input vectors $\mat{v}_i$, $\mat{v}_j$ are learned as their one-to-one correspondent $\mathbf{v}^{\prime}_{i}=\mathbf{O}^{\mathrm{T}} \mathbf{v}_i$, $\mathbf{v}^{\prime}_{j} = \mathbf{O}^{\mathrm{T}} \mathbf{v}_j$, assuming that the unitary matrix $\mathbf{O}$ is learned implicitly.
Note that to simultaneously diagonalize the normal matrices whose real parts are the weight matrices $\mat{U}^{[1]},\mat{U}^{[2]},\dots,\mat{U}^{[L]}$ in the scoring functions for labels, we have to assume that they form a commuting family as with the discussion in \ref{sim-diag}.
\section{Related work}
\label{sect:related_works}
\subsection*{Dependency Parsing}

In recent years, various graph-based parsers with attention mechanisms have been proposed.

\cite{mlp-attention} incorporated the attention mechanism used in machine translation~\citep{attention} into their graph-based parser. Their model receives vectors which are made by concatenating LSTM outputs corresponding to each word and its head candidates. Similarly, \cite{bilinear-attention} proposed a graph-based parser where they substitute the MLP-based classifier in~\citep{mlp-attention} with the bilinear one in their multi-task neural model, although they still use the MLP-based one in prediction of labels. Accordingly,~\cite{deepbiaffine} modified the model by~\cite{mlp-attention} using a biaffine classifier instead of an MLP-based one which enables the model to express not only the probability of a word receiving a particular word as dependent but also the prior probability of a word being a head.

Likewise, in transition-based parsing literature, the state-of-the-art parser on the English Penn TreeBank by~\cite{stackpointer} uses an attention mechanism based on a biaffne classifier which calculates the probability distribution of the next word which comes into the stack at each time step, with LSTM outputs corresponding to each word in the input sentence. The models proposed in this paper can be incorporated into these models.

\subsection*{Parameter Reduction in Neural Networks}
Recently, numerous methods toward parameter reduction of neural networks have been proposed.

As a similar approach to proposed methods, there is a method where a projection matrix is decomposed into smaller matrices by lower-rank approximation~\citep{lu:16}. In addition, \cite{ishihara:18} introduced eigenvector decomposition into neural tensor networks ~\citep{ntn} and analyzed the effects of parameter reduction. Although the paper~\citep{ishihara:18} is similar to the present paper in that their methods address parameter reduction in the bilinear term, our work is different in that we apply it to deep biaffine parser.

There are some methods to reduce parameters in a projection matrix by sharing them. \cite{yu:15} used circulant matrices in the fully connected layers. Our models are different from theirs in that we use circulant matrices in the bilinear terms. ~\cite{chen:2015} perform parameter reduction with a hash kernel and ~\cite{vikas:15} with special matrices like toeplitz matrices. While these can be also used for bilinear terms, the methods based on real diagonalization and circulant matrices are superior to them in computational efficiency.

~\cite{hinton:15} proposed a model called distillation and were able to train a model which was more compact than the original one. However, it needs a lot of time for training since it needs to be trained again for distillation. ~\cite{itay:16} achieved a significant reduction of parameters by the quantization, but the reported accuracy is inferior to the original model. Theoretically, these methods can be combined with the proposed methods.

\section{Experiments}
\label{sec:experiments}

\begin{table*}[t!]
\centering
\begin{tabular}{lcccccc}
\toprule
\bf Treebank & \multicolumn{2}{c}{\bf Baseline} & \multicolumn{2}{c}{\bf Symmetry Matrix} &     \multicolumn{2}{c}{\bf Circulant Matrix} \\
\bf             & \bf UAS          & \bf LAS  & \bf UAS      & \bf LAS  & \bf UAS     & \bf LAS \\\cmidrule(lr){1-7}
UD\_Czech & 93.72 & 91.89 & 93.45(-0.27) & 91.50(-0.39) & \bf 93.87(+0.15) & \bf 92.01(+0.12) \\
UD\_German & 87.57 & 84.27 & 87.21(-0.36) & 83.91(-0.36) & \bf 87.61(+0.04) & \bf 84.39(+0.12) \\
UD\_English & \bf 91.05 & \bf 89.42 & 90.95(-0.1) & 89.22(-0.2) & 91.04(-0.01)& 89.31(-0.11) \\
UD\_Chinese & 87.67 & 85.58 & 87.55(-0.12) & 85.41(-0.17) & \bf 87.96(+0.29) & \bf 85.65(+0.07) \\
UD\_Slovenian-SST & 75.63 & 69.53 & 74.94(-0.69) & 68.58(-0.95) & \bf 75.90(+0.27) & \bf 70.24(+0.71) \\
UD\_Latin & 70.90 & 64.53 & 70.27(-0.63) & 63.18(-1.35) & \bf 72.38(+1.48) & \bf 66.00(+1.47)\\
UD\_French-ParTUT & 91.82 & 89.78 & \bf 92.09(+0.27) & 89.94(+0.16) & 91.94(+0.12) & \bf 90.09(+0.31) \\
UD\_Galician-TreeGal & 80.10 & 74.77 & 80.05(-0.05) & 74.78 (+0.01) & \bf 80.24(+0.14) & \bf 75.68(+0.91) \\

\bottomrule
\end{tabular}
\caption{Main results on CoNLL 2017 dataset.}
\label{results-ud}
\end{table*}

\begin{table*}[t!]
\centering
\begin{tabular}{lcccccc}
\toprule
\bf Reduced Samples & \multicolumn{2}{c}{\bf Baseline} & \multicolumn{2}{c}{\bf Symmetry Matrix} &     \multicolumn{2}{c}{\bf Circulant Matrix} \\
\bf             & \bf UAS          & \bf LAS  & \bf UAS      & \bf LAS  & \bf UAS     & \bf LAS \\\cmidrule(lr){1-7}
0 / 4 & \bf 91.05 & \bf 89.42 & 90.95(-0.1) & 89.22(-0.2) & 91.04(-0.01) & 89.31(-0.11) \\
1 / 4 & \bf 90.32 & \bf 88.57 & 90.05(-0.27) & 88.30(-0.27) & 90.29(-0.03) & 88.49(-0.08) \\
2 / 4 & 88.98 & 87.08 & \bf 89.15(+0.17) & \bf 87.17(+0.09) & 88.72(-0.26) & 86.63(-0.45) \\
3 / 4 & 87.24 & 85.08 & 87.27(+0.03) & 85.06(-0.02) & \bf 87.59(+0.35) & \bf 85.38(+0.3) \\
\bottomrule
\end{tabular}
\caption{Results in UD\_English with fewer training samples.}
\label{results-reduced}
\end{table*}

\begin{table*}[t!]
\centering
\begin{tabular}{lcccccc}
\toprule

\bf Treebank & \bf UAS & \bf LAS \\ \cmidrule(lr){1-3}

UD\_Slovenian-SST & 74.98(-0.65) & 69.07(-0.46) \\
UD\_Latin & 71.96(+1.06) & 65.68(+1.15) \\
UD\_French-ParTUT &  92.17(+0.35) & 90.44(+0.66) \\
UD\_Galician-TreeGal & 79.68(-0.42) & 74.82(+0.05) \\

\bottomrule
\end{tabular}
\caption{Baseline model with reduced dimensions.
The numbers in parentheses are the differences from the baseline model with full dimensions.}
\label{results-halfdim}
\end{table*}

\subsection{Dataset and Implementation}

We compared the models described above on several languages in the CoNLL 2017 shared task for Universal Dependency Parsing dataset. We chose four languages which have relatively abundant training examples: UD\_Chinese, UD\_Czech, UD\_English and UD\_German. And we also selected four languages which have fewer training examples: UD\_French-ParTUT, UD\_Galician-TreeGal, UD\_Latin and UD\_Slovenian-SST.

As a baseline model, we used the dependency parser by Timothy Dozat\footnote{https://github.com/tdozat/Parser-v2} which achieved the highest accuracy on the shared task. The structures of the proposed models are based on the baseline model; we changed only the classifier part.

We only modified two hyper-parameters: we used no pretrained embeddings and initialized word embeddings with a uniform distribution. These settings remain the same throughout all experiments unless otherwise stated.

We use gold word segmentation and gold POS tags while word segmentation and POS tagging are included in the shared task. We excluded these two tasks because the objective of this research is to show the effects of the proposed methods on biaffine classifiers which are not used for those tasks.

\subsection{Results}
We show the results of the baseline and two proposed methods in Table~\ref{results-ud}. The method based on circulant matrices outperformed the others on almost all languages except for English where the baseline model achieved the best accuracy and French-ParTUT where the method based on symmetric matrix did so in UAS. Interestingly, the method based on symmetric matrices underperformed the baseline on most languages. This might be because of the restricted expressiveness of a symmetric weight matrix in comparison to a circulant one especially in that the former is not appropriate for expressing asymmetric relations.

\begin{table}[t!]
\centering
\begin{tabular}{lrrr}
\toprule
& \bf \shortstack{Number of \\ Parameters} & \bf Percentage\\ \cmidrule(lr){1-3}
Character LSTM & 241200 & 7.64\% \\
Bidirectional LSTM & 1927200 & 61.03\% \\
Arc MLP & 320800 & 10.16\% \\
Label MLP & 80200 & 2.54\% \\
Arc Classifier & 160400 & 5.08\% \\
Label Classifier & 377437 & 11.95\% \\
Others & 50400 & 1.60\% \\ \hline
TOTAL & 3157637 & 100\% \\
\bottomrule
\end{tabular}
\caption{\label{results} Percentage of baseline model parameters accounted for by each component.}
\label{table:ratio}
\end{table}

\begin{table*}[t!]
\centering
\begin{tabular}{lrrr}
\toprule
\bf  & \bf Baseline &  \bf Symmetry Matrix & \bf Circulant Matrix \\ \cmidrule(lr){1-4}
Arc Classifier & 160400& 1200& 1600\\
Label Classifier & 377437& 11100& 14800\\
Sum with shared parts & 3157637 & 2632100 & 2636200 \\ \cmidrule(lr){1-4}
Difference from the baseline & 0.0\%& \bf -16.64\%& \bf -16.51\%\\
\bottomrule
\end{tabular}
\caption{\label{results} Comparison of parameter sizes.}
 \label{table:comparison}
\end{table*}

\section{Analysis}

In this section, we examine the robustness to overfitting of the proposed methods.

\subsection{Relaxation of Overfitting}

First, to further examine how numbers of parameters affect the models, we conducted experiments on UD\_English treebank reducing the number of training examples by a quarter at a time.

The results of this experiment are shown in Table \ref{results-reduced}. Our methods performed better on smaller datasets where the number of training examples are less than or equal to a half of the original number of examples. This indicates that our methods not only became robust to overfitting through parameter reduction but also achieved high generalizability.

Second, we ran the baseline model with classifier dimensions reduced from 400 to 200 for the arc classifier and 100 to 50 for the label classifier and compared it with the proposed methods on languages from small treebanks. As shown in Table \ref{results-halfdim}, simply reducing the number of dimensions hindered the accuracies in some languages while the method based on circulant matrix consistently outperformed the baseline model with the original number of dimensions in all of these languages from small treebanks as shown in Table \ref{results-halfdim}. This result indicates the effectiveness of the proposed method on smaller datasets.

\subsection{Parameter Reduction}
Table \ref{table:ratio} shows the proportion of the total parameters which each part of the baseline model account for. While LSTMs occupy the largest portion of the model, the second largest part is the classifiers which account for about 17\% of the parameters. Table \ref{table:comparison} indicates that the proposed methods were able to reduce the number of parameters by more than 16\%.

\subsection{Parsing Speed}
To test the parsing speed, we used an NVidia GTX1080 GPU and parsed the test dataset of UD\_English.  As mentioned in Section~\ref{sec:proposed}, both proposed methods are superior to the baseline model in terms of time complexity. Actually, the methods based on symmetry  took 13.86 seconds, circularity 15.06 and the baseline 17.76 seconds. These results are in accordance with theoretical time complexity.

\section{Conclusion}
In this paper, we reduced the number of parameters in the weight matrices in biaffine classifiers based on the assumption of symmetry or circularity and examined their effects on the CoNLL 2017 shared task for Universal Dependency Parsing dataset. As a result, the method based on circulant matrices outperformed the baseline model in most of languages with about 16\% parameter reduction. As future work, the L1 regularization method for CompLex \citep{Trouillon:2016:CES:3045390.3045609} proposed by \citep{DBLP:conf/aaai/ManabeHS18}  may be integrated into our methods to further reduce the number of parameters in the bilinear function. The script for the bilinear functions proposed here is provided on the GitHub page below\footnote{https://github.com/TomokiMatsuno/PACLIC32/blob/master/my\_linalg.py}.

\if0
\begin{table}[t!]
\begin{center}
\begin{tabular}{|l|r|}
\hline  & 時間(秒) \\
\hline ベースライン & 14.99 \\
対称行列& 15.22 \\
巡回行列 & 15.70 \\
\hline
\end{tabular}
\end{center}
\caption{\label{results} UD\_English 評価データの解析にかかった時間}
 \label{table:parsetime}
\end{table}
\fi

\section{Acknowledgments}
We are grateful to Michael Wentao Li for proofreading the present paper. We thank the anonymous reviewers. This work was supported by JSPS KAKENHI Grant Number JP18K11457 and JST CREST Grant Number JPMJCR1513, Japan.

\bibliographystyle{plainnat}
{\footnotesize