Reconstruction of Bremsstrahlung γ -rays Spectrum in Heavy Ion Reactions with Richardson-Lucy Algorithm

,


I. INTRODUCTION
The strong interaction between nucleon pairs at short distance with high relative momentum in nucleus is called short-range correlation (SRC).Because of the SRC effect, the high momentum tail (HMT) in the momentum distribution of the nucleons in nuclei above the Fermi surface exists [1] and can be measured experimentally.So far much of our knowledge on the SRC, which gives suggestively rise to the EMC effect [2], is obtained from highenergy electron-nucleus scatterings [3][4][5][6] and from proton induced reactions [7,8].Very recently, it was pointed out that the existence of the HMT nucleons may harden the energy spectrum of the energetic photons emitted in np Bremsstrahlung process in heavy ion reactions at intermediate energies [9].
Because of the merit that the energetic photons experience rare final interactions with the surrounding nuclear medium once they are produced, experimental measurement of the Bremsstrahlung photons follows up.In our previously work [10], the energy spectrum up to 80 MeV of the Bremsstrahlung γ-rays has been measured and analyzed using the isospin-and momentum-dependent Boltzmann-Uehling-Uhlenbeck (IBUU) simulations.Despite of the insufficient statistics and the high energy border of E γ < 80 MeV, the experimental results suggest the existence of HMT at a confidential level of about 90% and showcase the feasibility to study SRC by the Bremsstrahlung photons.To improve the confidential level, further experimental efforts are required.
However, there is a hindrance.When comparing the experimental results to the theoretical calculation incorporating the HMT, the detector filtering effect must be taken into account.This is not trivial because the high energy γ-rays undergo complicated transport process in the detection system.Usually, it is convenient to solve the forward problem, i.e., to convert the theoretic prediction to the expected spectrum by the simulation of the detector response.From a practical point of view, it then hinders further comparative studies unless the accompanying response matrix of the detection system is always provided.Thus, it is favorable to solve the inverse problem, i.e., to reconstruct the original spectrum from the measured one.
The Richardson-Lucy (RL) algorithm [11,12] based on Bayesian theorem has been widely applied in optical imaging deblurring problems [13], mainly to repair the blurred images with known point distributions.In recent years, it has seen increasing applications in some nuclear physics problems.For instance, the velocity distribution of emitted particles measured experimentally can be used to infer the unknown velocity distribution of the center of mass system, thereby to determine the original momentum distribution of HIC products [14].The RL algorithm has also been applied in the analysis of particle correlation functions, from which one can infer the distribution of the source functions [15].In the field of energy spectrum measurement and analysis, the RL algorithm has been proven practical in obtaining the shell structure of unbound particle systems from the measured decay energy spectrum [16].In the data analysis of HIC experiments, statistical fluctuations can affect the stability of results, and RL algorithm is still an effective tool for optimizing the stability of solutions [17].In addition, during the iterative process on which the RL algorithm arXiv:2405.06711v2[nucl-th] 15 Aug 2024 is running, the output of each iteration maintains a positive value.For the problem of solving γ spectra, the RL algorithm avoids directly inverting the detector response matrix and efficiently solves the original γ spectra.
Therefore, we are motivated to introduce the RL algorithm in the reconstruction of the original γ spectrum, which can be compared directly to the model prediction.The paper is arranged as following.In section II, we describe the RL algorithm for solving the inverse problem in γ measurement.In section III we test and verify the accuracy of the RL algorithm.In section IV, the RL method is applied to reconstruct the original γ spectrum from the measured one in 25 MeV/u 86 Kr + 124 Sn [10], and the comparison between the reconstructed original γ energy spectrum and the simulation of IBUU model is discussed.Section V is the summary.

