Glauber Description of Reaction Cross Sections with Bogoliubov Many-Body Wave Functions

Introduction

The rapid advancements in accelerator and detection technologies have enabled the study of quantum objects governed by strong interactions through a wide range of observables in nucleus-nucleus collisions. Among the various theoretical frameworks, the Glauber model (GM) remains one of the most widely used tools for interpreting these data [1, 2]. For example, it serves as a foundation for establishing scaling in $pp$, $pA$ and $AA$ collisions and provides realistic initial-state conditions for hydrodynamic simulations of the quark-gluon plasma [3]. Furthermore, the GM is extensively employed to analyze inclusive cross sections, such as reaction, single-nucleon removal, and charge-changing cross sections, revealing exotic nuclear structures like halos, skins, and deformations.

The standard Glauber formalism is based on the eikonal approximation, assuming that projectile nucleons travel along straight-line trajectories and undergo multiple, independent sub-collisions with nucleons in the target [3]. However, the complex multiple-scattering operator is typically simplified into a one-body operator, known as the optical limit approximation (OLA). While convenient, the OLA often washes out crucial structural features, leading to an overestimation of the total reaction cross sections (RCS) [4].

To partially account for the multiple-scattering effects beyond the OLA, several approaches have been developed, including the nucleon-target formalism [5], the modified OLA [6], and the use of modified nucleon-nucleon profile functions [7]. Additionally, higher-order corrections have been investigated through expansions in

terms of density matrices [8, 9]. A particularly effective method for analyzing halo nuclei is the few-body Glauber approach, in which the core and valence nucleons are treated separately. In this framework, the OLA is applied to the integration over the core-target coordinates, while the integration over the valence-nucleon coordinate is performed without any approximation [10, 11]. This method has been extensively applied to analyze experimental RCS and momentum distributions of fragments, allowing the deduction of matter density distributions in halo candidates [12, 13, 14, 15, 16, 17]. To fully incorporate the multiple-scattering processes, in nucleus-nucleus collisions involving projectiles and targets with mass numbers $A_{\mathrm{P}}$ and $A_{\mathrm{T}}$, respectively, a $3 \times(A_{\mathrm{P}} + A_{\mathrm{T}} - 2)$-dimensional integral must be computed. In such cases, Monte Carlo integration (MCI) becomes the only viable option [18, 19].

A key ingredient for the MCI is the many-body wavefunction of the colliding nuclei. Previous applications have relied on phenomenological product states consisting of the harmonic-oscillator single-particle wave functions [19, 20] or the ab initio wave function obtained from the quantum Monte Carlo method [18, 30, 31]. At the same time, modern nuclear density functional theories (DFTs) have achieved remarkable success in describing ground-state and excited-state properties across the nuclear landscape. Yet their connection to reaction observables remains largely restricted to one-body ground-state densities [21, 22, 23, 24, 25, 26]. A consistent treatment that preserves the full many-body structure of DFT wave functions in reaction calculations has been missing.

In the present work, I establish such a microscopic connection. By projecting relativistic Hartree-Bogoliubov (RHB) wave functions onto fixed particle numbers and expressing them in Pfaffian form in coordinate space, I embed correlated many-body states directly into the MCI of the GM. This formulation allows pairing correlations and multiple-scattering effects to be treated on an equal footing for the first time. I first benchmark the method using the nearly spherical $^{16-24}\mathrm{O}$ isotopes and then extend the calculations to the deformed $^{18-28}\mathrm{Ne}$ and $^{20-34}\mathrm{Mg}$ isotopic chains. For the latter systems, I explicitly examine the orientation dependence of the RCS and restore rotational symmetry to obtain the physical $J = 0$ laboratory-frame results. The overall agreement with the experimental data demonstrates the reliability of the present framework for describing high-energy reaction cross sections. The theoretical framework and method are introduced in § 2. The results are described and discussed in Sec. III. Sec. IV contains a summary of the principal results.

Theoretical Framework

The high-energy collision of a projectile nucleus with mass number $A_P$ and a target nucleus with mass number $A_T$ is described within the Glauber model [27]. Under the adiabatic and eikonal approximations, the final state wave function of a projectile and target system $\Phi_f$ is described by the product of the ground-state wave functions of the projectile $\Phi_0^P$ and the target $\Phi_0^T$ nuclei, and the product of the phase-shift functions of a nucleon-nucleon collision $e^{i\chi_{NN}}$ as

