\section{Introduction}
\noindent {\it Introduction.$-$}
The interaction between vector mesons and the nucleon ($N$) is one of the most fundamental quantities to study the meson properties in nuclear matter and the modification of QCD condensates at finite baryon density ~\cite{Hayano-Hatsuda2010}.
Among others, the $\phi$-meson attracts particular interest, since the interactions of the $s\bar{s}$ pair with the nucleon and nucleus are not well understood at low energies. 
 
Experimentally, the production and absorption of the $\phi$-meson  in nuclear matter has been actively studied through photon and proton($p$)-induced reactions with various nuclear targets
(Spring8-LEPS~\cite{Ishikawa2005}, KEK-PS E325~\cite{KEK-PS2007}, JLab-CLAS~\cite{CLAS2010}, and COSY-ANKE~\cite{Polyanskiy2011}) as summarized in Ref.~\cite{Tolos2020}.
Recently, the ALICE Collaboration at LHC presented a new result on the $p$-$\phi$ interaction  through the femtoscopic analysis of the $p$-$\phi$ pair produced in the $pp$ reaction:
They reported a spin-averaged scattering length $a_{p\phi} = -0.85(34)_{\rm stat.}(14)_{\rm syst.}$ fm~\cite{ALICE2021}, an order of magnitude larger than that obtained from the photo-production data by the CLAS Collaboration at JLab~\cite{Dey:2014tfa} combined with the vector meson dominance \cite{Strakovsky2020}.
            
Theoretically, the properties of the $\phi$-meson in nuclear matter and its relation to the in-medium $\bar{s}s$ condensate have been extensively discussed~\cite{Hatsuda:1991ez,Asakawa:1994tp,Gubler:2016itj,Gubler:2018ctz}.  
The microscopic origin of the force between $\phi$ and $N$ is also an important open problem.
In particular, the two-pion exchange (TPE) as the major long-range contribution between a color-dipole and the nucleon~\cite{TarrusCastella:2018php} is of strong relevance to this problem and will shed new light on the important role of pion dynamics in hadron interactions.
  

 Under these experimental and theoretical circumstances, it is most desirable to carry out realistic lattice QCD simulations of the $N$-$\phi$ interaction.
In this Letter, we report a first simulation of the $N$-$\phi$ system in a large lattice volume $\simeq$ (8.1 fm)$^3$  with light dynamical quarks near the physical point.
We focus on the highest-spin $N$-$\phi$ system, the $^4S_{3/2}$ channel,  with the notation $^{2s+1}L_J$ ($s$ = total spin, $L$ = orbital angular momentum, $J$ = total angular momentum). This is because  its coupling to two-body open channels $\Lambda K(^2D_{3/2})$ and $\Sigma K(^2D_{3/2})$ are kinematically suppressed at low energies  due to their $D$-wave nature, similar to the situation for the highest-spin   $N$-$\Omega$ system (the $^5S_{2}$ channel)  \cite{Iritani2019PLB}. 
Moreover, decay processes with three or more particles in the final state such as $\Sigma \pi K$, $\Lambda \pi K$, and $\Lambda \pi\pi K$ are expected to be suppressed due to the small phase space.
Indeed, we will show below that visible effects from such decays cannot be seen from our lattice data in the $^4S_{3/2}$  channel.
By using the HAL QCD method~\cite{Ishii2007, Ishii2012, Aoki2020}, which converts the spacetime correlation of hadrons to the physical observables,  
we obtain the $N$-$\phi$ scattering phase shift, the $S$-wave scattering length $a_{0}^{(3/2)}$ and the effective range $r_{\rm eff}^{(3/2)}$ in the $^4S_{3/2}$ channel. (For recent applications of the HAL QCD method to the baryon-baryon interactions, see Refs.~\cite{Gongyo2018, Iritani2019PLB, Sasaki2020, Lyu2021}.)

 

