\section{Introduction}

The \gaia latest release of Non-Single Star catalogs
\citep[][hereafter {\it NSS}]{NSS}
includes the orbits of $181\,327$ single-lined spectroscopic binaries (SB1),
based on the radial velocities (RVs) obtained by the space-mission RVS spectrograph
\citep[see also][]{RVS_I_22,RVS_II_22,DR3_Katz}. The sample (hereafter {\it NSS-SB1}) is much larger than any previously-known catalog of SB1s. For example, the SB9 catalog \citep[e.g.,][]{SB9} lists $\sim 4\,000$ orbits, while in a recent work based on the APOGEE project, \cite{Price-Whelan20} have identified $\sim 1\,000$ orbits.

By increasing the known SB1s by two orders of magnitudes, the new catalog is a gold-mine candidate for learning about the statistical features of short-period binaries, like the eccentricity-period relation \citep[e.g.,][]{Mazeh2008,jorissen09}, the frequency of binaries as a function of the  primary mass \citep[][]{raghavan10,troup16, moe17} and mass-ratio distributions \citep[e.g.,][]{mg92,boffin12, boffin15,shahaf17}. These features have profound implication on our understanding of binary formation and evolution \citep[e.g.,][]{verbunt95,bate97,bate02,harada21}, and therefore were intensively discussed in the past \citep[e.g.,][]{duchene13,shahaf19}, based on the relatively small samples available then. The field is now open for new detailed studies based on the much-larger \gaia-SB1 catalog.

The \gaia catalog includes the Keplerian elements of the orbits, but not the \gaia RVs and their epochs themselves.
Instead, the catalog includes a few parameters that characterize the robustness of each orbit, like the \texttt{significance}, defined as the ratio between the primary semi amplitude $K$ and its uncertainty, and the \texttt{rv\_expected\_sig\_to\_noise}, defined as the signal-to-noise-ratio (SNR) of the \gaia RVS-spectra used.
As detailed by {\it NSS}, these parameters can be used to identify the more robust orbits and construct a sub-sample with a lower degree of false-orbit contamination.

Such an exercise was performed by the authors of {\it NSS} themselves, when they considered the eccentricity-period diagram. The full sample includes very short-period binaries with high eccentricity (VSPHE), expected to be rare
\citep[e.g.,][]{meibom05,Mazeh2008,winn16,price-whelan18,terquem21,barker22}.
When {\it NSS} considered only SB1s with
\texttt{significance}$>40$, these VSPHE systems disappeared.
Such a rejection of suspected binaries always comes with a price tag. For example,  the \texttt{significance}$>40$ threshold left a restricted sub-sample of only $28\,930$, smaller by a factor of $\sim6$ than the original one.