where $\boldsymbol{b}$ is the impact parameter. Taking the $z$ direction as the beam direction, $\hat{\boldsymbol{s}}_j^P$ ($\hat{\boldsymbol{s}}_k^T$) denotes the two-dimensional single-particle coordinate operator projected onto the $x$-$y$ plane of the $j$th ($k$th) nucleon from the center of mass of the projectile (target). With this approximation, we only need to evaluate the optical phase-shift function

where the profile function is defined as

The total RCS is evaluated by

where $P(\boldsymbol{b})$ is the reaction probability

The profile function is parametrized as

where $\sigma_{NN}^{\mathrm{tot}}$, $\alpha_{NN}$, and $\beta_{NN}$ are the total nucleon-nucleon cross section, the ratio between the real and imaginary parts of the scattering amplitude at the forward angle, and the slope parameter, respectively. The parameter sets for various incident energies listed in Ref. [28] for proton-proton (pp) and proton-neutron (pn) are employed. The nn (np) are taken to be the same as pp (pn).

Introducing the $A$-body density

where $\boldsymbol{R}=(\boldsymbol{r}_1,\cdots,\boldsymbol{r}_A)$ and $\hat{\boldsymbol{r}}_i$ is the single-particle coordinate operator of the $i$th nucleon from the origin. The complete Glauber amplitude can be written as

where $\tilde{\boldsymbol{s}}_j^P$ ($\tilde{\boldsymbol{s}}_k^T$) denotes the $xy$ component of the $j$th ($k$th) single-particle coordinate from the center-of-mass coordinate of the projectile (target) ($\tilde{\boldsymbol{r}}_j=\boldsymbol{r}_j-\boldsymbol{X}$ with $\boldsymbol{X}=\frac{1}{A}\sum_{j=1}^{A}\boldsymbol{r}_j$). Eq. 8 involves $3(A_P+A_T)$-dimensional integration. Following Refs. [18, 19, 29, 30, 31] I employ MCI to evaluate this multidimensional integration and $\rho_{A_P}^P\rho_{A_T}^T$ is used as the guiding function of the Metropolis iterative algorithm.

For the ground-state wave functions of the projectile $\Phi_0^P$ and target $\Phi_0^T$, products of single-particle functions generated from phenomenological potential [19] or wave functions obtained by the QMC method [18, 29, 30, 31] are used in previous works. In the present work, I use the many-body Bogoliubov wave function generated from the RHB model where the RHB equation

is solved by expanding the Dirac spinors in a harmonic oscillator (HO) basis. $E_k$ is the quasiparticle energy, $(U_k,V_k)^{\mathrm{T}}$ is the quasiparticle wave function, $\lambda$ is the chemical potential and $h_D$ is the single-particle Dirac Hamiltonian, $\Delta$ is the pairing potential. The corresponding Bogoliubov many-body wave functions, in the Thouless representation [32], are given by

where the skew-symmetric matrices $Z=(VU^{-1})^{*}$ and $\lvert0\rangle$ is the true vacuum. $c_i^\dagger$ and $c_k$ represent the single particle creation and annihilation operators. The coordinate-space amplitude for an even number of $A=2n$ nucleons can be expressed as

where

is the field operator. The Bogoliubov vacuum does not have fixed particle number. By projecting the wave functions Eq. (10) onto the $A$ particle sector, I obtain

Applying the fermionic Wick theorem yields

where

and $\phi_k(r)$ is the HO single-particle wave function. This representation encodes two-body correlations beyond a simple Slater determinant. For a deformed nucleus, the intrinsic Bogoliubov vacuum does not have definite angular momentum. The laboratory-frame state with angular momentum $J$ and projection $M$ is [33]

where

and $\hat{R}(\Omega)$ is the rotation operator. $\Omega\equiv(\alpha,\beta,\gamma)$ is the Euler angle. Note $\hat{R}(\Omega)\lvert\Phi_0(Z)\rangle=\lvert\Phi_0[Z^{(\Omega)}]\rangle$ with $Z^{(\Omega)}=DZD^{T}$ and $D_{ji}(\Omega)=\langle j\lvert\hat{R}(\Omega)\rvert i\rangle$, I can write

with

Because $[\hat{N},\hat{R}(\Omega)]=0$, following the same derivation of Eq. (14), I obtain

with

For the axially symmetric and reflection-symmetric even-even nuclei considered in the present work, the integration over the Euler angle $\gamma$ can be carried out analytically. The exact projected coordinate-space wave function consequently involves a two-dimensional (2D) integration over $\alpha$ and $\beta$. Direct evaluation of this 2D calculation is computationally demanding. I therefore employ a one-dimensional (1D) axial projection approximation in which only rotations about the $y$ axis are retained. The approximate projected coordinate-space amplitude is written as

