Neutron scattering is a direct, Coulomb-insensitive tool for probing the nuclear environment. The simplest neutronnucleus interaction quantity is the neutron total cross section, σ tot , which provides information about nuclear size and the ratio of elastic to inelastic components of nucleon scattering. Additionally, σ tot data are thought to be tightly correlated with a variety of structural nuclear properties of great interest including the neutron skin of neutron-rich nuclei [1] and thus the density dependence of the symmetry energy L, an essential equation-of-state input for neutron-star structure calculations [2][3][4].
In the crude "strongly absorbing sphere" (SAS) approximation, where a target nucleus absorbs incident neutrons passing within a nuclear radius, σ tot depends solely on the target nucleus size and the energy of the incident neutron:
where R = r 0 A 1 3 and λ is the reduced wavelength of the incident neutron with energy E in the center of mass [5,6]. While on average, experimental σ tot data comport with this naive model, the most prominent feature of experimental σ tot data is the oscillatory behavior centered about the average * Present address: Lawrence Livermore National Laboratory, Livermore, CA 94550, USA; pruitt9@llnl.gov † Present address: Department of Physics, University of Massachusetts-Lowell, Lowell, MA 01854.
‡ Present address: TRIUMF, Vancouver, BC V6T 2A3, Canada.
of Eq. ( 1), visible in Fig. 1. Peterson [7] interpreted these oscillations as the result of a phase shift between neutron partial waves passing around the nucleus (thus undergoing no phase shift) and waves passing through the nuclear potential, where they are refracted and exhibit a retardation of phase (an illustration is available in [6]). This explanation was termed the "nuclear Ramsauer effect" by Carpenter and Wilson [13] based on the analogous effect seen in electron scattering on noble gases. Following Angeli and Csikai [14], this explanation can be incorporated by imbuing the strongly absorbing sphere relations with a sinusoidal term:
where ρ = e -Im( ) and δ = Re( ), being the phase difference between a partial wave traveling around and traveling through the nucleus. The large amplitude of the oscillations suggests that elastic scattering accounts for a significant fraction of the total cross section, in turn implying a larger meanfree path for neutrons through the nucleus than might otherwise be expected in the absence of Pauli blocking [15,16]. If we approximate the nucleus with a real spherical potential of radius R and depth U , the total phase shift δ is
where C = 4 3 R is the average chord length through the sphere [14]. Rearranging Eq. ( 3) in terms of A and E and discarding [8][9][10][11][12]. Predictions for σ tot given by the "strongly absorbing sphere" (SAS) model [Eq. (1)] are shown as thin dashed lines for each nucleus. Regular oscillations about the SAS model are visible as is the trend for the oscillation maxima and minima to shift to higher energies as A is increased. leading constants yields
This form reveals an important relation: as A is increased, to maintain constant phase δ, E must also increase [6,7]. This is contrary to a typical resonance condition where an integer number of wavelengths are fitted inside a potential; in that case, to maintain constant phase as A is increased, E must be decreased. Thus these σ tot oscillations have been referred to as "antiresonances" or "echoes" [6,17]. Other authors [18] have exposed weaknesses in Angeli and Csikai's interpretation of Eq. ( 2) and have provided a more general semiempirical equation for σ tot . However, Eq. ( 2) is a valuable starting point for connecting σ tot with the depth and shape of the nuclear potential as experienced by neutrons. By including additional surface, spin-orbit, and other terms, optical models (OMs) have been used to successfully reproduce the general features of all manner of single-nucleon scattering data across the chart of nuclides up to several hundred MeV [19][20][21]. However, despite the excellent agreement with experiment, OMs involve the interaction of many partial waves with many sometimes-opaque terms in the potential, complicating intuitive understanding of the underlying physics at play. In particular, the isovector components of optical potentials are quite difficult to constrain as they depend on both proton and neutron scattering data, one or both of which are often unavailable. For example, when Dietrich et al. conducted an analysis of neutron total cross section differences between W isotopes, including standard isovector terms in their optical potential worsened the reproduction of experimental relative differences, an illustration of how poorly these isovector components are known [22].
With these considerations in mind, our present goal is twofold: first, to provide new isotopically resolved σ tot data useful for identifying the dependence of optical potential terms on nuclear asymmetry; and second, to conduct a dispersive optical model (DOM) analysis of these new σ tot data along with a large corpus of scattering and bound-state data to extract veiled structural quantities (e.g., neutron skin thicknesses and spectroscopic factors, or SFs) for several cornerstone, closed-proton-shell nuclei. Key findings of this DOM analysis are presented in the companion Letter [23].
By scattering secondary radioactive beams off of hydrogen targets in inverse kinematics, proton-scattering experiments are possible even on highly unstable nuclides. Because neutrons themselves must be generated as a secondary radioactive beam, neutron-scattering experiments are restricted to normal kinematics and σ tot measurements are possible only for relatively stable nuclides that can be formed into a target. At present, σ tot measurements above the resonance region on nuclides with short half-lives (shorter than the timescale of days) are technically infeasible for this reason, though a handful have been carried out on samples with half-lives in the tens to thousands of years [10,24,25].
Traditionally, σ tot measurements have relied on analogelectronics techniques for recording events, techniques that suffer from a large per-event dead time of up to several µs. For a typical analog intermediate-energy σ tot measurement with dozens or hundreds of energy bins, achieving statistical uncertainty at the level of 1% requires a thick sample to attenuate a sizable fraction of the incident neutron flux. If cross sections are in the 1-10 barn range, this means sample masses of tens of grams [8,12]. Producing an isotopically enriched sample of this size is often prohibitively expensive. As a result, there is a dearth of σ tot data on isotopically resolved targets from 1-300 MeV, even for closed-shell isotopes of special importance such as 3,4 He, 18 O, 64 Ni, 112,124 Sn, and 204,206 Pb (see Fig. 1.3 in [26]).
Recent developments in wave-form-digitizer technology have made it possible to reduce the per-event dead time by an order of magnitude or more, enabling a corresponding reduction in the necessary sample size. In 2008, we embarked on a campaign of σ tot measurements on isotopically enriched samples using these new technical capabilities, starting with 40,48 Ca from 15 E lab 300 MeV [27]. The data from that measurement were incorporated into several DOM analyses [28][29][30] that yielded proton and neutron SFs, charge radii, and initial estimates of the neutron skins [1] for these nuclei.
Here we significantly expand on that effort by providing σ tot results for the important closed-shell nuclides 16,18 O, 58,64 Ni, and 112,124 Sn. We also present a measurement on a very thin sample of the naturally monoisotopic 103 Rh to demonstrate that σ tot experiments over a broad energy range using only a few grams of material are feasible.
All σ tot measurements were carried out at the 15R beamline at the Weapons Neutron Research (WNR) facility of the Los Alamos Neutron Science Center during the 2016 and 2017 run cycles. Our experiment was modeled on previous σ tot measurements at WNR [8,12,27]. At WNR, broad-spectrum neutrons up to ≈700 MeV are generated by impinging proton pulses onto a water-cooled, 7.5-cm-long tungsten target (Fig. 2). Before the beam enters the experimental area, a permanent magnet deflects all charged particles generated by the proton pulses, allowing only neutrons and γ rays to reach the experimental area. At the entrance to the experimental area, the beam was collimated to 0.200 inches using steel donuts with a total thickness of 24 inches. In addition, the γ -ray content of the beam was suppressed using a plug of Hevimet (90% W, 6% Ni, 4% Cu by weight) at the upstream entrance of the collimation stack. After collimation, the beam passed successively through a flux monitor, the sample of interest, a veto detector, and finally the time-of-flight (TOF) detector approximately 25 meters from the neutron source. All detectors consisted of BC-400 fast scintillating plastic mated with photomultiplier tubes (PMTs) and encased in either a plastic or an aluminum housing. The flux monitor and veto detector each had scintillator thicknesses of 0.25 inches and the TOF detector had a scintillator thickness of 1 inch. Signals from all detectors and the target changer were relayed to a 500-MHz CAEN DT-5730 wave-form digitizer running custom software. To increase light collection and thus lower employable thresholds and also to improve time resolution, the TOF detector used two PMTs (one left, one right) mated to the same plastic scintillator and the PMTs' signals were summed before digitization.
The particular neutron beam structure at WNR dictates the energy range achievable for σ tot measurements (Fig. 3). Proton pulse trains, called "macropulses," are delivered to the tungsten target at 120 Hz. Each macropulse consists of ≈350 individual proton pulses, called "micropulses," spaced 1.8 µs apart. Each micropulse consists of a single proton packet that generates γ rays and neutrons within a tight temporal-spatial range. As neutrons from this micropulse travel along the beam path, high-energy neutrons separate in time from lowerenergy neutrons so that neutron energy can be determined by standard TOF techniques (see [31] for details). Because the γ rays and high-energy neutrons from later micropulses can overtake slower neutrons from an earlier micropulse, the distance of the TOF detector from the neutron source determines both the minimum neutron energy that can be unambiguously resolved and the maximum instantaneous neutron flux, critical to correcting for per-event dead time.
A programmable sample changer with six positions was used to cycle each sample into the beam at a regular interval of 150 seconds per sample. Once per macropulse, an analog signal from the sample changer was recorded to indicate its current position. The flux monitor was used to correct for variations in beam flux between macropulses. The veto detector suppressed events from charged-particle production in the samples and in air along the flight path.
Custom digitizer software was used to run the digitizer in two complementary modes, referred to as "DPP mode" and "wave-form mode." In DPP mode, triggers were initiated by the digitizer's onboard peak-sensing firmware. For each trigger, several quantities were recorded: the trigger time stamp, two charge integrals over the detected peak with different integration ranges (32 ns for the short integral, 100 ns for the long integral), and a 96-ns portion of the raw digitized wave form, referred to as a "wavelet." DPP mode was used for the vast majority of the experiment and accounts for ≈99% of the total data volume. In wave-form mode, the digitizer performs no peak sensing and was externally triggered. Upon triggering, the trigger time stamp and a very long wavelet (60 µs) were recorded. While wave-form mode data account for only ≈1% of the total data, the instantaneous data rate is much higher than in DPP mode because hundreds of µs of consecutive wave-form samples are stored. Roughly once every three seconds, the digitizer was switched to wave-form mode for one macropulse, then switched back to DPP mode as quickly as possible (10-40 ms, depending on run configuration). Except for the O and Rh samples, all samples were prepared as right cylinders 8.25 mm in diameter and ranging from 10-27 mm in length (see Table I for sample characteristics). For each element studied, a natural-abundance sample was also prepared as were two natural C samples and a natural Pb sample, useful for benchmarking against literature data. The samples were inserted into styrofoam sleeves and seated in the cradles of the sample changer. This design minimizes the amount of nontarget mass proximate to the neutron beam path. Our samples were generally much smaller than those used in previous measurements; for example, the Ni and Sn samples used in [8,12] had areal densities of 1.515 and 0.5475 mol/cm 2 , respectively, 12.7 and 6.5 times larger than for our Ni and Sn samples.
The O isotopes were prepared as water samples to increase the areal density of atoms and for ease of handling. Each water sample was contained by a cylindrical brass vessel with thin brass end caps (0.002 inches), and an empty brass vessel served as the blank. 16,18 O cross sections were calculated by subtracting the well-known H cross section from the raw H 2 O results. We used H σ tot data sets from Clement et al. [32] and FIG. 4. The effects of timing corrections on the γ -ray peak of a typical run are shown. The uncorrected spectrum is shown in black, the spectrum after correction with our software CFD is shown in blue, and the spectrum after correction with both our software CFD and γ averaging is shown in magenta. For this run, the final γ -ray peak FWHM after both corrections is 0.866 ns, comparable to the precision we achieved in our Ca study [27], which also employed γ averaging.
Abfalterer et al. [12], which together cover the range 0.5 E n 500 MeV and are in excellent agreement where their energy ranges overlap. In light of the additional uncertainty inherent to this subtractive σ tot determination, we prepared a deuterated water sample, from which the literature σ tot for D 2 could be subtracted, to serve as an additional cross-check. Due to the poor machining properties of Rh, the 103 Rh sample was prepared by purchasing and stacking a series of thin disks rather than by manufacturing a fused cylinder. These disks were held in place by a cylindrical plastic case with open ends.
The quantity of interest, σ tot , is related to the flux loss through a sample by
or, equivalently,
where I 0 is the neutron flux entering the sample, I t is the neutron flux transmitted through the sample without interaction, ρ A is the number density of nuclei in the sample, and is the sample length. For thin or low-density samples, flux attenuation through the sample will be small (e.g., 13% for our Ni samples at 100 MeV) and a large number of counts will be required to determine the cross section to high precision.
Two postprocessing steps were used to improve TOFdetector timing resolution (see Fig. 4). First, the wave form for each TOF-detector event was passed through a software constant-fraction discriminator (CFD) logic, improving precision by a factor of two. Second, a γ -ray-averaging procedure (cf. [27]) was used to improve the precision of each micropulse start time. The final corrected TOF resolution (taken as the FWHM of the γ -ray peak in the TOF spectra) ranged from 0.60-0.90 ns over the series of σ tot measurements. This is comparable to the resolution from our digitizer-mediated σ tot measurement on Ca isotopes in 2008 [27]. For context, for a 100-MeV neutron and a TOF detector distance of 25 meters, a TOF uncertainty of 0.80 ns translates to an energy resolution of ≈900 keV. For neutrons below ≈20 MeV, the TOF time resolution worsens because the traversal time through the 1-inch thickness of the TOF detector becomes non-negligible. However, because the TOF of these neutrons is already very long (several hundred ns or longer) the relative energy resolution ( E E ) is superior at low energies. As an example from one of our runs, a 5-MeV neutron with a 0.82-ns detector-traversal time and an inherent TOF resolution of 0.80 ns has an energy uncertainty of 13 keV. These energy uncertainties have been propagated through subsequent analysis into our σ tot results below.
Calculating the neutron energy requires knowledge of the flight-path distance to high precision. We determined this distance by calculating putative σ tot data for nat C from 3-15 MeV from our measurement and comparing the resonance peaks in this region with high-precision literature data sets. From this study, the mean TOF distance was determined as 2709 ± 1 cm for the Ni and Rh run configuration and 2554 ± 1 cm for the Sn and O run configuration.
Before cross sections could be tabulated, the per-event dead time had to be modeled and corrected for. Because events are not processed instantaneously, there is a brief period after each trigger during which the digitizer is busy processing that trigger. Any newly arriving events in this period will be ignored, privileging events arriving earlier and thus distorting TOF spectra and resulting cross sections. This busy period is referred to as the "analytic" or "per-event" dead time and can be corrected for according to standard techniques [31]. An additional complication is the possibility of flux variation between micropulses. If there is no variation, the fraction of time that the digitizer is dead for a given time bin i can be calculated [31]:
where N is the number of time bins in the micropulse, R x is the rate of detected events per micropulse in bin x, and P j is the probability that the digitizer is still busy from a trigger j bins ago. If the variation in beam flux is significant, a more advanced formula can be used; however, an examination of our flux-per-micropulse data showed very little flux variance across macropulses, except during the first 10% of the micropulses within each macropulse. In the final analysis we discarded these first 10% and used the simpler Eq. ( 7) to calculate the dead-time fraction.
To model the experimentally observed probability-dead, P j , we fitted a logistic function to the observed spectrum for time differences between consecutive events (Fig. 5). For a given bin i, the fraction of time that the digitizer is dead, F i , is a discrete convolution of the measured TOF spectrum with P j . Note that except for the first and last micropulses in a macropulse, all micropulses are consecutive, so dead-time effects can "wrap around" from the end of one micropulse to the next. For these wrap-around contributions (that is, j > i), the (mod N) term ensures that the bin referred to by ij is non-negative.
Because trigger processing is done in firmware onboard the digitizer, the per-event dead-times affecting our measurement were reduced to between 150-230 ns. After we calculated the average probability-dead for each time bin, the total number of events detected in that bin, N d [i], could be corrected to recover the true number of events that would have been detected in the absence of a per-event dead time:
where M is the total number of micropulse periods. At large TOFs (low energies) the correction is as low as a few percent, but at small TOFs (high energies), the digitizer is often still dead from the γ -ray flash and high-energy neutrons. In this regime the correction can be quite large (≈20% for our Ni/Rh runs, and ≈40% for our Sn/O runs). Still, the corrections needed for our measurement are far smaller than the typical analytic dead-time corrections required with the dead-time mitigation scheme of previous analog measurements [8,12]. In addition to analytic dead time, there is an additional dead-time effect associated with digitizer readout to the data acquisition computer (DAQ). During data collection, each pair of digitizer channels shares a common buffer for storing events. After several seconds of acquisition, the digitizer begins readout at which time the acquisition is paused and buffer contents are read out to the DAQ. However, because each buffer is independently read out to the DAQ, it is possible that buffers could be emptied and readied for new acquisition at slightly different times (10-40 ms apart), and a mismatch could develop between the number of macropulses seen on FIG. 6. TOF spectra after the analytic dead-time correction and the veto and integrated charge gating for the blank sample (in red) and the nat C sample (in blue), from the Ni/Rh experiment. The γ -ray peak is visible as a sharp spike at 90 ns, followed by the highestenergy neutrons at 130 ns. different channels. Such run-time interactions between the firmware and USB traffic of the DAQ were difficult to characterize, but we estimate that they might cause a systematic error of a few tenths of one percent in the number of macropulses seen by different channels, depending on the user-defined threshold and the buffer size. This effect could contribute to the discrepancy at the highest energies (>100 MeV) between our results and past analog-enabled measurements.
During analysis, it was noted that occasionally (1 in 400 macropulses), one or two adjacent macropulses would have an abnormally small number of events. The frequency of these "data dropouts" was similar to the rate of switching between DPP and wave-form modes; we suspect it is related to edge case behavior right before or after a mode switch. To mitigate this issue, we threw out any macropulse that had less than 50% of the average event rate in either the flux monitor or TOF detector channel.
After applying these corrections, the veto and integrated charge gates were applied to all events and surviving events were populated into TOF spectra (Fig. 6). Next, room background was subtracted (responsible for 0.1% to 1% of event rate, depending on energy) and spectra were mapped to the energy domain.
From these energy spectra, the raw cross sections were calculated, binwise, as follows:
where I 0 /I s is the ratio of counts in the energy spectra between the blank and sample, and M s /M 0 is the ratio of counts in the monitor detector between the sample and blank (for flux normalization). Finally, two isotope-dependent corrections were applied to the raw cross sections. First, because the blank sample contains air and not vacuum, the cross section of air must be added to each sample's cross section. Second, the cross section for 64 Ni was corrected for the isotopic enrichment of our sample (92.2%) using our measured nat Ni cross section. All other isotopes were sufficiently pure such that the impurity correction was negligible.
To validate our analysis, we first benchmarked our σ tot measurements of natural samples ( nat C, nat Ni, nat Sn, and nat Pb) against the high-precision data sets on natural samples from [8,12] (Fig. 7). Our natural sample results are in excellent agreement with these previous results from 3-100 MeV and show slight deviation above 100 MeV (a relative difference of up to 5% at 300 MeV), suggesting a small systematic error at high energies in one or both approaches when the instantaneous neutron flux is highest. As an additional diagnostic, we compared σ tot results from our long and short natural carbon targets and found excellent agreement, within 1% throughout the measured energy domain.
Extracting the 16,18 O σ tot required subtraction of the wellmeasured σ tot for H. To better characterize the additional systematic uncertainty associated with this subtractive analysis, we subtracted our measured values for 16 O neutron σ tot from our raw D 2 O and H 2 O data and calculated the D-to-H relative difference. A comparison of our D-to-H relative difference with that of [33] is shown in Fig. 8. Our results differ systematically from the previous (analog) measurement by 2%-3% throughout the energy range, comparable to the 2% systematic difference between our final 16 O neutron σ tot results and those of [12]. The size and uniformity of these Except for the already well-measured 16 O, our new data significantly extend knowledge of the neutron σ tot for each sample. In the cases of 18 O, 58 Ni, 103 Rh, and 124 Sn, almost no previous data were available above 20 MeV. Our new data are in good agreement with the previous measurements where available. In the cases of the rare isotopes 64 Ni and 112 Sn, data were available at only one energy, 14.1 MeV, from a study from more than 50 years ago [34] and our measurement is in excellent agreement, within 2%-3%.
Our results for relative differences between isotopic pairs 16,18 O, 58,64 Ni, and 112,124 Sn are shown in Fig. 11. For 16,18 O [Fig. 11(a)], the purely isoscalar SAS model [Eq. ( 1)] grossly reproduces the relative difference below 100 MeV, but fails completely above 100 MeV. Near 200 MeV, the 18 O σ tot crosses over that of 16 O resulting in a negative relative difference, in keeping with the Ramsauer-logic expectation of Eq. ( 2) that σ tot oscillation minima shift to higher energies as A is increased. In the relative-difference subfigures for 58,64 trend in Sn isotope-shift data [35] is also shown for reference and underpredicts the relative differences. In the DOM analyses presented below, we fitted only absolute σ tot data and did not directly fit these relative differences. Still, the relative differences between our individual DOM fits for 58,64 Ni and 112,124 Sn (black dasheddotted lines) show overall agreement with the experimental FIG. 9. Neutron σ tot for 16,18 O, 58,64 Ni, and 112,124 Sn: Our results and literature data. In the upper three panels, our digitizer-measured isotopic results are shown in red and corresponding analog-measured literature data [8,34,[36][37][38][39][40][41][42] are shown in blue. The data for 18 O have been shifted up by 1 barn for visibility. The lower three panels show residuals between our data and the literature data shown in the upper panels. relative differences, especially for the Sn relative difference. For the 16,18 O relative difference, there is an obvious phase mismatch between the oscillations of DOM calculation and the experimental data. This mismatch is symptomatic of a slight DOM overestimation of the 16 O radius (0.02 fm), which nudges the DOM-calculated 16 O σ tot rightward so that the 18 O crossover occurs at too low an energy. As was noted by Dietrich et al. in their study of σ tot relative differences in W isotopes, a simultaneous OM analysis along the entire isotopic chain, as in [28], may be required to realize the full isovector-constraining power latent in the relative differences.
The DOM is a phenomenological Green's function framework enabling a simultaneous and self-consistent analysis of nuclear structure and reaction data. An essential feature of the DOM is the enforcement of a dispersion relation between the complex components of the self-energy across the entire energy domain, allowing structural data from below the Fermi energy (e.g., charge densities, bound levels) to help constrain the potential above, and data from above the Fermi energy (e.g., elastic, reaction, and total cross sections) to help constrain the potential below. Using our new σ tot data for 16,18 O, 58,64 Ni, and 112,124 Sn, we performed a simultaneous fit on each isotopic pair and also revisited 40,48 Ca and 208 Pb. Compared to previous DOM analyses [1,28,29,43], we employ an updated version of the DOM that has been generalized for use with any combination of near-spherical even-even nuclei. Partial occupation of neutron open shells, as for the neutron d 5/2 valence shell in 18 O, is accommodated using the level's energy E and the pairing parameter :
where B(N, Z ) is the binding energy of the nucleus with N neutrons and Z protons. Occupation for the level is split into upper (n + ) and lower (n -) components:
where χ ≡ E -F , s ≡ (χ 2 + 2 ) 1 2 . Only the lower (occupied) component is included in calculations of bound-state quantities (e.g., total particle number, binding energy).
In the appendices, we provide the functional forms used to define the potential (Appendix A), optimized parameter values with uncertainties (Appendix B), and figures showing the quality of the DOM reproduction to each experimental data set (Appendix C). The other major methodological difference is the use of Markov chain Monte Carlo (MCMC) for parameter optimization, discussed below.
For additional details on the underlying DOM formalism, see [44,45]. To calculate cross sections from the self-energy, the standard R-matrix approach was used [46]. Except where indicated, experimental data used for fitting are the same as in [26]. To situate the reader, we describe the corpus of experimental data and DOM results for 16,18 O in full detail. The experimental data used and fit quality for 40,48 Ca, 58,64 Ni, 112,124 Sn, and 208 Pb are similar in quantity and quality and only key differences are noted. For systematics of neutron skins and binding energies, see the companion Letter [23].
A. 16
For protons, twenty-eight differential elastic cross section data sets and twenty analyzing power data sets from 10-200 MeV were incorporated. Only three proton reaction cross section data sets, ranging from 20-65 MeV, were available. As an added constraint, we used systematic trends from the comprehensive proton σ rxn review of Carlson [47] to generate proton σ rxn pseudodata from 70-200 MeV, which were included in the fit. These pseudodata are shown as gray open symbols in the proton σ rxn figures in Appendix C. For neutrons, ten differential elastic cross section data sets from 10 MeV to 95 MeV, a single neutron reaction cross section data point at 14 MeV, and our newly measured σ tot results for 16 O were included. In all, over sixty experimental nucleon scattering data sets were used to constrain the 16 O parameters.
In addition to nucleon scattering data, several sectors of bound-state data were included in the fit. Neutron (proton) 0p 1/2 and 0d 5/2 single-particle level energies were assigned according to the nucleon separation energies of 16 O and 17 O isotopes ( 16 O, 17 F isotopes) [48]. Charge density distributions were taken from the compilation of [49]. Since the time of that compilation, new experiments (particularly muonic-atom measurements) have improved the precision of many rootmean-square (rms) charge radii by roughly an order of magnitude [50]. To account for these improved data, we rescaled the distributions from [49] to recover the updated rms charge radii while still conserving particle number. We also fitted directly to the updated rms charge radii of [50]. Because the DOM self-energy does not necessarily conserve particle number, we included the "experimental" proton and neutron numbers of FIG. 11. 16,18 O, 58,64 Ni, 112,124 Sn neutron σ tot relative differences from our measurement. In each panel, the colored bands indicate regions of 1σ uncertainty due to target thickness imprecision (blue) and from both target thickness and statistics (red). The gray dashed lines show the prediction for the σ tot relative difference per the strongly absorbing sphere (SAS) model of Eq. ( 1), which assumes a simple A eight as part of the fit. Lastly, the total binding energy of 16 O from [48] was included as a constraint.
Numerous proton elastic scattering data for 18 O were available from the EXFOR database. Twenty-eight proton elastic differential cross sections were included ranging from 10-200 MeV. Unfortunately, no proton reaction cross section data were available at all in the relevant range of 10-200 MeV. As with 16 O, we generated proton reaction cross section pseudodata from systematic trends in [47] from 70-200 MeV. On the neutron side, two differential elastic cross section data sets were included, at 14 and 24 MeV, but no analyzing powers were available. One datum for the neutron reaction cross section, at 14.1 MeV, was incorporated as well. Our σ tot results for 18 O were the sole neutron total cross section data used in the fit. The energies of the proton and neutron 0p 1/2 and 0d 5/2 single-particle levels were assigned according to the same procedure used for 16 O.
Unlike 16 O, for 18 O, no charge density distribution was available from [49]. To approximate it, we rescaled the charge density distribution used for 16 O to give the 18 O rms charge radius of [50] while preserving eight units of charge. As with 16 O, we also fitted to the experimental rms charge radius directly, to the particle numbers N and Z, and the total binding energy.
Several aspects of the DOM potential make optimization challenging. Even with the reduced number of potential parameters used in this work (42 for 208 Pb and 43 for all other pairwise fits) compared to past DOM studies (for example, 60 or more in [1]), we found that classical gradient-descent methods were inappropriate for reliably searching the parameter space. A recent study [51] systematically compared Bayesian optical model optimization techniques to frequentist ones, the type almost universally used in previous analyses, and found that traditional algorithms may be overconfident in their parameter estimation. To avoid these problems, we used the affine-invariant MCMC library, EMCEE [52], for optimization and uncertainty characterization. For an in-depth introduction to applied MCMC, see [53].
In the ensemble-sampling approach, several hundred "walkers" are first randomly initialized in parameter space for each isotopic system to be fitted. At each subsequent step t during the random walk, each walker's position is updated from x t → x t+1 either by accepting a new position x with probability
or by remaining in the same position x with probability 1 -p( x → x ). New positions are proposed according to the stretch-move proposal distribution of [54] (for our stretch move scaling, we used α = 1.3 instead of the default α = 2.0, which improved the typical acceptance fraction from around 5% to 15%). In Eq. ( 12), the utility of a parameter vector conditional on the experimental data U ( x|D) was defined according to Bayes's rule (omitting the evidence term):
where D is the full set of constraining experimental data.
The parameter prior distribution P( x) was specified as uniform over a physically reasonable range for each parameter. For example, the diffusenesses of all Woods-Saxon potential geometry terms were restricted to 0.4-1.0 fm. Other more sophisticated choices for the prior distribution (e.g., broad truncated Gaussians) were tested and had little impact on the resulting posterior distributions. The likelihood function was defined as a least-squares function over all data sectors d: the model error increased linearly with respect to the scattering angle in the center-of-mass frame with units of % per degree. nlj are the single-particle energies for valence nucleons as calculated from separation energies in [48]. r rms is the root-meansquare charge radius and ρ q is the charge density distribution. Appendix A shows the parameter definitions and prior distributions used in the present analysis.
Due to the choice of functional form and finite model basis size, DOM predictions for nuclear observables suffer from inherent model error. For example, many previous OM analyses tend to easily reproduce low-angle experimental dσ d data taken at lower scattering energies but are increasingly discrepant with the data at high energies and at backward angles, where the predicted cross sections may differ from experimental results by an order of magnitude or more. This discrepancy indicates a deficiency in the potential form of the OM; ignoring it can lead to drastic underestimation of variances of extracted quantities. In this investigation, we found that the inclusion of reasonable model discrepancy terms in our utility function improved the visual fit to experimental data while broadening parameter uncertainties, in keeping with the methodological findings of [55]. Table II shows the model error terms we used for each data sector. We assigned model error for each data set according to how well preliminary fits could reproduce differing regions of each data sector, the flexibility of the functional forms, and intuition from the successes and failures of past OM analyses. In principle, the form of these model error terms could also be treated as random variables to be sampled over during MCMC, but due to computational limitations and the already-challenging size of the DOM parameter space, we elected to fix the model error terms. After N samples have been taken from the posterior distribution, a subset can be used to estimate the true parameter distributions, and physics results calculated for each sample. Ensuring that this subset is representative of the true posterior is discussed in the next section.
Following [52] we attempted an autocorrelation analysis to test for convergence and estimate the number of independent samples we had collected for each nucleus. Because of computational limitations on the number of walkers and steps used to approximate the posteriors, posterior estimation involves a finite MCMC sampling error. The integrated autocorrelation time for a physics feature f , denoted τ f , represents the number of steps required for a walker to produce a new, decorrelated posterior sample for the feature that is independent of the previous independent sample. In an ideal MCMC analysis, τ f could be accurately computed for each physics quantity and the MCMC sampling error could be robustly estimated. In practice, we found this to be computationally infeasible for the DOM parameter space. For example, in preliminary analysis of 18 O, we were able to perform N = 31 000 steps for each of 336 walkers (more than 100 000 CPU hours in total). Over this domain, we calculated the integrated autocorrelation time for each potential parameter p, denoted τ p , to be roughly 2800 steps. Assuming an N > 100τ p rule-of-thumb condition for convergence of the τ p estimate near its true value, the decorrelation time appears to be extremely long. In other words, from τ p alone, we could not exclude the possibility that the parameters had not yet fully "settled" in the region of their optimal values and begun independent sampling of the parameter posteriors. We note that the true τ f could be considerably smaller than τ p due to the highly correlated nature of DOM parameter space.
To proceed, we applied several commonsense tests to judge whether our parameter and extracted-quantity estimates were accurate. First, we sampled as long as possible and used as many parallel walkers as possible, given our computational resources. From time to time during sampling, we analyzed the mean walker positions and the mean walker position likelihood as a function of sampling step. Encouragingly, for all nuclei walkers quickly converged on a common region (within 1000 samples) and their mean parameter values stabilized soon afterward (within 10 000 samples), suggesting that walkers were sampling a reasonably optimal subspace. At this point, we considered the chain tentatively converged. As an additional test, we restarted sampling from a different (uniformly random) initial position for each nucleus and found that a similar optimal subspace was reached, again within roughly 1000 samples, indicating that our results are independent of the initial walker positions. Finally, for a "converged" chain, we calculated extracted physics quantities (e.g., neutron skins, scattering cross sections) for all walkers at several intervals to confirm that their mean values were stable. Again using 16 O and 18 O as an example, we found their mean neutron skin values varied by less than 0.001 and 0.01 fm, respectively, over several thousand sampling steps late in sampling. Out of caution (and given our expectation of very large autocorrelation times) we used only the terminal sample for each walker chain to produce the results presented here and in the companion Letter [23]. In the end, we expect that additional sampling could slightly reduce the estimated variance of each extracted quantity but have a negligible effect on the mean values. For all quantities derived from MCMC analysis, the estimated 16th, 50th, and 84th posterior percentile values are denoted as 50 84 16 . The range between the 16th and 84th percentiles corresponds to a 1σ -uncertainty range if the posteriors are assumed to be Gaussian. The median values and ranges for each parameter for each isotope system are listed in Appendix B.
Figure 12 in Appendix C shows the DOM fit of 16 O and experimental data. The experimental proton σ rxn , neutron dσ d , σ tot , and σ rxn charge density distribution, binding energy per nucleon, and p 1/2 and d 5/2 single-particle energy data are all well reproduced, suggesting that the DOM is effective for modeling nuclei as light as A = 16. Almost all experimental proton dσ d data are accurately reproduced by the DOM calculations with the exception of an overprediction of cross sections at backward angles and high energies, a regime known to be challenging from past OM analyses. In addition, the median DOM-generated rms charge radius, 2.72 fm, slightly exceeds the experimental value of 2.70 fm. Taken together with the 16,18 O relative difference results in panel (a) of Fig. 11, these overestimations indicate that the traditional OM assumption of radial proportionality with A 1/3 must be tweaked for a better description of 16 O.
To reproduce the 16 O proton σ rxn pseudodata generated from [47], a larger volume imaginary term was required above 100 MeV, which in turn reduced the spectroscopic strength for the valence π and ν p 1/2 nucleons by roughly 0.05. We also note the importance of the charge density distribution for determining the magnitude of the imaginary strength below the Fermi energy. For example, in test fits where the charge density was not included as a constraint, most of the negative imaginary strength was concentrated in the surface term between -30 < E < F MeV, and the tail of the charge density was overpredicted. With the charge density included as a constraint, the imaginary surface magnitude shrank by a factor of two and the volume term grew to compensate, pushing nucleon density deeper in energy space and increasing the binding energy closer to the experimental value.
While all data sectors contributed at least some information not fully captured by any other sector, the proton σ rxn , neutron σ tot , and charge density provided the most stringent constraints on the self-energy. The analyzing powers were the most difficult sector of experimental data to reproduce, with moderate deviations visible from 10-15 MeV for both protons and neutrons and above 100 MeV for protons [Figs. 12(b) and 12(d)]. Some of the difficulty with the analyzing powers is attributable to our neglecting of an imaginary spin-orbit term in the DOM potential used in this work, a choice made due to the unreasonable unbounded growth of the imaginary spin-orbit term as grows in the traditional • σ definition used in [21]. In a future analysis we intend to quantitatively investigate the importance of the imaginary spin-orbit term and to compare different options for its functional form.
Figure 13 in Appendix C shows the 18 O experimental data and the DOM fit. The paucity of 18 O experimental data presented a challenge for our analysis. To constrain the negative-energy domain of the potential, the only unambiguous experimental data were the neutron and proton separation energies and the overall binding energy. As with 16 O, broad agreement with experimental data was achieved for experimental proton and neutron dσ d data, the neutron σ tot , rms charge radius, binding energy per nucleon, and p 1/2 and d 5/2 single-particle energy data. The artificially scaled charge density and proton σ rxn data were also easily reproduced. Due to the deterioration of systematic trends from [47] below 70 MeV, we did not generate proton σ rxn pseudodata for lower energies, so the positive-energy surface term of the potential was largely unconstrained.
In symmetric 16 O, the proton and neutron potentials were identical except for the Coulomb interaction, so the neutron σ tot data provided information about both the proton and neutron imaginary strength at positive energies. For 18 O, this expectation of symmetric potentials was inapplicable, making proton σ rxn data essential for fixing the positive-energy imaginary strength for protons. In principle, 18 O proton and neutron differential elastic scattering cross sections about 100 MeV could jointly yield some information about the asymmetry dependence of the imaginary strength for 18 O, but no neutron elastic scattering data were available above 24 MeV. For a better characterization of this nucleus, even a single proton σ rxn datum between 10 and 50 MeV would be valuable. 40,48 Ca, 58,64 Ni, 112,124 Sn, and 208 Pb in Appendix C show 40,48 Ca, 64 Ni, 112,124 Sn, and 208 Pb experimental data and the DOM fits. The availability of single-nucleon scattering data for 40,48 Ca, 58,64 Ni, 112,124 Sn, and 208 Pb followed the same trends as that for 16,18 O: plentiful proton differential elastic scattering data, moderate coverage for neutron differential elastic cross sections and proton reaction cross sections on abundant isotopes ( 40 Ca, 58 Ni, and 208 Pb), with little to no coverage for neutron scattering or proton reaction cross section data on rare isotopes ( 48 Ca, 64 Ni, 112 Sn, 124 Sn). For 112 Sn and 124 Sn, however, even proton elastic scattering data sets were sparse and no data above 50 MeV were available, making our newly collected neutron σ tot data especially valuable in constraining the potential. For 40 Ca and 208 Pb, experimental proton reaction cross section data were available up to 200 MeV; for the other isotopes, proton reaction cross section pseudodata (discussed in the 16,18 O subsections) were used as a constraint. As for 18 O, no charge density parametrization was available for 112 Sn in [49], so we rescaled the available 124 Sn distribution to reproduce the 112 Sn charge radius.
Generally, all sectors of experimental data were well reproduced; exceptions include the large-angle (above 120 • ) proton elastic scattering data for 40 Ca and 208 Pb, where data sets were available up to 200 MeV, and the single-particle energies for neutron open shells in 112,124 Sn (see Figs. 18 and 19), where several levels are partially filled and clustered near the Fermi surface. Achieving more accurate single-particle energies while preserving particle number accuracy may require a more sophisticated treatment of pairing. Our new neutron σ tot data were well reproduced across the board, typically within 2% of the experimental value, by the DOM fits, suggesting that our Lane-like parametrization of the potential's asymmetry dependence [Eqs. (A15)-(A18)] is a promising starting point for extrapolation away from stability. We note that because 208 Pb was fitted on its own without an isotopic partner, initial fits showed that the asymmetry-dependence of the Hartree-Fock radius term was too poorly constrained to yield reliable neutron skin results; in the final treatment, this term was disabled for 208 Pb.
Table III shows DOM-calculated SFs for valence proton and neutron levels for all nine systems. Significant depletion from the mean-field expectation appears even in the light systems 16,18 O. In the present study, the extracted proton SFs show only a very weak dependence on neutron-richness within each isotopic pair, in keeping with the weak depen- dence extracted in (e, e p) and transfer reaction studies and at odds with knockout-reaction analyses that recover a strong asymmetry dependence [45,56]. The recent DOM analyses of [43,57] identified proton reaction cross sections above roughly 100 MeV as important for their successful reproduction of 40,48 Ca (e, e p) cross sections without arbitrary SF rescaling. Compared to the present work, these analyses found a much larger reduction of valence proton SFs in 48 Ca with respect to 40 Ca, indicative of an SF asymmetry dependence somewhere between the weak dependence deduced from transfer reactions and the very strong dependence from knockout reactions.
To understand the differences between these analyses, we conducted several diagnostic runs with artificially scaled Carlson pseudodata in 48 Ca. These diagnostic runs confirmed that fitting to appropriate high-energy proton reaction cross sections leads to larger 48 Ca proton imaginary strength both far above and far below the Fermi energy, an effect already seen in previous DOM work. However, the growth we observed in the imaginary potential was more modest compared to previous treatments, potentially explaining the weaker asymmetry-dependent SF reduction. We also note that in the present work, the high-energy neutron total cross sections and proton reaction cross sections appeared to have little impact on other extracted quantities such as neutron skins, as had been previously hypothesized for the neutron skin of 48 Ca [1]. We conclude that the different methodological choices, especially the focus of this work on simultaneous fitting of isotope pairs, is responsible for the differences in these asymmetry-dependent quantities. To further clarify the situation, the potentials of the present work should be used to generate (e, e p) cross sections that can be compared to the previous findings of [57].
Surprisingly, despite the extensive proton and neutron elastic scattering data for 16 O, 40 Ca, and 208 Pb, the extracted spectroscopic factor distributions and parameter uncertainties for these isotopes are just as wide as for those systems with barely any available elastic scattering data, such as 64 Ni. We tentatively conclude that the elastic scattering data we used are very weak constraints on the all-important imaginary terms of the optical potential, at least for the stable, spherical systems discussed here. Unfortunately, this suggests that elastic scattering measurements in inverse kinematics on radioactive beams are of diminishing utility for extrapolating optical potentials away from β stability. A program of proton reaction cross section and neutron total cross section measurements on radioactive targets could be useful for understanding the potential's near-Fermi-level asymmetry dependence but is experimentally daunting. Instead, a two-pronged approach may be required. On the experimental side, proton reaction and neutron total cross section measurements on stable isotopic chains can help identify which asymmetry-dependence forms are justifiable for increasingly asymmetric systems. On the theoretical side, sensitivity studies are needed to clarify how bound-state data on highly asymmetric systems connect to scattering cross sections.
Lastly, a few systematics in optical potential parameter values are worth mention. For most of the parameters, there was minimal variation with nuclear size or asymmetry, suggesting that a global DOM treatment using the functional forms we have selected is achievable. The radial term for the real central potential (r 1 ) and for the positive-energy imaginary volume and surface (r + 4 , r + 5 ) are nearly constant among 40,48 Ca, 58,64 Ni, 112,124 Sn, and 208 Pb, but the values for 16,18 O show moderate deviations, another indication that the geometric form of the potential is insufficient for light systems. As a consequence of the limited negative-energy data available for fitting, the negative-energy geometric terms (r - 4 , r - 5 , a - 4 , a - 5 ) show large variation. The nonlocalities for the negative imaginary components are systematically larger than those for the positive imaginary components. This suggests that while traditional OMs have been able to successfully reproduce positive-energy scattering data with strictly local potentials, description of hole properties requires true nonlocal character in the negative-energy potential. In practice, we found it impossible to simultaneously reproduce charge density distributions, binding energies, and scattering data unless the central potential and at least the volume imaginary terms were equipped with a nonlocality. In the end, for simplicity and generality, each element of the potential (except Coulomb) was treated nonlocally, but it is unclear which particular data are most important for constraining these several nonlocalities. As one moves further from stability to systems with even less (or no) scattering data available, the risk of overfitting will loom until this issue is resolved.
In preliminary fits, the imaginary volume magnitude (A - 4 ) component of the potential was shown to be strongly sensitive to the inclusion of the binding energy as a constraint during fitting. We expect the asymmetry dependence of this term (A - vol,asym ) to impact DOM-based predictions of the Ca, Ni, and Sn neutron drip lines (as in [28]), though in this work, this dependence was very poorly constrained due to the absence of experimental asymmetry-dependent data probing the most deeply bound nucleons. Because they encode information about how protons and neutrons share energy throughout the nucleus, experimental neutron-skin thicknesses could provide this kind of valuable information.
For the Ca, Ni, Sn, and Pb fits, the median positive-energy surface imaginary magnitude (A + sur,asym ) is positive, indicating enhancement in proton surface imaginary strength with increasing neutron richness and a corresponding decrease for neutron surface imaginary strength. Of course, the nuclei under study in the present work are stable; the trend for nuclei with large asymmetries, relevant for the r-process neutroncapture rate, is unknown.
By adopting a digitizer-driven approach, we measured σ tot on the important closed-shell nuclides 16,18 O, 58,64 Ni, and 112,124 Sn across more than two orders of magnitude in energy (3-450 MeV). Except at the highest energies, our results on natural targets are in good agreement with previous analogmediated measurements that required an order of magnitude more target material.
Using these new data and a suite of scattering and boundstate literature data on 16,18 O, 58,64 Ni, and 112,124 Sn, we extracted DOM potentials capable of reproducing a diverse range of scattering and structural data for both neutrons and protons, validating the use of the DOM away from doubly closed shells from A = 16 to A = 208, though with indications that the traditional A 1/3 radial dependence may require modification for light systems. These analyses further indicate that simultaneous fits of isotopically resolved neutron σ tot , proton σ rxn , and charge density distribution data on isotopic partners provide a more stringent constraint on the asymmetry dependence of both real and imaginary components. Experimental data are plotted as black points and pseudodata generated from [47] are plotted as gray open circles. Panels (f) show the neutron σ tot and σ rxn . The charge distributions of panels (g) are derived from the compilation of [49] (see comments in the DOM Analysis section), and are displayed with an arbitrary 1% uncertainty band in black. In panels (h), singleparticle energies nlj are shown as horizontal lines. In the "calc" column, DOM-calculated single-particle energies are plotted; the height of each rectangle spans the 1σ calculated uncertainty for that level. Panels (i) show DOM-calculated charge radii; the experimental charge radius is displayed using dark gray and light gray bands representing 1σ and 2σ uncertainties, respectively. Panels (j) show the DOM-calculated binding energy per nucleon; the experimental value is shown with a thin gray band.