Quality assessment and validation using other data sets were also performed by {\it NSS}
\citep[see \gaia DR3 documentation][]{rimoldini22}.\footnote{https://gea.esac.esa.int/archive/documentation/GDR3/pdf}
They presented a comparison between the \gaia DR3 orbital parameters and those from external catalogues --- SB9 \citep{SB9};
APOGEE \citep{Price-Whelan20} and WISE \citep{Petrosky21},
with success rates of $\sim 80\%$.
As expected, when considered
systems with orbital periods shorter than the RVS-observations time span, the recovery rate exceeds $90$\%.

In this work, we
use two external sources of information --- the LAMOST DR6
RV\footnote{https://dr6.lamost.org/} \citep{cui12}
and the Galactic Archaeology
with HERMES (GALAH) DR3\footnote{https://www.galah-survey.org/dr3/the\_catalogues/}
\citep{Buder21} surveys.
Instead of comparing orbital elements, we compare single RVs of $17\,563$ and $6\,018$ LAMOST and GALAH stars with the velocities expected by the \gaia orbits.
\cite{DR3_Katz} used a similar approach to compare the whole \gaia single RVS catalog with a few RV external catalogs.

We could validate the \gaia orbit for most cases, but encountered about $10$--$20$\% inconsistent RVs. Using these data, we constructed a function that estimates the probability of each of the {\it NSS-SB1} orbit to be a false solution, based on the \gaia period and a few of the published robustness parameters. Choosing a working point that allows for $\sim10$\% contamination, we were able to select a 'clean' catalog, consisting of \Nclean orbits.

Section~\ref{sec:comparison} compares the \gaia orbits with LAMOST and GALAH RVs. Section~\ref{sec:validation} suggests a clean SB1 sub-sample and presents some of its statistical aspects, and Section~\ref{sec:statistics} discusses two possible features of the sub-sample. We then briefly summarize and discuss our results in Section~\ref{sec:summary}.

\section{Comparing GAIA SB1 orbits
with LAMOST and GALAH RVs}
\label{sec:comparison}

\subsection{\gaia SB1s Cross match with LAMOST and GALAH}

We began by cross matching the LAMOST DR6 sample of $5\,776\,260$ RVs with the \texttt{gaia\_DR3\_source} catalog, using a sky match of < $1$ arcsec. Next, we removed sources with missing estimates of the RV uncertainty and low signal-to-noise ratio (SNR) in LAMOST g
(\textit{snrg < 10})
and i filters (\textit{snri < 10}), leaving us with $5\,146\,355$ LAMOST RVs associated with \gaia identified stars.

To avoid cases where the LAMOST pipeline has mistakenly grouped different observed sources as the same source, we performed a few quality cuts for sources having more than one observation ($2\,091\,842$ RVs), using the LAMOST parameters $T_{\mathrm{eff}}$, $\log{g}$ and $\mathrm{[Fe/H]}$ which were derived for each exposure.
We removed RVs that one or more of their LAMOST parameters deviated by more than $3\sigma$  from the distribution of this parameter, for the corresponding star.
This left us with  $5\,074\,002$ LAMOST RVs of $3\,850\,858$ different \gaia sources.

We then cross matched this sample, based on \gaia source ID, with the \texttt{nss\_two\_body\_orbit} table of SB1's ($181\,327$ sources with \texttt{nss\_solution\_type = SB1}).
We found $17\,436$ systems in common, of which $12\,230$ objects had a single LAMOST RV measurement, as shown in Fig.~\ref{fig:Nrvs LAMOST-SB1 Hist}.
Three sources
had $20$, $23$ and $31$ observations, for which
we compared the independently-derived LAMOST orbits with those of {\it NSS-SB1} (see subsection \ref{three-systems}).

\begin{figure}
\includegraphics[width=\columnwidth]{graphics/Nrvs_Hist_LAMOST.pdf}
\includegraphics[width=\columnwidth]{graphics/Nrvs_Hist_GALAH.pdf}
    \caption{Histogram of RV measurements per star for LAMOST DR6 (upper panel) and GALAH DR3 (lower panel) sources, and for those that were found with \textit{Gaia} SB1 orbits (blue and green).}
    \label{fig:Nrvs LAMOST-SB1 Hist}
\end{figure}

A similar cross match was performed with the GALAH DR3 \citep{Buder21} RVs.
We considered only measurements with available non gravitational-redshift-corrected RV (\textit{rv\_nogr\_obst}),
resulting in a sample of $562\,887$ RV measurements of $519\,668$ distinct sources.
We found $6\,018$ stars in common with the \gaia orbits,
$5\,552$ of which had a single GALAH-RV measurement,
as seen in Fig.~\ref{fig:Nrvs LAMOST-SB1 Hist}.

\begin{figure*}
	\includegraphics[width=14 cm]{graphics/SKY_SB1_LAMOST_GALAH.pdf}
\caption{Sky projection of $17\,436$
LAMOST (blue) and  $6\,018$ GALAH (green) observed stars with \gaia orbits.}
    \label{fig:SKY}
\end{figure*}

A density distributions of the LAMOST and GALAH systems in common with the \gaia SB1 binaries are plotted in equatorial coordinates in Fig.~\ref{fig:SKY}.
One can see that the LAMOST sources (marked in blue) are spread over the northern hemisphere while the GALAH sources (marked in green) are located in the south. $357$ sources are common to the two samples.

\subsection{Three LAMOST orbits}
\label{three-systems}

The best way of validating the \gaia orbits, as was done by {\it NSS} themselves, is to compare the \gaia Keplerian orbits with other-surveys' orbits. Unfortunately, this can not be done with the LAMOST or GALAH data sets as these surveys, by their nature, did not accumulate enough RVs for almost all of their observed stars.
Nevertheless, as can be seen in
Fig.~\ref{fig:Nrvs LAMOST-SB1 Hist}, there are three LAMOST exceptions for which we could derive independent orbital solutions, as detailed in Table~\ref{tab:LASMOT+20sources} of Appendix~\ref{sec:three-solutions}. Furthermore, we plotted the LAMOST RVs and \gaia orbits, phase folded by the \gaia period, by calculating the  LAMOST phases (including uncertainties) at the \gaia orbits (Figs.~\ref{fig:GaiaModel_3425964192178375680}--\ref{fig:GaiaModel_3376949338201658112})

We could verify the \gaia orbit of only one of the three systems --- \textit{Gaia} DR3 $3425964192178375680$. The LAMOST RVs of another system, \textit{Gaia} DR3 $3376949338201658112$, (see Fig.~\ref{fig:GaiaModel_3376949338201658112}) did not show any variation, and the third system,
\textit{Gaia} DR3 $3425556586900263424$,
yielded a similar but different-period solution.

Although a very small number, the three examples demonstrate two major points.
First, some of the \gaia orbits are probably false with no real variation. However, even for the case of \textit{Gaia} DR3 $3376949338201658112$, a few LAMOST RVs are consistent with the \gaia predicted velocities. Those are close to $\gamma$, the center-of-mass velocity of the orbit. Consequently, one can not verify the \gaia orbit with one LAMOST RV, unless that velocity is consistent with the orbit and significantly different from the \gaia $\gamma$ (see below).
Second, the orbital period of \textit{Gaia} DR3 $3425556586900263424$ is wrong. Nevertheless, most, if not all, the LAMOST RVs are within $1$--$2\sigma$ away from the \gaia expected velocities. Again, one LAMOST RV consistent with the \gaia orbit is no proof that all elements of the orbit are valid.

\newpage
\subsection{Comparing the \textit{Gaia} orbits with LAMOST and GALAH RVs}

To compare the LAMOST velocities with the \gaia orbits, we used the reported \textit{Gaia} Keplerian parameters ($P$, $K$, $e$, $t_p$, $\omega$, $\gamma$) and derived the expected RV --- $RV_{\rm SB1}$, at the LAMOST epoch, MJD $t_L$, and compared it
to the LAMOST observed RV, $RV_{\rm L}$, with the difference
\begin{equation}
   \mathcal{D}_{\rm SB1,L} = RV_{\rm SB1} - RV_{\rm L} \ .
\end{equation}

To derive the uncertainty of $RV_{\rm SB1}$, we used the reported \gaia-Keplerian-parameters correlation matrix, \texttt{corr\_vec}, to produce
$1\,000$ random draws of Keplerian parameters, excluding cases with negative or larger-than-one eccentricity value.
The standard deviation (STD) of these $1\,000$ velocities was considered as the $RV_{\rm SB1}$ uncertainty --- $\sigma[RV_{\rm SB1}]$.  The uncertainty of the difference was then derived by quadrature
of the uncertainties of \gaia-expected and LAMOST velocities:

\begin{equation}
    \sigma[\mathcal{D}_{\rm SB1,L}] =\sqrt{\sigma [RV_{SB1}]^2 + \sigma[ RV_{L}]^2} \ .
\end{equation}
We denote the difference $\mathcal{D}_{\rm SB1,L}$ in units of its own uncertainty as
\begin{equation}
    \mathcal{D}_{\rm SB1,L}^{\sigma}=
    \frac{\mathcal{D}_{\rm SB1,L}}
         {\sigma[\mathcal{D}_{\rm SB1,L}]}\ .
\end{equation}
Similar procedure was performed with the GALAH RVs,
obtaining $\mathcal{D}_{\rm SB1,G}$, $\sigma[\mathcal{D}_{\rm SB1,G}]$ and
$\mathcal{D}_{\rm SB1,G}^{\sigma}$.

\begin{figure*}
\includegraphics[width=18cm]{graphics/SB1_LAMOST_rv_diff.pdf}
\includegraphics[width=18cm]{graphics/SB1_GALAH_rv_diff.pdf}
    \caption{
    Differences between the LAMOST and GALAH RVs and the \gaia orbits.
    The left two panels present the RVs reported by LAMOST (upper panel; blue) and GALAH (bottom panel; green) as a function of the \gaia expected RVs.
    The middle panels show those velocities after subtracting \gaia $\gamma$, the estimated center-of-mass velocity of each binary.
    The right panels show histograms of the difference between \textit{Gaia} RV model and the LAMOST DR6 (upper) and GALAH DR3 (lower) measurements. Red dashed vertical lines mark the median values, suggesting offsets of $4.7$ km/s and $-0.28$ km/s of LAMOST and GALAH RVs, respectively, relative to the \gaia ones. }
    \label{fig:SB1_LAMOST_GALAH_rv_diff}
\end{figure*}

The differences between the LAMOST and GALAH RVs and the \gaia orbits are plotted in Fig.~\ref{fig:SB1_LAMOST_GALAH_rv_diff}.
The left two panels present the RVs reported by LAMOST (upper panel; blue) and GALAH (bottom panel; green) as a function of the \gaia expected RVs. These panels show a general agreement between the orbits and the observed RVs. The middle panels show those velocities after subtracting the \gaia $\gamma$ velocity of each binary, displaying large differences for some of the systems. Histograms of the RV differences are plotted in the right two panels.
The figure does not include the largest differences, in order to present the behaviour around zero difference. The median of the LAMOST-\gaia differences, denoted by a red dashed line, is $4.74$ km/s, suggesting a LAMOST RV zero-point offset. Similar shifts were noticed by \cite[][based on LAMOST DR3]{Anguiano18} and by \cite[][based on LAMOST DR7]{DR3_Katz}. For GALAH, a $-0.28$ km/s shift is emerging.

Fig.~\ref{fig:SB1_LAMOST_rv_dif_fix} shows in its two left panels histograms of the adopted uncertainties of the LAMOST and GALAH RVs, together with the uncertainties of the corresponding \textit{Gaia}-expected velocities.
Based on \citet{Buder21}, we opted to add in quadrature $0.5$,km/s to the reported uncertainty of the GALAH RVs.
The right two panels of the figure show histograms of the RV differences (after correcting for the LAMOST and GALAH zero points), in units of their uncertainties, $\sigma[\mathcal{D}_{\rm SB1,L}]$ and  $\sigma[\mathcal{D}_{\rm SB1,G}]$. We also plot a standard normal distribution expected for the same number of measurements.

The histograms present large wings that indicate a substantial number of non-Gaussian differences, probably originating from systematic errors. The fact that the significant differences appear in both the LAMOST and GALAH histograms suggests that most of them stem from false \gaia solutions, which is the main interest of this work, although we cannot rule out some cases for which the LAMOST or GALAH RVs are wrong. Fig.~\ref{fig:SB1_CornerRatios} shows that the significant RV differences occur in specific \gaia orbits, corroborating the notion that the differences are due to mistaken orbits.

\begin{figure*}
	\includegraphics[width=14cm]{graphics/SB1_LAMOST_rv_dif_fix.pdf}
		\includegraphics[width=14cm]{graphics/SB1_GALAH_rv_dif_fix.pdf}
    \caption{Left panels show histograms of RV uncertainties of the LAMOST (upper left, blue) and GALAH (bottom left) measurements. Also plotted histograms of the \gaia uncertainties of those orbits (grey). Right panels present histograms of the differences, in their uncertainty units, between \textit{Gaia} models and $24\,757$ LAMOST (upper right, blue) and $6\,628$ GALAH (lower right, green) measurements.
    Right-panels grey curves show simulated normal distributions with a mean of $0$ and STD of $1$. The outliers-immune, Median-Absolute-Deviation (MAD) based, STD of the differences distribution for LAMOST is $1.07$ while for GALAH is $1.20$.}
    \label{fig:SB1_LAMOST_rv_dif_fix}
\end{figure*}

\section{Validation the \gaia orbits}
\label{sec:validation}

\subsection{Identification of validated orbits}

To distinguish between true and false \gaia orbits,
we compare the LAMOST and GALAH RVs with two competing models:

(i) The  \gaia orbits, where we use  $\mathcal{D}_{\rm SB1,L}^{\sigma}$ and  $\mathcal{D}_{\rm SB1,G}^{\sigma}$ to quantify the agreement between the RVs and the expected \gaia RVs in LAMOST and GALAH epoch.

(ii) No RV modulation at all; the \gaia stellar velocity is simply the $\gamma$ velocity of the orbit.
For this model, we construct the differences for LAMOST RVs as
\begin{equation}
   \mathcal{D}_{\rm \gamma,L} = \gamma - RV_{\rm L} \ ,
\end{equation}
\begin{equation}
    \sigma[\mathcal{D}_{\rm \gamma,L}] =
    \sqrt{\sigma[\gamma]^2 + \sigma[RV_{L}]^2} \ ,
\end{equation}
and
\begin{equation}
    \mathcal{D}_{\rm \gamma,L}^{\sigma}=
    \frac{\mathcal{D}_{\rm \gamma,L}}
    {\sigma[\mathcal{D}_{\rm \gamma,L}]}\ .
\end{equation}
Similar differences were derived for the GALAH RVs:\\ $\mathcal{D}_{\rm \gamma,G}^{\sigma}=
\mathcal{D}_{\rm \gamma,G}/\sigma[\mathcal{D}_{\rm \gamma,G}]$.

Fig.~\ref{fig:SB1_Sigmas} displays the differences relative to the two models, in units of the corresponding uncertainties, for both the LAMOST and GALAH samples.

Accordingly, we divided our sample of \gaia SB1-LAMOST/GALAH sources into three groups:
\begin{itemize}
    \item $2\,445$ and $2\,777$ systems for which the radial velocity predicted by the \gaia orbits fits well with the LAMOST and GALAH RVs, respectively, and are {\it not} consistent (by more than $1\sigma$) with the $\gamma$ velocity.

    We consider these orbits as validated by LAMOST and GALAH RVs.
    \item $975$ and $666$ orbits for which at least one of the correspsonding LAMOST and GALAH RVs are  more than $3\sigma$ away from the \gaia orbits,

    considered as refuted by LAMOST and GALAH.
    \item All the rest are systems for which the LAMOST and GALAH RVs are within $3\sigma$ from both models or between $1$--$3\sigma$ away from the \gaia orbits,

    and therefore their RVs are inconclusive.
   \end{itemize}

The separations between the different groups are marked in the figure.

\begin{figure*}
	\includegraphics[width=8.5cm]{graphics/SB1_LAMOST_Sigmas.pdf}
		\includegraphics[width=8.5cm]{graphics/SB1_GALAH_Sigmas.pdf}
\caption{Differences (in units of their certainties) between the LAMOST (left) or GALAH (right) RVs and two competing \gaia models --- constant RV and {\it NSS-SB1} model.  Points inside the blue or green rectangle ($2\,445$ or $2\,777$) are sources where the LAMOST or GALAH RVs agree, by up to $1 \sigma$, with the {\it NSS-SB1} model and disagree, by over $3 \sigma$, with the constant model. Points right to the red dashed lines ($975$ or $666$) are sources where the LAMOST or GALAH RVs are inconsistent, by over $3 \sigma$, with the {\it NSS-SB1} model. The rest of the black points in the background are inconclusive sources.}
    \label{fig:SB1_Sigmas}
\end{figure*}

\begin{figure*}
	\includegraphics[width=19cm]{graphics/SB1_CornerRatios.pdf}
    \caption{Corner plot of the relative number of refuted orbits as a function of six published  parameters: Period, N$_{\rm RV}$, nGoF,
    Range/2K, GoF, and eSNR. Each 2D plot is divided into bins, with colours representing the fraction of refuted orbits in that bin. White bins have less than $5$ orbits, considered as not having enough information.}
    \label{fig:SB1_CornerRatios}
\end{figure*}

\subsection{The robust parameters and the validated orbits}

We now try identifying the \gaia robust parameters of the SB1 orbits that can predict the correctness of the solutions. We consider a few of the parameters listed at the \texttt{nss\_two\_body\_orbit} and the \texttt{gaia\_source} tables:

(i) \texttt{Period} --- The SB1 period (Period).

(ii) \texttt{rv\_n\_good\_obs\_primary} --- Total number of epoch RVs actually used for the primary in solving the SB1 model (N$_{\rm RV}$).

(iii) \texttt{rv\_renormalised\_gof} --- An empirical value defined by the \gaia team, that compares the scatter of the RVs relative to the obtained solution with the typical RV uncertainty of stars with similar  \texttt{rv\_template\_teff} and \texttt{grvs\_mag}. As part of the compilation of the NSS table, only sources with \texttt{rv\_renormalised\_gof} > $4$ where considered (nGoF).

(iv) \texttt{rv\_amplitude\_robust}/$2\times$\texttt{semi\_amplitude\_primary} --- The total range of the observed RVs, given by the {\it NSS-SB1} catalog, divided by twice the semi amplitude of the SB1 model. This parameter estimates the coverage of observed RVs of the orbital RV range (Range/2K).

(v) \texttt{goodness\_of\_fit} --- Goodness-of-fit statistic of the solution. This is the ‘gaussianized chi-square’ \citep{rimoldini22} cube root transformation), which for good fits should approximately follow a normal distribution with zero mean value and unit standard deviation (GoF).

(vi) \texttt{rv\_expected\_sig\_to\_noise} --- Expected signal to noise ratio \citep{rimoldini22} in the spectra used to obtain the radial velocity (eSNR).

Other parameters such as the \texttt{conf\_spectro\_period} (i.e. the probability of the period not being due to Gaussian white noise) or more complicated correlation between the orbit parameters were also been tested, but did not show any clear trend.

For each pair of these parameters, we derived a 2D binned map of the ratio of false orbits to the sum of the false and true orbits, presented in Fig.~\ref{fig:SB1_CornerRatios} in a corner-plot form. Obviously, the figure that presents the fraction of the refuted orbits depends on the separation between the false and valid orbits we adopted for Fig.~\ref{fig:SB1_Sigmas}.

Five out of the six robust parameters were found as having a significant predictive power ---
Period, N$_{\rm RV}$,
nGoF, Range/2K and eSNR.
In particular, the relative number of refuted orbits is high for short-period binaries, either with small nGoF, eSNR or N$_{\rm RV}$ parameters.

We note that the tendency for false orbits revealed by the figure is common to the LAMOST and GALAH RVs (see above) in all regions of the parameter space. One apparent exception is the long-period  orbits, for which the LAMOST fraction of RVs inconsistent with the \gaia orbits is larger than the corresponding GALAH RVs. The cause  for this small difference is not clear and might be due to the larger fraction of giants in the long-period orbits.

\begin{table}
\caption{
Best-fitted coefficients of the
Logistic-regression classification.
\label{tab:LR_coeff}}
\begin{center}
\begin{tabular}{l l r }
\hline
Coefficient & Parameter name & Value \\
\hline
$b_0$ &  & $3.8791$ \\
$b_1$ & log(Period) & $-0.6062$  \\
$b_2$ & N$_{\rm RV}$& $-0.041$\\
$b_3$ & log(nGoF) & $-2.9939$  \\
$b_4$ & Range/2K& $0.0153$  \\
$b_5$ & GoF & $ 0.4038$ \\
$b_6$ & eSNR & $0.9734$\\

\hline
\end{tabular}
\end{center}
\end{table}

\begin{figure}
	\includegraphics[width=9cm]{graphics/per_nGOF_allscores.pdf}
    \caption{Mean score as a function of the period and nGoF, presented as 2D binned diagram, color coded by the mean logistic-function score at each bin.}
    \label{fig:per_nGOF_allscores}
\end{figure}

\begin{figure}
	\includegraphics[width=9cm]{graphics/Score_hist.pdf}
    \caption{Histogram of logistic-regression scores of the full {\it NSS-SB1} sample. Low-score systems are probably false.}
    \label{fig:scores_hist}
\end{figure}

\subsection{Probability Function to identify the real \gaia orbits}

Following the qualitative trends presented in the figure, we built a statistical model using logistic regression
\citep{Bishop2006}\footnote{sklearn.linear\_model.LogisticRegression module in python} to estimate the probability that a \gaia SB1 solution is real, using the two groups of
Fig.~\ref{fig:SB1_Sigmas}.
For orbits common to both LAMOST and GALAH, we chose to list them as refuted orbits when one or more of the RVs met our $3\sigma$ away from the \gaia orbit model criteria. We then built a training set sample, composed of $70\%$ of the sample, and a testing set (the other $30\%$ of the sample), and used a \texttt{sigmoid} function, commonly used in logistic-regression framework, to find the best $N$-dimensional plane that separates the two groups.
Specifically, the function derived is:
\begin{equation}
S(z) = \frac{1}{1+e^{\textbf{z}}} ,
\label{eq:logit}
\end{equation}
where $\textbf{z} = b_0 + b_1z_1 + b_2z_2 + b_3z_3 + b_4z_4 + b_5z_4 + b_6z_6$ are the classification best fitted coefficients, listed in Table~\ref{tab:LR_coeff}.

The derived logistic function assigns a score value between $0$ and $1$ that estimates the validity for all orbits of {\it NSS-SB1}. These values are released in the supplement of this paper, for the community to use.
Again, the derived function  depends on the separation between the false and valid orbits we adopted.

Fig.~\ref{fig:per_nGOF_allscores} shows the distribution of the score in the period-nGoF plane, with color code based on mean scores in each bin, while
Fig.~\ref{fig:scores_hist} presents the distribution of the derived scores. A clear trend is evident in
Fig.~\ref{fig:per_nGOF_allscores} that resembles the shape seen in the associate diagram of Fig.~\ref{fig:SB1_CornerRatios}.

\newpage
\subsection{The Clean Sample}

To characterize the diagnostic ability of our binary
classification we used a ROC (Receiver Operating Characteristic) curve \citep{fawcett06} presented in Fig.~\ref{fig:ROC}. A ROC graph is a technique for visualizing, organizing and selecting classifiers based on their performance.  ROC curve is created by plotting the estimated True Positive Rate (TPR) against the estimated False Positive Rate (FPR) at various threshold settings.  The TPR (sensitivity) is a measure of the fraction of false orbits that are identified as such.

\begin{figure}
	\includegraphics[width=8cm]{graphics/ROC.pdf}
    \caption{ROC curve and our selected score limit of $0.587$ (red point), aimed to set a True Positive Rate (TPR) of $0.80$. The Area Under the Curve (AUC) is $0.78$, clearly above the $0.5$ value expected for a random-model classifier, marked by the red dashed line.}
    \label{fig:ROC}
\end{figure}

Adopting a score limit of $0.587$, with a sensitivity of $80\%$, implies that our 'clean' sample, with \Nclean orbits, might include $20\%$ of the false orbits.
Examining the testing sample, we estimate that the contamination of the clean sample is at an averaged level of $\sim10$\%, with a slightly higher level for long-period orbits.

Table~\ref{tab:ROC} lists a few representative threshold values and their corresponding TPR and FPR. These values can be used to select a different SB1 sample with desired sensitivity levels.

\begin{figure*}
	\includegraphics[width=17cm]{graphics/ParamsHistComparison.pdf}
    \caption{Histogram of the general {\it NSS-SB1} sample ($180\,281$ sources; gray) , the restricted sample of  \texttt{significance > 40} ($28\,774$ sources; red), and our clean \Nclean sources (blue) as function of selected \gaia parameters: Period, \texttt{significance}, N$_{\rm RV}$, \texttt{significance} $\times$ N$_{\rm RV}$, nGoF, G magnitude, GoF, eSNR and distance ($1/\mathrm{parallax}$ [Kpc]). }
    \label{fig:ParamsHistComparison}
\end{figure*}

\begin{table}
\caption{TPR and FPR values for different threshold scores of our logistic-regression clustering algorithm.
\label{tab:ROC}}
\begin{center}
\begin{tabular}{r r r }
\hline
TPR & FPR & Threshold \\
\hline
$0.11$ & $0.01$ & $0.146$ \\
$0.16$ & $0.01$ & $0.158$ \\
$0.21$ & $0.02$ & $0.185$ \\
$0.25$ & $0.02$ & $0.201$ \\
$0.31$ & $0.03$ & $0.228$ \\
$0.35$ & $0.03$ & $0.243$ \\
$0.4$ & $0.05$ & $0.267$ \\
$0.45$ & $0.06$ & $0.3$ \\
$0.5$ & $0.09$ & $0.337$ \\
$0.55$ & $0.14$ & $0.386$ \\
$0.6$ & $0.19$ & $0.422$ \\
$0.65$ & $0.24$ & $0.458$ \\
$0.7$ & $0.29$ & $0.493$ \\
$0.75$ & $0.34$ & $0.526$ \\
$0.8$ & $0.45$ & $0.587$ \\
$0.85$ & $0.55$ & $0.643$ \\
$0.9$ & $0.65$ & $0.695$ \\
$0.95$ & $0.76$ & $0.76$ \\

\hline
\end{tabular}
\end{center}
\end{table}

To demonstrate the realization of our method, Fig.~\ref{fig:ParamsHistComparison} presents the distribution of some of the robust parameters listed for the SB1 solutions for
three samples:
\begin{itemize}
    \item the {\it NSS-SB1} sample with $180\,281$ orbits with all six robust parameters released (gray curve),
    \item the restricted sample of  \texttt{significance > 40} of $28\,774$ orbits (red curve), considered by {\it NSS}, and
    \item our clean sample of \Nclean orbits (blue curve).
\end{itemize}

As can be seen, our clean sample is twice smaller than the general {\it NSS-SB1} one, yet larger by a factor of  $\sim 3$ than the \texttt{significance > 40} sample.
Most notably is that the clean sample includes orbits with lower \texttt{significance} and nGoF, as opposed to the  \texttt{significance > 40} sample. In addition, our clean sample seems to be centered around GoF of $0$, as expected by the SB1 solutions \citep{rimoldini22}.

\subsection{Eccentricity-Period Diagram}
\label{e-p}

To demonstrate the efficiency of our classification, we consider the eccentricity-period diagram, originally discussed by {\it NSS} and pointed out in the introduction of this paper. As in {\it NSS},  Fig.~\ref{fig:ecc_per} shows that the full sample (upper panel; black) does show many very short-period binaries with high eccentricity (VSPHE), while the restricted sample of \texttt{significance > 40} (lower panel; red) does not include any of those (lower panel). The middle panel (blue) with our clean sample, although much larger than the restricted one, is also clean of these apparently false orbits.

Out of all binaries of the clean sample, only the eccentricity of one system is significantly higher than the \citet{Mazeh2008} upper-envelope line. The system --- \gaia~DR3~$6727359218822543872$, with reported orbital period $P=9.13712\pm 0.00046$ days and eccentricity $e=0.862\pm0.017$, might be an interesting system to follow up, assuming its orbit is real.

\begin{figure*}
\includegraphics[width=14cm]{graphics/ecc_per.pdf}
\caption{Eccentricity vs.~period of the general {\it NSS-SB1} sample ($180\,281$ orbits; upper panel), our clean sample (\Nclean orbits; middle panel) and the restricted sample of  \texttt{significance > 40} ($28\,774$ orbits; bottom panel). Black dotted line marks the observed \citet{Mazeh2008} upper envelope.}
    \label{fig:ecc_per}
\end{figure*}

\section{Two statistical features of the clean SB1 sample}
\label{sec:statistics}

We turn now to plot the corner plot of the period, eccentricity and primary mass of the clean sample.
We wish to derive this plot separately for the  main-sequence (MS) and the evolved stars, as the correlation between these parameters might be different for the two populations.
To separate the MS from the giant stars we
plot in Fig.~\ref{fig:CMD} the positions of the binaries of the whole clean sample on the CMD, where the two populations are clearly resolved by:
\begin{equation}
  \mathrm{MS}:\begin{cases}
    -7.5+10(G_{\mathrm{BP}}- G_{\mathrm{RP}}) < G, & \mathrm{if~} G_{\mathrm{BP}}- G_{\mathrm{RP}} \leq 1.32 \ ,\\
    3.1+2(G_{\mathrm{BP}}- G_{\mathrm{RP}}) < G, & \mathrm{if~} G_{\mathrm{BP}}- G_{\mathrm{RP}} > 1.32 \ ,
  \end{cases}
\end{equation}
and the evolved stars from the long-period variables range by
\begin{equation}
  \mathrm{LPV}:
    -3.4+2(G_{\mathrm{BP}}- G_{\mathrm{RP}}) > G, \\ .
\end{equation}
We get $17\,603$ and $4\,237$ MS and evolved stars, respectively.

\begin{figure*}
	\includegraphics[width=12cm]{graphics/CMD.pdf}
    \caption{CMD of our clean-sample systems for which a primary-mass estimate available.}
    \label{fig:CMD}
\end{figure*}

We then plotted in Fig.~\ref{fig:CornerMS_EV}
the corner plots of the distributions of the period, eccentricity and primary mass of the MS and evolved stars.
Two features are emerging:
\begin{itemize}
    \item The MS period-eccentricity diagram suggests two populations: small-eccentricity binaries, in a period range of $2$--$14$ days, and another population with an eccentricity range that fits the \citet{Mazeh2008} curve. The evolved diagram is similar to the second MS population. One can argue that the small eccentricity population is composed of binaries that went through tidal circularization.
    \item There is a paucity of short-period (up to $5$ day) low-mass (up to $\sim 1.3$ $M_{\odot}$) binaries.
\end{itemize}

These two features, like many other statistical aspects, require further study.

\begin{figure*}
	\includegraphics[width=11cm]{graphics/Corner_CMD_Main.pdf}
		\includegraphics[width=11cm]{graphics/Corner_CMD_EV.pdf}
    \caption{Corner plots of our clean sample of $17\,603$ MS (upper panel) and $4\,237$  evolved (lower panel) stars.}
    \label{fig:CornerMS_EV}
\end{figure*}

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

We used the LAMOST DR6 and GALAH DR3 RVs
to validate the orbits of the {\it NSS} \gaia SB1 catalog, finding that $10$--$20$\% of the orbits are probably false.

Our analysis suggests that many of the short-period orbits are false, as was indicated by the eccentricity-period diagram of the original {\it NSS} paper. The reason for this might be the fact that the short-period range of orbits is densely populated by many possible independent false solutions, some of which are consistent with the small number of \gaia RVs. This is especially true for short-period eccentric orbits. Despite cautionary measures by {\it NSS}, some false orbits apparently made their way to the catalog.
Another cause for false orbits could be the low resolution of \gaia RVS, which did not allow resolving the two components of some double-lined spectroscopic binaries (SB2), introducing wrong solutions. This effect can be stronger in long-period SB2s, for which the RV modulation is smaller, and therefore the two components are not enough separated.

Based on the comparison with LAMOST and GALAH data sets we construct a function that estimates the probability of each of the orbits to be correct, using their published robust parameters.
The function and its values for each of the SB1 are given. Choosing a working point that allows for false-orbit contamination on the order of $10$\%, we were able to put up a selected catalog, consisting of
\Nclean orbits, for the community to use.
We have shown that for our particular choice, there are no very-short binaries with high eccentricity in the clean sample.

For three systems with more than $20$ LAMOST RVs measurements, we independently solved for their orbits.
The three newly derived orbits display different levels of consistency with the \gaia solutions, in agreement with the score derived by our
classifier for those \gaia orbits.

Obviously,
the value of the working point, as well as the border lines between false and validated orbits used to construct the probability function, are somewhat arbitrary. Other border lines and working points can be selected, resulting in different clean catalogs, with varying levels of contamination.

The clean sample offers many opportunities for studying the short-period binary population. For example, the period-eccentricity diagram of the MS stars of the clean sample shows a sub-sample of small-eccentricity binaries at a period range of $2$--$14$ days, probably composed of binaries that went through tidal circularization. This feature appears in the binaries of evolved stars too. The period-primary mass diagram shows
a paucity of short-period (up to $5$ day) low-mass (up to $\sim 1.3$ $M_{\odot}$) binaries.

These two examples, like many other statistical features hidden in the sample, require further extensive study. In particular, one has to consider observational selection effects that might bias the results, including the specific choices of the clean sample. Furthermore, as shown by {\it NSS}, there is an overlap between the eclipsing and the spectroscopic binaries for the short-period systems, and another overlap between the SB1 systems and the astrometric binaries. A global view of the binary population at hand requires a global analysis of the three \gaia
{\it NSS} catalogs \citep[see also][]{rimoldini22}, together with other studies, eclipsing binaries found by large photometric surveys \citep[e.g.,][]{Kepler-EB_16,ASAS_EB_06, ZTF_EB_20, ASAS_EB_22} in particular.

\newpage
\section*{Acknowledgements}

This research was supported by Grant No. 2016069 of the United States-Israel Binational Science Foundation (BSF) and by Grant No. I-1498-303.7/2019 of the German-Israeli Foundation for Scientific Research and Development (GIF) to TM and HWR.

This work has made use of data from the European Space Agency (ESA) mission \gaia (http://www.cosmos.esa.int/gaia), processed
by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the \gaia Multilateral Agreement.

Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission.
LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.

This work made use of the Third Data Release of the GALAH Survey \citep{Buder21}. The GALAH Survey is based on data acquired through the Australian Astronomical Observatory, under programs: A/2013B/13 (The GALAH pilot survey); A/2014A/25, A/2015A/19, A2017A/18 (The GALAH survey phase 1); A2018A/18 (Open clusters with HERMES); A2019A/1 (Hierarchical star formation in Ori OB1); A2019A/15 (The GALAH survey phase 2); A/2015B/19, A/2016A/22, A/2016B/10, A/2017B/16, A/2018B/15 (The HERMES-TESS program); and A/2015A/3, A/2015B/1, A/2015B/19, A/2016A/22, A/2016B/12, A/2017A/14 (The HERMES K2-follow-up program). We acknowledge the traditional owners of the land on which the AAT stands, the Gamilaraay people, and pay our respects to elders past and present. This paper includes data that has been provided by AAO Data Central (datacentral.org.au).

This research made use of ASTROPY,\footnote{http://www.astropy.org} a community-developed core Python package for Astronomy \citep{astropy13, astropy18},
TOPCAT tool, described in \citep{Taylor05},
and the cross-match service provided by CDS, Strasbourg.

\section*{Data Availability}

Data used in this study are available upon request from the corresponding
author.

\bibliographystyle{mnras}


\section{Supplementary}
\subsection{1.Computational details of density functional calculation}
For this prototypical case, we obtain the most relevant Hilbert space within $\pm10$ eV around the Fermi energy from the spin polarized LDA+$U$ ~\cite{Anisimov1993, Liechtenstein1995} electronic structure of the parent compound NdNiO${_2}$, using the linearized augmented plane wave~\cite{Singh} implementation~\cite{Blaha1990} of the density functional theory (DFT)~\cite{DFT1, DFT2}.
We take from Ref.~\cite{Li2019} the lattice structure with lattice constant $a=b=3.92\mathring{A},c=3.37\mathring{A}$, and the space group $P4/mmm$.
We set $U-J$ = $0.6 - 0.06$ Ry = 7.344 eV for Ni $d$-orbitals,$0.20 - 0.00$ Ry = 2.72 eV,  $0.60 - 0.00$ Ry = 8.16 eV for Nd $d$-,$f$-orbitals separately.

\vspace{-0.2cm}
\subsection{2.Extraction of many-body Hamiltonian}
We aim at obtaining a realistic SU(2)-symmetric many-body Hamiltonian $H^{\text{eff}}$ in Eq.(1).
To ensure a DFT-like level of accuracy, we further demand that when under a similar approximation, our many-body Hamiltonian needs to reproduce the self-consistent DFT Hamiltonian (equivalently the electronic band structure).
Since in LDA+$U$ the strong local electron-electron interaction is approximated in an effective Hartree-Fock manner~\cite{Anisimov1997}, a proper connection can be naturally made on the Wannier states basis by matching the self-consistent Hartree-Fock solution of our SU(2)-symmetric $H^{\text{eff}}$ with the self-consistent $H^{\text{LDA}+U}$~\cite{Wei2006}.
The requirement that the Hartree-Fock approximation of $H^{\text{eff}}$ needs to reproduce the LDA+$U$ band structure within the subspace of the active orbitals results in a rather unique set of parameters in $H^{\text{eff}}$~\cite{Wei2006}.

Specifically, since we demand that the self-consistent Hartree-Fock solution reproduce the self-consistent solution of the LDA+$U$ solution, the corresponding density matrix $\rho_{i\nu\nu^\prime}$ must be identical in both cases when represented in the same set of Wannier orbitals.
One can thus take $\rho_{i\nu\nu^\prime}$ directly from the self-consistent LDA+$U$ solution in the Wannier basis.
Furthermore, if one assumes that the structure of $U_{m m^{\prime\prime}m^{\prime\prime\prime} m^\prime}$ follows the Slater integral~\cite{Slater,Liechtenstein1995}, the entire $U_{m m^{\prime\prime}m^{\prime\prime\prime} m^\prime}$ can be fixed by just two parameters $U_{\text{eff}}$ and $J_{\text{eff}}$.
Combining $\rho_{i\nu\nu^\prime}$ and $U_{m m^{\prime\prime}m^{\prime\prime\prime} m^\prime}$, the effective Hartree-Fock potential $V^\text{HF}$ can then obtained.

In addition, orbital dependent $J_m$ in $H^{\text{eff}}$ can be straightforwardly obtained from the local site energy difference of Nd $d$-orbitals between spin-up and spin-down, when the $\frac{3}{2}$-spin $\mathbf{S}_i^f$ of the $f$-orbitals are all set to be along the spin-up direction.
Finally, the remaining intra-atomic hopping parameters $t_{ii mm^\prime}$ and site energy $\epsilon_m$ is then obtained by subtracting $V^\text{HF}$ from the intra-atomic part of $H^{\text{LDA}+U}$.

\begin{figure}[h]

    \vspace{-0.3cm}
	\begin{center}
	\includegraphics[width=0.9\columnwidth]{Reproduce.pdf}
	\end{center}
	\vspace{-0.6cm}
	\caption{
	Confirmation of the quality of $H^{^\mathbf{eff}}$ by comparing the electronic band structures with those from LDA+$U$ calculations under different magnetic configurations: (a)(b) spin majority of ferromagnetic order, (c)(d) spin minority of ferromagnetic order,  and (e)(f) antiferromagnetic order.
	The excellent agreement establishes a robust foundation for configurations with non-collinear spin correlations without long-range order.
	}
	\label{figs1}
	\vspace{-0.3cm}
\end{figure}

A simple criterion to check the quality of the resulting $H^{\text{eff}}$ (or the accuracy of the chosen $U_{\text{eff}}$ and $J_{\text{eff}}$) is the degree of the spin-independence of the intra-atomic $t_{ii mm^\prime}$.
A reasonable value of $U_{\text{eff}}$ and $J_{\text{eff}}$ should encode all the spin-dependence of the self-consistent $H^\text{HF}$.
Therefore, a non-negligible spin-dependence of the resulting $t_{ii mm^\prime}$ or site energy $\epsilon_m$ indicates clearly a need to improve the value of $U_{\text{eff}}$ and $J_{\text{eff}}$.
In practice, we find this criterion quite sufficient to pin down a rather narrow range of acceptable values of $U_{\text{eff}}$ and $J_{\text{eff}}$.
This then allows us to fix all the intra-atomic parameters in $H^{\text{eff}}$.

On the other hand, since the LDA+$U$ only included atomically local Hartree-Fock approximation, the inter-atomic hopping parameters $t_{ii^\prime mm^\prime}$ are unaffected in the approximation.
We therefore can simply take it from the spin-averaged $H^{\text{LDA}+U}$.
The next section lists some of the leading parameters.

The above procedure, if performed properly, should generate a SU(2)-symmetric many-body Hamiltonian $H^{\text{eff}}$ that can be further studied with any many-body solver one prefers, not just the Hartree-Fock approximation.
Furthermore, since the Hamiltonian does not require or guarantee a magnetic order, it can in principle be applied to metallic systems (or spin liquid insulators) that do not contain long-range order but show signs of the existence of local moments, for example Curie-Weiss susceptibility at high-temperature.

For the purpose of this work, in which we study the impact of non-collinear spin correlations, we verify that $H^{\text{eff}}$ is able to reproduce the $H^{\text{LDA}+U}$ solution under various magnetic orders.
Figure~\ref{figs1} demonstrates an excellent correspondence between our resulting band structure under the same Hartree-Fock treatment of local interactions and the LDA+$U$ band structure within the subspace of the selected orbitals.
Notice especially that the agreement occurs under both ferromagnetic and antiferromagnetic order with the \emph{same} set of parameters, despite the significantly different band structure under these two orders.
These results thus establish the high quality of these parameters and the validity of our effective Hamiltonian in describing various magnetic structures, including the paramagnetic parent compound.

\vspace{-0.2cm}
\subsection{3. Leading parameters for $H^{\text{eff}}$}
Table \ref{tabs1} lists some of the leading coefficients $\epsilon_m$ and $J_m$ of the resulting SU(2) symmetric many-body Hamiltonian.
The full interaction kernel $U_{m m^{\prime\prime}m^{\prime\prime\prime} m^\prime}$ is approximated by the Slater integral~\cite{Slater,Liechtenstein1995}, with $U_{\text{eff}} = 7.385$ eV and $J_{\text{eff}}=1.197$ eV obtained from the above procedure.
The leading terms of the resulting hopping parameters $t_{ii^\prime mm ^\prime}$ are $1.321$ eV between Ni $d_{x^2-y^2}$ and O $p_{x}$, $1.158$ eV between Nd $d_{xz}$/$d_{yz}$ and O $p_{z}$, and $1.024$ eV between Nd $d_{xy}$ and ``I'' orbitals.
The full $t_{ii^\prime mm ^\prime}$ parameters are available upon request.

\begin{table}[t]
    \caption{Leading parameters $\epsilon_m$ and $J_m$ in unit of eV.}
    \begin{ruledtabular}
        \begin{tabular}{llllll}
Nd & $d_{3z^2-r^2}$ & $d_{x^2-y^2}$ & $d_{yz}$ & $d_{xz}$ & $d_{xy}$ \\ \hline
$J_m$  & 0.272 & 0.379  & 0.268 & 0.268 & 0.217 \\
$\epsilon_m$ & 4.12 & 5.43  & 5.779 & 5.779 & 4.045 \\ \hline \hline
Ni & $d_{3z^2-r^2}$ & $d_{x^2-y^2}$ & $d_{yz}$ & $d_{xz}$ & $d_{xy}$ \\ \hline
$\epsilon_m$ & -55.448 & -55.101  & -55.616 & -55.616 & -55.284 \\\hline\hline
O & $p_{x}$ & $p_{y}$ & $p_{z}$ & ``I'' \\ \hline
$\epsilon_m$ & -3.488 & -2.542  & -2.602 & 4.93 \\
\end{tabular}
    \end{ruledtabular}
    \label{tabs1}
    \vspace{-0.3cm}
\end{table}

\vspace{-0.2cm}
\subsection{4. Implementation of non-collinear magnetic configuration}

In this study, the fluctuation of non-collinear spin directions is incorporated in our simulation by averaging large number of configurations with similar nearest short-range correlations $\chi = \sum_{\langle i, i^\prime\rangle}\mathbf{S}_{i} \cdot \mathbf{S}_{i^\prime}$,
but with negligible long-range order $\mathbf{M}_{\mathbf{q}} = \sum_{i}\mathbf{S}_{i} e^{i\mathbf{q}\cdot \mathbf{x}_i}$.
Obviously, a larger supercell is necessary in the search for configurations of various short-range correlation, but without long-range order.
Furthermore, averaging over large supercells with various sizes, shapes and orientations is an efficient way to avoid fictitious gap openings and shadow band foldings related to the artificial new spatial periodicity.

The enforcement of the spin directions according to the proposed magnetic configurations can be easily achieved within the Hartree-Fock approximation by constraining the atomically local one-body density matrix $\rho_{i\nu\nu^\prime}$ to be diagonal in the spin channel $\nu=\nu^\prime$ along the assigned spin direction.
Specifically, using the Euler angle we represent (rotate) $\rho_{i\nu\nu^\prime}$ in a local spin basis whose $z$-axis is along the assigned spin direction.
We then zero out the off-diagonal elements of $\rho_{i\nu\nu^\prime}$ in this local basis, and then rotate the representation back to the global one with $z$-axis of the spin along that of the lattice.
This constrained $\rho_{i\nu\nu^\prime}$ is then combined with  $U_{m m^{\prime\prime}m^{\prime\prime\prime} m^\prime}$ in Eq.(1) to evaluate the effective orbital-dependent potential within the Hartree-Fock approximation.
Naturally, the same constraint needs to be applied in each iteration of the self-consistent cycle, until a spin-density is converged.

\section{Supplementary}
\subsection{1.Computational details of density functional calculation}
For this prototypical case, we obtain the most relevant Hilbert space within $\pm10$ eV around the Fermi energy from the spin polarized LDA+$U$ ~\cite{Anisimov1993, Liechtenstein1995} electronic structure of the parent compound NdNiO${_2}$, using the linearized augmented plane wave~\cite{Singh} implementation~\cite{Blaha1990} of the density functional theory (DFT)~\cite{DFT1, DFT2}.
We take from Ref.~\cite{Li2019} the lattice structure with lattice constant $a=b=3.92\mathring{A},c=3.37\mathring{A}$, and the space group $P4/mmm$.
We set $U-J$ = $0.6 - 0.06$ Ry = 7.344 eV for Ni $d$-orbitals,$0.20 - 0.00$ Ry = 2.72 eV,  $0.60 - 0.00$ Ry = 8.16 eV for Nd $d$-,$f$-orbitals separately.

\vspace{-0.2cm}
\subsection{2.Extraction of many-body Hamiltonian}
We aim at obtaining a realistic SU(2)-symmetric many-body Hamiltonian $H^{\text{eff}}$ in Eq.(1).
To ensure a DFT-like level of accuracy, we further demand that when under a similar approximation, our many-body Hamiltonian needs to reproduce the self-consistent DFT Hamiltonian (equivalently the electronic band structure).
Since in LDA+$U$ the strong local electron-electron interaction is approximated in an effective Hartree-Fock manner~\cite{Anisimov1997}, a proper connection can be naturally made on the Wannier states basis by matching the self-consistent Hartree-Fock solution of our SU(2)-symmetric $H^{\text{eff}}$ with the self-consistent $H^{\text{LDA}+U}$~\cite{Wei2006}.
The requirement that the Hartree-Fock approximation of $H^{\text{eff}}$ needs to reproduce the LDA+$U$ band structure within the subspace of the active orbitals results in a rather unique set of parameters in $H^{\text{eff}}$~\cite{Wei2006}.

Specifically, since we demand that the self-consistent Hartree-Fock solution reproduce the self-consistent solution of the LDA+$U$ solution, the corresponding density matrix $\rho_{i\nu\nu^\prime}$ must be identical in both cases when represented in the same set of Wannier orbitals.
One can thus take $\rho_{i\nu\nu^\prime}$ directly from the self-consistent LDA+$U$ solution in the Wannier basis.
Furthermore, if one assumes that the structure of $U_{m m^{\prime\prime}m^{\prime\prime\prime} m^\prime}$ follows the Slater integral~\cite{Slater,Liechtenstein1995}, the entire $U_{m m^{\prime\prime}m^{\prime\prime\prime} m^\prime}$ can be fixed by just two parameters $U_{\text{eff}}$ and $J_{\text{eff}}$.
Combining $\rho_{i\nu\nu^\prime}$ and $U_{m m^{\prime\prime}m^{\prime\prime\prime} m^\prime}$, the effective Hartree-Fock potential $V^\text{HF}$ can then obtained.

In addition, orbital dependent $J_m$ in $H^{\text{eff}}$ can be straightforwardly obtained from the local site energy difference of Nd $d$-orbitals between spin-up and spin-down, when the $\frac{3}{2}$-spin $\mathbf{S}_i^f$ of the $f$-orbitals are all set to be along the spin-up direction.
Finally, the remaining intra-atomic hopping parameters $t_{ii mm^\prime}$ and site energy $\epsilon_m$ is then obtained by subtracting $V^\text{HF}$ from the intra-atomic part of $H^{\text{LDA}+U}$.

\begin{figure}[h]

    \vspace{-0.3cm}
	\begin{center}
	\includegraphics[width=0.9\columnwidth]{Reproduce.pdf}
	\end{center}
	\vspace{-0.6cm}
	\caption{
	Confirmation of the quality of $H^{^\mathbf{eff}}$ by comparing the electronic band structures with those from LDA+$U$ calculations under different magnetic configurations: (a)(b) spin majority of ferromagnetic order, (c)(d) spin minority of ferromagnetic order,  and (e)(f) antiferromagnetic order.
	The excellent agreement establishes a robust foundation for configurations with non-collinear spin correlations without long-range order.
	}
	\label{figs1}
	\vspace{-0.3cm}
\end{figure}

A simple criterion to check the quality of the resulting $H^{\text{eff}}$ (or the accuracy of the chosen $U_{\text{eff}}$ and $J_{\text{eff}}$) is the degree of the spin-independence of the intra-atomic $t_{ii mm^\prime}$.
A reasonable value of $U_{\text{eff}}$ and $J_{\text{eff}}$ should encode all the spin-dependence of the self-consistent $H^\text{HF}$.
Therefore, a non-negligible spin-dependence of the resulting $t_{ii mm^\prime}$ or site energy $\epsilon_m$ indicates clearly a need to improve the value of $U_{\text{eff}}$ and $J_{\text{eff}}$.
In practice, we find this criterion quite sufficient to pin down a rather narrow range of acceptable values of $U_{\text{eff}}$ and $J_{\text{eff}}$.
This then allows us to fix all the intra-atomic parameters in $H^{\text{eff}}$.

On the other hand, since the LDA+$U$ only included atomically local Hartree-Fock approximation, the inter-atomic hopping parameters $t_{ii^\prime mm^\prime}$ are unaffected in the approximation.
We therefore can simply take it from the spin-averaged $H^{\text{LDA}+U}$.
The next section lists some of the leading parameters.

The above procedure, if performed properly, should generate a SU(2)-symmetric many-body Hamiltonian $H^{\text{eff}}$ that can be further studied with any many-body solver one prefers, not just the Hartree-Fock approximation.
Furthermore, since the Hamiltonian does not require or guarantee a magnetic order, it can in principle be applied to metallic systems (or spin liquid insulators) that do not contain long-range order but show signs of the existence of local moments, for example Curie-Weiss susceptibility at high-temperature.

For the purpose of this work, in which we study the impact of non-collinear spin correlations, we verify that $H^{\text{eff}}$ is able to reproduce the $H^{\text{LDA}+U}$ solution under various magnetic orders.
Figure~\ref{figs1} demonstrates an excellent correspondence between our resulting band structure under the same Hartree-Fock treatment of local interactions and the LDA+$U$ band structure within the subspace of the selected orbitals.
Notice especially that the agreement occurs under both ferromagnetic and antiferromagnetic order with the \emph{same} set of parameters, despite the significantly different band structure under these two orders.
These results thus establish the high quality of these parameters and the validity of our effective Hamiltonian in describing various magnetic structures, including the paramagnetic parent compound.

\vspace{-0.2cm}
\subsection{3. Leading parameters for $H^{\text{eff}}$}
Table \ref{tabs1} lists some of the leading coefficients $\epsilon_m$ and $J_m$ of the resulting SU(2) symmetric many-body Hamiltonian.
The full interaction kernel $U_{m m^{\prime\prime}m^{\prime\prime\prime} m^\prime}$ is approximated by the Slater integral~\cite{Slater,Liechtenstein1995}, with $U_{\text{eff}} = 7.385$ eV and $J_{\text{eff}}=1.197$ eV obtained from the above procedure.
The leading terms of the resulting hopping parameters $t_{ii^\prime mm ^\prime}$ are $1.321$ eV between Ni $d_{x^2-y^2}$ and O $p_{x}$, $1.158$ eV between Nd $d_{xz}$/$d_{yz}$ and O $p_{z}$, and $1.024$ eV between Nd $d_{xy}$ and ``I'' orbitals.
The full $t_{ii^\prime mm ^\prime}$ parameters are available upon request.

\begin{table}[t]
    \caption{Leading parameters $\epsilon_m$ and $J_m$ in unit of eV.}
    \begin{ruledtabular}
        \begin{tabular}{llllll}
Nd & $d_{3z^2-r^2}$ & $d_{x^2-y^2}$ & $d_{yz}$ & $d_{xz}$ & $d_{xy}$ \\ \hline
$J_m$  & 0.272 & 0.379  & 0.268 & 0.268 & 0.217 \\
$\epsilon_m$ & 4.12 & 5.43  & 5.779 & 5.779 & 4.045 \\ \hline \hline
Ni & $d_{3z^2-r^2}$ & $d_{x^2-y^2}$ & $d_{yz}$ & $d_{xz}$ & $d_{xy}$ \\ \hline
$\epsilon_m$ & -55.448 & -55.101  & -55.616 & -55.616 & -55.284 \\\hline\hline
O & $p_{x}$ & $p_{y}$ & $p_{z}$ & ``I'' \\ \hline
$\epsilon_m$ & -3.488 & -2.542  & -2.602 & 4.93 \\
\end{tabular}
    \end{ruledtabular}
    \label{tabs1}
    \vspace{-0.3cm}
\end{table}

\vspace{-0.2cm}
\subsection{4. Implementation of non-collinear magnetic configuration}

In this study, the fluctuation of non-collinear spin directions is incorporated in our simulation by averaging large number of configurations with similar nearest short-range correlations $\chi = \sum_{\langle i, i^\prime\rangle}\mathbf{S}_{i} \cdot \mathbf{S}_{i^\prime}$,
but with negligible long-range order $\mathbf{M}_{\mathbf{q}} = \sum_{i}\mathbf{S}_{i} e^{i\mathbf{q}\cdot \mathbf{x}_i}$.
Obviously, a larger supercell is necessary in the search for configurations of various short-range correlation, but without long-range order.
Furthermore, averaging over large supercells with various sizes, shapes and orientations is an efficient way to avoid fictitious gap openings and shadow band foldings related to the artificial new spatial periodicity.

The enforcement of the spin directions according to the proposed magnetic configurations can be easily achieved within the Hartree-Fock approximation by constraining the atomically local one-body density matrix $\rho_{i\nu\nu^\prime}$ to be diagonal in the spin channel $\nu=\nu^\prime$ along the assigned spin direction.
Specifically, using the Euler angle we represent (rotate) $\rho_{i\nu\nu^\prime}$ in a local spin basis whose $z$-axis is along the assigned spin direction.
We then zero out the off-diagonal elements of $\rho_{i\nu\nu^\prime}$ in this local basis, and then rotate the representation back to the global one with $z$-axis of the spin along that of the lattice.
This constrained $\rho_{i\nu\nu^\prime}$ is then combined with  $U_{m m^{\prime\prime}m^{\prime\prime\prime} m^\prime}$ in Eq.(1) to evaluate the effective orbital-dependent potential within the Hartree-Fock approximation.
Naturally, the same constraint needs to be applied in each iteration of the self-consistent cycle, until a spin-density is converged.