PARTICLE CORRELATIONS IN HEAVY-ION COLLISIONS By Pierre Nzabahimana A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics—Doctor of Philosophy 2023 ABSTRACT In this thesis, our focus is on studying particle correlations in heavy-ion collisions to gain insights into nuclear systems in the final state reaction. Understanding these correlations is crucial for accessing the geometry, phase-space features, and time development of the collision’s final stages. We present two approaches to extract the relative distribution of particles at the last moment of the collision: the Gaussian parametrization source (GPS) and the deblurring method. In the GPS approach, we assume a Gaussian form source to approximate two-particle correlation functions using the Koonin-Pratt (KP) convolutional formula. This formula convolves the relative emission source with the squared two-particle relative wave function. We apply the approach to study the correlations of low-velocity alphas. We start by constructing the scattering wave function for the alpha-alpha pair by solving the Schrödinger equation, incorporating a potential tailored to match the measured phase shifts of the system. With this wave function, we interpret available data on alpha-alpha correlations in terms of emitting sources. In the deblurring approach, we propose using the Richardson-Lucy (RL) optical deblurring algorithm to deduce a source from the correlation function. The RL algorithm, derived from prob- abilistic Bayesian considerations, requires the optical object, image distributions, and convolution kernel to be positive definite. Fortunately, these conditions are satisfied by the corresponding quantities of interest within the KP formula. We demonstrate the success of the RL algorithm in restoring emitting sources from measured d–𝛼 correlations. Furthermore, we extend the deblurring approach to another field of nuclear physics by utilizing the RL algorithm on experimental nuclear physics data, leveraging only the observed energy spectrum and the detector’s response matrix (also known as the transfer matrix). This technique provides access to information regarding the shell structure of particle-unbound systems through the measured decay energy spectrum, which is not readily attainable through traditional approaches like chi-square fitting. In pursuit of the same objective, we develop a machine learning model that employs a deep neural network (DNN) classifier to identify resonance states from the measured decay energy spectrum. We evaluate the performance of both methods using simulated data and experimental measurements. Subsequently, we apply both algorithms to analyze the decay energy spectrum of 26O, as measured via invariant mass spectroscopy. Both the deblurring and DNN approaches indicate the presence of three peaks in the raw decay energy spectrum of 26O. Finally, we employ the transport model in this thesis to analyze two-proton (p-p) correlations in heavy-ion collisions at low incident energies per nucleon (E/A). Specifically, we utilize the Boltzmann-Uehling-Uhlenbeck (BUU) transport model to simulate the p–p source. Subsequently, we employ the source and the p-p kernel within the KP formula to calculate the correlations. Through a comparison between the correlations obtained from the BUU simulation and the RL algorithm, we gain a better understanding of the influence of fast and slow emissions on the measured correlations. We compute the angle-averaged and quadrupole components of the p-p source for the Ar + Sc and Xe + Au reactions at 𝐸/𝐴 =80 MeV. These sources are computed considering both momentum-independent and momentum-dependent nuclear equation of states (EOS), enabling us to observe the effect of the momentum-dependent EOS on the quadrupole component source. Copyright by PIERRE NZABAHIMANA 2023 ACKNOWLEDGEMENTS I would first like to thank my advisor, Pawel Danielewicz, not only for his academic guidance but also for his support, patience, and liberty. He is the kind of liberal and laid-back person who would not complain even when my progress is too little. I also appreciate that he has provided me with great support and opportunities to travel and attend conferences, and that he would go the extra mile to hunt for job opportunities and connect me to people. I cannot forget the battle he, Scott Pratt, and Kirsten Tollefson fought for me to get into Michigan State University. I am very grateful to them. I am also grateful to others of my committee members, Betty Tsang, Carlo Piermarocchi, Filomena Nunes, and Scott Pratt, for their guidance, whether through critical questions or heartfelt advice and of course, for putting up with my talks every year. I would also like to thank my collaborators: Thomas Redpath, Thomas Baumann, Pablo Giuliani, Paul Guèye, and Kin Tam. Working with them was fruitful. I want to extend my thanks to the supporting staff at the Department of Physics and Astronomy and NSCL/FRIB, especially Kim Crosslan, Elizabeth Deliyski, and Gillian Olsen. I would like to thank Rachel Younger and Mornetka Guèye for checking on me. I also want to give special gratitude to Paul Guèye for always being there for me. I want to thank the professors and instructors I encountered during my PhD coursework, such as Alex Brown, Michael Murillo, Filomena Nunes, Luke Robert, Wade Fisher, Mohammad Maghrebi, Danny Caballero, Carlo Piermarocchi, Vladimir Zelevinsky, and Yamazaki Yoshishige. The spiritual support and companionship of my peers from the Chi-Alpha Group and MSU chapter of NSBP mean a great deal to me. I take this opportunity to send my gratitude and best wishes to my cohort, including office mates, friends, and others, such as Joshua Wylie, Jordan Purcell, Sylvester Agbemava, Jacob Davison, Mengzhi Chen, Xian-Gai Deng, Linda H’lophe, Tong Li, Hossein Mahzoon, Avik Sarkar, Boyao Zhu, Kuan Zhu, Felix Ndayisabye, Jean Pierre Twagirayezu, Tracy Eduards, Timilehin Ogunbeku, Yani Udiani, Hao Lin, Letrell Harris,..., and whoever is reading this! v The support from my family, especially my beloved wife Dative Mariza, and my children Ella E.I. Gahozo and Eloi S. Cyizere, has been my everyday motivation, and I am grateful for them. To all the friends I met during this journey, their company means a lot to me. The work in this thesis is supported by the U.S. Department of Energy Office of Science under Grant No. DE-SC0019209, the U.S. National Science Foundation Grant No. PHY-2012040. vi TABLE OF CONTENTS CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Overview on particle femtoscopy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Transport Model 1.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Thesis organization . . 1 3 6 7 CHAPTER 2 Introduction . LOW-MOMENTUM CORRELATIONS OF 𝛼 PARTICLES IN HEAVY- ION COLLISIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 9 2.1 2.2 Two-Particle Correlation Functions . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 Results: 𝛼 − 𝛼 and 𝑑 – 𝛼 CF . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.4 Conclusion . . . . . . . . . CHAPTER 3 DEBLURRING APPROACH . . . . . . . . . . . . . . . . . . . . . . . 31 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1 . Introduction . 3.2 Deblurring method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.3 Deblurring to restore source function from two-particle correlation . . . . . . . 38 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.4 Summary . . . . . . . . . . . . CHAPTER 4 . . . . DEBLURRING EXPERIMENTAL DECAY ENERGY SPECTRA: THE 26O CASE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Introduction . . 46 4.1 Interplay of the Experiment and the Analysis Methods . . . . . . . . . . . . . . 47 4.2 4.3 Richardson-Lucy Deblurring Algorithm . . . . . . . . . . . . . . . . . . . . . 52 4.4 Tests of deblurring algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.5 Deblurring 26O decay energy spectrum . . . . . . . . . . . . . . . . . . . . . 59 . 62 4.6 Conclusions and Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 5 Introduction . MACHINE LEARNING MODELS . . . . . . . . . . . . . . . . . . . . 64 . 64 . 65 5.1 5.2 Deep neural network Model 5.3 DNN architecture to discover resonance states from the measured decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . energy spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.4 DNN Classification model for resonance peaks detection . . . . . . . . . . . . 75 5.5 Comparison of the DNN method with deblurring method . . . . . . . . . . . . 77 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.6 Summary . . . . . . CHAPTER 6 TRANSPORT MODEL SIMULATION FOR TWO-PARTICLE SOURCE FUNCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 6.1 BUU model 6.2 Two-proton correlation function . . . . . . . . . . . . . . . . . . . . . . . . . 85 6.3 Quadrupole moment of Source function . . . . . . . . . . . . . . . . . . . . . 97 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.4 Summary and outlook . . . CHAPTER 7 CONCLUSION AND OUTLOOK . . . . . . . . . . . . . . . . . . . . 101 vii 7.1 Conclusion . . . 7.2 Outlook . . . . . . . BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 . 102 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 APPENDIX A DERIVATION OF SCATTERING PHASE SHIFT . . . . . . . . . . . . 114 APPENDIX B THERMODYNAMIC MODEL FOR PARTICLE EMISSION . . . . . 117 APPENDIX C QUADRUPOLE COMPONENTS OF SOURCE FUNCTION FROM BUU SIMULATION . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 viii CHAPTER 1 INTRODUCTION In this thesis, we will focus on heavy-ion collisions at low energy. These are nuclear processes in which a heavy projectile, typically with a mass number 𝐴 > 4, impinges on another heavy target at an incident energy less than 100 MeV per nucleon. As illustrated in Fig. 1.1, initially, the two nuclei come into contact and compress against each other. Due to nucleon-nucleon collisions, both the temperature and entropy in the reaction will rise, leading to the formation of a finite system constituted of hot and dense nuclear matter. At this stage, some particles will leave the system; this is known as fast emission. Over time, the system undergoes expansion, thereby reducing temperature and density. The density continues to fall until it reaches the freeze-out moment, i.e., the nucleons in the system are no longer interacting. Most of the particles detected are those emitted at the freeze-out moment [1, 2]. Figure 1.1 The figure shows a schematic illustration of particle emission during a heavy-ion collision. Some particles leave the system early in the process (i.e., fast emission), while others leave later as a result of the decay of intermediate unstable nuclei (i.e., slow emission). This figure is a modified version of the figure in [3]. It is also possible that the emitted fragments (i.e., unstable intermediate nuclei) decay into daughters at a later time during the collision; this is known as secondary decays (see Fig. 1.1). These processes provide insights into the properties of the compound nucleus, the underlying nuclear interactions, and the correlations between the final-state particles. The aim of this thesis is to investigate the correlations between final-state particles to gain insights into the nuclear reaction system. These correlations are measured in heavy-ion collision 1 experiments and are defined as the ratio of coincidence yield (symbolized by 𝐴) to the product of single yield (symbolized by 𝐵): 𝐴(P1, P2) = [𝑅(q) + 1]𝐵(P1, P2), (1.1) In this equation, 𝑃1 and 𝑃2 represent the momenta of the first and second particles, respectively. When 𝐴 = 𝐵 (i.e., 𝑅(q) = 0, particularly notable at high values of 𝑞), the particles are considered uncorrelated; otherwise, they are correlated. During the collision, as the particles exit the high-density region of the system, they can interact and establish correlations (or anti-correlations) between their momenta, provided their separation distance is sufficiently small [1, 4]. When the particles are emitted due to the decay of unstable intermediate nuclei, the correlations are predominantly characterized by resonance peaks resulting from nuclear interactions between the particles. All of these aspects will be considered in this thesis. Additionally, the Coulomb interaction effect manifests in the correlation as a valley around zero relative momentum [5]. Consequently, our particular focus lies on charged particles such as 𝛼 − 𝛼 and d−𝛼, as they are likely to be detected in the heavy-ion collision experiments [6, 7, 8, 9]. Moreover, these correlations can be utilized to gain insights into the geometry of the emitting system at a later stage of the collision by calculating the relative distribution function of particle pairs in the emitting system, referred to as the source function [4, 10, 11, 12, 5, 13]. To extract a reliable source function, we employ trusted methodologies that utilize the measured correlations and the relative kernel of the pairs to reconstruct the source function. This thesis discusses the different methods we have developed to obtain a desirable source. Furthermore, in this dissertation, we employ a transport model simulation, introduced in Sec- tion 1.2, to calculate proton-proton correlations. Two-proton correlations serve as important tools for testing transport theory, as they are sensitive to the space-time properties of nuclear reactions [14] (Refs within). Additionally, we examine the effect of the momentum-(in)dependent nuclear equation of state (EOS); the EOS represents the behavior of energy or pressure of nuclear matter under variations in density, temperature, and neutron-proton asymmetry [15], on the two-proton source. 2 1.1 Overview on particle femtoscopy In the 1950s, Hanbury-Brown and Twiss developed a technique called intensity interferometry [16, 17], which is often referred to as the Hanbury-Brown and Twiss (HBT) effect, to measure the size of celestial objects. The method works similarly to a two-slit interferometer experiment [18]: two photons emitted from two distant points on the celestial object interfere and produce interference patterns on their way to detectors (see Fig.1.2). Due to this interference, there is a correlation in two-photon intensities, which can be used to measure the size and shape of the photons emitter. Specifically, the technique was used to examine radio-wave sources in the galaxies Cygnus and Cassiopeia [16]. Additionally, the method was utilized to measure the angular diameter of Sirius [17]. The advantage of intensity interferometry over other methods for measuring stellar size was that it reduces the effect of phase shifts from photons’ wave function in the measurements [19]. Figure 1.2 Illustration of HBT interferometry experiment [20]. In intensity interferometry, a correlation function, denoted as 𝑅(𝑟1, 𝑟2) + 1, is obtained from the counts registered at detectors 𝐷1 and 𝐷2 [10]. It is given by: 𝑅(𝑟1, 𝑟2) + 1 = < 𝑁12 > < 𝑁1 >< 𝑁2 > , (1.2) Here, 𝑁12 represents the number of counts at which particles are observed at both 𝐷1 and 𝐷2, and 𝑁𝑖+1,2 represents the number of counts of particles registered at 𝐷𝑖=1,2. The HBT technique was later adopted in high-energy physics, where in Ref.[21], pion-pion correlation was used to determine the radius and lifetime of the hadronic matter produced during 3 a proton-proton collision. Subsequently, in nuclear heavy-ion collisions, it was demonstrated through proton-proton correlation that the two-particle correlation can be used as a tool to study the space-time evolution characteristics of emitting systems[5, 13]. The authors in Ref. [21] and the references within showed that the interaction between emitters can be neglected in the correlation for high-energy heavy-ion collisions. However, in nuclear reactions, the two-particle correlation is sensitive to the interaction between particles. Here, we refer to the technique as the two-particle correlation function (CF). Interestingly, the similarity between 𝑅 + 1 in Eq. (1.2) and the one in (1.1) is apparent, with the only difference being that 𝑅 for intensity interferometry ((1.2)) is a function of the particles’ position, while it is a function of relative momentum in (1.1). As an example, in Fig. 1.3, panel (A) illustrates the correlation function calculated for the 𝑑 − 𝛼 pair, while panel (B) displays the energy levels of 6Li. The observed peaks in panel (A) correspond to the decay of 6Li and are influenced by the nuclear interaction incorporated into the correlation function through the wave function obtained from solving the Schrödinger equation for the d-𝛼 scattering problem. The two-particle correlation function is theoretically approximated as the convolution between the squared relative wave function and the source function. In Fig. 1.3, a Gaussian source function with a half-width of 6 fm was employed. The relationship between panels (A) and (B) can be understood as follows: the presence of peaks in panel (A) corresponds to the energy levels of 6Li in panel (B). Specifically, the peaks at 𝑞 ≈ 42 and 80 − 100 MeV/c correspond to the energy levels of 6Li at 2.18 and 4.31 MeV, respectively. 4 Figure 1.3 Panel (A) represents d−𝛼 correlation function calculated using source function of radius of 6 fm. The peaks are due to nuclear interaction between the particles. Panel (B) shows the energy levels of different states in 6Li. The first peak in panel (A) corresponds to the decay of the 3+ state of 6Li at E= 2.18 MeV into d and 𝛼. While the second peak corresponds to the decay of 2+ state of 6Li at E= 4.31 MeV. 1.1.1 The Two-Body Scattering Problem In order to calculate the wave function for CF, we first need to study the scattering phase shift. This is achieved by solving the two-body scattering problem, which involves solving the radial Shrödinger equation in space coordinates: (𝐻𝑙 𝑗 − 𝐸)𝑢𝑙 𝑗 (𝑟) = 0, (1.3) where the Hamiltonian operator 𝐻 includes both nuclear and Coulomb interactions, as well as a centrifugal term. The radial wave function is represented by u(r), while E is the center of mass energy. We integrate (1.3) with the boundary conditions 𝑢(𝑟 → 0) = 0 and 𝑑𝑢 𝑑𝑟 |𝑟→0 ≠ 0, and 𝑢(𝑟 → ∞) → 𝐶𝐹 + 𝐷𝐺, where 𝐹 and 𝐺 are regular and irregular Coulomb wave functions, respectively ( more details about coulomb wave function, we refer readers to the following references: [22, 23], and the factors 𝐶 and 𝐷 are functions of phase shifts. We can then solve for the phase shifts using a relation known as the logarithmic derivative, which compares the derivative of the logarithm of the solution to the radial Shrødinger equation inside the nuclear potential and the asymptotic solution 𝑢𝑎𝑠𝑦 ∝ 𝐶𝐹 + 𝐷𝐺 at a point 𝑟 = 𝑅𝑎 where the strong interaction approaches zero. The logarithmic derivative is written as: 𝑑 log(𝑢) 𝑑𝑟 |𝑟=𝑅𝑎 = 𝑑 log(𝑢𝑎𝑠𝑦) 𝑑𝑟 |𝑟=𝑅𝑎 . (1.4) 5 Using this relation, we can solve for 𝐷 𝐶 = tan(𝛿), and from there, it is straightforward to extract the phase shifts (𝛿). In Chapter 2, we present the 𝛼 − 𝛼 and d–𝛼 phase shifts in Fig. 2.2 and 2.3, respectively. These phase shifts are calculated via this procedure. 1.2 Transport Model Transport models, such as the Boltzmann-Uehling-Uhlenbeck (BUU) model, are essential tools for obtaining physics information on the nuclear equation of state (EOS). When the EOS is well understood, it is possible to relate nuclear heavy-ion collision experiments (above saturation density 𝜌0) to astrophysical observations and gain insight into astronomical phenomena, such as the evolution of the universe after the Big Bang, supernovas, and the structure of neutron stars [15]. The BUU describes the time evolution of a single-particle phase-space distribution and nucleon- nucleon collisions [1, 14, 15]. The particles follow classical trajectories under the influence of a mean-field potential and subsequently undergo binary collisions. The model requires several parameters for the scattering cross-section and mean field, but only those that fit the elementary data and microscopic theory are used. This implies that the model is a phenomenological model, which allows for flexible tuning of some input parameters until the simulations match the experimental data. By adjusting these parameters, we can constrain the EOS. For example, in Ref. [14], the BUU transport model was successfully applied to calculate the two-proton correlation function, and it was found that the two-proton correlation is significantly sensitive to the in-medium cross-section. Additionally, in Ref. [24], the transport model was used to study the effect of the nuclear symmetry energy, a measure of the difference between the binding energy of an asymmetric and symmetric nucleus (equal numbers of protons and neutrons), on two-nucleon correlations in heavy-ion collisions. The density dependence of the nuclear symmetry energy strongly affects the nucleon emission times in these collisions, leading to larger values of two- nucleon correlation functions for those symmetry energies. Therefore, two-nucleon correlations can be used as tools to extract information about nuclear symmetry energies. Moreover, knowledge about the density dependence of the nuclear symmetry energy is an important tool to understand the structure of radioactive nuclei and the nucleosynthesis during the presupernova evolution of 6 massive stars and the mechanisms for supernova explosions, to name a few. Furthermore, simulating heavy-ion collisions with BUU can provide access to several observ- ables that can be used to extract information about the EOS. In this thesis, specifically in Chapter 6, we focus on using BUU simulations to test the sensitivity of the source function to momentum- dependent stiff and soft EOS for the p-p correlation function measured in low-energy heavy-ion collisions. We calculated angle-averaged and quadrupole component source functions and observed that the quadrupole component of the source function is more sensitive to the momentum-dependent EOS. 1.3 Thesis organization In Chapter 2, we discuss a framework for calculating a two-particle correlation function using a parameterized Gaussian source. We first calculate the pair’s relative wave function and subsequently deduce the correlation as a convolution of the squared wave function and the source function. The framework is applied to the case of low-velocity alphas. In Chapter 3, we develop a deblurring method that utilizes the inverse Richardson-Lucy (RL) algorithm, originally developed in optics for optical image restoration. We apply it to restore the two-particle source function for the case of d–𝛼 correlation. Additionally, we demonstrate that considering three pairs simultaneously, such as 𝛼-𝛼, d–𝛼, and d–d, measured simultaneously, would help better constrain the geometry of the emitting system. In Chapter 4, we apply the deblurring technique to restore the decay energy spectrum measured in the three-body decay of the 26O nucleus. In Chapter 5, we develop a machine learning framework, a deep neural network (DNN), for classification. We specifically demonstrate that the DNN classifier is an important tool that can be used to identify resonance peaks/states in the measured spectrum. In Chapter 6, we carry out BUU simulation for low-energy heavy ion reactions to compute single proton distributions, which are used to simulate two-proton source functions. Subsequently, we use the source function and proton-proton kernel to construct a two-proton correlation function. We demonstrate that when BUU and deblurring source functions are compared, we can test the 7 sensitivity of the source to the EOS. The thesis is concluded with a general summary and outlook in Chapter 7." 8 CHAPTER 2 LOW-MOMENTUM CORRELATIONS OF 𝛼 PARTICLES IN HEAVY-ION COLLISIONS This chapter discusses low-momentum alphas’ correlation in heavy-ion collisions. We use Koonin- Pratt theoretical approximation to analyse those correlations for Gaussian parametrized emitting source. Section 2.2 provides an introduction to the two-particle correlation function (CF) and an overview of the calculations for the two-particle wave function and the two-component source model. The results of 𝛼-𝛼 and d–𝛼 CF calculations are discussed in Section 2.3. We conclude in Section 2.4. 2.1 Introduction Particle correlations at low relative velocities in heavy-ion collision have provided important insights into geometry and time development of the final stages of reactions [11, 12, 5, 13], as well as into their phase-space [25] and thermodynamic characteristics [26]. The more abundant the particle species, the easier it is to measure the correlations for the species and more relevant such correlations for overall understanding of the reaction. In low energy collisions, alpha particles are quite abundant, compared to other particles, so it would be natural to exploit correlations involving these particles, in particular 𝛼–𝛼, to learn about the collisions. However, the basic theoretic infrastructure for the purpose has been lacking and is delivered here. Information content in the correlations on a specific reaction could be augmented by simultaneously analysing correlations within any pair in a set of measured species, or correlation matrix. We discuss that, in the context of the 𝛼 particles, for the combination of 𝛼 particles and deuterons. The features of the geometry and time development of a reaction, that the correlations may provide access to, can be summarized [27] in terms of relative distribution of emission points for the two particles, in their center-of-mass velocity frame, 𝑆(r), where r is the relative distance. For convenience, 𝑆 can be normalized to 1 under integration over space. Depending on the physics in the two-particle system and measurement details, the correlation may be still blind to some aspects of the source function 𝑆 and these often include the large 𝑟 region. The particles emitted 9 far away from each other tend not to impact each other much. Accordingly, source functions are often represented in terms of a combination of the short-range and long-range contributions, 𝑆(r) = 𝜆 𝑆short(r) + (1 − 𝜆) 𝑆long(r) [28]. An effort is made to determine the functional form of 𝑆short, but for 𝑆long(r) contributes only in an integral form, in that the overall strength of the source at short distances is incomplete, 𝜆 < 1. For historical reasons, specifically the Hanbury-Brown and Twiss (HBT) interpretation of pion-pion correlations, the parameter 𝜆 is often termed the incoherence factor [28]. We will dedicate some attention to 𝜆 in this chapter. To connect source functions (SF), 𝑆 in the reactions, to correlation functions (CF) 𝐶, one needs to understand relative wave functions within the studied particle pairs. The 𝛼 − 𝛼 correlations were measured in the past [29], but the basic infrastructure to interpret them was missing. To infer the wave functions, and develop such an infrastructure, scattering phase shifts need to be studied. These studies can inform on the impact of strong interactions on the wave functions and the probability density in the relative distance between the particles. Ultimately, deviations of the probability density from an asymptotic average serve as a lens that maps features of SF onto CF, as schematically indicated in Fig. 2.1. Any resonances are associated with strong modifications of the probability density at short relative distances. Coulomb interactions play a role too [5], depleting the wave function at short relative distances for low asymptotic relative velocities. In the end, the correlations of charged particles can be used to learn about SF even when strong-interaction modifications of the relative wave functions are minor. Needless to say, the charged particles are also easier to measure precisely. As to the resonances, there might be some narrow states in heavier products from a reaction, that include narrowly constrained pair of investigated particles among decay products [30], confusing the interpretation of two-particle SF-CF relation interpretation. These need to be separately understood and simulated. 10 Figure 2.1 Interactions and symmetrization, affecting the relative state of two particles on their way to the detectors, and encoded in the relative wave function, allow to infer characteristics of particle emission from correlation functions at low relative velocities. When correlations within one type of particle species are measured, so are correlations with other species emitted from about the same reaction region, Fig. 2.1, obviously more information can be gained about that region. We illustrate the principal potential of a joint analysis within the example of 𝛼 − 𝛼 and d – 𝛼 pairs measured in very similar reactions. 2.2 Two-Particle Correlation Functions In this section, we discuss the experimental and theoretical sides of two-particle CF in reactions with multiparticle final states. CF modeling involves folding of two-particle wavefunctions with relative sources for the particles. 2.2.1 Correlation Experiment and Theory The two-particle correlation function that is sought in an experiment is the ratio [31, 32, 29] of two-particle coincidence observation probability density, 𝑃(p1, p2), to the product of single particle observation probability densities, 𝑃1(p1,2), 𝐶 (q, K) = 𝑁 𝑃(p1, p2) 𝑃1(p1) 𝑃1(p2) . (2.1) Here, p1,2 are the momenta of the two particles, 1 and 2, K = p1 + p2 is the total pair momentum and q is the relative momentum in the pair center of mass. The factor 𝑁 is a normalization constant chosen in such a way that 𝐶 (q, K) approaches unity for large 𝑞. With the normalization to be set in (2.1), the probability densities may be taken unnormalized and equal to the invariant differential yields, 𝑃1(p) = 𝐸 d3 𝑁 d𝑝3 and 𝑃(p1, p2) = 𝐸1 𝐸2 [31]. d6 𝑁 1 d𝑝3 2 d𝑝3 11 On the theoretical side, using reduction formulas and the assumption of the relative velocity in the pair being small compared to the characteristic scales for the system reminder, the correlation function may be expressed as a convolution of the square of the relative scattering wave function Ψ(r, q) with a relative source 𝑆(r) [19, 5, 33, 34]: ∫ 𝐶 (q) = 𝑑3𝑟 𝑆(r) |Ψ(r, q)|2 . (2.2) The above relation is called the Koonin equation. The wave function Ψ is one with outgoing boundary conditions. The vector q is the asymptotic relative momentum, 𝑞 = 𝜇(v1 − v2), where 𝜇 is the reduced mass and v1,2 are velocities of the two particles in the pair. For identical particles, (cid:1). The vector r is the or in the center of mass, the relative momentum reduces to q = 1 2 (cid:0)p1 − p2 relative position within the pair, when the last member of the pair decouples from the rest of the system, and 𝑆 is the distribution of those positions. At large 𝑟, the average value of the probability density |Ψ|2 is 1. Deviations of the modulus square from 1 at shorter distances, regulated by q, allow to probe the shape of 𝑆. The particular features include resonances that correspond to peaks in the probability density at short distances, Coulomb repulsion yielding depletion in the probability density at short 𝑟 and low 𝑞 and symmetrization for identical particles, yielding oscillations in |Ψ|2. The interplay of the features in |Ψ|2, dependent on q, with features of 𝑆, yields structures in 𝐶 (q) that can be used to learn on 𝑆(r) [32, 27]. From 𝛼 − 𝛼 Phase Shifts to Wave Function In order to access 𝑆(r) using the Koonin-Pratt equation, the square of the relative wavefunction, |Ψ|2, is needed. The wave functions are constrained by any measured phase shifts. The latter determines the large-𝑟 behavior of partial wave functions, but ambiguities can remain for shorter-𝑟 restrained only by the 𝑞-dependence of the phase shifts. Fortunately the typical sources are wide enough for the ambiguities to have quite tolerable impact on the 𝑆-inference. In fact, at times, the purely asymptotic form of the wave functions may be used in coarse analyses [35, 36]. Notably, the asymptotic conditions need to be switched between the wavefunctions describing the scattering experiment, where the conditions are of an outgoing spherical wave, to those for wavefunctions in 12 the Koonin equation, where they are of an incoming wave. To generate changes in the wavefunction, to match the measures phase shifts, strong-interaction potentials are adjusted at short distances. The steps that we undertake for the 𝛼 − 𝛼 pairs are similar to those for the d – 𝛼 pairs later on, so we keep some options open even when they are not needed for 𝛼 − 𝛼. The experimental 𝛼 − 𝛼 phase shifts for 𝑠 and 𝑑 waves, from measurements of Afzal et al. [37], are shown in Fig. 2.2. We exploit them to constrain the relative 𝛼 − 𝛼 wave function by solving the Schrödinger Equation (SE) (cid:20) − ℏ2 2𝜇 ∇2 + 𝑈 − 𝐸 (cid:21) Ψ = 0 (2.3) and matching the asymptotic behavior in the partial waves to the phase shifts. Here, 𝜇 is the reduced mass for the system, 𝐸 = ℏ2𝑞2/2𝜇 and 𝑈 is the interaction potential consisting of the Coulomb potential and adjusted nuclear potential. Benefiting from angular momentum conservation and Figure 2.2 𝑆-wave (top panel) and 𝐷-wave (bottom panel) phase shifts for an 𝛼 − 𝛼 system as a function of the kinetic energy in the center of mass. The symbols represent the data by Afzal et al. [37] and lines represent our calculations with a Wood-Saxon form for the strong-interaction potential. Resonance behaviors can be observed at ∼ 0.091 and ∼ 3.04 MeV, in the 𝑠- and 𝑑-wave, respectively. azimuthal symmetry around the asymptotic relative momentum 𝑞, we can decompose the wave function into partial waves, Ψ(r, q) = 1 𝑞𝑟 ∑︁ ℓ 𝑖ℓ (2ℓ + 1) 𝑢𝑞 ℓ (𝑟) 𝑃ℓ (𝑐𝑜𝑠𝜃) , (2.4) 13 05100100S-wave (deg)Our calculationAfzal et al., 19690510ECM (MeV)050100D-wave (deg) Table 2.1 Wood-Saxon potential parameters (c.f. Eq. (2.6)), arrived at by fitting experimental phase shifts for the 𝛼 − 𝛼 system [37]. For describing 𝑠-wave data a single attractive Wood-Saxon term was sufficient, while the 𝑑-wave data required a combination of an attractive and a repulsive term. Partial Wave 𝑠 𝑑 Potential Parameters Attractive Term 𝑅0 (fm) 3.636 2.523 𝑉0 (MeV) -23.398 -38.838 𝑎0 (fm) 0.676 0.738 Repulsive Term 𝑅1 (fm) 𝑉1 (MeV) 𝑎1 (fm) 20.000 2.040 0.547 where 𝜃 is the angle between r and q and 𝑃ℓ are Legendre polynomials. The radial wave functions 𝑢𝑞 ℓ satisfy the radial equations − ℏ2 2𝜇 d2 d𝑟 2 𝑢𝑞 ℓ (𝑟) = (cid:18) 𝐸 − ℏ2 ℓ(ℓ + 1) 2𝜇𝑟 2 (cid:19) − 𝑈ℓ (𝑟) 𝑢𝑞 ℓ (𝑟) , where we assume the potential 𝑈ℓ to be local in each partial wave. The potential in each partial wave, with ℓ suppressed, is of the form 𝑈 (𝑟) = 𝑉𝑐 (𝑟) + 𝑉0 𝑓 (𝑟, 𝑅0, 𝑎0) + 𝑉1 𝑓 (𝑟, 𝑅1, 𝑎1) +𝑖𝑊0(𝐸) 𝑓 (𝑟, 𝑅𝑤, 𝑎𝑤) . (2.5) (2.6) Here, 𝑉𝑐 is a Coulomb potential regularized at short distances where the strong-interaction potential becomes dominant. The Woods-Saxon form factor is 𝑓 (𝑟, 𝑅, 𝑎) = 1 1 + exp 𝑟−𝑅 𝑎 . (2.7) The energy-dependent imaginary potential is only used for the d – 𝛼 system. The parameters of the potentials from fitting the phase shifts in the 𝛼 − 𝛼 system are given in Table 2.1 and the description of the phase shifts with these parameters is illustrated in Fig. 2.2. For the relative velocities that are of interest to us and higher angular momenta, the 𝛼 − 𝛼 pair is not going to penetrate the centrifugal to reach distances where they can interact strongly, so we can use pure relative Coulomb wave functions there. In constructing the wave function for the Koonin equation, we make use of the Coulomb scattering wave function with outgoing boundary 14 conditions, Ψ𝑐, and write the scattering wave function as Ψ = Ψ𝑐 + (Ψ − Ψ𝑐) . (2.8) Then we only decompose the second difference term in the angular momentum states and end up with summing over only two terms there. 2.2.2 d – 𝛼 Phase Shifts and Wave Function Deuteron spin couples to spatial angular momentum, complicating the analysis of the scattering wave function. Still that analysis is simplified by the net angular momentum and parity conservation. At low relative velocities, coupling between partial waves tends to be weak and vanishes exactly for the minimally nonlocal potential in Eq. (2.6). Under those conditions, the experimental analysis of cross sections for the phase shifts and their theoretical description largely parallels the analysis with no spin. With the coupling the wave function components are Ψ𝑀 ′ 𝑀 (r, q) = √ 4𝜋 𝑞𝑟 √ ∑︁ 𝐽ℓ 2ℓ + 1 𝑖ℓ ⟨ℓ01𝑀 |ℓ1𝐽 𝑀⟩ ×⟨ℓ 𝑀 − 𝑀′ 1 𝑀′|ℓ1𝐽 𝑀⟩ ×𝑢𝑞ℓ𝐽 (𝑟) 𝑌ℓ 𝑀−𝑀 ′ (𝜃, 𝜙) , (2.9) with angles relative to q. In modeling the correlation function, modulus squared of the pair wave function is employed, c.f. (2.2), and averaging and summing over spin directions is employed. This principally can add complexity to the angular momentum structure (see Ref. [38] for insights in another context). However, in this work we consider only angle-averaged correlations, both on experimental and theoretical side. Under such circumstances, only angle-averaged modulus square of the wave function is needed [27], which greatly simplifies the structure that needs to be convoluted with the source, itself angle-averaged: ∫ dΩr |Ψ𝑀 ′ 𝑀 (r, q)|2 1 2𝑆 + 1 = ∑︁ 1 4𝜋 𝑀 𝑀 ′ 1 3(𝑞𝑟)2 (2𝐽 + 1) |𝑢𝑞ℓ𝐽 (𝑟)|2 . ∑︁ ℓ 𝐽 (2.10) Wood-Saxon potential parameters for ℓ = 0, 1, 2 waves, have been determined by Boal and Shillcock [19], by fitting low-energy experimental phase shifts by McIntyre and Haeberli [39]. We modified 15 Table 2.2 D-Wave Wood-Saxon potential parameters (c.f. Eq. (2.6)), arrived at by fitting low-energy experimental phase shifts for the d−𝛼 system [39]. For describing 3D3-wave data a real Wood- Saxon was sufficient, and for describing the 3D2 and 3D1-wave data it was necessary to add a weak imaginary term. The potential parameters for other low-ℓ partial waves can be found in [19]. Partial Wave 3D1 3D2 3D3 Potential Parameters real Term 𝑅0 (fm) 4.082 2.916 2.765 𝑉0 (MeV) -13.019 -31.000 -43.045 𝑎0 (fm) 0.401 0.638 0.730 Imaginary Term 𝑊0 (MeV) -8.364 -3.040 𝑅𝑤 (fm) 3.082 2.082 𝑎𝑤 (fm) 0.312 0.612 𝐶 (MeV2) 0.85 0.85 their parametrizations for 3D2 and 3D1 waves, see Table 2.2, to better fit the data, by including weak imaginary contributions, with 𝑊0(𝐸) = 𝑊0 (𝐸 − 𝐸th)2/[𝐶 + (𝐸 − 𝐸th)2] in Eq. (2.6). Here, 𝐸th is threshold energy for deuteron breakup. The real values of the phase shifts in 𝑑 waves, from measurements and calculations, are shown in Fig. 2.3. With nuclear and short-distance regularization of the Coulomb potential impacting only the few lowest partial waves, the series in (2.10) can be cast again as the summed up series for the pure Coulomb interaction with a correction for the lowest waves added. Figure 2.3 Real part of 𝑑 − 𝛼 scattering phase shifts in ℓ = 2 states. The symbols represent the measurements of Ref. [39] and the lines represent our calculations. 2.2.3 Source Function The source function 𝑆 in Eq. (2.2) is the distribution of the particles in relative displacement r at the time when the last particle in the pair interacts with the exterior, abruptly changing momentum [34]. The function can be deduced from transport calculations of the reactions [40] or parametrized, 16 051015Elab (MeV)0123 (rad)3d33d23d1 such as in the Gaussian form [19, 32, 41]. The normalized Gaussian source function, describing emission localized in space and time, is given by: 𝑆(𝑟) = √︃ 1 − 𝑟2 2𝑟2 0 𝑒 (2.11) (2𝜋𝑟 2 0 Here, 𝑟0 is the source radius (i.e., the source size). When particles are emitted from a localized )3 source over an extended period, the resulting source function becomes anisotropic and elongated in the direction of the pair velocity [32] and references within. Angle-averaged correlations only measure the angle-averaged sources [27], and in these, the prolonged emission gives rise to enhanced tails. For a source emitting particles with a finite lifetime 𝜏, the Gaussian source can be written([32] and references within): 𝑆(𝑟, 𝑡) = √︃ 1 − 𝑟2 2𝑟2 0 𝑒 − 𝑡2 2𝜏2 (2𝜋𝑟 2 0 𝜏)3 (2.12) Equation (2.12) provides information about the space-time evolution of the particle-emitting region, and 𝑡 is the time between emissions. Other source functions, such as the thermal source function, have also been used in the past. Pratt (1984) [13] studied the emissions of pions using the Wigner function 𝑓 (𝑟, 𝑡, 𝑝) = 𝛿(𝑟 − 𝑟0)𝑒(−𝑡/𝜏−𝐸 (𝑟,𝑃)/𝑇) by assuming the pion to be emitted from a spherical shell of radius 𝑟0, temperature 𝑇, and lifetime 𝜏, where 𝐸 is the energy of the pions. Figure 2.4 Dashed line superposed with crosses (green) represents normalized Gaussian source function, the solid line normalized source function of the form of reciprocal of hyperbolic cosine (i.e., ∝ 𝑐𝑜𝑠ℎ(𝑟 ′) ). Both functions are significant for low values of r and fall quickly at large values of r. But clearly, the Gaussian tail falls to zero much quicker than the tail of 𝑆(r) ∝ 1 1 𝑐𝑜𝑠ℎ(𝑟 ′) . 17 402002040r [fm]0.00000.00250.00500.00750.01000.01250.0150r2S(r) [fm1]GaussianS(r) 1/coch(r/rc) In this chapter, we utilize a simple time-independent and isotropic Gaussian form such as in Eq. (2.11). We cross-checked our calculation for the 𝛼 − 𝛼 CF using a reciprocal hyperbolic cosine source, denoted as 𝑆(𝑟) ∝ 1 cosh(𝑟 ′) , where 𝑟′ = 𝑟 𝑟𝑐 (see Fig. 2.7), where 𝑟𝑐 is a parameter. In Fig. 2.4, we compare the Gaussian source with the reciprocal hyperbolic cosine source as a function of relative distance. The latter converges to zero more slowly compared to the former. This reciprocal hyperbolic cosine function could be a good candidate for describing the emission source of alpha particles, which are abundant in low-energy reaction systems. 2.2.4 Two parameters source model (TPSM) As shown in Fig. 2.4, the model defined in (2.11) is significant at short relative distances and falls quickly at large relative distances r. Such source correspond to short-range emissions, but in some cases, long-distance emissions also matter. To consider both emissions, we represent the source function as a combination of both emissions in terms of short-range and long-range contributions (see also Section. 2.2 ), given by the equation [28]: 𝑆𝜆,𝑟0 (𝜆, 𝑟) = 𝜆𝑆𝑠ℎ𝑜𝑟𝑡 (𝑟, 𝑟0) + (1 − 𝜆)𝑆𝑤𝑖𝑑𝑒 (𝑟, 𝑟0), (2.13) where 𝜆 is the probability of detecting particles emitted at short relative distance (0 ≤ 𝜆 ≤ 1). 1 − 𝜆 is the probability of detecting particles emitted at wide relative distances. Regarding the source function at short and large relative distances, 𝑆𝑠ℎ𝑜𝑟𝑡 (𝑟, 𝑟0) and 𝑆𝑤𝑖𝑑𝑒 (𝑟, 𝑟0) we used Gaussian function. Thus, the source function becomes a two-parameter model (𝜆 and 𝑅0). Substituting this source in (2.2), we get; ∫ 𝐶 (𝑞) = 𝑑3𝑟𝜆𝑆𝑠ℎ𝑜𝑟𝑡 (𝑟; 𝑟0) (cid:16) |Ψ(r, q)|2 − 1 (cid:17) + 1. (2.14) The parameters 𝜆 and 𝑟0 play important roles in the calculation of the correlation function (CF). The parameter 𝜆 ensures the inclusion of short-range sources (i.e., fast emissions) and long-range sources (i.e., slow emissions) in the calculation. Notably, the effects of both parameters on the correlations are observed, as shown in Fig.2.5. Decreasing 𝜆 reduces both the height and width of the correlation, while decreasing 𝑟0 results in an increase in both aspects; this is unsurprising because the source function is proportional to 𝜆, and inversely proportional to the radius,𝑟0. Our 18 (a) (b) Figure 2.5 The plots illustrate the CF as a function of relative momentum (q) obtained from Eq. (2.14). In row (a): the top left panel represents 𝛼 − 𝛼 correlations for a fixed source radius, 𝑅0 =6 fm, with varying values of 𝜆. The top right panel shows different curves corresponding to different radii while keeping 𝜆 constant. The bottom row (b), left and right panels represent 𝛼 − 𝑑 correlations and follow a similar pattern to the top row (a) for varying 𝜆 and radii, respectively. The legends or titles of the plots provide the values of the parameters, 𝜆 and 𝑅0 used. Additionally, the insert zooms in on the second peak to offer a closer view. goal is to determine the combination of these parameters that reproduces the measured correlations. Additionally, it is easy to observe that Eq.(2.14) reverts to Eq. (2.11) when 𝜆 is set to unity. The results shown in Fig. 2.5 were obtained using Eq. (2.14). The left panel in Fig. 2.5a shows the 𝛼-𝛼 correlation function (CF) calculated for a fixed value of 𝑟0, and the different curves 19 50100150200q (MeV/c)2468101214d CFR0=5.9 fm=0.2=.25=0.3150100150200q (MeV/c)=0.31R0=6 fmR0=7 fmR0=7.5 fm correspond to different values of 𝜆. In the right panel, 𝜆 is fixed, and the curves represent different values of 𝑟0. Similarly, in Fig. (2.5b), the left panel shows the d-𝛼 CF for 𝑟0 = 5.9 fm, and the different curves correspond to different values of 𝜆. The right panel represents the d-𝛼 correlation for 𝜆 = 0.3 with different values of 𝑟0. In both cases, we observe resonance peaks from two-particle nuclear interactions; for the 𝛼 − 𝛼 CF (Fig. 2.5a), the first resonance peak at 𝑞 = 18 MeV/c is due to 𝛼 − 𝛼 S-wave interaction, while the second resonance peak at 𝑞 = 108 MeV/c is due to d-wave interaction. For the d–𝛼 CF (Fig. 2.5b), the first sharp peak at 𝑞 ≈ 42 MeV/c corresponds to the state 3D3 of 𝑑 − −𝛼, while the second peak at 𝑞 ≈ 84 MeV/c is due to the overlap of 3D2 and 3D1 states of d−𝛼 [42, 43]. Here, spectroscopic notation 2𝑆+1L𝐽 is used to denote the orbital angular momentum L, where 𝐿 = 2 stands for d-wave, while the upper-script number is 2𝑆 + 1 and the subscript number is the total angular momentum 𝐽. These quantum numbers were determined by employing the triangle selection rule, i.e., |𝐿 − 𝑆| ≤ 𝐽 ≤ 𝐿 + 𝑆, where 𝑆 is the spin number and 𝑆 = 1 for d – 𝛼 and 0 for 𝛼 − 𝛼 systems. In addition, the effect of long-range interaction (i.e., Coulomb interaction) is observed from the suppression in the CF at 𝑞 ≈ 0 MeV/c [5]. This effect is only observed for charged particle systems. 2.3 Results: 𝛼 − 𝛼 and 𝑑 – 𝛼 CF In this section, we discuss the 𝛼 − 𝛼 and d – 𝛼 CF measured in reactions 40Ar+27Al at E=44 MeV/A [29] and 36Ar+197Au at E=35 MeV/A [44]. We compare those correlations with theoretical calculation obtained utilizing Eq. (2.14). With this comparison, we extract the parameters for the emitting sources. In addition, we demonstrate that considering different particle pairs measured in the same reaction would be key for better constraining the geometry of the final state emission. It is important to mention that the correlations measured in an experiment, especially the structures in the low relative momentum range, are affected by detector resolutions. So far, these resolutions have not been used in C(q). However, in the rest of this chapter, we will discuss the measured CF. Therefore, it is necessary to apply the energy resolution in the calculated CF, so that we can reproduce all the features in the measured CF accurately. To account for that resolution in 20 the calculated correlation, we follow the smearing in q proposed in Ref. [42] by folding C(q) with a Gaussian function in q of width 𝜎𝐸 , adjusted to the energy resolution. We provide the value of 𝜎𝐸 used for each experimental correlation data that we analyze in this section and we will refer to it as the energy resolution parameter. 2.3.1 𝛼 − 𝛼 CF The correlations between alpha particles in low-energy reactions have not been well studied before, to the best of our knowledge. Here, we use the theoretical framework developed to analyze experimental alpha-alpha correlations. Fig. 2.7 shows the experimental and theoretical alphas correlations. The experimental data from Ref. [29] is shown as stars, while the theoretical CFs are shown as a dashed line ( green) for the Gaussian function and as a dashed line (blue) for the reciprocal hyperbolic cosine function 1 (i.e., 𝑆(𝑟) ∝ cosh(𝑟/𝑟𝑐) ), the energy resolution parameter,𝜎𝐸 = 1.9 MeV/c is used. Both source models closely reproduce the experimental data, and the corresponding parameters extracted from the fits are 𝜆 = 0.3 and 𝑟0 = 5.8 fm for the Gaussian function, and 𝜆 = 0.224 and 𝑟𝑐 = 3.5 fm for the reciprocal hyperbolic cosine source function. Those sources are plotted together with squared relative wave function in Fig. 2.6. 21 Figure 2.6 We plot the squared 𝛼 − 𝛼 wave function (|Ψ(r, q)|2) as a solid red line, along with the reciprocal hyperbolic cosine and Gaussian source functions represented by double-dashed and dash-dotted lines, respectively. These functions are plotted as a function of the relative distance 𝑟. The squared wave function, |Ψ(r, q)|2, corresponds to the 𝛼 − 𝛼 pair at relative momentum of 𝑞 ≈ 108 MeV/c. The shaded region indicates the range in the relative distance where both the wave function and source function are significant, i.e., the region of short relative distances. The parameters used to construct the sources can be found in the text. As mentioned earlier, in Fig. 2.6, we can observe that both the wave function and the emitting source function are important for small values of relative distances. At large distances, |Ψ|2 (shown as a solid red line) approaches unity, while the Gaussian and reciprocal hyperbolic cosine sources (represented by dash-dotted and double-dashed lines, respectively) decay to zero. At these distances, the interplay between |Ψ|2, which is dependent on q, and the source results in a correlation of 𝐶 (𝑞) = 1. In Fig. 2.7, the resonance peaks correspond to the decay of 8Be states are seen, where the first peak around 𝑞 = 18 MeV/c corresponds to the decay of the 𝐸 = 0.092 MeV ground state of 8Be [7], while the peak located at 𝑞 ≈ 108 MeV/c corresponds to the decay of the 𝐸 = 3.04 MeV first excited state of 8Be [7]. The energies of these decay states agree with the resonance energies observed in the phase shifts shown in Fig. 2.2. In addition, a peak around q=50 MeV/c is observed in the data, we discuss this peak in the next subsection. 2.3.1.1 Resonance peak at 𝑞 = 50 MeV/c In the experimental data, another peak appears at 𝑞 = 50 MeV/c; this peak can not be easily produced in the theoretical 𝛼 − 𝛼 CF such as shown in Fig. 2.7 because it is not directly related 22 010203040 r [fm]012345678||2 ||2 at q=104 MeV/c0.00000.00010.00020.00030.00040.00050.00060.0007S(r) [fm3]S(r) 1/coch(r/rc) S(r) (Gaussian) Figure 2.7 Figure compares calculated 𝛼 − 𝛼 CF and experimental data. Dashed blue line and dashed green are calculated CF for Gaussian and 1/Cosh source functions respectively in Fig. 2.6. The dots are experimental data from 40Ar+27Al reaction at 𝐸 = 44 MeV/A [29]. The insert shows more detail around the peak at 104 MeV/c. to any of the states of 8Be. This peak has been attributed to the decay of the 𝐸 = 2.42 MeV state of 9Be (𝐽 𝜋 = 5/2−) into 2𝛼 and an undetected neutron [45, 7, 29]. To obtain the correlation around this peak, we use particle number density distributions from the thermodynamics model as inputs to Eq. (2.15). First, we determine the proton and neutron chemical potentials at thermal equilibrium by solving coupled equations for the conservation number of protons and the number of neutrons in the system [46]. Then, we generate alpha particles through the sequential decay of 9Be→ 𝛼+5He, 5He→ n+𝛼. Subsequently, we use the generated alphas’ kinematics and proton and neutron chemical potentials to simulate particle number distributions, which are needed to construct 𝛼-𝛼 correlations, as shown in Fig. 2.8a. 23 (a) (b) Figure 2.8 Panel (a) represents alphas CF for the peak at around q=50 MeV/c, simulated using the thermodynamic model and Eq. (2.15). Panel (b) is the same as Fig. 2.7, but now includes the peak at 𝑞 = 50 MeV/c. 𝑑𝑁𝛼1−𝛼2 𝑑3𝑃𝛼1 𝑑3𝑃𝛼2 = 𝑅(𝑞) 𝑑𝑁𝛼1 𝑑3𝑃𝛼1 1 − 𝛼 2 𝑑3𝑃𝛼 2 𝑑𝑁 𝛼 2 𝑑3𝑃𝛼 2 𝑑𝑁 𝛼 𝑑3𝑃𝛼 1 𝑑𝑁 𝛼 1 𝑑3𝑃𝛼 1 𝑅(𝑞) = 𝑑𝑁𝛼2 𝑑3𝑃𝛼2 , . (2.15) where 𝑑𝑁 𝛼 𝑑3𝑃𝛼 1 1 − 𝛼 2 𝑑3𝑃𝛼 2 , 𝑑𝑁 𝛼 1 𝑑3𝑃𝛼 1 , and 𝑑𝑁 𝛼 2 𝑑3𝑃𝛼 2 are particle number distributions, obtained from thermody- namics model (see Appendix B). We can see the analogy between Eq. (2.15) and Eq. (2.1); in Eq. (2.15), the numerator and denominator correspond to the probability of coincidence emission and probability of single particle emission respectively. The thermodynamics parameters used in the calculation are temperature, 𝑇 = 12 MeV, mass, 𝐴 = 67, normal density, 𝜌0 = 0.15, neutron chemical potential, 𝜇𝑛 = −15.98 MeV and proton chemical potential, 𝜇𝑧 = −20.553 MeV. Fig. 2.8a shows the resulting correlations for the alphas pair produced from the sequential decay of the 𝐸 = 2.42 MeV state of 9Be into 2𝛼 +5 He via the resonance of 5Heg.s, where the neutron was not detected. Finally, the results in Fig. 2.8a are combined with 𝛼 − 𝛼 CF such as in Fig. 2.7 to give a full theoretical calculation of the 𝛼 − 𝛼 CF which reproduce well the experimental data as shown in the Fig. 2.8. This is done by replacing the portion of the 𝛼 − 𝛼 correlation function spanning 30-60 24 050100150q [MeV/c]010C(q)data (Ghetti et al. (2004)our calculation MeV/c with the correlation shown in Fig. 2.8a, scaled by a factor of .85 to match the calculated correlation function shown in Fig. 2.7. 2.3.2 d – 𝛼 CF In this section, we discuss the measured d–𝛼 correlation functions in 40Ar+26Al reactions at 𝐸/𝐴 = 44 MeV [8] and 36Ar+197Au reactions at 𝐸/𝐴 = 35 MeV [44]. For this system, we obtained the theoretical correlations by calculating the wave function using potential parameters from Ref. [32], which we revised to include an imaginary potential (refer to Eq. (2.10) and Table. 2.2 for the specific parameters we used). Studying the d–𝛼 correlation function provides insight into the geometry of the final state emission and serves as a valuable tool for understanding the decay of 6Li states. Fig. 2.11a and Fig. 2.11b depict the calculated d–𝛼 correlation functions plotted alongside the experimental data for the respective reactions. The best-fit parameters of the model source were 𝑟0 = 5.5 fm and 𝜆 = 0.6 for the plot in Fig. 2.11a (i.e., 40Ar+26Al reactions at 𝐸/𝐴 = 44 MeV). In Fig. 2.11b, there are two datasets: the triangles represent data from peripheral collision, while the stars represent data from a central collision in 36Ar+197Au reactions at 𝐸/𝐴 = 35 MeV. The best parameters for the peripheral collision are 𝑟0 = 5.00 fm and 𝜆 = 0.61, and for the central collision, they are 𝑟0 = 6.20 fm and 𝜆 = 0.61. From Fig. 2.9, it is clear that the calculated correlation functions reproduce both peaks in the experimental data. The first peak at 𝑞 = 42 MeV/c corresponds to the decay of the 𝐽 𝜋 = 3+ 6Li excited state at 𝐸 = 2.186 MeV [19, 44, 8], and the second peak between 80 − 100 MeV/c relative momentum corresponds to the overlap of the 𝐸 = 4.31 MeV (3D2) and 𝐸 = 5.65 MeV (3𝐷1) excited states of 6Li [19, 44, 8]. In the past, reproducing the shape of the d–𝛼 correlation function has been a very challenging task. Most of the time, the first resonance peak was successfully predicted, while the second one was overpredicted [19, 42, 44]. In Ref. [47], the authors attempted to solve this issue by applying both source geometry and collective motion in the d–𝛼 correlation function in Xe+Au collisions. However, the method still overpredicted the second peak. Here, to address this issue, we investigated the d–𝛼 scattering phase shift with close attention 25 paid to the 3D2 and 3D1 states, as they are responsible for the second peak in the d–𝛼 CF. We introduced an energy-dependent imaginary Wood-Saxon potential for the 3D2 and 3D1 states of the d–𝛼 system and chose the parameters of the potentials (see Table. 2.2) to accurately reproduce the experimental d–𝛼 scattering phase shifts, as shown in Fig. (2.3). (a) (b) Figure 2.9 The plots represent d – 𝛼 CF as a function of relative momentum; in panel (a) the dashed line is theoretical calculation and points is experimental data from 40Ar+27Al reaction at 𝐸 =44 MeV/A [8], and panel (b) is the same as (a) with experimental data (stars and triangles) from 36Ar+197Au reaction at 𝐸/𝐴 =35 MeV [44]. The energy resolution parameters used are; 3.5 MeV/c and 1.9 MeV/c for left and right panels, respectively. 2.3.3 𝛼 − 𝛼 and d −𝛼 measured simultaneously In the previous section, we discussed the correlation functions, namely 𝛼 − 𝛼 and 𝑑 − 𝛼, and extracted the corresponding source functions. In this section, our focus is on comparing these source functions for pairs measured in the same reaction. Additionally, we demonstrate that the d–𝛼 source function could be directly estimated from 𝛼 − 𝛼 and d–d sources if the correlation function three pairs were measured in one reaction. For demonstration, we plotted in Fig. 2.11b the source functions for the three pairs utilizing the experimental data from the 40Ar+27Al reaction for 𝛼 − 𝛼 and d–𝛼 correlations, and from the 197Au +40 Ar reaction at 𝐸/𝐴 = 60 MeV for the d-d correlation. When the source functions from three pairs measured in the same reaction are studied, it may play a crucial role in constraining the geometry of the final state emission. Fig. (2.10), (C) presents the 𝛼-𝛼 and d–𝛼 source functions multiplied by the squared relative distance for both pairs measured in the 40Ar+27Al reaction at 𝐸 =44 MeV/A. The dots represent the 26 50100150q [Mev/c]24C(q)theorydata (Ghetti et al.,2004)50100150q [Mev/c]024681012C(q)R0=6.20 fm, =0.61R0=5.00 fm, =0.6central (Zhu,(1992))peripheral (Zhu,(1992)) 𝛼-𝛼 source, while the solid line represents the d–𝛼 source. As depicted in the figure, the parameters for the 𝛼-𝛼 source are 𝑟0 = 6 fm and 𝜆 = 0.3, and for the d–𝛼 source, 𝑟0 = 5.5 fm and 𝜆 = 0.55. This results in the d–𝛼 and 𝛼-𝛼 source functions being different. However, the radii of both sources are approximately equal. Panel (a) and (b) illustrates the correlation function for 𝛼 − 𝛼 and d–𝛼 pairs, respectively, that correspond to the source functions in (c). Figure 2.10 Figure illustrates 𝛼 − 𝛼 correlation function, panel (a), and d–𝛼 correlation, panel (b), both measured in the same reaction, 40Ar+27Al reaction at 𝐸 =44 MeV/A [8]. Panel (c) displays the source functions: solid line displays d–𝛼 (red) source function and dots (green) shows 𝛼-𝛼 source function. To calculate the 𝑑 − 𝑑 correlation function (CF) depicted in Fig. 2.11a, we utilize Eq. (2.2) and construct the wave function using the Wooden Saxon potential parameters from Ref. [32]. The recalculated 𝑑 − 𝑑 correlations are then compared with experimental data obtained from the 197Au +40 Ar reaction at 𝐸/𝐴 = 60 MeV, where 𝐸1 + 𝐸2 = 75 − 125 MeV, where 𝐸1 and 𝐸2 stand for energies of the first and second deteuron particle, respectively. This allows us to determine the 𝑑 − 𝑑 source radius, denoted as 𝑅𝑑𝑑, which has been found to be 5 fm. Notably, this value is consistent with those reported in Ref. [32]. Our calculations (solid line) in Fig. 2.11a demonstrate 27 50100150200q [MeV/c]51015C(q)(a)data (Ghetti, et al., (2006))theory, pair50100150q [MeV/c]24C(q)(b)data (Ghetti, et al., (2006))theory, d pair402002040r [fm]0.000.010.020.03r2S(r) [fm3](c)d good agreement with the corresponding experimental data points (squares). Fig. 2.11b displays the source functions as a function of relative distance, clearly illustrating that the width of the 𝑑 − 𝑑 pair (dashed line) is smaller than that of the d–𝛼 pair (solid line) and the 𝛼 − 𝛼 pair (dots). We assume that the distribution of each of the two particles in the system is described by a Gaussian source function with widths 𝑅1 and 𝑅2 for particles 1 and 2, respectively. To obtain the two-particle relative distribution (relative source function) in the system, we perform a convolution of the two Gaussian functions, resulting in another Gaussian function with a width of 𝑅12 = √︃ 𝑅2 1 + 𝑅2 2, also referred to as the source radius. Thus, by substituting the indices 1 and 2 with the particle symbols (i.e., 𝛼 or 𝑑), we can express the widths (i.e., radii) of 𝑑–𝑑, 𝛼–𝛼, and 𝑑–𝛼 relative sources as √︃ 𝑅𝑑−𝑑 = 𝑑 + 𝑅2 𝑅2 𝑑, 𝑅𝛼−𝛼 = √︁𝑅2 𝛼 + 𝑅2 𝛼, and 𝑅𝑑−𝛼 = √︃ 𝑅2 𝛼 + 𝑅2 𝑑, respectively. Therefore, we deduce the mathematical relationship that connects the radii of the three pairs: 𝑅2 𝑑−𝛼 = 1 2 (cid:2)𝑅2 𝛼−𝛼 + 𝑅2 𝑑−𝑑 (cid:3) . (2.16) For example, if 𝑅𝛼−𝛼 = 6 fm and 𝑅𝑑−𝑑 = 5 fm, then 𝑅𝑑−𝛼 ≈5.5 fm. It would certainly be interesting to further validate this concept through simultaneous measurements of the pairs. In addition, it is important to note the difference between the measured 𝛼 − 𝛼 correlations presented in Fig. (2.10(a)) and 𝛼 − 𝛼 correlations shown in Fig. (2.8(b)). The correlation in Fig. (2.10(a)) corresponds to alphas’ total energy greater than 140 MeV. As can be observed, this correlation misses the low relative momentum part, including the resonance peak at q=18 MeV/c, which may affect the fit for the source function. We believe that this might explain the difference in source functions (see the parameters 𝑟0 and 𝜆 of the two sources discussed earlier in this section) for the two measured 𝛼 correlations. 28 (a) (b) Figure 2.11 Panel (a) shows the d–d correlation function; the solid line represents the calculated CF, and the squares represent d–d correlations measured in 197Au +40 Ar collision at 𝐸/𝐴 = 60 MeV for 𝐸1 + 𝐸2 = 75 − 125 MeV [32] (and references within). The source function used to calculate the d–d correlation is displayed as a dashed line in panel (b). The remaining curves in the panel represent source functions for d–𝛼 (dashed line) and 𝛼 − 𝛼 (dots) pairs measured in the 40Ar +27 Al reaction at 𝐸/𝐴 = 44 MeV [8]. 2.4 Conclusion We solved the Schr¨dinger equation (SE) for 𝛼 − 𝛼 scattering to obtain the scattering phase shifts, which were compared with experimental phase shifts to determine the best nuclear potential. Subsequently, the potential was used to determine the relative wave function of alphas. The squared relative wave function, along with the parameterized Gaussian source function, was employed in the Koonin-Pratt equation (Eq. (2.1)) to compute alphas correlations in low-energy reactions. We calculated the 𝛼 − 𝛼 correlation function (CF) and compared it to correlations measured in the 40Ar + 27Al reaction at 𝐸/𝐴 = 44 MeV [8] to estimate the parameters of the Gaussian source for this system. In that reaction, the parameters, radius, and 𝜆 for the 𝛼 − 𝛼 pair were found to be 5.80 fm and 0.30, respectively. Similarly, the d-𝛼 CF was analyzed using a two-parameter source model. Similar analyses have been carried out in Ref. [47] using different approaches, such as a two-parameter model and collective motion, and in [44] using a single-parameter Gaussian source model. However, our calculation reproduces the features in the measured d-𝛼 CF much better. Furthermore, we demonstrated that studying correlations between different particle pairs measured in one reaction helps to better constrain the geometry of the emission in the final state. We assumed the example of 𝛼 − 𝛼 and d−𝛼 measured in the same reaction (Fig. 2.10). Under the assumption 29 50100150q [Mev/c]0.40.60.81.0C(q)theorydata (Boal et al., 1990)05101520r [fm]0.0000.0010.0020.003S(r) [fm3]d dd that the 𝛼 − 𝛼, d- d, and d-𝛼 correlations are measured in the same reaction, we can easily predict the size of d−𝛼 if we know the sizes of 𝛼 − 𝛼 and d-d, but data is needed to accomplish this. This chapter primarily focuses on extracting one-dimensional source function from two-particle correlation functions, we are also interested to extend this work to three-dimensional source function information by analyzing three-dimensional particle correlations. Additionally, we may explore the deformation of the source in emission for distinguishable particles by expanding the source into angular components. While this phenomenon has been studied in high-energy reactions [48], it has received less attention in lower-energy cases. Once again, data is necessary to further investigate these aspects. 30 CHAPTER 3 DEBLURRING APPROACH In this chapter, we propose the use of the Richardson-Lucy (RL) optical deblurring algorithm for imaging source from the correlation function. We focus on the cases of deuteron-alpha (d–𝛼) correlations measured in a heavy-ion collision experiment. The main results of this chapter have been published and are accessible online in Ref. [49]. The chapter begins with a brief introduction, and in Section 3.2, we discuss the deblurring method. In Section 3.3, we apply the method to restore source function from two-particle correla- tion. Finally, the chapter concludes with a summary in Section 3.4. 3.1 Introduction In the previous chapter, the two-particle source functions have been parametrized in a Gaussian form and fitted to the correlation data. Alternatively, computing source function from the correlation measurements is an imaging problem, that principally, like elsewhere for imaging, invokes inversion and thus may suffer from instabilities. In this chapter, rather than applying an inversion directly, we take inspiration from optical deblurring which is an imaging problem too. One successful strategy there, that has been already ported into nuclear physics to restore decay energy spectra measurements [50], and cope with detector inefficiencies and reaction-plane uncertainties [51], is the Richardson-Lucy (RL) method [52, 53] that relies on the Bayes theorem and it follows an iterative procedure to find a self-consistent solution. Here, the algorithm only uses the correlation measurements and discretized transfer function, or kernel matrix (𝐾), as inputs. The RL method largely owes its success to the fact that it operates with strictly positive definite quantities, the probabilities. Conveniently, the corresponding quantities of interest within the KP formula (i.e., correlations and kernel matrix) are positive definite, even though the overall meaning of the KP formula differs from that providing context for the RL method. In maintaining the restored source function positive throughout the iteration procedure and by avoiding a direct 𝐾 inversion, serious singularity problems plaguing inverse problems are avoided. In addition, we test the RL method on the source restoration from the d-𝛼 correlation function calculated with a Gaussian source 31 function. We take both correlation functions without and with uncertainties. Those uncertainties are introduced in the correlation by adding Gaussian noise to the KP formula. In this chapter, we would like to recall the KP formula since it is important for us to discuss the deblurring method for restoring the source function from correlations, at 𝑞 −→ 0, the correlation function may be represented in terms of the KP formula: ∫ 𝐶 (q) = 𝑑3𝑟 |Ψ(−) q (r)|2 𝑆(r) ≡ ∫ 𝑑3𝑟 𝐾 (q, r) 𝑆(r) . (3.1) Here, Ψ(−) q is a 2-particle scattering wave function specified with incoming wave boundary con- ditions and asymptotically representing the center-of-mass relative momentum q. Possible spin indices are suppressed at this stage. The wave function normalization is such that the kernel in the KP relation, 𝐾 ≡ |Ψ|2, averages to 1 in the asymptotic zone of large 𝑟. The function 𝑆(r) is the probability distribution of particles 1 and 2 in their separation r in their center of mass, for the instant when they separate from the rest of the system and leave for the detectors. The ability to learn from the low relative-velocity correlations is further emphasized by subtracting unity from both sides of Eq. (3.1) and arriving at an equation for the 𝑅 correlation function [54]: 𝑅(q) = ∫ (cid:16) 𝑑3𝑟 |Ψ(−) q (r)|2 − 1 (cid:17) 𝑆(r) . (3.2) Provided the interaction within the particle pair is constrained, so that |Ψ| can be faithfully assessed for 𝑞 and 𝑟 of interest, then 𝑆 may be inferred from any structures in 𝑅. If the interaction is unknown, but generic assumptions on 𝑆 can be made, then the KP relation may be used to constrain the interaction between the particles. As illustrated in Fig. 3.1, the left panel shows the contour of the logarithm of d−𝛼 squared relative wave function, |Ψ|2. The features (bright spots) observed around 𝑞 = 42 and 80 MeV/c correspond to 3D3, and overlap of 2D1 and 2D2 interactions (we refer the readers to Chapter 2 for more details about those interaction). The right panel shows d−𝛼 correlations, obtained using the KP formula with a Gaussian source. The peaks around 42 and 84 MeV/c result from the interactions between deuteron and alpha particles, as shown in the left panel, and correspond to the decay of 6Li at E=2.18 MeV and the overlap states of 6Li at 4.31 and 5.6 MeV. The correlation function averaged over directions of q is related to the source function 32 Figure 3.1 The left panel displays a contour plot of the logarithm of the squared relative wave function |Ψ|2 of deuteron-alpha (𝑑 − 𝛼) pairs. The vertical axis shows the relative distance (𝑟) between the pairs, and the horizontal axis shows the relative momentum (𝑞). While, the right panel represents the 𝑑 − 𝛼 correlation function, obtained using the KP formula, in Eq. (3.3), also discussed in Chapter 2, for a Gaussian source function ( shown as stars in Fig. 3.2 ). The first and second peaks in the correlation represent the decay of 6Li at E=2.18 MeV and the overlap states of 6Li at 4.31 MeV (2+) and 5.6 MeV (1+1). These peaks correspond, respectively, to the bright spots observed in the left panel around 40 MeV/c and 84 MeV/c. averaged over the directions of the relative separation r: and 𝐶 (𝑞) = 4𝜋 ∫ ∞ 0 𝑑𝑟 𝑟 2 𝐾 (𝑞, 𝑟) 𝑆(𝑟) , 𝐾 (𝑞, 𝑟) = |Ψ(−) q (r)|2 , (3.3) (3.4) where the r.h.s. is the squared wave function is averaged over orientation of r relative to q. The squared wave function is given as [55]; |Ψ(−) q (r)|2 = ( 1 2𝑠 + 1 ) ∑︁ 𝐿 (2𝐿 + 1) 𝐴𝐿 (q, r)𝑃𝐿 (cos 𝜃), (3.5) where 𝐴𝐿 (q, r) = 𝑙𝑚𝑎𝑥∑︁ 𝑙′ 𝑚𝑎𝑥∑︁ (2𝑙 + 1)(2𝑙′ + 1) 𝑙=0 𝑙′=0 𝑅𝑒(𝑖𝑙−𝑙′𝑈∗ 𝑗 ′𝑙′ (q, r)𝑈 𝑗𝑙 (q, r)). 𝑗=𝑙+𝑠 ∑︁ 𝑗 ′=𝑙′+𝑠 ∑︁ 𝑗=𝑙−𝑠 𝑗 ′=𝑙′−𝑠 (2 𝑗 + 1) (2 𝑗 ′ + 1) (cid:169) (cid:173) (cid:173) (cid:171) 𝑙 𝑙′ 𝐿 0 0 0 2 (cid:170) (cid:174) (cid:174) (cid:172)    𝑙 𝑗 ′ 𝑙′ 𝐿 𝑗 𝑆 2    The symbols in () and {} are 3-j and 6-j symbol respectively. And 𝑈 𝑗𝑙 is the radial wave function, obtained by solving the Schrodinger equation with parameters of optical potential taken from 33 Chapter 2. For L=0 we reach to angle averaged squared wave function, |Ψ(−) q (r)|2 ( or refer to Eq. (2.10)). As may be apparent in Eqs. (3.1) and (3.2), the inference of 𝑆 from 𝐶 represents an imaging problem. In fact, for neutral pion pairs, with weak strong-interaction effects within the pair ignored, the kernel in (3.1) becomes 𝐾 (q, r) = 1 + cos 2q · r, so that the correlation 𝑅 in (3.2) becomes a Fourier transform of the source 𝑆. For the general task of imaging, in this work we reach for the Bayesian RL method that was originally developed for deblurring optical images, but has been by now invoked for nuclear problems bearing similarity to the optical deblurring. Here, we will carry that method to the application even farther from the method’s origins. 3.2 Deblurring method It is common for the differential results of nuclear measurements to be blurred due to the limited acceptance and resolution of the detectors. Hence, to be able to compare those measurements with the theory, we need good methodologies that incorporate such limitations. This served as our motivation to develop deblurring techniques with the primary goal of restoring the source function from the two-particle correlation function. In Chapter 4, we will demonstrate the use of the method in restoring the measured decay energy spectra from the 26O decay experiment. In this section, we introduce the terminologies and derive the deblurring method, but later in the chapter, we apply the method to determine the source function from d–𝛼 correlation function in heavy ion collision. 3.2.1 Blurring relation Let’s start with the blurring relation between a measured function, 𝑔, and the unknown function, G, can be stated as [56, 50, 51] ∫ 𝑔(𝑡′) = 𝑑𝑡 𝐴(𝑡′, 𝑡) G(𝑡) , (3.6) where, 𝑡′ is the measured variable, 𝑡 is the true variable, and 𝐴(𝑡′, 𝑡) is the conditional probability that the particles emitted at 𝑡 is measured at 𝑡′, also known as response function. The task at hand is to estimate the true function G, we achieve this using the RL algorithm that is discussed next. 34 3.2.2 Bayesian theory for RL algorithm In this section, we use the Bayes theorem to derive the RL algorithm, Eq. (3.12). For an ideal detector where all produced particles are perfectly measured, (cid:205) 𝑗 𝐴𝑖 𝑗 = 1, implies (cid:205) 𝑗 𝑔 𝑗 = (cid:205) 𝑗 G𝑗 = 𝑔𝑖 𝑁 and the probability of G𝑖 to occur is 𝑝(G𝑖) ≈ G𝑖 𝑁 , 𝑁, so the probability of measuring 𝑔𝑖, 𝑝(𝑔𝑖) ≈ Hence the conditional probability 𝑃(G𝑘 |𝑔𝑖) of G𝑘 to occur given 𝑔𝑖 is given by Bayes theorem 𝑃(G𝑘 |𝑔𝑖) = 𝑃(𝑔𝑘 |G𝑖)𝑃(G𝑖) (cid:205) 𝑗 𝑃(𝑔𝑘 |G𝑗 )𝑃(G𝑗 ) , where 𝑃(𝑔𝑘 |G𝑖) is a complement of 𝑃(G𝑘 |𝑔𝑖). We can then write 𝑃(G𝑘 ) as; 𝑃(G𝑖) = = ∑︁ 𝑘 ∑︁ 𝑘 𝑃(G𝑘 |𝑔𝑖) 𝑝(𝑔𝑘 ), 𝑃(𝑔𝑘 |G𝑖)𝑃(G𝑖) 𝑝(𝑔𝑘 ) (cid:205) 𝑗 𝑃(𝑔𝑘 |G𝑗 )𝑃(G𝑗 ) , Using 𝑃(G𝑖) ≈ G𝑖 𝑁 , 𝑃(𝑔𝑘 |G𝑖) = 𝐴𝑘𝑖 and (cid:205) 𝑗 𝑔 𝑗 = (cid:205) 𝑗 G𝑗 𝐴𝑘𝑖𝑔𝑘 (cid:205) 𝑗 𝐴𝑘 𝑗 𝑔 𝑗 G𝑖 = G𝑖 ∑︁ . 𝑘 (3.7) (3.8) (3.9) (3.10) Therefore, we can solve G, iteratively as discussed next. 3.2.3 Deblurring relation As mentioned before, the deblurring approach is based on the RL algorithm. In the optical blurring problem, a photon is measured with a property 𝑡′, while its true property is 𝑡. When the properties are discretized, such as in attributing the photon to a particular pixel, the relation in (3.6) becomes one in the matrix form between the distribution vectors: ∑︁ 𝑔𝑖 = 𝐴𝑖 𝑗 G𝑗 . 𝑖 (3.11) A deblurring method, such as RL, seeks to determine the distribution G, when knowing 𝑔 and 𝐴. To arrive at the RL strategy, a backward relation between 𝑔 and G is invoked, that involves a conditional probability 𝑄 that is complementary to 𝐴. Requiring the fulfillment of a Bayesian relation involving 𝐴 and 𝑄, G is searched through iterations G (𝔫+1) 𝑗 = G (𝔫) 𝑗 𝐴 𝑗𝑖 (cid:205)𝑖 𝛼𝑖 𝑔𝑖 𝑔 (𝔫) 𝑖 (cid:205)𝑖 𝛼𝑖 𝐴 𝑗𝑖 35 ≡ G (𝔫) 𝑗 𝑀 (𝔫) 𝑗 . (3.12) Here, 𝔫 is the iteration index, 𝑀 (𝔫) is an amplification factor, and 𝑔(𝔫) is prediction for the observation at 𝔫th iteration: 𝑔(𝔫) 𝑖 = ∑︁ 𝑖 𝐴𝑖 𝑗 G (𝔫) 𝑗 . (3.13) If we compare Eqs. (3.1) or (3.3) to (3.6), we can see a connection in the analogous mathematical structure. Moreover, each of the quantities in (3.1) has a probabilistic interpretation, though only 𝑆 ties directly to G in (3.6). We will primarily rely on the analogous mathematical structure in (3.1) or (3.3) and (3.6) and attempt to use the R-L method to deduce 𝑆. The weights 𝛼 in (3.12) can serve to focus attention on the region of relative momenta in the correlation function dominated by the interplay of the particles with each other. The iteration, 𝑛, in Eq. (3.12) stops when the ratio, 𝑔𝑖/𝑔(𝑛) 𝑖 , approaches one (see the inserted figure in Fig. 3.2). The figure illustrates the restoration of the Gaussian source function; the stars are the original Gaussian source function, used to construct d−𝛼 correlation function displayed in Fig. 3.1 right panel. The figure was obtained by establishing an analogous relationship between the blur function, Eq. (3.6), and discretized correlation, discussed in the next section. In Eq. (3.6), the quantity 𝑔 is analogous to the correlation function (𝐶 (𝑞)), 𝐴 to the kernel matrix (𝐾 (𝑟, 𝑞)), and G to the source function (𝑆). The dashed lines show the restored source function, we show the restoration for different numbers of iterations (n), where it is clear that when 𝑛 is small, 𝑛 =100 the restoration is not good and it improved as you increase 𝑛. In this example, convergence (or better restoration) is achieved at 𝑛 =50000. It is important to note that in the case of noisy distributions (i.e., Section 3.3 and Chapter 4) it is worth to stop the iteration as soon as a good restoration is observed because increasing the number of iterations beyond that may lead to instabilities in the restoration In our discussion, we consider the RL algorithm in one dimension case, but the multi-dimension case should be also straightforward and will be in consideration for our future works. 3.2.4 Regularization and error consideration In a typical experiment, the quantities 𝑔𝑖 in Eq. (3.11) represent measured values, such as decay energy spectra measurements (see Ref. [50, 57]). These measurements are often affected 36 Figure 3.2 The figure display the restoration of the source function from d−𝛼 correlation calculated assuming Gaussian using RL method. RL algorithm employed d–𝛼 correlation and Kernel matrix constructed from d–𝛼 squared wave function shown respectively on the right and left panel in Fig. 3.1. Different dashed lines correspond to the restoration for n number of iterations. The insert represents the ratio, 𝐴𝑛 = 𝑔/𝑔𝑛, as it is seen in Eq. (3.12). The algorithm converges when 𝐴𝑛 approaches one (see red dashed). The stars show the original source function. by finite statistics due to acceptance and resolution effects of the detectors, leading to fluctuations. In Ref. [50] as well as in the next chapter, we considered the counts in bins of measured energy spectra to fluctuate following the Poisson probability distribution function. This implies that the error of the bin’s counts is given as the square root of the bin content. However, in the case of correlations, we construct a noisy correlation by adding a relative momentum-dependent Gaussian noise of width, 𝜎. In this case, the uncertainty in the inference of G may be assessed by resampling 𝑔 for restoration, within the errors of the measurement. This can be done by constructing a set of alternative measurements that occur simultaneously, consistent with the best mean and variance, 37 through sampling from a normal probability distribution of 𝑔∗ 𝑖 , N (𝑔∗ 𝑖 |𝑔𝑖, 𝜎2) ∝ 𝑒− (𝑔∗ 𝑖 −𝑔𝑖 )2 2𝜎2 . (3.14) Then we perform RL restoration (Eqs. (3.12)-(3.13)) by replacing 𝑔𝑖 with 𝑔∗ 𝑖 to obtain G∗. However, estimating errors in the restored distribution is very challenging because the statistical error model in this distribution is not well known. In Chapter 4, we discuss uncertainty quantifi- cation in the restored decay energy distribution. But the RL algorithm suffers high-frequency instability due to noise amplification after a few iterations. These instabilities are suppressed by applying regularization in the algorithm. The first level of regularization is binning in the measured distributions and the response matrix. If the bin size is small compared to the resolution, it will result in oscillations in the restored samples. On the other hand, if the binning size is large, the resulting restored sample will be smooth. The other level of regularization is applied during the iteration process and works as follows: when instability appears, the regularization reduces it by increasing or decreasing the height of the oscillation. The regularization remains in effect until the number of iterations 𝑛 when the scaling ratio 𝑔𝑖/𝑔(𝑛) approaches unity. We discuss regularization 𝑖 in more detail in Chapter 4 (also see Refs. [50, 51]). 3.3 Deblurring to restore source function from two-particle correlation As illustration in this chapter, we choose correlations between deuteron and alpha particles. For this particle combination, scattering phase shifts have been measured and phenomenological potentials were developed allowing for calculations of scattering wave functions. Also correlation functions between those particles have been measured. Besides feasibility of the source inference with a deblurring algorithm, we will consider practicalities of the inference, such as the binning decisions for 𝐶 and 𝑆, errors of inference and impact of detector resolution. Typical correlation function from measurements is shown in the panel (a) of Fig. 3.3. The binning in relative momentum and scatter of points is quite typical for the measured light particle correlations in heavy-ion collisions. The pronounced resonance peak at 𝑞 ∼ 40 MeV/𝑐 represents the formation and decay of a 𝑑-wave resonance in 6Li and illustrates interest in the correlations as a 38 Figure 3.3 (a) Deuteron-alpha correlation vs magnitude of relative momentum q = 𝜇(v1 − v2) in the center of mass of the pair. Points represent the correlation measured in the 𝐸/𝐴 = 44 MeV 40Ar+27Al reaction by Ghetti et al. [8]. Line and shaded regions represent results from RL source restoration. (b) Source inferred from the measured correlation function in (a). The solid line and shaded areas represent results from RL source restoration. For comparison, a Gaussian source function with radius 𝑅0 = 4.5 fm is shown too. source of information about interactions in different channels. We use the features of the particular measurement as a general guidance in testing the capabilities of the deblurring algorithm in source restoration. Ahead of the restoration from data, we carry out tests where we first apply a forward relation between assumed source and correlation function. For the sake of deblurring, the source is discretized. To fix attention we take the source distance range of (0-50) fm and divide it into 𝑀 even bins and represent an isotropic 𝑆 in the form 𝑆(𝑟) = 𝑁 ∑︁ 𝑗=1 𝑆 𝑗 𝑔 𝑗 (𝑟) , (3.15) 39 50100150q [MeV/c] 24C(q)(a)Ghetti et. al., (2006)RL restoration 01020304050r [fm] 0.00000.00010.00020.00030.0004S(r) [fm3](b)Gaussian (R0=4.5 fm) Restored source (RL)2 uncertainty1 uncertainty where 𝑔 𝑗 is a characteristic function for the 𝑗’th bin, 𝑔 𝑗 (𝑟) = 1 , 0 ,    if 𝑟 𝑗−1/2 < 𝑟 < 𝑟 𝑗+1/2 , otherwise. (3.16) We compute the wave functions for the relation (3.1) of the source with the correlation using the interaction potentials developed by McIntyre and Haeberli to fit the 𝑑-𝛼 phase shifts from measurements. With the correlation function determined at momenta 𝑞𝑖, 𝑖 = 1, . . . , 𝑁, the mapping of the correlation onto the blurring problem amounts to 𝐶 (𝑞𝑖) ≡ 𝐶𝑖 ↔ 𝑔𝑖, 𝑆(𝑟 𝑗 ) ≡ 𝑆 𝑗 ↔ G𝑗 and 4𝜋 ∫ 𝑟 𝑗+1/2 𝑟 𝑗 −1/2 𝑑𝑟 𝑟 2 |Ψ(−) q𝑖 (r)|2 ↔ 𝐴𝑖 𝑗 . We generally use fewer points for the source than in the correlation function, 𝑀 ≤ 𝑁. For the selection of 𝑞 in Fig. 3.3(a) (𝑁 = 34), our source resolution ends up at 𝑟max/𝑀 ∼ 4 fm (𝑀 ∼ 14 for 𝑟max = 56 fm). Details on that will be provided later. Given the relative success of Gaussian sources in correlation analyses and the features of the particular measured function, we take the source for testing our restoration in the form 𝑆𝐺 (r) ∝ exp (− r2 2𝑅2 0 ). That source with 𝑅0 = 4.5 fm and normalized to 1 is shown with points in the panel (b) of Fig. 3.3. The correlation function generated with the forward source-correlation relation (3.3) is shown with points in the panel (a) of Fig. 3.4. Figure 3.4 Test of source restoration from a correlation function without noise. Panel (a) displays the discretized 𝑑-𝛼 correlation function (points) generated from the discretized Gaussian source function displayed in panel (b) (points). Panel (b) displays further the source restored with RL algorithm (solid line) from the discretized correlation function in (a). Finally, as a cross-check, panel (a) displays the correlation function (solid line) produced from the restored source. 40 Gaussian no noiseRL restoration We test restoration with the RL method both for a smooth and noisy input correlation functions 𝐶 (𝑞). The source restored from the smooth function in the panel (a) of Fig. 3.4 is illustrated with lines in the panel (b). It may be observed that the input and restored source cannot be distinguished within the resolution of the figure. To test the case of a noisy 𝐶, see [50], we add model fluctuations to the smooth 𝐶. Figure 3.5 Test of source restoration from a correlation function with noise. Points in panel (a) represent the correlation function from the model 𝑅 = 5.5 fm Gaussian source, with a sample Gaussian noise added. The Gaussian source itself is represented with a dashed line in panel (b). The shaded areas and the solid line in panel (b) represent results of restoration from an ensemble of correlation functions such as in (a) where the Gaussian noise was repeatedly sampled. The solid line shows the average restored source and the darker and lighter shaded areas show the extent of the 1-𝜎 and 2-𝜎 range in the distribution of restored source function values. Finally, the solid line and the darker and lighter shaded regions in panel (a) show similar information for the correlation functions from the ensemble of restored source functions. The narrow width of the 1-𝜎 region in (a), compared to the original assumed uncertainties in 𝐶, stems from the constraint of 𝑆 being nonnegative, built into the restoration, and 𝐶 dependent effectively on values of 𝑆 at just few points in 𝑟. Specifically, we observe that we approximate the scatter of points in the experimental 𝑑-𝛼 correlation function in Fig. 3.3(a) around a smooth function 𝐶sm with a Gaussian characterized by a 𝑞-dependent width approximately equal to 0.15√︁𝐶sm(𝑞). With this we sample noisy cor- √︁𝐶𝐺 (𝑞) N (0, 1), where 𝐶𝐺 (𝑞) is the relation functions for our tests from 𝐶 (𝑞) ∼ 𝐶𝐺 (𝑞) + 0.15 smooth function generated with the Gaussian source. A sample correlation function with noise is illustrated with stars in Fig. 3.5(a). We generate an ensemble of such correlation functions and the corresponding ensemble of restored sources. The average values of the restored sources at different 𝑟 are illustrated with a solid line in Fig. 3.5(b). The extent of the 1-𝜎 and 2-𝜎 ranges in the value distributions for restored sources at different 𝑟 are illustrated as dark and light shaded 41 50100150q [MeV/c] 024C(q)(a)Gaussian with noiseRL restoration 02040r [fm] 0.00000.00010.00020.00030.0004S(r) [fm3](b)Restored source (RL)Original source 2 uncertainty1 uncertainty areas, respectively. It can be observed that the restored values generally agree within 1-𝜎 with the original Gaussian source. It is important to note that the RL algorithm can suffer from noise amplification after a modest number of iterations, see Refs. [51, 50] and references within. In our calculation, we suppressed potential instability by applying a regularization in the algorithm. The first level of regularization is the binning choice in the source function (see Eq. (3.15)); too many bins lead to oscillations in restoration, and too few bins lead to the loss of information, and we discuss binning choice in detail later in the section. The second level of regularization, we use here, is the one developed in Ref. [51], where the parameter 𝜆 = 0.015 was chosen. The main goal of this section remains the restoration of a source from data following the RL algorithm. The 𝑑-𝛼 pairs yielding the correlation function in Fig. 3.3 have been measured by Ghetti et al. [8] at forward angles 0.7◦ < 𝜃 < 7◦ in 44 MeV/nucl 40Ar+27Al collisions. When narrow structures are measured in an experiment, such as the 𝑞 ∼ 40 MeV/𝑐 peak in the correlation function, then detector resolution needs to be considered. The impact of the resolution, as far as the forward relation between the source and correlation is concerned, is in the modification of the kernel in Eq. (3.1), where the original kernel 𝐾 (q, r) gets convoluted with an appropriate detector resolution function pertaining to q. Relative to the measured vector, the resolution can modify the magnitude of the vector q in the wave function, as well as its direction, especially for low 𝑞. However, for low values of the 𝑞𝑟 product, only low ℓ will matter in the wavefunction squared in the kernel, so the sensitivity to the q direction will be weak. On the other hand, in the presence of resonances the sensitivity to the pair c.m. energy, tied to the detector energy resolution, can be quite strong. In Ref. [42] it has been proposed to account for the smearing in 𝑞 by folding the original kernel with a Gaussian in 𝑞 of width 𝜎𝑞 adjusted to the energy resolution. For an angle-averaged correlation function, this yields 1 2𝜋𝜎𝑞 in Eq. (3.3). In the limit of low relative velocity 𝑣 for the pair, as compared to the velocity of the q′ (r)|2 , 𝐾 (𝑞, 𝑟) = |Ψ(−) (3.17) 𝑑𝑞′ √ − (𝑞−𝑞′ )2 2𝜎2𝑞 e ∫ pair c.m. 𝑉, simple kinematic considerations yield a relation between the resolution in energy 𝜎𝐸 42 𝜎𝐸 𝐸 𝑉. With the energy resolution in and that in relative velocity 𝜎𝑣, for angle-averaged q: 𝜎𝑣 = 1√ 6 𝜎𝐸 𝐸 ∼ 2% [58] and 𝑉 corresponding to the projectile-like fragments [8], the particular experiment we get 𝜎𝑞 = 𝜇 𝜎𝑣 ≃ 3 MeV/𝑐. An alternative to estimating 𝜎𝑞 using previously inspected 𝜎𝐸 /𝐸 is to assess 𝜎𝑞 directly from the correlation measurement when a narrow resonant state is present such as for 𝑑-𝛼. In the source inference from data it is also necessary to decide on the source discretization, i.e., 𝑟max and 𝑀 in our scheme. For this, we compare the correlation 𝐶 𝑅𝐿 from the source 𝑆𝑅𝐿 inferred through RL deblurring to the measured correlation 𝐶 and construct 𝜒2: 𝜒2 = 𝑁 ∑︁ 𝑖=1 (cid:16) 𝐶 𝑅𝐿 𝑖 − 𝐶𝑖 𝜖𝑖 (cid:17) 2 . (3.18) Here, the uncertainties are estimated as 𝜖𝑖 ≈ 0.15 𝐶𝑖. At fixed 𝑟max, 𝜎𝑞 and 𝑀, we carry out the √ RL deblurring and then minimize 𝜒2 under variation of 𝜎𝑞 and 𝑀. With the number of degrees of freedom (DOF) calculated as 𝑁–𝑀–2, where 𝑀 is for the number of source bins, and 2 is for 𝑀 and 𝜎𝑞 adjustments, we show in Fig. 3.6(a) 𝜒2/DOF obtained in this manner as a function of 𝑟max. We find a flat behavior of 𝜒2/DOF at 𝑟max ≳ 50 fm, close to 1 for the particular parametrization of 𝜖. In Fig. 3.6(b) we show a contour plot of 𝜒2/DOF at fixed 𝑟max = 56 fm when 𝜎𝑞 and 𝑀 are varied. We typically find the minimum at 𝜎𝑞 ∼ 3.8 MeV/𝑐 and 𝑟max/𝑀 ∼ 4 fm at different 𝑟max. In the source restoration illustrated in Fig. 3.3, we use 𝑟max = 56 fm, 𝑀 = 14 and 𝜎𝑞 = 3.7 MeV/𝑐. The proximity of 𝜎𝑞 from the fit to that from the resolution estimate should be noted. Compared to the Gaussian source there, that restored approaches faster low values by 𝑟 ∼ 10 fm, but then it has higher values above ∼ 15 fm. Notably, when an unnormalized Gaussian source parameterization is used to describe correlations it gets combined with Eq. 3.2 or angle-averaged version thereof. In the latter case, it is assumed that the strength completing source normalization is located at large 𝑟. 3.4 Summary In this chapter, we have demonstrated the use of a deblurring method, successful in optics, to infer the emission source from low relative-velocity correlation function. As the example, we chose the 𝑑-𝛼 correlation that features a narrow and overlapping broad resonances, as well as Coulomb 43 Figure 3.6 𝜒2 per degree-of-freedom results in describing the measured 𝑑-𝛼 correlation function in terms of the correlation from an RL restored source. Panel (a) shows minimal 𝜒2/DOF for 𝜎𝑞 and 𝑀 optimized, while 𝑟max is set. Panel (b) shows a contour plot of 𝜒2/DOF when 𝑟max = 56 fm and 𝜎𝑣 and 𝑀 are varied. The circle marks the minimum at 𝜎𝑞 = 3.7 MeV/𝑐 and 𝑀 = 14. depletion at low 𝑞. The source inference involves determination of relative wave function in order to generate the kernel for the KP relation. In parallel to the binning of the correlation function typical for experiment, the source and, correspondingly, the kernel gets discretized, yielding a transfer matrix. Impact of detector resolution can be accounted for in the matrix, in parallel to the physics connecting the source and correlation function. The source restoration progresses through RL iterations until source stabilization. Uncertainties in source determination can be assessed by resampling the experimental correlation function with experimental uncertainties. We have tested the source restoration from a 𝑑-𝛼 correlation function with a Gaussian source, both for an idealized function without uncertainties and with uncertainties. In both cases, an appli- cation of the KP relation followed by the RL deblurring returned a source information consistent with the input. In analyzing the measured 𝑑-𝛼 correlation, we demonstrated that, for sharp reso- nances, the impact of detector resolution may be read off from the correlation itself. The source restored from the data through RL deblurring is close, within restoration resolution, to a Gaussian source in the central part, but it first approaches low values more abruptly, to then exhibit a tail that the Gaussian source lacks. We hope that the RL or other optical deblurring algorithms, applied as here, may turn out being 44 2040M123456q [MeV/c](b)1234567892/DOF30405060rmax [fm]1.11.21.31.41.52/DOF(a) useful in inferring the sources from correlation 45 CHAPTER 4 DEBLURRING EXPERIMENTAL DECAY ENERGY SPECTRA: THE 26O CASE The deblurring method discussed in Chapter 3 is applied in this chapter to restore the decay energy spectrum for the three-body decay of 26O, as measured by the MoNA Collaboration [59]. In the next chapter, we will discuss a deep neural network classifier method to identify number of states in the measured decay energy, both method test and complete each other. The material discussed in this chapter was previously published in Ref. [50]. The first two sections provide an introductory description of the MoNa Collaboration experi- ments and an overview of the deblurring method. In Section 4.3, we discuss the practicalities of deblurring, specifically how the procedure is expanded to handle noisy data. The methodology is then tested on simulated data in Section 4.4. Moving on to Section 4.5, we delve into the deblurring of the measured 26O decay energy spectrum. Finally, in Section 4.6, we present our conclusions and outlook. 4.1 Introduction Invariant mass spectroscopy allows experimental access to unbound states. However, inter- preting and extracting physics from the measured decay energy spectra are often challenged by limited resolution and distortions caused by experimental acceptance effects. This is particularly true in investigations of neutron-unbound states, since they involve the measurement of neutrons and charged decay fragments in coincidence. In a decay experiment of this type, the neutron- unbound state is populated through a nuclear reaction induced by a rare isotope beam, typically proton-removal. The unbound state decays immediately, and by measuring the momentum vectors of the decay products, the invariant mass of the unbound system can be calculated. The measured decay energy spectrum can then be reconstructed by subtracting the masses of all constituents of the system. The squared invariant mass of the initial system, 𝑀 is given as [59]: 𝑀 2 = 𝑀 2 𝐴 + 𝑛 ∑︁ 𝑖 𝑚2 𝑁 + 2 𝑛 ∑︁ 𝑛+1 ∑︁ 𝑖 𝑗+1 𝐸𝑖𝐸 𝑗 − Pi.Pj, (4.1) 46 where 𝐸 and P are energies and momenta of the daughter particles, and 𝑀𝐴 and 𝑚𝑁 respectively represent the rest mass of daughters, fragments and neutron. The total 𝑛 + 1 terms correspond to the n neutrons plus fragment 𝑀𝐴 and 𝑖 and 𝑗 in the double sum index over all of the daughter particles. The decay energy, 𝐸𝑑 is defined as the difference between the invariant masses of the initial nucleus and decay products 𝐸𝑑 = 𝑀 − 𝑀𝐴 − 𝑛 ∑︁ 𝑖 𝑚 (𝑖) 𝑁 . (4.2) In this work, we are focusing on the two-neutron emission decay energy spectrum of 26O. This unbound nucleus was recently measured by the MoNA Collaboration [59], with the setup illustrated in Fig. 4.1. The exploration typifies efforts to learn about the structure of nuclei towards the neutron-drip line [60]. In general, measuring neutron momenta implies the use of a neutron detector array that usually has limited position resolution and detection efficiency. Similarly, measuring the momenta of charged particles involves tracking the particle trajectories back through a magnetic field and part of the reaction target to determine the angle and energy at the point of the breakup reaction, which is not accessible to direct measurements. The procedures introduce variations and uncertainties in such a way that the measured decay energy distribution is only a distorted and blurred image of the true decay energy spectrum of the unbound system. In the present work, we will utilize two novel methods of inferring features of the true decay spectrum: a deblurring algorithm and a deep neural network approach. There is much potential for these strategies outside of the particular problem. 4.2 Interplay of the Experiment and the Analysis Methods 4.2.1 Experiment and Construction of Transfer Matrix In the experiment considered here, the strongest distortion of the spectrum stems from the acceptance and resolution effects of the Modular Neutron Array and Large multi-Institutional Scintillator Array (MoNA-LISA), and the fact that it is hard to detect neutrons with good efficiency and determine their location with good precision. Detailed simulations of the detector setup allow to quantify the impact of the detection process on a decay energy spectrum and cast it in the form 47 Figure 4.1 (Color online) The MoNA experimental setup for invariant mass measurements in search of neutron-unbound states includes the Sweeper magnet, charged particle detector suite, and neutron detector array. The rare isotope beam (orange arrow) impinges on a reaction target where the unbound state is populated in a nuclear reaction. The charged breakup fragments (red shaded area) are directed by a magnetic dipole field into the charged particle detector suite, while the neutrons (green shaded area) travel along the beam direction to the neutron detector array. of a response matrix or transfer matrix, cf. Figs. 4.2 and 4.3. The matrix folded with any input decay spectrum and no detector distortions produces the spectrum expected to be measured in the experiment with those distortions imposed. In constructing the matrix, decays are simulated by randomly drawing the decay energy, 𝐸𝑑, from a uniform distribution and randomly selecting the orientation of the decay event in the 26O frame. Each decay is processed through a simulation of the detector response. The decay energy spectrum is then constructed from that response in the same fashion as for the measured data. By selecting a narrow range of input decay energies, the resulting ‘resolution-folded’ spectrum, 𝐸′ 𝑑, for a given 𝐸𝑑 is produced (cf. Fig. 4.3). The 𝐸𝑑-values used as examples are shown as thick red lines in Fig. 4.3 and red arrows in Fig. 4.2. The full response/transfer matrix is built from the resolution-folded spectra 𝐸′ 𝑑. A difference in normalization for the matrix shown here compared to Fig. 4.3 should be noted. In Fig. 4.3, the normalization is for practical purposes to illustrate the difference in response for two example 𝐸𝑑, and in Fig. 4.2 it is appropriate for the matrix in continuum limit, representing conditional probability density. As a further note, integration over the measured energy yields the probability of the event at a given input energy getting accepted, ∫ 𝑑𝐸′ 𝑑 = 𝐸𝑑 diagonal is marked in the figure to guide the eye. The 𝑑 |𝐸𝑑) = 𝑃(𝐸𝑑). The 𝐸′ 𝑑 𝑃(𝐸′ 48 Figure 4.2 (Color online) The response matrix 𝑃(𝐸′ 𝑑 |𝐸𝑑) of the MoNA experimental setup depicted in Fig 4.1 used in measuring the decay energy of the three-particle decay 26O →24 O + 2𝑛. The diagonal line corresponds to 𝐸𝑑 = 𝐸′𝑑 and is plotted on the figure to guide the eye, and the arrows point to the values of 𝐸𝑑 that correspond to the observed resonance states of 26O. See text for details. rapid decrease in the probability at high 𝐸′ 𝑑 indicates that an event at high 𝐸′ 𝑑 has a low chance to get recorded. 4.2.2 Accessing resonance properties It is common practice to assess the original, undistorted decay energy spectra with parameter estimation techniques. For example, neutron-unbound resonances are often [61, 62, 63, 60] modeled using energy-dependent Breit-Wigner line shapes (See Eq. (4.11)) [64]. Parameter estimation methods, such as 𝜒2 minimization, are used to extract the resonance energy, width and angular momentum for each resonance state. For the remainder of this chapter we refer to such methods as 49 0.00.51.01.52.02.53.03.54.0Ed (MeV)0.00.51.01.52.02.53.03.54.0E0d (MeV)0.00010.00020.00030.00040.00050.00060.00070.0008P(E0d|Ed) (MeV1) Figure 4.3 (Color online) Construction of individual columns in the response/transfer matrix (TM). A single bin (0.2 MeV width) in input energy 𝐸𝑑 is uniformly populated with events, as illustrated by the solid (red) histograms. Processing of the events through a simulation of the detection system yields corresponding event partition across bins in 𝐸′ 𝑑 illustrated by the open (green) histograms. traditional fit methods. The traditional approaches require decisions on the number of parameters to fit for the original spectrum, such as the choice of the number of resonances present in the explored energy range. The proposed deblurring method aims at restoring the features of the original spectrum without assuming how many states it contains. We complement the results from the novel approach by carrying out the standard chi-square minimization and assuming different numbers of resonance peaks in the data. 50 0.81.0(a)input at 500 KeVoutput from TM0246810Ed, E0d (MeV)0.000.050.10Normalized yield0.81.0(b)input at 2 MeVoutput from TM0246810Ed, E0d (MeV)0.000.050.10Normalized yield 4.2.2.1 The Richardson-Lucy deblurring procedure Our deblurring procedure employs the Richardson-Lucy (RL) algorithm initially developed to restore blurred images in optics [65, 66]. Over time the algorithm found use in astronomy [67] and medicine for medical images analysis [68], to list a few. In high-energy physics, analogous develop- ments progressed [69] without realization of the prior work elsewhere. Recently, Danielewicz and Kurata-Nishimura [51] have demonstrated that a nonlinear extension of the algorithm could be used to determine three-dimensional (3D) momentum distributions of products in intermediate-energy heavy-ion collisions. The RL algorithm derivation relies on the Bayes’ theorem and it follows an iterative procedure to find a self consistent solution. The algorithm only uses the distorted spec- trum and discretized response function of the apparatus, or Transfer Matrix (TM), as inputs. The spectrum entries and matrix elements are positive definite and carry probabilistic interpretation. The restoration of the original spectrum is an inverse problem, but it progresses in the deblurring without directly inverting the TM, an uncommon approach for inverse problems [70]. In maintain- ing the restored spectrum positive throughout the iteration procedure and by avoiding a direct TM inversion, serious singularity problems plaguing inverse problems are avoided. In Ref. [51], the RL was implemented without consideration of noise. In this thesis and the associated paper [50], we expand the utility of the algorithm by considering measurement statistics and improve on the assessment of what is actually learned from the data. However, in other fields, it has been demonstrated that the RL algorithm suffers from short-wavelength instability due to noise amplification after a limited number of iterations [71, 72, 73]. To overcome this challenge, we introduce a regularization in the algorithm that tames the short wavelength component in the deblurring solution. There are several options for such regularization, the Gaussian function smoothing being one such example. The smoothing requires considerations of a function width and boundary conditions [71]. Another regularization option is the use of denoising algorithms that invoke nonlinear combinations of derivatives of restored spectra [72, 74], commonly termed Total Variation (TV). We use a simple version of TV regularization employed in Ref. [51], but we make its strength 51 increase with energy, as the impact of noise on a restored spectrum increases at higher energy. With this approach we are able to arrive at stable deblurring solutions after just few RL iterations. 4.3 Richardson-Lucy Deblurring Algorithm 4.3.1 Setting In a nuclear decay experiment, the decaying nucleus, characterized by a total four momentum p = (𝐸, 𝑝𝑥, 𝑝𝑦, 𝑝𝑧), can be thought of as an emitter of particles that fly off towards the detector. The detector records the particles with some efficiency and allows to determine their four-momenta with some accuracy. From the combination of those four-momenta, the invariant mass of the decaying nucleus is determined, 𝑀 = √︁p2, and, over many events, the particle decay energy spectrum is established [59, 75]. Structures in that spectrum can tell us about the resonance states of the decaying nucleus. Limitation in the detector resolution makes the measured spectrum 𝑓 blurred compared to the true decay spectrum F of the nucleus. The blurring relation between 𝑓 and F can be written as ∫ 𝑓 (𝐸′ 𝑑) = 𝑑𝐸𝑑 𝑃(𝐸′ 𝑑 |𝐸𝑑) F (𝐸𝑑) . (4.3) Here, 𝐸′ 𝑑 is the measured energy, 𝐸𝑑 is the true energy and 𝑃(𝐸′ 𝑑 |𝐸𝑑) is the conditional probability that products for a nucleus decaying at 𝐸𝑑 are registered, the event is accepted and determined to represent the decay energy 𝐸′ 𝑑. In the context of an experiment, 𝑃(𝐸′ 𝑑 |𝐸𝑑) represents the response function of the apparatus, but in the context of blurring analyses it may be called a blurring or transfer function. As an extreme example, 𝑃(𝐸′ 𝑑 |𝐸𝑑) = 𝛿(𝐸′ 𝑑 − 𝐸𝑑) represents an ideal detector. Eq. (4.3) invokes the spectra 𝑓 and F in the limit of infinite measurement statistics. In practice, the spectra get discretized, most often simply binned. Moreover, in an experiment, 𝑓 only gets determined with some accuracy, and even 𝑃 gets established with some resolution. Under discretization, the blurring relation (4.3) acquires the matrix form ∑︁ 𝑓𝑖 = 𝑃𝑖 𝑗 F𝑗 , 𝑗 (4.4) where 1 ≤ 𝑖 ≤ 𝑁, 1 ≤ 𝑗 ≤ 𝑀 and 𝑃𝑖 𝑗 represents the conditional probability density integrated over a discretization form factor (typically Δ𝐸 bin) in 𝐸𝑑 and averaged over one in 𝐸′ 𝑑. As such, 52 the matrix elements 𝑃𝑖 𝑗 are positive and 𝑃𝑖 = (cid:205) 𝑗 𝑃 𝑗𝑖 represents probability than an event at decay energy near 𝐸𝑖 𝑑 is analysed. In our analysis of decay-energy spectra, we most often employ Δ𝐸𝑑 = 0.2 MeV binning. To construct the transfer matrix (TM), 𝑃(𝐸′ 𝑑 |𝐸𝑑), for the three-particle decay 26𝑂 →24 𝑂 + 2𝑛 experiment [59], we randomly draw the decay energy, 𝐸𝑑, from a uniform distribution and draw the orientation of the decay event in the frame of 26O. Each decay is then processed through the simulated response [59, 75] of the detector setup schematically illustrated in Fig. 4.1. The outcomes are sorted by bins in 𝐸𝑑 and 𝐸′ 𝑑 and their counts per 𝐸𝑑 bin entry become TM elements. The constructed matrix is illustrated in Fig. 4.2. The TM construction is additionally illustrated in Fig. 4.3 for individual 𝐸𝑑 bins. An 𝐸𝑑 bin is uniformly populated with events, as indicated by the solid (red) histograms shown in Fig. 4.3. Those events are processed through the simulation of the detector response and sorted according to 𝐸′ 𝑑 bins, as indicated by the open (green) histograms. After renormalization, the open (green) distributions in Fig. 4.3 become columns in the TM normalized as probability density 𝑃(𝐸′ 𝑑 |𝐸𝑑), or as contributions to the probability 𝑃𝑖 𝑗 in practical calculations with discretized spectra. 4.3.2 Deblurring The goal of deblurring is to estimate F when only 𝑓 and 𝑃 are known. The Richardson-Lucy (RL) algorithm [65, 66, 69, 51] relies on the conditional probability 𝑄(𝐸𝑑 |𝐸′ 𝑑) complimentary to 𝑑 |𝐸𝑑). Bayesian theorem linking the two probability densities yields a set of equations [51] 𝑃(𝐸′ that can be solved for F (𝐸𝑑) by iteration: 𝑓 (𝑛) 𝑗 = ∑︁ 𝑃 𝑗𝑖 F (𝑛) 𝑖 , 𝑖 = F (𝑛) 𝑖 F (𝑛+1) 𝑖 ∑︁ 𝑗 𝑓 𝑗 𝑓 (𝑛) 𝑗 𝑃 𝑗𝑖 𝑃𝑖 . (4.5) (4.6) Here, 𝑛 is the iteration index. We have chosen to start RL iterations with a rough guess for F (0), such as scaled up 𝑓 . The iterations is stopped once F (𝑛) ceases to change with 𝑛. For distributions that quickly change with their arguments, such as 𝐸𝑑 here, the long-term convergence may be slow and for large 𝑛 53 numerical seesaw instabilities in the arguments may set in. That instability can be tamed with a renormalization factor [72, 51] 𝐼 (𝑛) applied to the r.h.s. of (4.6): 𝐼 (𝑛) = 1 1 − 𝜆 D · ∇ . (cid:17) (cid:16) ∇𝐹 (𝑛) |∇𝐹 (𝑛) | (4.7) Here, D is a vector with components that are intervals over which F is discretized in its arguments (bin sizes), the divergence is approximated in low order based on that discretization and 𝜆 is a small positive number. In a one-dimensional case, such as here, the factor becomes simply 𝐼 (𝑛) 𝑖 = 1 1−𝜆 , 1 1+𝜆 , if F (𝑛) 𝑖 < F (𝑛) 𝑖−1,𝑖+1 if F (𝑛) 𝑖 > F (𝑛) 𝑖−1,𝑖+1 , , 1 , otherwise .    (4.8) This factor suppresses any patterns of maximae and minimae emerging on the discretization scale. However, when wider-scale maximae or minimae arise, the factor will be impacting them too. As uncertainties in the restored F will be of interest here, the use of the above regulation factor will introduce a relative error of the order of 𝜆 around the extrema of the restored F . 4.3.3 Fluctuations and other practicalities The blurring relation (4.3) invokes spectra in the limit of infinite statistics. However, the spectra are measured at finite statistics and its characteristics are expected to fluctuate compared to those at high statistics. Let 𝑓 represent the average event numbers registered in different bins of decay energy for measurement series carried out over a specific measurement time. If we carry out repeated measurement series over that time, event counts for individual bins will fluctuate in a Poisson-like manner. If we carry out just one measurement series, then the event count in the bin 𝑖, 𝑓𝑖, is our best estimate for the mean count and the best estimate for the mean squared deviation from that mean over repeated series [76]. When assessing uncertainties in the restored spectra, we build up an ensemble of alternative measurement results over the same time, consistent with the best estimates of the mean values for 54 decay energy bins and dispersion, by sampling the Poisson probability distribution for content 𝑓 ∗ 𝑖 , P ( 𝑓 ∗ 𝑖 | 𝑓𝑖) = 𝑖 𝑒− 𝑓𝑖 𝑓 𝑓 ∗ 𝑖 𝑓 ∗ 𝑖 ! . (4.9) We then carry out the RL restoration, Eqs. (4.5)-(4.6), with 𝑓𝑖 replaced by 𝑓 ∗ 𝑖 , arriving at F ∗ and we study the distribution of the latter within the ensemble. The algorithm requires F (0) ≥ 0 to start. However, we have not seen any significant sensitivity of the results to the fine details of F (0). In practice, the important factor is the number of iterations, a few hundreds is sufficient in our case, and the smoothing factor (see Eq. (4.10)). Within the higher end of the decay energy window in which we operate, usually up to 10 MeV, the counts tend to be low, fluctuating with energy and these fluctuations tend to be amplified in the restoration. Correspondingly, we make the parameter 𝜆 in the factor 𝐼, Eq. (4.8), increase with energy: 𝜆 = 𝜆0 (cid:16) 1 + (cid:17) 2(cid:17) , (cid:16) 𝐸 𝐸0 (4.10) and we typically use 𝜆0 = 0.035 and 𝐸0 = 6 MeV. The form and parameter values have been adjusted through experimentation. Notably, 𝜆 increases the bin-to-bin correlation, which is illustrated in Fig. 4.7. The values for 𝜆0 and 𝐸0 are chosen to reduce the noise oscillations in the restoration; this could also be done by increasing bin sizes. We choose 𝜆 to depend on E to suppress the oscillations in the restoration in the high E range, which is due to the finite statistics. 4.4 Tests of deblurring algorithm In this section we carry out tests of our deblurring procedures when applied to simulated data. We first consider data with negligible errors and then data with statistical errors comparable to those for the investigated decay energy measurements [59, 75]. Following physical expectations regarding the forms of input decay-energy spectrum, the spec- trum F (𝐸𝑑) in the tests is modeled as a superposition of Breit-Wigner distributions: F (𝐸𝑑) ≈ ∑︁ 𝐴𝑖 𝑖 0.5 Γ𝑖 (𝐸𝑑 − 𝐸𝑖)2 + (0.5 Γ𝑖)2 . (4.11) 55 Figure 4.4 (Color online) Restoration of decay-energy spectrum in the absence of noise. The dots (green) represent the original event distribution modelled with Eq. (4.11). Three wide peaks were assumed for the spectrum. The dashed (blue) line represents the blurred distribution, at adjusted normalization, and it has been obtained by folding the original distribution with the TM, cf. Eq. (4.4). The solid (red) line represents the distribution obtained by subjecting the blurred spectrum to deblurring with the RL algorithm, Eqs. (4.5) and (4.6). The restored and original distributions lie practically on top of each another. A binning in energy of 0.2 MeV was employed in generating these spectra. We are generally interested in the decay energy region extending up to 10 MeV, though we have also considered energies up to 14 MeV. Within such regions we have experimented with distributions containing (1–5) Breit-Wigner peaks at different energies 𝐸𝑖 and of different widths Γ𝑖 and amplitudes 𝐴𝑖. In the case we will use here for illustration, we take three peaks at 0.3, 2, and 4.5 MeV, with respective widths of 0.3, 0.85, and 1.3 MeV, see Fig. 4.4. For simplicity, we take 𝐴𝑖 ≡ 1. At first, we take the modeled input distribution F and multiply it by the TM to get 𝑓 . Up to some joint normalizing factor for both, these distributions stand for those in the limit of a very large statistics. The simulated input and measured distributions are illustrated in Fig. 4.4. To the simulated measured distribution we apply the RL algorithm, Eqs. (4.5) and (4.6). The restored distribution from the iteration is also shown in Fig. 4.4 and it lies practically on top of the original. This has been our typical finding for the limit of large statistics, no matter what input. In the restorations for large statistics, we usually can drop the smoothing factor (4.8). Next, we turn to simulations of ensembles of events, such as for real data. Specifically, we sample the shape of F within the energy range 𝐸𝑑 < 10 MeV, to get 𝐸𝑑 for a single event. Then we sample the probability density 𝑃(𝐸′ 𝑑 |𝐸𝑑) from TM to decide whether this event is accepted for 56 0246810Ed, E'd(MeV)01234Counts (104)restored spectrumoriginal spectrumblurred spectrum (× 400) Figure 4.5 (Color online) Restoration for four different examples, panel (a),(b), (c), and (d), of simulated data sets when Poisson noise is active. The sampled original spectrum is the same for each set and event statistics behind each set is similar to that believed to be behind the real data analysed in this work. The points represent the individually sampled sets with counts scaled up by a factor of 400. The dark blue and light blue bands illustrate the 𝜎 and 2𝜎 uncertainties resulting from spectra restoration with error sampling. In each panel, the dashed blue line represents the mean in the restoration ensemble for the set. The original spectrum (solid orange line) has three resonance peaks located at 0.3, 2, and 4.5 MeV with respective widths of 0.3, 0.85, and 1.3 MeV. Generally, we succeed in restoring the structures in the original spectrum using the RL algorithm. Binning of 0.2 MeV was used for the processed spectra. analysis and what the measured 𝐸′ 𝑑 is. We repeat the process until the number of analysed events is similar to that in the experiment. The needed number of input events provides a normalization for F . In Fig. 4.5, we show results from such four separate data simulations. Both the simulated measured 𝑓 (𝐸′ 𝑑) and underlying F (𝐸𝑑) are shown there. A measurement carried out over a specific beam time, with finite statistics, can be viewed as a member of an ensemble of measurements ran over the same time. We next attempt to simulate such an ensemble using only information in an individual generated data set, following the Poisson 57 012345Counts (104) (a)(a)(a)(a)(b)(b)(b)(b)restored spectrumoriginal spectrumsimulated data set (×400)95% uncertainty68% uncertainty0246810Ed, E'd(MeV)012345Counts (104)(c)(c)(c)(c)0246810Ed, E'd(MeV)(d)(d)(d)(d) distribution sampling discussed earlier, Eq. (4.9), to get 𝑓 ∗(𝐸′ 𝑑). To the individual 𝑓 ∗, we apply the RL deblurring algorithm to get an estimate of F . With this, we arrive at an ensemble of restored F that reflects uncertainties inherent in 𝑓 , within the methodology we adopt. In Fig. 4.5, we further show the characteristics of the ensemble of restored F , for each simulated data set, specifically the average values for the bins and 68% and 95% uncertainty ranges. It can be observed that the distributions of the restored values are generally consistent with the input F . Even if the dataset is noisy, the three peaks emerge in each case. Figure 4.6 (Color online) Outcome of averaging over restored distributions from 24 such simulations as in Fig. 4.5. The overall mean (dots) and the original spectrum (solid) are largely on top of each other. We complement the above resampling results by showing in Fig. 4.6 a distribution of restored F resulting from averaging over the distributions of restored F from a number of individual data simulations such as in Fig. 4.5. It can be seen that the average over a large number of ensembles 58 begins to approach the input F suggesting a faithful nature of the restored F even for finite statistics at the level of smoothness expected for decay spectra and accuracy that may be aimed at currently. Figure 4.7 (Color online) Contour plot of bin to bin Pearson correlation matrix for the restored spectrum when carrying out resampling for the case of the simulated spectrum in Fig. 4.5(d). The solid line guides the eye to indicate average behavior of the width for the main peak in the correlation – the finer details with energy can depend on the assumed original spectrum and even particular simulation. On average, the width grows with energy. The TM with binning for the measured decay energy as well as the RL algorithm with smoothing will generate correlations in results for different bins in the restored energy spectrum. Such correlations can limit the resolution that one can aim at for the restored spectrum. In resampling, we can test the emergence of the inter-bin correlations. This is demonstrated in Fig. 4.7 which shows bin to bin Pearson correlation matrix built from the restored spectrum shown in Fig. 4.5(d). A solid line in the figure guides the eye to show the average behavior for the width of the main peak in the correlation. Beyond variation tied to specific assumptions on the underlying spectrum, the width generally increases with the decay energy, starting at about 0.4 MeV at low 𝐸𝑑 and rising to 1.1 MeV at 𝐸𝑑 ∼ 10 MeV. 4.5 Deblurring 26O decay energy spectrum In the experiment [59, 75], two-neutron unbound 26O was produced via one-proton knockout from a 27F beam. The 26O nucleus decayed to 24O + n +n, and position and time-of-flight measurements of the daughter products were carried out in order to assess their momenta. The 59 0246810Ed2 (MeV)0246810Ed1 (MeV)0.40.20.00.20.40.60.81.0 Figure 4.8 (Color online) Analysis of the measured three-body decay energy spectrum for 26O → 24O + 2𝑛. Panel (a) shows the spectrum measured using invariant mass spectroscopy [59, 75]. Panel (b) shows the deblurred spectrum, as well as the peaks identified for the spectrum with the Deep Neural Network (DNN). Resonances behind the peaks in the spectrum near 0 and 1.3 MeV were also identified for 26O in Ref. [77] (0+ and 2+ states, respectively). Indications of a third peak at about 4 MeV were reported by Caesar et al. [60]. Panel (c) displays contributions from the three peaks identified by DNN, and shown in (b), to the measured spectrum, i.e., after blurring caused by the apparatus. Combination of those contributions (dashed line) matches closely the data (points). The width of the energy bin in processing the spectra is 0.2 MeV. momenta for 24O and two neutrons, measured in coincidence, were used to reconstruct the decay energy spectrum for 26O using the invariant mass technique. Previous invariant mass measurements have observed the ground state of 26O decaying directly into 24O and neutrons very near threshold [78, 77, 60, 59]. A recent experiment measured the half- life of this state to be 5 ps [75]. An excited state, 26O(2+), was also measured with a decay energy of 1.28 MeV above threshold [77]. Indications of a high-lying excited state, at around 4 MeV, were reported in Ref. [60], but Ref. [77] found no evidence of that state. Panel (a) of Fig. 4.8 shows the energy spectrum of the three-body decay of 26O as recorded in the experiment performed at NSCL 60 0246810E0d (MeV)0204060Counts(a)26O24O+2n T. Redpath (2019)0246810Ed (MeV)010203040 Counts (104)(b)restored (deblurring)peaks from DNN95% uncertainty (deblurring)68% uncertainty (deblurring)0246810E0d (MeV)0204060 Counts(c)peak1 from DNNPeak2 from DNNpeak3 from DNNtotalmeasured spectrum Figure 4.9 (Color online) Chi-square per degree of freedom versus the number of peaks included in the fit to the experimental decay energy spectrum (shown in Fig. 4.8(a)). The horizontal axis represents the number of peaks, 𝑛, and each peak is described by three parameters (see the text for details). Increasing the number of peaks is equivalent to increasing the number of parameters for fitting. It may be seen that the three-peak case yields the minimal chi-square per degree of freedom. [59]. Only the first peak, from those mentioned above, is easily seen. The deblurring technique discussed in the previous section helps to extract more information from the measured decay energy spectrum. Panel (b) of Fig. 4.8 displays the spectrum restored from the measured spectrum of the 26O system, using the deblurring method, Eq. (4.6) with an energy-dependent smoothing parameter of Eq. (4.10). The bumps evident in the restored spectrum near 0 MeV and 1.3 MeV, respectively, can be recognized as the 𝐽 = 0+ and 2+ states of 26O nucleus identified in Ref. [77]. We associate the broad peak between 4 and 6 MeV with the third 26O state observed in [60]. The panel (b) in Fig. 4.8 includes peaks that DNN attributed to the original spectrum and corresponding contributions of those peaks to the observed spectrum. We discuss the DNN analysis of decay energy spectrum in Chapter 5. In comparing our method with traditional methods, we have performed chi-square minimization by fitting the measured decay energy spectrum with the resolution-folded BW distribution (see Section. IV). We started with one peak BW function, and gradually increased the number of peaks to five. Each peak is described by three parameters, i.e., amplitude, peak position, and peak width, which implies that the number of fit parameters is three times the number of peaks. In Fig. 4.9, we present the values of chi-square per degree of freedom, 𝜒2/𝑛𝑑𝑓 , versus the number of peaks, 𝑛. A decrease in 𝜒2/𝑛𝑑𝑓 may be observed from 𝑛 = 1 to 𝑛 = 3 and then an increase from 𝑛 = 3 to 5, 61 12345n23452/ndf which implies that three peaks are sufficient to describe the data. It is important to emphasize that the deblurring method does not require any assumption about the number of peaks in the spectrum in order to carry out the restoration, whereas, in the chi-square approach as well as DNN (to be discussed in the next section), one needs to invoke some peaks explicitly (or parameters) in the model. From its side, the deblurring method can suggest the type and number and type of parameters needed in the chi-square fitting or DNN. 4.6 Conclusions and Outlook We applied the deblurring method, successful in optics and employing the RL algorithm, to the restoration of the energy spectrum from the three-body decay of 26O. As presented here, the algorithm requires only the measured distribution in energy and the TM, with elements only labeled by energy, to operate. Two-dimensional distributions of photons are typically employed in optics and such and higher dimensions in nuclear applications can be envisioned. The inversion implicit in the algorithm is largely stabilized by the positive-definite probabilistic nature of the measured and restored distributions and of the TM elements. When significant noise is present in the deblurred distribution, though, a short wavelength instability may develop in the restored distribution in the limit of many restoration iterations. With the relative noise growing with energy, due to fewer counts there, we stabilize that instability with an energy-dependent regularization in the individual restoration steps. Ahead of the data, we tested the method in the restoration of a simulated decay energy spectrum without and with significant noise, as was illustrated in Figs. 4.4 and 4.5. Then, we applied the method to the measured energy spectrum of the three-body decay of 26O. Three peaks were observed in restored spectrum. Two of those were found in the low energy region, at about 0.15 and 1.5 MeV, which may be tied to the previously identified (0+) and (2+) states of 26O. The third peak is located between 4 and 6 MeV in the restored spectrum, and such a peak was previously reported in Ref. [60]. Possible nuclear applications of the deblurring method we described, besides recovering de- cay energy spectra, and the aforementioned 3D distributions in heavy-ion collisions, can include 62 restoration of emitting source distribution from particle correlations in heavy-ion collisions. The emitting source function gives information about spatial geometry and time development of the final stages of reactions [11, 12, 5, 13], as well as their phase-space [79] and thermodynamic characteristics [26]. 63 CHAPTER 5 MACHINE LEARNING MODELS This chapter presents our machine learning framework designed for identifying resonance states from noisy measured decay energy spectra. As an alternative to the deblurring approach discussed in Chapter 4, in this chapter, we employ a deep neural network (DNN) classifier method for resonance state identification in the measured three body decay energy spectrum of 26O system. The comparison between deblurring and DNN classification methods is discussed in Sec. 5.5. The main results discussed in this chapter were recently published in Ref. [50]. The chapter is organized as follows: a brief introduction is given in Section 5.1. In Section 5.2, we introduce the deep neural network (DNN) models. We apply the DNN classification model to the simulated dataset in Section 5.4. In Section 5.3, we apply DNN to experimental data. Finally, we provide a summary in Section 5.6. 5.1 Introduction Machine learning (ML) is a significant research field in modern science. It involves techniques that allow computers to learn from data and make predictions. By using large datasets, ML enables us to extract valuable information about different physical processes and underlying fundamental scientific principles. Ideally, in ML, the models are built to perform tasks, i.e., to learn and make predictions from data without explicit instructions being given. There are numerous machine learning (ML) algorithms, which can be broadly categorized into two groups: supervised and unsupervised learning. Supervised learning involves using algorithms like deep neural networks, decision trees, and random forests, among others. In supervised learning, the training data is labeled, and the algorithms identify patterns within the datasets. This enables the algorithms to make predictions about future events or data that were not part of the training set. On the other hand, unsupervised learning, such as clustering and the Boltzmann machine algorithm, is a method for finding patterns and relationships in the datasets without prior knowledge of the system. In this chapter, we focus on supervised learning, specifically the DNN classification algorithm, and apply the algorithm to build the model for analyzing measured decay energy spectra. For 64 readers interested in learning more about supervised and unsupervised algorithms, we refer them to Ref. [80]. In addition to the RL-based deblurring algorithm, discussed in Chap. 4 and in Ref. [50], we implemented a Deep Neural Network (DNN) classification algorithm to identify the number of resonance states in the decaying nucleus (i.e., 26O →24 O + n + n ) from the measured decay energy distribution. The DNN methods have been popular in face [81] and speech [82] recognition. In the field of particle and nuclear physics, the methods have been applied to particle identification and event selection [83, 84, 85, 86, 87, 88, 89]. In the present work, the DNN uses a training dataset generated from a Breit-Wigner (BW) resonance distribution, folded with the experimental response matrix and sampled according to a Poisson distribution [50]. This process yields a dataset which resembles experimental data. The dataset is labeled and grouped into classes based on the number of resonance peaks introduced in the BW distribution. 5.2 Deep neural network Model A Deep Neural Network (DNN) is a machine learning algorithm that belongs to the family of neural networks. It is characterized by having multiple hidden layers between the input and output layers. Inspired by the information processing of the brain [90], DNNs are designed to learn and represent complex patterns and hierarchical features from data. A key theoretical result for DNNs is the Universal Approximation Theorem, which states that a feedforward neural network with a single hidden layer can approximate any continuous function on a closed and bounded input domain, given a sufficient number of hidden neurons. However, the depth of DNNs (i.e., having more than one hidden layer) significantly enhances their capabilities to model and learn complex relationships between inputs features, and outputs. There are two types of DNN models: DNN for regression problems and DNN for classification. The main difference between them is that the former uses continuous datasets for prediction; for example, the DNN is applied to predict house pricing [91]. The advantage of DNN for regression problems over linear regression is that DNN can learn the complex non-linear relationship between input and output via the activation function introduced in each layer. The latter is used with categorical datasets to decide a class y∗ of an 65 unknown data object x∗ based on the class y of a known data object x [92]. A popular example is the MNIST dataset [93], where DNN is used to classify handwritten digits. In this chapter, we use DNN classification to analyze decay energy spectrum measurements. In this section, we discuss DNN terminologies such as training and testing of the DNN model, feed-forward, and back-propagation, and present the mathematical relationships that make up the DNN in general. 5.2.1 Designing DNN Designing an optimum DNN that does not under-fit or overfit dataset is not an easy task since training a DNN model involves many parameters that need to be optimized. Among them, a few are listed: the number of hidden layers, the number of neurons in each layer, hyperparameters, number of epochs. Here, we briefly describe constructing a DNN that better serves our purpose. Each DNN consists of three main layers; the input layer, hidden layers, and output. The depth of a network is determined by the number of hidden layers it possesses, it is called a deep neuron network when it has more than one hidden layers. It is simple to determine the number of neurons in the input and output layer when you know the problem you are trying to solve. For regression problems, the output layer has just one neuron, and the number of neurons in the input layer is fixed to be the dimensions of the inputs. Whereas, the number of output neurons in the classification equals the number of categorical labels (classes) in your problem. In contrast, determining the number of hidden layers and the corresponding neurons number is complicated. To the author’s knowledge, there is no fixed mathematical theory to determine the number of neurons in the hidden layer, but more often the rule of thumb is used. With the rule of thumb, you try different choices of the number of neurons and hidden layers and select the one which gives the optimum solution. An example of fully connected DNN architecture is illustrated in Fig. 5.1. In the diagram, 𝑎 (0) 𝑚 represents m𝑡ℎ neuron in the input layer (0). The a(𝑙) 𝑗 is the 𝑗 𝑡ℎ neuron in layer 𝑙; 𝑏 (𝑙) is the bias in layer l and 𝑤 𝑗 𝑘 denotes the weight matrix of layer l. While ˜𝑦 is the approximation prediction (solution) of the true function 𝑦. The outputs of any neuron in a layer are given by a non-linear function known as the activation 66 Figure 5.1 Fully connected neural network for classification. function; also, this output serves as an input in the next layer. The activation function is a differentiable function and the most commonly used are: Sigmoid/Logistic hyperbolic tangent Rectifier linear unit (ReLu) function Softmax activation function 𝜎(𝑥) = 1 1 + 𝑒−𝑥 ; 𝜎(𝑥) = 𝑒𝑥 − 𝑒−𝑥 𝑒𝑥 + 𝑒−𝑥 ; 𝜎(𝑥) = max(0, 𝑥); 𝜎(𝑥𝑖) = 𝑒𝑥𝑖 (cid:205)𝑖 𝑒𝑥𝑖 ; (5.1) (5.2) (5.3) (5.4) softmax is used as an activation function of the last layer of a neural network. The softmax activation function is mostly used in multi-class classification problems because it calculates the probability of each class. 67 5.2.2 Training a DNN Training a DNN model involves teaching the neural network to map the input set to the output set. This mapping can be expressed as the following relation: 𝐻 (𝐼1, . . . , 𝐼𝑛) → (𝑦1, . . . , 𝑦𝑛), (5.5) where 𝐻 represents the non-linear function that takes the input set 𝐼1, . . . , 𝐼𝑛 and produces the outputs (the solution of the DNN model), 𝑦1, . . . , 𝑦𝑛. The goal of training the model is to minimize the cost function. The most commonly used cost function in regression problems is sum of square errors (MSE), 𝑀𝑆𝐸 = 1 𝑁 𝑁 ∑︁ 𝑖=1 (𝑦𝑖 − ˜𝑦𝑖)2. In the classification problem, cross-entropy (CEL) loss is the commonly used cost function, 𝐶𝐸 𝐿 = −𝑦𝑖 log( ˜𝑦𝑖) + (1 − 𝑦𝑖) log(1 − ˜𝑦𝑖), (5.6) (5.7) where 𝑦 = 𝑦1, 𝑦2, ..., 𝑦𝑛 be the set of targeted values from the data, and ˜𝑦 = ˜𝑦1, ˜𝑦2, ..., ˜𝑦𝑛 be DNN output. Training a DNN consists of two stages; a feed-forward network that consists of passing a data set from the inputs throughout to the outputs layer, and back-propagation, which is essentially a process to find optimum parameters of the DNN. 5.2.3 Feed forward Feed-forward is one direction flow of information from the input layer to the output. At the output, the network approximates a function ˜𝑦 = 𝑓 (𝑎, 𝑊) which is comparable to the target variable in the data, 𝑦. The measure of how well ˜𝑦 reproduces 𝑦 is given the cost function (see Eq. (5.6) and (5.7)). Mathematically, feed-forward proceeds as the following: let 𝑎 (𝑙) be the set of the inputs in l𝑡ℎ, 𝑎 (𝑙) =𝑎 (𝑙) 0 𝑛 . Then, j𝑡ℎ neuron in the 𝑙𝑡ℎ layer is defined as a weighted sum: , . . . 𝑎(𝑙) , 𝑎 (𝑙) 1 𝑍 𝑙 𝑗 (𝑎) = 𝑛 ∑︁ 𝑘=1 𝑤 𝑘 𝑗 𝑎 (𝑙−1) 𝑗 + 𝑏 (𝑙) 𝑗 , (5.8) where 𝑊 𝑘 𝑗 weights matrix,and 𝑏 is a bias term. In matrix notation, Eq. (5.8) can be simply written: Z = wa + b. 68 (5.9) The output of the neuron is obtained by applying an activation function on 𝑍. Here, we use the Rectified Linear Unit (ReLU) activation function [94] for the hidden layers, and the softmax for the output layer (See those activation functions in the previous subsection). 𝐴𝑙 (𝑍) = 𝜎(𝑍 𝑙 (a, w))), using ReLu activation, 𝐴𝑙 (𝑍 (a, W)) = max(0, 𝑍 𝑙 (a, W)), (5.10) where 𝜎 stands for activation function. It is important to note that the inputs of the current layer are the outputs of the previous layer. Mathematically, the process from the first layer to the output is a recurrent process: The first hidden layer is represented as 𝐴1 = 𝜎(𝑍 1), the second layer as 𝐴2 = 𝜎(𝑍 2) = 𝜎(𝜎(𝑍 1)), and the 𝑙𝑡ℎ layer as 𝐴𝑙 = 𝜎( 𝐴𝑙−1). The outputs in the 𝑙𝑡ℎ layer are summarized as: 𝐴𝑙 = 𝜎(𝑍 𝑙 (𝜎𝑙−1(𝑍 𝑙−1(...𝜎1(𝑍 1)))). (5.11) The goal of training the DNN is to be able to obtain parameters such as w orb that minimize the cost function. The procedure to obtain these optimal parameters is called back-propagation and is discussed in the next subsection. 5.2.4 Back-propagation Back-propagation is a procedure to update the weights and/or bias from the output back to the input, as illustrated in Fig. 5.2 until the optimal values are reached. The procedure starts by computing the chain derivative of the cost function with respect to the weight and bias, also known as the gradient of the error function, which is subsequently used to update the weight, 𝑊, and bias, 𝑏 𝑊 ∗ = 𝑊 − 𝛼∇𝑤 𝐸 , 𝑏∗ = 𝑏 − 𝛼∇𝑏𝐸 ,    69 (5.12) Figure 5.2 The figure shows a schematic illustration of the back-propagation process. The process starts from outputs as shown by red dashed arrows, and propagation through the network to inputs. The propagation refers to as chain derivative of a cost function with respect to weights and bias. where, 𝑊 ∗ and 𝑏∗ represent the updated weight and bias values, respectively. Moreover 𝛼 denotes the learning rate and ∇𝑤 𝐸 = 𝜕𝐸 (𝑊,𝑏) 𝜕𝑊 and ∇𝑏𝐸 = 𝜕𝐸 (𝑊,𝑏) 𝜕𝑏 are the partial derivatives (or gradients) of cost function (i.e., mean squared error, Eq. (5.6)), denoted as 𝐸, with respect to the weights and biases, respectively. These derivatives play a crucial role in the process of finding the optimum weight and bias values and can be expressed as follows: ∇𝑤 𝐸 = 𝑑𝐸 𝑑𝑊 𝑊 𝑙−1 𝜕 𝐴𝑙 𝜕 𝐴𝑙−1 𝑊 𝑙−2 𝜕 𝐴𝑙−1 𝜕 𝐴𝑙−2 · · · = 𝑊 𝑙 𝜕𝐸 𝜕 𝐴𝑙 𝜕 𝐴2 𝜕 𝐴1 𝑊 𝑖(cid:17) (cid:16) 𝑙 (cid:214) 𝑊 2 𝑋, (cid:16) = 𝑖=2 𝑖 (𝑍𝑖) 𝜎′ (cid:17) 𝑋, 𝑙 (cid:214) 𝑖=1 (5.13) In Equation (5.13), 𝜎′ 𝑖 (𝑍𝑖) represents the derivative of the activation function. The differentiation process outlined here pertains to the chain derivative of the cost function 𝐸 with respect to the weight 𝑊. A similar procedure can be carried out for the chain derivative of 𝐸 with respect to the bias 𝑏. It is essential to note that the choice of the ReLU activation function, 𝜎(𝑍) = max(0, 𝑍), is crucial. Other activation functions such as sigmoid and hyperbolic tangent have small derivatives when 𝑍 is large. If several small values of these derivatives are multiplied together, the product will be even smaller. Consequently, the partial derivative with respect to the weight will be nearly zero, making it challenging for the algorithm to update the weight. This issue in deep learning is referred to as the vanishing gradient problem [95]. Furthermore, Eq. (5.13) shows how changing the weights and biases in a network affect the cost 70 function. Different optimization algorithms commonly used in DNN include gradient stochastic descent (GSD) algorithm and adaptive moment estimation (Adam) algorithm to list a few. We adopted the ADAM algorithm, which is currently the most effective and popular DNN. The algorithm is a first-order gradient-based optimization of a stochastic function, and the algorithm requires a little computational memory [96]. For more details about Adam algorithm, we refer the reader to Ref. [96]. When the entire training dataset passes completely through the feed-forward network and back-propagation, it is said to have completed one cycle, which is called an epoch. Each epoch corresponds to updated weights and biases. To check if the optimum values of weights and biases have been reached, we monitor MSE for each epoch. If the MSE is still improving, it implies that the optimal weights and biases have not been reached yet. In this case, we continue to increase the number of epochs until the MSE no longer improves. In the present work, the MSE was found to cease to improve after 50 epochs (see the training/testing accuracy curves in Section 5.3). 5.3 DNN architecture to discover resonance states from the measured decay energy spectrum To complement the results from the deblurring method discussed in Chapter 4, we built a machine learning (ML) tool to classify the number of peaks in the observed decay energy spectrum. A fully connected DNN, schematically illustrated in Fig. 5.3, is defined with the equations: 𝐴(𝑙+1) 𝑖 = 𝑏 (𝑙+1) 𝑖 + 𝑊 (𝑙+1) 𝑖 𝑗 𝑎 (𝑙) 𝑗 , ∑︁ 𝑗+1 𝑎 (𝑙) = 𝑍 ( 𝐴(𝑙)). (5.14) (5.15) where 𝑎 (𝑙) and 𝐴(𝑙+1) are the input and output layers and 𝑊 (𝑙+1) bias of the (𝑙 + 1)𝑡ℎ layer. The non-linear activation function is 𝑍 (𝑥) = ReLu = max(0, 𝑥). The and 𝑏 (𝑙+1) are the weights and 𝑖 𝑗 Relu [94] is commonly used as the activation function in neural network models. The function 𝑓 (𝑥)𝑖 = Softmax = 𝑒 𝑥𝑖 𝑖 𝑒 𝑥𝑖 is used in the output layer to normalize or scale the output so that it may (cid:205)𝑁 be interpreted as a probability [97]. We implement the network using the categorical cross-entropy loss function, 𝐿 = − (cid:205)𝑁 𝑖 𝑦𝑖 log( ˜𝑦𝑖), that is suitable for a multi-class classification problem [98]. Here, 𝑦𝑖 is the 𝑖𝑡ℎ actual value and ˜𝑦𝑖 is the 𝑖𝑡ℎ predicted value (output of the DNN). Then, the Adam algorithm [96], a popular optimizer in DNN models, is used to solve for the optimal weights 𝑊𝑖 𝑗 . 71 Figure 5.3 (Color online) The figure shows a schematic illustration of deep neuron network archi- tecture designed for the classification model. In the present work the input (at the input layer) is the decay energy spectrum, and at the output layer, is a labeled (class) value that tells the number of states in the spectrum. The parameters used to train and test the model are displayed in Table 5.1. The architecture and training specifications of the DNN model are displayed in Table 5.1 and the network design is shown in Fig. 5.3. The DNN classifier is trained using simulated datasets to learn plausible patterns in the decay energy spectra. The dataset is simulated by by folding Breit-Wigner line shapes, Eq. (4.11), with the TM in order to resemble the experimental spectra. We then distort the folded distribution according to Poisson noise to produce a noisy distribution similar to experimental measurements. Note that we utilize the bins, which consist of 50 bins spanning from 0 to 10 MeV, each with a width of 0.2 MeV, as inputs to the DNN. The parameters 𝐸𝑖 and Γ𝑖 in Eq. (4.11), with 𝑖 = 1, . . . , 5, are randomly drawn from a uniform distribution. In this work, we consider the parameters to stem from the range of values displayed in Table 5.1. We divided the training data set into five classes of spectra according to the number of resonances contributing to the decay energy spectrum. The first class, C1, assumed two resonance states with energies 𝐸1 and 𝐸2. The second class, C2, assumed three resonances at 𝐸1, 𝐸2, and 𝐸3. The third class, 𝐶3, assumed three resonances at 𝐸1, 𝐸2, and 𝐸5. The fourth class, C4, assumed four resonances at 𝐸1, 𝐸2, 𝐸3, and 𝐸5. The class C5 contained any other spectrum that does not belong in the first four 72 Table 5.1 The table shows the hyper-parameters used to design the DNN classification model and our assumptions to generate the training data set. The first part of the table displays parameters that made the DNN architecture (i.e., numbers of layers and neurons in each layer). The second part (middle) shows other hyper-parameters and also shows the positions of the peaks (𝐸𝑖) used in Breit-Wigner distribution, Eq. (4.11) which was multiplied with TM to obtain the training set. The last part consists of classes, C1, C2, C3, C4, and C5 created in such a way that each class has distributions with a number of peaks and/or features different from other classes. DNN architecture Layers Input Layer 1𝑠𝑡 Hidden layer 2𝑛𝑑 Hidden layer Output layer Other hyper-parameters Optimizer (Adam) Epoch number (200) Batch number (20) Learning rate (0.004) Learning rate (0.004) prediction of DNN classifier on experimental data Number of neurons Activation function 50 300 500 5 Peak location (MeV) 0.00 ≤ 𝐸1 ≤ 0.30 1.10 < 𝐸2 ≤ 1.80 1.90 < 𝐸3 ≤ 2.70 2.70 < 𝐸4 ≤ 4.00 4.00 < 𝐸5 ≤ 6.00 ReLU ReLU ReLU Softmax Peak width (MeV) 0.08 ≤ Γ1 ≤ 0.50 0.50 < Γ2 ≤ 1.10 1.10 < Γ3 ≤ 1.30 1.30 < Γ4 ≤ 1.80 1.80 < Γ5 ≤ 2.10 Class label 1. C1 2. C2 3. C3 4. C4 5. C5 Description Peak Position (resonances states) 𝐸1 and 𝐸2 𝐸1, 𝐸2, and 𝐸3 𝐸1, 𝐸2, and 𝐸5 𝐸1, 𝐸2, 𝐸3 and 𝐸5 elsewhere classes. For convenience, we assign to class C5 four kinds of spectra: spectra with one peak at E1, spectra with two peaks at E1 and E3, spectra with three peaks at E1, E4, and E5 and spectra with four peaks at E1, E2, E3, and E4. It is important to note that, in choosing values for 𝐸1,2,3,5, we made sure to include all the 26O states that were previously reported (see Refs. [77, 60]). The mean values of 𝐸1,2,3,4,5 have been equal to about 0.15, 1.50, 2.40, 3.35 and 5.00 MeV, respectively. We generated 6,000 spectra for each class, producing a data set containing 30,000 simulated spectra to train and test the model. From these, 60% of the data set was used for training, and 40% was used for testing. The optimal model was achieved for the values of the parameters displayed in Table 5.1. The model’s performance was evaluated based on the training/testing accuracy curves illustrated in the panel (a) of Fig. 5.4. The performance of the DNN classifier, as shown in panel (a) of Fig. 5.4, was assessed in terms of accuracy. The accuracy, as the metric used to evaluate the classification model, is the number of correct predictions out of the total number of predictions. An accuracy equal to 1 stands for 73 Figure 5.4 (Color online) DNN model to identify resonance states from measured decay energy spectrum. Panel (a) indicates the training and testing accuracy of the model. The curves converge at ≈ 0.75 on both data sets, which means the model predicts 75% of the data set correctly. Panel (b) represents the confusion matrix which tells how well the DNN classifier was able to classify spectra: C1, C2, C3, C4, and C5 are five classes we used to train the model, and the detail about each class is discussed in the text and Table 5.1. Panel (c) displays the ratio of distributions predicted to belong in a given category (C𝑖, 𝑖 =1, 2, 3, 4, 5) over the total number of distributions used in the DNN model prediction. The ratio helps to estimate the class where the measured spectrum fits. We found that C3 has the highest fraction, which suggests that there is a high chance for the experimental decay energy spectrum of 26O system decaying 24O+ 2n from invariant mass spectroscopy measurements belongs to C3. the perfect performance of a model, and 0 stands for complete failure. As shown in the figure, the model achieves an accuracy between 0.7 and 0.75 after training for 40 epochs. Panel (b) in Fig. 5.4 displays the confusion matrix, which gives information about the classifier’s performance in assigning each simulated spectrum to the correct class. The elements on the diagonal represent a normalized number of ideally classified spectra, and the off-diagonal elements represent the misclassified spectra. The first and the second class show a high number of misclassified spectra because those two classes have similar peaks in the low energy regime (< 2.5 MeV). Hence, it is harder for the network to distinguish them, especially when the data set fluctuates significantly. Each element of the confusion matrix is estimated from representative test sets, counting the number of distributions assigned to each class and normalized by the number of spectra in that class. After the DNN is trained, we use it to in classify the experimental spectrum. For example, if the assigned class is C1, this means that the spectrum is perceived to have two states around the 74 positions described before. If it is C4, the spectrum is recognized as having four peaks. However, we only have one measured decay energy spectrum from the experiment investigating the three-body decay of 26O into 24O and 2n. A fair prediction is expected, when enough samples are passed to the DNN model. For that reason, we carry out error resampling for the measured spectrum to obtain a dataset that one can use in the model prediction. With this process, we have generated 10,000 samples of distributions and estimated the classes to which each spectrum from the resampling belongs. We evaluated the number of distributions predicted to be in a given class as a fraction of the distributions in that class per the total number of distributions, see panel (c) of Fig. 5.4. With this, the number of resonance states most likely there in the measured decay energy spectrum corresponds to the class with the highest fractional value. As evident in Fig. 5.4(c), more than 75% of the total distributions used in the prediction belong to class C3, which suggests within our statistical framework the presence of three peaks in the measured spectrum. The locations of those peaks are approximately equal to the mean values of parameters 𝐸𝑖 of C3 at 0.15, 1.50, and 5.00 MeV. The mean values of half-widths with which the classifier sorts the three peaks, are, respectively, 0.29, 0.80, and 1.85 MeV. Finally, these peaks correspond to the resonance states of 26O reported in Refs. [77, 60]. 5.4 DNN Classification model for resonance peaks detection In this section, we test the DNN classification model developed in the previous section (Sec- tion 5.3) on the distribution simulated in such a way that it resembles the spectrum measured by MoNA-Lisa detectors. This distribution is simulated assuming that the resonance state follows the Breit-Wigner (BW) line shape of peak position (decay energy) 𝐸𝑖 and width Γ𝑖. The expression of the BW distribution is exhibited in Eq. (4.11). The decay energy distributions generated from BW are multiplied with the simulated response matrix of the MoNA experimental setup [50] that incorporates the detector’s acceptance and resolution and then distorted with Poisson noise to end up with spectra which resemble the measured ones. It is worth mentioning that the response matrix used here is built based on three body decay of 26O experiment [99]. In generating the distribution used to test the DNN classification model, we assume three peaks at 0.3, 1.4, and 4.5 MeV, with 75 Figure 5.5 The figure displays simulated decay energy distributions. The solid line corresponding to the energy, 𝐸𝑑, in the horizontal, represents the distribution obtained using the Breit-Wigner distribution (see Eq. (4.11)). The dots represent energy distribution ( corresponding to the energy, 𝐸′ 𝑑, in the horizontal ) was obtained by multiplying the distribution shown as a solid line by the transfer matrix shown in Fig. 4.2 and adding Poisson noise to produce fluctuations in the distribution. widths of 0.3, 1, and 1.5 MeV, respectively. This distribution is shown in Fig. 5.5: the solid line (yellow) represents the original BW distribution and dots (blue) show the distribution obtained by multiplying the original BW distribution (i.e. solid line) with the transfer matrix and subsequently fluctuating it with Poison noise. Here, we aim to use the latter in DNN to request the number of peaks in the original distribution. To do that, we use the DNN classification model developed in the previous section (Section 5.3) and follow the procedure used for the measured decay energy spectrum. The corresponding results are displayed in Fig. 5.6, where the bar shows the percentage of the distributions predicted to members of each of the five classes. As distributed in the figure, around 73% of all distributions belong to class C3; this tells us that the considered distribution has 3 peaks and the mean value of their parameters, peak positions and widths are 0.3, 1.4, and 4.5 MeV, peak position and 0.3, 1, and 1.5, peak width. Those parameters agree with the parameters we used to create the original distribution in Fig. 5.5. 76 0246810Ed, E0d (MeV)0.00.51.01.52.02.53.0Counts (104) original spectrumdataset(x400) Figure 5.6 The bar plots displays the ratio of distributions predicted to belong in a given category (C𝑖, 𝑖 =1, 2, 3, 4, 5) over the total number of distributions used in the DNN model prediction. The class with a ratio high than 0.6 is estimated as the class where the measured spectrum fits ( i.e., here, 𝐶3 fulfills that condition). 5.5 Comparison of the DNN method with deblurring method In Chapter 4, we employed the deblurring method discussed in Chapter 3 to restore the measured decay energy spectrum of the 260 system. The restored spectrum provides insightful information regarding the energy levels of 260. In this chapter, we introduce a DNN classification model for identifying the number of resonance peaks in the measured decay energy. The primary distinction between the two methods lies in their approach. The proposed deblur- ring method focuses on restoring the original spectrum’s features without assuming the number of states it contains. Conversely, in the DNN classifier, we assign probabilities to hypotheses with varying numbers of states in the original spectrum. Referring to the results of both methods, as de- picted in Fig. 4.8, when applied to the same dataset, these two approaches both test and complement each other. The deblurring process allows us to distinctly observe the resonance peaks, whereas the DNN enables us to identify the peak count in the measured spectrum, along with estimating 77 associated parameters such as peak position (𝐸𝑖) and width (Γ𝑖). To summarize, using deblurring to restore distorted decay energy results in a smooth spectrum with clearly defined peaks corresponding to the resonance states of 260, all achieved at a low cost. 5.6 Summary We constructed a deep neural network (DNN) classification model with the same objective as the deblurring technique: to identify resonance states of 26O from the measured decay energy spectrum. Initially, the DNN was trained and tested using a simulated dataset. Subsequently, we employed the DNN model to estimate the presence of three peaks in the spectrum, located approximately at mean positions of 0.15 MeV and 1.50 MeV for the first and second peaks, and at around 5.00 MeV for the third peak. The half-widths of these three peaks were determined to be approximately 0.29 MeV, 0.80 MeV, and 1.85 MeV, respectively. The concurrence observed between the deblurring and DNN classification employed in our analyses suggests that there might be three resonance states of 26O influencing the measured decay energy spectrum. In the future, it will be interesting to extend the DNN method by building a DNN that solve an inverse problem to estimate the parameters of the resonance states. 78 CHAPTER 6 TRANSPORT MODEL SIMULATION FOR TWO-PARTICLE SOURCE FUNCTION The dynamics of the nuclear collision system can be understood by solving appropriate transport equations. In the literature, various transport models such as Boltzmann-Uehling-Uhlenbeck (BUU) and Quantum Molecular Dynamics (QMD) were developed to study heavy ion collision dynam- ics [15] and the reference within. In this chapter, we will focus on the BUU transport approach developed by Danielewicz [100]. The BUU code allows us to calculate the single-particle phase- space distribution, which is then utilized to determine the source function. However, since BUU involves several adjustable parameters, one can likely learn correlation data by comparing the source function obtained from BUU with the one obtained from the deblurring method discussed in Chapter 3. The parameters yielding the correct source function, which reproduces the mea- sured correlations, can then be used to infer the nuclear equation of state (EOS). Equation of State (EOS) is the mathematical relationship that connects thermodynamic quantities such as energy, temperature, density, and isospin, which is defined as the ratio between the number of protons and neutrons in nuclear matter. For nuclear matter, the EOS plays a crucial role in determining the macroscopic properties of laboratory nuclei and neutron-star crusts. The saturation density and energy of symmetric nuclear matter, where the number of neutrons and protons is equal, are reasonably well determined from the masses and radii of stable nuclei, as the differences in neutron and proton numbers are not significant. In this chapter, we employ the BUU transport model to calculate the proton-proton (p-p) source function and subsequently utilize it to compute p-p correlations in heavy-ion collisions. We calculate p-p correlation in the low-energy heavy-ion reactions (E/A = 50-100 MeV), i.e., 36Ar +45 Sc reaction at 80 MeV/nucleon. Additionally, we test the model on heavier systems, namely the 129Xe + 197Au reaction at 80 MeV/nucleon. However, for this system, we only simulate the source function as experimental measurements are currently unavailable. We then compare the BUU-calculated source function with the deblurring calculation for proton-proton correlation and examine the sensitivity of the source function to the momentum-dependent soft and stiff EOS. 79 Section 6.1 discusses the BUU transport model including testing the time evolution of the reaction from the initial condition at t = 0 fm/c up to t ≥ 80 fm/c. We systematically examine the effect of EOS on the reaction density as a function of time. In Section 6.2, we delve into the calculation of p-p correlations and compare them with the deblurring calculations. Section 6.3 focuses on the quadrupole angular component of the BUU source function and investigates its sensitivity to the EOS. Finally, we summarize the chapter in Section 6.4. 6.1 BUU model The BUU model is a one-body transport model that has been extensively used to describe the dynamics of the heavy ion collision at the intermediate energies in the semi-classical limit. In this transport model, the system is characterized by a single particle density at any time. In the semi-classical limit, the single particle phase space distribution, 𝑓 = 𝑓 (p, r, 𝑡), can completely specify the state of the system if the distribution is given for all particles. We use the BUU code to solve for the one body phase space distribution 𝑓 that satisfies the Boltzmann equation [100, 101], (cid:16) 𝜕 𝜕𝑡 + P 𝑚 .∇𝑟 − ∇𝑟𝑈 (r).∇𝑝 (cid:17) 𝑓 = 𝐼𝑐𝑜𝑙 (𝜎𝑖𝑛, 𝑓 ). (6.1) The left-hand side of the Eq. (6.1) accounts for the motion of the particle in the mean-field potential velocity v = P/𝑚 = ∇𝑝𝐸 and force ∇𝑟𝑈 (r); where 𝑚, P, and 𝐸 are mass, momentum, and total energy. On the right hand side, 𝐼𝐶𝑜𝑙 is the collision term that accounts for Nucleon-Nucleon collisions 𝐼𝑐𝑜𝑙 = ∫ 𝑔 ℎ3 ∫ 𝑑3P2 𝑑Ω 𝑑𝜎12 𝑑Ω (cid:104) ˆ𝑓1 ˆ𝑓2 𝑓1′ 𝑓2′ − 𝑓1 𝑓2 ˆ𝑓1′ ˆ𝑓2′ (cid:105) , 𝑣12 (6.2) where 𝑓𝑖=1,2 is single particle phase-space distribution and ˆ𝑓 = 1 + 𝜏 𝑓 ; 𝜏 = −1, 0, 1 correspond Fermi-Dirac, Boltzmann, and Bose-Einstein statistic, 𝑔 is spin degeneracy, and 𝑑𝜎12 𝑑Ω is differential cross section for scattering angle Ω . The primed indices represent the final momenta while the initial momenta are P1 and P2 with relative velocity 𝑣12, relative momentum q12 = P1 − P2 for identical masses and q12 = 𝜇(v1 − v2) for distinguish masses, where 𝜇 is a reduced mass. The total momentum P12 = P1 + P2. The final relative momentum is determined by the scattering angle, Ω = (𝜃, 𝜙) so that the final relative momentum is 𝑞1′2′ = 𝑞12(sin 𝜃 cos 𝜙, sin 𝜃 sin 𝜙, 𝑐𝑜𝑠𝜃). 80 Notably, for the elastic collision, the total momentum remains constant while the final and initial relative momenta are equal (i.e., q12 = q1′2′, hence q1′2′ is determined by scattering angle. The initial locations and momenta of the nucleons in the colliding nuclei are assumed to populate uniformly a nucleon Fermi sphere with a Fermi momentum determined from local density using Thomas-Fermi approximation. The local density from Thomas Fermi approximation is utilized as initial conditions to the BUU transport equation. This density is obtained by solving the coupled Thomas-Fermi equations, as detailed in Ref. [102]: ˜𝐸𝑛 (𝐾 𝐹 (𝜌𝑛)) + 𝑎1∇2(𝜌/𝜌0) + 𝑎𝑇 4 𝛿(𝑟) − 𝜇𝑛 = 0, ˜𝐸 𝑝(𝐾 𝐹 (𝜌 𝑝)) + 𝑎1∇2(𝜌/𝜌0) + 𝑎𝑇 4 𝛿(𝑟) + 𝑈Coul − 𝜇 𝑝 = 0,    (6.3) These equations are obtained by requiring that the energy is at a minimum in the ground state of the nucleus with a definite number of protons and neutrons, i.e., 𝜌𝑖 = 𝜌0. Here, 𝜇𝑖 (not to be confused with reduced mass) is the chemical potential, ˜𝐸𝑖 is the single-particle energy as a function of the ˜𝐸𝑖 = √︁𝐾 𝐹 + 𝑚𝑖 (𝜌𝑖). In Eq. (6.1), 𝑈 represents the mean-field potential Fermi momentum 𝐾 𝐹, with the isospin asymmetry 𝛿(𝑟) = (𝜌𝑛 − 𝜌 𝑝)/(𝜌) and 𝜌 = 𝜌𝑛 + 𝜌 𝑝. The index 𝑖 stands for the particle species (i.e., neutron,𝑛 or proton, 𝑝). The values of parameters, such as 𝑎1, 𝑎𝑇 , and other parameters we used, can be found in Ref. [102]. The Thomas-Fermi equations are solved by employing the boundary condition of neutron and proton densities: 𝜌𝑛,𝑝 (𝑟 → ∞) = 0 and 𝑑𝜌𝑛, 𝑝 𝑑𝑟 |𝑟 = 0 = 0. These equations are numerically solved by starting with initial values of 𝜌𝑛/𝑝 and iteratively adjusting the chemical potentials [101]. The momentum dependent is introduced in BUU model, Eq. (6.1) through the potential. The momentum-dependent potential that we used here was discussed in Ref [102] and is shown in Fig. 6.1, each curve corresponds to the potential parameterizations, but here we only show the value of the ratio, 𝑚∗/𝑚 for each curve but the complete set of parameters can be found in Ref. [102]. In the ratio. 𝑚∗/𝑚, 𝑚∗ = 𝑃𝐹 𝑣𝐹 , where 𝑃𝐹 and 𝑣𝐹 are fermi momentum and velocity, respectively and 𝑚 mass of nucleon. In Fig. 6.1, the linear dependence of the potentials on energy in the low-energy regime is interpreted as an indication that the effective mass provides a good characterization of those potentials in that region. 81 Figure 6.1 Real part of the optical potential, as a function of nucleon energy, at different densities, 𝜌. Different curves represent different parametrizations of the mean-field (MF) potential, the dashed line represents the potential for the ratio, 𝑚∗/𝑚 = 0.65, and the solid line represent the potential for 𝑚∗/𝑚 = 0.70. The parameters of the potential can be found in Ref. [102]. For the ratio 𝑚∗/𝑚 = 0.65, the first set of parameters was used. Stars represent the MF potential without It is notable that the potential grows as 𝜌 increases. The figure was momentum dependence. extracted from Ref. [102]. Fig. 6.2 shows the density profile of proton (dashed line) and neutron (dash-dotted line) for intermedium mass nuclei, 40Ca, 48Ca, 40Ar, and 45Sc obtained from TF equation; such nuclei are suitable in the low energy heavy ion collision experiments. The solid lines represent the total density profile 𝜌. As expected from the TF equation, the density decreases slowly in the range of low 𝑟, but the tail falls off rapidly. Notably, the effect of coulomb interaction is seen in the difference between proton and neutron density in the case of a symmetric nucleus (i.e., 40Ca). The BUU model is a complex model which requires many input parameters such as symmetry energy, momentum dependence, impact parameter, and in medium cross-section to list a few. Changing one parameter can affect several observables. Previously, in Ref. [14], It was reported that the two-proton source function depends on in-medium cross-section. In fact, in a coarse estimate, the size of the source function is given by, 𝑅 ∝ √ 𝜎 𝐴 [14]. Also, the sensitivity of soft and stiff symmetry energy proton-proton correlation was explored for 52Ca+48Ca reaction in Ref. [24]. Here, we explore the effect of different factors such as impact parameter, and momentum-dependent soft and stiff symmetry energy EOS on the two-proton source functions. For soft and stiff EOS, we consider both momentum-independent and momentum-dependent cases in the BUU model. Testing the effects of these parameters on the reaction dynamics can be done at different stages 82 0200400m [MeV]10050050Uopt [MeV]=0.2 fm30200400m [MeV]10050050=0.15 fm30200400m [MeV]10050050=0.1 fm3m*/m=0.70m*/m=0.65 Figure 6.2 The figure displays the density profile versus distance for nuclei: 40Ca, 48Ca, 40Ar, and 45Sc obtained from TF equation. The dashed line represents proton density, 𝜌 𝑝, the dash-dotted line shows the neutron density, 𝜌𝑛, and the solid line shows the total density, 𝜌 = 𝜌𝑛 + 𝜌 𝑝. during the collision, starting from the initial condition to the freeze-out moment (i.e., the moment the last particle leaves the system). Additionally, one could include different properties of the equation of state (EoS), such as symmetry energy and momentum dependence, in the mean-field potential. The potential symmetry energy we use has density dependence, 𝐸sym(𝜌) 𝑝𝑜𝑡 = 𝑆0(𝜌/𝜌0)𝛾, where 𝜌0 = 0.16 fm−3 is the density of nuclear matter, 𝛾 is a constant, and typically 𝛾 < 1 corresponds to a soft EoS, while 𝛾 ≥ 1 corresponds to a stiff EoS and 𝑆0 is a parameter. We assessed the effect of different EOSs on proton and neutron density profile, Fig. 6.3 (a) and (b) illustrate the proton and neutron densities of 40Ar and 45Sc nuclei obtained from the solution of TF equation for soft and stiff EOS. In general, for both soft and stiff EOS, in the momentum- dependence case, the proton(neutron) density profile is decreased in the range of low values of 𝑟, and then in the intermediate values, the density profile falls faster for the case of momentum- 83 r (fm)0.000.050.100.15(r)40Car (fm)48Capn0246r (fm)0.000.050.100.15(r)40Ar0246r (fm)45Sc (a) (b) Caption Figure 6.3 Proton and neutron density for Ar, and Sc nuclei. The densities were obtained from the TF equation where different parameters of the nuclear equation of state were varied. Figure (a) represents soft EOS for both momentum-dependent and independent cases and Figure (b) is the same as (a) but for stiff EOS. independent soft and stiff EOS. We do not see much difference in the densities profile of stiff and soft EoS, however, the gap between densities resulting from the momentum-dependent and momentum-independent models is small for stiff compared to soft EOS. Progressing in the BUU simulation, we evaluated the time evolution of the reaction during heavy- ion collisions and tested the effect of input parameters mentioned earlier (i.e., impact parameter and momentum-dependent soft and stiff EOS of symmetric matter). Figure 6.4 shows the contour plots for density in the 129Xe+197Au reaction at 50 MeV/A in the reaction plane at 0, 50, 100, and 150 fm/c times, obtained from the BUU simulation with a stiff EOS. The rows represent the different impact parameters: 𝑏 = 1.3, 3.1, and 6.2 fm, respectively, from top to bottom. At t = 0 fm/c, the 84 012345r (fm)0.000.020.040.060.08(r) (fm3)40Arp, softn, softp, soft mo-depn, soft mo-dep012345r (fm)45Sc012345r (fm)0.000.020.040.060.08(r) (fm3)40Arp, stiffn, stiffp stiff mo-depn stiff mo-dep012345r (fm)45Sc three systems (in three rows) are distinguished by their impact parameters. Around 50 fm/c, the projectiles and targets emerge completely, and it is clear that the particles have started leaving the system at around 100 fm/c. The neck formation is visible for the case of 𝑏 = 6.2 fm. Finally, around 150 fm/c, the formation of rings is present for 𝑏 = 1.3 fm and 𝑏 = 3.1 fm. At this point, the system is separated and forms two fragments for 𝑏 = 6.2 fm. Furthermore, we investigated the effect of switching to momentum-dependent behavior in BUU for stiff and soft symmetry energy EOSs on the reaction during the collision. As an example, we considered two systems: a heavy system, 129Xe+197Au, and a light system, 36Ar+45 Sc, at 80 A MeV. The density contour plots of these reactions are shown in Fig. 6.5 and 6.6. As displayed in Figs. 6.5a, 6.5b, 6.6a, and 6.6b, at a time equal to 80 fm/c, the difference between momentum-dependent (mo- dep.) and momentum-independent (mo-indep.) EoSs become apparent. In addition, we can make a rough comparison between Fig. 6.4 and Fig. 6.5. For example, the parameters used to simulate the second row of Fig. 6.4 and the first row of Fig.6.5 (b) are the same, except for the beam energy of 50 AMeV for the former and 80 AMeV for the latter. However, it is clear that after 80 fm/c, the system shown in Fig. 6.5 (b) has expanded more than the system shown in Fig. 6.4 has after 100 fm/c. Moreover, the effect of the momentum (in)dependent stiff and soft EOS on the dynamics of the system is also notable. For example, in Fig. 6.5 and 6.6, at a time equal to 80 fm/c, we can observe the difference between the density contour plots for the momentum-independent and momentum-dependent EoSs. However, for the small system (Fig. 6.6), this difference is not significant. 6.2 Two-proton correlation function Two-particle correlation can be measured in heavy-ion collisions for the final state emission. By utilizing the BUU transport simulation discussed in the previous section, one can compute the two-particle source function which is used together with the proton-proton kernel to calculate the correlation function. In this section, we discuss the two-proton correlation function for the source function obtained using the BUU model. Then we compare these correlations with the ones obtained using the deblurring method, discussed in Chapter 3. We explored the source function for 85 Figure 6.4 Density contour plots for the reaction 129Xe+197Au at 50 MeV in the reaction plane, at selected times during the reaction. Results displayed in each row correspond to a value impact parameter, 𝑏 = 1.3 fm, 3.1 fm, and 6.2 fm for the first, second, and third low respectively. We considered the momentum-independent stiff equation of state in the BUU simulation. two systems that undergo central collision; one light,36Ar +45 Sc and a heavy system, 129Xe +197 Au at 80 AMeV but we calculated the 𝑝 − 𝑝 correlation function for light system (i.e., 36Ar +45 Sc reaction), since for this system, we have access to experimental data of our interests (i.e., low energy reactions). Furthermore, we test the sensitivity of the source function on the equation of state, we consider soft and stiff both momentum-dependent and momentum-independent cases. The p-p correlation we discussing here is theoretically defined according to the Koonin-Pratt equation and can be written as, ∫ 𝑑3𝑟 [∫ 𝑑3𝑅 𝑓 (P/2, R + r/2, 𝑡′) 𝑓 (P/2, 𝑅 − 𝑟/2, 𝑡′)] |Ψ𝑝 𝑝 (q, r)|2 | ∫ 𝑑3𝑟 𝑓 (P/2, r, 𝑡′)|2 , where, ∫ 𝑑3𝑅 𝑓 (P/2, R + r/2, 𝑡′) 𝑓 (P/2, 𝑅 − 𝑟/2, 𝑡′) | ∫ 𝑑3𝑟 𝑓 (P/2, r, 𝑡′)|2 = 𝑆𝐵𝑈𝑈 (𝑟), (6.4) (6.5) is the source function calculated directly from the relative position (r) and total momenta (P) of the particles at the time when the last particle is emitted [14] and references within. Here, R = 0.5(r1 + r2) is the center of mass coordinate and |Ψ𝑝 𝑝 (q, r)|2 p-p kernel at relative momentum, 86 b=6.2 fmb=3.1 fmb=1.3 fm (a) (b) Figure 6.5 Density contour plots for the reaction 129Xe+197Au at 80 MeV in the reaction plane, at selected times during the reaction. Panel (a) represents densities for a soft EOS, where the first row shows the results for momentum-independent and the second row for momentum-dependent soft EOS. Panel (b) is similar to panel (a) for the case of stiff equation of state results. Results in the first and second rows correspond to densities for the soft momentum independent and dependent EoS, respectively. 87 10010soft mo-indep.EOS0 fm/c40 fm/c80 fm/c1001010010soft mo-dep.EOS1001010010Z (fm)X (fm)10010stiff mo-indep.EOS0 fm/c40 fm/c80 fm/c1001010010stiff mo-dep.EOS1001010010Z (fm)X (fm) (a) (b) Figure 6.6 The similar as Fig. 6.5 for 36Ar +45 Sc reaction. 88 10010soft mo-indep.EOS0 fm/c40 fm/c80 fm/c1001010010soft mo-dep.EOS1001010010Z (fm)X (fm)10010stiff mo-indep.EOS0 fm/c40 fm/c80 fm/c1001010010stiff mo-dep.EOS1001010010Z (fm)X (fm) q and r, Ψ𝑝 𝑝 (q, r) is calculated as discussed in Refs. [27, 54, 103]. The phase-space distribution of particles with momentum P, position r = r1 − r2, and at time 𝑡′ after both particles have been emitted is described by a Wigner distribution 𝑓 (P, r, 𝑡′) [14] and reference within. This Wigner distribution can be expressed as a function of the single-particle emission function 𝑔(P, r, 𝑡) for the last emission points of particle at time 𝑡. Therefore, the function 𝑓 (P, r, 𝑡′) is given by 𝑓 (P, r, 𝑡′) = ∫ 𝑡′ −∞ 𝑑𝑡.𝑔(P, r − P(𝑡′ − 𝑡)/𝑚, 𝑡), (6.6) where 𝑚 is the two-particle reduced mass. The dependence of position and time in Eq. (6.5) shows that the source function describes the space-time evolution of the reaction (i.e., see Fig. 6.4-6.3). The BUU transport model has predicted two types of emission: pre-equilibrium emission and late emission [104]. At pre-equilibrium emission (i.e., emission after a few fm/c), particles are emitted with high velocity, and the resulting source function is dominated by spatial properties. On the other hand, the particles emitted at a later time, also known as slow emission, originate from the sequential decay of unstable nuclei produced during the multifragmentation reaction. Both the pre- equilibrium and slow emissions contribute to the shape of the source function at different distances. The pre-equilibrium emissions primarily affect the shape of the source function at short distances, while the slow emissions contribute to the tail of the source function, which is approximated as an exponential decay [14]. 89 Figure 6.7 The two-proton source function, obtained from BUU simulation for 129Xe+197Au col- lision at 80 MeV/nucleon and b=3.1 fm, is shown as a function of relative position, r for 200-400 MeV/c cut in total momentum. The top row displays stiff and soft cases: solid line for momentum- dependent EoSs and dash-dotted line for no momentum-dependent EoSs. The bottom row compares soft (dash-dotted line) and soft (solid line) cases, with the left and the right panels showing the source functions for momentum-dependent and no momentum-dependent EoSs, respectively. As observed in Fig. 6.7 and 6.8, the tails of the source function presented in the figures decay rapidly after a few fm. This suggests that the contribution from slow emissions may be missing, which is not surprising, we discuss the contribution of the slow emissions to the p-p correlations later in this section. Concerning the sensitivity of source function on different choices of EoSs, in Fig. 6.7, we observe the difference between soft momentum-dependent and momentum- independent EoS as well as between soft and stiff momentum-dependent EoSs. But there is no difference observed between EoSs in. Fig. 6.8. 90 0.00.10.20.30.4stiff mo-dep EOSstiff mo-indep EOSsoft mo-dep EOSsoft mo-indep EOS02468100.00.10.20.30.4stiff mo-dep EOSsoft mo-dep EOS0246810soft mo-indep EOSstiff mo-indep EOSr (fm)SBUU(r)[103 fm3] Figure 6.8 Similar as Fig. 6.7 but for 36Ar +45 Sc reaction at 80 MeV/nucleon with 𝑏 =1.9 fm. We employ the source function obtained from a BUU simulation for a central36Ar +45 Sc collision at 80 A MeV and impact parameter 𝑏 = 2 fm with a 𝑝 − 𝑝 Kernel in the Koonin-Pratt equation to calculate the correlation functions depicted in Fig. 6.9. In this figure, the triangles represent correlations measured in the central36Ar +45 Sc reaction at 80 MeV/A, with proton pairs’ total momenta restricted to the range of 200-400 MeV/c those data can be found in Ref. [14] and the references within. The lines in the figure correspond to correlations computed using sources from BUU simulations. In panel (a), the (dashed) solid lines correspond to the momentum (in)dependent stiff EOS, while in panel (b), the (dashed) solid lines represent the momentum (in)dependent soft EOS. Panel (c) presents the comparison between the momentum-dependent stiff (dashed line) and soft (solid line) EoS, while panel (d) displays the comparison between the momentum-independent stiff (dashed line) and soft (solid line) EoS. 91 0.00.51.01.52.0stiff mo-dep EOSstiff mo-indep EOSsoft mo-dep EOSsoft mo-indep EOS02468100.00.51.01.52.0stiff mo-dep EOSsoft mo-dep EOS0246810soft mo-indep EOSstiff mo-indep EOSr (fm)SBUU(r)[103 fm3] Figure 6.9 The 𝑝 − 𝑝 correlation functions obtained through BUU simulation for36Ar +45 Sc collision at 𝐸/𝐴 = 80 MeV with b=1.9 fm impact parameter are displayed as a function of relative momentum. The lines represent the correlations function that correspond to different equations of state (EoSs) in separate panels, with mo-indep. and mo-dep standing for momentum independent and momentum dependent, respectively (refer to the Legends). Those correlations are obtained using the source function defined in Eq. (6.7). Experimental data, extracted from Ref. [14], is indicated by triangles. However, the BUU simulation correlations differ from the data in all panels. The shape of these correlation functions reflects the interplay among several factors, including the short-range nuclear interaction, the antisymmetrization effect, and the long-range Coulomb interaction between the emitted protons. The attractive S-wave nuclear interaction plays a crucial role in generating the observed prominent peak at a relative momentum of 𝑞 = 20 MeV/c. This peak is incorporated in the correlation function, R(q), through the kernel, which is computed from the radial proton-proton relative wave function. The kernel’s contribution is particularly significant at short relative distances, where the S-wave nuclear interaction dominates. On the other hand, the 92 0.00.51.01.5R(q)+1(a)data stiff mo-indepstiff mo-dep(b)soft mo-indepsoft mo-dep020406080q (MeV/c) 0.00.51.01.5R(q)+1(c)stiff mo-depsoft mo-dep020406080q (MeV/c) (d)stiff mo-indepsoft mo-indep minimum at 𝑞 = 0 MeV/c in R(q) arises from the long-range repulsive Coulomb interaction. As depicted in Fig. 6.9, there is no notable distinction between the equations of state (EoSs). However, it can be observed that the momentum-dependent EoS tends to yield slightly higher values at the peak for both the stiff and soft EoS. Regarding the comparison between the results for soft and stiff EoSs, it appears that for the soft EoS the correlations exhibit a slightly higher value than the ones for the stiff EoS at the peak position. Nevertheless, due to the insignificant difference between the EoSs, solely manipulating different EoSs is insufficient to reproduce the observed data accurately. The difference observed between correlations simulated from BUU and the data in Fig. 6.9, is expected because BUU describes early emissions but not late emissions and both emissions contribute to the data. In this section, we propose a way to address this issue. In the existing literature, researchers have made attempts to correct the proton-proton correlation function obtained from BUU simulations by multiplying it with a factor, 𝜆, resulting in 𝜆𝑅(𝑞) (see Refs. [14, 104] or the two-parameter source model in Chapter (2) for more details). We attempted a similar approach by using 𝜆 = 0.5 to reproduce the experimental p-p correlations measured in the Ar+Sc reaction at E/A = 80 MeV. We observed that this could successfully reproduce the correlations in the high relative momentum range (𝑞 > 20 MeV/c), but our model struggled to predict the shape of the correlations in the low 𝑞 regime (𝑞 < 20 MeV/c), where particle contributions from secondary decay played a significant role. Furthermore, it is crucial to emphasize that secondary decay contributes to the tail of the source function. However, our chosen transport model, the BUU model, does not generate particles through secondary decay. As a result, the correlations simulated using the BUU model do not include contributions from these decays. As previously mentioned, attempting to address this issue by simply multiplying the parameter 𝜆 to the correlations obtained from the BUU model proved inadequate for fitting these correlation functions. Consequently, we propose a novel approach to incorporate the effects of secondary decay into the source function. This involves adding a parameterized exponential tail to the source function derived from BUU simulations, as defined in Eq. (6.5); 𝑆(𝑟) = 𝑆𝐵𝑈𝑈 (𝑟) + 𝐶 (𝑟/𝐵)2 exp(−𝑟/𝐵), (6.7) 93 Figure 6.10 The 𝑝-𝑝 source function for 36Ar +45 Sc at a beam energy of 80 MeV/nucleon and an impact parameter of 𝑏 = 1.9 fm is presented. The dash-dotted (orange) line represents the source obtained from the BUU simulation using a stiff momentum-dependent EOS. The dashed line corresponds to the source function obtained using Eq. (6.7), and the solid line shows a function, 𝐶 (𝑟/𝐵)2 exp(−𝑟/𝐵), with 𝐶 = 0.63 × 10−4 fm−3 and 𝐵 = 5 fm. This function is incorporated into the tail of the source function from the BUU simulation to account for the contribution from slow emissions (i.e., the source function for particles resulting from secondary decay, which the BUU simulation does not produce). where C and B are parameters. This tail represents the source function of emissions from secondary decays, which the BUU code does not provide. In Fig. 6.10, we have plotted the P-P source function obtained from the BUU model simulation for the Ar+Sc reaction at E/A=80 MeV. We have chosen to illustrate the source corresponding to a stiff momentum-dependent EOS as an example. The solid line represents the tail function, 𝐶 (𝑟/𝐵)2 exp(−𝑟/𝐵) with the parameters 𝐶 = 0.63 × 10−4 fm−3 and 𝐵 = 5fm (we will use these parameters in the rest of the section). The dash-dotted line shows the BUU source function (𝑆𝐵𝑈𝑈 (𝑟) as defined in Eq. (6.5)), while the dashed line represents the source function obtained from Eq. (6.7). We utilize the source function defined in Eq. (6.7) to generate the correlation functions displayed as lines in Fig. 6.11. Various line styles correspond to results derived from different equations of state (as specified in the legends). Triangles within Fig. 6.11 represent experimental data. In the same figure, we also present results obtained from the deblurring method, extensively discussed in Refs. [50, 49], represented by the blue band. It’s worth noting that the experimental data falls 94 0510152025r [fm]0.00.51.01.52.0S(r)[103fm3]stiff mo-dep EOS +tailstiff mo-dep EOStail010200.000.050.10 within this band. Furthermore, the correlation obtained using the corrected BUU source function from Eq. (6.7) fits the data well. Figure 6.11 The 𝑝 − 𝑝 correlation function is plotted as a function of relative momentum (𝑞). The triangles represent data from Ref. [14] for 36Ar +45 Sc at 80 MeV/nucleon, with a 200-400 MeV/c cut in the pair’s total momentum. Different lines show correlations obtained using the source function defined in Eq. (6.7). Experimental data, extracted from Ref. [14], are indicated by the triangles. The blue band represents the correlation function obtained using the RL algorithm (deblurring) source function. Overall, both methods, deblurring and the corrected BUU source function, reproduce the data. In panel (a), the solid line depicts outcomes for soft momentum-dependent equations of state (EoS), while the dash-dotted line represents soft momentum-independent EoS results. In panel (b), solid lines correspond to stiff momentum-dependent EoS, and dash-dotted lines show stiff momentum-independent EoS results. Panel (c) displays the correlation results: a solid line for stiff momentum-dependent EoS and a dash-dotted line for soft momentum-dependent EoS. Finally, 95 0.250.500.751.001.251.50R(q)+1(a)data soft mo-depsoft mo-indepRL(b)stiff mo-depstiff mo-indep020406080q (MeV/c) 0.250.500.751.001.251.50R(q)+1(c)stiff mo-depsoft mo-dep020406080q (MeV/c) (d)soft mo-indepstiff mo-indep panel (d) exhibits solid lines for soft momentum-independent EoS and dash-dotted lines for stiff momentum-independent EoS. We anticipate that expanding the correlation function into angular components may yield similar results. The subsequent section will delve into the angular component expansion of the source function, but applying this source to calculate p-p correlations remains a topic for future work. Figure 6.12 𝑝-𝑝 source function for36Ar +45 Sc at a beam energy of 80 MeV/nucleon and an impact parameter of b=1.9 fm. The solid line, and the blue bands correspond to the source restored using the RL algorithm. The different lines represent source functions from the BUU simulation corrected as defined in (6.7). These sources were used to produce the correlation function displayed in Fig. 6.11. The source functions corresponding to the corrected correlations are illustrated in Fig. 6.12. The solid line superposed with stars represents the source function obtained using the deblurring method (RL algorithm), with the dark and light blue band denoting the 𝜎 and 2𝜎 uncertainties associated with error re-sampling [50]. The different lines of various colors correspond to sources obtained from BUU simulations for different equations of state, as indicated in the legend of the figure. 96 05101520r [fm] 0.00.51.01.52.0S(r) [103fm3] soft mom-indep. stiff mom-indep. soft mom-dep. STIFF mom-dep.RL restoration2 uncertainty uncertainty 6.3 Quadrupole moment of Source function Expanding the source function into angular components is important because it helps to examine the source function from different directions relative to the beam direction. This source function can be expressed by summing over all angular components: 𝑆(𝑟) = ∑︁ ∑︁ 𝑙=0 𝑚 𝑆𝑙𝑚 (𝑟). (6.8) Here, 𝑙 represents the orbital angular momentum and 𝑚 represents its projection (−𝑙 ≤ 𝑚 ≤ 𝑙). In low-energy reactions, where the anisotropies are weak, only a few terms in the expansion are sufficient. In the case of a pair of identical particles, even values of 𝑙 (𝑙 = 0, 2, 4, ...) are significant due to symmetrization in the emission of particles. The moments of 𝑆(𝑟) can be obtained using the following expression: 𝑆𝑙𝑚 (𝑟) = ∫ 1 √ 4𝜋 𝑑Ω 𝑌 ∗ 𝑙𝑚 (Ω)𝑆(r). (6.9) If the function, 𝑆(𝑟) is expanded in Taylor expansion around 𝑟 = 0, the lowest order terms of the Taylor expansion that contribute to the coefficients of 𝑆 will involve 𝑙-th order derivatives and increase as 𝑟 𝑙 [105]. The projection 𝑙𝑚 in 𝑆 originates from the spherical harmonic components 𝑌𝑙𝑚, meaning that 𝑆𝑙𝑚 is related to 𝑌𝑙𝑚 as 𝑆𝑙𝑚 ∝ 𝐴(𝑟)𝑌𝑙𝑚, where 𝐴(𝑟) represents the relative distance-dependent factor of the source function. The spherical harmonic is a function of spherical coordinates Ω, where Ω = (𝜃, 𝜙), and 𝜃 and 𝜙 represent the polar and azimuth angles, respectively. We can express 𝑌𝑙𝑚 in terms of unnormalized spherical harmonics denoted as 𝑌 𝑚 𝑙 [106]: 𝑌𝑙𝑚 = (−1)𝑚 (cid:104) (2𝑙 + 1) (𝑙 − 𝑚)! 4𝜋(𝑙 + 𝑚)! (cid:105) 1/2 𝑌 𝑚 𝑙 (𝜃, 𝜙). (6.10) Here, 𝑌 𝑚 𝑙 (𝜃, 𝜙) = 𝑃𝑚 𝑙 (cos 𝜃)𝑒𝑖𝑚𝜙, where 𝑃𝑚 𝑙 represents the associated Legendre function. For simplicity, we will consider the cases where 𝑚 = 0 and examine the cases for 𝑙 = 0 and 𝑙 = 2, also known as monopole and quadrupole, respectively. For 𝑙 = 0, we have 𝑌00 = 1, which corresponds to the 𝑆00 component of the source function, equivalent to the angle-averaged source function. The angle-averaged source functions are illustrated in Fig. 6.8 and 6.7 for the 36Xe+45Sc and 129Xe+197Au reactions at 𝐸/𝐴 = 80 MeV, respectively. However, as shown in Fig. 6.8 for 97 the 36Xe+45Sc reaction, this component does not show a significant difference between equations of state (i.e., stiff momentum-dependent and momentum-independent EOS, and soft momentum- dependent and momentum-independent EOS). For the 129Xe+197Au reactions, we can see some differences between the soft momentum-dependent and momentum-independent EOS in the pair total momentum spanning 200 − 400 MeV/c. However, in the same total momentum range, the stiff EOSs results are indistinguishable. On the other hand, for high total momentum spanning 600 − 800 MeV/c, we can see a difference between the stiff momentum-dependent and momentum- independent EOS results, similarly for the soft momentum-dependent and momentum-independent EOS results in the low relative distance regime. Figure 6.13 The 𝑝 − 𝑝 𝑆20 source function is plotted for a 0-400 MeV/c cut in the total momentum of the pair using BUU simulation for 36Xe+45Sc at 80 MeV/nucleon and 𝑏 =1.9 fm. The top row displays stiff and soft cases, with the solid line representing momentum-dependent EoSs and the dots representing no momentum-dependent EoSs. The bottom row compares soft (dots) and stiff (solid line), with the left panel showing momentum-dependent EoSs and the right panel showing no momentum-dependent EoSs. 98 0.000.020.040.06S20(r)[103]StiffStiff_momSoftSoft_mom051015r (fm)0.000.020.040.06S20(r)[103]Stiff-momSoft-mom051015r (fm)StiffSoft Figure 6.14 Similar to Fig. 6.13 but for129Xe+197Au reaction at 80 MeV/nucleon with 𝑏 =3.1 fm. We extend the calculation to examine 𝑆20, which is related to the spherical harmonic component (cid:0)3(cos 𝜃)2 − 1(cid:1), where 𝜃 is the angle between the pair’s total momentum and the 𝑌20 = 1 4 √︃ 5 𝜋 relative distance between the two particles. For the two aforementioned reactions, we utilized BUU simulation to calculate 𝑆20 and analyzed the cuts in total momentum as well as other momentum directions relative to the beam direction (i.e., transverse momentum (𝑃𝑇 ), momentum in the X- direction (𝑃𝑥), and momentum in the Y-direction (𝑃𝑦), See appendix. B). In this section, We discussed quadrupole components of the proton-proton relative source function, 𝑆20 for different equations of state (EoSs). Fig. 6.13 displays quadrupole (i.e., l=2 for m=0) moment of p-p source function for 36Xe+45Sc at a beam energy of 80 MeV/nucleon with total momentum of the pair spanning in 0-400 MeV/c. The top left panel represents stiff momentum- dependent (dots) and no momentum-dependent (solid line) EoSs and the top right displays soft momentum-dependent (dots) and no momentum-dependent(solid line) EoSs. In the bottom row, we compared stiff and soft EoSs, where the left bottom shows stiff (solid line) and soft (dots) both momentum-dependent, while the bottom right panel represents stiff (solid line) and soft (dots) both no momentum-dependent. Fig. 6.14 is similar to 6.13 but for 129Xe+197Au reaction. 99 0.000.020.04S20(r)[103]StiffStiff_momSoftSoft_mom05101520r (fm)0.000.020.04S20(r)[103]Stiff-momSoft-mom05101520r (fm)StiffSoft 6.4 Summary and outlook In this chapter, we discussed proton-proton the source function for different symmetric matter equations of state using BUU transport model. We investigated both momentum-dependent and non-momentum-dependent EoSs and found that, in the case of the 36Xe+45Sc reaction, there was no significant difference in the angle-averaged source function between the two EoSs. By utilizing this source function, we proceeded to calculate the p-p correlation function and compared it with the correlation function obtained through the deblurring method. Remarkably, the deblurring method accurately reproduced the measured correlations, while the BUU correlations failed to do so at low relative momentum (𝑞 < 20 MeV/c). This is not surprising as deblurring doesn’t rely on assumptions about dynamics and is obligated to reproduce the data, unlike BUU. To gain further insights, we expanded the source functions into angular components in both reactions. We specifically focused on the quadripole components, as the angle-averaged source function is equivalent to the monopole component of the source. Through this analysis, we discerned a clear distinction between momentum-dependent and momentum-independent soft EoSs, as well as between momentum-dependent stiff and soft EoSs. Moving forward, it would be intriguing to calculate the correlation function using expanded source functions to investigate if the disparity between the BUU model’s correlation and the data at low q values is linked to the higher-order angular components of the source functions. This analysis should be applied to reaction systems with both higher and intermediate masses, offering an intriguing understanding of the nuclear equation of states. 100 CHAPTER 7 CONCLUSION AND OUTLOOK 7.1 Conclusion In this thesis, we studied the correlations between particles in low-energy heavy-ion collision to learn about nuclear systems at final state reactions. These correlations are an important tool to study the space-time information of the system during a collision. We developed frameworks that allowed us to access that information. In Chapter 2, we discussed two-particle correlation in heavy-ion collisions at low relative velocities and demonstrate that to have access to space-time information at the final stages of reactions. one needs to construct the scattering wave function for a pair at low relative velocities, by solving the Schrödinger equation with potential constructed to match measured phase shifts for the system. Using the wavefunction we interpret available data on 𝛼–𝛼 correlation in terms of emitting sources, to our knowledge for the first time in the literature. Besides incoherent pair emission, the data contain the effect of 9Be resonances decaying into the 𝛼 pairs. Information on emission geometry can be enhanced by combining pairs of particles, ultimately considering matrices of correlations and sources. As an example involving 𝛼 particles, we consider the possibility of simultaneous measurements OF 𝛼 – 𝛼, d – 𝛼, and d – d correlations. In Chapter 3, we developed a novel deblurring method, that uses the inverse Richardson-Lucy algorithm which was originally developed in optics for optical image restoration. We demonstrated that the method can be used for two-particle source function imaging from correlation function. The method only needs the kernel of the pair and measured correlation as inputs. We applied the method on d–𝛼 correlation measured in low energy heavy ion collision and showed its success on imaging source function for a positively defined kernel and correlations. Additionally, in Chapter 4, we utilized the deblurring method to restore decay energy spectra for the 26O case. Our technique enables accessing information about the shell structure of particle-unbound systems from the measured decay energy spectrum, which is not readily attainable using traditional approaches such as chi-square fitting. 101 Furthermore, in Chapter 5 we were interested in machine learning techniques where we devel- oped a deep neural network (DNN) classifier as a machine learning model to identify resonance states from the measured decay energy spectrum of 26O →24 O + n + n, obtained through invariant mass spectroscopy. The findings from the deblurring and DNN classifier indicated tree resonance peaks in the measured decay energy spectrum of 26O. Finally, in Chapter 6, we utilized the BUU transport model to simulate the proton-proton source function. We compared this source with the one obtained through the deblurring technique for the p-p correlation function measured in 36Ar+45Sc at 80 A MeV. Our objective was to examine the sensitivity of the source function to the nuclear equation of state (EOS). However, we did not observe significant sensitivity of the angle- averaged p-p source function to the EOS. Additionally, this source function failed to reproduce the low relative momentum in the measured correlation. This discrepancy may be attributed to the insufficient constraint on the calculation of the correlation function’s slow emission contributions. On the other hand, the deblurring method effectively reproduces the measured p-p correlation function. Furthermore, in search of the sensitivity of the source function on EOS, we expand the source function into angular components where we explored the sensitivity of quadrupole moment of the function and as shown in Figs. 6.13 and 6.14 we can observe the difference in EOSs but the corresponding correlations have to be yet determined. 7.2 Outlook The present thesis has yielded novel insights that lay the foundation for future investigations. Our primary focus has been on angle-averaged two-particle correlations, encompassing the angle- averaged source function and wave function. However, expanding the source function to incorporate angular components and examining the system from various orientations could provide enhanced understanding of the deformed source function and its impact on particle emission in the final state. Thus far, we have successfully employed the deblurring method to address one-dimensional problems, achieving satisfactory outcomes. Our forthcoming endeavors involve extending this method to three-dimensional problems. For instance, we intend to apply the method to a three- dimensional source function, enabling us to access the requisite 3D information and infer 3D 102 correlation functions. Furthermore, in relation to the BUU transport model simulation, we have investigated the quadrupole moment of the p-p source function. The subsequent step entails calculating the cor- responding correlation function for this source. To this end, we will initially determine the p-p quadrupole kernel and subsequently employ it in conjunction with the quadrupole source to de- duce the correlation. The total correlation function will then be obtained as the summation of all correlation components, encompassing the monopole and quadrupole. Additionally, we possess a keen interest in expanding our proficiency in machine learning to construct a Deep Neural Network (DNN) model for solving inverse problems. We plan to test this model by utilizing it to compute two-particle correlations and comparing its results with those obtained from the deblurring method. 103 BIBLIOGRAPHY [1] [2] [3] G. Bertsch and S. Das Gupta, “A guide to microscopic models for intermediate energy heavy ion collisions,” Physics Reports, vol. 160, no. 4, pp. 189–233, 1988. URL: https://www.sciencedirect.com/science/article/pii/0370157388901706 N. Herrmann, J. P. Wessels, and T. Wienold, “Collective flow in heavy-ion collisions,” Annual Review of Nuclear and Particle Science, vol. 49, no. 1, pp. 581–632, 1999. URL: https://doi.org/10.1146/annurev.nucl.49.1.581 “Probing heavy ion collisions by using d-𝛼 correlation,” accessed: 2023-06-25. URL: https://groups.nscl.msu.edu/hira/projects/chan/chan05.htm [4] W. G. Gong, W. Bauer, C. K. Gelbke, N. Carlin, R. T. de Souza, Y. D. Kim, W. G. Lynch, T. Murakami, G. Poggi, D. P. Sanderson, M. B. Tsang, H. M. Xu, S. Pratt, D. E. Fields, K. Kwiatkowski, R. Płaneta, V. E. Viola, and S. J. Yennello, “Intensity-interferometric test of nuclear collision geometries obtained from the boltzmann-uehling-uhlenbeck equation,” Phys. Rev. Lett., vol. 65, pp. 2114–2117, Oct 1990. URL: https://link.aps.org/doi/10.1103/PhysRevLett.65.2114 [5] [6] [7] [8] S. E. Koonin, “Proton pictures of high-energy nuclear collisions,” Physics Letters B, vol. 70, no. 1, pp. 43–47, 1977. P. A. DeYoung, M. S. Gordon, X. q. Lu, R. L. McGrath, J. M. Alexander, D. M. de Castro Rizzo, and L. C. Vaz, “Emission times in low-energy heavy-ion reactions by particle-particle correlations,” Phys. Rev. C, vol. 39, pp. 128–131, Jan 1989. URL: https://link.aps.org/doi/10.1103/PhysRevC.39.128 R. Ghetti, E. Pollacco, J. Helgesson, R. Dayras, S. Cavallaro, A. Pagano, M. Geraci, E. De Fil- ippo, C. Volant, J. Charvet et al., “Correlation functions and emission time sequence of light charged particles from projectile-like fragment source in 𝐸/𝐴 = 44 MeV and 77 MeV 40Ar+ 27Al collisions,” Nucl. Phys., vol. 765, no. nucl-ex/0412038, pp. 307–318, 2004. R. Ghetti, J. Helgesson, G. Lanzanò, E. De Filippo, M. Geraci, S. Aiello, S. Cavallaro, A. Pagano, G. Politi, J. Charvet, R. Dayras, E. Pollacco, C. Volant, C. Beck, D. Mahboub, and R. Nouicer, “Correlation functions of light charged particles from projectile-like fragment source in 𝐸/𝐴 = 44 and77 MeV 40Ar + 27Al collisions,” Nuclear Physics A, vol. 765, no. 3, pp. 307–318, 2006. URL: https://www.sciencedirect.com/science/article/pii/S0375947405012224 [9] G. Verde, A. Chbihi, R. Ghetti, and J. Helgesson, “Correlations and characterization of emitting sources,” The European Physical Journal A-Hadrons and Nuclei, vol. 30, no. 1, pp. 81–108, 2006. [10] D. H. Boal, C.-K. Gelbke, and B. K. Jennings, subatomic physics,” Rev. Mod. Phys., vol. 62, pp. 553–602, https://link.aps.org/doi/10.1103/RevModPhys.62.553 “Intensity interferometry in Jul 1990. URL: 104 [11] W. Bauer, “Particle-particle correlations and the space-time structure of heavy ion collisions,” Progress in Particle and Nuclear Physics, vol. 30, pp. 45–64, 1993. [12] U. Heinz, B. Tomášik, U. A. Wiedemann, and Y. F. Wu, “Lifetimes and sizes from two-particle correlation functions,” Physics Letters B, vol. 382, no. 1, pp. 181–188, 1996. URL: https://www.sciencedirect.com/science/article/pii/0370269396006570 [13] S. Pratt, “Pion interferometry for exploding sources,” Physical Review Letters, vol. 53, no. 13, p. 1219, 1984. [14] G. Verde, P. Danielewicz, D. A. Brown, W. G. Lynch, C. K. Gelbke, and M. B. Tsang, “Probing transport theories via two-proton source imaging,” Phys. Rev. C, vol. 67, p. 034606, Mar 2003. URL: https://link.aps.org/doi/10.1103/PhysRevC.67.034606 [15] H. Wolter, M. Colonna, D. Cozma, P. Danielewicz, C. M. Ko, R. Kumar, A. Ono, M. B. Tsang, J. Xu, Y.-X. Zhang, E. Bratkovskaya, Z.-Q. Feng, T. Gaitanos, A. Le Fèvre, N. Ikeno, Y. Kim, S. Mallik, P. Napolitani, D. Oliinychenko, T. Ogawa, M. Papa, J. Su, R. Wang, Y.-J. Wang, J. Weil, F.-S. Zhang, G.-Q. Zhang, Z. Zhang, J. Aichelin, W. Cassing, L.-W. Chen, H.-G. Cheng, H. Elfner, K. Gallmeister, C. Hartnack, S. Hashimoto, S. Jeon, K. Kim, M. Kim, B.-A. Li, C.-H. Lee, Q.-F. Li, Z.-X. Li, U. Mosel, Y. Nara, K. Niita, A. Ohnishi, T. Sato, T. Song, A. Sorensen, N. Wang, and W.-J. Xie, “Transport model comparison studies of intermediate-energy heavy-ion collisions,” Progress in Particle and Nuclear Physics, vol. 125, p. 103962, 2022. URL: https://www.sciencedirect.com/science/article/pii/S0146641022000230 [16] R. H. Brown, R. Jennison, and M. Gupta, “Apparent angular sizes of discrete radio sources: Observations at Jodrell Bank, Manchester.” Nature, vol. 170, no. 4338, pp. 1061–1063, 1952. [17] R. H. Brown and R. Q. Twiss, 2. A Test of a New Type of Stellar Interferometer on Sirius. Cambridge, MA and London, England: Harvard University Press, 2013, pp. 8–12. ISBN 9780674366688. URL: https://doi.org/10.4159/harvard.9780674366688.c4 [18] D. G. C. Jones, “Two slit interference-classical and quantum pictures,” European Journal of Physics, vol. 15, no. 4, p. 170, Jul. 1994. URL: https://dx.doi.org/10.1088/0143-0807/15/ 4/003 [19] D. H. Boal and J. C. Shillcock, “Nuclear interferometry and thermal freeze out,” Phys. Rev. C, vol. 33, pp. 549–556, Feb 1986. URL: https://link.aps.org/doi/10.1103/PhysRevC.33.549 [20] P. Samuelsson, I. Neder, and M. Büttiker, “Entanglement at finite temperatures in the electronic two-particle interferometer,” Physica Scripta, vol. 2009, no. T137, p. 014023, dec 2009. URL: https://dx.doi.org/10.1088/0031-8949/2009/T137/014023 [21] C. Ezell, L. J. Gutay, A. T. Laasanen, F. T. Dao, P. Schübelin, and F. Turkot, “Determination of the space-time extension of centrally produced hadronic matter using the intensity interferometer technique in pp interactions at 28.5 GeV/c,” Phys. Rev. Lett., vol. 38, pp. 873–876, Apr 1977. URL: https://link.aps.org/doi/10.1103/PhysRevLett.38.873 105 [22] M. Abramowitz and I. A. Stegun, Handbook of mathematical functions with formulas, graphs, and mathematical tables. US Government Printing Office, 1948, vol. 55. [23] I. J. Thompson and F. M. Nunes, Nuclear Reactions for Astrophysics: Principles, Calculation and Applications of Low-Energy Reactions. Cambridge University Press, 2009. [24] L.-W. Chen, V. Greco, C. M. Ko, symmetry induced by energy on two-nucleon correlation functions neutron-rich nuclei,” Phys. Rev. Lett., vol. 90, p. 162701, Apr 2003. URL: https://link.aps.org/doi/10.1103/PhysRevLett.90.162701 in heavy-ion collisions and B.-A. Li, “Effects of [25] G. F. Bertsch, “Meson phase space density in heavy ion collisions from interferometrfoy,” Phys. Rev. Lett., vol. 72, no. 15, pp. 2349–2350, Apr. 1994, corrected in [79]. URL: https://link.aps.org/doi/10.1103/PhysRevLett.72.2349 [26] S. Chapman, P. Scotto, and U. Heinz, “Model independent features of the two-particle correlation function,” Acta Physica Hungarica New Series Heavy Ion Physics, vol. 1, no. 1, pp. 1–31, 1995. [27] D. A. Brown and P. Danielewicz, “Imaging of sources in heavy-ion reactions,” Physics Letters B, vol. 398, no. 3-4, pp. 252–258, 1997. [28] M. A. Lisa, S. Pratt, R. Soltz, and U. Wiedemann, “Femtoscopy in relativistic heavy ion collisions: two decades of progress,” Annu. Rev. Nucl. Part. Sci., vol. 55, pp. 357–402, 2005. [29] R. Ghetti, J. Helgesson, G. Lanzanò, E. De Filippo, and M. Geraci, “Two-particle correlation functions in projectilelike emission experiments at very forward angles,” Phys. Rev. C, vol. 70, p. 027601, Aug 2004. URL: https://link.aps.org/doi/10.1103/PhysRevC.70.027601 [30] T. K. Nayak, T. Murakami, W. G. Lynch, K. Swartz, D. J. Fields, C. K. Gelbke, Y. D. Kim, J. Pochodzalla, M. B. Tsang, H. M. Xu, F. Zhu, and K. Kwiatkowski, “Emission temperatures from the decay of particle unstable complex nuclei,” Phys. Rev. C, vol. 45, pp. 132–161, Jan 1992. URL: https://link.aps.org/doi/10.1103/PhysRevC.45.132 [31] D. Anchishkin, U. Heinz, and P. Renk, “Final state interactions in two-particle pp. 1428–1439, Mar 1998. URL: interferometry,” Phys. Rev. C, vol. 57, https://link.aps.org/doi/10.1103/PhysRevC.57.1428 [32] D. H. Boal, C.-K. Gelbke, and B. K. Jennings, subatomic physics,” Rev. Mod. Phys., vol. 62, pp. 553–602, https://link.aps.org/doi/10.1103/RevModPhys.62.553 “Intensity interferometry in Jul 1990. URL: [33] S. Pratt and M. B. Tsang, “Viewing the liquid-gas phase transition by measuring proton correlations,” Phys. Rev. C, vol. 36, pp. 2390–2395, Dec 1987. URL: https://link.aps.org/doi/10.1103/PhysRevC.36.2390 [34] P. Danielewicz and P. Schuck, “Formulation of particle correlation and cluster production in heavy-ion-induced reactions,” Physics Letters B, vol. 274, no. 3, pp. 268–274, Jan. 1992. URL: http://www.sciencedirect.com/science/article/pii/037026939291985I 106 [35] R. Lednicky and V. L. Lyuboshits, “Final State Interaction Effect on Pairing Correlations Between Particles with Small Relative Momenta,” Yad. Fiz., vol. 35, p. 1316, 1981. [36] S. Mrówczyński, “On the neutron—proton correlations and deuteron production,” Feb. 1992. URL: https: Physics Letters B, vol. 277, //www.sciencedirect.com/science/article/pii/0370269392909543 pp. 43–48, no. 1, [37] S. A. AFZAL, A. A. Z. AHMAD, and S. ALI, “Systematic survey of 𝛼 − 𝛼 interaction,” Rev. Mod. Phys., vol. 41, pp. 247–273, https://link.aps.org/doi/10.1103/RevModPhys.41.247 the Jan 1969. URL: [38] P. Danielewicz, P. Singh, and J. Lee, “Symmetry energy III: Isovector skins,” Nucl. Phys. A, vol. 958, pp. 147–186, Feb. 2017. URL: http://www.sciencedirect.com/science/article/pii/ S0375947416302895 [39] L. McIntyre and W. Haeberli, “Phase shift analysis of d-𝛼 scattering,” Nuclear Physics A, vol. 91, no. 2, pp. 382–398, 1967. [40] W. G. Gong, P. Danielewicz, C. K. Gelbke, N. Carlin, R. T. de Souza, Y. D. Kim, W. G. Lynch, T. Murakami, G. Poggi, M. B. Tsang, H. M. Xu, S. Pratt, K. Kwiatkowski, V. E. Viola, S. J. Yennello, and J. C. Shillcock, “Two-deuteron correlation functions in 14N+27al collisions at e/a=75 mev,” Phys. Rev. C, vol. 47, pp. R429–R432, Feb 1993. URL: https://link.aps.org/doi/10.1103/PhysRevC.47.R429 [41] D. L. Mihaylov, V. Mantovani Sarti, O. W. Arnold, L. Fabbietti, B. Hohlweger, and A. M. Mathis, “A femtoscopic correlation analysis tool using the schrödinger equation (cats),” The European Physical Journal C, vol. 78, no. 5, p. 394, 2018. [42] Z. Chen, C. K. Gelbke, W. G. Gong, Y. D. Kim, W. G. Lynch, M. R. Maier, J. Pochodzalla, M. B. Tsang, F. Saint-Laurent, D. Ardouin, H. Delagrange, H. Doubre, J. Kasagi, A. Kyanowski, A. Péghaire, J. Péter, E. Rosato, G. Bizard, F. Lefèbvres, B. Tamain, J. Québert, and Y. P. Viyogi, “Inclusive two-particle correlations for 16induced reactions on 197Au at E/A=94 MeV,” Phys. Rev. C, vol. 36, pp. 2297–2308, Dec 1987. URL: https://link.aps.org/doi/10.1103/PhysRevC.36.2297 [43] P. E. Shanley, “Three-body model of 6Li and deuteron-alpha-particle scattering,” Physical Review, vol. 187, no. 4, p. 1328, 1969. [44] Z. Fan, “Two particles correlation and excited state population in nucleus -nucleus collision,” Ph.D. dissertation, Michigan State University/Department of Physics and Astronomy, 1992. [45] T. Brown, P. Papka, B. Fulton, D. Watson, S. Fox, D. Groombridge, M. Freer, N. Clarke, N. Ashwood, N. Curtis et al., “Decay studies for states in be 9 up to 11 mev: Insights into the n+ be 8 and 𝛼+ he 5 cluster structure,” Physical Review C, vol. 76, no. 5, p. 054605, 2007. [46] J. I. Kapusta, “Particle production in the nuclear fireball model,” Physical Review C, vol. 16, no. 4, p. 1493, 1977. 107 [47] G. Verde, P. Danielewicz, W. Lynch, C. Chan, C. Gelbke, L. Kwong, T. Liu, X. Liu, D. Seymour, R. Shomin et al., “d–𝛼 correlation functions and collective motion in xe+ au collisions at e/a= 50 mev,” Physics Letters B, vol. 653, no. 1, pp. 12–17, 2007. [48] P. Danielewicz and S. Pratt, “Analyzing correlation functions with tesseral and cartesian spherical harmonics,” Physical Review C, vol. 75, no. 3, p. 034907, 2007. [49] P. Nzabahimana and P. Danielewicz, “Source function from two-particle correlation 2023. URL: through deblurring,” Physics Letters B, vol. 846, https://www.sciencedirect.com/science/article/pii/S0370269323005816 p. 138247, [50] P. Nzabahimana, T. Redpath, T. Baumann, P. Danielewicz, P. Giuliani, and P. Guèye, “Deconvoluting experimental decay energy spectra: The 26O case,” Phys. Rev. C, vol. 107, p. 064315, Jun 2023. URL: https://link.aps.org/doi/10.1103/PhysRevC.107.064315 [51] P. Danielewicz and M. Kurata-Nishimura, “Deblurring for nuclei: 3d characteristics of heavy-ion collisions,” Phys. Rev. C, vol. 105, no. 3, p. 034608, 2022. [52] W. H. Richardson, “Bayesian-Based Iterative Method of Image Restoration*,” Journal the Optical Society of America, vol. 62, no. 1, pp. 55–59, Jan. 1972. URL: of http://www.osapublishing.org/josa/abstract.cfm?uri=josa-62-1-55 [53] L. B. Lucy, “An iterative technique for the rectification of observed distributions,” The Astronomical Journal, vol. 79, p. 745, Jun. 1974. URL: https://ui.adsabs.harvard.edu/abs/ 1974AJ.....79..745L [54] D. A. Brown and P. Danielewicz, “Imaging of sources in heavy-ion reactions,” Physics Letters B, vol. 398, no. 3, pp. 252–258, 1997. URL: https://www.sciencedirect.com/science/ article/pii/S0370269397002517 [55] P. Danielewicz, P. Singh, and J. Lee, “Symmetry energy iii: Isovector skins,” Nuclear Physics A, vol. 958, pp. 147–186, 2017. URL: https://www.sciencedirect.com/science/ article/pii/S0375947416302895 [56] G. Zech, “Iterative unfolding with algorithm,” Nuclear the Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 2013. URL: https: //www.sciencedirect.com/science/article/pii/S0168900213003148 richardson–lucy vol. 716, pp. 1–9, [57] P. Nzabahimana, T. Redpath, T. Baumann, P. Danielewicz, P. Giuliani, and P. Guèye, “Deblurring decay-energy distribution from an invariant-mass measurement,” Journal of Physics: Conference Series, vol. 2586, no. 1, p. 012042, sep 2023. URL: https://dx.doi.org/10.1088/1742-6596/2586/1/012042 [58] G. Lanzanó, A. Pagano, S. Urso, E. De Filippo, B. Berthier, J. L. Charvet, R. Dayras, R. Legrain, R. Lucas, C. Mazur, E. Pollacco, J. E. Sauvestre, C. Volant, C. Beck, B. Djerroud, and B. Heusch, “Using BaF2 crystals as detectors of light charged particles at intermediate energies,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, 108 Spectrometers, Detectors and Associated Equipment, vol. 312, no. 3, pp. 515–520, Feb. 1992. URL: https://www.sciencedirect.com/science/article/pii/016890029290199E [59] T. Redpath, “Measuring dissertation, ISBN Michigan 9781392603222 2019. URL: http://ezproxy.msu.edu/login?url=https://www.proquest.com/ dissertations-theses/measuring-half-life-o-26/docview/2345890497/se-2 Ph.D. and Astronomy. State University/Department half-life of of Physics o-26,” the [60] C. Caesar, J. Simonis, T. Adachi, Y. Aksyutina, J. Alcantara, S. Altstadt, H. Alvarez-Pol, N. Ashwood, T. Aumann, V. Avdeichikov, M. Barr et al., “Beyond the neutron drip line: The unbound oxygen isotopes 25o and 26o,” Phys. Rev. C, vol. 88, p. 034313, Sep 2013. URL: https://link.aps.org/doi/10.1103/PhysRevC.88.034313 [61] D. Chrisman, A. N. Kuchera, T. Baumann, A. Blake, B. A. Brown, J. Brown, C. Cochran, P. A. DeYoung, J. E. Finck, N. Frank, P. Guèye, H. Karrick, H. Liu, J. McDonaugh, T. Mix, B. Monteagudo, T. H. Redpath, W. F. Rogers, R. Seaton-Todd, A. Spyrou, K. Stiefel, M. Thoennessen, J. A. Tostevin, and D. Votaw, “Neutron-unbound states in 31Ne,” Phys. Rev. C, vol. 104, p. 034313, Sep 2021. URL: https://link.aps.org/doi/10.1103/PhysRevC.104.034313 [62] A. Revel, O. Sorlin, F. M. Marqués, Y. Kondo, J. Kahlbow, T. Nakamura, N. A. Orr, F. Nowacki, J. A. Tostevin, C. X. Yuan, N. L. Achouri, H. Al Falou, L. Atar, T. Aumann, H. Baba, K. Boretzky, C. Caesar, D. Calvet, H. Chae, N. Chiga, A. Corsi, H. L. Crawford, F. Delaunay, A. Delbart, Q. Deshayes, Z. Dombrádi, C. A. Douma, Z. Elekes, P. Fallon, I. Gašparić, J.-M. Gheller, J. Gibelin, A. Gillibert, M. N. Harakeh, W. He, A. Hirayama, C. R. Hoffman, M. Holl, A. Horvat, A. Horváth, J. W. Hwang, T. Isobe, N. Kalantar-Nayestanaki, S. Kawase, S. Kim, K. Kisamori, T. Kobayashi, D. Körper, S. Koyama, I. Kuti, V. Lapoux, S. Lindberg, S. Masuoka, J. Mayer, K. Miki, T. Murakami, M. Najafi, K. Nakano, N. Nakatsuka, T. Nilsson, A. Obertelli, F. de Oliveira Santos, H. Otsu, T. Ozaki, V. Panin, S. Paschalis, D. Rossi, A. T. Saito, T. Saito, M. Sasano, H. Sato, Y. Satou, H. Scheit, F. Schindler, P. Schrock, M. Shikata, Y. Shimizu, H. Simon, D. Sohler, L. Stuhl, S. Takeuchi, M. Tanaka, M. Thoennessen, H. Törnqvist, Y. Togano, T. Tomai, J. Tscheuschner, J. Tsubota, T. Uesaka, Z. Yang, M. Yasuda, and K. Yoneda, “Extending the southern shore of the island of inversion to 28F,” Phys. Rev. Lett., vol. 124, p. 152502, Apr 2020. URL: https://link.aps.org/doi/10.1103/PhysRevLett.124.152502 [63] S. Leblond, F. M. Marqués, J. Gibelin, N. A. Orr, Y. Kondo, T. Nakamura, J. Bonnard, N. Michel, N. L. Achouri, T. Aumann, H. Baba, F. Delaunay, Q. Deshayes, P. Doornenbal, N. Fukuda, J. W. Hwang, N. Inabe, T. Isobe, D. Kameda, D. Kanno, S. Kim, N. Kobayashi, T. Kobayashi, T. Kubo, J. Lee, R. Minakata, T. Motobayashi, D. Murai, T. Murakami, K. Muto, T. Nakashima, N. Nakatsuka, A. Navin, S. Nishi, S. Ogoshi, H. Otsu, H. Sato, Y. Satou, Y. Shimizu, H. Suzuki, K. Takahashi, H. Takeda, S. Takeuchi, R. Tanaka, Y. Togano, A. G. Tuff, M. Vandebrouck, and K. Yoneda, “First observation of 20B and 21B,” Phys. Rev. Lett., vol. 121, p. 262502, Dec 2018. URL: https://link.aps.org/doi/10.1103/PhysRevLett.121.262502 109 [64] A. M. Lane and R. G. Thomas, “R-matrix theory of nuclear reactions,” Rev. Mod. Phys., vol. 30, pp. 257–353, Apr 1958. URL: https://link.aps.org/doi/10.1103/RevModPhys.30.257 [65] W. H. Richardson, “Bayesian-based iterative method of image restoration,” JoSA, vol. 62, no. 1, pp. 55–59, 1972. [66] L. B. Lucy, “An iterative technique for the rectification of observed distributions,” Astron. J., vol. 79, p. 745, 1974. [67] É. Thiébaut, L. Dénis, F. Soulez, and R. Mourya, “Spatially variant psf modeling and image International Society for Optics and deblurring,” in Adaptive Optics Systems V, vol. 9909. Photonics, 2016, p. 99097N. [68] Z. Al-Ameen and G. Sulong, “Deblurring computed tomography medical images using a novel amended landweber algorithm,” Interdiscipl. Sci.: Comput. Life Sci.,, vol. 7, no. 3, pp. 319–325, 2015. [69] G. D’Agostini, “A multidimensional unfolding method based on Bayes’ theorem,” Nucl. Instrum. Methods. Phys. Res. A, vol. 362, no. 2, pp. 487–498, Aug. 1995. URL: https://www.sciencedirect.com/science/article/pii/016890029500274X [70] R. Grech, T. Cassar, J. Muscat, K. P. Camilleri, S. G. Fabri, M. Zervakis, P. Xanthopoulos, V. Sakkalis, and B. Vanrumste, “Review on solving the inverse problem in eeg source analysis,” J. Neuroengineering.Rehabil., vol. 5, no. 1, pp. 1–33, 2008. URL: https://doi.org/10.1186/1743-0003-5-25 [71] T. T. Fister, G. T. Seidler, J. O. Cross, and K. P. Nagle, “Deconvolving instrumental and intrinsic broadening in core-shell x-ray spectroscopies,” Phys. Rev. B, vol. 75, p. 174106, May 2007. URL: https://link.aps.org/doi/10.1103/PhysRevB.75.174106 J. J. Kas, W. T. Elam, J. J. Rehr, [72] N. Dey, L. Blanc-Feraud, C. Zimmer, P. Roux, Z. Kam, J.-C. Olivo-Marin, and J. Zerubia, “Richardson–lucy algorithm with total variation regularization for 3d confocal microscope deconvolution,” Microsc. Res. Tech., vol. 69, no. 4, pp. 260–266, 2006. [73] J. Benlliure, and M. Caamaño, “Unfolding the response of a zero- J. Vargas, the Δ resonance,” Nucl. degree magnetic spectrometer Instrum. Methodes. Phys. Res. Sect., A:, vol. 707, pp. 16–25, Apr. 2013. URL: https://www.sciencedirect.com/science/article/pii/S016890021201635X from measurements of [74] L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algo- rithms,” Physica D, vol. 60, no. 1-4, pp. 259–268, 1992. [75] T. Redpath, T. Baumann, J. Brown, D. Chrisman, P. DeYoung, N. Frank, P. Guèye, A. Kuchera, H. Liu, C. Persch et al., “New segmented target for studies of neutron unbound systems,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 977, p. 164284, 2020. 110 [76] G. Bohm and G. Zech, Introduction to statistics and data analysis for physicists. DESY, 2010. ISBN 978-3-935702-41-6 [77] Y. Kondo, T. Nakamura, R. Tanaka, R. Minakata, S. Ogoshi, N. Orr, N. Achouri, T. Aumann, H. Baba, F. Delaunay et al., “Nucleus o 26: A barely unbound system beyond the drip line,” Physical review letters, vol. 116, no. 10, p. 102503, 2016. [78] Z. Kohley, T. Baumann, D. Bazin, G. Christian, P. A. DeYoung, J. Finck, N. Frank, M. Jones, E. Lunderberg, B. Luther et al., “Study of two-neutron radioactivity in the decay of o 26,” Physical review letters, vol. 110, no. 15, p. 152501, 2013. [79] G. F. Bertsch, from Interferometry,” Phys. Rev. Lett., vol. 77, no. 4, pp. 789–789, Jul. 1996. URL: https://link.aps.org/doi/10.1103/PhysRevLett.77.789 “Meson Phase Space Density in Heavy Ion Collisions [80] A. Géron, Hands-on machine learning with Scikit-Learn, Keras, and TensorFlow. " O’Reilly Media, Inc.", 2022. [81] S. Lawrence, C. Giles, A. C. Tsoi, and A. Back, “Face recognition: a convolutional neural- network approach,” IEEE Trans. Neural Networks, vol. 8, no. 1, pp. 98–113, 1997. doi: 10.1109/72.554195. [82] A. B. Nassif, I. Shahin, I. Attili, M. Azzeh, and K. Shaalan, “Speech recognition using deep neural networks: A systematic review,” IEEE Access, vol. 7, pp. 19 143–19 165, 2019. doi: 10.1109/ACCESS.2019.2896880. [83] D. Guest, K. Cranmer, and D. Whiteson, “Deep learning and its application to lhc physics,” Annual Review of Nuclear and Particle Science, vol. 68, pp. 161–181, 2018. [84] K. T. Matchev and P. Shyamsundar, “Thickbrick: optimal event selection and categorization in high energy physics. part i. signal discovery,” Journal of High Energy Physics, vol. 2021, no. 3, pp. 1–53, 2021. [85] G. Carleo and M. Troyer, “Solving the quantum many-body problem with artificial neural networks,” Science, vol. 355, no. 6325, pp. 602–606, 2017. [86] S. Whiteson and D. Whiteson, “Machine learning for event selection in high energy physics,” Engineering Applications of Artificial Intelligence, vol. 22, no. 8, pp. 1203–1217, 2009. [87] Y. Fujimoto, K. Fukushima, and K. Murase, “Mapping neutron star data to the equation of state using the deep neural network,” Phys. Rev. D, vol. 101, no. 5, p. 054016, 2020. [88] P. Bedaque, A. Boehnlein, M. Cromaz, M. Diefenthaler, L. Elouadrhiri, T. Horn, M. Kuchera, D. Lawrence, D. Lee, S. Lidia et al., “Ai for nuclear physics,” Eur. Phys. J. A, vol. 57, no. 3, pp. 1–27, 2021. [89] A. Boehnlein, M. Diefenthaler, N. Sato, M. Schram, V. Ziegler, C. Fanelli, M. Hjorth-Jensen, T. Horn, M. P. Kuchera, D. Lee, W. Nazarewicz, P. Ostroumov, K. Orginos, A. Poon, X.-N. Wang, A. Scheinker, M. S. Smith, and L.-G. Pang, “Colloquium: Machine 111 learning in nuclear physics,” Rev. Mod. Phys., vol. 94, p. 031003, Sep 2022. URL: https://link.aps.org/doi/10.1103/RevModPhys.94.031003 [90] J.-W. Lin, “Artificial neural network related to biological neuron network: A review,” Ad- vanced Studies in Medical Sciences, vol. 5, no. 1, pp. 55–62, 2017. [91] S. Lu, Z. Li, Z. Qin, X. Yang, and R. S. M. Goh, “A hybrid regression technique for house prices prediction,” in 2017 IEEE International Conference on Industrial Engineering and Engineering Management (IEEM), 2017. doi: 10.1109/IEEM.2017.8289904. pp. 319–323. [92] S. Dreiseitl and L. Ohno-Machado, “Logistic regression and artificial neural network a methodology review,” Journal of Biomedical Informatics, classification models: vol. 35, no. 5, pp. 352–359, 2002. URL: https://www.sciencedirect.com/science/article/pii/ S1532046403000340 [93] Y. LeCun, L. D. Jackel, L. Bottou, C. Cortes, J. S. Denker, H. Drucker, I. Guyon, U. A. Muller, E. Sackinger, P. Simard et al., “Learning algorithms for classification: A comparison on handwritten digit recognition,” Neural networks: the statistical mechanics perspective, vol. 261, no. 276, p. 2, 1995. [94] I. Daubechies, R. DeVore, S. Foucart, B. Hanin, and G. Petrova, “Nonlinear approximation and (deep) relu networks,” arXiv preprint arXiv:1905.02199, 2019. [95] H. H. Tan and K. H. Lim, “Vanishing gradient mitigation with deep learning neural network optimization,” pp. 1–4, 2019. doi: 10.1109/ICSCC.2019.8843652. [96] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014. [97] S. Sharma, S. Sharma, neural networks,” Towards Data Science, vol. 6, no. 12, pp. 310–316, 2017. URL: https://towardsdatascience.com/activation-functions-neural-networks-1cbd9f8d91d6 and A. Athaiya, “Activation functions in [98] A. Rusiecki, “Trimmed categorical cross-entropy for deep learning with label noise,” Electron. Lett., vol. 55, no. 6, pp. 319–320, 2019. URL: https://doi.org/10.1049/el.2018.7980 [99] T. Redpath, Measuring the Half-life of O-26. Michigan State University. Physics, 2019. [100] P. Danielewicz, “Nonrelativistic and relativistic landau/fokker-planck equation for arbitrary statistics,” Physica A: Statistical Mechanics and its Applications, vol. 100, no. 1, pp. 167– 182, 1980. URL: https://www.sciencedirect.com/science/article/pii/0378437180901570 [101] H. Lin and P. Danielewicz, “One-body langevin dynamics in heavy-ion collisions intermediate energies,” Phys. Rev. C, vol. 99, p. 024612, Feb 2019. URL: at https://link.aps.org/doi/10.1103/PhysRevC.99.024612 [102] P. Danielewicz, the mean-field momentum-dependence using elliptic flow,” Nuclear Physics A, vol. 673, no. 1, pp. 375–410, 2000. URL: https://www.sciencedirect.com/science/article/pii/S037594740000083X “Determination of 112 [103] D. A. Brown and P. Danielewicz, “Optimized discretization of sources imaged in heavy-ion reactions,” Physical Review C, vol. 57, no. 5, p. 2474, 1998. [104] G. Verde, D. A. Brown, P. Danielewicz, C. K. Gelbke, W. G. Lynch, and M. B. Tsang, “Imaging sources with fast and slow emission components,” Phys. Rev. C, vol. 65, p. 054609, May 2002. URL: https://link.aps.org/doi/10.1103/PhysRevC.65.054609 [105] P. Danielewicz and S. Pratt, “Analysis of low-momentum correlations with cartesian 2005. URL: harmonics,” Physics Letters B, vol. 618, https://www.sciencedirect.com/science/article/pii/S0370269305006404 pp. 60–67, no. 1, [106] J. Applequist, “Maxwell-cartesian spherical harmonics in multipole potentials and atomic orbitals,” THEORETICAL CHEMISTRY ACCOUNTS, vol. 107, no. 2, pp. 103–115, JAN 2002. [107] A. Messiah, Quantum mechanics: volume II. North-Holland Publishing Company Ams- terdam, 1962. [108] I. J. Thompson and F. M. Nunes, Nuclear reactions for astrophysics: principles, calculation and applications of low-energy reactions. Cambridge University Press, 2009. [109] A. Mekjian, “Thermodynamic model for composite-particle emission in relativistic heavy- ion collisions,” Physical Review Letters, vol. 38, no. 12, p. 640, 1977. 113 APPENDIX A DERIVATION OF SCATTERING PHASE SHIFT A.1 Phase shift In this section, the expression for the phase shift is derived. The calculated phase shift, along with experimental data, is presented in Fig. 2.3 for the 𝑑 − 𝛼 system and in Fig. 2.2 for the 𝛼 − 𝛼 system. Let’s start with the two-particle Schrödinger Equation (SE) defined as: (𝐻 − 𝐸)Ψ(q, r) = 0 where the two-particle Hamiltonian 𝐻 is given by: 𝐻 = − ℏ2 2𝑚1 ∇2 r1 − ℏ2 2𝑚2 ∇2 r2 + 𝑉ef(r1 − r2) (A.1) (A.2) Here, the effective potential 𝑉ef(r1 − r2) depends only on the relative distance r = r1 − r2. The potential 𝑉ef(r) includes the nuclear central potential, the long-range Coulomb interaction, and the centrifugal term ℏ2𝑙 (𝑙+1) 2𝜇𝑟2 . The Hamiltonian in Eq. (A.2) can be transformed and written as a separable expression: 𝐻 = − ℏ2 2𝑀 ∇2 R − ℏ2 2𝜇 ∇2 r + 𝑉ef(r) (A.3) In Eq. (A.3), the first term on the right-hand side represents the Hamiltonian of the center-of- mass motion and represents the motion of a free particle. In contrast, the second and third terms represent the Hamiltonian of the relative motion of the two particles, incorporating the interaction between them. We denote the reduced mass of the two-particle system, 𝜇 = 𝑚1𝑚2 𝑚1+𝑚2 , and total mass, 𝑀 = 𝑚1 + 𝑚2. For our purpose, we consider the relative motion party of Hamiltonian and write a radial Schr¨dinger equation as the following: 𝑑2 𝑑𝑟 2 (𝑢) = 2𝜇 ℏ2 (cid:18) 𝑉 (𝑟) + ℏ2𝑙 (𝑙 + 1) 2𝜇𝑟 2 (cid:19) 𝑢, − 𝐸 (A.4) where 𝑢 is the radial wave function, 𝑉 (𝑟) = 𝑉𝑛 (𝑟) + 𝑉𝑐 (𝑟) is the potential, and 𝑉𝑛 and 𝑉𝑐 represent the nuclear and Coulomb potentials, respectively. In our case, the nuclear potential is modeled by 114 the Wood-Saxon potential of the form: 1 1 + 𝑒 where 𝑅0 is the potential range and 𝑎 is the diffuseness. The nuclear potential, 𝑉𝑛 may contain one 𝑉𝑛 (𝑟) ∝ (A.5) 𝑟 −𝑅 𝑎 , 0 than one term, be real or complex (i.e., optical potential) depending on your purpose. The radial wave function inside the nuclear potential, 𝑅𝑖𝑛 = 𝑈/𝑟 is numerically calculated through solving Eq. (A.4). Outside the potential, we get an asymptotic radial wave function which goes as combines the regular (𝐹𝑙) and irregular (𝐺𝑙) Coulomb wave functions. Please refer to the references for more details [107, 108]. The asymptotic radial wave function can be expressed as: 𝑅𝑎𝑠 (𝑟)𝑟→∞ = 𝑒𝑖𝛿𝑙 (cos(𝛿𝑙)𝐹𝑙 (𝜂, 𝑘𝑟) + sin(𝛿𝑙)𝐺𝑙 (𝜂, 𝑘𝑟)) , ∝ 𝐹𝑙 (𝜂, 𝑘𝑟) + tan(𝛿𝑙)𝐺𝑙 (𝜂, 𝑘𝑟). (A.6) where 𝜂 = 𝑍1𝑍2𝛼√︁ 𝑚 scattering bodies and 𝛼 = 1/137 is the fine-structure constant. 2𝐸 is the Sommerfeld factor, where 𝑍1 and 𝑍2 are the charges of the two The phase shift tan(𝛿𝑙) may directly obtained from the logarithmic derivative of the inner solution 𝑅𝑖𝑛 at the matching point 𝑟 = 𝑟0, i.e., when 𝑉 (𝑟) → 0. From, 𝛽 = 𝑑 𝑑𝑟 (ln(𝑅𝑖𝑛 (𝑟)))|𝑟=𝑟0 , and comparing with asymptotic solution, we get; tan(𝛿𝑙) = 𝐹𝑙 − 𝛽𝐹′ 𝑙 𝐺𝑙 − 𝛽𝑙𝐺′ 𝑙 . (A.7) (A.8) where 𝐹′ 𝑙 and 𝐺′ 𝑙 are derivetive of 𝐹𝑙 and 𝐺𝑙 with respect to 𝑟 at 𝑟 = 𝑟0. A.2 Resolution to correlation function Experimental correlation functions are often affected by the energy resolution of the detectors. In order for theoretical correlations to reproduce the features observed in the measured correlation, we need to incorporate the resolution in the theoretical correlation. This is achieved by folding the calculated correlation with a Gaussian function of an appropriate width (𝜎𝑞). 115 The folded correlation can be expressed as in the equation below: ∫ 𝐶 (𝑞′) = 𝐺 (q, q′)𝐶 (q)𝑑3q, (A.9) where, 𝐺 (q, q′) is assumed to be the Gaussian function of momenta vectors q and q′, 𝐺 (q, q′) ∝ − (q−q′ )2 2𝜎2𝑞 and 𝐶 (𝑞) represents the original calculation of the correlation function. 𝑒 𝐶 (𝑞′) ∝ ∫ − (q−q′ )2 2𝜎2𝑞 𝑑3𝑞, 𝐶 (𝑞)𝑒 (A.10) ∫ ∝ 2𝜋 𝐶 (𝑞)𝑒 − 𝑞2+𝑞′2 −2𝑞𝑞′ 𝑐𝑜𝑠 𝜃 2𝜎2𝑞 𝑞2𝑑𝑞𝑑 (𝑐𝑜𝑠𝜃). By integrating over 𝑐𝑜𝑠𝜃 and performing some manipulations, we obtain a simplified expression; 𝐶 (𝑞′) ∝ ∫ 2𝜋𝜎2 𝑞 𝑞′ (cid:34) 𝑑𝑞𝑞𝐶 (𝑞) 𝑒 − (𝑞+𝑞′ )2 2𝜎2𝑞 − 𝑒 − (𝑞−𝑞′ )2 2𝜎2𝑞 (cid:35) . Thus, the correlation with the correction can be written as; 𝐶 (𝑞′) ∝ ∫ 2𝜋𝜎2 𝑞 𝑞′ (cid:34) 𝑑𝑞𝑞𝐶 (𝑞) 𝑒 − (𝑞+𝑞′ )2 2𝜎2𝑞 − 𝑒 − (𝑞−𝑞′ )2 2𝜎2𝑞 (cid:35) . (A.11) (A.12) The integral in Equation (A.12) is taken over the entire range of relative momentum to include all resonance peaks. 116 APPENDIX B THERMODYNAMIC MODEL FOR PARTICLE EMISSION The thermodynamic model, such as the fireball thermal model, allows us to study various properties of particle emissions in nuclear reactions. This model is employed to compute the density distribu- tion of particles emitted in heavy-ion collisions. Subsequently, the density can be utilized to study the correlation between particles in the final state emissions. In Chapter 2, we used the density distribution of alpha particles produced from the decay of the 9Be nucleus to calculate correlations between them. To determine the particle density distribution, one first needs to determine the thermodynamic parameters that characterize our system, such as temperature, proton and neutron chemical potentials, and then use these parameters in Eq. (B.1). Let us consider a non-interacting gas. In momentum space, the density distribution of particles of type 𝑖 is defined as [46]; 𝑑𝑁𝑖 𝑑3𝑃𝑖 = ≈ (2𝑆𝑖 + 1)𝑉 (2𝜋ℏ)3 (cid:16) 𝛽 exp (2𝑆𝑖 + 1)𝑉 (2𝜋ℏ)3 (cid:18) exp −𝛽 1 (cid:104) 𝑃2 2𝑚 − 𝐵.𝐸 (cid:20) 𝑃2 2𝑚𝑖 (cid:105) (cid:17) − 𝛽𝜇𝑖 (cid:21) ± 1 (cid:19) − 𝐵.𝐸 + 𝛽𝜇𝑖 . (B.1) The total average number of particles of species 𝑖 is obtained by integrating the equation above (Eq. B.1) over momentum space: < 𝑛𝑖 >= 𝑉 (2𝑠𝑖 + 1) ℏ3 (cid:18) 𝑚𝑖𝑇 2𝜋 (cid:19) 3/2 ( 𝜇𝑖 −𝐵.𝐸 ) 𝑇 , 𝑒 (B.2) where 𝜇𝑖 = 𝜇𝑛𝑁𝑖 + 𝜇𝑧 𝑍𝑖. Here, 𝜇𝑛 and 𝜇𝑧 are the neutron and proton chemical potentials, respectively. 𝑇 is the temperature, 𝑉 is the volume of the system, 𝐵.𝐸 is the binding energy, and 𝑠𝑖 is the spin statistics of particle 𝑖. In Equation B.2, we fix temperature and solve for are 𝜇𝑛 and 𝜇𝑧. To find them, we need to solve the following system of equations: 𝑍𝑡𝑜𝑡 = (cid:205) 𝑍𝑖 < 𝑛𝑖 >, 𝑁𝑡𝑜𝑡 = (cid:205) 𝑁𝑖 < 𝑛𝑖 > .    117 Here, 𝑁𝑖 and 𝑍𝑖 are the number of neutrons and protons of species 𝑖, and 𝑁tot and 𝑍tot are the total number of neutrons and protons, respectively, in the original system. The proton and neutron chemical potentials were calculated for a system of volume 𝑉 = 𝐴/𝜌 and 𝜌 = 𝑐𝜌0, for a fixed temperature (T). To verify our calculations, we calculated the ratio of alpha to proton and compared it with the results published in reference [109]. The comparison is shown in Fig. (B.1), which demonstrates good agreement. Figure B.1 Ratio of Alpha to proton numbers versus system volume: The dashed line represents our calculation (we used parameters 𝑇 = 50 MeV, 𝜌0 = 0.15, 𝑍tot = 30, and 𝑁tot = 30 from Refs. [109]). The blue points display the calculation from Ref. [109]. 118 APPENDIX C QUADRUPOLE COMPONENTS OF SOURCE FUNCTION FROM BUU SIMULATION Here, we represent the quadrupole components of p-p source functions when a cut is applied in different directions of momenta. Fig. C.1 illustrates the p-p quadrupole component of the source for the BUU simulation of 131Xe+197Au reaction at E/A=80 MeV. The first row represents the soft nuclear equation of state (EOS), while the second row represents the stiff EOS. In both rows, the solid line represents the momentum-independent EOS, and the dots represent the momentum- dependent EOS. The columns represent cuts in momentum, from left to right: the cut in total momentum, 𝑃, 200-400 MeV/c; the cut in transverse momentum, 𝑃𝑇 = √︃ 𝑃2 𝑥 + 𝑃2 𝑦, 0-400 MeV/c; the cut in the x-component of the momentum, 𝑃𝑥 = 𝑃𝑇 cos(𝜙), 0-400 MeV/c; and the cut in the y-component of the momentum, 𝑃𝑦 = 𝑃𝑇 sin(𝜙), 0-400 MeV/c, where 𝜙 is azumithal angle. Notably, there is a significant difference between the soft and momentum-dependent soft EOSs, especially in the case of the PT cut. The similar description is applied to Fig. C.2 for 36Ar+45Sc reaction at a beam energy of 80 MeV. Figure C.1 Two-proton source function as a function of relative position obtained using the BUU simulation code for 131Xe+197Au at a beam energy of 80 MeV/nucleon and an impact parameter of b=3.1 fm. The solid and dashed-dot curves correspond to the stiff and soft momentum-independent EoS, respectively. The dashed and solid curves superposed with stars represent the soft and stiff momentum-dependent EoS, respectively. 119 r (fm)0.00000.00010.00020.0003S20 (r)P: 200400 MeV/cr (fm)PT: 0400 MeV/cr (fm)Px: 0400 MeV/cr (fm)Py: 0400 MeV/cSoftSoft_mom0510152025r (fm)0.00000.00010.00020.0003S20 (r)0510152025r (fm)0510152025r (fm)0510152025r (fm)StiffStiff_mom Figure C.2 Two-proton 𝑆20 component of the source function as a function of relative position obtained using the BUU simulation code for 36Ar+45Sc at a beam energy of 80 MeV/nucleon and an impact parameter of b=3.1 fm. The solid and dashed-dot curves correspond to the stiff and soft momentum-independent EoS, respectively. The dashed and solid curves superposed with stars represent the soft and stiff momentum-dependent EoS, respectively. 120 r (fm)0.00000.00020.00040.00060.0008S20 (r)P: 200400 MeV/cr (fm)PT: 0400 MeV/cr (fm)Px: 0400 MeV/cr (fm)Py: 0400 MeV/cSoftSoft_mom02468101214r (fm)0.00000.00020.00040.00060.0008S20 (r)02468101214r (fm)02468101214r (fm)02468101214r (fm)StiffStiff_mom