A. Model Description
In optical deblurring problems, a certain real property µ of a photon is measured as ν in experiments.The distribution function F(µ) of the real property µ and the distribution function f (ν) of the measured property ν is connected by the following equation, Here, P (ν|µ) is the conditional probability of measuring photons with real properties µ as ν.
Defining Q(µ|ν) to represent the probability density of the real value being µ if ν is measured, which is the complementary probability density of P , one can connect P (ν|µ) and Q(µ|ν) by the Bayesian law.Assuming that the measured value is within the range of dν, the probability within the range of dµ satisfies the following equation, Equivalently, Q(µ|ν) can be expressed as, and F(µ) can be expressed as where W (ν) is the weight of the relative importance of specific data point at ν to infer F. The weight of each measurement in the calculation can be manually specified to reduce the impact of unreliable data on the calculation.Meanwhile, by varying the W one can test the stability of the final solution obtained through calculation.
Use the RL method to iteratively solve the above two equations, one writes and Here the superscript r is the iteration index.It is seen that starting with a given initial distribution F (0) (µ), one can obtain the final distribution after certain iteration times achieving convergence.The implicit dependency of P on F can be handled in the iteration of the equation, and the non-negative distribution remains unchanged during the iteration process.By combining the above two equations [14], we can obtain Here, A (r) (µ) is called the amplification factor, and f (r) (ν) is the distribution function of ν obtained after the r th iteration, that is B. Application to the γ Energy Spectrum Reconstruction In the experiment of HIC, γ with a certain energy µ is measured to have energy ν through the response of the detector.Similar to optical deblurring problems, the RL algorithm can also be used to iteratively solve the original energy spectrum of γ.
The relationship between the distribution function E(µ) of the discret real energy spectrum and the distribution function e(ν) of the measured energy spectrum is written in matrix form, Here, D ij is the detector response matrix, which can also be understood as the conditional probability that the γ measurement with real energy µ is ν.The values of matrix elements D ij can be obtained generally by Geant4 simulation if the experimental setup and data analysis scheme are provided.
According to the RL formulation, one can reconstruct the original energy spectrum using the following iteration, where r in the above equation represents the iteration index, and the amplification factor A (r) i after each iteration can be written as where e (r) is the value of the predicted measurement spectrum obtained in the r th iteration, i.e. e (r) T ji is the weighted energy spectrum transformation matrix, which has been automatically normalized as, In HIC experiments, the manifestation of measuring energy spectrum is the counts of detector signals in different energy bins.The count distribution in each bin satisfies the Poisson distribution.Accordingly, the weight W of each measurement energy point used in iteration can be set to √ e.

