11institutetext: ICTEAM Institute, UCLouvain, Louvain-la-Neuve, Belgium
11email<EMAIL_ADDRESS>22institutetext:
Center of Artificial Intelligence Technology, Skoltech, Moscow, Russia
22email<EMAIL_ADDRESS>33institutetext: Dept. of Mathematics and
Operational Research, University of Mons, Mons, Belgium
33email<EMAIL_ADDRESS>
# Direct Exoplanet Detection Using
L1 Norm Low-Rank Approximation††thanks: This work was supported by the Fonds
de la Recherche Scientifique-FNRS under Grant no T.0001.23. Simon Vary is a
beneficiary of the FSR Incoming Post-doctoral Fellowship.

Hazan Daglayan 11 0009-0006-4843-6913 Simon Vary 11 0000-0002-9929-4357
Valentin Leplat 22 0000-0002-3313-1547 Nicolas Gillis 33 0000-0001-6423-6897
P.-A. Absil 11 0000-0003-2946-4178

###### Abstract

We propose to use low-rank matrix approximation using the component-wise
L1-norm for direct imaging of exoplanets. Exoplanet detection by direct
imaging is a challenging task for three main reasons: (1) the host star is
several orders of magnitude brighter than exoplanets, (2) the angular distance
between exoplanets and star is usually very small, and (3) the images are
affected by the noises called speckles that are very similar to the exoplanet
signal both in shape and intensity. We first empirically examine the
statistical noise assumptions of the L1 and L2 models, and then we evaluate
the performance of the proposed L1 low-rank approximation (L1-LRA) algorithm
based on visual comparisons and receiver operating characteristic (ROC)
curves. We compare the results of the L1-LRA with the widely used truncated
singular value decomposition (SVD) based on the L2 norm in two different
annuli, one close to the star and one far away.

###### Keywords:

L1 norm low-rank approximation Laplace distribution direct imaging exoplanet
detection

##  1 Introduction

In the field of exoplanet detection, the vast majority of planets (about 99%)
have been detected by indirect methods. Over the past decade, we have observed
the rapid development of high-contrast imaging as a promising technique for
the detection of exoplanets. Although very challenging, direct imaging
provides two main advantages compared to indirect methods. First, it gives
access to the photons of the planets, so we can obtain information about the
atmospheric composition of the planets [2]. Second, direct imaging allows us
to capture images of planets in a matter of hours, while indirect methods can
take years, as they often depend on observing the orbital motion of planet
around its host star. This temporal efficiency becomes especially critical
when searching for planets with wider orbits. Indirect methods may require an
extended duration to detect planets in these regions, while direct imaging
opens up opportunities for their discovery.

Due to the small angular separation between the host star and potential
exoplanets, direct imaging requires high-resolution images and, therefore,
large and high resolution ground-based telescopes. This causes the light to be
diffracted by atmospheric turbulence as it passes through the atmosphere.
Despite the use of coronagraphs, such as the well-known Lyot coronagraph or,
more recently, vortex coronagraphs [16], to reduce the high contrast caused by
the brightness of the star and adaptive optics techniques to avoid the
aberrations caused by the refraction of the light, standard signal processing
methods may still be unable to detect the planet in the images, due to the
residual aberrations in the form of quasi-static speckles that are often
brighter than the planet and resemble the planet in shape.

Angular differential imaging (ADI) is a widely used technique in astronomy to
reduce the effects of speckle noise in the images [14]. This technique is
based on observations made in pupil tracking mode, in which the star is fixed
in the center of the images as the Earth rotates in a night, which causes the
exoplanet to rotate around the host with time. ADI aims at building a
reference point spread function (PSF) that reproduces a model of the speckle
field to be subtracted from the images and aligned (in some methods, also
combined) with the signal of potential exoplanets.

Several algorithms are used to build the reference PSF in ADI. The most
popular ones are based on low-rank approximations to build the reference PSF:
principal component analysis (PCA) [19, 1], its annular version (AnnPCA) [12,
4], non-negative matrix factorization [18], the local low-rank plus sparse
plus Gaussian decomposition (LLSG) [11], and more recently low-rank plus
sparse trajectory [20]. Regime-switching model [7] also combines the
advantages of numerous PSF subtraction techniques using low-rank
approximations.

Methods based on the low-rank assumption are obtained by transforming the data
cube into a matrix such that each frame corresponds to a row of the matrix. In
the cases where we fit the low-rank approximation by minimizing the Frobenius
norm, i.e., by PCA (or equivalently, truncated SVD), this corresponds to the
maximum likelihood estimator under the i.i.d. white Gaussian noise assumption.
However, several recent lines of work observed that the residual datacube,
which is obtained by subtracting the low-rank part from the original data, is
more compatible with the Laplacian distribution that has heavier tails instead
of Gaussian [17, 6].

In this paper, we propose to perform the low-rank background approximation
using the component-wise L1 norm and to apply a data-dependent approximation.
This approach leverages two distinct annuli in the frames, one close and the
other far from the star. Firstly, it involves an analysis of the data to
identify whether the data is better fitted by a Gaussian or Laplacian
distribution. Subsequently, we compare the results both visually and through
ROC curve analysis.

The rest of this paper is structured as follows. In Section 2, we propose an
alternative to PCA in the context of exoplanet detection, namely a component-
wise L1 norm low-rank matrix approximation. We investigate different
statistical assumptions on the data and analyze the performances, then apply
the appropriate low-rank approximation to the data and present the
experimental results in Section 3. Finally, we conclude in Section 4 and
discuss potential future works.

##  2 Models and Methods

###  2.1 Planet flux estimation

Let $M\in\mathbb{R}^{T\times N^{2}}$ be a matrix of observations that consists
of $T$ unfolded frames with size $N\times N$, i.e., each row of the matrix
represents a single vectorized frame. The model for $M$ proposed in [6],
assuming a single planet located at position $g\in[N]\times[N]$ in the first
frame, is expressed as

$M=L+a{\color[rgb]{0,0,0}{}_{g}}P_{g}+E,\quad\mathrm{rank}(L)\leq{\color[rgb]{0,0,0}k},$
(1)

where $L$ is the low-rank background, $E$ is the noise,
$a{\color[rgb]{0,0,0}{}_{g}}$ is the intensity of the planet referred to as
the _flux_ , $P_{g}$ is the planet signature along the trajectory.

For such models based on low-rank approximations, the choice for the rank
value is crucial. Indeed, if the rank is too small, the signal of the speckles
will remain in the residual $M-L$, making it difficult to separate the signal
of the planet from the speckles. Conversely, if it is too large, the signal of
the planet will be captured by the low-rank matrix, and it will be challenging
to find the signal of the planets in the residual.

When the error $E$ is Gaussian distributed, the maximum likelihood estimator
for $L$ is obtained by the minimization of the Frobenius norm. Classical
methods, such as AnnPCA and LLSG, fit the low-rank component as follows:

$\displaystyle\hat{L}=\underset{L}{\arg\min}\|M-L\|_{F}\quad\text{{\color[rgb]{0,0,0}subject
to}}\quad\mathrm{rank}(L)\leq{\color[rgb]{0,0,0}k},$ (2)

where $\|A\|_{F}$ denotes the entry-wise L2-norm of $A$ (the Frobenius norm),
which can be solved using the truncated SVD. This step is followed by
subtracting the low-rank component and identifying the planet
$a{\color[rgb]{0,0,0}{}_{g}}P_{g}$ in the residual matrix by solving the
minimization problem

$\displaystyle\hat{a}_{g}=\arg\min_{a_{g}>0}\left\|M-\hat{L}-a_{g}P_{g}\right\|_{2},$
(3)

for all possible $P_{g}\in\Lambda$.

In order to detect the presence of a planet, (3) has to be solved for all
planet positions $g$ of potential interest, which are all the pixels except
for the ones that are too close to the star and those located at the corners
since they cannot adequately capture the rotational motion of the planet; see
Section 2.2 for the detection step.

Recently, it has been observed that the error term $E$ has heavy tails and
more closely follows the Laplacian distribution [17]. Consequently, it was
proposed in [6] to identify the planet using L1 minimization

$\displaystyle\hat{a}_{g}=\arg\min_{a_{g}>0}\left\|M-\hat{L}-a_{g}P_{g}\right\|_{1},$
(4)

where $\|\cdot\|_{1}$ is the component-wise L1 norm, that is,
$\|M\|_{1}=\sum_{i,j}|M[i,j]|$.

However, this approach introduces an inconsistency in the noise assumption
between the speckle subtraction (2) and planet identification (4). While L1
minimization proves effective for planet identification, it raises the
question of why not apply the L1 norm to background subtraction as well.

Consequently, we propose to fit the low-rank component using the component-
wise L1 norm

$\displaystyle\hat{L}=\underset{L}{\arg\min}\|M-L\|_{1}\quad\mathrm{s.t.}\quad\mathrm{rank}(L)\leq
k,$ (5)

followed by the planet estimation in L1 norm as stated in (4). This approach
allows us to maintain a consistent noise assumption across both the low-rank
speckle subtraction and the planet estimation phases. By doing so, we aim to
enhance the overall integrity and reliability of the analysis.