with

To avoid singularity when calculate $Z$ for the case of vanishing pairing, a regularization factor $\eta=10^{-3}$ is added to $\overline{\Delta}_{i\bar{i}}$. Since the Glauber theory is formulated within a non-relativistic kinematic framework and the resulting profile function Eq. (6) is therefore non-relativistic, as a first attempt to apply a relativistic microscopic approach to reaction observables, I neglect the contributions from the small components and retain only the large components to construct the coordinate-space many-body wave function Eq. (14). One should emphasize that, although several previous works have attempted to connect relativistic DFT to reaction observables via the Glauber model, these studies are limited to the OLA and require only the one-body density from relativistic DFT [22, 23, 34].

The relativistic energy density functionals (REDFs) are rooted in quantum hadrodynamics (QHD), in which a nucleus is described as a system of Dirac nucleons interacting via meson exchange through an effective Lagrangian (see, e.g., Ref. [35]). In the mean-field approximation, the theory can be reformulated in terms of equivalent local interactions with density-dependent couplings. The parameters of a given REDF are constrained by nuclear matter properties and further fine-tuned to experimental data of finite nuclei. In the present work, I employ the DD-PC1 parameter set [36]. For the pairing channel, I adopt a separable finite-range interaction adjusted to reproduce the density dependence of the pairing gap in symmetric nuclear matter obtained with the Gogny force [37, 38, 39]. A fully consistent derivation of both the nuclear structure input and the reaction profile function within a unified relativistic microscopic framework remains highly nontrivial. Therefore, the nucleon-nucleon profile function in Eq. (6) is taken from empirical nucleon-nucleon scattering data [28]. Calculations are performed in a basis with 10 oscillator shells. As in Ref. [19], Coulomb breakup contributions are neglected, as they are expected to be small for light systems. I emphasize that the present approach is not a fully microscopic or ab initio framework in which the nuclear structure and reaction ingredients are derived from a common underlying interaction. The term “microscopic” is used here to indicate that the self-consistent many-body Bogoliubov wave functions are incorporated directly into the complete Glauber multiple-scattering calculation, rather than entering the reaction model only through their one-body densities.

Results and Discussions

Spatial configurations in the $A$-body Cartesian coordinate space $(\boldsymbol{r}_1,\cdots,\boldsymbol{r}_A)$ are generated via a Metropolis random walk with a step size of $\delta= 1.0$ fm. The resulting $10^6$ MC configurations, sampled according to the probability distribution of the guiding function are used to perform the multidimensional integration over projectile and target coordinates. The 2D integration over $\boldsymbol{b}$ is performed using a 21st-order Gauss-Kronrod quadrature. To evaluate Eq. (22) for axially deformed nuclei, the integration over $\beta$ is performed using a 14-point Gauss-Legendre quadrature.

As mentioned in Ref. [19], the root-mean-square (RMS) radii can be used to evaluate the precision of the MCI. As a benchmark, the RMS radii of $^{16-24}\mathrm{O}$ isotopes calculated using the MCI with Pfaffian wave functions

are presented in Figure 1 (a). The center-of-mass (C.M.) uncorrected RMS radii agree well with those obtained from RHB calculations, i.e., the deviation is at most $\sim0.2\%$ and can be reduced further by increasing the number of MC configurations. The C.M. corrected RMS radii obtained by MCI of the multidimensional integration

Vector (PDF)

with $\rho_A(\mathbf{R})$ as the guiding function are also plotted in Figure 1(a). The C.M. corrected RMS radii are typically smaller than that of the corresponding C.M. uncorrected value by $\sim0.05$ fm, consistent with the trend reported in Ref. [19]. For the RCS of $^{12}$C+$^{16}$O, Figure 1(b) displays the energy dependence from 40 to 1000 MeV/nucleon. The Pfaffian-based MCI calculation shows excellent agreement with experimental data, particularly at higher energies where the eikonal approximation is most robust. In addition, since the RHB calculations predict vanishing pairing, the one calculated with Pfaffian-based MCI and Slater determinant (consists of single-particle wave functions from RHB calculations but without pairing correlations) based MCI agrees.

The interaction cross sections (ICS) $\sigma_I$ for the oxygen isotopic chain $^{16\text{--}24}$O on a $^{12}$C target at 1.0 GeV/nucleon are analyzed in Figure 2. It is evident that the OLA, which utilizes one-body densities,

Vector (PDF)