III. ALGORITHM TEST
In order to verify the effectiveness of the RL algorithm in solving the γ original energy spectrum, the IBUU model was applied to calculate the energy spectrum of the bremsstrahlung γ emission from HICs.The spectrum is then filtered by the detector response obtained by Geant4 [18] simulations.Once the energy spectrum from the detector output is generated, one can reconstruct inversely the input energy spectrum using the RL algorithm and compare to the original one.The consistency is a proof of the effectiveness of the RL algorithm.
For the test we take the 25 MeV/u 86 Kr + 124 Sn reactions performed on the compact spectrometer for heavy ion experiment (CSHINE) [10,19].A 4 × 4 CsI(Tl) hodoscope was mounted on CSHINE to measure the high energy γ rays [19,20].The central 2 × 2 units in the hodoscope were used to record the event center where the highest energy deposit was found.The total γ energy was then collected from all neighbouring units fired with a given time window.For the details of the experiment, one can refer to [19].With this detector setup, one can simulate the detector response with Geant4 packages.As an example, Fig. 1 shows the response of the detector caused by E γ = 60 MeV single energy photons.The red histogram presents the probability P i (E j ) of a single energy E j photon (here E j = 60 MeV indicated by the blue line) producing a response E i in the hodoscope, equivalently to say, D ij = P i (E j ).
The BUU transport model [21][22][23] effectively describes the dynamics of nucleus-nucleus collisions, with its main equation given by where f (r, p, t) is the probability of finding a particle at time t with momentum p at position r, U represents the mean-field potential, and the evolution of f (r, p, t) by elastic and inelastic two-body collisions is governed by the collision term I coll .The original γ spectrum is calculated by the IBUU transport model, in which the SRC and HMT effects have been incorporated [9,24].The simulations utilize the initial single-nucleon momentum distribution n(k) that includes HMT.The behavior of HMT is remarkably consistent from the deuteron to infinite nuclear matter.The n(k) for the HMT induced by SRC [25] is given by where k represents the single nucleon momentum, k F is Fermi momentum and λ is the high-momentum cutoff.Due to the low probability of producing Bremsstrahlung photons in the reaction, the impact of bremsstrahlung on nucleon kinematics is also negligible, and hence, perturbation methods can be used to calculate the probability of photon production in each np collision.Based on the single boson exchange (OBE) method, a well fitted expression for the probability of elementary double differential photon production can be applied in IBUU simulations [26], Here, E γ represents the energy of produced photons, and E max represents the total available energy in the center of mass frame.The coefficient α is α = 0.7319−0.5898βi , where β i represents the nucleon velocity.The total photon production probability per event is the sum of the probabilities of all np collisions producing photons throughout the entire process.From the total probability, one can obtain the original γ spectrum with given collision times N evt .
In our test, the stiffness parameter x of the symmetry energy is taken as x = −1.The triangular distribution of b ≤ 5 fm is adopted to describe the collision centrality.Three ratios of the high momentum tail (R HMT ) are adopted in the calculation, including R HMT = 0% as a free Fermi gas (FFG), R HMT = 15%, and 30%, respectively [9,24].Taking the theoretic spectrum with R HMT = 0% (not shown) as the initial spectrum in the first iteration E (0) , one can solve reversely the original spectrum using the RL algorithm with the formulae (9) to (13), as shown by the solid symbols.Indeed the choice of E (0) does not affect the final output.For the high statistics test in panel (a), it is shown the reconstructed spectrum (solid symbols) is in agreement with input original one (red dashed curves).To estimate the uncertainty in the gamma spectrum obtained through the RL algorithm, we used multinomial sampling on the central energy spectrum derived from the experimental data.The sampled spectra were individually processed using the RL iterative algorithm, yielding a set of reconstructed spectra.For each bin in these spectra, the standard deviation across the set was calculated and used as the uncertainty in the final extracted spectrum.For the low statistics test in panel (b), the reconstructed data points situates also on top of the original curve, while the statistical uncertainty are significant.It is demonstrated that the RL algorithm is very effective in solving the γ original energy spectrum.Clearly, when the counts of measurements increase, the smoothness of inverse solution will become better.It is then clearly favorable to increase the statistics in order to reduce the error of the inverse solution results and to improve the confidence level of the conclusion in real applications.