In fact, the proposed approach consists of addressing the optimization problem

$\min_{L,a_{g}}\|M-L-a_{g}P_{g}\|_{1}\text{ subject to }\mathrm{rank}(L)\leq
k$

by first minimizing over $L$, then minimizing over $a_{g}$. (An interesting
direction for future work would be to consider alternating minimization.
Preliminary results can be found in [5].)

The L1 low-rank approximation in (5) is an NP-hard problem, even in the rank-
one case [10]. Hence, most algorithms to tackle (5), such as alternating
convex optimization [13], the Wiberg algorithm [8], and augmented Lagrangian
approaches [21], do not guarantee to find a global optimal solution, unlike in
the case of PCA. Moreover, the computed solutions are sensitive to the
initialization of the algorithms.

We use Algorithm 1 (L1-LRA) suggested by [10] to solve (5). It solves the
problem using an exact block cyclic coordinate descent method, where the
blocks of variables are the columns of $\hat{U}$ and the rows of $\hat{V}$ of
the low-rank approximation $\hat{L}=\hat{U}\hat{V}$. We initialized the
algorithm with the truncated SVD solution, denoted by
${\color[rgb]{0,0,0}\mathrm{H}^{\mathrm{SVD}}_{k}({\cdot})}$ in Algorithm 1.
In our experiments, we apply an annular version, similar to annular PCA
(AnnPCA) [12, 4], that selects only the pixels of $M$ in a certain annulus.
Indeed, as the values of pixels decrease away from the star, it is usually
better to calculate the low-rank approximation of each annulus separately.

Input: Image sequence $M\in\mathbb{R}^{t\times n}$, rank
${\color[rgb]{0,0,0}k}$, maximum number of iteration maxiter

1

2$\dot{U},\dot{S},\dot{V}^{T}=$
${\color[rgb]{0,0,0}\mathrm{H}^{\mathrm{SVD}}_{k}(M)}$

3 $\hat{U}=\dot{U}\dot{S}$; $\quad\hat{V}=\dot{V}^{T}$

4

5for  _$i=$ 1: maxiter_ do

6 $R=M-\hat{U}\hat{V}$

7 for  _${\color[rgb]{0,0,0}j}=1:{\color[rgb]{0,0,0}k}$_ do

8 $R=R+\hat{U}[:,{\color[rgb]{0,0,0}j}]\hat{V}[{\color[rgb]{0,0,0}j},:]$

9
$\hat{U}[:,{\color[rgb]{0,0,0}j}]=\underset{u}{\min}\|R-u\hat{V}[{\color[rgb]{0,0,0}j},:]\|_{1}$

10
$\hat{V}[{\color[rgb]{0,0,0}j},:]=\underset{v}{\min}\|R^{T}-v\hat{U}[:,{\color[rgb]{0,0,0}j}]\|_{1}^{T}$

11 $R=R-\hat{U}[:,{\color[rgb]{0,0,0}j}]\hat{V}[{\color[rgb]{0,0,0}j},:]$

12

13 end for

14

15 end for

Algorithm 1 L1-LRA [10]

To solve the minimization problem in steps 7-8 of Algorithm 1, we use the
exact method from [9]; these subproblems are weighted median problems which
can be solved in closed form.

###  2.2 Planet detection

The detection procedure consists of declaring positive the positions $g$ where
a detection metric—e.g., an SNR or a likelihood ratio—exceeds a given
threshold. A likelihood ratio map $\Lambda$ was proposed in [6], which
consists of L1 norm likelihood ratios $\Lambda_{g}(R)$ based on maximizing
log-likelihood of Laplace distribution using the solution of (3) or (4)
because it has been shown to provide better results in practice [6]

$\displaystyle\log\Lambda_{g}(R)=-\\!\\!\\!\\!\\!\sum_{(t,r)\in{\Omega}_{g}}\\!\\!\\!\\!\frac{|R(t,r)-\hat{a}_{g}P_{g}(t,r)|-|{R(t,r)}|}{\sigma_{R(r)}},$
(6)

where $R=M-\hat{L}$, $\sigma_{R}$ is the standard deviation of $R$ computed
along the time dimension, and $\Omega_{g}$ is the set of indices $(t,r)$ of
pixels whose distance from the trajectory is smaller than half the diffraction
limit, for more details see [6]. We will also use the L2 norm $\Lambda$ map
which consists of L2 norm likelihood ratios $\Lambda_{g}(R)$ using the
solution of (3) or (4)

$\displaystyle\log\Lambda_{g}(R)=-\frac{1}{2}\sum_{\mathclap{(t,r)\in{\Omega}_{g}}}\frac{|R(t,r)-\hat{a}_{g}P_{g}(t,r)|^{2}-|{R(t,r)}|^{2}}{\sigma^{2}_{R(r)}}.$
(7)

##  3 Numerical Experiments

In order to analyze the performance of the L1-LRA algorithm for background
subtraction in exoplanet detection, we compare the algorithms visually by
plotting their log-likelihood detection maps $\Lambda$ and also using receiver
operating characteristic (ROC) curves. Additionally, we empirically
investigate the methods in terms of fitting the data to Gaussian and Laplacian
distributions. The Python codes of the implementations are publicly available
from https://github.com/hazandaglayan/l1lra_for_exoplanets.

We tested the algorithms using the publicly available dataset _sph3_ for the
exoplanet data challenge [3]. The ADI cube obtained with the VLT/SPHERE-IRDIS
instrument has 228 frames with size $160\times 160$,  and it has a total field
rotation of 80.5 degrees. It is a real dataset, but it has 5 synthetically
injected planets at different distances from the star.

###  3.1 Empirical estimation of the noise distributions

In order to analyze the suitability of different noise assumptions, we fit the
Gaussian and the Laplacian distribution to the residual data, i.e., the data
after subtracting the low-rank component using PCA or L1-LRA. We look at two
different annuli separately, one that is close to the star at $4\lambda/D$
separation and one more distant from the star at $10\lambda/D$ and measure the
goodness of fit visually and by the coefficient of determination $\rho^{2}$.

In Figure 1, we observe that the residual data follows a Gaussian distribution
after applying PCA and a Laplace distribution after applying L1-LRA in
general. However, the Laplacian distribution provides a better fit for the
residual cube distribution in the tails, regardless of whether PCA or L1-LRA
is used, for both small and large separations.

(a) PCA at 4$\lambda/D$

(b) L1-LRA at 4$\lambda/D$

(c) PCA at 10$\lambda/D$

(d) L1-LRA at 10$\lambda/D$

Figure 1: Residual cube after low-rank approximation applied for small and
large separations

To assess the quality of fit of a distribution to the residual data, we use
the correlation of determination as a metric. The correlation of determination
$\rho^{2}$ quantifies the relationship between two sample sets, denoted as
$X=(x_{1},\dots,x_{n})$ and $Y=(y_{1},\dots,y_{n})$, and is computed as
follows:

$\displaystyle\rho^{2}=\frac{\left(\sum_{i=1}^{n}(x_{i}-\bar{x})(y_{i}-\bar{y})\right)^{2}}{\left(\sum_{i=1}^{n}(x_{i}-\bar{x})^{2}\right)\left(\sum_{i=1}^{n}(y_{i}-\bar{y})^{2}\right)}.$
(8)

where $\bar{x}$ and $\bar{y}$ represent the means of the variables in samples
$X$ and $Y$. This metric aids in comparing the goodness of fit of different
distributions. A $\rho^{2}$ value close to 1 indicates that the data is more
consistent with the distribution. In our experiments, $X$ represents the
height of the bins of the data histogram, while $Y$ corresponds to the values
of the probability density function within these bins. In Table 1 and 2, we
compute $\rho^{2}$ value between the tails of the data distribution and the
Gaussian or Laplacian distributions. The highest $\rho^{2}$ is obtained with
L1-LRA for both small and large separations.

Table 1: The coefficient of determination $\rho^{2}$ for small separation 4$\lambda/D$ Rank | PCA | L1-LRA
---|---|---
Gaussian | Laplacian | Gaussian | Laplacian
5 | 0.9876 | 0.9918 | 0.9878 | 0.9920
10 | 0.9866 | 0.9890 | 0.9878 | 0.9938
15 | 0.9841 | 0.9948 | 0.9916 | 0.9972
Table 2: The coefficient of determination $\rho^{2}$ for large separation 10$\lambda/D$ Rank | PCA | L1-LRA
---|---|---
Gaussian | Laplacian | Gaussian | Laplacian
5 | 0.9866 | 0.9920 | 0.9894 | 0.9940
10 | 0.9859 | 0.9936 | 0.9872 | 0.9948
15 | 0.9912 | 0.9954 | 0.9922 | 0.9960