\noindent {\it HAL QCD Method.$-$}
Thanks to the Haag-Nishijima-Zimmermann reduction formula for composite particles~\cite{Zimmermann1987},
the interaction between hadrons can be extracted as the energy-independent non-local potential $U(\bm r,\bm r')$ ~\cite{Ishii2007}.
Such a potential can be defined from the normalized four-point function,
 \begin{eqnarray}\label{eq-R}
    R(\bm{r},t)&=&
   \frac{\sum_{\bm{x}}\langle 0| N(\bm{r}+\bm{x},t)\phi(\bm{x},t)\overline{\mathcal{J}}(0)|0\rangle} 
    {\sqrt{Z_{\phi}Z_N}  e^{-(m_{N}+m_{\phi})t}} \nonumber \\
    &=&\sum_n a_n \psi_n(\bm{r})e^{-(\Delta E_n)t}+O(e^{-(\Delta E^*)t}).
\end{eqnarray}
Here the energy eigenvalue and the energy shift of each elastic scattering state are given by $E_n=\sqrt{m_N^2+\bm{k}_n^2}+\sqrt{m_\phi^2+\bm{k}_n^2}$ 
and   $\Delta E_n=E_n-(m_N+m_\phi)$, respectively, with $\bm{k}_n$ being the relative momentum  in the center of mass frame.  
The equal-time Nambu-Bethe-Salpeter wave function for each eigenstate is denoted by $\psi_n(\bm r)$.  
The wave function renormalization constant for $\phi$ ($N$) is  given by $Z_{\phi}$ ($Z_N$).
 The contributions from the inelastic scattering states are exponentially suppressed as $O(e^{-(\Delta E^*)t})$ with $\Delta E^*$ being the energy of the inelastic threshold relative to $m_N+m_{\phi}$. 
The overlapping factor between the $N$-$\phi$   source operator  $\overline{\mathcal{J}}(0)$
 and the $n$-th eigenstate $|n\rangle$ is given by  $a_n=\langle{n}|\overline{\mathcal{J}}(0)|0\rangle$.
The wall-type  operator  with the Coulomb gauge fixing is used for the source, while the local operators are adopted for the sink,
\begin{equation}\label{eq-op}
 \begin{split}
  & N_\alpha(x)=\epsilon_{abc}\left[{{u^a}^T}(x)C\gamma_5d^b(x)\right]u^c_\alpha(x), \\
  & \phi_i(x)=\delta_{ab}\bar{s}^a(x)\gamma_is^b(x).
 \end{split}
\end{equation}
 Here $a$, $b$, and $c$ are the color indices, $i$ is the vector index, $\alpha$ is the Dirac index restricted to the upper two components, $C=\gamma_4\gamma_2$ is the charge conjugation, and $T$ denotes the transposition for the Dirac index.
Misner's method of approximate partial wave decomposition on a cubic grid, originally developed in the field of numerical relativity~\cite{Misner:1999ab}, is used to make the $S$-wave projection on the lattice~\cite{Miyamoto2020}.
The normalized four-point functions $R(\bm{r},t)$ with the spin projection onto the states with $J=3/2$ are used to extract the potential.

It has been shown that $R(\bm{r},t)$ satisfies the integro-differential equation ~\cite{Ishii2012},
\begin{equation}\label{eq-TDHAL}
 \begin{split}
  &\left[ \frac{1+3\delta^2}{8\mu}\frac{\partial^2}{\partial t^2}-\frac{\partial}{\partial t}-H_0 +O(\delta^2\partial_t^3) \right] R(\bm{r},t) \\
 &= \int d\bm{r}'U(\bm{r},\bm{r}')R(\bm{r}',t),
 \end{split}
\end{equation}
with $H_0=-\nabla^2/(2\mu)$, $\mu=m_N m_\phi/(m_N+m_\phi)$ being the reduced mass,  $\delta=(m_N-m_\phi)/(m_N+m_\phi)$ being the
 mass asymmetry.  We have neglected the term of $O(\delta^2\partial_t^3)$ in the above, since it is found to be consistent with zero within the statistical error in our simulation.  In practical calculations, a derivative expansion of the non-local potential is employed; $U(\bm{r},\bm{r}')=V(r)\delta(\bm{r}-\bm{r}')+ \sum\limits_{n=1}V_{n}(\bm{r})\nabla^{n}\delta(\bm{r}-\bm{r}')$.
 Then, the leading-order  central potential $V(r)$ reads,
\begin{equation}\label{eq-V}
 V(r)= R^{-1}(\bm{r},t)\left[ \frac{1+3\delta^2}{8\mu}\frac{\partial^2}{\partial t^2}-\frac{\partial}{\partial t}-H_0 \right] R(\bm{r},t).
\end{equation}
The truncation error of the derivative expansion is found to be small at low energies by using the method of finite volume spectra analysis~\cite{Iritani2019Jhep, Lyu2022} (Fig. S1 in the Supplemental Material). 


\noindent{\it Lattice setup.$-$}  
The ($2+1$)-flavor gauge configurations are generated with the Iwasaki gauge action at $\beta=1.82$ and the nonperturbatively $O(a)$-improved Wilson quark action with stout smearing at nearly physical quark masses~\cite{Ishikawa2016}.
The lattice spacing is $a\simeq0.0846$ fm ($a^{-1}\simeq2333$ MeV), and the lattice volume is $L^4=96^4$, corresponding to $La\simeq8.1$ fm.
The isospin-averaged masses of $\pi$, $K$, $\phi$, and $N$ 
 from the single-state fitting  in the interval $t/a=15$-25 for the mesons and 11-18 for the nucleon are 
 listed in Table ~\ref{tab-mass}  with statistical errors in the parentheses.
Corresponding experimental values~\cite{PDG2020} are shown in the third column.
The lattice results are about 6\%  larger than the experimental values for the pseudo-scalar mesons, while about 3\% (2\%) larger for $\phi$ ($N$).
 With our light quark masses slightly heavier than the physical masses, $\phi$ is below the $K\bar{K}$ threshold, 
which prohibits a fall-apart decay of  $\phi$.

\begin{table}[htbp]
\caption{ Isospin-averaged hadron masses with statistical errors obtained from (2+1)-flavor lattice QCD simulations together with the experimental values.
}
\begin{tabular}{ccc}
  \hline\hline
    Hadron &~~Lattice [MeV]  &~~Expt. [MeV] \\
  \hline
    $\pi$ &\ \ 146.4(4)  & \ \ 138.0\\
    $K$  & \ \ 524.7(2) & \ \ 495.6\\
    $\phi$ &\ \ 1048.0(4) &\ \ 1019.5\\
    $N$  &\ \ 954.0(2.9)  & \ \ 938.9 \\
  \hline\hline
\end{tabular} \label{tab-mass}
\end{table}

We use 200 configurations separated by 10 trajectories.
For each configuration, the forward and backward propagations are averaged, the hypercubic symmetry on the lattice (4 rotations) is utilized, and 80 measurements are performed by shifting the source position in a temporal direction in order to reduce the statistical fluctuation.
In total, the number of measurements is 128000.
With the  periodic boundary condition for all directions, the quark propagators are calculated by the domain-decomposed solver~\cite{Ishikawa2021},
and the hadron correlation functions are obtained by the unified contraction algorithm~\cite{Doi2013}.
The contraction diagram with the pair annihilation of $s\bar{s}$  is not considered, 
since it is known to be small, e.g. less than 0.4\% correction to the 
 $\phi$-meson mass \cite{Chakraborty:2017hry}.
The statistical errors are evaluated by the jackknife method with a bin size of 20 configurations throughout this Letter.
A comparison with a bin size of 40 configurations shows that the bin size dependence is small.
\begin{figure}[htbp]
  \centering
  \includegraphics[width=8cm]{Fig1}
  \caption{(Color online).
  The $N$-$\phi$ potential $V(r)$ in the ${^4S_{3/2}}$ channel as a function of separation $r$ at
 Euclidean time $t/a= 12$ (red squares), $13$ (green circles) and $14$ (blue triangles).
  } \label{Fig1}
\end{figure}

\noindent{\it Numerical results.$-$} 
The $N$-$\phi$ potential $V(r)$ in the $^4S_{3/2}$ channel defined in Eq.~(\ref{eq-V}) with the
lattice measurement of $R(\bm r, t)$ is shown in Fig.~\ref{Fig1} for Euclidean times, $t/a=12$, $13$, and $14$. (See Fig. S2 in the Supplemental Material for the $t$-dependence of $V(r)$ in a wider range of $t$.)
These Euclidean times are chosen such that they are large enough to suppress contaminations from excited states in the single-hadron correlator and simultaneously small enough to avoid exponentially
increasing statistical errors.
The potentials from different $t$ agree with each other within statistical errors, which indicates that the $N$-$\phi$ correlation function is dominated by the elastic scattering states in the $^4S_{3/2}$ channel; we do not find the effects of the two-body open channels  ($\Lambda K(^2D_{3/2})$ and $\Sigma K(^2D_{3/2})$) as well as the  three-body open channels including $N\phi \rightarrow \{\Sigma^* K, \Lambda(1405) K\} \rightarrow \{\Lambda \pi K, \Sigma \pi K$\}.
This is in sharp contrast to the  $^2S_{1/2}$ case where we found that the $N$-$\phi$ potential shows a clear $t$-dependence, as expected from the $S$-wave fall-apart decay into $\Lambda K(^2S_{1/2})$ and $\Sigma K(^2S_{1/2})$.

The potential $V(r)$ in the $^4S_{3/2}$ channel shown in Fig.~\ref{Fig1} is attractive for all distances and has a characteristic two-component structure, the 
attractive core at short distance and the attractive tail at long distance,  similar to the case of the $N\Omega(^5S_{2})$ potential~\cite{Iritani2019PLB}.  We note that the Pauli exclusion principle between quarks, which partially gives rise to the repulsive core in the $NN$ interaction~\cite{Inoue2012,Oka2000},  does not operate in the present case, since $N$ and $\phi$ have no common valence quarks.  
 
As has been discussed  for the interaction between color dipoles~\cite{Bhanot:1979vb,Fujii:1999xn,Brambilla:2015rqa},
  non-perturbative gluon exchange is  expected to appear in the form of  the TPE at long distance.
  The idea was generalized to the interaction between a color-dipole and the nucleon with the result, 
 $V(r \gg (2m_{\pi})^{-1} ) = -\alpha  \frac{\exp (-2m_\pi r)}{r^{2}}$, where $\alpha$ is proportional to $m_{\pi}^4$ \cite{TarrusCastella:2018php}.  
To check such a long distance behavior of $V(r)$,  we
 show in Fig.~\ref{Fig2} the spatial effective energy, $E_{\rm eff}(r)=-\frac{\ln[-V(r)r^{2}/\alpha]}{r}$,  
as a function of $r$ with $\alpha \simeq 91~{\rm MeV}\cdot{\rm fm}^{2}$ determined by fitting the lattice data of $V(r)$ at long distance.  
We find that $E_{\rm eff}(r)$ has a plateau at  $2m_\pi = 292.8$ MeV  for $r > 1.0$ fm, which indicates that the long-range part of the $N$-$\phi$ potential is indeed dominated by the TPE.

    
\begin{figure}[htbp]
  \centering
  \includegraphics[width=8cm]{Fig2}
  \caption{(Color online).
  The spatial effective energy $E_{\rm eff}(r)$ as a function of separation $r$ at
 Euclidean time $t/a= 12$ (red squares), $13$ (green circles) and $14$ (blue triangles).
The orange dashed line corresponds to $2m_\pi$ with lattice pion mass $m_\pi=146.4$ MeV.
  } \label{Fig2}
\end{figure}

In order to convert the potential to physical observables, we perform an uncorrelated fit of the lattice QCD potential 
by using two different functional forms,
\begin{eqnarray}\label{eq-fit}
 \! \! {\rm{A}}: V_{\rm{fit}}(r) &=& \sum_{i=1,2}a_i e^{-(r/b_i)^2} + a_3 m_\pi^4 f(r;b_3)  \frac{e^{-2m_\pi r}}{r^{2}},  \nonumber \\
\! \! {\rm{B}}:  V_{\rm{fit}}(r)&=&\sum_{i=1,2,3}a_i e^{-(r/b_i)^2}.
\end{eqnarray}

The fit A is motivated by the TPE tail at long distance with an overall strength proportional to $m_{\pi}^4$~\cite{TarrusCastella:2018php},
  while the fit B is  a purely  phenomenological Gaussian form  for comparison.
  In fit A, we consider two types of  form factors  commonly used in the $NN$ potentials: (i) The Nijmegen-type form factor \cite{Stoks:1994wp}, 
  $f_{\rm erfc}(r;b_3)= \left[  {\rm erfc}\left( \frac{m_{\pi}}{\Lambda}-\frac{\Lambda r}{2}  \right) - e^{2m_{\pi}r} {\rm erfc}\left(  \frac{m_{\pi}}{\Lambda}+\frac{\Lambda r}{2}   \right)  \right]^2/4$
with $\Lambda = 2/b_3$ and the complementary error function ${\rm erfc}(x)= \frac{2}{\sqrt{\pi}} \int_x^{\infty} e^{-z^2} dz$. 
This is motivated by the exponential-type regularization of the pion propagator in the momentum space, $1/(k^2+m_{\pi}^2)  \rightarrow e^{-(k/\Lambda)^2 } /(k^2+m_{\pi}^2 )$. 
 (ii) The Argonne-type  form factor~\cite{Wiringa1995}, $ f_{\rm exp}(r;b_3) =  \left(1-e^{-(r/b_3)^{2}}\right)^2$. 
We refer to fit A with (i) and (ii) as fit A$_{\rm erfc}$ and fit A$_{\rm exp}$, respectively.
The pion mass in fit A$_{\rm erfc, exp}$ is taken to be  $m_\pi=146.4$ MeV, and the fit range is  chosen as  $0 < r < 3.0$ fm.
We found that the fit A$_{\rm erfc, exp}$ and the fit B provide an equally good fit 
 ($\chi^2_{\rm d.o.f}=$0.3-0.4) and are stable against the choice of $t$.
 In Table~\ref{tab-prm}, we show the fit results for $t/a=14$, which are expected to have least contamination from the inelastic states. 
Changing the fit range of the potential to  $0.1<r<2.5$ fm does not affect the results within statistical errors.
Also we found that the simple fitting functions  such as the 
Yukawa  form $\sim  - \frac{\exp (-\mu r)}{r}$ ~\cite{Brodsky1990,Gao2001} and  the van der Waals (Casimir-Polder) form $\sim - \frac{1}{r^k}$ with $k=6$ ($7$)~\cite{Appelquist1978}  cannot reproduce the lattice data.
 
\begin{table}[hbtp]
\caption{The fit parameters in Eq.~(\ref{eq-fit}) with statistical errors quoted in the parentheses at $t/a=14$.
The fit range is $0<r<3.0$ fm. In $a_3 m_\pi^{4n}$, we take  $n=1$ and $n=0$ for fit A and B, respectively.
A$_{\rm erfc}$ (A$_{\rm exp}$) denotes fit A with the Nijmegen-type (Argonne-type) form factor.}
\begin{tabular}{cccc}
  \hline\hline
    fit~~&A$_{\rm erfc}$~~~&A$_{\rm exp}$~~~&B\\
  \hline
    $a_1$ [MeV] ~~~&-376(20)~~~&-371(27)~~~&-371(19)\\
    $b_1$ [fm] ~~~&0.14(1)~~~&0.13(1)~~~&0.15(3)\\
    $a_2$ [MeV] ~~~&306(122)~~~&-119(39)~~~&-50(35) \\
    $b_2$ [fm] ~~~&0.46(4)~~~&0.30(5)~~~ &0.66(61)\\
    $a_3 m_\pi^{4n}$ [MeV$\cdot{\rm fm}^{2n}$]~~~&-95(13)~~~&-97(14)~~~&-31(53)\\
    $b_3$ [fm] ~~&0.41(7)~~~&0.63(4)~~~&1.09(41) \\
  \hline\hline
\end{tabular} \label{tab-prm}
\end{table}


\begin{figure}[hbtp]
  \centering
  \includegraphics[width=8cm]{Fig3}
  \caption{
    The $N$-$\phi$ scattering phase shifts $\delta^{(3/2)}_0$ in the $^4S_{3/2}$ channel obtained from $V_{\rm fit}(r)$ at $t/a=12$  (red squares),  $13$ (green circles), and $14$ (blue triangles).
  } \label{Fig3}
\end{figure}



Figure~\ref{Fig3} shows the $N$-$\phi$ scattering phase shifts $\delta^{(3/2)}_0$ in the $^4S_{3/2}$ channel as a function of the center of mass kinetic energy $E_{\mathrm{CM}}=\sqrt{m_N^2+{k}^2}+\sqrt{m_\phi^2+{k}^2}-(m_N+m_\phi)$ obtained by using $V_{\rm fit}(r)$ with the fit A$_{\rm erfc}$. 
The scattering phase shifts from different $t$ are consistent with each other within statistical errors.
 The scattering length $a^{(3/2)}_0$ and the effective range $r^{(3/2)}_{\rm eff}$ can be extracted from the 
 effective range expansion for small $k$ as $k\cot\delta^{(3/2)}_0(k)=-\frac{1}{a^{(3/2)}_0}+\frac12 r^{(3/2)}_{\rm eff}k^2 +O(k^4)$.
 In Table~\ref{tab-scattering}, $a^{(3/2)}_0$ and $r^{(3/2)}_{\rm eff}$ are shown for the present  pion mass $m_{\pi}=146.4$ MeV; 
 the central values and the statistical errors are obtained from the data at $t/a=14$ with the  fit A$_{\rm erfc}$,  while 
  the systematic errors in the second parentheses are estimated by comparing results for $t/a=12$-14 with A$_{\rm erfc, exp}$ and B.
Other possible systematic errors are as follows; (i) the finite volume effect of the potential, which is expected to be small as  $\exp (-2m_{\pi} (L/2)a) \lesssim 0.3\%$ due to the large volume, and (ii) the finite cutoff effect, which is also expected to be small as $O((a\Lambda_{\rm QCD})^2)\sim O(1)\%$ due to the non-perturbative $O(a)$-improvement.
As an alternative estimate, we study how much scattering parameters are sensitive to the potential in short range, where the finite cutoff effect should be most prominent.
It is found that $a^{(3/2)}_0$ and $r^{(3/2)}_{\rm eff}$ change only $\sim 2\%$ even when we completely cut the short-range potential at $r < 0.1$ fm.

\begin{table}[tbhp]
\caption{The scattering length $a^{(3/2)}_0$ and the effective range $r^{(3/2)}_{\rm eff}$ obtained by using $V_{\rm fit}(r)$ at $m_\pi=146.4$ MeV
 with statistical and systematic errors. Estimated central values using  a model-dependent extrapolation  of $V_{\rm fit}(r)$ to $m_\pi=138.0$ MeV are 
  also shown for comparison.}
\begin{tabular}{ccc}
  \hline\hline
    $m_\pi$ [MeV]~&$a^{(3/2)}_0$ [fm]~~~&$r^{(3/2)}_{\rm eff}$ [fm]\\
  \hline
    146.4 ~~~&$-1.43(23)_{\rm stat.}\left(^{+36}_{-06}\right)_{\rm syst.}$~~~&$2.36(10)_{\rm stat.}\left(^{+02}_{-48}\right)_{\rm syst.}$\\
      138.0 ~~~&$\simeq -1.25$~~~&$\simeq 2.49$\\

 \hline\hline
\end{tabular} \label{tab-scattering}
\end{table}

To estimate how the scattering parameters change toward the physical quark mass,
we keep $a_{1,2,3}$ and $b_{1,2,3}$ in $V_{\rm fit}(r)$ fixed in fit A$_{\rm erfc, exp}$ and smoothly change the long-range potential by taking the isospin-averaged physical pion mass $m_{\pi} =138.0$ MeV  in the region where the TPE is dominated ($r>1.0$ fm from Fig.\ref{Fig2}).
By calculating the scattering phase shifts with such a potential with the physical masses of  $\phi$ and $N$, we obtain estimated values of  $a^{(3/2)}_0$ and $r^{(3/2)}_{\rm eff}$ for  $m_{\pi} =138.0$ MeV in Table~\ref{tab-scattering}. 
 Although the range of the TPE is increased by the smaller pion mass, the characteristic $m_{\pi}^4$ behavior  of the TPE strength makes the overall attraction weaker.
Note that this is only a model-dependent qualitative estimate and needs to be confirmed by future physical-point simulations.

\begin{figure}[tbhp]
  \centering
  \includegraphics[width=8cm]{Fig4}
  \caption{
    The scattering length $a^{(3/2)}_0$ obtained from $V(r;r_{\rm c})=\theta (r_{\rm c}-r)  V_{\rm fit}(r)$ as a function of  cutoff length $r_{\rm c}$ at $t/a=12$  (red squares),  $13$ (green circles), and $14$ (blue triangles).
  } \label{Fig4}
\end{figure}

Since we do not have reliable information on the $N\phi(^2S_{1/2})$ potential from lattice QCD at the moment due to the effect of the open channels,
a comparison of our results with spin-averaged scattering parameters should be made with caution. With this reservation in mind, our $a^{(3/2)}_0$ is found to be one or two orders of magnitude larger than the previous theoretical results in QCD sum rules~\cite{Koike1997,Klingl1997kf}.
Such a discrepancy may be due to the difficulty of obtaining the long-range TPE contribution from the low-order truncation of the operator product expansion in QCD sum rules.  In fact,  the magnitude of $a^{(3/2)}_0$ becomes considerably smaller when the long-range potential is cut off. 
Shown in Fig.~\ref{Fig4} is $a^{(3/2)}_0$ as a function of cutoff length $r_{\rm c}$ obtained by the potential, 
$  V(r;r_{\rm c}) = \theta (r_{\rm c}-r)  V_{\rm fit}(r)$ with the fit A$_{\rm erfc}$.
Considerable decrease of  $|a^{(3/2)}_0|$ from $1.43$ fm at $r_{\rm c}=\infty$ to about $0.1$ fm at $r_{\rm c}=0.5$ fm can be seen.


\noindent{\it Summary.$-$}
In this Letter, we present a first lattice QCD calculation on the interaction of the  $N$-$\phi$ system in the $^4S_{3/2}$ 
channel based on the  ($2+1$)-flavor simulations with nearly physical quark masses.
The interaction potential in the  $N\phi(^4S_{3/2})$ channel is extracted from lattice data of the 
hadronic spacetime correlation using the HAL QCD method.
The potential is found to be attractive for all distances and appears to be a combination of an attractive core at short distances and a two-pion exchange (TPE) tail at long distances ($r >$ 1 fm).
The latter is well fitted by the characteristic form of the TPE obtained by the interaction of a color-dipole and the nucleon.
The scattering parameters obtained from our potential at $m_{\pi}=146.4$ MeV is summarized in Table \ref{tab-scattering}.
By examining the potential fitted to the lattice data, we find that the scattering length $a^{(3/2)}_0$ is sensitive to the length scale of $r > 0.5$ fm.  
Also, we suggest that the $N$-$\phi$ attraction could be weaker at the physical pion mass due to the characteristic $m^4_{\pi}$ dependence of the strength of the TPE.
 
Our $a^{(3/2)}_0$ is substantially larger in magnitude than the previous calculations of the spin-averaged $a_0$ using QCD sum rules but is comparable to the spin-averaged $a_0$  by ALICE Collaboration within the error bar ~\cite{ALICE2021}.
Also, our $r^{(3/2)}_{\rm eff}$ is about three times smaller than the spin-averaged $r_{\rm eff}$ by ALICE Collaboration.
To make a firm comparison between theory and experiment, we need to perform a $\Lambda K$-$\Sigma K$-$N\phi$ coupled-channel analysis to reveal the complex dynamics of the $^2S_{1/2}$ and $^4S_{3/2}$  
channels as well as to perform physical-point simulations.


The present lattice QCD study near the physical point provides a first step to exploring the interaction of $s\bar{s}$ with the nucleon
 from first principle.  The heavier system such as $c\bar{c}$ interacting with the nucleon pioneered in ~\cite{Kawanai:2010ev,Sugiura:2019pye}
  is also an important problem to be studied near the physical point. 

\begin{acknowledgments}
{\it Acknowledgments.$-$}
We thank members of the HAL QCD Collaboration for stimulating discussions. 
Y. L. thanks Xu Feng for the helpful discussions.  We thank members of the PACS Collaboration for the gauge configuration generation conducted on the K computer at RIKEN.
The lattice QCD measurements have been performed on Fugaku and HOKUSAI supercomputers at RIKEN.
We thank ILDG/JLDG \cite{ldg}, which serves as essential infrastructure in this study.
This work was partially supported by HPCI System Research Project (hp120281, hp130023, hp140209, hp150223, hp150262, hp160211, hp170230, hp170170, hp180117, hp190103, hp200130, hp210165, and hp210212),  the National Key R\&D Program of China (Contract Nos. 2017YFE0116700 and 2018YFA0404400), the National Natural Science Foundation of China (Grant Nos. 11935003, 11975031, 11875075, and 12070131001), the JSPS (Grant Nos. JP18H05236, JP16H03978, JP19K03879, JP18H05407, and JP21K03555), the MOST-RIKEN Joint Project ``Ab initio investigation in nuclear physics'', ``Priority Issue on Post-K computer'' (Elucidation of the Fundamental Laws and Evolution of the Universe), ``Program for Promoting Researches on the Supercomputer Fugaku'' (Simulation for basic science: from fundamental laws of particles to creation of nuclei), and Joint Institute for Computational Fundamental Science (JICFuS).

\end{acknowledgments}



\