consistently overestimates the measured $\sigma_I$ values, in agreement with the results presented in Ref. [22]. In contrast, the MCI results, including both Slater determinant and Pfaffian wave functions-provide a significantly better description of the experimental data. This improvement demonstrates that multiple-scattering effects beyond the OLA are essential for a precise description of reaction observables in these isotopes, as already pointed out in Refs. [5, 6, 7, 8, 9, 15, 43]. For ${}^{16-22}\mathrm{O}$, the differences between our predictions and the experimental data remain within approximately 20 mb, corresponding to at most 2% of the measured values. It should be noted that the predicted $\sigma_R$ values are systematically smaller than the measured $\sigma_I$. Strictly speaking, $\sigma_R$ should be larger than $\sigma_I$ because of contributions from inelastic channels, although the difference is expected to be negligibly small at incident energies around 1 GeV/nucleon [6, 44, 45]. Therefore, the actual uncertainties of the present model may be slightly larger than those inferred from a direct comparison with the experimental $\sigma_I$ data. As demonstrated in Ref. [20], the difference between $\sigma_I$ and $\sigma_R$ is typically no more than 2%–3% of $\sigma_R$. Taking this difference into account, our calculations underestimate the physical reaction cross sections by no more than approximately 5% for the isotopes considered here. Special attention should be paid to ${}^{24}\mathrm{O}$ where the present MCI underestimates the interaction cross sections while the OLA reproduces the experimental values fortuitously which may indicate the need to improve of the nucleon-nucleon effective interactions for the study of neutron-rich nuclei in the mean-field approximation [46].

The extension to Ne and Mg requires particular care because their self-consistent RHB ground states are generally deformed. Figure 3 shows the quadrupole deformations [(a)–(c)] and matter radii [(d)–(f)], calculated with and without pairing for the O, Ne, and Mg chains. The pairing energies obtained from RHB calculations are also shown in Figure 3 [(g)–(i)]. Whereas the oxygen ground states remain nearly spherical, several Ne and Mg isotopes exhibit sizable prolate or oblate deformation. Switching pairing on or off can also change the self-consistent equilibrium deformation. Consequently, a direct comparison between paired and unpaired intrinsic calculations does not isolate a pure pairing effect: it simultaneously changes the pairing field, intrinsic shape, and orientation-dependent reaction geometry.

Vector (PDF)

For an axially deformed intrinsic state, the transverse distribution seen by the target depends on the angle $\beta$ between the intrinsic symmetry axis and the beam direction. This dependence is illustrated in Figure 4 for a prolate configuration of ${}^{20}\mathrm{Ne}$ and an oblate configuration of ${}^{24}\mathrm{Ne}$. For prolate ${}^{20}\mathrm{Ne}$, the cross section is smallest for the tip-on configuration and largest when the symmetry axis is perpendicular to the beam. The trend is reversed for oblate ${}^{24}\mathrm{Ne}$. The intrinsic reaction cross section is therefore not a rotationally invariant laboratory observable. The upper panel of Figure 4 displays the angular Gauss–Legendre weights entering the projections onto $J=0$, 2, and 4. The projected cross sections are calculated using the approximate projected wave functions constructed as coherent superpositions of intrinsic configurations rotated about the $y$ axis, as given in Eq. (22). In the present application, the $J=0$ projection constructs the appropriate laboratory-frame ground state and removes the arbitrary choice of intrinsic orientation. Owing to the opposite orientation dependence of the reaction cross section for prolate and oblate shapes, the projected $J=0$, 2, and 4 states exhibit opposite spin-dependent trends. For the prolate nucleus ${}^{20}\mathrm{Ne}$, $\sigma_R$ decreases with increasing $J$, whereas the ordering is reversed for the oblate nucleus ${}^{24}\mathrm{Ne}$. Nevertheless, the differences among the projected cross sections remain small, at approximately 1%.

Figure 5 and Figure 6 show the RCS for the Ne and Mg chains, respectively. The unprojected Pfaffian calculations correspond to a fixed orientation of the intrinsic symmetry axis relative to the beam. They display pronounced fluctuations that follow the changes in equilibrium deformation shown in Figure 3. These fluctuations arise predominantly because differently deformed intrinsic states present different transverse geometries in the selected orientation. After projection onto $J=0$, the calculated cross sections exhibit much smoother isotopic behavior and generally follow the experimental trends. The projected Pfaffian and projected Slater-determinant results are close for most of the nuclei considered. This indicates that, for the present well-bound systems and the inclusive reaction cross section, the dominant differences between the original intrinsic calculations arise from deformation and orientation rather than from a strong direct sensitivity to pairing correlations. The comparison also illustrates an important requirement for microscopic structure-reaction calculations. Although an intrinsic mean-field state provides a convenient representation of nuclear deformation, it cannot be inserted directly into a laboratory-frame reaction calculation without specifying or restoring its orientation. In particular, differences between paired and unpaired self-consistent calculations cannot be identified with direct pairing effects when the two solutions have different intrinsic deformations. Symmetry restoration or a physically justified treatment of orientation is required before such structural effects can be interpreted.

