\section{Introduction}

The semi-leptonic and non-leptonic decays of $ B $ mesons are the frequently studied decays and also the dominant production channels of charmed mesons. During the last decades, for many important cases such as providing precise value of CKM element $ V_{cb} $,  the channels of $ B $ decays to S-wave ground states of $ D $ mesons have been extensively measured and studied by the ALEPH, CLEO, OPAL, BABAR and Belle Collaborations\cite{ALEPH1997.373,CLEO1999.3746,OPAL2000.15,Belle2002.258,BABAR2008.032002,BABAR2010.011802,Belle2010.112007,Belle2016.032006} besides the theoretical studies.

In recent years, many collaborations reported several charmed resonances including some orbitally excited $ D $ mesons, which attracts lots of attention.
The Belle and BABAR Collaborations reported the semi-leptonic $ B $ decays to P-wave $ D^* $ mesons by using fully reconstructed $ B $ tags \cite{BABAR2008,Belle2008} and the Belle, BABAR and LHCb Collaborations reported the non-leptonic decays $ B \to D_0^*\pi(K) $ \cite{Belle2004,Belle2007.012006,BABAR2009.112004,LHCb2015.032002,LHCb2015.092002}.
They inspired many theoretical studies on the excited charmed states using different models, for example, the light-front quark model \cite{Kang2018}, the constituent quark model\cite{Segovia2011.094029}, as well as the Bethe-Salpeter method \cite{FuHF2011}, etc.