IV. APPLICATION TO EXPERIMENTAL RESULTS ANALYSIS A. Experiment Data Solution
The RL algorithms discussed above are applied to reconstruct the original energy spectrum of the Bremsstrahlung γ emitted in 25 MeV/u 86 Kr + 124 Sn performed at CSHINE [10,19].Iteration times is an important factor to ensure the convergence, meanwhile to keep the computational efficiency.The terminating point of the iteration is identified by monitoring the standard deviation δ exp of the experimental spectrum and the predicted spectrum in detector response after each iteration.The quantity δ exp is defined FIG.4: (Color Online) The the output of the RL algorithm for the original spectrum (upper row) and the predicted spectrum in measurement (lower row) after 0 th (a,e), 5 th (b,f), 10 th (c,g) and 100 th (d,h) iteration, respectively.See text for details.
as the following, here, e i is the input measurement value of the corresponding i th energy point and e P i is the predicted value obtained in the current iterations.N p is the number of bins.If √ e i in the denominator equals zero, it is replaced by e P i to avoid a segment fault.The distribution of δ exp as function of iteration time is depicted in Fig. 3.It is shown that after a certain number of iterations, the value of δ exp gradually converges and remains constant.The final nonzero δ exp value mainly comes from the contribution of the null bins where the experimental statistics is zero.
In order to gain a more detailed view on the iteration process, we present in Fig. 4 the output of the RL algorithm for the original spectrum (upper row) and the predicted spectrum in measurement (lower row) after the 0 th , 5 th , 10 th and 100 th iteration, respectively.As the initial input, we adapt the theoretic spectrum with R HMT = 0% as the real original spectrum, as shown by the red dashed curve in panel (a).After the detector response filter, one expects the spectrum of the detector output, shown by the blue histogram in panel (e), it exhibits large disagreement with the experimentally measured spectrum shown by the red histogram.After 5 iterations, the low energy part is fast converged because in the low energy part the full energy can easily be collected by the detector and only the diagonal matrix element is nonzero.With the interation times increasing, the predicted spectrum of the detector output approaches gradually to the real one, as shown in pan-els (f-g).After 100 times of the iteration, as shown in (h) the predicted spectrum of detector output is in good agreement with the experimental spectrum.Meanwhile, as shown in panel (d), the original γ spectrum is situating in the vicinity of the IBUU calculation, except for some jumping points at high energy side, originated from the large fluctuations on the experimental spectrum.Thus, in our analysis, we set conservatively the iteration number to 500 in RL algorithm implementation.The inverse solution of high-energy γ spectrum in 25 MeV/u 86 Kr + 124 Sn is presented in Fig. 5 in comparison to IBUU predictions with R HMT = 0%, 15% and 30%, respectively.The central values of scatter points are reconstructed by RL algorithm from the input measured spectrum and the error bars are derived from the standard deviations when one samples the content in the corresponding bin using multinomial distribution to ensure that the total count of our pseudo-gamma spectrum matches that of the experimental statistics.Appendix A and B present the reconstructed original spectrum (in CM) using RL method and the experimental spectrum of detector output (in laboratory), respectively.With these results, one can readily analyze the most probable HMT ratio of nucleons in nuclei.

B. HMT Analysis
Next, we are to determine the favored ratio R HMT from the RL solution.Let's define χ 2 c to evaluate the closeness (overall deviation) of the theoretic spectrum to the inversely solved original γ spectrum in the experiment, where the two terms in the square bracket denote the possibility in the i th bin of energy for experiment and theoretic prediction, respectively.The super and subscript c represents the percentage of HMT in IBUU simulations, σ i in the denominator is the standard deviation obtained from the sampling simulation and n DOF refers to the number of degree of freedom (DOF), which is approximately taken as the number of data points in the corresponding energy range.In addition, the scales of different theoretical curves shall be adjusted to minimize the χ 2 c (the χ 2 c s mentioned in the following process are all minimized).This quantity is used to evaluate the closeness of the inversely solved energy spectrum to various spectra of IBUU model simulations.Energies below E γ < 30 MeV are largely influenced by collective resonance and statistical emissions, which are excluded from this analysis.At energies higher than 80 MeV, the impact of cosmic rays becomes significant, limiting the applicability of the IBUU model.Table .I lists the χ 2 c values using different R HMT in the energy range of interest (ROI) from 35 to 80 MeV.From the table, the minimum χ 2 c value is found at R HMT = 25%.It indicates that the HMT exists in nuclei, qualitatively consistent with our previous analysis.On the other hand, the most favored ratio of R HMT = 25% is somewhat high than the most probable value derived in [10], see next.
It is also necessary to conduct the hypothesis test in order to determine the probability to exclude the R HMT = 0%, i.e., to exclude the assumption of null high momentum tail of nucleons in nuclei.To achieve this goal, the theoretic curve of R HMT = 0% is multiplied by a factor ηN evt to have a hypothetical spectrum with certain counts in each bin without fluctuation.Here N evt = 51269108 is the number of events recorded in experiment, and η = 1.22 is a normalization factor to make the summing counts above 20 MeV equal to that of experimental data.With this hypothetical spectrum, one can generate a psudo-experimental spectrum, taking the fluctuations into account by sampling the content in each bin of E γ with multinomial distribution.Starting with the psudo-experimental spectra, one can repeatedly use the detector response filter and solve the RL algorithm to get the original γ spectra with the zero HMT assumptions.By computing the difference of χ 2 0 − χ 2 25 with ( 16) for the psudo-experimental spectra in the same way as done in the real experiment, one can infer the possibility to conclude a non-zero HMT (R HMT = 25% here) where the input condition is R HMT = 0%., the area on the left (right) side of the arrow adds up to 80% (20%), suggesting that there is a probability of 20% to observe a non-zero HMT behavior from a "gedanken experiment" without HMT on the nucleon momentum distribution.Equivalently, it is stated that the experimental results exclude the zero HMT hypothesis at 80% confidential level.
Finally, we check the look-elsewhere effect.Tab.II lists the values of χ 2 c at different R HMT and in different en- a Beyond 30%, the χ 2 c in these two columns decreases with R HMT due to the existence of the last data points, but the difference of the hardness of the spectrum by varying R HMT becomes less pronounced and smaller than the uncertainty.Thus the calculations is terminated at 30%.ergy ROI, together with the favored R HMT and the corresponding confidential level to exclude the zero HMT assumption.It is shown that the favored ratio of R HMT in different analysis ranges is changed.With the only exception in ROI from 30 to 80 MeV, all the other analyses favor a non zero HMT, with the best value depends on the ROI.The confidential level varies from 42% to 96%, and the average confidential level of excluding the zero HMT hypothesis is about 74%, suggesting the existence of SRC effect in nuclei.
The results from this new analysis is in qualitative agreement with our previous conclusion by analyzing directly the measured spectrum without introducing the RL method to solve the inverse problem [10].Considering the small difference of χ 2 c here as well as the small difference of the likelihood by varying R HMT in [10], it is understandable that the favored value of R HMT exhibits some variation.Meanwhile, due to the low experimental statistics, the confidential level to exclude zero R HMT is yet not high, i.e., in the vicinity of 2σ at most.In order to increase the resolution power of HMT ratio from the Bremsstrahlung γ emission, it is necessary to increase the statistics and to extend the high energy border because the discrepancy caused by different R HMT is increasingly pronounced if one goes to higher γ energy.Further efforts are accordingly required in the future experiments [27].