Figure 3: Left column: $\Lambda$ map results of the methods L1L1, L1L2, L2L1,
and L2L2 from top to bottom. The detected planets are circled in white and the
non-detected planets are circled in red. Right column: Intensity of the pixels
in the $\Lambda$ maps. Those maps have been flattened such that the horizontal
axis corresponds to the pixels of the $\Lambda$ map, and the vertical axis
indicates the intensity of the pixels in the $\Lambda$ map. The center pixels
of the planets and the pixels near the planet are shown in thick red and
orange dots, respectively. Other pixels are displayed as blue squares. The
black line corresponds to the highest blue square.

###  3.2 Performance comparison of L1-LRA and PCA

We test the performance advantage of using the L1-LRA instead of PCA
approximation in terms of visual quality and ROC curves tested on the sph3
test data [3]. We also provide ablation studies using combinations of mixed
models, where the low-rank $L$ is first subtracted from the original data $M$
using L1-LRA or PCA, then used in the planet detection. As such, we have four
possible algorithms to investigate; see Table 3.

Table 3: Descriptions of algorithms | Planet detection
---|---
Background subtraction | L1 (4) and (6) | L2 (3) and (7)
L1-LRA (5) |  L1L1 |  L1L2
PCA (2) |  L2L1 |  L2L2

We first compared the algorithms, visually in the likelihood ratio map
$\Lambda$ map using the dataset with five synthetic planets. We selected the
best performing rank according to the average likelihood ratio over the
locations of the injected planets. Based on this, we chose ranks 5, 6, 5, and
10 for L1L1, L1L2, L2L1, and L2L2, respectively.

Figure 3 shows the $\Lambda$ map and the intensity of pixels in the $\Lambda$
map obtained with the four algorithms. In each right-hand plot, the black line
shows the lowest threshold for which there is no false positive. At this
threshold, we observe that L1L1 and L2L1 detect the five injected planets,
whereas L1L2 and L2L2 miss one of them. Moreover, the likelihood ratios of the
planet pixels (the red dots in the figures) obtained by L1L1 are much larger
than the threshold.

To evaluate the performance of the four methods, we put them to the test using
synthetically generated data and examine their results using ROC curves in
Fig. 4. The detection map used as input for the ROC curve procedure of [6] is
the $\Lambda$ map. We deleted the five injected planets from the dataset using
VIP-HCI package [12, 4]. Then, we created  50 different datasets by injecting
two planets in each,  180 degrees apart, and placed at the separation
4$\lambda/D$. We set the intensity of each injected planet as  $1\sigma$ where
$\sigma$ is the standard deviation of the annulus. We applied the same
procedure to a larger separation of $10\lambda/D$ with the intensity
$0.6\sigma$. The four methods were tested for three different rank values
$k=\\{5,10,15\\}$.

In the ROC curve results, we focus on the number of true positives before the
first false positives are found, as done in [6]. In the ROC curves obtained by
injecting the planets at the small separation (left-hand column of Figure 4),
L1L1 always gives the best results. This is also the case for the large
separation with $k=5$ (top right plot of Figure 4). In the remaining two plots
(i.e., large separation and $k=\\{10,15\\}$), the four algorithms perform
rather similarly.  Moreover, we can see a relationship between the ROC curve
performances of the algorithms and the ${\color[rgb]{0,0,0}\rho^{2}}$ values.
L1L1, which corresponds to blue curves, has the highest
${\color[rgb]{0,0,0}\rho^{2}}$ value, while L2L2, which corresponds to red
curves, has the lowest value. In addition, L2L1 generally has the second-
highest ${\color[rgb]{0,0,0}\rho^{2}}$ value while generally performing the
second-best on ROC curves.

Figure 4:  ROC curves. Left column plots belong to data cubes where the
planets are injected in small separation (4$\lambda/D$). Right column plots
belong to data cubes where the planets are injected in large separation
(10$\lambda/D$).The rank is 5, 10, and 15 from top to bottom. Brown bars show
the coefficient of determination $\rho^{2}$ for each algorithm.

We also compared the L1-LRA and PCA algorithms using SNR map as a detection
map, which is a classical approach in direct imaging [15]. We use the
implementation of SNR available in the VIP package [12, 4]. Figure 5 shows the
ROC curves of SNR maps applied after PCA and L1-LRA. We use the same inputs
and procedure as described above. Except for the results of
${\color[rgb]{0,0,0}k}=5$, L1-LRA outperforms PCA in all ROC curves for both
separations.

Figure 5:  ROC curves of SNR maps. The left column plot belongs to data cubes
where the planets are injected in small separation (4$\lambda/D$), while the
right column plot belongs to data cubes where the planets are injected in
large separation (10$\lambda/D$).

##  4 Conclusion

In this study, we proposed a low-rank approximation algorithm in the sense of
the component-wise L1-norm, termed L1-LRA, for background subtraction to
detect exoplanets. A comparison with the commonly used PCA based on the L2
norm was conducted within two distinct annuli: one near the host star and the
other farther away. Through empirical analysis, we examined the statistical
noise assumptions of our models and evaluated the performance of the L1-LRA
algorithm using visual comparisons and ROC curves. In addition to the success
observed in these comparisons, our evaluation of performance metrics,
including the coefficient of determination $\rho^{2}$, indicates that the L1L1
approach consistently outperforms alternatives. The proposed L1-LRA algorithm
presents a promising solution for the challenging task of exoplanet detection.

There are several directions for future studies, such as using different
initializations for the L1 norm approximation or designing more efficient
algorithms (e.g., based on smoothing techniques). In concurrent work [5], we
have investigated a technique that consists in estimating both $L$ and $a_{g}$
simultaneously in (1).

## References

  * [1] Amara, A., Quanz, S.P.: Pynpoint : An image processing package for finding exoplanets: Pynpoint. Monthly Notices of the Royal Astronomical Society 427(2), 948–955 (Dec 2012). https://doi.org/10.1111/j.1365-2966.2012.21918.x
  * [2] Bowler, B.P., Nielsen, E.L.: Occurrence rates from direct imaging surveys. In: Handbook of Exoplanets, pp. 1967–1983. Springer International Publishing (2018). https://doi.org/10.1007/978-3-319-55333-7_155
  * [3] Cantalloube, F., et al.: Exoplanet imaging data challenge: benchmarking the various image processing methods for exoplanet detection. In: Adaptive Optics Systems VII. vol. 11448, pp. 1027–1062. SPIE (2020)
  * [4] Christiaens, V., Gonzalez, C., Farkas, R., Dahlqvist, C.H., Nasedkin, E., et al.: VIP: A Python package for high-contrast imaging. The Journal of Open Source Software 8(81), 4774 (Jan 2023). https://doi.org/10.21105/joss.04774
  * [5] Daglayan, H., Vary, S., Absil, P.A.: An alternating minimization algorithm with trajectory for direct exoplanet detection. In: ESANN 2023 proceedings, European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning (2023). https://doi.org/10.14428/esann/2023.ES2023-137
  * [6] Daglayan, H., Vary, S., Cantalloube, F., Absil, P.A., Absil, O.: Likelihood ratio map for direct exoplanet detection. In: 2022 IEEE 5th International Conference on Image Processing Applications and Systems (IPAS). vol. Five, pp. 1–5 (2022). https://doi.org/10.1109/IPAS55744.2022.10052997
  * [7] Dahlqvist, C.H., Cantalloube, F., Absil, O.: Regime-switching model detection map for direct exoplanet detection in adi sequences. Astronomy & Astrophysics 633, A95 (2020)
  * [8] Eriksson, A., Van Den Hengel, A.: Efficient computation of robust low-rank matrix approximations in the presence of missing data using the l 1 norm. In: 2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. pp. 771–778. IEEE (2010)
  * [9] Gillis, N., Plemmons, R.J.: Dimensionality reduction, classification, and spectral mixture analysis using non-negative underapproximation. Optical Engineering 50(2), 027001–027001 (2011)
  * [10] Gillis, N., Vavasis, S.A.: On the complexity of robust PCA and $\ell_{1}$-norm low-rank matrix approximation. Mathematics of Operations Research 43(4), 1072–1084 (2018)
  * [11] Gomez Gonzalez, C.A., et al.: Low-rank plus sparse decomposition for exoplanet detection in direct-imaging ADI sequences: The LLSG algorithm. Astronomy & Astrophysics 589, A54 (May 2016). https://doi.org/10.1051/0004-6361/201527387
  * [12] Gomez Gonzalez, C.A., et al.: VIP : Vortex Image Processing Package for High-contrast Direct Imaging. The Astronomical Journal 154(1), 7 (Jun 2017). https://doi.org/10.3847/1538-3881/aa73d7
  * [13] Ke, Q., Kanade, T.: Robust $l_{1}$ norm factorization in the presence of outliers and missing data by alternative convex programming. In: 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05). vol. 1, pp. 739–746 vol. 1 (2005). https://doi.org/10.1109/CVPR.2005.309
  * [14] Marois, C., Lafreniere, D., Doyon, R., Macintosh, B., Nadeau, D.: Angular Differential Imaging: A Powerful High-Contrast Imaging Technique. The Astrophysical Journal 641(1), 556–564 (Apr 2006). https://doi.org/10.1086/500401
  * [15] Mawet, D., Milli, J., Wahhaj, Z., Pelat, D., Absil, O., Delacroix, C., Boccaletti, A., Kasper, M., Kenworthy, M., Marois, C., Mennesson, B., Pueyo, L.: Fundamental limitations of high contrast imaging set by small sample statistics. The Astrophysical Journal 792(2), 97 (aug 2014). https://doi.org/10.1088/0004-637x/792/2/97
  * [16] Mawet, D., Riaud, P., Absil, O., Surdej, J.: Annular groove phase mask coronagraph. The Astrophysical Journal 632(2), 1191 (2005)
  * [17] Pairet, B., Cantalloube, F., Gomez Gonzalez, C.A., Absil, O., Jacques, L.: STIM map: Detection map for exoplanets imaging beyond asymptotic Gaussian residual speckle noise. Monthly Notices of the Royal Astronomical Society 487(2), 2262–2277 (Aug 2019). https://doi.org/10.1093/mnras/stz1350
  * [18] Ren, B., Pueyo, L., Zhu, G.B., Debes, J., Duchêne, G.: Non-negative matrix factorization: robust extraction of extended structures. The Astrophysical Journal 852(2), 104 (2018)
  * [19] Soummer, R., Pueyo, L., Larkin, J.: Detection and characterization of exoplanets and disks using projections on Karhunen-Loève eigenimages. The Astrophysical Journal 755(2), L28 (Aug 2012). https://doi.org/10.1088/2041-8205/755/2/L28
  * [20] Vary, S., Daglayan, H., Jacques, L., Absil, P.A.: Low-rank plus sparse trajectory decomposition for direct exoplanet imaging. In: ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). pp. 1–5. IEEE (2023)
  * [21] Zheng, Y., Liu, G., Sugimoto, S., Yan, S., Okutomi, M.: Practical low-rank matrix approximation under robust l1-norm. In: 2012 IEEE Conference on Computer Vision and Pattern Recognition. pp. 1410–1417. IEEE (2012)