In 2013, the LHCb Collaboration reported several resonances around 3000 MeV, $ D_J(3000)^0 $ and $ D_J^*(3000)^{+,0} $\cite{LHCb2013}. The $ D_J(3000) $ and $ D_J^*(3000) $ were observed in the $ D^*\pi $ and $D\pi$ invariant mass spectrum respectively. Their quantum numbers $J^P$ are still undetermined and many theoretical studies give different assignments \cite{Sun2013,Lu2014,YuGL2015,Godfrey2016}. In our previous works\cite{lsc2017,lsc2018}, we calculated the strong decays and the leptonic productions of $ D_J(3000)$ and we favoured it as the excited $2P(1^{+'})$ broad state . For $D_J^*(3000)$, we calculated its strong decays and our results favoured it as the excited scalar $2P(0^{+})$ state \cite{Tan2018}.

We notice that, in current experiments and theories, the knowledge of $ B $  semi-leptonic and non-leptonic decays to orbitally excited $ D^{*}_0 $ meson is still rather poor. Thus this work will focus on the leptonic decays $ B \to D_0^* \ell^- \overline{\nu}_{\ell} $ and non-leptonic decays $B \to D_0^* X$, where the initial state could be $B^-$ or $\overline{B}^0$, the final state $D_0^*$ is the excited scalar $D_0^{*}(nP)^0$ or $D_0^{*}(nP)^+$ ( $n=1,~2,~3$), and $X$ is a light meson. Currently, the $1P$ states $D_0^{*}(2400)^0$ and $D_0^{*}(2400)^+$ have been well studied, while $2P$ and $3P$ states haven't. The newly detected $D_J^*(3000)^{+,0}$ are treated as the $2P$ scalars in this paper and our results will help to determine their quantum numbers. The processes of $B$ leptonic decays to them could be their important production ways.

In our previous study \cite{gzk2019}, we found that large relativistic corrections exist in the processes where a heavy-light excited state is involved. We also found that the highly excited state has larger relativistic effect than its corresponding ground state. Thus when a process includes an excited state, a relativistic method or model is needed. In this paper, we use the Bethe-Salper (BS) method based on the relativistic BS equation. The relativistic effect is well concerned by solving the BS equation and applying the BS wave function.

The rest contents of this paper are organized as follows: in section 2, we present the formalism of semi-leptonic production process, including the leptonic and hadronic matrix elements by using the BS method. Then the factorization approach is used to derive the formalism of non-leptonic process in section 3. In section 4, we show our numerical results and comparison with the results of other model. Finally, discussions and short summary are given in section 4.

\section{Formalism of semi-leptonic decays}

We take $ B^{-}(\overline{B}^0) \to D_{0}^{*0}(D_{0}^{*+})\ell^- \overline{\nu}_{\ell}$ as an example to show the calculation details of semi-leptonic process. The Feynman diagram is shows in Fig. \ref{feynman-sl}.
\begin{figure}[htb!]
	\centering
	\includegraphics[scale=0.5]{Feynman_B_D_semi-leptonic-eps-converted-to.pdf}
	\caption[fig1]{Feynman diagram of the semi-leptonic decay $B^{-}(\overline{B}^0)\to D_0^{*0(+)} \ell^- \overline{\nu}_{\ell}$.}
	\label{feynman-sl}
\end{figure}

The transition amplitude $T  $ can be expressed as :
\begin{equation}
	T=\frac{G_F}{\sqrt{2}}V_{cb} l^{\xi} \left< D_0^*(P_f) | J_{\xi} | B^-(P) \right>,
\end{equation}
where $ G_F $ is the Fermi weak coupling constant, $ V_{cb} $ is the CKM matrix element, $J_{\xi}$ is the charged weak current and $ l^{\xi}  $ is the leptonic matrix.
\begin{align}
	& J_{\xi} = \bar{c} \gamma_{\xi}(1-\gamma_{5})b , \\
	& l^{\xi} = \bar{u}_{\ell}(p_l)\gamma^{\xi}(1-\gamma_5){v_{\nu}}(p_{{\nu}}) .
\end{align}

Then the square of amplitude can be expressed by the function of hadronic and leptonic tensor,
\begin{equation}
	|T|^2 = \frac{G_F^2}{2}|V_{cb}|^2 l^{\xi \xi '} h_{\xi \xi '},
\end{equation}
where the leptonic tensor can be written as following form:
\begin{equation}
\begin{split}
l^{\xi \xi '} = 8(p_{\nu}^{\xi} p_{\ell}^{\xi '} + p_{\ell}^{\xi}p_{\nu}^{\xi '} - p_{\ell} \cdot p_{\nu} g^{\xi \xi '} + i \varepsilon ^{\xi \xi' p_{\ell} p_{\nu}}).
\end{split}
\end{equation}

We derive the hadronic matrix element by using the relativistic BS method
\begin{equation}
	\begin{split}
		h_{\xi}&=\left< D_0^*(P_f) | J_{\xi} | B^-(P) \right> \\
		 &=  \int \frac{{\rm{d}}^3 q}{(2\pi)^3} \mathrm{Tr}\left[ \frac{\slashed{P}}{M} \overline{\varphi}_{P_f}(q_{f\perp}) \gamma_{\xi} (1-\gamma_{5}) \varphi_{P}(q_{\perp}) \right] = n_1 P_{\xi} + n_2 P_{f\xi} ,
	\end{split}
\end{equation}
where $q_{f}=q-\frac{m_u}{m_b+m_u}P_{f}$, $ \varphi_{P}  $ and $ \varphi_{P_f}  $ are the instantaneous BS wave functions of the initial state $ B $ meson and final state $ D_0^* $ mesons. They are obtained by completely solving the BS equation. The processes of solving the Salpeter equation and obtaining the wave functions are not shown here. More details can be found in our previous works \cite{WangGL2004,WangGL2006, WangGL2007}. We just give a brief review in the appendix.  $ n_i $ are the form factors whose results are shown in next section.

$ h_{\xi \xi '} $ is the hadronic tensor,
\begin{equation}
\label{htesor}
	\begin{split}
	h_{\xi \xi '} =& \left< D_0^*(P_f) | J_{\xi} | B^-(P) \right>^{\dagger} \left< D_0^*(P_f) | J_{\xi '} | B^-(P) \right> \\
	=& n_1^2 P_{\xi}P_{ \xi '} + n_2^2 P_{f \xi }P_{f \xi '} + n_1n_2 (P_{\xi}P_{f\xi '}+P_{f \xi} P_{\xi '}).
	\end{split}
\end{equation}

Then the decay width can be given by the phase-space integral
\begin{equation}
	\Gamma=\frac{1}{2M} \int \frac{{\rm{d}}^3 \bm{P}_f {\rm{d}}^3 \bm{p}_{\ell} {\rm{d}}^3 \bm{p}_{{\nu}_{\ell}}}{(2\pi)^9 2E_f 2E_{\ell} 2E_{{\nu}_{\ell}}} (2\pi)^4\delta^4(P-P_f-p_{\ell}-p_{{\nu}_{\ell}})  |T|^2.
\end{equation}
After the simplification, it can be rewritten as
\begin{equation}
	\Gamma=\frac{1}{64\pi M } \int   |T|^2  \frac{|\bm{P}_f|}{E_f} {\rm{d}} |\bm{P}_f |  \frac{|\bm{p}_{\ell}|}{E_{\ell}} {\rm{d}} |\bm{p}_{\ell} |.
\end{equation}

\section{Formalism of non-leptonic decays}

The Feynman diagram of non-leptonic decay $B^{-}\to D_0^{*0} X$ or $\overline{B}^{0}\to D_0^{*+} X$ is shown in Fig. \ref{feynman-nl}, where $ X $ could be $ \pi^-,~K^-, ~\rho^-$ or $ K^{*-} $.

\begin{figure}[htb!]
	\centering
	\includegraphics[scale=0.5]{Feynman_B_D_non-leptonic-eps-converted-to.pdf}
	\caption[fig2]{Feynman diagram of the non-leptonic decay $B^{-}(\overline{B}^0)\to D_0^{*0(+)} X$.}
	\label{feynman-nl}
\end{figure}

The effective Hamiltonian $\mathcal{H}_{eff}  $ can be expressed as \cite{Ali1998.094009,Choi2009.114003}:
\begin{equation}
	\begin{split}
	\mathcal{H}_{eff} = \frac{G_F}{2}\left[ V_{cb}V^*_{uq}(c_1(\mu) O_1 + c_2(\mu) O_2) + h.c.  \right],
	\end{split}
\end{equation}
where, $ O_i $ is the 4-quark operator containing the charged weak current $ (\bar{q}u)_{V-A} $ and $ (\bar{c}b)_{V-A} $; $ c_i $ is the Wilson coefficient which depends on the renormalization scale $ \mu $.

By using the factorization approach, the transition matrix elements$ \left< D_0^*X|\mathcal{H}_{eff}|B \right> $ involving the 4-quark operators can be split into the product of two matrix elements $\left< D_0^*|(\bar{c}b)_{V-A}|B \right> $ and $ \left< X|(\bar{q}u)_{V-A}|0 \right> $ \cite{Fakirov1978.315,Bauer1987.103,Ali1998.094009}.

Then the transition amplitude can be expressed as
\begin{equation}
	T=\frac{G_F}{\sqrt{2}} V_{cb} V^*_{uq} a_1(\mu) \left<D_0^{*0} | (\bar{c}b)_{V-A} | B^- \right> \left< X|(\bar{q}u)_{V-A}|0 \right>,
\end{equation}
where, $ q $ denotes the $ d $ or $ s $ quark; $ a_1=c_1+\frac{1}{N_c}c_2 $ is the effective Wilson coefficient, where $ N_c=3 $ is the number of colors. We choose the scale $ \mu \approx m_b$ for $ B $ decays and adopt the effective Wilson coefficient $ a_1=1.14 $ \cite{Ivanov2006.054024}. The annihilation matrix element can be written as
\begin{equation}
	\begin{split}
	\left< X_P | (\bar{q}u)^{\mu}_{V-A} |0 \right> &= -i f_{X} P_{X}^{\mu}, \\
	\left< X_V | (\bar{q}u)^{\mu}_{V-A} |0 \right> &=  f_{X} M_{X}\varepsilon_{X}^{\mu},
	\end{split}
\end{equation}
where $ X_P $ means pseudoscalar ($ \pi , K $) and $ X_V $ means vector mesons ($ \rho , K^* $);  $ f_X $ is the corresponding decay constant; $ \varepsilon^{\mu} $ is the polarization vector of $ X_V $ and it satisfies the completeness relation $  \sum\limits \varepsilon_{\lambda}^{\mu} \varepsilon_{\lambda}^{\nu}=\frac{P_X^{\mu}P_X^{\nu}}{M_X^2}-g^{\mu\nu}$.

Same as the semi-leptonic case, we write the square of amplitude by the hadronic and light meson tensor
\begin{equation}
	|T|^2 = \frac{G_F^2}{2}|V_{ub}|^2 |V_{uq}|^2 a_1^2 h_{\mu\nu} X^{\mu\nu},
\end{equation}
where $ h_{\mu\nu} $ is same as Eq. \ref{htesor}, and the light meson tensor is
\begin{equation}
	\begin{split}
	X_P^{\mu\nu} &= P_X^{\mu}P_X^{\nu}f_X^2,	\\
	X_V^{\mu\nu} &= (P_X^{\mu}P_X^{\nu} - M_X^2 g^{\mu\nu}) f_X^2.
	\end{split}
\end{equation}

Then the decay width can be obtained by
\begin{equation}
	\Gamma = \int \frac{|\bm{P}_f|}{32\pi^2 M^2} \sum |T|^2 {\rm{d}} \Omega.
\end{equation}

\section{RESULTS AND DISCUSSION}

In our calculations, we adopt the same parameters as what we used before\cite{Tan2018}: $m_u=0.305{\ \rm{GeV}},\ m_d=0.311{\ \rm{GeV}},\ m_s=0.50{\ \rm{GeV}},\   m_c=1.62{\ \rm{GeV}}$,  $\alpha=0.060{\ \rm{GeV}}$, $\lambda=0.210{\ \rm{GeV}}^2$ and $\Lambda_{QCD}=0.270 {\ \rm{GeV}}$. The involved mesons' masses are: $M_{D_0^*(2400)^0}=2.318{\ \rm{GeV}} $, $M_{D_0^*(2400)^+}=2.351 {\ \rm{GeV}}$,
 ${ M_{D_J^*(3000)^{(0,+)}}=3.008 {\ \rm{GeV}}}$, $ M_{D_0^*(3P)^{0,+}}=3.183{\ \rm{GeV}}$, $M_{B^0}=5.2796{\ \rm{GeV}}$ and  $ M_{B^{\pm}}=5.2793{\ \rm{GeV}}$.

The CKM matrix elements\cite{PDG2018} and the involved mesons' decay constants are\cite{PDG2018,LHCb2013,LiQ2017}: $ |V_{cb}|=0.0422 $, $ |V_{ud}|=0.9742 $, $ |V_{us}|=0.2243 $,
$f_\rho=205{\ \rm{MeV}} $, $f_\pi=130.4{\ \rm{MeV}}$, $   f_\eta=130{\ \rm{MeV}}$ and $  f_k=156.2{\ \rm{MeV}}$.

\subsection{Semi-leptonic decays}

The form factors relevant to the hadronic transition matrix elements of $B\to D_0^*(1P-3P) \ell \nu_{\ell}$ are shown in Fig. \ref{ff}, where $ t=(P-P_f)^2 $ and  $ t_m $ is the momentum transfer at the zero recoil(the maximum value of $t$).

\begin{figure}[htb!]
	\centering
	\subfigure[The form factors for $D_0^*(2400),\ l=e $ or $\mu $]{
		\label{ff2400e}
		\includegraphics[width=0.45\textwidth]{formfactor_2400_0_e_zxjs_May_2-eps-converted-to.pdf}}
	\hspace{0.2cm}
	\subfigure[The form factors for $D_0^*(2400),\ l=\tau $]{
		\label{ff2400tau}
		\includegraphics[width=0.45\textwidth]{formfactor_2400_0_tau_zxjs_May_2-eps-converted-to.pdf}}

	\subfigure[The form factors for $D_J^*(3000),\ l=e$ or $\mu $]{
		\label{ff3000e}
		\includegraphics[width=0.45\textwidth]{formfactor_3000_0_e_zxjs_May_2-eps-converted-to.pdf}}
	\hspace{0.2cm}
	\subfigure[The form factors for $D_J^*(3000),\ l=\tau $]{
		\label{ff3000tau}
		\includegraphics[width=0.45\textwidth]{formfactor_3000_0_tau_zxjs_May_2-eps-converted-to.pdf}}

	\subfigure[The form factors for $D_0^*(3P),\ l=e$ or $\mu $]{
		\label{ff3183e}
		\includegraphics[width=0.45\textwidth]{formfactor_3183_0_e_zxjs_May_2-eps-converted-to.pdf}}
	\hspace{0.2cm}
	\subfigure[The form factors for $D_0^*(3P),\ l=\tau $]{
		\label{ff3183tau}
		\includegraphics[width=0.45\textwidth]{formfactor_3183_0_tau_zxjs_May_2-eps-converted-to.pdf}}
	\caption{The form factors for the semi-leptonic production process of $D_0^*(2400)$ , $D_J^*(3000)$ and $ D_0^*(3P) $ , where $ S^+=(n_1+n_2)/2,\ S^-=(n_1-n_2)/2 $.}
	\label{ff}
\end{figure}

Table \ref{slresuluts2400} shows the decay widths and  branching ratios of semi-leptonic production of the ground state $ D_0^*(2400) $.  With varying the parameters by $ \pm 5 \% $, we can obtain the uncertainty of the results. Because there is almost no difference between the results of $ l=e $ and $ l=\mu $, only the values of $ l=e $ are given below.

\begin{table}[htb!]
	\renewcommand\arraystretch{0.99}
	\caption[results]{The results of semi-leptonic production of $ 1P $ state $ D_0^*(2400) $ }
	\label{slresuluts2400}
	\vspace{0.5em}\centering
	\begin{tabular}{lcc}
		\toprule[1.5pt]
		\multicolumn{1}{c}{Channels} & Widths (GeV) & BR  \\
		\midrule[1pt]
		$ B^{-} \to D_0^*(2400)^0 e^- \overline{\nu}_e $ & $1.03 \sim 1.70  \times 10^{-15} $ &  $ 2.57 \sim 4.24 \times 10^{-3} $ \\
		 $ \overline{B}^{0} \to D_0^*(2400)^+ e^- \overline{\nu}_e $ & $1.04 \sim 1.69  \times 10^{-15} $ &  $ 2.39 \sim 3.91 \times 10^{-3} $  \\
		$ B^{-} \to D_0^*(2400)^0 \tau^- \overline{\nu}_{\tau} $ & $ 1.20 \sim 1.72 \times 10^{-16} $ & $ 2.98 \sim 4.28 \times 10^{-4} $ \\
		 $ \overline{B}^{0} \to D_0^*(2400)^+ \tau^- \overline{\nu}_{\tau} $ & $ 1.13 \sim 1.60 \times 10^{-16} $ & $ 2.62 \sim 3.70 \times 10^{-4} $ \\
		\bottomrule[1.5pt]
	\end{tabular}
\end{table}

For comparison, we also give the results of cascade decays to $ D\pi  $ which are shown in Table \ref{cascadesl}.  Recently, Ref.\cite{Kang2018}  used covariant light-front quark model to calculate the channel $ B^{+} \to D_0^*(2400)^0 l^+ \overline{\nu}_l   \to D^- \pi^+ $ , whose result of branching ratio is $ 2.31\pm 0.25 \times 10^{-3} $. And Ref.\cite{Segovia2011.094029} shows that $ \mathcal{B}(B^{+} \rightarrow \overline{D}_{0}^{* 0} l^{+} \nu_{l}) \times \mathcal{B}(\overline{D}_{0}^{* 0} \rightarrow D^{-} \pi^{+})=2.15 \times 10^{-3}  $ and $ \mathcal{B}\left(B^{0} \rightarrow D_{0}^{*-} l^{+} \nu_{l}\right) \times \mathcal{B}\left(D_{0}^{*-} \rightarrow D^{0} \pi^{-}\right)=1.80 \times 10^{-3} $ with the constituent quark model.
Considering the uncertainty, our results are consistent with the experimental and other models' results.

\begin{table}[htb!]
	\renewcommand\arraystretch{0.90}
	\caption[results]{The branching ratios($ \times 10^{-3} $)  of cascade decays to $ D \pi $(with $ l=e $).  The branching ratio $\mathcal{B}( D_0^*(2400) \to D \pi ) $ is from our previous work\cite{Tan2018}. }
	\label{cascadesl}
	\vspace{0.5em}\centering
	\begin{tabular}{lcccc}
		\toprule[1.5pt]
		\multicolumn{1}{c}{Channels} &  Ours  &  BABAR\cite{BABAR2008}   &   Belle\cite{Belle2008}   & PDG\cite{PDG2018}  \\
		\midrule[1pt]
		$ \mathcal{B}  (B^{-} \to D_0^*(2400)^0 e^- \overline{\nu}_l  \to D^0 \pi^0 )$ &  $ 0.85 \sim 1.39  $ & $ $ &  $  $   \\
		$ \mathcal{B}( B^{-} \to D_0^*(2400)^0 e^- \overline{\nu}_l   \to D^+ \pi^- )$ &  $ 1.72 \sim 2.81   $ & $ 2.6 \pm 0.9  $ &  $ 2.4 \pm 1.0  $ & $ 2.5\pm 0.5  $  \\
		$ \mathcal{B}( \overline{B}^{0} \to D_0^*(2400)^+ e^- \overline{\nu}_l   \to D^0 \pi^+ )$ &  $ 1.60 \sim 2.61  $ &  $ 4.4 \pm 1.4  $ &  $ 2.0 \pm 1.2  $  & $ 3.0\pm 1.2 $ \\
		$ \mathcal{B}( \overline{B}^{0} \to D_0^*(2400)^+ e^- \overline{\nu}_l   \to D^+ \pi^0 )$ &  $ 0.79 \sim 1.29  $  & $ $ &  $  $  \\
		\bottomrule[1.5pt]
	\end{tabular}
\end{table}

Then we use the same method to calculate the $2P$ state $D_J(3000)$, which is shown in Table \ref{slresuluts3000}. Unlike the ground state, the results of $ 2P $ state are much lower and have large uncertainty by varying the input parameters, while the predicted results of $ 3P $ state are given in Table \ref{slresuluts3183} and get smaller uncertainty.

\begin{table}[htb!]
	\renewcommand\arraystretch{1.01}
	\caption[results]{The results of semi-leptonic production of  $ 2P $ state $ D_J^*(3000) $ }
	\label{slresuluts3000}
	\vspace{0.5em}\centering
	\begin{tabular}{lcc}
		\toprule[1.5pt]
		\multicolumn{1}{c}{Channels} & Widths (GeV) & BR   \\
		\midrule[1pt]
		$ B^{-} \to D_J^*(3000)^0 e^- \overline{\nu}_e $ & $ 1.6	\sim 80.2 \times 10^{-18} $ & $ 0.4	\sim 20.0 \times 10^{-5} $ \\
		 $ \overline{B}^{0} \to D_J^*(3000)^+ e^- \overline{\nu}_e $ & $ 1.2 \sim 83.0 \times 10^{-18} $ & $ 0.3 \sim 19.2 \times 10^{-5} $  \\
		$ B^{-} \to D_J^*(3000)^0 \tau^- \overline{\nu}_{\tau} $ & $ 0.4 \sim 69.8 \times10^{-20} $ & $ 0.1 \sim 17.3 \times 10^{-7} $ \\
		 $ \overline{B}^{0} \to D_J^*(3000)^+ \tau^- \overline{\nu}_{\tau} $ & $ 0.8 \sim 72.9\times10^{-20} $ & $ 0.2 \sim 16.8 \times 10^{-7} $ \\
		\bottomrule[1.5pt]
	\end{tabular}
\end{table}

\begin{table}[htb!]
	\renewcommand\arraystretch{0.99}
	\caption[results]{The results of semi-leptonic production of $ 3P $ state $ D_0^*$ }
	\label{slresuluts3183}
	\vspace{0.5em}\centering
		\begin{tabular}{lcc}
			\toprule[1.5pt]
			\multicolumn{1}{c}{Channels} & Widths (GeV) & BR   \\
			\midrule[1pt]
			$ B^{-} \to D_0^*(3P)^0 e^- \overline{\nu}_e $ & $3.10 \sim 7.37  \times 10^{-18} $ &  $ 0.77 \sim 1.83\times 10^{-5} $ \\
			 $ \overline{B}^{0} \to D_0^*(3P)^+ e^- \overline{\nu}_e $ & $ 3.24 \sim 7.88  \times 10^{-18} $ &  $ 0.75 \sim 1.82 \times 10^{-5} $  \\
			$ B^{-} \to D_0^*(3P)^0 \tau^- \overline{\nu}_{\tau} $ & $ 2.49 \sim 3.67 \times 10^{-20} $ & $ 6.19 \sim 9.13 \times 10^{-8} $ \\
			 $ \overline{B}^{0} \to D_0^*(3P)^+ \tau^- \overline{\nu}_{\tau} $ & $ 2.45 \sim 3.78 \times 10^{-20} $ & $ 5.65 \sim 8.73 \times 10^{-8} $ \\
			\bottomrule[1.5pt]
	\end{tabular}
\end{table}

Why the same parameters varying leads to the abnormal results of excited states? We consider that the different structures of BS wave function possibly play an important role here.   Fig. \ref{wfBD} shows the wave function values changing with the relative momentum$ |q| $. The wave functions of $ B^- $ and $ 1P $ state are all positive without nodes. When we varying the input parameters, the curve will have some small shift. And the shift could cause the small uncertainty in the overlapping integral. For excited states, we can find that the wave functions have nodes . For $ 2P $ states, the wave function changes from positive to negative after the nodes. In the overlapping integral, it causes the cancellation and the final results will be highly suppressed. Then if we vary the input parameters, a small shift of the wave function could cause a large uncertainty.

For $ 3P $ states, the wave function has two nodes, and the value change from negative to positive after the second node. The cancellation gets smaller than the case of the $ 2P $ state. Thus, the final branching ratio seems to be fine and the uncertainty is not very large.

\begin{figure}[htb!]
	\centering
	\subfigure[Wave functions of  $ B^- $ ]{
		\label{0-wf}
		\includegraphics[width=0.45\textwidth]{wf_0minus_1S0_B-_May_2-eps-converted-to.pdf}}
	\hspace{0.2cm}
	\subfigure[Wave functions of $ 1P $ state $ D_0^*(2400) $]{
		\label{1Pwf}
		\includegraphics[width=0.45\textwidth]{wf_0plus_1P_3P0_D0star2400_0_May_2-eps-converted-to.pdf}}

	\subfigure[Wave functions of $ 2P $ state $ D_J^*(3000) $]{
		\label{2Pwf}
		\includegraphics[width=0.45\textwidth]{wf_0plus_2P_3P0_D0star3000_0_May_2-eps-converted-to.pdf}}
	\hspace{0.2cm}
	\subfigure[Wave functions of $ 3P $ state $ D_0^*(3P) $]{
		\label{3Pwf}
		\includegraphics[width=0.45\textwidth]{wf_0plus_3P_3P0_D0star3183_0_May_2-eps-converted-to.pdf}}
	\caption{Wave functions of $ B^- $ and $ D_0^* $.}
	\label{wfBD}
\end{figure}

Considering the masses of these states have errors, the branching ratio of their semi-leptonic production changing with their masses are given in Fig. \ref{brmass}. For $ 1P $ state $ D_0^*(2400) $, the results have small changes. But the branching ratio of $ 2P $ state $ D_J^*(3000)$ dramatically decrease to nearly zero with increase of mass. The mass changing will also cause the wave function shift. That means the overlapping integral cancellation will increase as the mass increasing. For $ 3P $ state, it can be seen that the curves of $ l=e$ have minimum points around $ m=3.175 {\ \rm{GeV}} $, which means the overlapping integrals have the maximum cancellation at that mass value. After that, the values increase again. For $ l= \tau $, because of the small phase space, the branching ratios have the downtrend from the beginning to the end.

\begin{figure}[htb!]
	\centering
	\subfigure[Branching ratio vers the mass of $D_0^*(2400), l=e$]{
		\label{br2400emu}
		\includegraphics[width=0.45\textwidth]{BR_mass_B-D0star2400_e_0_+_HPC_May_2-eps-converted-to.pdf}}
	\hspace{0.2cm}
	\subfigure[Braching ratio vers the mass of $D_0^*(2400), l=\tau$]{
		\label{br2400tau}
		\includegraphics[width=0.45\textwidth]{BR_mass_B-D0star2400_tau_0_+_HPC_May_2-eps-converted-to.pdf}}

	\subfigure[Braching ratio vers the mass of $D_J^*(3000), l=e$]{
		\label{br3000emu}
		\includegraphics[width=0.45\textwidth]{BR_mass_B-D0star3000_e_0_+_HPC_May_2-eps-converted-to.pdf}}
	\hspace{0.2cm}
	\subfigure[Braching ratio vers the mass of $D_J^*(3000), l=\tau$]{
		\label{br3000tau}
		\includegraphics[width=0.45\textwidth]{BR_mass_B-D0star3000_tau_0_+_HPC_May_2-eps-converted-to.pdf}}

	\subfigure[Braching ratio vers the mass of $D_J^*(3P), l=e$]{
		\label{br3183emu}
		\includegraphics[width=0.45\textwidth]{BR_mass_B-D0star3183_e_0_+_HPC_May_2-eps-converted-to.pdf}}
	\hspace{0.2cm}
	\subfigure[Braching ratio vers the mass of $D_J^*(3P), l=\tau$]{
		\label{br3183tau}
		\includegraphics[width=0.45\textwidth]{BR_mass_B-D0star3183_tau_0_+_HPC_May_2-eps-converted-to.pdf}}
	\caption{The branching ratios of  semi-leptonic production change with the mass of $ D_0^* $.}
	\label{brmass}
\end{figure}

Then, the normalized lepton spectra of the semi-leptonic production are presented in Fig. \ref{ls}. Because there are almost no difference between $ l=e $ and $ l=\mu $, only the channels of $ B^- \to D_0^{*0} e^- \overline{\nu}_e $ and $ B^- \to D_0^{*0} \tau^- \overline{\nu}_{\tau}  $ are given here. The spectrum peaks of $ 2P $ and $ 3P $ states move left because phase space decreases, especially for $ l=\tau $.

\begin{figure}[htb!]
	\centering
	\subfigure[Lepton spectra when $ l=e $]{
		\label{lse}
		\includegraphics[width=0.45\textwidth]{lepton_spectra_e_0_2-eps-converted-to.pdf}}
	\hspace{0.2cm}
	\subfigure[Lepton spectra when $ l=\tau $]{
		\label{lstau}
		\includegraphics[width=0.45\textwidth]{lepton_spectra_tau_0_2-eps-converted-to.pdf}}
	\caption{The spectra of differential decay width vs $ |\bm{p}_l| $(with
		normalization).}
	\label{ls}
\end{figure}

\subsection{Non-leptonic production}

The non-leptonic production results of $ 1P $ state $ D_0^*(2400) $ are shown in Table \ref{nonleptonic2400}. The $ D_0^*\pi $ and $ D_0^*\rho $ channels get the order of $ 10^{-3} $, while $ D_0^*K $ and $ D_0^*K^* $ modes are in the order of $ 10^{-5} $ and $ 10^{-4} $, respectively.

\begin{table}[htb!]
	\renewcommand\arraystretch{1.01}
	\caption[results]{Non-leptonic production results of the $ 1P $ state $ D_0^*(2400) $. }
	\label{nonleptonic2400}
	\vspace{0.5em}\centering
	\begin{tabular}{lcc}
		\toprule[1.5pt]
		\multicolumn{1}{c}{Channels} & Widths (GeV) & BR   \\
		\midrule[1pt]
		$ B^{-} \to D_0^*(2400)^0 \pi^-  $ & $ 3.35 \sim 5.78 \times 10^{-16} $ & $ 0.83 \sim 1.44 \times 10^{-3} $ \\
		 $ \overline{B}^{0} \to D_0^*(2400)^+ \pi^-  $ & $ 3.46 \sim 5.80 \times 10^{-16} $ & $ 0.80 \sim 1.34 \times 10^{-3} $  \\
		$ B^{-} \to D_0^*(2400)^0  K^- $ & $ 2.48 \sim 4.26 \times 10^{-17} $ & $ 0.62 \sim 1.06 \times 10^{-5} $ \\
		 $ \overline{B}^{0} \to D_0^*(2400)^+  K^- $ & $ 2.56 \sim 4.27 \times 10^{-17} $ & $ 5.91 \sim 9.85 \times 10^{-6} $ \\
		$ B^{-} \to D_0^*(2400)^0  \rho^- $ & $ 0.74 \sim 1.30\times10^{-15} $ & $ 1.85 \sim 3.23 \times 10^{-3} $ \\
		 $ \overline{B}^{0} \to D_0^*(2400)^+  \rho^- $ & $ 0.76 \sim 1.30 \times10^{-15} $ & $ 1.77 \sim 2.99 \times 10^{-3} $ \\
		$ B^{-} \to D_0^*(2400)^0 K^{*-}$ & $ 4.25 \sim 7.46 \times 10^{-17} $ & $ 1.06 \sim 1.86 \times 10^{-4} $ \\
		 $ \overline{B}^{0} \to D_0^*(2400)^+ K^{*-}$ & $ 4.37 \sim 7.44 \times 10^{-17} $ & $ 1.01 \sim 1.72 \times 10^{-4} $ \\
		\bottomrule[1.5pt]
	\end{tabular}
\end{table}

There have been some experimental results of their cascade decays. For comparison, our results of the cascade decays are shown in Table \ref{cascadenl}.  The Belle and BABAR Collaborations show that the branching ratios $ \mathcal{B}\left(B^{-} \rightarrow D_{0}^{* 0} \pi^{-}\right) \times \mathcal{B}\left(D_{0}^{* 0} \rightarrow D^{+} \pi^{-}\right)$ are $ (6.1 \pm 0.6 \pm 0.9 \pm 1.6) \times 10^{-4} $\cite{Belle2004} and $ (6.8 \pm 0.3 \pm 0.4 \pm 2.0) \times 10^{-4} $\cite{BABAR2009.112004}, respectively. Our results are consistent with them. For the charged $ D_0^{*\pm} $ meson, the Belle and LHCb Collaborations give the results $ \mathcal{B}(\overline{B}^{0} \rightarrow D_{0}^{*+} \pi^{-}) \times \mathcal{B}\left(D_{0}^{*+} \rightarrow D^{0} \pi^{+}\right) =(6.0 \pm 1.3 \pm 1.5 \pm 2.2) \times 10^{-5}$\cite{Belle2007.012006} and $ (7.7 \pm 0.5 \pm 0.3 \pm 0.3 \pm 0.4) \times 10^{-5}$\cite{LHCb2015.032002}, respectively, which are much lower than the branching ratio of the previous channel. For $ D_0^*K $ channel, the LHCb collaboration shows the result  $ \mathcal{B}\left(B^{-} \rightarrow D_0^{*0} K^{-}\right) \times \mathcal{B}\left(D_0^{*0} \rightarrow D^{+} \pi^{-}\right)=(6.1 \pm 1.9 \pm 0.5 \pm 1.4 \pm 0.4) \times 10^{-6} $\cite{LHCb2016.119901}, which is also lower than our calculation. From the perspective of symmetry, the non-leptonic results of these channels should be similar. But different experimental results show marked discrepancy, which need more experimental data accumulations and theoretical attentions.

\begin{table}[htb!]
	\renewcommand\arraystretch{0.90}
	\caption[results]{The branching ratios of cascade decays to $ D \pi $. }
	\label{cascadenl}
	\vspace{0.5em}\centering
	\resizebox{\textwidth}{!}{
		\begin{tabular}{llcllc}
			\toprule[1.5pt]
			\multicolumn{2}{c}{Channels} & BR & \multicolumn{2}{c}{Channels} & BR \\
			\midrule[1pt]
			\multirow{2}{*}{$ B^{-} \to D_0^*(2400)^0 \pi^- $}	& $  \to D^0 \pi^0 $ &  $ 2.76 \sim 4.76 \times 10^{-4} $ & \multirow{2}{*}{$ \overline{B}^{0} \to D_0^*(2400)^+ \pi^- $}	& $  \to D^0 \pi^+ $ &  $ 5.34 \sim 8.94 \times 10^{-4} $   \\
			& $ \to D^+ \pi^- $ &  $ 5.57 \sim 9.63  \times 10^{-4} $ & & $ \to D^+ \pi^0 $ &  $ 2.66 \sim 4.45  \times 10^{-4} $  \\

			\multirow{2}{*}{$ B^{-} \to D_0^*(2400)^0 K^- $}	& $  \to D^0 \pi^0 $ &  $ 2.04 \sim 3.51\times 10^{-5} $ & \multirow{2}{*}{$ \overline{B}^{0} \to D_0^*(2400)^+ K^- $}	& $  \to D^0 \pi^+ $ &  $ 3.95 \sim 6.58 \times 10^{-5} $   \\
			& $ \to D^+ \pi^- $ &  $ 4.12 \sim 7.09  \times 10^{-5} $ & & $ \to D^+ \pi^0 $ &  $ 1.96 \sim 3.27  \times 10^{-5} $  \\

			\multirow{2}{*}{$ B^{-} \to D_0^*(2400)^0 \rho^- $}	& $  \to D^0 \pi^0 $ &  $ 0.61 \sim 1.07\times 10^{-3} $ & \multirow{2}{*}{$ \overline{B}^{0} \to D_0^*(2400)^+ \rho^- $}	& $  \to D^0 \pi^+ $ &  $ 1.18 \sim 2.00 \times 10^{-3} $   \\
			& $ \to D^+ \pi^- $ &  $ 1.23 \sim 2.16  \times 10^{-3} $ & & $ \to D^+ \pi^0 $ &  $ 5.86 \sim 9.93  \times 10^{-4} $  \\

			\multirow{2}{*}{$ B^{-} \to D_0^*(2400)^0 K^{*-} $}	& $  \to D^0 \pi^0 $ &  $ 3.50 \sim 6.15 \times 10^{-5} $ & \multirow{2}{*}{$ \overline{B}^{0} \to D_0^*(2400)^+ K^{*-} $}	& $  \to D^0 \pi^+ $ &  $ 0.68 \sim 1.15 \times 10^{-4} $   \\
			& $ \to D^+ \pi^- $ &  $ 0.71 \sim 1.24  \times 10^{-4} $ & & $ \to D^+ \pi^0 $ &  $ 3.35 \sim 5.71  \times 10^{-5} $  \\

			\bottomrule[1.5pt]
	\end{tabular}}
\end{table}

Then, like the previous section, the non-leptonic productions of $ 2P $ and  $ 3P $ states are also considered and the results are shown in Table \ref{nonleptonic3000} and \ref{nonleptonic3183}, respectively.

\begin{table}[htb!]
	\renewcommand\arraystretch{1.01}
	\caption[results]{Non-leptonic production results of the $ 2P $ state $ D_0^*(3000) $. }
	\label{nonleptonic3000}
	\vspace{0.5em}\centering
	\begin{tabular}{lcc}
		\toprule[1.5pt]
		\multicolumn{1}{c}{Channels} & Widths (GeV) & BR  \\
		\midrule[1pt]
		 $ B^{-} \to D_0^*(3000)^0 \pi^-  $ & $ 0.07\sim 54.7 \times 10^{-18} $ & $ 0.02\sim 13.4 \times 10^{-5} $ \\ $ \overline{B}^{0} \to D_0^*(3000)^+ \pi^-  $ & $ 0.003\sim 55.3 \times 10^{-18} $ & $ 0.0006\sim 12.8 \times 10^{-5} $ \\
		 $ B^{-} \to D_0^*(3000)^0  K^- $ & $ 0.08\sim 37.9\times 10^{-19} $ & $0.02\sim 9.44 \times 10^{-6} $ \\ $ \overline{B}^{0} \to D_0^*(3000)^+  K^- $ & $ 0.02\sim 39.0\times 10^{-19} $ & $ 0.004\sim 9.02 \times 10^{-6} $ \\
		 $ B^{-} \to D_0^*(3000)^0  \rho^- $ & $ 0.04\sim 10.8\times10^{-17} $ & $ 0.10\sim 26.8\times 10^{-5} $ \\ $ \overline{B}^{0} \to D_0^*(3000)^+  \rho^- $ & $ 0.01\sim 11.1\times10^{-17} $ & $ 0.03\sim 25.7 \times 10^{-5} $ \\
		 $ B^{-} \to D_0^*(3000)^0 K^{*-}$ & $ 0.28\sim 59.0\times 10^{-19} $ & $ 0.07\sim 14.7 \times 10^{-6} $ \\ $ \overline{B}^{0} \to D_0^*(3000)^+ K^{*-}$ & $ 0.11\sim 61.2 \times 10^{-19} $ & $ 0.03\sim 14.1 \times 10^{-6} $ \\
		\bottomrule[1.5pt]
	\end{tabular}
\end{table}

Similar to the semi-leptonic occasion, the results of $ 2P $ state $ D_0^*(3000) $ have a large uncertainty , which can also be explained by the node structure of BS wave function. Because the phase spaces of $ 2P $ and $ 3P $ state are close, their branching ratios reach similar magnitude.

\begin{table}[htb!]
	\renewcommand\arraystretch{1.01}
	\caption[results]{Non-leptonic production results of the $ 3P $ state $ D_0^* $. }
	\label{nonleptonic3183}
	\vspace{0.5em}\centering
		\begin{tabular}{lcc}
			\toprule[1.5pt]
			\multicolumn{1}{c}{Channels} & Widths (GeV) & BR   \\
			\midrule[1pt]
			$ B^{-} \to D_0^*(3P)^0 \pi^-  $ & $ 1.87 \sim 5.59 \times 10^{-18} $ & $ 0.47 \sim 1.39 \times 10^{-5} $ \\ $ \overline{B}^{0} \to D_0^*(3P)^+ \pi^-  $ & $ 1.86 \sim 6.00 \times 10^{-18} $ & $ 0.43 \sim 1.39 \times 10^{-5} $  \\
			$ B^{-} \to D_0^*(3P)^0  K^- $ & $ 1.43 \sim 4.02 \times 10^{-19} $ & $ 0.36 \sim 1.00 \times 10^{-6} $ \\ $ \overline{B}^{0} \to D_0^*(3P)^+  K^- $ & $ 1.41 \sim 4.30 \times 10^{-19} $ & $ 0.33 \sim 0.99 \times 10^{-6} $ \\
			$ B^{-} \to D_0^*(3P)^0  \rho^- $ & $ 0.38 \sim 1.03 \times10^{-17} $ & $ 0.96 \sim 2.56 \times 10^{-5} $ \\ $ \overline{B}^{0} \to D_0^*(3P)^+  \rho^- $ & $ 0.38 \sim 1.11 \times10^{-17} $ & $ 0.89 \sim 2.55 \times 10^{-5} $ \\
			$ B^{-} \to D_0^*(3P)^0 K^{*-}$ & $ 2.10 \sim 5.47  \times 10^{-19} $ & $ 0.52 \sim 1.36 \times 10^{-6} $ \\ $ \overline{B}^{0} \to D_0^*(3P)^+ K^{*-}$ & $ 2.11 \sim 5.89 \times 10^{-19} $ & $ 0.49 \sim 1.36 \times 10^{-6} $ \\
			\bottomrule[1.5pt]
	\end{tabular}
\end{table}

We also draw the branching ratios changing with the mass of $ D_0^* $ in Fig. \ref{nonleptonic}. Like the semi-leptonic production case, the non-leptonic production results of  $ 2P $ state are sensitive to the mass. The curves of $ 1P $ and $ 3P $ states stay relatively stable when the mass values change. There are also minimum points of $ 3P $ states' curves at $ 3.175{\ \rm{GeV}} $, where the maximum cancellation of overlapping integral occurs.

\begin{figure}[htb!]]
	\centering
	\subfigure[Braching ratio vers the mass of $ 1P $ state $D_0^*(2400)^0$]{
		\label{nl2400zero}
		\includegraphics[width=0.45\textwidth]{Br_mass_non-leptonic_D0star2400_0_May_2-eps-converted-to.pdf}}
	\hspace{0.2cm}
	\subfigure[Braching ratio vers the mass of $ 1P $ state $D_0^*(2400)^+$]{
		\label{nl2400plusl}
		\includegraphics[width=0.45\textwidth]{Br_mass_non-leptonic_D0star2400_plus_May_2-eps-converted-to.pdf}}

	\subfigure[Braching ratio vers the mass of $ 2P $ state $D_J^*(3000)^0$]{
		\label{nl3000zero}
		\includegraphics[width=0.45\textwidth]{Br_mass_non-leptonic_D0star3000_0_May_2-eps-converted-to.pdf}}
	\hspace{0.2cm}
	\subfigure[Braching ratio vers the mass of $ 2P $ state $D_J^*(3000)^+$]{
		\label{nl3000plus}
		\includegraphics[width=0.45\textwidth]{Br_mass_non-leptonic_D0star3000_plus_May_2-eps-converted-to.pdf}}

	\subfigure[Braching ratio vers the mass of $ 3P $ state $D_0^*(3183)^0$]{
		\label{nl33183zero}
		\includegraphics[width=0.45\textwidth]{Br_mass_non-leptonic_D0star3183_0_May_2-eps-converted-to.pdf}}
	\hspace{0.2cm}
	\subfigure[Braching ratio vers the mass of $ 3P $ state $D_0^*(3183)^+$]{
		\label{nl3183plus}
		\includegraphics[width=0.45\textwidth]{Br_mass_non-leptonic_D0star3183_plus_May_2-eps-converted-to.pdf}}
	\caption{The branching ratios of   non-leptonic production of $ D_0^* $.}
	\label{nonleptonic}
\end{figure}

\section{SUMMARY}

Based on the instantaneous BS framework, we calculate the semi-leptonic and non-leptonic productions of several excited $ D_0^* $ states from $ B$ mesons. For $ 1P $ state $ D_0^*(2400) $, the branching ratios of $ B \to D_0^* e^- \overline{\nu}_{e} $ are in the order of $ 10^{-3} $, which is consistent with the results of present experiments and other models. For non-leptonic channels, the experiments didn't get quite consistent results while our calculating consists with parts of present experimental results.
For $ 2P $ states $ D_J^*(3000) $, we get suppressed branching ratios in the order of $ 10^{-5} \sim 10^{-6}$ and large uncertainty in both semi-leptonic and non-leptonic channels. The cancellation in overlapping integral, which is caused by its one-nodes structure of the BS wave functions, could explain the abnormal results. For $ 3P $ states, their ratios are of the same order of magnitude as $ 2P $ states' results because they have similar phase spaces. The two-nodes structure of $ 3P $ states wave functions makes the cancellation smaller than that of $ 2P $ states and get the minimum branching ratios if their masses are around 3.175 GeV.
Our work could give some inspiration to future experiment and we expect more attention on these production processes of the orbitally excited $ D $ mesons.

\section*{ACKNOWLEDGEMENTS}
This work was supported in part by the National Natural Science
Foundation of China (NSFC) under Grant No.~11575048, No.~11405037, No.~11505039.
We thank the HPC Studio at Physics Department of Harbin Institute of Technology for access to computing resources through INSPUR-HPC@PHY.HIT.


\section{Introduction}
In the last fifty years, nucleon knockout reactions induced by intermediate energy protons of the type $(p,pN)$ in normal (forward) kinematics have been one of the most successful tools for studying the single-particle nature of stable nuclei \cite{Jac66,Jac73,Kit85,Wak17}. These reactions are sometimes referred to as quasifree scattering due to the dominance of scattering processes between the proton and the knocked-out nucleon. The most widely used theoretical approach to analyze these reactions is the distorted-wave impulse approximation (DWIA) method of Chant and Roos \cite{Cha77,Cha83}. A detailed review of this method including its uncertainties and applications to existing forward-kinematics $(p,pN)$ data has been reported in Ref.~\cite{Wak17}.

With the availability of radioactive beams at energies up to 450 MeV/u, there has been a renewed interest to elucidate single-particle properties of unstable nuclei by using nucleon knockout reactions on the hydrogen target in inverse kinematics. The DWIA method has first been theoretically applied to inverse kinematics $(p,pN)$ reactions with the scattering wave functions treated with the eikonal approximation \cite{Aum13,Oga15} and later in (standard) partial-wave expansion form \cite{Yos18}. The latter method was then successfully applied in several experimental studies \cite{Oli17,Ele19,Tan19,Chen19,Kaw18} carried out at RIKEN, Japan. However, an extensive investigation of the sensitivity of calculated cross sections to various choices of inputs and corrections (as in Ref.~\cite{Wak17}) was not done in these studies.

Besides DWIA, other theoretical methods such as the three-body Faddeev equation in the Alt-Grassberger-Sandhas formulation (FAGS) \cite{Cre08,Cre19} and the transfer-to-the-continuum (TC) method \cite{Mor15} have also been used to analyze the $(p,pN)$ experimental data measured in inverse kinematics. Although there is a good consistency between DWIA and TC formalism for the specific $^{15}$C($p$,$pn$)$^{14}$C reaction \cite{Yos18}, a systematic comparison between the DWIA (eikonal and partial-wave form), TC, and FAGS analyses on the actual experimental data is essential to determine the range of applicability for $(p,pN)$ reactions.

Among the possible uses of proton-induced nucleon knockout reactions, the study of quenching single-particle strength and its proton-neutron asymmetry dependence is one of the most important subjects. A reduction of 30\%--40\% with respect to the independent-particle model (IPM) limit in spectroscopic factors (SFs) deduced from $(e,e'p)$ experiments was observed at NIKHEF \cite{Lap93}. This quenching of the SF, quantified as the reduction factor $R_s$, is due to the lack of short-range (including tensor) and long-range correlations in the IPM and standard shell model (SM) calculations \cite{Pan97,Dic04,Sick07,Bar09,Jen11,Cip15}.

Systematic analysis of nucleon removal reactions on light composite targets ($^{9}$Be and $^{12}$C) \cite{Gade08,Tos14} suggests a strong dependence on the proton-neutron asymmetry defined as ${\Delta S=S_p-S_n}$ $(\Delta S=S_n-S_p)$ for proton (neutron) removal, where $S_p$ $(S_n)$ is the proton (neutron) separation energies. However, such strong dependence cannot be observed in systematic nucleon-transfer studies \cite{Lee06,Kay13,Fla13,Fla18,Xu19}. A weak dependence of $R_s$ on $\Delta S$ is also supported by \textit{ab initio} coupled-cluster (CC) \cite{Jen11} and self-consistent Green's function (SCGF) \cite{Cip15} calculations.

Very recently, a series of $(p,pN)$ measurements for carbon-, nitrogen-, and oxygen-isotope beam with incident energy range of 300--450 MeV/u in inverse kinematics was performed at the R$^3$B/LAND setup at GSI Helmholtzzentrum f\"{u}r Schwerionenforschung in Darmstadt, Germany \cite{Pan16,Atar18,Fer18,Holl19}. While the measurements in Ref.~\cite{Pan16} are exclusive, those in Refs.~\cite{Atar18,Fer18,Holl19} provide the semi-inclusive cross sections and momentum distributions, in which various bound states of the residual nucleus are summed. The eikonal DWIA \cite{Pan16,Atar18,Holl19}, FAGS \cite{Fer18}, and TC \cite{Gom18} methods were applied to these experimental data and gave a generally similar conclusion that the reduction factors depend very weakly on the proton-neutron asymmetry, in contradiction with the much steeper asymmetry found in nucleon removal analysis. However, as was noted in Ref.~\cite{Gom18}, the $(p,pN)$ data analyzed by different reaction models exhibit some discrepancies due to choices of inputs and nonrelativistic treatments. These results also slightly underestimate the magnitude of $R_s$ given by \textit{ab initio} calculations \cite{Jen11,Cip15}. To clarify the inconsistencies between these models and give a reliable evaluation of the experimental data, it is of great interest to perform a careful DWIA analysis on the reduction factors with the GSI $(p,pN)$ data \cite{Pan16,Atar18,Fer18,Holl19}.

The content of this paper is as follows. In Sec.~\ref{secformalism} the formulation of the standard partial-wave DWIA formalism is given. In Sec.~\ref{secresult} the cross sections and reduction factors are calculated for all published GSI data. The impact of several corrections in the DWIA framework on these observables are investigated. The transverse momentum distributions for the specific $^{12}$C($p$,$2p$)$^{11}$B case is also discussed. Finally, the summary is given in Sec.~\ref{secsum}.

\section{Formalism}
\label{secformalism}
The A($p$,$pN$)B knockout reaction is analyzed with the same partial-wave DWIA framework as in Ref.~\cite{Yos18}. Observables with superscript A are evaluated in the A-rest frame while those without the superscript are in the three-body center-of-mass (c.m.) frame, also called the G frame. The transition amplitude for the A($p$,$pN$)B reaction is given by
\begin{align}
T_{\bm{K}_0\bm{K}_1\bm{K}_2}^{nljm}
&=
\Braket{
\chi_{1,\bm{K}_1}^{(-)}\chi_{2,\bm{K}_2}^{(-)}
|t_{pN}|
\chi_{0,\bm{K}_0}^{(+)}\varphi^{nljm}
},
\label{eqtrans}
\end{align}
where $\chi_{i,\bm{K}_i} (i=0,1,2)$ are the distorted scattering wave functions of the $p$-A, $p$-B, and $N$-B systems, respectively. $\bm{K}_i$ is the momentum of particle $i$ in the G frame. As has been shown in Eqs.(3.8)-(3.15) of Ref.~\cite{Wak17}, when the kinematic coupling term in the exit channel Hamiltonian is approximated, the three-body scattering wave function can be separated into two two-body distorted-wave functions in Eq.~(1), where $\bm{K}_i$ can now be interpreted as the relative momentum. The superscripts $(+)$ and $(-)$ specify the outgoing and the incoming boundary conditions of these scattering waves, respectively. The relative single-particle wave function of the $N$-B system bound inside A is denoted as $\varphi^{nljm}$ where $n$, $l$, $j$, and $m$ are the principal quantum number, the orbital angular momentum, the total angular momentum, and its third component, respectively. $t_{pN}$ is a transition operator for the $p$-$N$ scattering, which is sometimes called an effective interaction \cite{FL85}. The absolute square of its matrix element is proportional to the $pN$ elastic-scattering cross section ${d\sigma_{pN}}/{d\Omega_{pN}}$.

Following the same theoretical treatment as in Ref.~\cite{Yos18}, upon disregarding the spin-orbit distortion, the momentum distribution (MD) is given by
\begin{align}
\frac{d\sigma}{d\bm{K}_\mathrm{B}^\mathrm{A}}
=
&C_0 \int d\bm{K}_1^\mathrm{A}d\bm{K}_2^\mathrm{A}
\delta(E_f^\mathrm{A}-E_i^\mathrm{A})
\delta(\bm{K}_f^\mathrm{A}-\bm{K}_i^\mathrm{A}) \nonumber \\
&\times \frac{E_1 E_2 E_\mathrm{B}}{E_1^\mathrm{A} E_2^\mathrm{A} E_\mathrm{B}^\mathrm{A}}
\frac{d\sigma_{pN}}{d\Omega_{pN}}
\sum_{m}(2\pi)^2
\lvert
\bar{T}_{\bm{K}_0\bm{K}_1\bm{K}_2}^{nljm}
\rvert ^2, \label{eq.md}
\end{align}
where
\begin{align}
C_0
&\equiv
\frac{E_0^\mathrm{A}}{(\hbar c)^2 K_0^\mathrm{A}}
\frac{f_{pN}}{(2l+1)}\frac{\hbar^4}{(2\pi)^3 \mu_{pN}^2}.
\end{align}
The factor $f_{pN}$ equals to 1 for $(p,pn)$ and $1/2$ for $(p,2p)$ reactions. $E_i$ and $E^{\rm A}_i$ are the total (relativistic) energies of particle $i$ in the G and A-rest frames, respectively.

The G-frame $pN$ scattering cross section in Eq.~(\ref{eq.md}) is related to the one in the two-nucleon c.m. frame, which we term the t frame, through
\begin{align}
\frac{d\sigma_{pN}}{d\Omega_{pN}}=
\eta^2\frac{d\sigma_{pN}^\textrm{\,t}}{d\Omega_{pN}^\textrm{\,t}},
\end{align}
where
\begin{align}
\eta=\left(\frac{E_1^\mathrm{t} E_2^\mathrm{t} E_0^\mathrm{t} E_N^\mathrm{t}}{E_1 E_2 E_0 E_N}\right)^{1/2}, \label{eq.mol}
\end{align}
is the M\o ller factor required for the transformation of $t_{pN}$ from the t frame to the G frame in relativistic kinematics \cite{Mol45,Ker59}. The total energy of the struck nucleon $E_N$ is determined by the momentum conservation of the two colliding nucleons as \cite{Yos18}
\begin{align}
E_N=\dfrac{\hbar^2}{2\mu_{NB}}\left[\bm{K}_1+\bm{K}_2-\dfrac{(A+1)}{A}\bm{K}_0\right]^2.
\end{align}

The reduced transition amplitude in Eq.~(\ref{eq.md}) is given by
\begin{align}
\bar{T}_{\bm{K}_0\bm{K}_1\bm{K}_2}^{nljm}
&=
\int d\bm{R}\,\chi_{1,\bm{K}_1}^{*(-)}(\bm{R})\,
\chi_{2,\bm{K}_2}^{*(-)}(\bm{R})\,
\chi_{0,\bm{K}_0}^{(+)}(\bm{R})\,
\nonumber \\
&\times
\varphi^{nljm}(\bm{R}) e^{-i\bm{K}_0\cdot\bm{R}/A}.
\end{align}
For the bound-state wave function $\varphi^{nljm}$ generated by a local potential, the effect of nonlocality is taken into accounted in the interior region by multiplying $\varphi^{nljm}$ by the Perey factor \cite{Per63}
\begin{align}
F_\textrm{PR}(R)=C_\textrm{PR}\left[1-\dfrac{\mu_{N\mathrm{B}}}{2\hbar^2}\beta^2 V_{N\mathrm{B}}(R)\right]^{-1/2},
\end{align}
where the nonlocality range $\beta=0.85$ fm for nucleon \cite{Per62}, and $V_{N\mathrm{B}}$ is the single-particle binding potential. The factor $C_\textrm{PR}$ is chosen so that the modified bound-state wave function is normalized to unity. Similarly, for a scattering wave function obtained from a Dirac phenomenology optical potential (OP), the relativistic velocity-dependent term modifies the wave function by what we refer to as the Darwin factor \cite{Arn81,Ham90}
\begin{align}
F_\textrm{DW}(R)=\left[\dfrac{E_i+U_S(R)-U_V(R)}{E_i}\right]^{1/2},
\end{align}
where $U_S$ and $U_V$ are the scalar and vector potentials in the Dirac equation, respectively. This Darwin factor is regarded as a kind of nonlocality correction and has been well known to be very important in order to fully take into account relativistic effect in $(e,e'p)$ reactions \cite{Udi95}.

The cylindrical transverse momentum distributions (TMD) are obtained from the MD as
\begin{align}
\frac{d\sigma}{dK_{\mathrm{B}b}^\mathrm{A}}
&=
2\pi\int dK_{\mathrm{B}z}^\mathrm{A} K_{\mathrm{B}b}^\mathrm{A}
\frac{d\sigma}{d\bm{K}_\mathrm{B}^\mathrm{A}}.
\end{align}
The integrated single-particle cross section is then calculated from the TMD as
\begin{align}
\sigma_\textrm{sp}
&=
\int \frac{d\sigma}{dK_{\mathrm{B}b}^\mathrm{A}}dK_{\mathrm{B}b}^\mathrm{A}. \label{eq.sp}
\end{align}

\section{Results and discussion \label{secresult}}
In this section, we applied the DWIA framework described in the preceding section to the GSI $(p,pN)$ data \cite{Pan16,Atar18,Fer18,Holl19} to evaluate the single-particle cross section and deduce the reduction factor. The impact of nonlocality corrections, the M\o ller factor, and energy dependence of final-state OP on the $R_s$ are clarified.

\subsection{Numerical inputs}
\label{subsecinput}

We perform DWIA calculations for all 18 published $(p,pN)$ cases of the R$^3$B/LAND setup \cite{Pan16,Atar18,Fer18,Holl19}. In addition, we also perform calculations with several other choices of inputs to estimate the theoretical uncertainty on reduction factors. All calculations in this work adopt the relativistic treatment of the kinematics, which is essential to reproduce the correct MD \cite{Yos18}. The nonlocality correction is taken into account in our calculation through the use of the Perey factor for the single-particle bound state and Darwin factor for distorted waves.

For the distorting potential of the $p$-A, $p$-B, and $N$-B systems, we use the EDAD2 parameter set of the Dirac phenomenology \cite{Coo93}. Calculations using the EDAD1, EDAD3 \cite{Coo93}, and the ``democratic'' EDAD \cite{Coo09} Dirac OP sets give a difference of 10\% at most, which is consistent with the one observed in normal kinematics \cite{Wak17}.

The single-particle wave function of the struck nucleon is obtained from a Woods-Saxon potential with central and spin-orbit components defined in the same manner as in Refs.~\cite{Gade08,Gom18}. For both components, a diffuseness $a=0.7$~fm is used for all the cases. The radius parameter is adjusted following the prescription $\langle r^2 \rangle=[A/(A-1)]\langle r^2 \rangle_\textrm{HF}$, where $\langle r^2 \rangle_\textrm{HF}$ is the single-particle mean square radius of the Hartree-Fock calculation with the Skyrme SkX interaction \cite{Bro98}. The depth $V_\textrm{so}=6$ MeV is fixed for the spin-orbit term while the central one is adjusted to reproduce the experimental separation energies. This choice of binding potential gives a difference in the $(p,pN)$ cross section within 10\% compared with the one used in the $(e,e'p)$ analysis \cite{Kra01,Wak17}. Based on this result and the investigations of Refs.~\cite{Gade08,Gom18} on the ambiguity of different effective nucleon-nucleon (\textit{NN}) interactions used in the Hartree-Fock calculation, we adopt an uncertainty of 10\% for single-particle wave functions.

For the $pN$ elementary cross section, we employ the one generated by the $t$-matrix parametrization of Franey and Love \cite{FL85} with a final-energy prescription, which has been suggested to be the best approximation for the half-off-shell amplitude \cite{Red70}. Different choices of the on-shell approximation such as initial-energy and average-energy prescriptions give an uncertainty of 2\% for $(p,2p)$ and 8\% for $(p,pn)$ processes. The large discrepancies in neutron knockout processes are due to the asymmetric shape of the $pn$ scattering angular distribution that makes it more sensitive to the choice of energy prescription. Note that other choices of the \textit{NN} cross section such as those from the Reid93 potential \cite{Reid93} in TC \cite{Gom18}, the CD-Bonn potential \cite{Mac01} in FAGS \cite{Fer18}, or from experimental database fitting \cite{Arn03} in the eikonal DWIA \cite{Aum13} give essentially the same \textit{NN} cross sections up to 350 MeV.

The theoretical spectroscopic factor for each state of the bound residual core is the same as in Refs.~\cite{Gom18,Holl19}, which is computed by the standard SM with the WBT interaction \cite{War92} and includes the c.m. correction \cite{Die74}. Since many of the considered nuclei are weakly bound, the use of SF calculated from the SM, compared with the IPM limit, provides a more proper description of single-particle-strength fragmentation near the Fermi level.

\subsection{Reduction factors} \label{sec.redfac}

\begin{table*}[htbp]

	\caption{Experimental cross sections $\sigma_\textrm{exp}$ \cite{Atar18,Pan16,Fer18,Holl19}, calculated ones $\sigma_\textrm{th}$, and reduction factors $R_s$. See the text for details.}
	\begin{ruledtabular}
		\begin{tabular}{lcddc}
			Reaction & \multicolumn{1}{l}{$E_\textrm{beam}$ (MeV/u)} & \multicolumn{1}{c}{$\sigma_\textrm{th}$ (mb)} & \multicolumn{1}{c}{$\sigma_\textrm{exp}$ (mb)} & \multicolumn{1}{c}{$R_s$} \\ \hline
			$^{10}$C($p$,$pn$)$^{9}$C  & 386   & 12.95 & 16.3(22)[14]      & 1.26(29)   \\
			$^{11}$C($p$,$2p$)$^{10}$B & 325   & 15.68 & 18.2(9)[10]       & 1.16(19)   \\
			$^{11}$C($p$,$pn$)$^{10}$C & 325   & 14.07 & 17.0(15)[21]      & 1.21(27)   \\
			$^{12}$C($p$,$2p$)$^{11}$B & 398   & 22.04 & 19.2(18)[12]      & 0.87(16)   \\
			$^{12}$C($p$,$pn$)$^{11}$C & 398   & 27.43 & 30.0(32)[27]      & 1.09(23)   \\
			$^{13}$O($p$,$2p$)$^{12}$N & 401   & 5.77  & 5.78(91)[37]      & 1.00(22)   \\
			$^{14}$O($p$,$2p$)$^{13}$N & 351   & 13.28 & 10.23(80)[65]     & 0.77(13)   \\
			$^{15}$O($p$,$2p$)$^{14}$N & 310   & 18.07 & 18.92(182)[120]   & 1.05(19)   \\
			$^{16}$O($p$,$2p$)$^{15}$N & 451   & 27.78 & 26.84(90)[170]    & 0.97(15)   \\
			$^{17}$O($p$,$2p$)$^{16}$N & 406   & 9.16  & 7.90(26)[50]      & 0.86(14)   \\
			$^{18}$O($p$,$2p$)$^{17}$N & 368   & 20.01 & 17.80(104)[113]   & 0.89(15)   \\
			$^{21}$O($p$,$2p$)$^{20}$N & 449   & 5.58  & 5.31(23)[34]      & 0.95(15)   \\
			$^{21}$N($p$,$2p$)$^{20}$C & 417   & 3.25  & 2.27(34)          & 0.70(14)   \\
			$^{21}$N($p$,$pn$)$^{20}$N & 417   & 38.87 & 48.52(404)        & 1.25(23)   \\
			$^{22}$O($p$,$2p$)$^{21}$N & 414   & 6.90  & 6.01(41)          & 0.87(14)   \\
			$^{22}$O($p$,$pn$)$^{21}$O & 414   & 36.24 & 39.24(234)        & 1.08(19)   \\
			$^{23}$O($p$,$2p$)$^{22}$N & 445   & 4.97  & 4.93(96)          & 0.99(24)   \\
			$^{23}$O($p$,$pn$)$^{22}$O & 445   & 50.05 & 54.0(108)         & 1.08(28)   \\
		\end{tabular}%
		\label{tab.result}%
	\end{ruledtabular}
\end{table*}%

The results of our DWIA calculations are presented in Table ~\ref{tab.result}. Because of the semi-inclusive nature of the data concerned, the calculation results are the sum of the cross sections corresponding to several bound-state configurations of the residual nucleus. The beam energy in the middle of the target is shown in the second column. The third column indicates the theoretical cross sections $\sigma_\textrm{th}=\sum C^2S \times \sigma_\textrm{sp}$, where $\sigma_\textrm{sp}$ is the single-particle cross section of a specific configuration. The fourth column shows the experimental cross section, with statistical (round brackets) and systematic (square brackets) uncertainties, corresponding to the beam of $^{10-12}$C \cite{Pan16,Holl19}, $^{13-21}$O \cite{Atar18} and $^{21}$N,$^{22,23}$O \cite{Fer18}. The reduction factor $R_s=\sigma_\textrm{exp}/\sigma_\textrm{th}$ is given in the last column.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.48\textwidth]{Main-theory.eps}
	\caption{
		Reduction factors deduced from $(p,pN)$ reactions using DWIA compared with the reduced SFs calculated from CC \cite{Jen11} (black circles) and SCGF \cite{Cip15,Atar18} with different interactions (blue and green triangles). The black (red) texts in the header label the isotope whose proton (neutron) is knocked out. See the text for details.
	}
	\label{fig.main-theory}
\end{figure}

The reduction factors as a function of the proton-neutron asymmetry $\Delta S$ is shown in Fig.~\ref{fig.main-theory}. The value calculated from the present DWIA analysis is indicated by red squares with error bars propagated from the experimental uncertainties reported in Refs.~\cite{Pan16,Atar18,Fer18,Holl19} and the theoretical uncertainties estimated in Sec.~\ref{secformalism}, which is about 14\%--16\%. The total relative uncertainties for the extracted reduction factor are ranging from 15\% to 25\%, which is similar to those reported in the recent systematic $(p,d)$ analysis \cite{Xu19}. In general, a trend of the reduction factor about 0.9--1.0 with a very weak asymmetry dependence is observed. Properly taking the uncertainties into account, the reduction factors $R_s=0.87(16)$ for $^{12}$C and $R_s=0.97(15)$ for $^{16}$O are overestimated compared with the $(e,e'p)$ \cite{Lap93} results of $R_s=0.57(6)$ and $R_s=0.65(5)$, respectively. By performing a linear function fitting, the $\Delta S$ dependence of $R_s$ is obtained as $R_s=0.947(36)-2.6(27)\times 10^{-3}\Delta S$ with a reduced $\chi^2/N$ of 0.74. As discussed in Sec.~\ref{sec.reaction}, the close-to-unity reduction factor does not necessarily mean that the quenching effect observed in $(p,pN)$ reactions is weak but rather indicates a fundamental problem in current reaction models.

The reduction factors from DWIA are compared with the reduced SFs, which is the ratio of the SF to the IPM limit, from \textit{ab initio} self-consistent Green’s function (SCGF) \cite{Atar18,Cip15} and coupled-cluster (CC) \cite{Jen11} models. We note that all of these values are presented as a function of \textit{experimental} $\Delta S$. The CC calculation \cite{Jen11} uses the chiral \textit{NN} interaction at next-to-next-to-next-to-leading order (N3LO) \cite{Ent03} with the cutoff at $\Lambda_\textit{NN}=500$ MeV while the SCGF of Ref.~\cite{Cip15} uses the same \textit{NN} interaction in addition to an NNLO three-nucleon force \cite{Roth12} with $\Lambda_{3N}=400$ MeV. We compare also with the result of the more recent SCGF calculation in Ref.~\cite{Atar18} based on the NNLO-sat \cite{Eks15} for the \textit{NN} interaction, which is more optimized for the mass region in that study. The present DWIA calculation shows a reasonable agreement with \textit{ab initio} results, especially about the weak dependence of the trend on proton-neutron asymmetry.

\subsection{Comparison with $\bm{(e,e'p)}$ results and other reaction models} \label{sec.reaction}
Figure \ref{fig.main-reaction} shows the reduction factors from the present DWIA calculation in comparison with those analyzed by other theoretical reaction models. The slope of the reduction factor obtained by the present analysis is in excellent agreement with the value deduced by TC model with fixed-energy EDAD2 OP \cite{Gom18}. It is also consistent (within the uncertainty range) with the slopes observed in Refs.~\cite{Atar18,Fla18}.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.48\textwidth]{Main-reaction.eps}
	\caption{
		Same as Fig.~\ref{fig.main-theory} but compared with other $(p,pN)$ analysis performed with the eikonal DWIA \cite{Atar18} (green triangles), FAGS \cite{Fer18} (pink triangles), TC \cite{Gom18} (blue dashed line) for the GSI data, and similar partial-wave DWIA for the RIKEN/RCNP data \cite{Kaw18} (black circles). The results from Refs.~\cite{Gom18,Kaw18} are those performed with Dirac OP. The blue shaded band indicates the trend observed in the analysis of nucleon removal reaction with composite targets \cite{Gade08,Tos14}.
	}
	\label{fig.main-reaction}
\end{figure}

The values of the reduction factor given by DWIA overestimate those reported by $(e,e'p)$ analysis \cite{Lap93} and other reaction models such as TC \cite{Gom18} and the eikonal DWIA \cite{Atar18}. This means the present DWIA calculation gives a smaller cross section compared with these. However, the better agreement between the $(e,e'p)$ results and those given by other $(p,pN)$ analyses of the same GSI data will be, in fact, caused by the lack of several crucial corrections in the reaction models used in these analyses. We will return to this point below.

The overshooting of the $(e,e'p)$ results means that the DWIA cross section is smaller than the observed value. This indicates that some contributions from higher-order processes such as multistep scattering or channel coupling are included in the GSI data, especially in the large recoil momentum region. In fact, a smaller reduction factor around 0.7 is observed from the analysis of $(p,2p)$ data measured at RIKEN/RCNP \cite{Kaw18} using the same DWIA framework as in this study. An important feature of the measurement at RIKEN/RCNP is a very
constrained kinematics corresponding to the quasifree condition. This supports the conclusion that the lack of higher-order effects will be the main reason for the underestimation of the GSI data with the current DWIA calculation.

The discrepancies between the current DWIA results and those using other reaction models can be explained by several factors. First, the nonlocality corrections in the single-particle and scattering wave functions are not presented in the TC and eikonal DWIA models. Second, the TC calculation \cite{Gom18} uses the energy-independent optical potentials evaluated at $E_\textrm{beam}/2$. On the other hand, the large discrepancy with the FAGS method \cite{Fer18} is due to two reasons pointed out in Ref.~\cite{Gom18}, the lack of relativistic kinematics in the FAGS framework as shown in \cite{Yos18} and different choices of optical potentials used in Ref.~\cite{Fer18}. The latter has been confirmed by recent FAGS calculation with a more proper OP for the energy range considered, where the cross section has been reduced by almost 30\% \cite{Cre19}.

We further investigate the possible source of discrepancies between different reaction models used in $(p,pN)$ studies. Impacts of the lack of nonlocality corrections and energy-dependent potentials on the reduction factors are illustrated in Fig.~\ref{fig.deltaNL}. The relative difference $\Delta R$ is evaluated with respect to the ``reference'' DWIA result shown in Sec.~\ref{sec.redfac}. As seen from the lines in Fig.~\ref{fig.deltaNL}, these effects are strongly associated with the separation energy and affect the reduction factor though the cross section in the opposite ways. For the weakly bound nucleus, the single-particle wave function is more extended and since the nonlocality correction only affects the interior of the wave function, its effect gets weaker with decreasing separation energy. On the other hand, since the energy-independent potentials of the outgoing nucleons are evaluated at half the incident beam, the deviation from energy-dependent OP is minimized for the small separation energy, where the kinematics most resembles the quasifree condition. Some deviations from both lines are because of the cross section for each case is the sum of single-particle wave functions with different orbitals. In general, by neglecting both the nonlocality correction and the energy dependence of the OP, DWIA calculations are expected to have a similar magnitude to those from TC \cite{Yos18}.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.48\textwidth]{deltaNL.eps}
	\caption{
		The relative difference with respect to the reference DWIA results as a function of knocked-out-nucleon separation energy. The red squares and black circles represent the DWIA calculations with energy-dependent potentials and without nonlocality corrections, respectively.
	}
	\label{fig.deltaNL}
\end{figure}

Finally, the effect of the M\o ller factor is presented in Fig.~\ref{fig.deltaMol}. This factor, shown in Eq.~(\ref{eq.mol}), has a relativistic origin \cite{Mol45,Ker59,Udi95} and therefore is directly related to the incident energy. For the considered energy range of 300--450 MeV, the M\o ller factor can contribute about 18\%--26\% to the reduction factor, which also explains the small magnitude of the nonrelativistic FAGS model. We note that, although the TC model does not explicitly include the M\o ller factor, its consistency with the DWIA model as found in \cite{Yos18} suggests that the M\o ller-factor effect has been implicitly included in the relativistic treatment of TC \cite{Mor15}.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.48\textwidth]{deltaMol.eps}
	\caption{
		The relative difference between the reference DWIA results and those without the M\o ller factor as a function of beam energy.
	}
	\label{fig.deltaMol}
\end{figure}

Although the lack of nonlocality correction, energy dependence in OP, and the M\o ller factor may cancel with the lack of higher-order effects in some models and give a better agreement with $(e,e'p)$ results, they also obscure the true nature of the problem. As we have seen, the effects of these corrections are highly dependent on the separation energy of the single-particle wave function and beam energy, so more proper inclusion of these corrections will be very essential for future knockout studies.

\subsection{Transverse momentum distribution}

In inverse kinematics nucleon knockout reactions, a comparison between the measured and calculated momentum distribution of the residual nucleus can reveal a lot of information about the reaction mechanism as well as the validity of the theoretical model. We consider the cylindrical TMD of $^{12}$C($p$,$2p$)$^{11}$B \cite{Pan16}, which is some of the highest resolution data from the R$^3$B/LAND experiments.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.46\textwidth]{12C.eps}
	\caption{
		Cylindrical transverse momentum distribution of the $^{12}$C($p$,$2p$)$^{11}$B reaction. The experimental data are taken from Ref.~\cite{Pan16}. The DWIA calculations with and without the nonlocality corrections are presented as the black solid and red dashed lines, respectively. The blue dotted line represents the result with nonlocality corrections scaled down by a factor of 0.655.
	}
	\label{fig.12C}
\end{figure}

The possible absence of higher-order effects is pronounced in the comparison between the DWIA result (black solid line) and the experimental data \cite{Pan16} shown in Fig.~\ref{fig.12C}. The calculation results are the sum of the theoretical TMD corresponding to the bound states of the $^{11}$B core multiplied by the reduction factor $R_s$ of 0.87. One sees that the calculation result including the nonlocality corrections underestimates the experimental data in the momentum region about 150--300 MeV/$c$. We note that this undershooting of theoretical prediction is also observed in the FAGS calculation for the same data using OP fit to elastic-scattering data at the proper energy range \cite{Cre19}. Moreover, a similar discrepancy with the data in the high-momentum region found in the Cartesian TMD calculated with TC \cite{Gom18} implies that some contributions are also not properly accounted for in other reaction models. Because of the success of the DWIA framework for the same system in the quasifree condition \cite{Wak17}, the discrepancies reinforce our claim that they originate from the higher-order effects that take place when the recoil momentum becomes high.

The red dashed line represents the DWIA result without nonlocality correction, i.e., neglecting the Perey and Darwin factors. The corresponding $R_s$ is 0.77. Although this line may seem to improve the result, we emphasize that such a prescription is inappropriate and inconsistent with the known properties of $(p,pN)$ scattering as reviewed in
Ref.~\cite{Wak17}. The blue dotted line is the DWIA result with nonlocality corrections multiplied by 0.66, which is the ratio of the $R_\textrm{s}$ in Table \ref{tab.result} to the one determined by $(e,e’p)$ \cite{Lap93}. If the $(e,e’p)$ $R_s$ is correct, the blue dotted line should agree with the experimental data. Therefore, it is clear that some additional contributions are necessary to reproduce the data.

\section{Summary}
\label{secsum}
We have performed an analysis on all the published data to date for the $(p,pN)$ reaction in inverse kinematics by the R$^3$B collaboration by using the standard partial-wave DWIA formalism. Our study focuses on evaluating the source of ambiguity in DWIA calculations of the reduction factor and investigating the discrepancies between various reaction models currently used for inverse kinematics $(p,pN)$ data.

Our study suggests a very weak dependence of the reduction factor on the proton-neutron asymmetry $\Delta S$. This result is consistent with previous analyses on $(p,pN)$ using TC, eikonal and partial-wave DWIA \cite{Gom18,Atar18,Kaw18,Holl19} and \textit{ab initio} calculations \cite{Cip15,Atar18,Jen11}.

More importantly, the present study suggests that the lack of a proper treatment of higher-order effects may considerably affect the cross section at kinematics far from the quasifree condition. That effect, which does not manifest in experiments performed around the recoilless condition, becomes more crucial in the semi-inclusive type of integrated cross-section measurements like those carried out at the GSI R$^3$B/LAND setup. However, it is mostly hindered by the lack of essential corrections such as nonlocality, relativistic, M\o ller factor, and energy-dependent OP as well as the considerably large uncertainty from various choices of distorting potentials and single-particle wave functions. The proper inclusion of these corrections in the future $(p,pN)$-reactions analyses will be required. Furthermore, unless measurements with restricted kinematics similar to that in Ref.~\cite{Kaw18} are performed, a reaction model that takes into account higher-order processes will be necessary.

Recently, a consistent description of the nonlocality in bound and scattering wave functions for $(e,e'p)$ has been done with the nonlocal dispersive optical model approach \cite{Atk18}. The incorporation of such treatment in the DWIA framework for proton-induced nucleon knockout reactions is in progress and will be reported elsewhere.

\begin{acknowledgments}
We thank Mario G\'{o}mez-Ramos and Matthias Holl for providing essential details of their calculations. We also thank Tetsuo Noro for valuable discussions on the DWIA calculation. One of the authors (N.T.T.P) would like to thank the Research Center for Nuclear Physics, Osaka University for their hospitality during his stay in which this work was initiated. He also wants to thank the Vietnam MOST for its support through the Physics Development Program Grant No. {\DJ}T{\DJ}LCN.25/18. The computation was carried out with the computer facilities at the Research Center for Nuclear Physics, Osaka University. This work was supported in part by Grants-in-Aid of the Japan Society for the Promotion of Science (Grants No. JP16K05352), and the RCNP Young Foreign Scientist Promotion Program.
\end{acknowledgments}


\section{Introduction}

Given positive integers $m, n$ and a real number $0 < \gamma < mn$, it is define the multilinear fractional operator $\mathcal{I}_{\gamma}$ by

\[
\mathcal{I}_{\gamma} (f_1, ..., f_m)(x) = \int_{(\mathbb{R}^{n})^{m}} \frac{f_1(y_1) \cdot \cdot \cdot f_m(y_m)}{(|x - y_1| + \cdot \cdot \cdot +|x-y_m|)^{nm - \gamma}} dy_1 \cdot \cdot \cdot dy_m, \,\,\,\,\, x \in \mathbb{R}^{n}.
\]

\

Y. Lin and S. Lu in \cite{Lin} proved Hardy space estimates for the multilinear fractional operator $\mathcal{I}_{\gamma}$. More precisely, they proved that if $0 < \gamma < n$, $0 < p_1, . . . , p_m, q \leq 1$, and $q$ such that $\frac{1}{q}= \frac{1}{p_1} + \cdot \cdot \cdot + \frac{1}{p_m} - \frac{\gamma}{n} > 0$, then
\[
\| \mathcal{I}_{\gamma}(f_1, ..., f_m) \|_{L^{q}} \leq C \| f_1 \|_{H^{p_1}} \cdot \cdot \cdot \| f_m \|_{H^{p_m}}.
\]
Recently, D. Cruz-Uribe, K. Moen and H. van Nguyen in \cite{Cruz-Uribe} generalized the result of Lin and Lu to weighted Hardy
spaces on the full range $0 < \gamma < nm$.

The purpose of this note is to give a counter-example to show that the multilinear fractional operator $\mathcal{I}_{\gamma}$ is not bounded from a product of Hardy spaces into a Hardy space. For them, we consider $n=1$, $m=2$, $\gamma = \alpha +1$ with $0 < \alpha < 1$, so
$1 < \gamma < 2$, $2 - \gamma = 1 - \alpha$ and the multilinear fractional operator $\mathcal{I}_{\alpha +1}$ in this case is given by
\[
\mathcal{I}_{\alpha +1}(f_1,f_2)(x) = \iint_{\mathbb{R}^{2}} \frac{f_1(s) f_2(t)}{(|x-s|+|x-t|)^{1-\alpha}} ds dt, \,\,\,\ x \in \mathbb{R}.
\]
We will prove that the operator $\mathcal{I}_{\alpha +1}$ is not bounded from $H^{1}(\mathbb{R}) \times H^{p}(\mathbb{R})$ into $H^{q}(\mathbb{R})$, for $0 < p \leq (\alpha + 1)^{-1}$ and $\frac{1}{q}= \frac{1}{p} - \alpha$.

\

We briefly recall the definition of Hardy space on $\mathbb{R}^{n}$. The Hardy space $H^p(\mathbb{R}^n)$ (for $0 < p < \infty$) consists of tempered distributions $f \in \mathcal{S}'(\mathbb{R}^n)$ such that for some Schwartz function $\varphi$ with $\int \varphi = 1$, the maximal operator
$$
{\displaystyle (\mathcal{M}_{\varphi }f)(x)=\sup _{t>0}|(\varphi_{t}*f)(x)|}$$
is in $L^p(\mathbb{R}^n)$, where $\varphi_t(x):=\frac{1}{t^n}\varphi(\frac{x}{t})$.
In this case we define $\left\Vert  f \right\Vert_{H^{p}}:=\left\Vert \mathcal{M}_\varphi f \right\Vert _{p}$ as the $H^{p}$ ``norm". It can be shown that this definition does not depend on the choice of the function $\varphi$.
For $1 < p < \infty$, it is well known that $H^{p}(\mathbb{R}^{n}) \cong L^{p}(\mathbb{R}^{n})$, $H^{1}(\mathbb{R}^{n}) \subset L^{1}(\mathbb{R}^{n})$ strictly, and for $0 < p < 1$ the spaces $H^{p}(\mathbb{R}^{n})$ and $L^{p}(\mathbb{R}^{n})$ are not comparable.

The following sentences hold in Hardy spaces $H^{p}(\mathbb{R}^{n})$ for $0 < p \leq 1$ (see pp. 128-129 in \cite{St}):

\

(S1) A bounded compactly supported function $f$ belongs to $H^{p}(\mathbb{R}^{n})$ if and only if it satisfies the moment conditions $\int x^{\beta} f(x) dx = 0$ for all $|\beta| \leq n (p^{-1} -1)$.

\

(S2) If $f \in L^{1}(\mathbb{R}^{n}) \cap H^{p}(\mathbb{R}^{n})$ then $\int  x^{\beta} f(x) dx = 0$, whenever $|\beta| \leq n (p^{-1} -1)$ and the function $x^{\beta} f(x)$ is in $L^{1}(\mathbb{R}^{n})$.

\

To obtain our result we will compute explicitly in Section 2 the Fourier transform of the kernel $\left(  |x-s|+|x-t| \right)^{\alpha - 1}$ in the $x$ variable, this allows us to get the following identity
\[
\int_{\mathbb{R}} \mathcal{I}_{\alpha +1}(a_1, a_2)(x) dx = \frac{\alpha - 1}{\alpha} \iint_{\mathbb{R}^{2}} a_{1}(s) \, a_2(t) \, |t-s|^{\alpha} \, ds dt,
\]
valid for bounded functions $a_1$ and $a_2$ having compact support with $\int a_1 =0$ or $\int a_2 =0$. Then, from (S2), the counter-example will follow to consider $a_1 \in H^{1}(\mathbb{R})$ and $a_2 \in H^{p}(\mathbb{R})$ such that $\iint_{\mathbb{R}^{2}} a_{1}(s) \, a_2(t) \, |t-s|^{\alpha} \, ds dt \neq 0$.

\

\textbf{Notation:} We use the following convention for the Fourier transform in $\mathbb{R}$ $\widehat{f}(\xi) = \int f(x) e^{-ix \xi} \, dx$. As usual we denote with $\mathcal{S}(\mathbb{R})$ the Schwartz space on $\mathbb{R}$.

\section{Preliminaries}

We start with the following lemma.

\begin{lemma} For $0 < \alpha < 1$ and $s \neq t \in \mathbb{R}$ fixed, let $K_{s,t}^{\alpha}$ be the function defined in $\mathbb{R}$ by
\[
K_{s,t}^{\alpha}(x)= \left(  |x-s|+|x-t| \right)^{\alpha - 1}, \,\,\,\, x \in \mathbb{R}.
\]
Then
\[
\widehat{K_{s,t}^{\alpha}}(\xi) = - 2^{\alpha} \, \Gamma(\alpha) \sin \left(\frac{(\alpha -1)\pi}{2} \right) e^{-i \xi(\frac{s+t}{2})}|\xi|^{-\alpha} + |t-s|^{\alpha - 1} \sgn{(t-s)} \int_{s}^{t} e^{-i x \xi} dx
\]
\[
- \, \frac{|t-s|^{\alpha}}{\alpha} e^{-i \frac{(s+t)}{2} \xi} \cos \left(\frac{|t-s| \xi}{2} \right) \, + \,
\frac{i \, 2^{\alpha} \, \xi e^{-i \xi(\frac{s+t}{2})}}{\alpha}
\int_{0}^{\frac{|t-s|}{2}} x^{\alpha} \sin(x \xi) dx,
\]
in the distributional sense.
\end{lemma}

\begin{proof} First we assume that $s < t$. Then for each $\phi \in \mathcal{S}(\mathbb{R})$ fixed, we have
\[
\left(\widehat{K_{s,t}^{\alpha}}, \phi \right)= \left(K_{s,t}^{\alpha}, \widehat{\phi} \, \right) = \int_{\mathbb{R}} K_{s,t}^{\alpha}(x)
\widehat{\phi}(x) dx
\]
\[
= \int_{t}^{+\infty}  K_{s,t}^{\alpha}(x) \widehat{\phi}(x) dx + \int_{s}^{t}  K_{s,t}^{\alpha}(x) \widehat{\phi}(x) dx +
\int_{-\infty}^{s}  K_{s,t}^{\alpha}(x) \widehat{\phi}(x) dx = I + II + III.
\]
Let us now proceed to compute each one of these integrals,
\[
I = \int_{\mathbb{R}} \chi_{(t, +\infty)}(x) \left( 2x -(s+t) \right)^{\alpha -1}
\, \widehat{\phi}(x) dx
\]
\[
= 2^{\alpha - 1} \int_{\mathbb{R}}  x^{\alpha -1} \chi_{(\frac{t-s}{2}, +\infty)}(x)
\, (e^{-i(\cdot) \frac{(s+t)}{2}}\phi) \,\, \widehat{} \,\, (x) dx
\]
\[
= 2^{\alpha - 1} \int_{\mathbb{R}} x_{+}^{\alpha - 1} (e^{-i(\cdot) \frac{(s+t)}{2}}\phi) \,\, \widehat{} \,\, (x) dx - 2^{\alpha - 1} \int_{\mathbb{R}} x^{\alpha - 1}
\chi_{(0, \frac{t-s}{2})}(x) (e^{-i(\cdot) \frac{(s+t)}{2}}\phi) \,\, \widehat{} \,\, (x) dx
\]
\[
=2^{\alpha - 1} \int_{\mathbb{R}} x_{+}^{\alpha - 1} (e^{-i(\cdot) \frac{(s+t)}{2}}\phi) \,\, \widehat{} \,\, (x) dx
\]
\[
- 2^{\alpha - 1} \int_{\mathbb{R}} \left( \frac{(t-s)^{\alpha} e^{-i \xi \frac{(t-s)}{2}}}{2^{\alpha} \alpha} + \frac{i \xi}{\alpha}
\int_{0}^{\frac{t-s}{2}} x^{\alpha} e^{-ix \xi} dx \right) e^{-i \xi \frac{(s+t)}{2}}\phi(\xi) d\xi,
\]
to compute $III$ we proceed as in $I$, thus
\[
III = 2^{\alpha - 1} \int_{\mathbb{R}} x_{-}^{\alpha - 1} (e^{-i(\cdot) \frac{(s+t)}{2}}\phi) \,\, \widehat{} \,\, (x) dx
\]
\[
- 2^{\alpha - 1} \int_{\mathbb{R}} \left( \frac{(t-s)^{\alpha}}{2^{\alpha}\alpha}e^{i \xi (\frac{t-s}{2})} - \frac{i \xi}{\alpha}
\int_{0}^{\frac{t-s}{2}} x^{\alpha} e^{ix \xi} dx \right) e^{-i \xi \frac{(s+t)}{2}}\phi(\xi) d\xi,
\]
so
\[
I + III= 2^{\alpha - 1} \int_{\mathbb{R}} |x|^{\alpha - 1} (e^{-i(\cdot) \frac{(s+t)}{2}}\phi) \,\, \widehat{} \,\, (x) dx
\]
\[
- \, \int_{\mathbb{R}} \frac{(t-s)^{\alpha}}{\alpha} e^{-i \frac{(s+t)}{2} \xi} \cos \left(\frac{(t-s) \xi}{2} \right) \phi(\xi) d\xi
\]
\[
+ \, \int_{\mathbb{R}} \left( \frac{i \, 2^{\alpha} \, \xi e^{-i \xi(\frac{s+t}{2})}}{\alpha}
\int_{0}^{\frac{(t-s)}{2}} x^{\alpha} \sin(x \xi) dx \right) \phi(\xi) d\xi.
\]
Now $II$ is easy, indeed
\[
II = \int_{\mathbb{R}} \chi_{(s,t)}(x) (t-s)^{\alpha - 1} \, \widehat{\phi}(x) dx
= \int_{\mathbb{R}} \left( (t-s)^{\alpha - 1} \int_{s}^{t} e^{-ix \xi} dx \right) \phi(\xi) d\xi.
\]
Since
\[
\widehat{|x|^{\alpha -1}}(\xi) = -2 \Gamma(\alpha) \sin \left(\frac{(\alpha -1)\pi}{2} \right) |\xi|^{-\alpha}
\]
(see equation (12), pp. 173, in \cite{Gelfand}), the lemma follows for the case $s < t$. Finally, exchanging the roles of $s$ and $t$ we obtain the statement of the lemma.
\end{proof}

\begin{corollary} If $a_1$ and $a_2$ are two bounded functions on $\mathbb{R}$ with compact support and such that $\int a_1 = 0$ or $\int a_2 = 0$, then
\[
\int_{\mathbb{R}} \mathcal{I}_{\alpha +1}(a_1, a_2)(x) dx = \frac{\alpha - 1}{\alpha} \iint_{\mathbb{R}^{2}} a_{1}(s) \, a_2(t) \, |t-s|^{\alpha} \, ds dt.
\]
\end{corollary}

\begin{proof} It is easy to check that $\mathcal{I}_{\alpha +1}(a_1, a_2)(\cdot) \in L^{1}(\mathbb{R})$. Let $\varphi \in \mathcal{S}(\mathbb{R})$ be an {\bf even} function such that $\varphi(0)=1$ and for $\epsilon > 0$ let $\varphi_{\epsilon}(x) = \varphi(\epsilon x)$. Since
\[
\int_{\mathbb{R}} \mathcal{I}_{\alpha +1}(a_1, a_2)(x) dx =
\lim_{\epsilon \to 0^{+}} \int_{\mathbb{R}} \mathcal{I}_{\alpha +1}(a_1, a_2)(x) \varphi_{\epsilon}(x)dx,
\]
we will proceed to compute this limit.
\[
\lim_{\epsilon \to 0^{+}} \int_{\mathbb{R}} \mathcal{I}_{\alpha +1}(a_1, a_2)(x) \varphi_{\epsilon}(x)dx =
\lim_{\epsilon \to 0^{+}} \iint_{\mathbb{R}^{2}} a_1(s) a_2(t) \left(\int_{\mathbb{R}} K_{s,t}^{\alpha}(x) \varphi_{\epsilon}(x)dx \right)dsdt
\]
\[
= \lim_{\epsilon \to 0^{+}} \iint_{\mathbb{R}^{2}} a_1(s) a_2(t) \left(\int_{\mathbb{R}} \widehat{K_{s,t}^{\alpha}}(\xi)
\widehat{\varphi_{\epsilon}}(\xi)d\xi \right)dsdt
\]
\[
= \iint_{\mathbb{R}^{2}} a_1(s) a_2(t) \lim_{\epsilon \to 0^{+}} \left(\int_{\mathbb{R}} \left( \widehat{K_{s,t}^{\alpha}}(\epsilon \, \xi) +
2^{\alpha} \, \Gamma(\alpha) \sin \left(\frac{(\alpha -1)\pi}{2} \right) |\epsilon \xi|^{-\alpha} \right) \widehat{\varphi}(\xi)d\xi \right)dsdt
\]
\[
= \frac{\alpha - 1}{\alpha} \iint_{\mathbb{R}^{2}} a_{1}(s) \, a_2(t) \, |t-s|^{\alpha} \, ds dt,
\]
where the third equality follows from the moment condition of $a_1$ (or $a_2$) and the last one from Lemma 1 and that $\varphi(0)=1$.
\end{proof}

\section{A Counter-example}

We take $a_1(s) = \chi_{(-1,0)}(s) - \chi_{(0,1)}(s)$ and $a_2(t) = a_1(t-2)$. From (S1) it follows that $a_1 \in H^{1}(\mathbb{R})$ and
$a_2 \in H^{(\alpha +1)^{-1}}(\mathbb{R})$ for each $0 < \alpha <1$. A computation gives

\[
\iint a_1(s) a_2(t) |t-s|^{\alpha} ds dt = \frac{4 \cdot 3^{\alpha + 2} - 4^{\alpha +2} - 6 \cdot 2^{\alpha +2} +4}{(\alpha +1) (\alpha +2)}
\neq 0.
\]
${}$\\
From (S2) and corollary 2 it obtains that $\mathcal{I}_{\alpha +1}(a_1, a_2)(\cdot) \notin H^{1}(\mathbb{R})$, for each $0 < \alpha <1$.
For $0 < p < (\alpha+1)^{-1}$ and $\frac{1}{q} = \frac{1}{p} - \alpha$, we take $N$ as any fixed integer with $N > p^{-1}-1$, then the set of all bounded, compactly supported functions for which $\int x^{\beta} f(x) dx =0$, for all $|\beta| \leq N$ is dense in $H^{r}(\mathbb{R})$ for each $p \leq r \leq 1$ (see 5.2 b), pp. 128, in \cite{St}). In particular, there exists $b \in H^{p}(\mathbb{R})$ such that $\|a_1 \|_{H^{1}} \| a_2 - b \|_{H^{(\alpha+1)^{-1}}} < \left| \int_{\mathbb{R}} \mathcal{I}_{\alpha+1}(a_1, a_2) (x) dx \right|/2C$. Then
\[
\left| \int_{\mathbb{R}} \mathcal{I}_{\alpha+1}(a_1, b) (x) dx \right| \geq \left| \int_{\mathbb{R}} \mathcal{I}_{\alpha+1}(a_1, a_2) (x) dx \right| -  \int_{\mathbb{R}} | \mathcal{I}_{\alpha+1}(a_1, a_2 - b) (x) | dx
\]
\[
\geq \left| \int_{\mathbb{R}} \mathcal{I}_{\alpha+1}(a_1, a_2) (x) dx \right| - C \|a_1 \|_{H^{1}} \| a_2 - b \|_{H^{(\alpha+1)^{-1}}}
\]
\[
> \frac{1}{2}\left| \int_{\mathbb{R}} \mathcal{I}_{\alpha+1}(a_1, a_2) (x) dx \right| > 0
\]
where the second inequality follows from Theorem 1.1 in \cite{Cruz-Uribe} with $p_1 = 1$, $p_2 = (\alpha+1)^{-1}$ and $q=1$. But then the operator $\mathcal{I}_{\alpha + 1}$ is not bounded from $H^{1}(\mathbb{R}) \times H^{p}(\mathbb{R})$ into $H^{q}(\mathbb{R})$ for $0 < p \leq (\alpha + 1)^{-1}$ and $\frac{1}{q}= \frac{1}{p} - \alpha$, since $\int_{\mathbb{R}} \mathcal{I}_{\alpha+1}(a_1, b) (x) dx \neq 0$.

\qquad

We conclude this note by summarizing our main result in the following theorem.

\begin{theorem} For $1 < \gamma < 2$, let $\mathcal{I}_{\gamma}$ be the multilinear fractional integral operator given by
\[
\mathcal{I}_{\gamma}(f_1,f_2)(x) = \iint_{\mathbb{R}^{2}} \frac{f_1(s) f_2(t)}{(|x-s|+|x-t|)^{2-\gamma}} ds dt, \,\,\,\, x \in \mathbb{R}.
\]
Then the operator $\mathcal{I}_{\gamma}$ is not bounded from $H^{1}(\mathbb{R}) \times H^{p}(\mathbb{R})$ into $H^{q}(\mathbb{R})$ for \,\,\,\,\,\,\,\,\,\, $0 < p \leq \gamma^{-1}$ and $\frac{1}{q}= 1 +\frac{1}{p} - \gamma$.
\end{theorem}


\section{Introduction}

Breast cancer is the second leading cause of cancer-related death among women in the United States. Screening mammography has been highly effective in reducing mortality rate but suffers from a high rate of false positive recalls, leading to added costs and stress. Better computer-aided diagnostic tools could improve patient outcomes by helping radiologists. Our goal is to build a system that provides highly accurate as well as interpretable predictions, in order to assist radiologists in cancer detection. To this end, we train an object detection network to predict the location of suspicious lesions and classify them.

\section{Data}

Our dataset \cite{NYU_dataset} comprises 229,426 screening mammography exams (1,001,093 images, with each exam containing the four standard views) from 141,473 unique patients. Among these, 5,832 exams had at least one biopsy performed within 120 days of a screening mammogram. Within this set of exams, 985 breasts had malignant findings, 5,556 had benign findings and 234 had both. For these exams, radiologists retrospectively annotated the locations of the biopsied lesions at a pixel level.

\section{Methods}

\label{sec:methods}

We used bounding boxes on the annotations of biopsied lesions to train a Feature-Pyramidal Network \cite{fpn}-based Faster-RCNN model \cite{faster_rcnn}. In most experiments, we used a ResNet-50 \citep[R-50;][]{resnet} backbone pretrained on ImageNet \cite{imagenet} to initialize our models. We also conducted experiments with ResNet-101 (R-101) and ResNeXt-101 \citep[X-101;][]{resnext} backbones. The final classifier predicts three classes: benign lesions, malignant lesions, or nothing. We used \texttt{maskrcnn-benchmark} \cite{mrcnn-benchmark} to train our models. We evaluated models after every 5,000 mini-batches and select the best models according to the score computed on the validation set. We used the standard setup of the framework, aside from the following exceptions.

\vspace{-1mm}
\paragraph{Use of non-annotated images} We used images without annotated lesions to train our model. These are treated as negatives for both the region-proposal network (RPN) and the classifier branch. The goal of this approach is to reduce overfitting on the comparatively small set of annotations. For our base experiment, during training, half of the images are sampled from an exam with a biopsy and half without.

\vspace{-1mm}
\paragraph{Resolution} We used a maximum resolution of $2200 \times 3000$ for the R-50 backbone, by resizing the mammograms with bilinear interpolation ($1700 \times 2700$ for R-101, $1300 \times 2100$ for X-101). Indeed, prior work \cite{multi-view} highlighted the importance of operating at a high resolution. Consequently, we also reduced the batch size to 4 (2 per GPU) due to our model's large memory requirements.

\vspace{-1mm}
\paragraph{IoU thresholds} We also relaxed the IoU threshold for foreground objects in the RPN from 0.7 to 0.5, as in \cite{detecting_lesions, deformable}. Compared to natural images, our data set has (i) noisier annotations with less precise boundaries, and (ii) fewer annotations per image. To circumvent (i), we also tried isotropic bounding box rescaling by a factor drawn uniformly at random in $\left[0.8, 1.2\right]$. We set the IoU threshold of the final non-maximum suppression to 0.1, following the rationale of \citet{detecting_lesions} that for mammograms, ``overlapping detections are expected to happen less often than in usual object detection".

\vspace{-1mm}
\paragraph{Optimization} We generally did not tune the learning rate schedules according to the advice of \citet{lr_schedule} (``recommended lr schedule'') despite lowering the batch size. Indeed, found that it did not strongly impact results while considerably increasing training time. We used gradient clipping (for norm $>3$) in all our experiments. We experimented with using batch normalization but found it not to help in our experiments.

\vspace{-1mm}
\paragraph{Inference}
During inference, we take the maximum malignant prediction among the bounding boxes of each view and then average over both views to get the predicted probability for the whole breast. We found this setup to perform better than others, including taking the maximum over views or the mean over boxes. In addition, we found it useful to decrease the score threshold at inference time from 0.05 to 0.001, as the tail-like behavior is important when computing AUC.

\section{Results}

\begin{table}[t]
\centering \small
\begin{tabular}{lrlrl}
\toprule
\textbf{model} & \multicolumn{2}{c}{\textbf{test set}} & \multicolumn{2}{c}{\textbf{reader study}}  \\
\midrule
\textbf{\citet{pnas_paper} single model} &  0.886 & $\pm$ 0.003 & - & \\
\textbf{\citet{pnas_paper} ensemble} & 0.895 & & 0.876 \\
\hline
\textbf{Base setup} & 0.891 & $\pm$ 0.005 & 0.845 & $\pm$ 0.007 \\
\textbf{+ biopsy ratio 0.75} & 0.895 & $\pm$ 0.004 & 0.855 & $\pm$ 0.012 \\
\textbf{+ biopsy ratio 1} & 0.887 & $\pm$ 0.011 & 0.855 & $\pm$ 0.013  \\
\textbf{+ bounding box scaling} & 0.890 & $\pm$ 0.006 & 0.841 & $\pm$ 0.010 \\
\textbf{+ classifier*5} & 0.903 & $\pm$ 0.007 & 0.859 & $\pm$ 0.013 \\
\textbf{+ bb scaling, classifier*5, ratio 0.75} & 0.888 & $\pm$ 0.006 & 0.855 & $\pm$ 0.006 \\
\textbf{+ recommended lr schedule} & 0.897 & $\pm$ 0.007 & 0.858 & $\pm$ 0.007 \\
\textbf{+ R-101 backbone} & 0.887 & $\pm$ 0.011 & 0.834 & $\pm$ 0.011 \\
\textbf{+ X-101 backbone} & \textbf{0.908} & $\pm$ 0.014 & \textbf{0.866} & $\pm$ 0.020 \\
\hline
\textbf{Ensemble} & \textbf{0.919} & & \textbf{0.879} \\
\textbf{Ensemble + \citet{pnas_paper} ensemble} & \textbf{0.930} & & \textbf{0.895} \\
\bottomrule
\end{tabular}
\caption{AUC comparison against earlier results on the test set and reader study. Means and standard deviations for our models are computed over 3 random initializations. Base refers to the R-50 setup. For biopsy ratio, we increased the proportion of images from exams with biopsies in training. For classifier*5, we quintuple the weight on the classifier loss.}
\label{tab:main_results}

\end{table}

The results are shown in \autoref{tab:main_results}. Our methods compare favorably against the models of \citet{pnas_paper}. Indeed, we achieve a 23\% relative error reduction in AUC on the test set. Note that both models use pixel-level annotations. Due to high variance, we ensemble all runs together for our "Ensemble" setup. On the subpopulation of the test set used for the reader study in \citet{pnas_paper}, our results also compare favorably to both earlier models and radiologists in terms of ROC, although the improvements are smaller. The reader study setup contains more exams from patients who underwent biopsy. Thus we believe this model distinguishes better between negative cases and cases that need a biopsy but is not as good in distinguishing between benign and malignant cases within the biopsied population. Ensembling our method with \citet{pnas_paper} improves the reader study AUC further to 0.895, compared to 0.778 for the average radiologist.

\section{Analysis}

\begin{wrapfigure}{r}{0.5\textwidth}
  \captionsetup{format=plain}
  \begin{center}
  \vspace{-8mm}
  \begin{tabular}{c c}
  \includegraphics[width=0.22\textwidth,trim=16mm 12mm 0mm 0mm]{354_mask.png} & \includegraphics[width=0.22\textwidth,trim=16mm 12mm 0mm 0mm]{354_preds.png}
  \end{tabular}
  \end{center}
  \vspace{-4mm}
  \caption{An example of an annotation of a malignant lesion (left) and predictions (right) from a X-101 model on a test set breast.}
  \label{fig:pred}
\end{wrapfigure}

Our validation metrics varied significantly between checkpoints. We believe that this is due to (i) the interaction of the components of the training loss, (ii) components of the training loss being only loosely related to the final metric, (iii) a small batch size rendering optimization unstable. This motivates ensembling all runs to produce our final model.

We notice a clear trade-off between using larger networks and using a larger resolution, with the R-50, R-101 and X-101 setups performing competitively. Unfortunately, due to computational constraints, we were unable to test R-101 and X-101 on the same resolution as R-50. Overall, our ensemble combining different backbones and resolutions performed best.

In Figure~\ref{fig:pred}, we show an annotation and the corresponding prediction using an X-101 model. The model accurately predicts a malignant lesion with high probability (0.99). It also predicts a benign lesion with low probability (0.36) for which there is no ground-truth annotation. These bounding boxes can highlight suspicious regions and help radiologists understand predictions from our models.

\midlacknowledgments{We would like to thank Yiqiu Shen and Jungkyu Park for comments on earlier versions of this manuscript. We also gratefully acknowledge the support of Nvidia Corporation with the donation of some of the GPUs used in this research. This work was supported in part by grants from the National Institutes of Health (R21CA225175 and P41EB017183).}