Vector (PDF)

Vector (PDF)

The calculated RCS are compared with the measured ICS. For the Ne isotopes, the present model underestimates the measured $\sigma_I$ values for ${}^{18}\mathrm{Ne}$ and ${}^{20}\mathrm{Ne}$. The largest discrepancy occurs for ${}^{18}\mathrm{Ne}$ and amounts to approximately 4%. Because the RCS is expected to exceed the ICS by approximately 2%–3%, the calculated $\sigma_R$ for ${}^{18}\mathrm{Ne}$ may underestimate the corresponding physical RCS by as much as approximately 6%–7%. By contrast, the calculations overestimate the measured $\sigma_I$ values for ${}^{26}\mathrm{Ne}$ and ${}^{28}\mathrm{Ne}$ by at most approximately 4%. Once the expected 2%–3% difference between $\sigma_I$ and $\sigma_R$ is taken into account, these results are compatible with the physical RCS at the 1%–2% level.

For the Mg isotopes, the calculated RCS generally reproduce the isotopic trend of the experimental data. The result for ${}^{20}\mathrm{Mg}$ lies below the measured $\sigma_I$ by approximately 5%; after accounting for the expected difference between $\sigma_I$ and $\sigma_R$, the physical RCS may therefore be underestimated by approximately 7%–8%. For the heavier Mg isotopes, the calculated $\sigma_R$ values are generally larger than the measured $\sigma_I$ values. An overestimate of $\sigma_I$ by approximately 2%–4% is largely consistent with the expected inequality $\sigma_R > \sigma_I$. The largest deviation appears for ${}^{34}\mathrm{Mg}$, for which the calculated value exceeds the measured interaction cross section by approximately 6%. Even after allowing for the expected $\sigma_R-\sigma_I$ difference, a residual discrepancy of approximately 3%–4% remains. Nevertheless, considering the experimental uncertainties and the distinction between reaction and interaction cross sections, the present calculations provide an overall reasonable description of the measured Ne and Mg systematics.

To assess the accuracy of the 1D axial projection approximation [Eq. (22)], we have also performed a calculation using the full 2D angular integration for ${}^{20}\mathrm{Ne}$ at an incident energy of 1 GeV/nucleon. The resulting RCS is 1092 mb, compared with 1109 mb obtained using the 1D approximation at the same energy, corresponding to a relative difference of approximately 1.6%. This small difference indicates that the 1D approximation provides a reliable description of the total RCS for the system considered here and does not affect the main conclusions of the present study.

Summary

In summary, I have developed a microscopic framework that embeds particle-number-projected Bogoliubov many-body wave functions directly into the MCI of the complete Glauber multiple-scattering operator. The coordinate-space Pfaffian representation makes it possible to sample the many-body probability distribution of an RHB state without reducing the nuclear input to its one-body density. Applications to O isotopes demonstrate that complete multiple-scattering calculations reduce the calculated RCS relative to the OLA and generally improves agreement with the measured ICS. Applications to the deformed Ne and Mg chains show that RCS evaluated from fixed-orientation intrinsic states can exhibit large geometrical variations. These variations are strongly correlated with the self-consistent deformation. Restoring rotational symmetry produces smooth $J=0$ laboratory-frame systematics and provides a consistent basis for comparing microscopic calculations with unpolarized reaction measurements.

The present work establishes the Pfaffian-based Glauber MCI as a general interface between Bogoliubov nuclear structure calculations and high-energy reaction observables. Future applications to weakly bound systems and more exclusive observables may provide direct access to correlations encoded in the many-body wave function beyond the one-body density. The significant discrepancies between the calculations and the experimental data for nuclei at or near the drip-line, such as ${}^{24}\mathrm{O}$, ${}^{18}\mathrm{Ne}$, and ${}^{20}\mathrm{Mg}$, may indicate limitations of the mean-field description in weakly bound systems, including deficiencies in the effective nucleon-nucleon interaction and the treatment of continuum coupling.

Acknowledgments

This work has been supported by the National Natural Science Foundation of China under Grant No. 12005107. Calculations were performed on the Guangxi Key Laboratory of Nuclear Physics and Technology High Performance Computing Platform.