∎

e1e-mail<EMAIL_ADDRESS>
11institutetext: Institut für Kernphysik, Technische Universität Darmstadt,
D-64289 Darmstadt, Germany 22institutetext: Departamento de Física Atómica,
Molecular y Nuclear, Facultad de Física, Universidad de Sevilla, Apartado
1065, E-41080 Sevilla, Spain 33institutetext: Instituto Interuniversitario
Carlos I de Física Teórica y Computacional (iC1), Apdo. 1065, E-41080 Sevilla,
Spain 44institutetext: Centro Nacional de Aceleradores, U. Sevilla, J.
Andalucía, CSIC, Avda Thomas A Edison, 7, E-41092 Sevilla, Spain
55institutetext: School of Physics Science and Engineering, Tongji University,
Shanghai 200092, China 66institutetext: Institute for Advanced Study of Tongji
University, Shanghai 200092, China

# The Hussein-McVoy formula for inclusive breakup revisited

A Tribute to Mahir Hussein

M. Gómez Ramos addr1  J. Gómez-Camacho addr2,addr4  Jin Lei addr5,addr6  A. M.
Moroe1,addr2,addr3

(Received: date / Accepted: date)

###### Abstract

In 1985, Hussein and McVoy [Nuc. Phys. A445 (1985) 124] elucidated a formula
for the evaluation of the nonelastic breakup (“stripping”) contribution in
inclusive breakup reactions. The formula, based on the spectator core model,
acquires a particularly simple and appealing form in the eikonal limit, to the
extent that it has become the standard procedure to analyze single-nucleon
knockout reactions at intermediate energies. In this contribution, a critical
assessment of this formula is presented and its connection with other,
noneikonal expressions discussed. Some calculations comparing the different
formulae are also presented for the one-nucleon removal of 14O+9Be reaction at
several incident energies.

###### Keywords:

Breakup reactions Inclusive breakup Hussein-McVoy model

###### pacs:

PACS code1 PACS code2 more

##  1 Introduction

Breakup reactions have been extensively used to extract nuclear structure
information (binding energies, spectroscopic factors, electric response to the
continuum, etc) and have also permitted to improve our understanding of the
dynamics of reactions among composite systems. When the projectile dissociates
into two fragments, the process can be described as an effective three-body
problem, which can be schematically represented as $a+A\rightarrow b+x+A$,
where $a$ represents the projectile which eventually dissociates into $b+x$.
Even in this simplified three-body picture, the theoretical description of the
process is not straightforward due to the presence of three particles in the
final state.

In some applications, one is interested in the inclusive process in which only
one of the fragments (say, $b$) is measured experimentally, that we represent
schematically as $A(a,b)X$. These inclusive cross sections are needed, for
example, in the application of the surrogate method Esc12  and in
spectroscopic studies by means of intermediate-energy knockout reactions Han03
; Tos01 ; Tos14 .

The evaluation of inclusive breakup reactions poses a challenging theoretical
problem because many processes can in principle contribute to the $b$ singles
cross section. When the two fragments $b$ and $x$ “survive” and the target
remains in its ground state, the process is referred to as elastic breakup
(denoted EBU hereafter), also called diffraction dissociation.

The remaining part of the inclusive breakup cross section, that we denote
globally as nonelastic breakup (NEB), includes those processes in which the
$x$ particle interacts nonelastically with the target nucleus. This involves,
for example, the transfer of $x$ to a bound state of the residual system
$B=A+x$, the fusion of $x$ forming a compound nucleus (incomplete fusion) or
simply the target excitation by $x$. If $x$ is a composite system, it also
includes any process in which the latter is broken or excited in any way. The
explicit evaluation of all these processes is not possible in general so
several authors proposed closed-form formulae which avoid the sum over the
final states. Interestingly, all these formulae display a common structure,
given by

$\left.\frac{d^{2}\sigma}{dE_{b}d\Omega_{b}}\right|_{\mathrm{NEB}}=-\frac{2}{\hbar
v_{a}}\rho_{b}(E_{b})\langle\varphi_{x}|W_{x}|\varphi_{x}\rangle,$ (1)

where $\rho_{b}(E_{b})=k_{b}\mu_{b}/[(2\pi)^{3}\hbar^{2}]$ is the density of
states (with $\mu_{b}$ the reduced mass of $b+B$ and $k_{b}$ their relative
wave number), $W_{x}$ is the imaginary part of the optical potential $U_{x}$,
which describes $x+A$ elastic scattering. Expression (1) offers a intuitively
appealing interpretation of nonelastic breakup. As particle $b$ is scattered,
fragments $x$ and $A$ interact with each other. In the original Hamiltonian,
the interaction between $x$ and $A$ will be represented by a real operator
depending on the internal degrees of freedom of the target nucleus. After
application of a Feshbach reduction, this interaction is replaced by the
complex potential $U_{xA}$, describing $x+A$ elastic scattering and whose
imaginary part accounts for nonelastic breakup events. Indeed, the difficulty
in this interpretation is to provide a proper description of the state
$|\varphi_{x}\rangle$ (the $x$-channel wave function hereafter), which should
describe the relative motion of $x$ and $A$, compatible with the incoming
boundary conditions, and with the fact that fragment $b$ will be finally
detected with a given momentum $\vec{k}_{b}$.

One of the first of these expressions was due to Hussein and McVoy HM85 ,
given explicitly in the next section. In the same work, they also derived an
approximate expression obtained by treating the distorted waves appearing in
(1) in the Glauber (also referred to as eikonal) approximation. Since this
approximation is valid at high energies, the HM formula so obtained is
expected to be accurate also at high energies. In fact, this eikonal formula
is the key tool to evaluate the NEB part of the inclusive cross section in
nucleon removal knockout reactions used to study spectroscopy of nucleon hole-
states. In these experiments, measured cross sections are compared with
theoretical calculations for EBU and NEB, with the latter being evaluated with
the Glauber version of the HM formula. These theoretical cross sections are
commonly evaluated assuming single-particle wavefunctions for the removed
nucleon and they are later multiplied by the required spectroscopic factors
derived, for example, from shell-model calculations. Then, the ratio
$R_{s}=\sigma^{\mathrm{exp}}/\sigma^{\mathrm{theo}}$ is computed. Typically,
one obtains $R_{s}<1$, which has been interpreted as an effect of additional
correlations not present in small-scale shell-model calculations, presumably
leading to a larger fragmentation of single-particle strengths (and a
subsequent reduction of spectroscopic factors). Moreover, these studies have
found a systematic dependence of this ratio on the separation energy of the
removed nucleon, with $R_{s}$ becoming smaller and smaller as the separation
energy becomes larger Tos14 . Some authors have interpreted this result as an
indication of additional correlations (coming from tensor and short-range
components of the nucleon-nucleon interaction). However, this interpretation
has been recently put into question by other authors, because this trend is
apparently not observed in other reactions, such as transfer Fla13  and
$(p,pN)$ reactions Ata18 ; Kaw18 ; Gom18 .

