\section{Introduction}
Lattice QCD simulations have become an indispensable tool to explore the non-perturbative baryon-baryon interactions, in particular, hyperon-nucleon and hyperon-hyperon interactions~\cite{Aoki:2020bew,Illa:2020nsi,Wagman:2017tmp,Inoue:2010es,Beane:2010hg,Ishii:2006ec}, which have so far remain difficult to be studied experimentally. Nonetheless, lattice QCD simulations of coupled channel baryon-baryon interactions are very challenging as well (see, e.g., ~Refs.~\cite{Briceno:2013bda,Sasaki:2019qnh}). In this sense,  reliable theoretical studies of lattice QCD simulations are in urgent need, particularly for those simulations performed for unphysical pion masses, which remain unconstrained. One of such theoretical tools is chiral perturbation theory~\cite{Weinberg:1978kz,Gasser:1983yg,Gasser:1987rb}(see Ref.~\cite{Scherer:2012xha} for a review) and the resulting chiral nuclear forces~\cite{Weinberg:1990rz,Weinberg:1991um,Ordonez:1995rz,Bedaque:2002mn,Epelbaum:2008ga,Machleidt:2011zz}. It should be noted that the conventional chiral nuclear forces are based on the heavy baryon (HB) version of chiral perturbation theory~\cite{Jenkins:1990jv,Bernard:1992qa}.    In recent years, it has been proposed that one can build a covariant verion of chiral nuclear forces with the covariant baryon chiral perturbation theory~\cite{Fuchs:2003qc,Geng:2013xn}. Such proposals have been explored mainly by two groups, the Bochum group~\cite{Epelbaum:2012ua,Baru:2019ndr,Ren:2019qow} and the BUAA group~\cite{Ren:2016jna,Li:2016mln,Li:2018tbt,Xiao:2018jot,Song:2018qqm,Xiao:2020ozd,Bai:2020yml,wang:2020myr,Liu:2020uxi}.

In addition to being covariant, the newly proposed scheme has another interesting feature, i.e., some low energy constants contribute to different partial wave phase shifts. This is in contrast with the conventional chiral nuclear force, where a single LEC only contributes to a particular partial wave. The particular feature raises two intriguing questions: 1) whether such correlations are physical and supported by data; 2) how they can be utilized to better understand lattice QCD simulations for unphysical pion masses or baryon-baryon interactions where there are little or no constraints. It is the main purpose of the present work to address these two questions.

To address the first question, we turn to the physical nucleon-nucleon phase shifts, which have been measured very precisely. We show that it is possible to fix the unknown LECs by fitting to either $\delta_{3S1}$ and then predict $\delta_{3D1}$ and $\varepsilon_1$ up to  $T_\mathrm{lab.}\approx 100$ MeV. Similarly, we can also fit to $\delta_{3D1}$ and $\varepsilon_1$ and predict $\delta_{3S1}$. These demonstrate that the correlation contained in the covariant chiral EFT is indeed correct.

To address the second question, we apply the same approach to study the lattice QCD nucleon-nucleon phase shifts for $m_\pi=469$ MeV~\cite{Inoue:2011ai}. Although the covariant approach works more or less fine, but the largest $T_\mathrm{lab.}$ for which a reasonable description can be achieved is reduced to about 10 MeV. This seems to indicate that  for some unknown reasons, the application region of the covariant chiral EFT for unphysical pion masses is much smaller.

To check whether this is particular to the covariant chiral nuclear force, we turn to the next-to-leading order HB chiral effective field theory, because in this case, only one more LEC is needed~\cite{Epelbaum:2014efa}, while in the covariant case, more LECs are needed and therefore they cannot be properly determined by the limited lattice QCD data~\cite{Inoue:2011ai}.  We show that the NLO HB chiral nuclear force can describe the lattice QCD data up to $T_\mathrm{lab.}\approx40$ MeV, but the description of the physical nucleon-nucleon phase shifts remain similar to the leading order covariant ChEFT. We can trace the improvement to the decoupling of the different partial waves and to the fact that one independent LEC appears for the mixing angle.

This work is organized as follows. In Sec.II, we briefly introduce the leading order covariant ChEFT and the next-to-leading order HB ChEFT, with particular emphasis on the description of light quark mass dependence.  In Sec.III, we study the predictive power of the covariant ChEFT for both the physical nucleon-nucleon phase shifts and the lattice QCD nucleon-nucleon phase shifts. Then, we explore how the next-to-leading order HB chiral nuclear force can describe the lattice QCD data. A short summary and outlook is given in Sec.IV.