V. CONCLUSIONS
For the first time, the Richardson-Lucy algorithm has been applied to solve inversely the original energy spectrum of Bremsstrahlung γ-rays in heavy ion reactions.The γ spectrum in 25 MeV/u 86 Kr+ 124 Sn measured at CSHINE are reconstructed and compared to the isospin-and momentum-dependent Boltzmann-Uehling-Uhlenbeck transport model simulations.It is inferred that a ratio of high momentum tail R HMT = 25% describes the experimental data points with the least χ 2 c in the ROI from 35 to 80 MeV and the confidential level to exclude the zero R HMT is averagely 74%, in qualitative agreement with our previous conclusion without introducing the RL inverse solution process.
Although the application of the RL method has its advantage that the detector response has already been folded in solving the inverse problem to reconstruct the original γ spectrum, the resolution power to quantify the HMT ratio is indeed determined by the statistics and the upper limit of the γ energy covered by the experiment.So, to explore further the SRC effect in atomic nuclei with enhanced accuracy from the Bremsstrahlung γ emission in heavy ion reactions, further experimental efforts are called for.
Photon Energy, E

FIG. 2 :
FIG. 2: (Color Online) Comparison of the input and output spectra (in CM reference frame) of RL algorithm for a high-statistics (a) and low-statistics (b) test.The dashed curves are the input from IBUU simulations with R HMT = 15%.The solid squares with the grey bands are the output of the RL algorithm.Shown in the inset are the expected spectra in laboratory with the detector response filtering.

FIG. 3 :
FIG.3:(Color Online) The variation of the standard deviation δ exp between the experimental spectrum and the predicted spectrum in measurement as a function of iteration time in solving the high-energy γ spectrum reversely.

TABLE I :
χ 2 c s of direct solution high-energy γ spectrum compared with different HMT percentages in the energy ROI from 35 to 80 MeV.

TABLE II :
The χ 2 c s between experimental solution results and theoretical models at different HMT ratios in different energy ROI, the HMT ratio preferences, the confidence level of supporting the favored R HMT are also listed.