Clearly, the conclusions strongly rely on the validity of the formulae used in
the evaluation of the inclusive cross sections. Whereas there is a general
consensus on the evaluation of the EBU cross sections, with different
approaches leading to consistent results (as, for instance, the distorted-wave
Born approximation (DWBA) Bau83 , the continuum-discretized coupled-channels
(CDCC) method Aus87  and a variety of semiclassical approaches Typ94 ; Esb96 ;
Kid94 ; Cap04 ), the reliability of the NEB calculations has been more
controversial. One of the criticisms concern the validity of the Glauber HM
formula at the energies commonly used in these experiments (several tens of
MeV per nucleon).

In this paper, we revisit the HM formula and its Glauber limit, and discuss
its connection with other inclusive breakup formulae proposed by other
authors. We present also preliminary numerical calculations comparing some of
these formulae, which show that the Glauber limit of the HM formula is a
reasonable approximation at the energies at which the nucleon-knockout
reactions have been measured.

Figure 1: Relevant coordinates for the three-body models described in the
text.

##  2 Review of inclusive breakup formulae

###  2.1 The Hussein-McVoy (HM) formula

One of the first attempts to provide a closed-form formula for the inclusive
breakup cross section was due to Hussein and McVoy in their seminal 1985 paper
HM85 . The HM derivation makes use of the spectator assumption for the
detected fragment $b$, which means that this fragment does not participate
directly in the breakup, so that the breakup is produced by any nonelastic
scattering of the participant $x$ with the target. By summing over all $x+A$
states which leave $A$ in an excited state, they arrived to the following
formula for the double differential cross section for NEB:

$\left.\frac{d^{2}\sigma}{dE_{b}d\Omega_{b}}\right|^{\mathrm{HM}}_{\mathrm{NEB}}=-\frac{2}{\hbar
v_{a}}\rho_{b}(E_{b})\langle\varphi^{\mathrm{HM}}_{x}|W_{x}|\varphi^{\mathrm{HM}}_{x}\rangle,$
(2)

where $\varphi^{\mathrm{HM}}_{x}$ is defined by the so-called non-
orthogonality overlap, which is given by:

$\varphi^{\mathrm{HM}}_{x}(\vec{r}_{x})=\langle{\vec{r}}_{x}|\varphi^{\mathrm{HM}}_{x}\rangle=\langle{\vec{r}}_{x}\chi^{(-)}_{b}|\chi_{a}^{(+)}\phi_{a}\rangle.$
(3)

Here, $|\phi_{a}\rangle$ is the projectile ground state,
$|\chi^{(+)}_{a}\rangle$ and $|\chi^{(-)}_{b}\rangle$ are distorted scattering
states for the $a+A$ and $b+B$ systems, respectively. $|\vec{r}_{x}\rangle$ is
a state with given separation of fragment and target. The relevant coordinates
are shown in Fig. 1.

###  2.2 The Eikonal Hussein-McVoy formula (EHM)

Hussein and McVoy obtained further insight on their formula (2) by treating
the distorted waves $\chi^{(+)}_{a}(\vec{r}_{a})$ and
$\chi^{(-)}_{b}(\vec{r}_{b})$ in the Glauber (also known as eikonal)
approximation. The Glauber approximation to an elastic-scattering distorted
wave is:

$\chi_{\vec{k}}^{(+)}({\vec{r}})=\mathrm{e}^{i\vec{k}\cdot\vec{r}}\exp\left[+i\int_{-\infty}^{z}\Delta
k\left(z^{\prime},b\right)\mathrm{d}z^{\prime}\right]$ (4)

where the incident momentum $\vec{k}$ points along the positive $z$-axis, and
$\vec{b}$ is the component of $\vec{r}$ perpendicular to $z$. The exponent in
the second factor is

$\Delta k\left(z^{\prime},b\right)\equiv-\frac{k}{2E}U(z,b),$ (5)

which, when integrated along the entire trajectory, gives the optical phase
shift

$2\delta(b)=\int_{-\infty}^{\infty}\Delta
k\left(z^{\prime},b\right)\mathrm{d}z^{\prime}=2\int_{0}^{\infty}\Delta
k\left(z^{\prime},b\right)\mathrm{d}z^{\prime},$ (6)

that is related to the partial-wave optical S-matrix, i.e.,
$S(b)=\exp[{2i\delta(b)}]$.

A key step in the HM derivation is the choice of the $U_{a}$ potential,
distorting the incident wave, which is taken as the sum of the corresponding
fragment-target potentials:

$U_{a}=U_{bA}+U_{xA}.$ (7)

As we shall see, this is a particularly fortunate choice, which takes into
account breakup effects in the entrance channel. Usual DWBA approaches would
use a distorting potential depending only on the $a-A$ coordinate, which does
not induce breakup. Indeed, they could get away with this complicated
distorting potential because the eikonal approximations were assumed: the $x$
and $b$ fragments move with the same average velocity as the projectile and
hence their momenta are given by

$\vec{k}_{x}=\left(m_{x}/m_{\mathrm{a}}\right)\vec{k}_{\mathrm{a}},\quad\vec{k}_{b}=\left(m_{b}/m_{a}\right)\vec{k}_{a}.$
(8)

With the particular choice (7) and the assumption (8) one obtains the
following result

$\displaystyle\varphi^{\mathrm{EHM}}_{x}({\vec{r}}_{x})=$
$\displaystyle\int\mathrm{d}^{3}{\vec{r}}_{b}\chi_{b}^{(-)*}({\vec{r}}_{b})\chi_{a}^{(+)}({\vec{r}}_{a})\phi_{a}({\vec{r}}_{bx})$
(9) $\displaystyle=$
$\displaystyle\mathrm{e}^{i\vec{k}_{x}\cdot{\vec{r}}_{x}}\exp\left[i\int_{-\infty}^{z_{x}}\Delta
k_{x}\left(z^{\prime},b_{x}\right)\mathrm{d}z^{\prime}\right]$
$\displaystyle\times\int\mathrm{d}^{3}{\vec{r}}_{b}\mathrm{e}^{i\vec{q}\cdot{\vec{r}}_{b}}S_{bA}\left(b_{b}\right)\phi_{\mathrm{a}}\left({\vec{r}}_{bx}\right)$

with $\vec{q}=\vec{k}_{b}-\vec{k}^{\prime}_{b}$, the average momentum
transferred in $b-A$ elastic scattering.

The last factor is then conveniently expressed as

$\int\mathrm{d}^{3}{\vec{r}}_{b}\mathrm{e}^{i\vec{q}\cdot{\vec{r}}_{b}}S_{bA}\left(b_{b}\right)\phi_{\mathrm{a}}\left({\vec{r}}_{bx}\right)\equiv\mathrm{e}^{iq_{\parallel}z_{x}}\tilde{\phi}_{a,b}(\vec{q},\vec{b}_{x}).$
(10)

Replacing this result into (2) one obtains for the double differential cross
section:

$\displaystyle\left.\frac{d^{2}\sigma}{dE_{b}d\Omega_{b}}\right|^{\mathrm{EHM}}_{\mathrm{NEB}}$
$\displaystyle=\frac{2}{\hbar v_{a}}\rho_{b}(E_{b})\frac{E_{x}}{k_{x}}$
$\displaystyle\times\int\mathrm{d}^{2}\vec{b}_{x}\left|\tilde{\phi}_{a,b}(\vec{q},\vec{b}_{x})\right|^{2}\left[1-\left|S_{xA}(b_{x})\right|^{2}\right].$
(11)

It should be noticed that the NEB depends only on the asymptotic properties,
this is, the $S$ matrices, of the interaction of $b$ and $x$ with the target.
There is no sensitivity on the wavefunctions in the interaction region. This
is a result of the eikonal approximation, plus the particular choice of the
distorted interaction, which included the imaginary potential $W_{xA}$ which
ultimately generates the NEB.

In many applications, one is interested in the total yield of fragment $b$,
which is obtained upon integration of the previous formula over the angular
and energy variables, resulting:

$\displaystyle\sigma_{\mathrm{NEB}}^{\mathrm{EHM}}=$
$\displaystyle\frac{2}{v_{\mathrm{a}}}(2\pi)^{3}\frac{E_{x}}{\hbar
k_{x}}\int\mathrm{d}^{3}\vec{r}_{b}\mathrm{d}^{3}\vec{r}_{x}\left|\phi_{a}\left(\vec{r}_{bx}\right)\right|^{2}$
$\displaystyle\times$
$\displaystyle\left|S_{bA}(b_{b})\right|^{2}\left[1-\left|S_{xA}\left(b_{x}\right)\right|^{2}\right].$
(12)

This equation has an appealing and intuitive form: the integrand contains the
product of the probabilities for the core being elastically scattered by the
target, $|S_{bA}(b_{b})|^{2}$, times the probability of the valence particle
being absorbed, $(1-|S_{xA}(b_{x})|)^{2}$. These probabilities are weighted by
the projectile wave function squared, and integrated over all possible impact
parameters. Because of the Glauber approximation, Eq. (2.2) is expected to be
accurate at high energies (above $\sim$100 MeV per nucleon). In fact, this
formula has been extensively employed in the analysis of intermediate-energy
knockout reactions (see e.g. Han03 ; Tos01 ; Tos14  and references therein)
mostly aimed at obtaining spectroscopic information of nucleon hole states.