\section{Theoretical framework}
The ChEFT we employ in the present work is described in detail in Refs.~\cite{Ren:2016jna,Xiao:2018jot,wang:2020myr}. Here
we only provide a concise introduction and relevant formulation needed for studying the light quark mass dependence of the NN interaction.

At leading order, the covariant chiral potential in the pion-full version can be written as.
\begin{equation}\label{1}
  V_{LO}=V_{CT}+V_{OPE}
\end{equation}
where $V_{CT}$ represents contact contributions and $V_{OPE}$ denotes
one-pion exchanges. The contact contributions are described by
four covariant four-fermion contact terms without derivatives~\cite{Xiao:2018jot},  namely,
\begin{eqnarray}
  V_{CT} &=& C_{S}\left[\overline{u}\left(\pmb{p'},s'_{1}\right)u\left(\pmb{p},s_{1}\right)\right]\left[\overline{u}\left(-\pmb{p'},s'_{2}\right)u\left(-\pmb{p},s_{2}\right)\right]\nonumber \\
   &+&  C_{V}\left[\overline{u}\left(\pmb{p'},s'_{1}\right)\gamma_{\mu}u\left(\pmb{p},s_{1}\right)\right]\left[\overline{u}\left(-\pmb{p'},s'_{2}\right)\gamma^{\mu}u\left(-\pmb{p},s_{2}\right)\right] \\
   &+&  C_{AV}\left[\overline{u}\left(\pmb{p'},s'_{1}\right)\gamma_{5}\gamma_{\mu}u\left(\pmb{p},s_{1}\right)\right]\left[\overline{u}\left(-\pmb{p'},s'_{2}\right)\gamma_{5}\gamma^{\mu}u\left(-\pmb{p},s_{2}\right)\right] \nonumber \\
   &+&   C_{T}\left[\overline{u}\left(\pmb{p'},s'_{1}\right)\sigma_{\mu\nu}u\left(\pmb{p},s_{1}\right)\right]\left[\overline{u}\left(-\pmb{p'},s'_{2}\right)\sigma^{\mu\nu}u\left(-\pmb{p},s_{2}\right)\right],\nonumber
\end{eqnarray}
where $\pmb{p}$ and $\pmb{p'}$ are initial and final three momentum, $s_1(s'_1)$, $s_2(s'_2)$ are spin projections, $C_{S,V,AV,T}$ are LECs, and $u(\overline{u})$ are Dirac spinors,
\begin{equation}\label{3}
  u(\pmb{p},s)=N_{P}\begin{pmatrix} 1 \\ \frac{\pmb{\sigma \cdot p}}{E_{P}+M} \end{pmatrix}\chi_{S}, \qquad N_{P}=\sqrt{\frac{E_{P}+M}{2M}}
\end{equation}
with $\chi_{S}$ being the Pauli spinor and $E_{P}(M)$ being the nucleon energy (mass). The one-pion-exchange potential in momentum space reads
\begin{eqnarray}
   V_{OPE}(\pmb{p',p}) &=& -\left(g_{A}^2/4f_{\pi}^{2}\right)\times \left[\overline{u}\left(\pmb{p'},s_{1}'\right)\pmb{\tau_{1}}\gamma^{\mu}\gamma_{5}q_{\mu}u(\pmb{p},s_{1})\right]\nonumber\\
    &\times& \frac{ \overline{u}\left(\pmb{-p'},s_{2}'\right)\pmb{\tau_{2}}\gamma^{\mu}\gamma_{5}q_{\nu}u\left(\pmb{-p},s_{2}\right)}{\left(E_{p'}-E_{p}\right)^{2}-\left(\pmb{p'-p}\right)^{2}-m_{\pi}^{2}},
\end{eqnarray}
where $m_{\pi}$ is the pion mass,  $\tau_{1,2}$ are the isospin matrcies, $g_A = 1.26$, and $f_{\pi} = 92.4$ MeV. Note that the leading order covariant potentials already contain all the six spin operators needed to describe nucleon-nucleon scattering.

Up to the next-to-leading order, the conventional (non-relativistic) chiral potential contains three terms~\cite{Epelbaum:2014efa}, i.e.,
\begin{equation}
  V_{NLO}=V_{CT}+V_{OPE}+V_{TPE},
\end{equation}
where the contact term, the one-pion exchange term, and the two-pion exchange term read, respectively,
\begin{equation}
  V_{CT}=V^{(0)}_{CT}+V^{(2)}_{CT},
\end{equation}
\begin{equation}
  V^{(0)}_{CT}=C_{S}+C_{T} \overrightarrow{\sigma_{1}} \cdot \overrightarrow{\sigma_{2}},
\end{equation}
\begin{eqnarray}
   V^{(2)}_{CT} &=& C_{1}q^2+C_{2}k^2+\left(C_{3}q^2+C_{4}k^2\right)\left(\overrightarrow{\sigma_{1}} \cdot \overrightarrow{\sigma_{2}}\right)\nonumber\\
   &+& \frac{i}{2}C_{5}\left(\overrightarrow{\sigma_{1}} + \overrightarrow{\sigma_{2}}\right)\cdot\left(\overrightarrow{k}\times\overrightarrow{q}\right)+C_{6}\left( \overrightarrow{q}\cdot \overrightarrow{\sigma_{1}}\right)\left( \overrightarrow{q}\cdot \overrightarrow{\sigma_{2}}\right)\nonumber \\
   &+&C_{7}\left ( \overrightarrow{k}\cdot \overrightarrow{\sigma_{1}}\right)\left( \overrightarrow{k}\cdot \overrightarrow{\sigma_{2}}\right),
\end{eqnarray}
\begin{equation}
  V_{OPE}=-\frac{g^2_{A}}{4F^2_{\pi}}\frac{\overrightarrow{\sigma_1}\cdot\overrightarrow{q}\overrightarrow{\sigma_2}\cdot\overrightarrow{q}}{q^2+m^2_{\pi}},
\end{equation}
\begin{eqnarray}
  V_{TPE} &=& V_{C}+\pmb{\tau_1\cdot\tau_2}W_C\nonumber \\
    &+& \left[V_S+\pmb{\tau_1\cdot\tau_2}W_S\right] \overrightarrow{\sigma_1}\cdot\overrightarrow{\sigma_2}\nonumber\\
   &+& \left[V_T+\pmb{\tau_1\cdot\tau_2}W_T\right]\overrightarrow{\sigma_1}\cdot\overrightarrow{q}\overrightarrow{\sigma_2}\cdot\overrightarrow{q} \nonumber\\
   &+&  \left[V_{LS}+\pmb{\tau_1\cdot\tau_2}W_{LS}\right] i\left(\overrightarrow{\sigma_1}+\overrightarrow{\sigma_2}\right)\cdot\left(\overrightarrow{q}\times\overrightarrow{k}\right),
\end{eqnarray}
where $\overrightarrow{k}=\left(\overrightarrow{p}+\overrightarrow{p}'\right)/2$, $\pmb{\tau_i}$ denote the isospin Pauli matrices associated with the nucleon i, while $V_{C,S,T,LS}$ and
$W_{C,S,T,LS}$ are scalar functions which depend on the nucleon momenta. The order-$Q^2$ contributions take the following form
\begin{eqnarray}
  W^{(2)}_{C} &=& -\frac{L(q)}{384\pi^2F^4_{\pi}}\big[4m^2_{\pi}\left(5g^4_A-4g^2_A-1\right)\nonumber\\
  &+& q^2\left(23g^4_A-10g^2_A-1\right)+\frac{48g^4_Am^4_{\pi}}{\left(4m^2_{\pi}+q^2\right)}\big],
\end{eqnarray}

\begin{equation}
  V_T^{(2)}=-\frac{1}{q^2}V^{(2)}_S=-\frac{3g^4_A}{64\pi^2F^4_{\pi}}L(q),
\end{equation}

\begin{equation}
  V^{(2)}_C= V^{(2)}_{LS}= W^{(2)}_S=W^{(2)}_T=W^{(2)}_{LS}=0.
\end{equation}
The loop function $L(q)$ is defined in dimensional regularization (DR) as
\begin{equation}
  L(q)=\frac{\sqrt{4m^2_{\pi}+q^2}}{q}In\frac{\sqrt{4m^2_{\pi}+q^2}+q}{2m_{\pi}}.
\end{equation}

The contact potentials can be projected into different partial waves in the $|LSJ\rangle$ basis. In the present work, we are only
interested in the $^3S_1$-$^3D_1$ coupled channel. In the covaraint ChEFT, the corresponding partial wave potentials read
\begin{eqnarray}
 V_{3S1}  &=& \frac{\xi_{N}}{9}\left[C_{3S1}\left(9+R^{2}_{\pmb{p}}R^{2}_{\pmb{p'}}\right)+\hat{C}_{3S1}\left(R^{2}_{\pmb{p}}+R^{2}_{\pmb{p'}}\right)\right],
\end{eqnarray}
\begin{eqnarray}
  V_{3D1} &=& \frac{8\xi_{N}}{9}C_{3S1}R^{2}_{\pmb{p}}R^{2}_{\pmb{p'}},
\end{eqnarray}
\begin{eqnarray}
 V_{3S1-3D1}  &=& \frac{2\sqrt{2}\xi_{N}}{9}\left[C_{3S1}R^{2}_{\pmb{p}}R^{2}_{\pmb{p'}}+\hat{C}_{3S1}R^{2}_{\pmb{p}}\right],
\end{eqnarray}
\begin{eqnarray}
 V_{3D1-3S1}  &=& \frac{2\sqrt{2}\xi_{N}}{9}\left[C_{3S1}R^{2}_{\pmb{p}}R^{2}_{\pmb{p'}}+\hat{C}_{3S1}R^{2}_{\pmb{p'}}\right],
\end{eqnarray}
where $ \xi_{N}=4 \pi N_{p}^{2}N_{p}^{'2}$, $R_{\pmb{p}}=|\pmb{p}|/(E_{p}+M)$, and $R_{\pmb{p'}}=|\pmb{p'}|/(E_{p'}+M)$.

In the conventional ChEPT up to NLO, the $^3S_1$-$^3D_1$ coupled channel potentials read
\begin{eqnarray}
 V_{3S1}&=&\widetilde{C}_{3S1}+ C_{3S1}\left(\pmb{p^2}+\pmb{p'^2}\right),
\end{eqnarray}

\begin{eqnarray}
 V_{3S1-3D1}  &=& C_{3DS1} \pmb{p^2},
\end{eqnarray}

\begin{eqnarray}
 V_{3D1-3S1}  &=& C_{3DS1} \pmb{p'^2},
\end{eqnarray}
where $C_{3S1}$, $\widetilde{C}_{3S1}$, and $C_{3DS1}$
are linear combinations of $C_S$, $C_T$, and $C_{1,2,3,4,5,6,7}$, which are not relevant to our purpose here. We note that there are no LECs contributing to $V_{3D1}$ up to NLO. As a result, only the OPE and TPE contribute to $\delta_{3D1}$.

To study pion mass dependent lattice QCD data,  in the covariant ChEFT, following Ref.~\cite{Bai:2020yml},    we  add a pion-mass dependent term in the  $^3S_1$ potential but keep the momentum expansion at the leading order, similar to, e.g., Ref.~\cite{Haidenbauer:2017dua}. The resulting LECs then read, explicitly,
\begin{eqnarray}
  C_{3S1}  &\rightarrow & C_{3S1}+C_{3S1}^{\pi} m_{\pi}^{2},\\
  \hat{C}_{3S1}  &\rightarrow  & \hat{C}_{3S1}+\hat{C}_{3S1}^{\pi} m_{\pi}^{2}.
\end{eqnarray}

Similarly, in the NLO conventional ChEFT, one also adds one pion mass dependent term
\begin{eqnarray}
  \widetilde{C}_{3S1} &\rightarrow & \widetilde{C}_{3S1}+C_{3S1}^{\pi} m_{\pi}^{2},\\
  C_{3S1}&\rightarrow&C_{3S1},\\
  C_{3DS1}&\rightarrow&C_{3DS1}.
\end{eqnarray}
 Comparing the leading order covariant ChEFT and the  next-to-leading order HB ChEFT, it is clear that one more LEC is needed in the latter, i.e., $C_{3DS1}$. In the covariant ChEFT, however, there are only two LECs for the $^3S_1$-$^3D_1$ coupled channel.  As a result, one can anticipate that the mixing angle can be better described in the latter.

\section{Results and discussions}
In this section, we first demonstrate the predictive power of the covariant ChEFT, in particular, its ability to connect $\delta_{3S1}$, $\delta_{3D1}$, and $\varepsilon_1$, both for physical $M_\pi$ and $M_\pi=469$ MeV. Then to confirm the observation that the lattice QCD data can only be described up to $T_\mathrm{lab.}\approx 10$ MeV, we turn to the next to leading order HBChEFT, where there is no correlation between the LECs for different partial waves, and therefore we expect that the lattice QCD simulations can be better described~\footnote{It has been shown that in Ref.~\cite{Ren:2016jna} the LO covariant ChEFT can describe the $J=0,1$ nucleon-nucleon partial waves at a level similar to the NLO HB ChPT, although with a smaller number of LECs. }. We show that even in such a case the lattice QCD data can only be qualitatively described  up to $T_\mathrm{lab.}\approx40$ MeV. The energy range for which the NLO HBChPT works is still much reduced compared to the physical nucleon-nucleon phase shifts.

\renewcommand{\arraystretch}{1.5}

\begin{figure*}[htpb]
\includegraphics[scale=0.29]{1f.pdf}
\includegraphics[scale=0.29]{2f.pdf}
\includegraphics[scale=0.28]{3f.pdf}
\caption{ Fitting to $\delta_{3S1}$ and predicting $\delta_{3D1}$ and $\varepsilon_1$ by the covariant ChEFT. Solid circles denote the data fitted and empty ones denote those not fitted. \label{fig:3s1to3d1}}
\end{figure*}

\begin{table}[htpb]
\centering
\caption{Values of the LECs of the best fit to the Nijmegen $\delta_{3S1}$.The LECs in the  first column are obtained by fitting to four experimental data up to $T_\mathrm{max}=100$ MeV and those in the second column by fitting to five experimental data up to $T_\mathrm{max}=150$ MeV.\label{tab:3s1to3d1}}
\begin{tabular}{c|  c |c}
 \hline\hline
  $C_{3S1}$  [$10^{-2}$ MeV$^{-2}$] & $\hat{C}_{3S1}$   [$10^{-2}$ MeV$^{-2}$]  &  $\Lambda$ [MeV]\\ [0.3ex]
 \hline
 $-0.141$ & $-0.686$ & 905 \\[0.3ex]
 \hline
 $-0.126$ & $-0.603$ & 920 \\[0.3ex]
 \hline\hline
\end{tabular}
\end{table}

\subsection{Covariant ChEFT}
First, we fix the  two  LECs $C_{3S1}$ and $\hat{C}_{3S1}$ and the cutoff by fitting to $\delta_{3S1}$ for $T_\mathrm{lab.}<100$ and 150 MeV, respectively. The resulting LECs and phase shifts are shown in Table~\ref{tab:3s1to3d1} and the left panel of Fig.~\ref{fig:3s1to3d1}. With the same LECs and cutoff, we predict $\delta_{3D1}$ and $\varepsilon_1$, which are shown in the middle and right panels of Fig.~\ref{fig:3s1to3d1}. Two things are noteworthy. First, the $^3S_1$ phase shifts can be described quite well, regardless of the $T_\mathrm{max}$ fitted. Second, the predicted
$\delta_{3D1}$ and $\varepsilon_1$ are in reasonable agreement with the data, particularly the former. The description of the mixing angle is not very satisfactory, but the qualitative feature is produced.

Next we fix the two LECs  $C_{3S1}$ and $\hat{C}_{3S1}$ and the cutoff by fitting to $\delta_{3D1}$ and $\varepsilon_1$. Then we predict $\delta_{3S1}$. The corresponding LECs and phase shifts are shown in Table \ref{tab:3d1to3s1} and Fig.~\ref{fig:3d1to3s1}. We note that  the predicted $\delta_{3S1}$ are also in reasonable agreement with data.
Contrary to the previous case, now one can better describe both $\delta_{3D1}$ and $\varepsilon_1$, particularly the latter.
This is because the mixing angle is much small in magnitude compared to $\delta_{3S1}$, while in our fitting, we assigned the same weights to $\delta_{3S1}$, $\delta_{3D1}$, and $\varepsilon_1$.

We stress that in both cases the correlation implicitly contained in the contact terms of the covariant framework is supported. Similar correlations have been noticed in the strangeness $S=-3/-4$ baryon-baryon~\cite{Liu:2020uxi} as well as the $\Lambda_c N$ systems~\cite{Song:2020isu,Song:2021war}, although for unphysical pion masses.

\begin{figure*}[htpb]
\includegraphics[scale=0.293]{11f.pdf}
\includegraphics[scale=0.29]{22f.pdf}
\includegraphics[scale=0.29]{33f.pdf}
\caption{ Fitting to $\delta_{3D1}$ and $\varepsilon_1$ and predicting $\delta_{3S1}$ by the covariant ChEFT. Solid circles denote the data fitted and empty ones denote those not fitted. \label{fig:3d1to3s1}}
\end{figure*}

\begin{table}[htpb]
\centering
\caption{Same as Table~\ref{tab:3d1to3s1}, but fitted to $\delta_{3D1}$ and $\varepsilon_1$.\label{tab:3d1to3s1}}
\begin{tabular}{c|  c |c}
 \hline\hline
  $C_{3S1}$  [$10^{-4}$ MeV$^{-2}$] & $\hat{C}_{3S1}$   [$10^{-2}$ MeV$^{-2}$]  &  $\Lambda$ [MeV]\\ [0.3ex]
 \hline
 0.268 & 0.123 & 1135 \\[0.3ex]
 \hline
 0.259 & 0.122 & 1130 \\[0.3ex]
 \hline\hline
\end{tabular}
\end{table}

Now we turn to the lattice QCD nucleon-nucleon phase shifts obtained for $M_\pi=469$ MeV, which have been studied in Refs.~\cite{Bai:2020yml,Hu:2020djy}. First, we notice that the energy range for which a simultaneous description of $\delta_{3S1}$, $\delta_{3D1}$, and $\varepsilon_1$ can be achieved is much reduced~\cite{Bai:2020yml}. Therefore, in the following, we check the largest kinetic energy up to which a reasonable description can be achieved. We again explore the two strategies adopted for studying the physical nucleon-nucleon phase shits.

First, we fit to $\delta_{3S1}$ up to $T_\mathrm{lab.}\approx6$ and 12 MeV, and predict $\delta_{3D1}$ and $\varepsilon_1$. The corresponding LECs~\footnote{The optimal cutoff values are a bit small compared to the preferred value of $0.5\sim1$ GeV, indicating that one should be less demanding in reproducing the lattice QCD data. On the other hand, as we can see in Table IV,  the optimal cutoff values obtained by fitting to $\delta_{3D1}$ and $\varepsilon_1$ look more natural.} are listed in Table~\ref{tab:lat3s1to3d1} and phase shifts given in Fig.~\ref{fig:lat3s1to3d1}. Clearly, although the central values of $\delta_{3S1}$ are fitted quite well, the predicted $\delta_{3D1}$ and $\varepsilon_1$ are a bit off. Two features of the lattice QCD data are responsible for this failure. First, the lattice QCD $\delta_{3D1}$ in the energy range studied are practically zero. Second, the mixing angle $\varepsilon_1$ are very small as well. As a matter of fact, as noted in Ref.~\cite{Bai:2020yml}, the mixing angle $\varepsilon_1$ for $M_\pi=469$ MeV is similar to that for $M_\pi=138$ MeV, exhibiting little light quark mass dependence.

\begin{figure*}[htpb]
\includegraphics[scale=0.288]{111f.pdf}
\includegraphics[scale=0.295]{222f.pdf}
\includegraphics[scale=0.288]{333f.pdf}
\caption{ Fitting to $\delta_{3S1}$ and predicting $\delta_{3D1}$ and $\varepsilon_1$ by the covariant ChEFT. Solid circles denote the data fitted and empty ones denote those not fitted.\label{fig:lat3s1to3d1}}
\end{figure*}

\begin{table}[htpb]
\centering
\caption{Values of the LECs of the best fit to the lattice QCD simulations of
$\delta_{3S1}$. The LECs in the first column are obtained by fitting to five lattice QCD data and those in the second column by fitting to six lattice QCD data.\label{tab:lat3s1to3d1}}
\begin{tabular}{c|  c |c}
 \hline\hline
  $C_{3S1}$  [$10^{-4}$ MeV$^{-2}$] & $\hat{C}_{3S1}$   [$10^{-1}$ MeV$^{-2}$]  &  $\Lambda$ [MeV]\\ [0.3ex]
 \hline
 0.669 & $-0.344$ & 243 \\[0.3ex]
 \hline
 0.544 & $-0.559$ & 231 \\[0.3ex]
 \hline\hline
\end{tabular}
\end{table}

Considering that the uncertainties of $\delta_{3S1}$ are much larger than those of $\delta_{3D1}$ and $\varepsilon_1$, which are almost negligible, we expect that the situation could be improved for strategy 2. Indeed, as one can see from Fig.~\ref{fig:lat3d1to3s1} with the corresponding LECs given in Table~\ref{tab:lat3d1to3s1}, now $\delta_{3D1}$ and $\varepsilon_1$ are described reasonably well. However, $\delta_{3S1}$ are only in reasonable agreement with the lattice QCD data within the relatively large uncertainties of the lattice QCD simulations.

A particularly puzzling observation is that from $M_\pi=138$ MeV to $M_\pi=469$ MeV, the largest energy  to which the leading order covariant ChEFT is applicable is reduced from about $T_\mathrm{lab.}\approx100$ MeV to about 10 MeV.  This needs to be better understood in the future.

\begin{figure*}[htpb]
\includegraphics[scale=0.295]{1111f.pdf}
\includegraphics[scale=0.288]{2222f.pdf}
\includegraphics[scale=0.288]{3333f.pdf}
\caption{ Fitting to $\delta_{3D1}$ and $\varepsilon_1$ and predicting $\delta_{3S1}$ by the covariant ChEFT. Solid circles denote the data fitted and empty ones denote those not fitted.\label{fig:lat3d1to3s1}}
\end{figure*}

\begin{table}[htpb]
\centering
\caption{Same as Table~\ref{tab:lat3s1to3d1}, but fitted to $\delta_{3D1}$ and $\varepsilon_1$.\label{tab:lat3d1to3s1}}
\begin{tabular}{c|  c |c}
 \hline\hline
  $C_{3S1}$  [$10^{-4}$ MeV$^{-2}$] & $\hat{C}_{3S1}$   [$10^{-2}$ MeV$^{-2}$]  &  $\Lambda$ [MeV]\\ [0.3ex]
 \hline
 0.207 & 0.780 & 548 \\[0.3ex]
 \hline
 0.260 & 0.825 & 513\\[0.3ex]
 \hline\hline
\end{tabular}
\end{table}

\subsection{Conventional ChEFT}

In this section, we check how the physical and the lattice QCD nucleon-nucleon phase shifts can be described by the NLO conventional ChEFT. As noted earlier, in the conventional ChEFT, there are two LECs for the $^3S_1$ partial wave, one LEC for the mixing angle, and no LEC for the $^3D_1$ partial wave. As a result, the contribution to $\delta_{3D1}$ solely comes  from the one-pion and two-pion exchanges. In this framework, the correlation between $\delta_{3S1}$, $\delta_{3D1}$, and $\varepsilon_1$ is induced by the one- and two-pion exchanges. We have checked that the two-pion exchange contributions are negligible.  It should be noted that because the decoupling of LECs, we can only fit to $\delta_{3S1}$, $\delta_{3D1}$, and $\varepsilon_1$ simultaneously.

In Fig.~\ref{fig:conphys}, we compare the NLO conventional ChEFT results with the corresponding empirical data, with the LECs given in Table ~\ref{tab:conphys}. As can be clearly seen, up to $T_\mathrm{lab.}\approx100$ MeV, the empirical data can be described quite well, as already noted in Ref.~\cite{Ren:2016jna}, which is similar to the covariant ChEFT~\cite{Ren:2016jna,wang:2020myr}. In addition, the extra LEC for the mixing angle also allows for a better quantitative description of this observable.

Now we turn to the lattice QCD data obtained for $M_\pi=469$ MeV. The corresponding results are shown in Fig.~\ref{fig:conlat} and Table~\ref{tab:conlat}. One can see that the lattice QCD data can be described reasonably well only up to $T_\mathrm{lab.}\approx40$ MeV. Compared to the physical case, the energy range for which the NLO ChEFT works is also reduced by a factor of  about three.  Therefore, we can conclude that as the pion mass increases, the validity region of ChEFT in terms of kinetic energies for the nucleon-nucleon scattering is reduced.

\begin{figure*}[htpb]
\includegraphics[scale=0.29]{1.pdf}
\includegraphics[scale=0.29]{2.pdf}
\includegraphics[scale=0.29]{3.pdf}
\caption{ Simultaneous description of the  PWA93   phase shifts by the NLO conventional ChEFT.\label{fig:conphys}}
\end{figure*}

\begin{table}[htpb]
\centering
\caption{Values of the LECs of the best fit to the Nijmegen $\delta_{3S1}$, $\delta_{3D1}$, and $\varepsilon_1$ by  the NLO conventional ChEFT.\label{tab:conphys}}
\begin{tabular}{c| c| c |c}
 \hline\hline
  $C_{3S1P}$  [$10^{-1}$ MeV$^{-2}$] & $C_{3S1}$   [$10^{-7}$ MeV$^{-4}$]  & $C_{3DS1}$  [$10^{-8}$ MeV$^{-4}$] &  $\Lambda$ [MeV]\\ [0.3ex]
 \hline
 0.210 & 0.242 & $-0.557$  & 662 \\[0.3ex]
 \hline
 0.251 & 0.254 & $-0.448$  & 670 \\[0.3ex]
 \hline\hline
\end{tabular}
\end{table}

\begin{figure*}[htpb]
\includegraphics[scale=0.29]{4.pdf}
\includegraphics[scale=0.29]{5.pdf}
\includegraphics[scale=0.29]{6.pdf}
\caption{ Simultaneous description of the HALQCD469  phase shifts by the NLO conventional ChEFT.\label{fig:conlat}}
\end{figure*}

\begin{table}[htpb]
\centering
\caption{Same as Table~\ref{tab:conphys}, but fitted to the  lattice QCD simulations by the NLO conventional ChEFT.\label{tab:conlat}}
\begin{tabular}{c| c| c |c}
 \hline\hline
  $C_{3S1P}$  [$10^{-1}$ MeV$^{-2}$] & $C_{3S1}$   [$10^{-7}$ MeV$^{-4}$]  & $C_{3DS1}$  [$10^{-7}$ MeV$^{-4}$] &  $\Lambda$ [MeV]\\ [0.3ex]
 \hline
 0.115 & 0.235 & $-0.191$  & 555 \\[0.3ex]
  \hline
 0.147 & 0.315 & $-0.195$  & 540 \\[0.3ex]
 \hline\hline
\end{tabular}
\end{table}

\section{Summary and outlook}
In this work, we explored the predictive power of the covariant chiral effective field theory for the nucleon-nucleon scattering in the $^3S_1$-$^3D_1$ coupled channel, because the same low energy constants are responsible for the three observables $\delta_{3S1}$, $\delta_{3D1}$, and $\varepsilon_1$. We checked the correlation using the physical $NN$ phase shifts and found that it works pretty well up to $T_\mathrm{lab.}\approx100$ MeV. However, once applied to study the lattice QCD nucleon-nucleon phase shifts obtained for $M_\pi=469$ MeV, we found that the covariant ChEFT only works up to $T_\mathrm{lab.}\approx 10$ MeV at a quantitative level. Nevertheless, at a qualtitative level, the correlation predicted by the covariant ChEFT seems to hold for a larger kinetic energy.

To confirm the observed reduction of the validity region of ChEFT, we performed studies using the next-to-leading order conventional ChEFT. We found that with an additional independent LEC for the mixing angle $\varepsilon_1$, the lattice QCD data can be reasonably described up to $T_\mathrm{lab.}\approx40$ MeV, which is larger than that of the leading order covariant ChEFT, but still roughly a factor of three smaller than the corresponding value achieved in describing the physical nucleon-nucleon phase shifts.

The correlation exhibited by the covariant ChEFT can be of relevance in understanding lattice QCD simulations of other baryon-baryon systems~\cite{Song:2018qqm,Li:2018tbt,Liu:2020uxi,Song:2020isu}, such as the strangeness $-3/-4$ systems~\cite{Liu:2020uxi}. However, as the present study demonstrated, the validity region of ChEFT for unphysical pion masses might be considerably smaller than that for physical pion masses. This should be kept in mind when studying lattice QCD simulations obtained for large unphysical pion masses. In a recent work, it was shown that for the $\Lambda_c N$ system the lattice QCD data can only be described at a quantitative level up to $E_\mathrm{c.m.}\approx 5$ MeV~\cite{Song:2021war}.

It should be noted that our current study is tied to the HALQCD simulations of the nucleon-nucleon interaction. It will be of great value if other lattice QCD collaborations can perform more simulations, which may help better understand the pion mass dependence of the nucleon-nucleon interaction, and baryon-baryon interactions in general. In recent years, a large amount of the-called exotic hadrons, which cannot easily fit into the conventional quark model, have been discovered at word-wide high energy facilities, such as LHC, KEKB, and BEPC. Because many of them are located close to the thresholds of two  conventional hadrons, they are believed to be hadronic molecules. Clearly, in such a picture, lattice QCD simulations of hadron-hadron interactions are very important to understand these exotic hadrons. Our current work could also of relevance to such studies.

\section{Acknowledgements}

This work is partly supported by the National Natural Science Foundation of China under Grants No.11735003, No.11975041,  and No.11961141004, and the fundamental Research Funds for the Central Universities.