###  2.3 The three-body (3B) model of Austern et al.

In Aus87 , Austern and collaborators derived a three-body formula for the
inclusive breakup cross section using as starting point the post-form
representation of the exact transition amplitude

$\displaystyle\frac{d^{2}\sigma}{d\Omega_{b}E_{b}}$
$\displaystyle=\frac{2\pi}{\hbar
v_{a}}\rho(E_{b})\sum_{c}|\langle\chi^{(-)}_{b}\Psi^{c,(-)}_{xA}|V_{\mathrm{post}}|\Psi^{(+)}\rangle|^{2}$
$\displaystyle\times\delta(E-E_{b}-E^{c}),$ (13)

where $V_{\mathrm{post}}\equiv V_{bx}+U_{bA}-U_{bB}$ is the post-form
transition operator, $\Psi^{(+)}$ is the system wavefunction with the incident
wave in the $a+A$ channel, and $\Psi^{c,(-)}_{xA}$ are the eigenstates of the
$x+A$ system, with $c=0$ denoting the $x$ and $A$ ground states. Thus, for
$c=0$ this expression gives the EBU part, whereas the terms $c\neq 0$ give the
NEB contribution.

In this model, excitations of the target are not explicitly considered
(although they are effectively taken into account by means of the optical
potentials $U_{xA}$ and $U_{bA}$). Thus, the total wavefunction is given by

$|\Psi^{(+)}\rangle\approx|\Psi^{\mathrm{3B(+)}}\phi^{0}_{A}\rangle,$ (14)

where $|\phi^{0}_{A}\rangle$ is the target ground-state and
$|\Psi^{\mathrm{3B(+)}}\rangle$ is the solution of the three-body equation:

$[\hat{T}_{aA}+\hat{T}_{bx}+V_{bx}+U_{bA}+U_{xA}-E]|\Psi^{\mathrm{3B(+)}}\rangle=0.$
(15)

Using the Feshbach projection formalism and the optical model reduction they
obtain a closed-form expression for the inclusive breakup cross section. The
latter can be further split into its EBU and NEB components Kas82 . The EBU is
given by

$\displaystyle\left.\frac{d^{2}\sigma}{dE_{b}d\Omega_{b}}\right|^{\mathrm{3B}}_{\mathrm{EBU}}$
$\displaystyle=\frac{2\pi}{\hbar
v_{a}}\rho(E_{b})\int{d\Omega_{x}}\,\rho(E_{x})$
$\displaystyle\times\left|\left\langle\chi_{b}^{(-)}\chi_{x}^{(-)}\left(\mathbf{k}_{x}\right)\left|V_{\mathrm{post}}\right|\Psi^{\mathrm{3B(+)}}\right\rangle\right|^{2}$
(16)

where $\rho(E_{x})$ is the density of states of particle $x$, whereas the NEB
part is given by.

$\left.\frac{d^{2}\sigma}{dE_{b}d\Omega_{b}}\right|^{\mathrm{3B}}_{\mathrm{NEB}}=-\frac{2}{\hbar
v_{a}}\rho_{b}(E_{b})\langle\varphi^{\mathrm{3B}}_{x}|W_{x}|\varphi^{\mathrm{3B}}_{x}\rangle.$
(17)

The three-body $x$-channel wavefunction
$\varphi^{\mathrm{3B}}_{x}({\vec{r}}_{x})$ is obtained by solving the
inhomogeneous equation

$(E^{+}_{x}-K_{x}-{U}_{x})\varphi^{\mathrm{3B}}_{x}({\vec{r}}_{x})=\langle{\vec{r}}_{x}\chi_{b}^{(-)}|V_{\mathrm{post}}|\Psi^{\mathrm{3B(+)}}\rangle$
(18)

where $E_{x}=E-E_{b}$. This equation can also be written in integral form as

$\displaystyle\varphi^{\mathrm{3B}}_{x}({\vec{r}}_{x})$
$\displaystyle=G^{(+)}_{x}\langle{\vec{r}}_{x}\,\chi_{b}^{(-)}|V_{\mathrm{post}}|\Psi^{\mathrm{3B(+)}}\rangle$
$\displaystyle=\int d^{3}\vec{r}\mkern
2.0mu\vphantom{r}^{\prime}_{x}\langle{\vec{r}}_{x}|G^{(+)}_{x}|{\vec{r}\mkern
2.0mu\vphantom{r}^{\prime}_{x}}\rangle\langle\vec{r}\mkern
2.0mu\vphantom{r}^{\prime}_{x}\,\chi_{b}^{(-)}|V_{\mathrm{post}}|\Psi^{\mathrm{3B(+)}}\rangle$
(19)

where $G^{(+)}_{x}=(E^{+}_{x}-K_{x}-{U}_{x}-i\epsilon)^{-1}$ is the optical
model Green’s function of particle $x$.

Austern et al. derived also an interesting alternative expression for
$\varphi^{\mathrm{3B}}_{x}({\vec{r}}_{x})$, namely,

$\varphi^{\mathrm{3B}}_{x}({\vec{r}}_{x})=\langle{\vec{r}}_{x}\chi_{b}^{(-)}|\Psi^{\mathrm{3B(+)}}\rangle.$
(20)

It is worth noting that the formulae (20) and (2.3) are formally equivalent
and, as such, they must provide identical results. Although (20) may seem
simpler, in practice, it requires that the three-body wavefunction
$\Psi^{\mathrm{3B(+)}}$ be accurate in the full configuration space, since
there is no natural cutoff in the integration variable ${\vec{r}}_{b}$. By
contrast, in Eqs. (18) and (2.3), the presence of the $V_{\mathrm{post}}$
operator will tend to emphasize small $b-x$ separations and hence one requires
only an approximate three-body wavefunction accurate within that range. This
can be achieved, for instance, expanding $\Psi^{\mathrm{3B(+)}}$ in terms of
$b-x$ eigenstates, as done in the CDCC method Aus87 , or in terms of Weinberg
states Pan13 . The implementation of the method with CDCC wavefunctions is
numerically challenging and the first calculation of this kind was only
recently reported Jin19 .

###  2.4 The Ichimura, Austern, Vincent (IAV) formula

Ichimura, Austern and Vincent IAV85  proposed a simpler DWBA version of the 3B
formula above. In DWBA, the exact wavefunction $\Psi^{(+)}$ is approximated by
the factorized form:

$|\Psi^{(+)}\rangle\approx|\chi_{a}^{(+)}\phi_{a}\phi^{0}_{A}\rangle.$ (21)

With this approximation, the NEB component of the $b$ singles cross section
becomes

$\left.\frac{d^{2}\sigma}{dE_{b}d\Omega_{b}}\right|^{\mathrm{IAV,post}}_{\mathrm{NEB}}=-\frac{2}{\hbar
v_{a}}\rho_{b}(E_{b})\langle\varphi^{\mathrm{post}}_{x}|W_{x}|\varphi^{\mathrm{post}}_{x}\rangle,$
(22)

that is formally identical to (17) but with the $x$-channel wavefunction given
now by

$\varphi^{\mathrm{post}}_{x}({\vec{r}}_{x})=G^{(+)}_{x}\langle{\vec{r}}_{x}\,\chi_{b}^{(-)}|V_{\mathrm{post}}|\chi_{a}^{(+)}\phi_{a}\rangle.$
(23)

The IAV model has been recently revisited by several groups Car16 ; Jin15 ;
Pot15  and its accuracy assessed against experimental data with considerable
success Jin17 ; Pot17 .

###  2.5 The Udagawa, Tamura (UT) formula

Udagawa and Tamura Uda81  derived a formula similar to that of IAV, but making
use of the prior form DWBA. Their final result can be again expressed in the
form (1),

$\left.\frac{d^{2}\sigma}{dE_{b}d\Omega_{b}}\right|^{\mathrm{UT}}_{\mathrm{NEB}}=-\frac{2}{\hbar
v_{i}}\rho_{b}(E_{b})\langle\varphi^{\mathrm{prior}}_{x}|W_{x}|\varphi^{\mathrm{prior}}_{x}\rangle,$
(24)

where $\varphi^{\mathrm{prior}}_{x}$ is the solution of a $x-A$ inhomogeneous
equation similar to Eq. (23) but replacing in the source term the post-form
transition operator, $V_{\mathrm{post}}$, by its prior form counterpart,
$V_{\mathrm{prior}}\equiv U_{xA}+U_{bA}-U_{a}$, i.e.:

$(E^{+}_{x}-K_{x}-{U}_{x})\varphi^{\mathrm{prior}}_{x}({\vec{r}}_{x})=\langle{\vec{r}}_{x}\,\chi_{b}^{(-)}|V_{\mathrm{prior}}|\chi^{(+)}_{a}\phi_{a}\rangle.$
(25)

##  3 Relation among theories

In this section we discuss the connection between the formalisms outlined in
the previous section, with emphasis in the HM formulae. Some of these
relations have already been discussed in previous works, most notably in the
comprehensive work of Ichimura Ich90 . We focus here on the HM formulae (2)
and (2.2) due their relevance in the analysis and interpretation of knockout
studies.

The HM formula (2) can be readily obtained starting from Austern’s identity
(20), and making the replacement
$|\Psi^{\mathrm{3B(+)}}\rangle\approx|\chi_{a}^{(+)}\phi_{a}\rangle$. It is
also enlightening to see the connection of the HM formula (2) with the IAV
model. For that, one needs to transform first the IAV formula, Eq. (22) into
its prior form. This can be done using the following relation due to Li,
Udagawa and Tamura Li84 :

$\varphi^{\mathrm{post}}_{x}(\vec{r}_{x})=\varphi^{\mathrm{prior}}_{x}(\vec{r}_{x})+\varphi^{\mathrm{HM}}_{x}(\vec{r}_{x}).$
(26)

Replacing (26) into Eq. (22) one gets

$\displaystyle\left.\frac{d^{2}\sigma}{dE_{b}d\Omega_{b}}\right|^{\mathrm{IAV,prior}}_{\mathrm{NEB}}$
$\displaystyle=\left.\frac{d^{2}\sigma}{dE_{b}d\Omega_{b}}\right|^{\mathrm{UT}}_{\mathrm{NEB}}+\left.\frac{d^{2}\sigma}{dE_{b}d\Omega_{b}}\right|^{\mathrm{HM}}_{\mathrm{NEB}}$
$\displaystyle+\left.\frac{d^{2}\sigma}{dE_{b}d\Omega_{b}}\right|^{\mathrm{IN}}_{\mathrm{NEB}},$
(27)

with the interference (IN) term

$\left.\frac{d^{2}\sigma}{dE_{b}d\Omega_{b}}\right|^{\mathrm{IN}}_{\mathrm{NEB}}=-\frac{4}{\hbar
v_{a}}\rho_{b}(E_{b})\mathrm{Re}\langle\varphi^{\mathrm{prior}}_{x}|W_{xA}|\varphi^{\mathrm{HM}}_{x}\rangle.$
(28)

Equation (27) represents the post-prior equivalence of the NEB cross sections
in the IAV model, with the RHS corresponding to the prior-form expression of
this model. The first term is just the NEB formula proposed by Udagawa and
Tamura Uda81 , which is formally analogous to the IAV post-form formula (22),
but with the $x$-channel wave function given by
$\varphi^{\mathrm{prior}}_{x}({\vec{r}}_{x})$. The second term in (27) is the
HM formula of Eq. (2) and the last term arises due the interference between
the UT and HM terms.

The HM formula [Eqs. (2) and (3)], is then recovered from (27) by neglecting
altogether the UT term. This result suggests that this HM formula is an
incomplete NEB theory, since the UT term can be very large, or even dominant
Jin15b 111This conclusion was in fact shared by M. Hussein himself who, at
least in recent years, was aware of the incompleteness of the HM formula
Hussein . .

The situation for the eikonal HM formula, Eq. (2.2), is qualitatively
different. In this case, we cannot start from the DWBA formula (27), since the
latter assumes that the auxiliary potential potential $U_{a}$ is a function of
the $a-A$ relative coordinate. This is not the case of the HM choice, Eq. (7).
If $U_{a}$ is allowed to be a function of both ${\vec{r}}_{a}$ and
${\vec{r}}_{bx}$ coordinates (26) becomes:

$\displaystyle\varphi^{\mathrm{post}}_{x}$
$\displaystyle=G_{x}\langle{\vec{r}}_{x}\,\chi_{b}^{(-)}|V_{\mathrm{prior}}|\psi^{3B(+)}_{xb}\rangle+\langle{\vec{r}}_{x}\chi^{(-)}_{b}|\psi^{3B(+)}_{xb}\rangle.$
(29)

But, for the HM choice of the $U_{a}$ potential, Eq. (7), $V_{\mathrm{prior}}$
vanishes identically, and one recovers Austern’s formula, Eq. (20). This shows
that the EHM approximation (2.2) incorporates three-body effects which go
beyond its non-eikonal form (2). In fact, Eq. (2.2) represents the Glauber
limit of the three-body formula by Austern al.. In this regard, one may expect
the eikonal HM formula to be more accurate than its original noneikonal
counterpart, whenever the Glauber approximation is justified (i.e., at high
energies).

##  4 Application to nucleon knockout from 14O

In this section, we present some preliminary numerical results comparing the
HM, EHM and IAV formulae. A more systematic study, including more detailed
observables, such as momentum distributions, is in progress and will be
presented elsewhere.

###  4.1 Practical considerations

Some remarks on the numerical implementation of the different models are in
order. The post-form IAV formula faces the problem of the marginal convergence
of the integrals appearing in the source term of Eq. (23), due to the
oscillatory character of the functions appearing in the initial and final
states. To overcome this problem, several regularization procedures have been
used in the literature. Huby and Mines Hub65  and Vincent Vin68  multiply the
source term by an exponential convergence factor, that damps the contribution
of the integral at large distances. Vincent and Fortune Vin70  use integration
in the complex plane to transform the oscillatory functions into decaying
exponentials. One may also replace the distorted waves
$\chi_{b}^{(-)}({\vec{r}}_{b})$ by wave packets constructed by averaging these
distorted waves over finite energy intervals Jin15 ; Jin15b . The resulting
averaged functions become square-integrable and the source term of Eq. (23)
vanishes at large distances.

Notice that this problem does not arise in the prior form version of this
formula because, in this case, the transition operator ($V_{\mathrm{prior}}$)
makes the source term short-ranged. In Jin15b , a numerical comparison between
the post and prior IAV formulae was performed for the 58Ni(d,p)X reaction and
they were found to be yield almost identical results, provided a
regularization procedure was applied to the post formula.

###  4.2 Numerical results

For these test calculations we consider the one-neutron and one-proton removal
processes taking place in the 14O+9Be reaction at two different incident
energies, namely, 53 MeV/u and 80 MeV/u. Experimental data for this reaction
were reported in Ref. Fla12  and analyzed in terms of the eikonal model as
well as the semiclassical “transfer the continuum” method of Bonaccorso and
Brink Bon88 ; Bon91 , which provides an alternative to the eikonal method for
the case of neutron removal.

In the calculations presented here, we focus on the NEB part of the cross
section. We compare the DWBA IAV [Eq. (22)], the HM [Eq. (2)] and the Eikonal
HM formulae [Eq. (2.2)]. Owing to the aforementioned convergence issues, the
calculations with the IAV method were performed using its prior form formula.
Intrinsic spins are ignored for simplicity. For a meaningful comparison, the
same structure inputs and optical potentials are considered in all these
calculations. In particular, the wavefunction of the removed nucleon is
generated with a Woods-Saxon potential with parameters $R_{0}=3.29$ fm,
$a=0.67$ fm and the depth adjusted to reproduce the proton ($S_{p}=4.6$ MeV)
or neutron ($S_{n}=23.2$ MeV) separation energy, as appropriate. We used
energy-independent neutron-target and core-target potentials derived from the
$t\rho$ and $t\rho\rho$ approximation, respectively, assuming for the nucleon-
nucleon t-matrices the parametrization of Hof78  with nucleon-nucleon cross
sections from Cha90  and imaginary-to-real ratios from Ray79  and densities
extracted from Hartree-Fock calculations with the SkX interaction Bro98 . For
the IAV and HM calculations, the potential between projectile and target is
also needed. This potential has been computed via a folding of the neutron-
target and core-target potentials with the square of the wave-function between
neutron and core in the projectile. For the EHM calculation, the optical limit
of the Glauber approximation was used to derive the required S-matrices from
these potentials.

Figure 2: Integrated nonelastic breakup cross sections for the 9Be(14O,13N)
(top) and 9Be(14O,13O) (bottom) reactions at 80 MeV/u and 53 MeV/u computed
with the IAV, HM and EHM formulae.

The results of these calculations are presented in Fig. 2, in the form of bar
diagrams. The upper and bottom panels correspond to the one-proton and one-
neutron removal reactions. For the neutron removal case, we find that the
three methods give close results. This suggests that the nonorthogonality
term, which is absent in the HM calculation, is small in the case of removal
of a tightly bound nucleon. For the proton removal, the HM result tends to
overestimate the IAV result. This indicates that, at lest in this case, the HM
and nonorthogonality terms interfere destructively (c.f. last term in Eq.
(27)]. Remarkably, the Glauber version of the HM formula is in better
agreement with the IAV result than the original (i.e. non-eikonal) HM formula.
As discussed in the previous section, this result can be interpreted
recognizing that the EHM can be regarded as an approximation to the full
three-body IAV theory and, as such, includes effectively contributions from
the three terms in Eq. (27).

Two main conclusions can be drawn from this preliminary analysis. First, the
non-eikonal HM formula provides an accurate approximation to the more
elaborated IAV formula for deeply bound nucleons, but fails for more weakly
bound nucleons. Second, the EHM formula represents a very good approximation
to the IAV formula both for well bound and weakly bound nucleons and for
energies as low as 50 MeV/u. This result seems to add support to the use of
the EHM formula in the analysis of knockout reactions.

We note that some other approximations are implicit in the derivation of all
discussed formulae, including the IAV one. For example, all these formulae are
based on the spectator assumption for the core fragment. This means that the
latter simply scatters elastically by the target nucleus, but does not
participate in the nucleon-removal dynamics. For example, possible
rescattering effects of the struck nucleon are not taken into account. These
rescattering effects, not considered by any of the presented theories, could
modify the cross section for nucleon removal, whose systematics are currently
under discussion Tos14 . A recent comprehensive review on the subject is given
in Ref. Aum20 .

##  5 Summary and conclusions

In this work, we have reexamined the nonelastic breakup formula devised by
Hussein and McVoy (HM) and extensively employed in knockout studies. We have
shown that this formula can be derived from the Ichimura, Austern, Vincent
(IAV) model, recently revisited and applied by several groups, as well as from
the three-body formula of Austern et al. Aus87 . We have also shown that,
owing to the particular choice of the auxiliary interaction $U_{aA}$, the
eikonal version of the HM formula (EHM) incorporates genuine three-body
effects. These effects are also present in the more general three-body formula
of Austern et al., but in a more complicated way.

Preliminary calculations for the one-nucleon removal in the 14O+9Be reaction
shows that the EHM formula reproduces accurately the results of the IAV model.
For the removal of the strongly bound neutron ($S_{n}=23.2$ MeV), the
noneikonal HM formula is also in good agreement with the IAV result. However,
for the one-proton removal ($S_{p}=4.6$ MeV), the noneikonal HM formula tends
to overestimate the IAV result. It will be interesting to extend these
calculations to other systems and energies to see whether these conclusions
remain.

###### Acknowledgements.

This work has been partially supported by the Spanish Ministerio de Ciencia,
Innovación y Universidades under projects FIS2017-88410-P and
RTI2018-098117-B-C21 and by the European Union’s Horizon 2020 research and
innovation program under Grant Agreement No. 654002. J.L. acknowledges support
from the National Natural Science Foundation of China (Grants No. 12035011,
No. 11975167, and No. 11535004), by the National Key R&D Program of China
(Contracts No. 2018YFA0404403 and 2016YFE0129300), and by the Fundamental
Research Funds for the Central Universities (Grant No. 22120200101). M. G.-R.
acknowledges support from the Alexander von Humboldt foundation.

## References

  * (1) J.E. Escher, J.T. Burke, F.S. Dietrich, N.D. Scielzo, I.J. Thompson, W. Younes, Rev. Mod. Phys. 84, 353 (2012).  DOI 10.1103/RevModPhys.84.353
  * (2) P. Hansen, J. Tostevin, Annu. Rev. of Nucl. and Part. Sci 53(1), 219 (2003).  DOI 10.1146/annurev.nucl.53.041002.110406
  * (3) J. Tostevin, Nucl. Phys. A 682(1), 320 (2001).  DOI 10.1016/S0375-9474(00)00656-4
  * (4) J.A. Tostevin, A. Gade, Phys. Rev. C 90, 057602 (2014).  DOI 10.1103/PhysRevC.90.057602
  * (5) M. Hussein, K. McVoy, Nucl. Phys. A 445(1), 124 (1985).  DOI 10.1016/0375-9474(85)90364-1
  * (6) F. Flavigny, et al., Phys. Rev. Lett. 110, 122503 (2013).  DOI 10.1103/PhysRevLett.110.122503
  * (7) L. Atar, et al., Phys. Rev. Lett. 120, 052501 (2018).  DOI 10.1103/PhysRevLett.120.052501
  * (8) S. Kawase, et al., Prog. Theor. Exp. Phys 2018, 021D01 (2018).  DOI 10.1093/ptep/pty011
  * (9) M. Gómez-Ramos, A. Moro, Phys. Lett. B 785, 511 (2018).  DOI 10.1016/j.physletb.2018.08.058
  * (10) G. Baur, R. Shyam, F. Rösel, D. Trautmann, Phys. Rev. C 28, 946 (1983).  DOI 10.1103/PhysRevC.28.946
  * (11) N. Austern, Y. Iseri, M. Kamimura, M. Kawai, G. Rawitscher, M. Yahiro, Phys. Rep. 154, 125 (1987).  DOI 10.1016/0370-1573(87)90094-9
  * (12) S. Typel, G. Baur, Phys. Rev. C50, 2104 (1994).  DOI 10.1103/PhysRevC.50.2104
  * (13) H. Esbensen, G.F. Bertsch, Nucl. Phys. A600, 37 (1996).  DOI 10.1016/0375-9474(96)00006-1
  * (14) T. Kido, K. Yabana, Y. Suzuki, Phys. Rev. C 50, R1276 (1994).  DOI 10.1103/PhysRevC.50.R1276
  * (15) P. Capel, G. Goldstein, D. Baye, Phys. Rev. C 70, 064605 (2004).  DOI 10.1103/PhysRevC.70.064605
  * (16) A. Kasano, M. Ichimura, Physics Letters B 115(2), 81 (1982).  DOI 10.1016/0370-2693(82)90800-0
  * (17) D.Y. Pang, N.K. Timofeyuk, R.C. Johnson, J.A. Tostevin, Phys. Rev. C 87, 064613 (2013).  DOI 10.1103/PhysRevC.87.064613
  * (18) J. Lei, A.M. Moro, Phys. Rev. Lett. 123, 232501 (2019).  DOI 10.1103/PhysRevLett.123.232501
  * (19) M. Ichimura, N. Austern, C.M. Vincent, Phys. Rev. C 32, 431 (1985).  DOI 10.1103/PhysRevC.32.431
  * (20) B.V. Carlson, R. Capote, M. Sin, Few-Body Syst. 57(5), 307 (2016).  DOI 10.1007/s00601-016-1054-8
  * (21) J. Lei, A.M. Moro, Phys. Rev. C 92, 044616 (2015).  DOI 10.1103/PhysRevC.92.044616
  * (22) G. Potel, F.M. Nunes, I.J. Thompson, Phys. Rev. C 92, 034611 (2015).  DOI 10.1103/PhysRevC.92.034611
  * (23) J. Lei, A.M. Moro, Phys. Rev. C 95, 044605 (2017).  DOI 10.1103/PhysRevC.95.044605
  * (24) G. Potel, G. Perdikakis, B.V. Carlson, M.C. Atkinson, W.H. Dickhoff, J.E. Escher, M.S. Hussein, J. Lei, W. Li, A.O. Macchiavelli, A.M. Moro, F.M. Nunes, S.D. Pain, J. Rotureau, The European Physical Journal A 53(9), 178 (2017).  DOI 10.1140/epja/i2017-12371-9
  * (25) T. Udagawa, T. Tamura, Phys. Rev. C 24, 1348 (1981).  DOI 10.1103/PhysRevC.24.1348
  * (26) M. Ichimura, Phys. Rev. C 41, 834 (1990).  DOI 10.1103/PhysRevC.41.834
  * (27) X.H. Li, T. Udagawa, T. Tamura, Phys. Rev. C 30, 1895 (1984).  DOI 10.1103/PhysRevC.30.1895
  * (28) J. Lei, A.M. Moro, Phys. Rev. C 92, 061602 (2015).  DOI 10.1103/PhysRevC.92.061602
  * (29) M. Hussein.  Private Communication
  * (30) R. Huby, J.R. Mines, Rev. Mod. Phys. 37, 406 (1965).  DOI 10.1103/RevModPhys.37.406
  * (31) C.M. Vincent, Phys. Rev. 175, 1309 (1968).  DOI 10.1103/PhysRev.175.1309
  * (32) C.M. Vincent, H.T. Fortune, Phys. Rev. C 2, 782 (1970).  DOI 10.1103/PhysRevC.2.782
  * (33) F. Flavigny, A. Obertelli, A. Bonaccorso, G.F. Grinyer, C. Louchart, L. Nalpas, A. Signoracci, Phys. Rev. Lett. 108, 252501 (2012).  DOI 10.1103/PhysRevLett.108.252501
  * (34) A. Bonaccorso, D.M. Brink, Phys. Rev. C 38, 1776 (1988).  DOI 10.1103/PhysRevC.38.1776
  * (35) A. Bonaccorso, D.M. Brink, Phys. Rev. C 44, 1559 (1991).  DOI 10.1103/PhysRevC.44.1559
  * (36) G. Hoffmann, et al., Physics Letters B 76(4), 383 (1978).  DOI 10.1016/0370-2693(78)90888-2
  * (37) S.K. Charagi, S.K. Gupta, Phys. Rev. C 41, 1610 (1990).  DOI 10.1103/PhysRevC.41.1610.  URL https://link.aps.org/doi/10.1103/PhysRevC.41.1610
  * (38) L. Ray, Phys. Rev. C 20, 1857 (1979).  DOI 10.1103/PhysRevC.20.1857
  * (39) B. Alex Brown, Phys. Rev. C 58, 220 (1998).  DOI 10.1103/PhysRevC.58.220
  * (40) T. Aumann, et al., Prog. Part. Nucl. Phys. (accepted), arXiv:2012.12553 (2020)