CONSTRAINING NUCLEON EFFECTIVE MASS SPLITTING USING NEUTRON AND PROTON OBSERVABLES FROM HEAVY-ION COLLISIONS By Chi En Teh A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics—Doctor of Philosophy 2023 ABSTRACT This dissertation focuses on the study of the nuclear equation of state, with an emphasis on the extraction of effective mass splitting between neutron and proton. A series of experiments was proposed and approved to use the neutron-to-proton (n/p) ratios from heavy-ion collisions as a probe to constrain the nuclear equation of state. The latest experiments were carried out at the National Superconducting Cyclotron Laboratory (NSCL) using three detection systems: the Microball, for extracting impact parameter information, the HiRA10 silicon array, for detecting light charged particles from hydrogen to helium isotopes, and the vLANA neutron wall, for measuring neutrons. 40 48 Each of the two calcium isotope beams, Ca and Ca, were delivered at two different energies, 56 MeV/u and 140 MeV/u, and targeted onto stationary targets. These targets 58 68 112 124 were enriched with either nickel isotopes, Ni and Ni, or tin isotopes, Sn and Sn, to produce a range of isospin asymmetries, (N − Z)/A, with values ranging from 0.02 to 0.14. A total of 16 beam-target configurations were explored in the experiment, but this dissertation 40 primarily analyzes beam-target pairs with the extreme asymmetry values, specifically, Ca + 58 48 Ni, Ca + 64Ni, 40 Ca + 112Sn, and 48 Ca + 124Sn, at both beam energies. Key findings include the successful reconstruction of light charged particle spectra with HiRA10 and the revelation of the isoscaling property across the entire measured range. The isoscaling of charged particles enables the construction of pseudo-neutron and corresponding coalescence-invariant (CI) neutron spectra. Due to large systematic errors, the CI n/p ratios from one single reaction may not be reliable. These limitations were greatly mitigated by constructing the double ratio of CI n/p ratios from the most and the least neutron-rich systems. The double ratio was then analyzed through a Bayesian framework informed by the Improved Quantum Molecular Dynamics (ImQMD) model, showing sensitivity to the nucleon effective mass splitting. The second part of this work dedicated significant effort towards the calibration and data analysis of vLANA, a neutron wall array equipped with a veto wall to discriminate against charged particles. The in-depth calibration resulted in improvements in accuracy, precision, and effective handling of large variety of data obtained over a period of two months when the experiments were performed. The intrinsic efficiency of neutron detection was evaluated through a sophisticated simulation using a newly developed neuSIM4 code. The study culminated with the presentation of neutron spectra for the Ca + Ni systems at incident energy of 140 MeV/u, confirming the isoscaling property of neutrons. Furthermore, the substitution of pseudo-neutrons with actual reconstructed neutrons in the CI n/p double ratio suggested ∆m∗np /δ = −0.10 ± 0.09, aligning with the findings from the light charged particle analysis and previous works. In summary, the methodology developed in this dissertation paves a way for future investigations of effective mass splitting, contributing to our understanding of the nuclear equation of state. It also underscores the need for ongoing collaboration with theorists and calls for further exploration using other theoretical models. ACKNOWLEDGEMENTS A deep sense of gratitude engulfs me as I pen down my acknowledgements to the numerous individuals who significantly contributed to the completion of my graduate studies. Their shared wisdom, encouragement, and companionship have immensely enriched this incredible journey. At the forefront of this distinguished group is my advisor, Prof. Betty Tsang. Her steadfast dedication to teaching, infectious passion for scientific discovery, and continuous drive for excellence were an ever-glowing lighthouse throughout my academic pursuit. My sincere thanks go to my committee members, Prof. William Lynch, Prof. Pawel Danielewicz, Prof. Laura Chomiuk, and Prof. Edward Brown. Their valuable perspectives, enduring support, and wise counsel were instrumental in honing my research and broadening my understanding of our field. Invaluable credit is due to the faculty and post-doctoral members of the HiRA group, particularly Prof. Zbigniew Chajecki, Prof. Kyle Brown, and Rensheng Wang, Som Paneru, Daniele Dell’Aquila, Chenyang Niu, Genie Jhang, Taras Lokotko, and Curtis Hunt. Their deep expertise, thoughtful instruction, and professional advice significantly shaped my research journey. I am profoundly grateful to my fellow students, some of whom have graduated and advanced to academia or industry, including Chun Yuen Tsang (Tommy), Kuan Zhu, Jeonghyeok Park, Sean Sweany, Om Khanal, Joseph Wieske, Hongyi Wu, Adam Anthony, Jonathan Barney, Justin Estee, Juan Manfredi, Ryan Yeung, Poulomi Dey, Jessica Ranshaw, Zac Benzerara, and Mira Ghazali. Their open-hearted exchange of knowledge and skills remarkably enriched my journey, and their dedication to science continually inspired me. I am particularly indebted to Prof. Jenny Lee, whose support and inspiration during my undergraduate studies in Hong Kong paved the way for my pursuit of a PhD in experimental nuclear physics at Michigan State University. My thanks also go to the diligent board members of the Artificial Intelligence (AI) iv club at Michigan State University. Their commitment to AI and education was a constant wellspring of inspiration. Likewise, I extend my heartfelt gratitude to the Malaysian Student Organization (MSO) at Michigan State University for their supportive companionship during my academic expedition. I would like to express my special appreciation to my friend, Alex Chandra. His cama- raderie throughout my academic adventure was a comforting source of strength. My family, whose love and support spanned oceans and time zones, deserve my utmost gratitude. Special mention goes to Ain Hanani and Cadence for their undying affection. As I recount my gratitude, I would like to appreciate the enlightening conversations on physics, life, and culinary explorations I have had with everyone mentioned. These discussions, frequently laced with humor, have been instructive, insightful, and will remain a cherished part of my academic journey. The financial support from the U.S. National Science Foundation, under Grant No. PHY-1565546, and the U.S. Department of Energy (Office of Science), under Grant No. DE-NA0003908, greatly contributed to making this work possible. Fanurs C.E. Teh East Lansing, MI July 7th, 2023 TABLE OF CONTENTS LIST OF ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Equation of state and nuclear matter . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Liquid-drop model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Nuclear equation of state . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.3 Symmetry energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.1.4 Effective mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2 Heavy-ion collisions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.2.1 Spectral yield ratios . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.3 Organization of dissertation . . . . . . . . . . . . . . . . . . . . . . . . . 16 CHAPTER 2 EXPERIMENT E15190-E14030 . . . . . . . . . . . . . . . . . . . 18 2.1 Impact parameter detection system – Microball . . . . . . . . . . . . . . 21 2.1.1 Adaptation for the experiment . . . . . . . . . . . . . . . . . . . . . . 23 2.1.2 Impact parameter determination . . . . . . . . . . . . . . . . . . . . 25 2.2 Charged-particle detection system – HiRA10 . . . . . . . . . . . . . . . . 30 2.2.1 Configuration for the experiment . . . . . . . . . . . . . . . . . . . . 31 2.2.2 Particle identification of charged particles . . . . . . . . . . . . . . . 32 2.3 Neutron detection system – vLANA . . . . . . . . . . . . . . . . . . . . . 35 2.3.1 Large area neutron array . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.3.2 Forward array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.3.3 Veto wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 CHAPTER 3 RESULTS FROM CHARGED PARTICLES . . . . . . . . . . . . 46 3.1 Charged-particle spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.1.1 Energy spectra in lab frame . . . . . . . . . . . . . . . . . . . . . . . 46 3.1.2 Transverse momentum spectra . . . . . . . . . . . . . . . . . . . . . . 48 3.2 Isoscaling from charged particles . . . . . . . . . . . . . . . . . . . . . . . 55 3.2.1 Ca + Ni systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.2.2 Ca + Sn systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.3 Neutron-proton double ratio . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.3.1 Coalescence-invariant spectra . . . . . . . . . . . . . . . . . . . . . . 66 3.4 Bayesian analysis of double ratio spectrum . . . . . . . . . . . . . . . . . 70 3.4.1 Gaussian Process emulation . . . . . . . . . . . . . . . . . . . . . . . 71 3.4.2 Bayesian inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 CHAPTER 4 CALIBRATIONS FOR THE NEUTRON DETECTION SYSTEM 83 4.1 Hit position calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.1.1 Coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.1.2 Calibrating hit positions using Veto wall bar shadows . . . . . . . . . 87 4.2 Time-of-flight calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.2.1 Using Forward array as a reference for event start time . . . . . . . . 92 vi 4.2.2 Calibrating ToF using prompt gamma events . . . . . . . . . . . . . . 94 4.3 Pre-processing of analog-to-digital converter values . . . . . . . . . . . . 97 4.3.1 Non-linearity correction . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.3.2 Saturation correction . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.4 Light output calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.4.1 Removal of position dependence . . . . . . . . . . . . . . . . . . . . . 106 4.4.2 Calibration of light output . . . . . . . . . . . . . . . . . . . . . . . . 108 4.5 Pulse shape discrimination . . . . . . . . . . . . . . . . . . . . . . . . . . 109 CHAPTER 5 RECONSTRUCTION OF NEUTRON SPECTRA . . . . . . . . . 116 5.1 Background analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.1.1 Scattered neutron background . . . . . . . . . . . . . . . . . . . . . . 116 5.1.2 Random background at radio frequency . . . . . . . . . . . . . . . . . 121 5.2 Neutron detection efficiencies . . . . . . . . . . . . . . . . . . . . . . . . 126 5.2.1 Modifications for calculating the neutron energy spectrum . . . . . . 128 5.2.2 Geometric efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 5.2.3 Intrinsic efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 5.3 Neutron spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 5.3.1 Transverse momentum spectra . . . . . . . . . . . . . . . . . . . . . . 146 5.3.2 Spectral yield ratios . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 CHAPTER 6 SUMMARY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 APPENDIX A SPECTRAL YIELDS FOR LIGHT CHARGED PARTICLES . . 184 APPENDIX B TRANSVERSE MOMENTUM VERSUS RAPIDITY . . . . . . . 193 vii LIST OF ABBREVIATIONS ADC analog-to-digital converter ANM asymmetric nuclear matter CFD constant fraction discriminator CI coalescence invariant DAQ data acquisition DSSD double-sided silicon strip detector EoS / EOS equation of state FA Forward array FRIB Facility for Rare Isotope Beams FWHM full width at half maximum GM geometric mean GP Gaussian process HIC heavy-ion collision HiRA High Resolution Array ImQMD Improved quantum molecular dynamics LANA Large area neutron array LHS Latin hypercube sampling MCMC Markov chain Monte Carlo NM neutron matter NSCL National Superconducting Cyclotron Laboratory NUTS No-U-Turn Sampler NW Neutron wall NWA Neutron wall A NWB Neutron wall B viii PID particle identification PMT photomultiplier tube PPSD position-corrected VPSD PSD pulse shape discrimination RF radio frequency QDC charge-to-digital converter QMD quantum molecular dynamics SEMF semi-empirical mass formula SNM symmetric nuclear matter TDC time-to-digital converter ToF / TOF time-of-flight vLANA Large area neutron array with Veto wall VPSD value-assigned pulse shape discrimination VW Veto wall ix CHAPTER 1 INTRODUCTION Navigating the intricacies of an atom’s nucleus reveals a rich, complex landscape of nuclear physics. This thesis embarks on a journey into these depths, focusing especially on experi- mental efforts to measure neutrons and charged particles resulting from heavy-ion collisions. The comparison of these measurements with theoretical predictions provides valuable insights into nucleon effective mass splitting. Our exploration begins with a comprehensive overview of the fundamental concepts and theories, building towards an understanding of nucleon effective mass and its implications. Subsequently, we turn our attention to heavy-ion collision, a terrestrial tool for probing the nature of nuclear matter. From these collisions, we derive observables that serve as sensitive probes for the nuclear equation of state and effective mass splitting. The chapter concludes with a detailed roadmap of the thesis. 1.1 Equation of state and nuclear matter An equation of state (EOS) is a function, often rooted in theory or empirical observation, which describes the state of matter under a specific set of physical conditions. It achieves this by establishing relationships between thermodynamic variables, such as pressure, volume, temperature, and internal energy. Examples in thermodynamics include the ideal gas law and the van der Waals equation, which describe the behavior of gases and liquids without the need to track the exact trajectories and positions of individual molecules. The concept of EOS is wide-ranging in its application, found in fields from thermodynamics to astrophysics, providing a theoretical foundation necessary to comprehend the properties and behavior of various forms of matter. Within the realm of nuclear physics, the EOS is also applied to describe nuclear matter. Like its counterparts in other fields, a nuclear EOS needs to be able to reproduce the observed properties of matter without having to replicate the kinematics of individual neutrons and protons. 1 1.1.1 Liquid-drop model Introduced by George Gamow in 1930, the liquid-drop model offers a means to understand the structure and behavior of atomic nuclei [1]. It emerged during a period of groundbreaking advancements in nuclear physics, drawing on the conceptual framework of classical physics and incorporating elements from the then-nascent field of quantum mechanics. By comparing an atomic nucleus to a droplet of liquid, the model simplifies the complexity of an atomic nucleus into a few state variables. It presents an intuitive understanding of collective motion in a nucleus, allowing for the use of familiar principles of classical mechanics. This model presupposes a spherical shape for nuclei, an assumption which, while not always exact, provides a first-order approximation facilitating analysis and prediction. Despite its oversimplifications and inherent limitations, the liquid-drop model is a valuable pedagogical tool in nuclear physics till this day. It is often employed as an initial step for understanding more intricate models, such as the nuclear shell model, and serves as a foundational base in the study of nuclear deformation, fission, and heavy-ion collisions (HIC). The state variables in this model include the proton number Z, the neutron number N , the mass number A = N + Z, and the binding energy BE(N, Z). To specify a nucleus, it is only necessary to know two of the three numbers N , Z, and A. Throughout this thesis, we will use variable pairs such as (N, Z) and (A, Z) interchangeably to refer to a nucleus, with the understanding that the third number is determined by the other two. Building on Gamow’s initial proposal, Carl Friedrich von Weizsäcker formulated the semi-empirical mass formula (SEMF) in 1935 [2], sometimes referred to as the Weizsäcker formula. This equation varies slightly in different literature, but generally takes the following form: Z(Z − 1) (N − Z)2 1 BE(N, Z) = αV A − αS A2/3 − αC 1/3 − αA + αP 1/2 , (1.1) A A A where αV , αS , αC , αS are empirical parameters obtained by fitting to existing nuclide mass. Each of these terms represents different physical phenomenon: • +αV A represents the volume term, which considers the nucleus as a homogeneous 2 sphere where each nucleon interacts only with its nearest neighbors. • −αS A2/3 is the surface energy term, accounting for the fact that nucleons on the surface have fewer neighboring nucleons. Z(Z − 1) • −αC represents the Coulomb term, resulting from the repulsive force between A1/3 protons. (N − Z)2 • −αA is the asymmetry term, expressing the energy cost associated with A unequal numbers of neutrons and protons as a consequence of the Pauli exclusion principle. 1 • +αP is the pairing term, which accounts for the spin coupling of nucleons. A1/2 It is crucial not to underestimate the limitations of this model. Liquid-drop model does not consider the shell structure, deformation, and it is unable to predict the binding energy of unstable nuclei accurately. It also performs inadequately for light nuclei, where the droplet metaphor becomes less applicable due to the small number of nucleons and the much larger contributions from the surface effects. coefficient value (MeV) aV 14.64 aS 14.08 aC 0.64 aA 21.07  +11.54 if (N, Z) is even-even  aP −11.54 if (N, Z) is odd-odd  0 if (N, Z) is even-odd or odd-even  Table 1.1 Values of coefficients in the semi-empirical mass formula (SEMF) as reported in [3]. The value of the aP coefficient in the pairing term varies based on the parities of N and Z. The SEMF coefficient values depend on the mass data used and the fitting method. Table 1.1 presents a set of coefficients obtained by fitting to 2166 isotopes in AME2016 3 [4, 5] with A ≥ 50 as reported in [3]. Isotopes with A < 50 are excluded, aligning with the acknowledged limitation of the SEMF for light nuclei. 1.1.2 Nuclear equation of state Nuclear matter is theoretical construct of an infinite system of nucleons with homogeneous neutron and proton densities. Despite being an idealized medium, it offers important insights into the dynamics of heavy-ion collisions (HICs) [6, 7, 8, 9, 10, 11] and the structure of neutron stars [12, 13, 14, 15, 16]. In this model, Coulomb interactions are assumed to be absent and the isospin formalism is employed. Nuclear matter is essentially a quantum many-body system composed of two types of fermions, the neutron and proton. Despite their differing electric charges, these fermions share similar masses (939.6 MeV/c2 for the neutron and 938.3 MeV/c2 for the proton) and exhibit indistinguishable behavior under the strong interaction. This property motivates the introduction of Wigner’s isotopic spin [17], more commonly known as isospin, which encapsulates the neutron-proton exchange symmetry. Under this formalism, the neutron and proton are treated as two states of the same particle, the nucleon, with an isospin of 1/2. The convention in nuclear physics denotes the neutron and proton with isospin projections of +1/2 and −1/2, respectively [18]. In nuclear matter, the energy per nucleon E/A = E = E(ρ, δ) depends on both the baryon density, ρ = ρn + ρp , (1.2) and the isospin asymmetry, ρn − ρp δ= , (1.3) ρ where ρn and ρp denote the neutron and proton number densities, respectively. E is often expressed as a Taylor expansion in δ about the symmetric nuclear matter (SNM) at δ = 0: E(ρ, δ) = ESNM (ρ) + S(ρ)δ 2 + O(δ 4 ) , (1.4) where ESNM (ρ) ≡ E(ρ, 0) and S(ρ) denotes the symmetry energy. Contrary to its name, the 4 symmetry energy characterizes the contribution due to the presence of asymmetry, drawing a parallel to the asymmetry term in the SEMF. It is worth noting the absence of odd powers of δ in Equation 1.4. This absence can be understood in terms of isospin symmetry, as the nuclear interaction does not significantly favor either protons or neutrons. Mathematically, this translates to the fact that the system’s energy would not change under a transformation that swaps protons and neutrons. In reality, protons, unlike neutrons, carry an electric charge and experience additional Coulomb repulsion. This condition breaks the symmetry between protons and neutrons, altering the behavior of nuclear matter. The Coulomb interaction is typically considered as a separate term or correction in nuclear physics calculations [19, 20, 21]. To better understand ESNM , it is common to consider its Taylor expansion in ρ about ρ0 : 1 1 ESNM = E0 + K0 x2 + Q0 x3 + O(x4 ) , (1.5) 2 6 where ρ0 is defined as the saturation density of SNM such that ∂ESNM =0, (1.6) ∂ρ ρ=ρ0 and x is a dimensionless parameter introduced as ρ − ρ0 x≡ . (1.7) 3ρ0 Here, the factor of 3 in the denominator arises from historical convention. Prior studies in the early 1960s [22, 23] adopted internucleonic separation r as the variable for expanding E/A rather than ρ, which is more commonly used today. Hence, for consistency, the second- ∂ 2E order coefficient is expressed as 32 ρ2 2 . The general form of the nth-order coefficient in ∂ρ Equation 1.5 is then given by ∂ n ESNM 3n ρn0 , (1.8) ∂ρn ρ=ρ0 which always has a unit of energy. The terms within the expansion have distinct physical interpretations. The leading term 5 E0 represents the energy per nucleon of SNM at ρ0 . Experimentally, E0 ≈ −16 MeV and ρ0 ≈ 0.16 fm−3 have long been established from the properties of stable nuclei. While the linear term is absent, it is worth noting its relation to the pressure of SNM, ∂ESNM ∂ESNM PSNM = − = ρ2 , (1.9) ∂V ∂ρ where V denotes volume and ESNM = AESNM is the total energy. At equilibrium (ρ = ρ0 ), the pressure of SNM is expected to vanish. The quadratic term in the expansion corresponds to the incompressibility coefficient, ∂ 2 ESNM K0 = 9ρ20 . (1.10) ∂ρ2 ρ=ρ0 Giant monopole resonance (GMR) measurements have led to estimates of K0 = 230 ± 30 MeV [23, 24, 25, 26, 27]. Lastly, the cubic term, ∂ 3 ESNM Q0 = 27ρ30 , (1.11) ∂ρ3 ρ=ρ0 is associated with skewness. Current experimental data provide insufficient constraint for its value [28]. 1.1.3 Symmetry energy The symmetry energy is a fundamental quantity in nuclear physics that characterizes the energy cost associated with creating neutron-proton asymmetry in nuclear matter. Analogously to the expansion of ESNM in Equation 1.5, the symmetry energy, 1 ∂ 2 E(ρ, δ) S(ρ) ≡ , (1.12) 2 ∂δ 2 δ=0 can be similarly expanded in ρ about ρ0 (x = 0): 1 1 S(ρ) = S0 + Lx + Ksym x2 + Qsym x3 + O(x4 ) , (1.13) 2 6 6 where the coefficients are given by S0 = S(ρ0 ) , (1.14) ∂S(ρ) L = 3ρ0 , (1.15) ∂ρ ρ=ρ0 ∂ 2 S(ρ) Ksym = 9ρ20 , (1.16) ∂ρ2 ρ=ρ0 ∂ 3 S(ρ) Qsym = 27ρ30 . (1.17) ∂ρ3 ρ=ρ0 Here, S0 denotes the magnitude of the symmetry energy at the saturation density ρ0 , sometimes also represented as J [28, 29]. Since the first-order term of ESNM vanishes, the linear term L is proportional to the pressure of nuclear matter at ρ0 . L is also known as the slope parameter. Ksym and Qsym are known as the curvature and skewness parameters, respectively. The higher-order terms beyond Qsym have been neglected in this expansion. Despite its fundamental role in diverse nuclear and astrophysical phenomena, symmetry energy remains one of the least constrained aspects of the nuclear equation of state. A plethora of laboratory experiments have strived to constrain S(ρ) at sub-saturation densities 2 (ρ ≲ ρ0 ), with ρ ≈ ρ0 being the most well-constrained density region based on mass 3 measurements. These include investigations of the excitation energies of Isobaric Analog States (IAS) [30], Pygmy Dipole Resonances (PDR) [31, 32], and electric dipole polarizability studies [33, 34]. These techniques, along with studies of isospin diffusion in heavy-ion collisions (HICs) [35, 36, 37], have provided valuable insights into the behavior of S(ρ) at sub-saturation densities. In the supra-saturation density regime (ρ ≳ ρ0 ), the symmetry energy becomes even more elusive. Despite the challenges, understanding the symmetry energy at supra-saturation densities is of paramount importance, particularly in the context of neutron stars, which exist in a high-density environment beyond normal nuclear densities. Neutron stars, partially supported against further gravitational collapse by the neutron degeneracy pressure due to the Pauli exclusion principle, require the repulsive nuclear forces from the symmetry energy 7 to exist beyond 0.7M⊙ [38, 39], a maximum mass limit first established by Oppenheimer and Volkoff in 1939 using a simple degenerate cold Fermi gas model [40]. Modern observations indicate that neutron stars typically have masses around 1.4M⊙ , and current observed neutron star masses range from ∼ 1.2M⊙ to ∼ 2.2M⊙ [41, 42, 43, 44], far exceeding the maximum mass limit supported by neutron degeneracy pressure alone. Therefore, S(ρ) significantly influences the structural properties of neutron stars [45], including the radius [44, 46], moment of inertia [47, 48], and cooling mechanisms [49, 50]. The landmark observation of a neutron star merger event, GW170817 [51], underscored the critical importance of the symmetry energy in accurately characterizing these extraordinary astrophysical laboratories [52, 53, 54]. Figure 1.1 The x-axis represents the baryon density relative to ρ0 , and the y-axis represents the symmetry energy S(ρ) in MeV. Various curves represent different Skyrme parametrization of symmetry energy used in [55]. The cyan shaded region is constrained from heavy-ion collision experiments. This figure is taken from [56]. Figure 1.1 offers an illustrative snapshot of the density dependence of symmetry energy based on a select set of parametrizations from [55]. The cyan shaded region represents constraints from HIC experiments. A band of constraints at sub-saturation densities is 8 emerging, with the most stringent constraint at ρ ∼ 0.7ρ0 indicating a symmetry energy of ∼ 27MeV. For densities far above ρ0 or below 0.7ρ0 , the symmetry energy remains largely unconstrained. Nonetheless, our understanding of the symmetry energy continues to be refined as new experimental data and theoretical insights become available. Despite being a cornerstone of nuclear physics and astrophysics, the symmetry energy still harbors numerous mysteries, particularly regarding its high-density behavior. The properties of neutron stars do give us insight into the total pressure of these high-energy systems, but a comprehensive understanding of nuclear physics is necessary to decipher what fraction of this pressure arises from the symmetry energy, and what fraction comes from the symmetry matter EoS. Continued theoretical and experimental investigations, along with developments in neutron star astrophysics and gravitational wave astronomy, are essential for unraveling these mysteries and thus deepening some of our understanding of the dense matter equation of state. 1.1.4 Effective mass Although nuclear EOS is useful in outlining the bulk or macroscopic properties of nuclear matter, it does not provide detailed description regarding the individual behavior of nucleons. Ab initio calculations on nuclear many-body systems present a potential alternative. However, such calculations are not only computationally intensive and often limited to smaller systems, but they are also influenced by our limited understanding of the high-density behavior of three and four-body interactions [57]. Even approximations such as the Chiral effective field theory, which applies a high momentum cutoff, are not exempt from these challenges. Thus, an approach that balances computational feasibility and accurate reflection of individual nucleonic behavior within nuclear matter is sought. One common approach that addresses this balance is the mean-field theory. This ap- proximation theory allows us to view nucleons as moving independently within an averaged potential field. Various theories, such as the Hartree-Fock (HF) theory [58, 59], its extension the Skyrme-Hartree-Fock (SHF) [60], and the more complex Dirac-Brueckner Hartree-Fock 9 (DBHF) theory [61], have been developed within this framework. A distinctive feature of these phenomenological theories is their flexibility. Their interactions can be adjusted to change their high-density behavior without significantly affecting their low-density predictions, making them versatile tools for studying nuclear matter. Inside a mean-field potential, the behavior of nucleons can be understood as though they possess a mass different from their free mass. This understanding led to the introduction of the concept of an effective mass, first proposed by K. A. Brueckner in 1955 [62]. The effective mass was introduced to explain the momentum dependence of the single-nucleon potential, Uτ (ρ, δ, k), in nuclear matter with density ρ and isospin asymmetry δ. Here, Uτ depends explicitly on the momentum k (wavenumber), and the subscript τ denotes the nucleon type, which can be either neutron (n) or proton (p). To motivate the definition of effective mass, we consider the Hamiltonian of a nucleon, k2 H= + Uτ (ρ, δ, k) , (1.18) 2mτ in which the natural unit ℏ = 1 is used, and mτ denotes the nucleon free mass. To write the Hamiltonian equation, we differentiate the Hamiltonian with respect to momentum: ∂H k dUτ (ρ, δ, k) = + (1.19) ∂k mτ dk   k mτ dUτ = 1+ . (1.20) mτ k dk This motivates the definition of effective mass as [63, 64, 65] !−1 mτ dUτ m∗τ ≡ mτ 1 + (τ ) , (1.21) kF dk (τ ) k=kF (τ ) which is conventionally evaluated at the Fermi momentum kF . From Equation 1.21, it can be seen that the effective mass m∗τ is a function of both density ρ and isospin asymmetry δ. When comparing the values of effective masses from empirical data and theoretical calculations, it is common to quote the values at the saturation density ρ0 ≈ 0.16 fm−3 . Throughout this thesis, it is understood that when m∗τ is quoted for a single 10 value, it is the value at ρ0 ; otherwise, m∗τ is assumed to be a function of ρ. The dependence on δ can be removed by introducing the isoscalar and isovector effective masses, m∗s and m∗v , respectively. In the limit of SNM (δ = 0), neutron-proton exchange symmetry implies m∗n = m∗p ≡ m∗s . For general values of δ, the conversion between (m∗s , m∗v ) and (m∗n , m∗p ) is given by [66]   1 2ρτ 1 2ρτ 1 = + 1 − , (1.22) m∗τ ρ m∗s ρ m∗v where ρτ denotes the baryon density of type τ = n, p. By definition, both m∗s and m∗v are independent of δ, which sometimes make them more convenient to work with when comparing nuclear systems with different isospin asymmetries. In recent times, a significant amount of scientific attention has been directed towards understanding the nuances of effective masses in neutron-rich matter, particularly focusing on the neutron-proton effective mass splitting: m∗n − m∗p ∆m∗np ≡ , (1.23) mN where mN is the free nucleon mass. The neutron-proton effective mass splitting bears significant relevance to a diverse range of phenomena in nuclear physics and astrophysics. Gerry Brown and his colleagues, half a century ago, noted with surprise, “It seems strange to us that people making calculations in nuclear matter do not worry about this point, however, since it remains in direct conflict with our conclusions, unless understood as a specific effect of the finiteness” [67]. Regrettably, despite its profound significance, the neutron-proton effective mass splitting remains elusive, particularly regarding its sign and magnitude, which have been long-standing subjects of contention and debate in the field. This enduring uncertainty has led to a variety of responses and interpretations, underscoring the complexity of the issue [68]. Certain microscopic calculations, such as the Landau-Fermi liquid theory [69] and the non-relativistic Brueckner-Hartree-Fock (BHF) method [70], have indicated that ∆m∗np > 0. On the contrary, calculations based on relativistic mean-field (RMF) theory and DBHF theory [71, 72] suggest ∆m∗np < 0. 11 An in-depth study of nucleon effective masses in neutron-rich matter can shed light on many unresolved aspects. For instance, the equilibrium neutron-proton ratio in primordial nucleosynthesis is influenced by ∆m∗np [73]. Furthermore, it is vital for comprehending the properties of neutron stars [74, 75, 76], properties of mirror nuclei [77], the determination of the neutron and proton drip-lines [78], and various isospin-sensitive observables in heavy-ion collisions (HICs) [79, 80, 81, 82, 83]. Given that ∆m∗np depends on the isospin asymmetry δ, it may sometimes be more convenient to work with [84] mN mN fI ≡ ∗ − ∗ , (1.24) ms mv which is independent of δ. The conversion between ∆m∗np and fI can be deduced from Equation 1.22 as "  ∗  #−1 ∗ 2 1 1 m − m ∆m∗np = −2δfI · 2 − δ2 s v (1.25) mN (m∗s )2 m∗s m∗v  ∗ 2 ms ≈ −2δfI · . (1.26) mN Here, an approximation has been made by neglecting the second term in the square brackets, which is justified for small values of δ < 1. This formulation allows us to quote the numerical value of ∆m∗np in units of δ, facilitating a more direct comparison with data extracted from nuclear systems with different isospin asymmetries. As we further navigate through this thesis, we will highlight experimental efforts to measure neutron-to-proton ratios (n/p) derived from heavy-ion collisions. These ratios provide crucial insights enabling us to constrain neutron-proton effective mass splitting by employing a Bayesian framework informed by the Improved Quantum Molecular Dynamics (ImQMD) model [85]. 1.2 Heavy-ion collisions Heavy-ion collisions (HICs) are instrumental in nuclear physics for exploring the properties of nuclear matter under extreme conditions of high density and temperature. These collisions 12 essentially create a miniature cosmic “laboratory”, enabling the investigation of nuclear matter properties beyond the saturation density at ρ0 ≈ 0.16 fm−3 , which is typical in the core of stable nuclei. Figure 1.2 Schematic representation of a heavy-ion collision between two nuclei, adopted from [86]. The figure illustrates the progression of the collision over time: from the initial trajectories of the nuclei (left) to the formation of a high-density region (middle) and ultimately the development of transverse pressure (right). The bottom and back panels respectively depict contours of constant density in the reaction plane (xz-plane) and the transverse plane (xy-plane). Studying the complex interplay between compression and expansion dynamics during a collision can provide crucial insights into the EOS. As illustrated in Figure 1.2, a significant aspect of these dynamics is the role played by spectators and participants in shaping the collision. Spectators are nucleons that largely bypass the primary collision zone without substantial interaction, continuing along their initial trajectory. They are represented by the magenta arrows in the three-dimensional schematic. Conversely, participants are the nucleons that actively engage in the collision, contributing to the formation of a high-density region near the collision vicinity. The participation of these nucleons establishes significant 13 pressure gradients, particularly in transverse directions, as depicted in the back panel and by the red arrows in the three-dimensional schematic. The dynamic interplay between participants and spectators during a collision results in the emergence of a strongly interacting system, which is characterized by the redistribution of nucleons and momentum transfer, eventually causes an outward flow of particles from the high-density region. Particles emitted from this region would encode information about the nuclear EOS and other nuclear matter properties at high densities, giving us valuable observables in HIC experiments. 1.2.1 Spectral yield ratios HICs serve as a potent tool to investigate the density dependence of the symmetry energy and the momentum dependence of the symmetry potential. One common HIC observables that is sensitive to these properties is the spectral yield ratios of particles emitted during the collision. The spectral yield of a particle, denoted as Y , is defined as the number of particles emitted per unit phase space volume, averaged over all collision events. This yield is typically represented as a function (spectrum) of energy and transverse momentum. In nuclear physics, there are several spectral yield ratios commonly used to probe the nuclear EOS, including the neutron-to-proton ratio n/p [87, 88, 89, 90], the pion ratio π − /π + [91, 92, 93], and the kaon ratio K 0 /K + [94]. To differentiate with the double ratio that will be introduced later in this section, we refer to these ratios as single ratios. In this thesis, we primarily focus on analyzing the the n/p single ratio, Yn SRn/p ≡ . (1.27) Yp This ratio serves as a direct measure of the relative production rates of neutrons and protons within various segments of the phase space. Nucleons with different isospin states, i.e., neutrons and protons, experience symmetry potentials of opposite signs. Consequently, a variation in the n/p ratio would offer valuable insights into the behavior of nuclear matter under different isospin asymmetries, effectively probing the symmetry energy. 14 The changes in the n/p single ratio can sometimes be qualitatively understood by vi- sualizing the nucleons in the collision as classical entities. For instance, when ∆m∗np < 0, high-momentum neutrons emitted from the compressed participant region experience a more repulsive potential, which results in a higher acceleration process for neutrons than for pro- tons with the same momentum which tends to enhance the n/p ratio at high momentum. Conversely, when ∆m∗np > 0, the acceleration of protons is enhanced, leading to a lower ratio [83]. However, interpreting the n/p single ratio is not always straightforward. Hence, model- ing the intricate dynamics of HICs and interpreting these observables necessitates the use of transport models. These models provide a semi-classical description of the dynamics of nuclear collisions by tracking nucleons and other particles. One such model is the Quantum Molecular Dynamics (QMD) model [95, 96]. In this work, we employ an improved version of the QMD model, the Improved-QMD (ImQMD), developed initially at the China Insti- tute of Atomic Energy (CIAE) by Prof. Zhuxia Li’s research group and later refined by Prof. Yingxun Zhang in the early 2000s [85]. The ImQMD model uses mean-field potentials derived from the Skyrme energy density functional and incorporates explicit momentum dependence. In addition to the single ratio, another valuable observable is the double ratio. This involves comparing the single ratios from two different reaction systems, labeled as system 1 and system 2. Conventionally, for nuclear-rich systems (δ > 0), the system with larger δ is denoted as system 2, and the other as system 1. The double ratio is defined as the single ratio of system 2 divided by the single ratio of system 1, (2) SRn/p DRn/p ≡ (1) . (1.28) SRn/p This quantity is particularly useful as it helps to minimize the impact of systematic errors and model uncertainties, thereby providing a more reliable way to interpret experimental data and constrain the nuclear EOS [89, 90]. Collectively, spectral yield ratios play a critical role in exploring HICs, offering key insights into symmetry energy and effective masses. When coupled with modern transport models 15 like ImQMD, they provide robust tools for understanding the complex dynamics inherent in HICs. 1.3 Organization of dissertation The dissertation is organized into six chapters, including this introductory chapter. Additionally, an extensive bibliography and two appendices are included at the end. Chapter 2 details the experiments conducted under the designations E15190-E14030. This chapter explores the three detection systems utilized in the experiments: the impact parameter detection system (Microball), the charged-particle detection system (HiRA10), and the neutron detection system (vLANA). Each subsection elaborates on the detector configurations as well as the adaptation strategies to achieve the physics goals of this experiment. Chapter 3 discusses the results from charged particle detection. It begins with the reconstruction of energy spectra and transverse momentum spectra. An isoscaling analysis is then performed on the transverse momentum spectra for Ca + Ni and Ca + Sn systems. The coalescence-invariant (CI) spectrum is introduced and used to probe the effective mass splitting using a Bayesian framework informed by the Improved Quantum Molecular Dynamics (ImQMD) model. Chapter 4 is dedicated to the processing and calibrations of raw data collected by the neutron detection system, vLANA, which include hit position calibration, time-of-flight calibration, preprocessing of analog-to-digital converter (ADC) values, light output calibration, and pulse shape discrimination. Chapter 5 focuses on the reconstruction of neutron spectra. An in-depth background analysis is presented, detailing scattered neutron background and random background at radio frequency (RF). It then discusses neutron detection efficiencies, including geometric and intrinsic efficiency. The newly developed neuSIM4 code is introduced and used to simulate the light response of vLANA for an accurate determination of intrinsic efficiency. The chapter concludes with the neutron spectra for the Ca + Ni systems at 140 MeV/u, and updates the 16 studies on isoscaling and Bayesian analysis from Chapter 3. Chapter 6 provides a summary of the research discussed in this dissertation work. It encapsulates the key findings, discusses their implications, and offers suggestions for future research in the area. Following these chapters, a comprehensive bibliography is provided that lists all the literature and resources referred to in the writing of this dissertation. Lastly, appendices include further details that complement the main chapters. 17 CHAPTER 2 EXPERIMENT E15190-E14030 Two experiment proposals, E15190 and E14030, were accepted to conduct studies at the S2 Vault of the National Superconducting Cyclotron Laboratory (NSCL) in Michigan, USA, now known as the Facility for Rare Isotope Beams (FRIB). This combined experiment, referred to as E15190-E14030 or simply the experiment, ran from early February to late March 2018. vLANA (neutron detection system) Neutron Neutron wall A Microball wall B (impact parameter Veto detection system) wall LANA 𝜃lab = 39.37∘ beam line Forward array target HiRA10 (charged-particle detection system) Figure 2.1 Schematic top view plot for the experimental setup. Not to scale. During the experiment, data were collected in runs. Each run is stored as a binary file, numbered by a four-digit run number. Each run typically lasts 30 minutes, and is manually started and stopped by the experimentalists on shift. Storing data in terms of runs reduces the risk of losing a significant amount of data if the binary files become corrupted due 18 to unexpected issues with the data acquisition (DAQ) system or other computer software. Between runs, the experimental configuration may change, allowing us to explore the physics under different conditions and beam intensities. Runs 1000-1999 and 3000-3999 were reserved for testing and calibration. Runs 2000-2999 and 4000-4999 were reserved for real beam data collection, although testing and calibration were still done intermittently in some of these runs. The electronic run log (Elog) contains 40 more detailed information [97]. Runs 2000-2999 took place in February and used the Ca 48 beam; runs 4000-4999 took place in March and used the Ca beam. Both beam isotopes were delivered at two energies, 56 MeV/u and 140 MeV/u, giving us a total of 2 × 2 = 4 58 64 beam configurations. On the other hand, four isotope-enriched targets, namely Ni, Ni, 112 124 Sn, and Sn, were prepared, with their thicknesses listed in Table 2.1. Every target was 40 48 bombarded by Ca and Ca at both energies, giving us a total of 4 × 2 × 2 = 16 beam-target configurations. Throughout this thesis, we denote the beam-target pairs using the notation, “B + T”. For 40 example, Ca + 58Ni implies that the 40 Ca beam was used to collide with the 58 Ni target. Table 2.2 lists all 4 × 2 = 8 beam-target pairs, along with their corresponding (N − Z)/A values to indicate the asymmetry of the system. These beam-target pairs are often grouped target thickness [mg/cm2 ] 58 Ni 5.0 64 Ni 5.3 112 Sn 6.1 124 Sn 6.47 Table 2.1 Thicknesses of the four isotope-enriched targets. Ca + Ni Ca + Sn beam target (N − Z)/A beam target (N − Z)/A 40 58 40 112 Ca Ni 0.020 Ca Sn 0.079 40 64 40 124 Ca Ni 0.077 Ca Sn 0.146 48 58 48 112 Ca Ni 0.094 Ca Sn 0.125 48 64 48 124 Ca Ni 0.143 Ca Sn 0.186 Table 2.2 Beam-target pairs in this experiment and their corresponding (N − Z)/A values. 19 into two categories, Ca + Ni and Ca + Sn, depending on the target element. In this work, we only analyze the beam-target pairs in each category with the lowest and the highest 40 (N − Z)/A values within the category, namely, Ca + 58Ni, 48 Ca + 64Ni, 40 Ca + 112Sn, and 48 Ca + 124Sn, each at both beam energies. Beam-target pairs with intermediate (N − Z)/A values will be explored in future works. Figure 2.2 Experimental setup overview, with the chamber opened. LANA, which consists of Neutron wall A (NWA) and Neutron wall B (NWB), is not visible in this photograph. It is installed right behind the Veto wall (VW). Three detection systems were designed and constructed, namely, the charged-particle detection system (HiRA10), the neutron detection system (vLANA), and the impact parameter detection system (Microball). All detection systems will be discussed in detail in the subsequent sections. The relative positions of the detection systems are shown in Figure 2.1 20 as a schematic top view plot. Notice that both Microball and HiRA10 are placed inside a thin-wall aluminum chamber. During the experiment, the chamber will be pumped down to a pressure of ∼10−6 to ∼10−5 Torr to allow charged particles reaching the two detection systems with negligible interactions with the air molecules. The vLANA is composed of four components: Forward array (FA), Veto wall (VW), Neutron wall A (NWA), and Neutron wall B (NWB). The Forward array is placed inside the chamber, whereas the other three components are placed outside the chamber. The FA is used to provide a start time from charged particles emitted from the beam projectile, the VW is used to reject charged particles that reach the neutron walls, and the NWA and NWB are used to detect neutrons. Collectively, the NWA and NWB are also known as the Large Area Neutron Array (LANA) [98]. Figure 2.2 shows a photograph of the experimental setup. 2.1 Impact parameter detection system – Microball Microball was originally designed and constructed to be a 4π array with 9 concentric rings of CsI(Tl) scintillators compactly arranged around the target to optimize the detection of light charged particles emitted from beam-target interactions [99]. In this experiment, the Microball is the detector system with the largest solid angle coverage. It is used to measure the charged-particle multiplicity, Nc , which is then used to quantify the centrality of collisions in terms of impact parameter b. Impact parameter is defined to be the transverse distance between the centers of two particles, the beam isotope and the target isotope. We shall discuss how the relation b = b(Nc ) is established in detail in subsection 2.1.2. The determination of b is crucial because the centrality of collisions is known to play a role in the determination of nuclear equation of state [86, 100, 101, 102]. Central collisions (b ≈ 0) tend to result in more violent emission than peripheral collisions (b ≈ Rbeam + Rtarget ). In other words, b is expected to have a negative correlation with the number of particles emitted N . In Figure 2.3, a schematic side view of the Microball is drawn. The beam travels from right to left, which aligns with the +z-axis of the Microball. The most forward ring is Ring 21 Figure 2.3 Schematic side view of the original Microball design, adapted from [99]. Some rings and CsI(Tl) scintillators are removed for this experiment. See subsection 2.1.1 for details. 1, spanning θ = 9.0 ± 5.0◦ ; the most backward ring is Ring 9, spanning θ = 159.0 ± 12.0◦ . The scintillators on each ring compactly span the entire azimuthal range of ϕ ∈ [0, 2π]. All scintillators are facing toward the center of the target, marked by the “+” icon in the figure. The gap between neighboring rings is 0.10 mm, whereas the gap between neighboring scintillators on the same ring is 0.12 mm. The front surface of each scintillator is covered with a thin layer of aluminum foil at a thickness of 0.5–0.9 mg/cm2 . This reflective layer serves to protect the scintillator from external light and to increase internal reflection. Additionally, due to the high number of delta electrons emitted in heavy-ion collisions [103], Sn/Pb foils are used to cover the scintillators in the forward rings. The number of scintillators is different for each ring. Table 2.3 tabulates the geometric parameters of the Microball, which is capable of mounting 96 scintillators in total. In each ring, the scintillators are arranged compactly, with the first scintillator numbered as 22 scintillator 1, centered at ϕ = 0◦ (top). More specifically, for a ring with n scintillators, each scintillator has a δϕ coverage of 360◦ /n, centered at 360◦ 360◦ ϕ1 = 0◦ , . . . , ϕi = (i − 1) · , . . . , ϕn = (n − 1) · . (2.1) n n Ring No. 1 2 3 4 5 6 7 8 9 No. of scintillator 6 10 12 12 14 14 12 10 6 Distance from center [mm] 110 80 60 50 50 50 45 47 50 Min. θ [deg] 4.0 14.0 28.0 44.0 60.0 80.0 100.0 123.0 147.0 Max. θ [deg] 14.0 28.0 44.0 60.0 80.0 100.0 123.0 147.0 171.0 Table 2.3 Geometric parameters of Microball. Adapted from Table 1 in [99]. 2.1.1 Adaptation for the experiment Some rings and scintillators are removed for this experiment. Ring 1 (most forward) and Ring 9 (most backward) are removed to avoid heavy-ion beams directly hitting the Microball. Ring 6 is removed to allow target ladder to be installed. Scintillators 3, 4, 5, 6 in Ring 3, scintillators 4, 5, 6 in Ring 4, and scintillators 4, 5, 6 in Ring 5 are removed to allow charged particles emitted from the target to hit HiRA10 directly. Scintillator 4 in Ring 8 is excluded in the analysis due to bad performance. The final configuration of Microball is shown in Figure 2.4, leaving us 59 CsI(Tl) scintillators in total. Figure 2.4 also reveals the experimental hit pattern of Microball from an actual data run in spherical coordinates. As expected from kinematics, there are more hits in the scintillators at forward angles (small θ) than at backward angles (large θ). Along each ring (constant θ), the number of hits stays roughly the same due to azimuthal symmetry. Notice that there is no way of differentiating the hit positions within the same scintillator. An opening was constructed on the Microball in the range of ϕ ∈ [70◦ , 150◦ ] and θ ∈ [20◦ , 80◦ ] to allow direct access of charged particles from the target to HiRA10. However, some regions on the HiRA10 were obscured by the Microball scintillators. As a result, data collected from these obscured regions of HiRA10 were omitted from the analysis. A detailed discussion 23 180 1e6 160 1.0 140 R8-01 R8-02 R8-03 R8-05 R8-06 R8-07 R8-08 R8-09 R8-10 0.8 120 R7-01 R7-02 R7-03 R7-04 R7-05 R7-06 R7-07 R7-08 R7-09 R7-10 R7-11 R7-12 100 0.6 θ (deg) 80 HiRA10 R5-01 R5-02 R5-03 R5-07 R5-08 R5-09 R5-10 R5-11 R5-12 R5-13 R5-14 0.4 60 R4-01 R4-02 R4-03 R4-07 R4-08 R4-09 R4-10 R4-11 R4-12 40 R3-01 R3-02 R3-07 R3-08 R3-09 R3-10 R3-11 R3-12 0.2 20 R2-01 R2-02 R2-03 R2-04 R2-05 R2-06 R2-07 R2-08 R2-09 R2-10 0 0 50 100 150 200 250 300 350 0.0 φ (deg) Figure 2.4 Microball hit pattern in run 4520. Each rectangle represents a scintillator. Ring number and scintillator number are labeled in the format of “RX -Y ”, e.g. “R2-10” means scintillator 10 of Ring 2. about HiRA10 is provided in section 2.2. Figure 2.6 and Figure 2.7 display photographs of the Microball used in this experiment, clearly showing the opening for HiRA10. The pulse-shape signals collected from a Microball scintillator consists of two components, fast and tail, as shown in Figure 2.8. The red dashed line indicates the alpha line before the punch-through at Tail ≈ 3400 ADC. A punch-through occurs when the energy of a particle is so high that it penetrates the entire scintillator. The fast component and the tail component refer to integrations of the pulse over short and long time intervals, respectively. It is a common technique to use two pulse-shape components from CsI(Tl) scintillators for particle identification (PID) [104, 105]. More discussion about punch-through and PID will given in subsection 2.2.2, where we show PID plots with better resolution from HiRA10. Before analyzing the Microball multiplicity count for impact parameter determination, a manual acceptance cut is made for every CsI(Tl) scintillator on the MicroBall and every beam-target configuration. The cut is drawn using a green solid line in Figure 2.8. Hits outside of the acceptance cuts are regarded as noise, and they have been excluded from the 24 Figure 2.5 Simulated hit pattern of Microball and HiRA10. Some overlapping regions between Microball and HiRA10 can be seen in R2-03, R2-04, and R5-07. Microball multiplicity count. A zoomed-in plot at the lower ADC region has been given in Figure 2.9 to better show the PID lines of light charged particle. 2.1.2 Impact parameter determination In heavy-ion collision studies, the impact parameter b plays a crucial role. Collisions occurring at small impact parameters (b ≈ 0) are referred to as central collisions, and are known to produce more fragmentation than peripheral collisions, which occur at large impact parameters (b ≈ Rbeam + Rtarget ). Therefore, the impact parameter b is expected to be negatively correlated with the number of particles emitted N , i.e. the smallest impact parameter correspond to the largest multiplicity. Denote Nc as the observed multiplicity of Microball. The subscript c is used to signify vast majority detected particles are charged particles. To establish the relation b = b(Nc ), there are several assumptions that need to be made. First, we ignore the fluctuation of N and Nc at any given b, so all quantities related to N and Nc are just expressed as the expected values from now on. Second, the experimentally observed multiplicity Nc is positively correlated 25 Figure 2.6 A photograph of the Microball taken after the teardown of the experiment. Some rings and scintillators are placed back to their original positions before storage, e.g. Ring 6. Figure 2.7 An opening is made on the Microball for HiRA10. Notice that Ring 6 is removed. to the actual N . This allows us to infer b as a function of Nc rather than N , which is not accessible experimentally. Third, b(Nc ) is a monotonic function. The first assumption is a bold one because the fluctuation of N or Nc at fixed b is non-negligible due to the stochastic nature of collisions, as demonstrated by many transport models [106, 107, 108, 109]. In the data analysis, we would only gate on a wide range of b values, which allows for more room of 26 Figure 2.8 Particle identification of Microball scintillator, generated from run 4660, Microball Ring 2, scintillator 2. The green solid line is the acceptance cut. The red dashed line indicates the alpha line. The lines for proton, deuteron, and triton are below the alpha line, and they are labeled in Figure 2.9. uncertainty. The Microball is used to measure the charged-particle multiplicity, which is assumed to decrease monotonically with impact parameter b. The probability of producing particles with multiplicity less than or equal to Nc is proportional to the geometrical area of a disk with radius b(Nc ), we can calculate the deduced impact parameter, b̂ = b/bmax , which is a dimensionless variable associated with the centrality of the collision. The probability of producing particles with multiplicity less than or equal to Nc is 58 64 112 124 Ni Ni Sn Sn 40 Ca at 56 MeV/u 7.15 7.22 8.95 8.60 40 Ca at 140 MeV/u 6.76 6.89 8.43 8.22 48 Ca at 56 MeV/u 6.93 7.16 8.72 8.36 48 Ca at 140 MeV/u 5.96 6.97 8.31 8.05 Table 2.4 The bmax for all 16 reaction systems in the unit of fm. The rows represent the target isotopes, whereas the columns represent the beam isotopes as well as their energies. 27 Figure 2.9 A zoomed-in version of Figure 2.8. The red dashed lines indicate the alpha, triton, deuteron, and proton lines, from top to bottom. proportional to the geometrical area of a disk with radius b(Nc ), σ(Nc ) = π(b(Nc ))2 . Denote P (Nc ) as the probability of detecting Nc charged-particles on the Microball, then simple geometry tells us v u ∞ uX b(Nc ) ∝ t P (N ) . (2.2) N =Nc The probability P (Nc ) can be estimated according to C(Nc ) P (Nc ) = P∞ , (2.3) Nc =1 C(Nc ) where C(Nc ) is the number of events with Microball multiplicity Nc in the experimental data. To finalize the impact parameters, we introduce the reduced impact parameter b̂, defined as b(Nc ) b̂ ≡ , (2.4) bmax where bmax is the largest impact parameter considered in data analysis, and its values for different reaction systems are tabulated in Table 2.4. Head-on collisions happen at b̂ = 0, and they are expected to produce the largest multiplicity; peripheral collisions happen at 28 b̂ ≈ 1, and they are expected to produce the smallest multiplicity. In this dissertation, we are interested in central collisions with b̂ ≤ 0.4. 8 8 Ebeam = 56 MeV/u Ebeam = 140 MeV/u 7 40 Ca + 58 Ni 7 40 Ca + 58 Ni Impact parameter b [fm] 6 40 Ca + 64 Ni 6 40 Ca + 64 Ni 48 Ca + 58 Ni 48 Ca + 58 Ni 5 48 Ca + 64 Ni 5 48 Ca + 64 Ni 4 4 3 3 2 2 1 1 0 5 10 15 20 25 0 5 10 15 20 25 Microball multiplicity Nc Microball multiplicity Nc 10 10 Ebeam = 56 MeV/u Ebeam = 140 MeV/u 40 Ca + 112 Sn 40 Ca + 112 Sn 8 8 Impact parameter b [fm] 40 Ca + 124 Sn 40 Ca + 124 Sn 48 Ca + 112 Sn 48 Ca + 112 Sn 6 48 Ca + 124 Sn 6 48 Ca + 124 Sn 4 4 2 2 0 5 10 15 20 25 0 5 10 15 20 25 Microball multiplicity Nc Microball multiplicity Nc Figure 2.10 Impact parameter b as a function of the charged-particle multiplicity Nc for all 16 reaction systems. In Figure 2.10, we show the impact parameter b as a function of Nc for all 16 reaction systems. The top two panels are Ca + Ni systems; the bottom two panels are Ca + Sn systems. The left two panels are for 56 MeV/u; the right two panels are for 140 MeV/u. This plot is replotted based on the thesis work by S. Sweany [110]. In the analysis, a selection criterion is applied to Nc , effectively creating a cut on the impact parameter values considered. 29 2.2 Charged-particle detection system – HiRA10 The High Resolution Array (HiRA) is an array of detector telescopes that was first designed and developed for charged-particle detection in rare isotope beam experiments [111] with the capability of particle identification (PID). In this experiment, the original design was upgraded to bring forth the “HiRA10” [110, 112]. The number “10” in the new name came from the fact that we replaced the 4-cm long CsI(Tl) crystals in the original HiRA telescope by some 10-cm long CsI(Tl) crystals, increasing its capability to fully absorb the energy of charged particles with higher energies (≳ 100 MeV/u). Figure 2.11 A schematic plot showing the internal components of the HiRA10 telescope. The charged particles would enter from the top, penetrating a double- sided silicon strip detector (not shown) before reaching the CsI(Tl) crystals. Adapted from [110]. Figure 2.11 shows a schematic plot of the HiRA10 telescope. The front surface is a 32 × 32 double-sided silicon strip detector (DSSD) with a thickness of 1.5 mm. Each strip on the DSSD, both front side and back side, has a width of 1.95 mm, offering an excellent position resolution of 0.347◦ ± 0.002◦ . Details about the performance of DSSD can be found in [111, 113]. After the DSSD, four CsI(Tl) crystals of length 10 cm are packed into four 30 quadrants. The front surface of the crystal has a dimension of 35 mm × 35 mm, whereas the back surface has a dimension of 44.6 mm × 44.6 mm. Details about the design and construction of HiRA10 can be found in [110]. 2.2.1 Configuration for the experiment HiRA10 telescope is designed to be used in a variety of experiments for charged-particle detection. While each telescope can be used independently as a single detector, multiple telescopes can be combined to form a detector array with a larger coverage and one coherent data acquisition system. This modular design allows experimenters to easily construct a detector array of various geometry and size that can satisfy the experimental constraints and needs. Figure 2.12 A photograph of HiRA10 array taken during the preparation of the experiment. There are a total of 12 telescopes, arranged in three vertical columns, called “towers”, each consists of four telescopes. In this experiment, 12 telescopes with identical designs are used to construct the HiRA10 array as shown in Figure 2.12. The telescopes are arranged into three vertical towers, each consists of four telescopes. The front surfaces of the telescopes are oriented to face the target normally, with a distance ranging from 32.9 cm to 33.3 cm. While the actual positions of the telescopes are not precisely planned in advance, the positions are fixed rigidly throughout 31 the entire experiment campaign. ROMER arm is then used to measure the exact positions for every telescope. Specifically, the ROMER arm measures coordinates for all the four precision-machined dimples at the front frame of each telescope. The coordinates for the pixels on the DSSD can then be calculated based on the dimple positions and the known geometry of the telescope. 160 tele-03 140 tele-02 tele-07 120 tele-06 tele-11 φ [deg] tele-01 tele-10 100 tele-05 tele-09 tele-00 tele-04 80 tele-08 6020 30 40 50 60 70 80 θ [deg] Figure 2.13 Coverage of HiRA10 array. The azimuthal angle ϕ starts from the top and increases clockwise when viewed along the beam direction; this is the same convention used in Figure 2.4 for the Microball. The polar angle coverage of the array starts from θ ≈ 20◦ to θ ≈ 80◦ with respect to the beam direction defined as the +z-axis. The azimuthal angle coverage ϕ starts from ϕ ≈ 70◦ to ϕ ≈ 160◦ ; ϕ = 0 is defined as the direction pointing upward in the lab (the +x-axis), and it increases clockwise when viewed along the beam direction. Figure 2.13 shows the coverage of the HiRA10 telescopes in spherical coordinate system. 2.2.2 Particle identification of charged particles The design of HiRA10 telescope places a DSSD in front of the CsI(Tl) crystals to allow for particle identification (PID) using the “∆E-E technique” [114, 115, 116]. The DSSD serves as the ∆E detector, while the CsI(Tl) crystals serves as the E detector. The technique assumes that a charged particle will penetrate the ∆E detector before stopping entirely in 32 the E detector. If a particle is so energetic that it also penetrates through the E detector, we call it a punch-through. In either case, both ∆E and E detectors will measure some energy losses of the particle, which depend on the mass, electric charge, and incident energy. This provides a basis for particle identification. A commonly used model to describe the energy losses of charged particles is the Bethe- Bloch equation [117]. For an incident particle with energy E and electric charge q, the energy loss per unit length when passing through a medium, i.e. the electronic stopping power, is given by 2   2me c2 β2    dE 4πρe αℏq 2 − = ln · −β . (2.5) dz me βe Eex 1 − β2 Here, β is the velocity of the particle in units of the speed of light c, and it is related to the particle mass m and energy E by s  2 mc2 β= 1− . (2.6) E The medium matter is parametrized by the electron number density ρe and the mean excitation energy Eex . The rest of the constants are the electron mass me , the fine structure constant α, the reduced Planck constant ℏ, and the elementary charge e. As shown in Figure 2.14, the stopping power of a medium increases with decreasing energy. In other words, when a particle enters a medium, it loses energy at an increasing rate. The particle will eventually stop if the medium is sufficiently thick, otherwise punch-through will occur. For given medium and kinetic energy, the stopping power is a function of the particle mass m and electric charge q. Particularly, the stopping power is proportional to q 2 , which results in two isotope groups in Figure 2.14 – one for hydrogen and one for helium. Within each group, the lines are separated with smaller gaps by different masses. This separation mechanism allows for particle identification. In this experiment, the 1.5 mm-thick DSSD (∆E detector) is thin enough to allow most particles to penetrate it and enter the CsI(Tl) crystals (E detector). The E detector has a thickness of 10 cm and a higher electron number density than the ∆E detector. As a result, 33 proton helium-3 Electronic stopping power −dE/dz [MeV/mm] 102 deuteron helium-4 triton helium-6 101 100 101 102 103 Energy E [MeV] Figure 2.14 Electronic stopping powers of light isotopes in silicon according to the Bethe-Bloch equation defined in Equation 2.5. Figure 2.15 Particle identification plot of HiRA10 array. Adapted from [112]. many charged particles would stop completely in the E detector, depositing all remaining energy into the crystals. In Figure 2.15, the vertical axis (∆E) and the horizontal axis (E) can be compared to the respective axes in Figure 2.14. In our data analysis, the measured values from both ∆E and E detectors have to be added together to obtain the total energy 34 loss of the particle in a HiRA10 telescope. If the particle does not punch-through the E detector, then this total energy loss is equal to the incident energy of the particle. particle 1.5 mm Si + 10 cm CsI(Tl) 1.5 mm Si only 10 cm CsI(Tl) only p 199.56 15.539 198.29 d 132.46 10.447 131.63 t 104.47 8.252 103.82 3 He 237.19 18.284 235.67 4 He 199.17 15.509 197.90 6 He 156.74 12.325 155.75 Table 2.5 Punch-through energies in MeV/A for light charged particles computed using LISE++ [118, 119]. The second column presents energies required to first penetrate a 1.5 mm-thick DSSD (∆E detector), then a 10 cm-thick CsI(Tl) crystals (E detector). The third and fourth columns respectively show energies necessary to fully penetrate solely a 1.5 mm-thick Si crystal and a 10 cm-thick CsI(Tl) crystal. Incident angles are set to be normal to the detector surfaces. If the particle does punch-through the E detector, then only a portion of the remaining energy after the ∆E detector is deposited into the E detector. As the particle energy increases beyond the punch-through energy, the stopping power of the E detector decreases. This explains the bending of PID lines after the punch-through energies that can be seen in Figure 2.15. The punch-through energies of six light charged particles are listed in Table 2.5. Much effort has been put into calibrating the energy measured by the DSSD and CsI(Tl) crystals of HiRA10, as well as correcting for punch-through particles and reaction losses as described in [110, 112, 120, 121], and will not be repeated here. 2.3 Neutron detection system – vLANA Neutron detection is performed by the vLANA (LANA with Veto Wall), which consists of Forward array (FA), Veto wall (VW), Neutron wall A (NWA), and Neutron wall (B). We refer readers to Figure 2.1 for a schematic plot of the setup. The two neutron walls, NWA and NWB, are collectively referred as the Large Area Neutron Array (LANA) in [98, 122]. LANA is used to collect scintillation lights due to neutrons, FA is used to measure the time-of-flight of neutrons, and VW is used to veto charged particles. The primary goal of this experiment is to construct neutron spectra over the range from 35 20 MeV to 200 MeV. To achieve this, we need to measure two things accurately: kinetic energies and neutron counts. To measure the kinetic energies accurately, we need to establish good measurement for the time-of-flight (ToF). Once ToF t is measured, the kinetic energy K can be calculated using relativistic kinematics,    1  K = mn c2  r  r 2 − 1 ,  (2.7) 1− ct where r is “distance-of-flight”, c is the speed of light, and mn is the neutron mass. The timing resolution in this experiment is about δt ≈ 0.2 ns [122, 123]. The energy resolution δK depends on both K and r. With everything else being fixed, δK increases with increasing K, but decreases with increasing r. To achieve the desired energy resolution of δK ≈ 5 MeV at K = 200 MeV, we need to place the neutron detector at least ≳ 365 cm away from the target according to " 2 #3/2 δt K r≥ · mn c3 1+ −1 . (2.8) δK mn c2 In this experiment, NWA and NWB are placed at 517.5 cm and 441.6 cm away from the target, with energy resolutions of 3.5 MeV and 4.5 MeV, respectively. The time-of-flight (ToF) is measured by comparing the time difference between the arrival times at LANA and FA for each collision event at the target, with FA serves as a reference time for the start time and LANA serves as the stop time. While FA is not the actual start time of the collision event, it is a good approximation for the start time because the FA is placed close to the target (∼ 12 cm) when compared to the distance between the target and NWB (∼ 440 cm). Furthermore, with some assumptions made about how the FA is triggered in each event, the ToF shift due to the difference between the actual start time and the FA start time can be corrected using the prompt gamma-ray signals. More discussions will be given in subsection 2.3.2 and section 4.2. 36 To measure the neutron counts accurately, we need to optimize the purity of neutrons and calculate the efficiency of detection. Figure 2.16 Top panel shows a two-dimensional time-of-flight spectrum of the neutron wall B (NWB) without vetoing, in which the charged particle lines for protons (p), deuterons (d), and tritons (t) are clearly visible; the triangle symbols in magenta mark their respective punch-through points. Bottom panel shows same spectrum but after charged-particle vetoing. Prompt gamma-ray is also visible in both panels at around 15 ns in time-of-flight. Purity of neutrons is achieved by rejecting the unwanted particles aside from neutrons. In heavy-ion collisions, emitted particles that have a chance to exit the vacuum chamber and travel through air until they are detected by LANA include neutrons, gamma rays, and charged particles such as protons, deuterons, tritons, etc. Electrons are excluded because they 37 are easily absorbed by the chamber wall as well as the atmosphere. To reject gamma rays, NE-213 organic liquid scintillation material was chosen to fill the LANA structure [124, 98]. To reject charged particles, VW is placed in between the target and LANA, allowing the use of anti-coincidence technique for tagging the unwanted charged particles. Specifically, if both LANA and VW are triggered for a particle, the particle is rejected as a charged particle. In Figure 2.16, we show the ToF spectra before and after charged-particle vetoing. Both the x-axis (time-of-flight) and y-axis (light output) have been calibrated. Detailed description of the calibration procedures will be given in chapter 4. The efficiency of detection is the ratio of the number of detected neutrons to the number of emitted neutrons. The determination of efficiency is often broken down into two parts: the geometric efficiency and the intrinsic efficiency. Geometric efficiency is associated with the geometric acceptance of the detector, whereas intrinsic efficiency is associated with the reaction processes in the scintillation material. They are discussed in more detail in subsection 5.2.2 and subsection 5.2.3, respectively. 2.3.1 Large area neutron array Large area neutron array (LANA) consists of two “neutron walls” (NWs), NWA and NWB. NWA and NWB are placed along the direction of 39.37◦ with respect to the beam direction, at distances of 517.5 cm and 441.6 cm, respectively, from the target. Both NWs face the target perpendicularly and share identical dimensions aside from the placement. In this dissertation work, we primarily analyze the data from NWB. Each NW is made up of 25 horizontal bars stacked one on top of the other, with a gap of 3/8 inch (0.9525 cm), as shown in Figure 2.17. Each bar is made up of a hollow Pyrex container with a thickness of 1/8 inch (0.3175 cm). Inside each container, liquid scintillator NE-213 is used to fill up the cuboid cavity which has a transverse height of 3.0 inches (7.62 cm), a longitudinal thickness of 2.5 inches (6.35 cm), and a length of 76.0 inches (193.04 cm). This setup gives a detection area of ∼ 2 m × 2 m. When a passing radiation such as gamma ray, charged particle or neutron excites the 38 Figure 2.17 A picture of a neutron wall (NW) with the cover removed. There are 25 bars in total, numbered from 0 to 24 inclusively from bottom to top. Each bar is attached with two photomultiplier tubes (PMTs) at both ends. The bluish glow of the NW bars results from the illumination of an ultraviolet light source. Photo courtesy of K. Zhu. scintillator, light output is produced depending on the radiation type and energy. To collect the light outputs for further analysis, two photomultiplier tubes (PMTs) are installed at the two ends of each bar. Emitted photons excited by the radiation will travel through the scintillation material, reflect off the inner surfaces of the Pyrex container, and then be collected by the PMTs. On the outer surfaces, a thick black paper are attached between neighboring bars to absorb any escaping photons, preventing them from entering the neighboring bars. Light output signals collected from each PMT are processed to give three quantities: the signal arrival time, the Fast charge integration and the Total charge integration. To achieve this, the output signals from every PMT are processed according to an electronics setup shown in Figure 2.18. Each PMT contains a photocathode, a series of dynodes, and an anode. The photocathode 39 Figure 2.18 Schematic of the electronics setup for the neutron wall (NW) for one PMT from one bar. Adapted from [123, 122]. absorbs photons and releases electrons. The released electrons are accelerated and multiplied in number by the dynodes. The anode collects electrons from the last dynode and converts them into a voltage pulse. Each PMT is equipped with a passive voltage divider circuit [98] to allow readouts from both the anode and the last dynode. The anode signal is used to quantify the light output, whereas the last-dynode signal is used to determine the arrival time of the light output. The last-dynode signal, or simply the dynode signal, is sent to a LeCroy 1806 16-channel constant-fraction discriminator (CFD) [125] to determine the timing from the leading edge of the pulse signal. The CFD output is duplicated into three copies: one to create a logical OR with of other CFD outputs from other PMTs, one to generate gate signals for the Fast charge integration, and one to the CAEN V1190A time-to-digital converter (TDC) [126]. The trigger signal is then generated by taking the logical OR signal from all PMTs in anti-coincidence with a global busy signal, which is generated as an OR between the computer busy and the 40 fast busy. This trigger signal is duplicated into multiple copies, and one of the copies is sent to the TDC as a reference time. The TDC takes in the individual CFD outputs and records their time differences with respect to the reference time from the trigger, giving an arrival time for each PMT. These arrival times need to be calibrated to infer hit position and time-of-flight (ToF), which are discussed in section 4.1 and section 4.2, respectively. The anode signal is sent to an in-house passive splitter at NSCL. The splitter splits the signal into two equal copies: one to the CAEN V862 charge-to-digital converter (QDC) [127] for Total charge integration, and one to the CAEN V792 QDC [128] for Fast charge integration. Both modules have 32 channels and 12-bit resolutions. The Total charge integrates over the entire pulse (340 ns), whereas the Fast charge integrates over the fast component (first 30 ns) only. The Total charge integration uses a common gate of width ∼ 3 µs generated by the trigger signal, whereas the Fast charge integration uses an individual gate generated by the individual CFD output in anti-coincidence with the global busy. In section 4.3, we discuss corrections for these analog-to-digital (ADC) values from the QDC modules to account for non-linearity and saturation. Collecting both the Fast and Total charge integrations allows us to perform neutron-gamma pulse shape discrimination (PSD) [129, 98, 122] that is discussed in section 4.5. The Total charge integration is also used to calibrate the light output in MeV electron equivalent (MeVee) units, and it is discussed in section 4.4. We have described the setup for the analog electronics. In this experiment, a digital electronics system [130, 129] is also used in parallel for benchmarking and development purposes. Instead of recording only the Fast and Total signals, the digital electronics system the entire pulse at a sampling rate of 500 MHz, i.e. 2 ns per sample. PSD using digitized pulses has been studied in [129]. The improved PSD algorithm has been adopted for analog signals in this experiment, which will be discussed in section 4.5. 2.3.2 Forward array The primary purpose of the forward array (FA) is to measure a reference time for the event start time. The event start time is then used as the zero point for the time-of-flight 41 (ToF) determination of LANA. Figure 2.19 A picture of the Forward array (FA). The FA consists of 18 annulus sectors, numbered using marker pen from 1 to 16; the remaining two wider sectors are not labeled. Photo courtesy of J. Barney. The forward array is made up of 18 annulus sectors, 16 of which have an angle of 18◦ and 2 of which have an angle of 36◦ , as shown in Figure 2.19. The 18◦ -sector has an inner diameter of 1 inch and an outer diameter of 4.5 inches; the 36◦ -sector has identical inner diameter but a shorter outer diameter of 3.0 inches. On each sector, an NE-110 plastic scintillator is mounted for charged-particle detection. Each scintillator is painted with Bicron BC-620 for internal reflection, wrapped with aluminized Mylar foil to reduce crosstalk, and attached to a Hamamatsu R5600U PMT with an E5780 base for light collection. During the experiment, the FA is placed ∼ 12 cm downstream from the target inside the chamber, with beam axis passing through the central hole. It is the only component of the vLANA system that is placed inside the vacuum chamber. When a heavy-ion collision event 42 occurs at the target, a few dozens of charged particles will be emitted. Due to the proximity of the FA to the target, many of the FA scintillators will be hit by these charged particles and produce light signals. The average number of triggered scintillators per event is ∼ 12.7, with over 90% of all events triggering at least five FA scintillators. The arrival time for each FA scintillator can then be determined. The FA start time is calculated as the average of all the arrival times of the FA scintillators after excluding the extreme values. The FA start time does not equal to the actual event start time that happens at the target. Instead, the FA start time tFA > 0 should happen after the event start time t ≡ 0. Nonetheless, it is reasonable to assume that the values of tFA are distributed with a small variance from event to event in a given run. In other words, we assume the FA start time is delayed by a constant amount from the event start time. This allows us to correct for the shift by studying the prompt gamma ray arrival time at the LANA. Prompt gamma ray are emitted promptly after a heavy-ion collision, and the difference between the emission time and the event start time is negligible in this experiment. section 4.2 will be dedicated to the discussion of the time-of-flight calibration, which includes the aforementioned correction. 2.3.3 Veto wall The Veto wall (VW) is designed to veto out the charged particles. Its construction is a collaborative effort between Western Michigan University (WMU) and Michigan State University (MSU). Section 2.5 in [122] provides a detailed description of prototyping and construction. 43 Figure 2.20 A picture of the Veto wall (VW). The VW consists of 25 vertical plastic scintillator bars. PMTs are installed at both ends of each bar. Photo courtesy of K. Zhu. The VW consists of 25 vertical plastic scintillator bars stacked side by side in a zigzag pattern, as shown in Figure 2.20. VW is placed along the direction of 39.37◦ with respect to the beam axis, at a distance of 393.3 cm from the target. This placement of the VW is between the target and the NWB. Furthermore, there is an overlapping width of 3 mm for every pair of adjacent VW bars. This design makes any particles emitted from the target to be blocked from directly hitting the LANA without passing through the VW. Notice that the VW bars are installed vertically, whereas the NW bars are installed horizontally. This resembles the design of a double-sided silicon strip detector (DSSD), which 44 allows the hit position of a particle to be determined in both the horizontal and vertical directions. In section 4.1, we will discuss how VW is also used to calibrate the hit position of the NWB. Each VW bar is made up of EJ-200 (1.023 g/cm2 ) plastic [131], with a length of 250 cm, a transverse width of 9.4 cm, and a longitudinal thickness of 1 cm. Both ends of a VW bar are coupled with PHOTONIS XP3462 PMTs for light collection. The thickness of 1 cm is chosen such that it is sufficiently thick for charged particles passing through the VW to always induce a detectable amount of light. To give an example, a 200 MeV proton would lose ∼ 4.6 MeV in the scintillator according to the Bethe-Bloch equation [117]. On the other hand, the thickness cannot be made arbitrarily thicker because the interactions between neutrons and the VW need to be kept to a minimum. 45 CHAPTER 3 RESULTS FROM CHARGED PARTICLES This chapter focuses on the outcomes derived from the charged-particle detection system, HiRA10, and the impact-parameter detection system, Microball. The subsequent chapters will delve into the data analysis and findings from the neutron detection system, vLANA. We will first discuss the construction of charged-particle spectra, including energy spectra in the lab frame and transverse momentum spectra. The subsequent section discusses the isoscaling results, which are examined using the transverse momentum spectra. Finally, we explore the construction of neutron-proton (n/p) ratios based on isoscaling, leveraging them to investigate the nuclear equation of state (EOS) parameters and nucleon effective masses through transport model calculations and Bayesian analysis. 3.1 Charged-particle spectra This thesis focuses on the spectral yields of light charged particles, including protons (p), deuterons (d), tritons (t), 3He, and 4He, resulting from the central collisions of 40 Ca + 58Ni, 48 Ca + 64Ni, 40 Ca + 112Sn, and 48 Sn + 124Ni at both 56 MeV/u and 140 MeV/u. Isotopes of higher mass numbers have been excluded due to their low statistics detected by HiRA10. To maintain brevity of this chapter, we only include plots for a selected combination of collision systems and beam energies directly in this section. The comprehensive set of plots is presented in Appendix A. 3.1.1 Energy spectra in lab frame We start our discussion with the energy spectra in the lab frame. The particles are identified employing the ∆E-E method, as delineated in subsection 2.2.2. Once identified, the particles are deployed to populate distinct two-dimensional histograms according to their types. Each histogram is characterized by the lab (kinetic) energy Elab as x-axis and lab theta θlab as y-axis. The energy of these particles is deduced from their respective energy losses in the detectors, whereas the theta angles are estimated from the hit positions on the 46 HiRA10 detectors. In order to account for fluctuations in beam intensity during the experiment, we normalize the lab energy spectrum to the number of analyzed events. In this study, we concentrate on central collisions, defined as events with reduced impact parameters b̂ ≤ 0.4. This corresponds to b ≲ 3 fm for all the reaction systems discussed in this chapter. The actual impact parameter cut is implemented in terms of Microball multiplicity, as elaborated in subsection 2.1.2. For each reaction system, we define Nevent as the total number of central events, which serves as an important normalization factor for various spectra in this thesis. The final step in obtaining the energy spectrum is to correct for the solid angle coverage and detector efficiency of the HiRA10 array. The energy spectrum, expressed as a double- differential multiplicity that depends on energy and angle, is estimated from a two-dimensional ′ ′ histogram. For a specific bin at (Elab , θlab ), we calculate the differential multiplicity as follows: d2 M(A,Z) N(A,Z) (i′ , j ′ )   ′ ′ 1 1 1 (E , θ ) ≈ · · · . (3.1) dElab dΩlab lab lab Nevent ϵ(θlab′ ′ ) 2π sin θlab ∆Elab · ∆θlab ′ Here, ϵ(θlab ) ∈ [0, 1] is the overall detector efficiency, and N(A,Z) (i′ , j ′ ) is the count of particles of type (A, Z) in the bin indexed by i′ and j ′ , which correspond to Elab ′ ′ and θlab , respectively. The histogram uses uniform bins with widths of ∆Elab = 1 MeV/A and ∆θlab = 0.1◦ . Owing to the limited coverage and dynamic range of the detectors, the resulting two- dimensional phase space is incomplete, leading to rectangular phase spaces as illustrated in Figure 3.1. Each panel in the figure corresponds to a specific particle type, including p, d, t, 3He, 4He, and 6He. Within this confined phase space, it is conspicuous that the yield diminishes as Elab or θlab increases. The phase space delineated by four boundaries is determined by the constraints of the experiment. The horizontal boundaries stem from the limited θ coverage from 30◦ to 75◦ by HiRA10. The vertical boundaries result from the cuts in kinetic energies of the particles, with the lower one representing the energy threshold and the upper one symbolizing the particle punch-through. These energy cuts vary for different particles, as enumerated in Table 3.1. Note that the punch-through values stated for p, d, and t are slightly conservative, being 47 Figure 3.1 Two-dimensional phase space of lab θ versus lab energy for light charged particles, including p, d, t, 3He, 4He, and 6He, from a 40Ca + 58Ni collision at 140 MeV/u. energy threshold maximum energy detectable particle [MeV/A] [MeV/A] p 20.0 198.0 d 15.0 131.5 t 12.0 104.0 3 He 20.0 200.0 4 He 18.0 200.0 6 He 13.0 200.0 Table 3.1 Energy ranges for light charged particles. lower than the precise punch-through energies in Table 2.5. For the helium isotopes, we have arbitrarily set the upper boundary to be 200 MeV/u as their yields deplete long before reaching their respective punch-through energies (see Appendix A). 3.1.2 Transverse momentum spectra In analysis of heavy-ion collisions (HICs), it is often useful to distinguish nucleons into participants and spectators, as discussed in section 1.2. The interaction dynamics of these 48 participant nucleons, actively involved in the collision, provide significant insights into nuclear EOS and nucleon effective masses under high density conditions. These nucleons engage in a squeeze-out motion, ejecting from the reaction plane, contributing notably to the transverse momentum of the resulting emission [132, 133, 134]. This motivates the construction of transverse momentum spectra. Similar to the lab θ-E histogram, a two-dimensional histogram detailing transverse momentum, pT = p sin θ , (3.2) and rapidity,   1 E + p cos θ y = ln , (3.3) 2 E − p cos θ offers a complete kinematics of the emitted particles. In these expressions, all quantities are in the lab frame. The variable p denotes the magnitude of the momentum vector, while θ represents the angle between the momentum vector and the beam direction. The total energy p of the particle, denoted by E, is calculated using the identity E = p2 + m2 , where m is the mass of the particle. The two-dimensional pT -y phase space, akin to the lab θ-E phase space, is confined due to the finite coverage and dynamic range of the detectors. As depicted in Figure 3.2, the phase space presents an annular sector, bound by four boundaries. Two boundaries converging at the origin emerge from the restricted θlab coverage of HiRA10, ranging from 30◦ to 75◦ . The other two boundaries are defined by our energy cuts, with the lower and upper boundaries denoting our energy threshold and the particle punch-through energy, respectively. Through the reaction kinematics, we can derive the analytic forms of these boundaries (Appendix B). The energy-cut boundary is given by q E = p2T + m2 · cosh y . (3.4) The θ-cut boundary is given as ! pT θ = arctan p . (3.5) p2T + m2 · sinh y 49 Figure 3.2 Two-dimensional phase space representation of deuterons from a 40Ca + 58 Ni collision at 140 MeV/u. Note that the horizontal axis has been normalized by the beam rapidity in the lab, and the vertical axis has been normalized by the mass number A of the deuteron. The shaded region indicates the mid-rapidity slice, i.e. 0.4 < ŷ (lab) < 0.6. These analytic boundaries enable us to draw precise phase space cuts, thereby correcting for the yield loss beyond the phase space boundaries. Our primary focus is on the fragments that have actively participated in the collisions. In the lab frame, the normalized rapidity, y (lab) ŷ (lab) ≡ (lab) , (3.6) ybeam can provide information about whether the fragments are beam-like (ŷ (lab) ≈ 1) or target-like (ŷ (lab) ≈ 0). Selection of these participant fragments is typically achieved by imposing a mid-rapidity cut [135]. To keep the notation simple, we will henceforth drop the superscript “(lab)” and use ŷ to denote the normalized rapidity in the lab frame. For constructing transverse momentum spectra, we only select fragments from the mid-rapidity region, i.e. 0.4 < ŷ < 0.6, as indicated by the shaded region in Figure 3.2. As the phase space is confined, the resulting transverse momentum spectrum by cutting 50 w/o correction w/ correction 10 2 d 2 M [(MeV/c) −1 ] 10 d(pT /A) dŷ 3 10 4 100 150 200 250 300 350 400 450 500 pT /A [MeV/c] Figure 3.3 One-dimensional pT /A spectra of deuterons from a 40Ca + 58Ni collision at 140 MeV/u. The black points are the raw spectra, whereas the green points are the spectra after the coverage correction. around ŷ = 0.5 ± 0.1 and projecting onto the pT /A axis is incomplete beyond the phase space boundaries. To account for this loss of yield, we apply a coverage correction to the transverse momentum spectrum by dividing the raw spectrum by the coverage fraction within the mid-rapidity slice at each pT /A bin. Figure 3.3 illustrates the raw and corrected transverse momentum spectra from Figure 3.2. Error bars indicate the statistical uncertainties. Note that for very high and very low pT /A, the coverage correction would inflate the uncertainties, as the coverage fraction approaches zero. To improve the reliability of the data points, we have chosen a coarse binning of 20 MeV/c for pT /A and remove data points with fractional errors exceeding 20%. This correction assumes that the rapidity dependence of the yield is negligible within the mid-rapidity slice. To check the validity of this assumption, we have examined the pT /A spectra with different widths of the mid-rapidity slice, ŷ ± δ ŷ, and the results show that the mean values of the spectra are consistent with each other after the coverage correction, up to pT /A = 450 MeV/c, as indicated by the vertical dashed lines. 51 10 2 d 2 M [(MeV/c) −1 ] 10 3 d(pT /A) dŷ δŷ = 0.04 δŷ = 0.16 10 4 δŷ = 0.08 δŷ = 0.20 δŷ = 0.12 δŷ = 0.10 100 150 200 250 300 350 400 450 500 pT /A [MeV/c] Figure 3.4 One-dimensional pT /A spectra of deuterons from a 40Ca + 58Ni collision at 140 MeV/u with varying widths of the mid-rapidity slice, δ ŷ. Therefore, the two-dimensional pT -spectrum of particle (A, Z) is given by d2 M(A,Z) N(A,Z) (i′ , j ′ )   1 1 (p′T /A, ŷ ′ ) ≈ · ′ · , (3.7) d(pT /A) dŷ Nevent ϵ(θlab ) ∆(pT /A) · ∆ŷ ′ ′ where ϵ(θlab ) is the detector at θlab , which can be calculated using Equation 3.5, and N(A,Z) (i′ , j ′ ) is the count of particles of type (A, Z) in the bin indexed by i′ and j ′ , which correspond to p′T /A and ŷ ′ , respectively. The histogram has a uniform bin widths of ∆(pT /A) = 1 MeV/c and ∆ŷ = 0.01. To obtain the one-dimensional transverse momen- tum spectrum at mid-rapidity region, we integrate the two-dimensional spectrum over the mid-rapidity slice, i.e. " Z # dM(A,Z) ′ 1 0.5+δ/2 d2 M(A,Z) (p /A) = lim dŷ (3.8) d(pT /A) T δ→0 δ 0.5−δ/2 d(pT /A) dŷ 5 X 1 N(A,Z) (i′ , j ′ ) ≈ · ′ , (3.9) Nevent η · ϵ(θlab ) ∆(pT /A) |ŷj′ −0.5|≤0.1 where η = η(p′T /A) is the coverage fraction at p′T /A over the mid-rapidity slice, the summation is over all the bins j ′ that fall within the mid-rapidity region, and i′ is the bin index of p′T /A. While the coverage correction is not needed when studying the yield ratios of the same 52 40 Ca+ 58 Ni at 56 MeV/u 48 Ca+ 64 Ni at 56 MeV/u 10 1 10 2 d 2M d(pT /A) dŷ [(MeV/c) ] 10 3 −1 10 4 10 5 10 6 10 7 40 Ca+ 58 Ni at 140 MeV/u 48 Ca+ 64 Ni at 140 MeV/u 10 1 10 2 d 2M d(pT /A) dŷ [(MeV/c) ] 10 3 −1 10 4 10 5 p He3 10 6 d He4 t 10 100 7 200 300 400 500 600 100 200 300 400 500 600 pT /A [MeV/c] pT /A [MeV/c] Figure 3.5 Transverse momentum spectra of light charged particles from the Ca + Ni systems. isotopes from different collision systems because the phase space boundaries are identical, it is essential when comparing the yield ratios of different isotopes, for example, Yd /Yp and Yt /Y 3He . For consistency, we apply the coverage correction to all transverse momentum spectra in this thesis, except when computing yield ratios for the isoscaling analysis in the subsequent section, section 3.2. In this specific case, we opt to leave the spectra uncorrected to avoid inflating the error bars of ratios, as both the numerator and denominator of the ratio would have their own inflated uncertainties from the coverage correction. In Figure 3.5 and Figure 3.6, we present the transverse momentum spectra of all light charged particles, up to 4He, for reactions on the nickel and tin targets, respectively. Each quadrant in the figures corresponds to a different beam configuration, specified by the beam 53 40 Ca+ 112 Sn at 56 MeV/u 48 Ca+ 124 Sn at 56 MeV/u 10 1 10 2 d 2M d(pT /A) dŷ [(MeV/c) ] 10 3 −1 10 4 10 5 10 6 10 7 40 Ca+ 112 Sn at 140 MeV/u 48 Ca+ 124 Sn at 140 MeV/u 10 1 10 2 d 2M d(pT /A) dŷ [(MeV/c) ] 10 3 −1 10 4 10 5 p He3 10 6 d He4 t 10 100 7 200 300 400 500 600 100 200 300 400 500 600 pT /A [MeV/c] pT /A [MeV/c] Figure 3.6 Transverse momentum spectra of light charged particles from the Ca + Sn systems. isotope and energy. Corrections for coverage have been applied to these spectra. The error bars in these figures represent statistical uncertainties. We have ensured the reliability of data by excluding data points with fractional errors exceeding 20% from the plots. An examination of these spectra reveals a common pattern: the yields of all isotopes decline exponentially or at a faster rate as pT /A increases. Spectral shapes and magnitudes from the same quadrant, hence the same beam configuration, in the two figures resemble each other. This suggests that the spectral characteristics are more influenced by the beam, especially the beam energy, rather than the target, whether as light as nickel or as heavy as tin. In the following section, section 3.2, this beam-energy influence will manifest itself in the isoscaling analysis. 54 3.2 Isoscaling from charged particles In this section, we discuss the isoscaling from light charged particles. A pair of collision systems are needed to define an isoscaling ratio, R21 . Conventionally for neutron-rich systems (δ > 0), the system with larger isospin asymmetry is labeled as “2”, whereas the other system is labeled as “1”. In both systems, an isotope fragment (N, Z) can selected to compute its isoscaling ratio, Y2 (N, Z) R21 (N, Z) = , (3.10) Y1 (N, Z) where Y is the isotope yield from the respective reaction as labeled by the subscript. Table 3.2 lists the isospin asymmetries for all the reaction systems used in this section. reaction system δ 40 Ca + 58Ni 0.02041 48 Ca + 64Ni 0.14286 40 Ca + 112Sn 0.07895 48 Ca + 124Sn 0.18605 Table 3.2 Values of isospin asymmetry δ for all the reaction systems used in isoscaling analysis. In multifragmentation processes, the isoscaling ratio, denoted as R21 (N, Z), exhibits an exponential relationship with neutron number N and proton number Z. This relationship, reported in literature [136, 137, 138, 139], can be mathematically represented as: R21 (N, Z) = C exp (αN + βZ) , (3.11) Here, C is a normalization factor, while α and β are empirical parameters inferred by fitting the isoscaling ratios derived from a range of isotopes. This exponential relationship emerges as matter expands from high-density conditions within the nucleus interior to lower densities, with isotopic clusters forming in the process, and it can be justified using the grand canonical ensemble approach [140, 141]. The parameters α and β correspond to the chemical potential differences of neutrons and protons, respectively, in the two reaction system. They are given by 1 (2) µn − µ(1)  α= n (3.12) T 55 and 1 (2) µp − µ(1)  β= p . (3.13) T (i) Here, T represents the equilibrium temperature after the reaction. The variables µn and (i) µp represent the chemical potentials of neutrons and protons for reaction i, respectively. By convention, reaction 2 is the system with larger isospin asymmetry δ, i.e. it is more neutron- rich than reaction 1. Given the chemical potential of neutrons in the more neutron-rich (2) (1) system (reaction 2) is higher, µn > µn , we expect the parameter α to be positive; similarly, we expect β to be negative. These chemical potentials represent the characteristic energy at which nucleon occupancies in phase space transition from a degenerate Fermi gas to a Boltzmann gas. This energy depends not only on the nucleon’s kinetic energy but also its potential energy, including the Coulomb potential. Isoscaling occurs under conditions where both reaction systems reach the same equilibrium temperature [137, 139]. To quantify this temperature, one can use the yields of isotopes d, t, 3 He, and 4He from experimental data. This method is based on the work of S. Albergo et al. [141] and such temperature is expressed as 14.29 TH−He =  . (3.14) Y (d) · Y (4 He)  ln 1.59 · Y (t) · Y (3 He) This equation will be referred to as the Albergo thermometer. 3.2.1 Ca + Ni systems 40 We have chosen Ca + 58Ni (δ ≈ 0.0204) and 48 Ca + 64Ni (δ ≈ 0.143) as the two isoscaling systems for the Ca + Ni systems. Both beam energies, 56 MeV/u and 140 MeV/u, will be discussed. In Figure 3.7, we present the isoscaling fits for the Ca + Ni systems at 56 MeV/u and 48 140 MeV/u. The isoscaling ratio R21 is computed by dividing the yield of Ca + 64Ni system by the yield of 40Ca + 58Ni system for each isotope, indicated by the vertical axis in logarithmic scale. For the 56 MeV/u systems, the yields are sums from 110 MeV/c to 250 MeV/c in pT /A; for the 140 MeV/u systems, the yields are sums from 130 MeV/c to 350 MeV/c in 56 Y( 48 Ca + 64 Ni) at 56 MeV/u Y( 48 Ca + 64 Ni) at 140 MeV/u 1.8 Y( 40 Ca + 58 Ni) 40 1.8 Y( Ca + Ni) 58 α = 0.46 α = 0.31 1.5 β = − 0.44 β = − 0.23 C = 1.08 1.5 C = 1.06 R21 (N, Z) 1.0 R21 (N, Z) 0.8 1 2 1.0 Z= 2 Z= Z= 1 Z= 0 1 2 0.8 0 1 2 N N Figure 3.7 The isoscaling fits for the Ca + Ni systems at 56 MeV/u (left) and 140 MeV/u (right). The points are experimental measurements, whereas the dashed lines are the isoscaling fits. pT /A. It can be seen that the isoscaling fit (dashed lines) is in good agreement with the experimental data (points) for both beam energies. To be more general, the isoscaling ratio can also be computed as a function of pT /A using the one-dimensional transverse momentum spectra in Figure 3.5. For each pT /A bin, the same isoscaling fit is performed, yielding a pair of α and β values, while sharing the same normalization factor C across all pT /A bins. More precisely, we perform least-squares minimization on the objective function n X X m 2 S= rij (3.15) i=1 j=1 by varying C, αi , and βi . Here, i is the bin index of pT /A, j is the index of the isotopes, n is the number of pT /A bins, m is the number of isotopes, and the residual rij is defined as (ij) R21 − C exp (αi Nj + βi Zj ) rij ≡ (ij) , (3.16) δR21 (ij) (ij) where R21 is the isoscaling ratio of particle j at the i-th pT /A bin and δR21 is the respective standard error. 57 Y( 48 Ca + 64 Ni)/Y( 40 Ca + 58 Ni) at 56 MeV/u 2.00 t 1.75 N−Z=1 d 4 He p 3 He 1.50 R21 1.25 N−Z=0 1.00 0.75 N − Z = − 1 0.50 0.5 α ∆µ/T 0.0 0.5 β 12 pT /A [MeV/c] 10 TH − He [MeV] 8 40 Ca + 58 Ni 6 48 Ca + 64 Ni 0 100 200 300 400 500 pT /A [MeV/c] Figure 3.8 Top: The isoscaling ratio spectra for the Ca + Ni systems at 56 MeV/u. The points are experimental measurements, whereas the dashed lines are the isoscal- ing fits. Middle: The fitted α and β values. Bottom: The TH−He temperatures for both systems. The findings are depicted in Figure 3.8 and Figure 3.9, presenting the Ca + Ni systems at 56 MeV/u and 140 MeV/u, respectively. Each figure shares a common horizontal axis, representing the transverse momentum per nucleon, pT /A. The uppermost panels display the isoscaling ratio spectra, where data points represent experimental measurements and dashed lines indicate the isoscaling fits. The term “isoscaling holds” is employed when the 58 Y( 48 Ca + 64 Ni)/Y( 40 Ca + 58 Ni) at 140 MeV/u 2.00 t 1.75 d 4 He p 3 He 1.50 N−Z=1 R21 1.25 N−Z=0 1.00 N−Z= −1 0.75 0.50 0.5 α ∆µ/T 0.0 β 0.5 12 pT /A [MeV/c] 10 TH − He [MeV] 8 40 Ca + 58 Ni 6 48 Ca + 64 Ni 0 100 200 300 400 500 pT /A [MeV/c] Figure 3.9 Top: The isoscaling ratio spectra for the Ca + Ni systems at 140 MeV/u. The points are experimental measurements, whereas the dashed lines are the isoscaling fits. Middle: The fitted α and β values. Bottom: The TH−He temperatures for both systems. dashed lines align with the experimental data points. For some isotopes, the dashed lines might stretch beyond the last data point at higher pT /A values. This extension is a result of the isoscaling fit, which necessitates data points from a minimum of three isotopes at each pT /A bin. In this analysis, the fit is executed as long as there are at least three isotopes with valid R21 (pT /A) values. The panels in the middle exhibit the fitted α and β as functions of 59 pT /A, both of which remain fairly constant throughout the available pT /A range. Although α and β manifest the opposite signs, their magnitudes are comparable. The lowest panels reveal the TH−He temperatures for both systems, which demonstrate a steadily increasing trend with pT /A. As anticipated, when isoscaling is observed, the temperatures are found to be similar for both systems. Upon comparing the isoscaling ratio spectra at 56 MeV/u (Figure 3.8) and 140 MeV/u (Figure 3.9), it is discernible that isoscaling persists over a broader pT /A range because more high-pT isotopes are available at increased beam energies. The magnitudes of α and β are smaller at 140 MeV/u than at 56 MeV/u due to the increased temperatures at higher beam energies. 3.2.2 Ca + Sn systems 40 The chosen pair of isoscaling systems for the Ca + Sn combination in our study are Ca + 112 48 Sn (δ ≈ 0.0790) and Ca + 124Sn (δ ≈ 0.186). For both systems, observations are carried out at two different beam energies, 56 MeV/u and 140 MeV/u. Y( 48 Ca + 124 Sn) at 56 MeV/u Y( 48 Ca + 124 Sn) at 140 MeV/u 40 112 40 112 1.8 Y( Ca + Sn) 1.8 Y( Ca + Sn) α = 0.46 α = 0.34 1.5 β = − 0.44 β = − 0.23 C = 1.05 1.5 C = 1.05 R21 (N, Z) 1.0 R21 (N, Z) 2 0.8 =2 1.0 Z= Z =1 Z 1 Z= 0 1 2 0.8 0 1 2 N N Figure 3.10 The isoscaling fits for the Ca + Sn systems at 56 MeV/u (left) and 140 MeV/u (right). The points are experimental measurements, whereas the dashed lines are the isoscaling fits. 60 In Figure 3.10, we present the isoscaling fits for the Ca + Sn systems at 56 MeV/u and 140 MeV/u. The data points include yields from 110 MeV/c to 250 MeV/c in pT /A for the 56 MeV/u systems and from 130 MeV/c to 350 MeV/c in pT /A for the 140 MeV/u systems. Similarly to the Ca + Ni systems, the isoscaling fit (dashed lines) is in good agreement with the experimental data (points) for both beam energies. Y( 48 Ca + 124 Sn)/Y( 40 Ca + 112 Sn) at 56 MeV/u 2.00 t 1.75 N−Z=1 d 4 He p 3 He 1.50 R21 1.25 1.00 N−Z=0 0.75 N − Z = − 1 0.50 0.5 α ∆µ/T 0.0 0.5 β 12 pT /A [MeV/c] 10 TH − He [MeV] 8 40 Ca + 112 Sn 6 48 Ca + 124 Sn 0 100 200 300 400 500 pT /A [MeV/c] Figure 3.11 Top: The isoscaling ratio spectra for the Ca + Sn systems at 56 MeV/u. The points are experimental measurements, whereas the dashed lines are the isoscaling fits. Middle: The fitted α and β values. Bottom: The TH−He temperatures for both systems. 61 Y( 48 Ca + 124 Sn)/Y( 40 Ca + 112 Sn) at 140 MeV/u 2.00 t 1.75 d 4 He p 3 He 1.50 N−Z=1 R21 1.25 N−Z=0 1.00 N−Z= −1 0.75 0.50 0.5 α ∆µ/T 0.0 β 0.5 12 pT /A [MeV/c] 10 TH − He [MeV] 8 40 Ca + 112 Sn 6 48 Ca + 124 Sn 0 100 200 300 400 500 pT /A [MeV/c] Figure 3.12 Top: The isoscaling ratio spectra for the Ca + Sn systems at 140 MeV/u. The points are experimental measurements, whereas the dashed lines are the isoscaling fits. Middle: The fitted α and β values. Bottom: The TH−He temperatures for both systems. The exact same fitting methodology is employed to this data set, and similar findings are obtained. These are graphically represented in Figure 3.11 and Figure 3.12, illustrating the outcomes for the 56 MeV/u and 140 MeV/u systems, respectively. It is worth noting that the isoscaling phenomena, as well as the behavior of α and β and the TH−He temperatures, show consistent trends with the previous Ca + Ni systems, reinforcing the observations made 62 in the previous section: • Isoscaling persists over a broader pT /A range at 140 MeV/u than at 56 MeV/u. • α and β stay relatively constant throughout the available pT /A range. • α and β are comparable in magnitude but opposite in sign. • Magnitudes of α and β are smaller at 140 MeV/u than at 56 MeV/u. • TH−He is higher in the 140 MeV/u system than in the 56 MeV/u system. • TH−He increases with pT /A. The behavior of α and β mimics that observed in the Ca + Ni system. For both systems, we find that α and β exhibit lower magnitudes at the higher beam energy of 140 MeV/u, compared to 56 MeV/u. This expected behavior arises as α and β are inversely proportional to temperature, so their magnitudes decrease as the beam energy, and thus the temperature, increases. Alongside this, the TH−He temperatures for both the 56 MeV/u and 140 MeV/u systems show a consistent increasing trend with pT /A. As previously discussed in subsection 3.1.2, the transverse momentum spectra of Ca + Ni and Ca + Sn demonstrate similarities for the same beam energy. These similarities are also mirrored in the isoscaling parameters α and β. In Figure 3.13, we present the fitted α and β values for all the isoscaling systems across each pT /A bin. The error bars indicate the standard errors derived from the fitting procedure, and the diagonal dashed line provides a visual reference denoting when |α| = |β|. The green arrows represent the directions of increasing pT /A for respective data points. As pT /A ascends, indicating emissions from earlier stages of the reaction, the magnitudes of α and β become more divergent. This divergence may be linked to the more extreme conditions prevalent during the initial stages of the reaction, including higher temperatures, densities, and net charges, thereby influencing the observed α and β values. In this study, we have shown that isoscaling holds true across the entire range of pT /A measured with ||α| − |β|| < 0.2, but J. W. Lee et al. [139] demonstrated 63 0.1 140 MeV/u 0.2 0.3 Ca + Ni Ca + Sn 56 MeV/u β 0.4 0.5 Ca + Ni Ca + Sn 0.60.1 0.2 0.3 0.4 0.5 0.6 α Figure 3.13 The fitted α and β values for all the isoscaling systems. The diagonal dashed line indicates when |α| = |β|. The green arrows indicate the direction of increasing pT /A. in an experiment with a higher beam energy of 270 MeV/u that the isoscaling property ceases beyond a specific pT /A value, producing ||α| − |β|| > 0.3. The phenomenon where |α| ≈ |β| can also be discerned directly in the isoscaling ratio spectra. If we assume that α = −β, then 3.11 can be rewritten as R21 (N, Z) ≈ Ceα(N −Z) . (3.17) This formulation implies that isotopes sharing the same neutron excess number, N − Z, should exhibit comparable isoscaling ratios. This aggregation of isotopes with identical N − Z values is most prominent in the systems with beam energy at 56 MeV/u, as depicted in 3.8 and 3.11. Within these figures, we can see that the p- 3He and d- 4He isoscaling pairs are grouped closely. The t spectrum stands solitary at present. But if we hypothesize that neutrons follow the same scaling pattern, R21 (n) should be similar to R21 (t). 64 This isoscaling ratio grouping becomes less evident but still remains visible in the systems at 140 MeV/u, as illustrated in 3.9 and 3.12. This decrease in grouping is attributed to the smaller magnitudes of |α| and |β| as a consequence of increased temperature. To conclude, any discernible divergence within the isoscaling ratios of the same N − Z groups could suggest that the magnitudes of α and β are not precisely identical. 3.3 Neutron-proton double ratio As previously discussed in subsection 1.2.1, the neutron to proton yield ratio (n/p) is a crucial metric in the investigation of heavy-ion collisions (HICs). It is employed as an instrument for probing the nuclear equation of state (EOS) and the effective masses of nucleons. This ratio serves as a statistical indicator, evaluating the relative abundance of neutrons and protons within specific energy ranges. An intuitive approach for calculating the single n/p ratio involves constructing the absolute neutron spectrum Yn and the absolute proton spectrum Yp , followed by dividing the former by the latter, SRn/p = Yn /Yp . However, this methodology presents certain challenges. Firstly, obtaining an accurate measure of the neutron spectrum Yn tends to be far more complex than measuring the proton spectrum Yp . The subsequent chapters of this thesis will delve deeper into this issue, focusing on neutron detection, analysis, and spectrum reconstruction to shed light on the specific difficulties inherent to this process. Secondly, neutron and proton detections are generally measured by two distinct detectors, each with their own unique detection mechanisms, calibrations, and efficiencies. Therefore, it often presents a non-trivial task to appropriately normalize the two spectra to the same absolute scale with a high degree of precision. Lastly, a significant challenge arises from the fact that many theoretical computations and models, including ImQMD that we employ in this work, struggle to accurately predict the abundance of isotope fragments emitted from HICs. This leads to difficulties when attempting to directly compare Yn and Yp from both simulations and experimental data. To bypass these impediments, the literature has proposed several adaptations to the 65 straightforward single n/p ratio [90, 89, 87]. The focus of this section will be a detailed discussion that eventually leads to the construction and application of the n/p double ratio. 3.3.1 Coalescence-invariant spectra The coalescence mechanism plays a key role in many transport models, facilitating the formation of light clusters from individual nucleons. Models from the Quantum Molecular Dynamics (QMD) family, including ImQMD, typically simulate the dynamical evolution of nucleons as localized Gaussian wave packets with fixed forms [142, 143]. As a result, these simulations yield raw outputs comprising nucleons, each distinguished by their respective isospin projections (n or p), positions, and momenta. To compare the simulation results with experimental data, we need a procedure to coalesce the nucleons into clusters. This is because many experiments are capable of measuring fragments above A = 1, such as d, t, 3 He, and 4 He, and it is generally desirable to compare the yields from as many isotopes as possible. Hence, a process known as coalescence is needed to coalesce individual nucleons at the end of the simulation into clusters. Although the specifics of the coalescence model may vary, the general method involves grouping nucleons into clusters based on their relative positions and momenta. For instance, nucleons residing close to one another in the six-dimensional (r, p) phase space are grouped, with the exact parameters that defines proximity established empirically. Nucleons grouped within the same cluster are then treated as a bound state, characterized by a specific mass number A and charge number Z, like a nucleus. Despite its practical utility, the coalescence model does not provide an entirely precise representation of light cluster formation. The derived cluster yields often diverge from experimental observations. To reconcile this disparity, the concept of the coalescence-invariant (CI) spectrum was introduced [90, 89]. The CI spectrum effectively decomposes clusters back into constituent neutrons and protons, allocating to each nucleon its respective position and momentum post-coalescence. This strategy minimizes the bias introduced by the coalescence model, facilitating a more meaningful comparison between experimental data and theoretical 66 calculations. Although it may appear counterintuitive to initially perform coalescence and subsequently disassemble the clusters, it is important to note that it is typically not feasible to measure isotopes with exceedingly high mass numbers A with high statistical accuracy in experi- ments. The yield of a particular isotope generally decreases as A increases. Moreover, other experimental constraints, such as detector acceptance and resolution, inhibit the reliable mea- surement of isotopes with high A. Hence, if we were to exclude the coalescence step and immediately compile all the neutrons and protons to construct the spectra, we would inad- vertently include a large quantity of nucleons originating from clusters with high A values based on the coalescence model, even though they are not experimentally accessible. In this section, we will establish the CI spectra for neutrons and protons, utilizing light charged particles such as p, d, t, 3He, and 4He. Specifically, the CI proton spectrum can be formulated as follows: X Y (p(CI) ) = Z · Y (A, Z) (3.18) A≤4, Z≤2 = Y (p) + Y (d) + Y (t) + 2Y ( 3He) + 2Y ( 4He) . (3.19) Analogously, the CI neutron spectrum can be expressed as: X Y (n(CI) ) = (A − Z) · Y (A, Z) (3.20) A≤4, Z≤2 = Y (n) + Y (d) + 2Y (t) + Y ( 3He) + 2Y ( 4He) . (3.21) However, this CI neutron spectrum still necessitates the absolute yield of the neutron, Y (n), which HiRA10 cannot measure. The concept of pseudo-neutron, ñ, which is about crafting the neutron spectrum from only charged-particle spectra, was introduced in an attempt to circumvent this issue [144]. In Figure 3.14, the pseudo-neutron spectra for the two Ca + Ni systems at 140 MeV/u are plotted alongside the proton spectra for comparison. The idea hinges on the assumption that neutrons, like charged particles, follow the isoscaling behavior discussed in section 3.2. Hence, 67 40 Ca + 58 Ni at 140 MeV/u 48 Ca + 64 Ni at 140 MeV/u 10 1 d 2 M [(MeV/c) −1 ] 10 2 d(pT /A) dŷ p pseudo-n 10 150 3 200 250 300 350 400 150 200 250 300 350 400 pT /A [MeV/c] pT /A [MeV/c] Figure 3.14 The pseudo-neutron spectra from the two Ca + Ni systems at 140 MeV/u. the isoscaling ratio R21 (N = 1, Z = 0) for neutrons should align with the simple exponential form R21 = C exp (αN + βZ), a model which has successfully described the isoscaling ratios of light charged particles. Under this premise, it can be demonstrated that: R21 (1, 0) R21 (2, 1) R21 (t) = ∴ R21 (n) = · R21 (p) . (3.22) R21 (0, 1) R21 (1, 2) R21 ( 3He) Consequently, the yield of pseudo-neutrons is approximated by [144, 145]: Y (t) Y (ñt/ 3He ) ≡ · Y (p) . (3.23) Y ( 3He) While this approach is appealing, the effects of the Coulomb interaction and secondary decay processes need to be considered. Particularly at lower energies, the Coulomb force continues to act on charged particles, altering the observed clusters depending on their charge-to-mass ratio. However, this effect is partially mitigated in R21 as both reactions 1 and 2 possess the same total charge. Further, the comparison of single and double ratios with transport model predictions implicitly assumes the neglect of secondary decay. While this may be approximately correct for the less favored decay channels, such as tritons and helium-3 due to their lower decay q-values, the dominant decay channels for the unstable emitted nuclei—protons, neutrons, and alpha particles—might introduce notable contributions. 68 Despite these considerations, the pseudo-neutron approximation provides valuable insights especially when the neutron spectrum is a lot more challenging to measure than the charged particle spectra. After this chapter, we will discuss our effort to measure the neutron spectrum directly using the vLANA neutron detection system. 1.4 1.2 1.0 SR(n (CI) /p (CI) ) 0.8 0.6 0.4 48 Ca + 64 Ni 0.2 40 Ca + 58 Ni 0.0100 150 200 250 300 350 400 pT /A [MeV/c] Figure 3.15 The coalescence-invariant single n/p ratio as a function of pT /A for the Ca + Ni systems at 140 MeV/u. Having pseudo-neutron enables the construction of the CI neutron spectrum utilizing exclusively the data from charged-particle spectra. In Figure 3.15, we showcase the CI single n/p ratio as a function of pT /A for the Ca + Ni systems at 140 MeV/u. These ratios are directly derived from the transverse momentum spectra in Figure 3.5, with the error bars representing the statistical uncertainties: Y (n(CI) ) SR(n(CI) /p(CI) ) = . (3.24) Y (p(CI) ) 40 Both reaction systems, Ca + 58Ni and 48 Ca + 64Ni, are neutron-rich (δ > 0), with the latter exhibiting a higher isospin asymmetry than the former. This disparity is reflected in the fact that the single ratio values for the latter system are persistently higher across the entire 40 pT /A range. We label Ca + 58Ni as system 1 and 48 Ca + 64Ni as system 2. To further suppress systematic biases instigated by the distinct construction procedures 69 1.8 1.7 1.6 1.5 DR(n (CI) /p (CI) ) 1.4 1.3 1.2 1.1 1.0100 150 200 250 300 350 400 pT /A [MeV/c] Figure 3.16 The coalescence-invariant double n/p ratio as a function of pT /A for the Ca + Ni systems at 140 MeV/u as demonstrated in Figure 3.15. of CI n and CI p, we introduce the CI double n/p ratio, which is calculated by the ratio of single CI n/p ratios, SR2 (n(CI) /p(CI) ) DR(n(CI) /p(CI) ) = . (3.25) SR1 (n(CI) /p(CI) ) The outcome is presented in Figure 3.16. It is evident that the CI double n/p ratio persistently surpasses unity, increasing across the entire pT /A range. We will now adopt this experimental double ratio from HiRA10 against the theoretical calculations given by ImQMD. 3.4 Bayesian analysis of double ratio spectrum In this section, we will discuss the use of Bayesian analysis on the double ratio spectrum to infer equation-of-state (EOS) parameters and effective masses of nucleons. We constrain ourselves to vary S0 ≡ S(ρ0 ) (3.26) and ∂S(ρ) L ≡ 3ρ0 , (3.27) ∂ρ ρ=ρ0 where S(ρ) is the symmetry energy and ρ0 is the nuclear saturation density at ≈ 0.16 fm−3 . As for the effective masses, we will vary the isoscalar effective mass m∗s and the isovector 70 effective mass m∗v , seeking to constrain these parameters by comparing the experimental double ratio spectrum presented in Figure 3.16 with the theoretical calculations from ImQMD. Typically, least-squares regression is employed to find the best-fit parameters within the input space of (S0 , L, m∗s , m∗v ) for a given observable, like the double ratio spectrum. However, least-squares regression faces several challenges. Although it can provide uncertainty estimates, these measures often assume Gaussian residuals. The lack of comprehensive uncertainty quantification limits our understanding of the inferred EOS parameters and effective masses. Furthermore, least-squares regression lacks a mechanism for systematically incorporating prior knowledge about parameters, which could be particularly beneficial given the existing uncertainties around these parameters. In contrast, Bayesian analysis provides a more comprehensive approach to uncertainty quantification through full posterior distributions and enables the incorporation of prior knowledge systematically. Thanks to advancements in computational hardware and software, Bayesian analysis has become increasingly popular in recent years. In addition, Bayesian analysis provides a measure of protection against overfitting because any tendencies toward overfitting are reflected in increased parameter uncertainty in the posterior distributions. Given these advantages, we opt for Bayesian analysis for our task. The remaining of this section will start with a brief introduction to Gaussian Process regression and Bayesian inference, followed by a discussion on the inference results from the double ratio spectrum. 3.4.1 Gaussian Process emulation Gaussian Process (GP) emulation [146], also known as GP regression, is a non-parametric regression technique applied across numerous domains, including nuclear physics [147, 148, 149]. It offers several advantages, such as flexibility and inherent uncertainty quantification, making it well-suited for Bayesian analysis. Of paramount importance to our work is the GP emulator’s ability to expedite the inference of the double ratio prediction by ImQMD across any point x(ImQMD) = (S0 , L, m∗s , m∗v ) ∈ R4 . (3.28) 71 This acceleration is significant, especially when considering that while simulating a single collision event with ImQMD takes approximately ∼ 1 s, and obtaining an accurate prediction of the double ratio spectrum necessitates at least ∼ 105 collision events. Conducting regression, be it least-squares or Bayesian, entails evaluating the double ratios at numerous points (∼ 103 ) within the parameter space. Given the computational expense, directly employing ImQMD for this task is unfeasible. In contrast, a well-trained GP emulator can predict the double ratio spectrum across any point in the parameter space swiftly (≲ 1 s) while simultaneously quantifying the uncertainty. For a more general language, we represent the input parameters as x ∈ Rn , in which n is the dimension of the input space, and the output space as y ∈ R. A Gaussian Process constitutes a collection of random variables, where any finite set possesses a joint Gaussian distribution. In the context of regression, it is a prior over functions, which is fully characterized by its mean function µ(x) and covariance k(x, x′ ): µ(x) = E[f (x)] , (3.29) k(x, x′ ) = E[(f (x) − µ(x))(f (x′ ) − µ(x′ ))] . (3.30) Here, f (x) is the function we aim to learn. Typically, we set µ(x) = 0 for simplicity as the GP is flexible enough to model the mean arbitrarily well. The covariance k(x, x′ ), also known as the kernel, is a similarity measurement between function values at two different points in the input space. In this work, we will always be using the radial basis function (RBF) kernel, also known as the squared exponential kernel: n ! 1 X (xi − x′i )2 k(x, x′ ) = exp − , (3.31) 2 i=1 ℓ2i where ℓi > 0 is the length scale of the i-th dimension that controls the smoothness of the function. In the training process of the GP emulator, we provide a set of training data, which consists of input points X ∈ Rm×n and output values y ∈ Rm . This data set is used to learn the hyperparameters of the GP, specifically the length scales ℓi . This is typically done 72 by maximizing the log marginal likelihood of the training data under the GP prior, which provides a balance between fit to the data and complexity of the function. Predicting at a new point x∗ ∈ Rn involves having a joint Gaussian distribution over the observed output values y ∈ Rm and the predicted function value f∗ = f (x∗ ) ∈ R:      2 y  K(X, X) + diag(σ ) K(X, x∗ )   ∼ N 0,   , (3.32) f∗ K(x∗ , X) k(x∗ , x∗ ) where K(X, X′ ) ∈ Rm×m has matrix elements k(x, x′ ), K(X, x∗ ) = K(x∗ , X)T ∈ Rm denotes the vector with elements k(x, x∗ ), and σ 2 ∈ Rm is the vector of noise variances in the training data. Conditioning this joint Gaussian on the observed values y ∈ Rm , we get the posterior distribution over f∗ which is also Gaussian: f∗ ∼ N (µ∗ , σ∗2 ) , (3.33) where µ∗ = K(x∗ , X)[K(X, X) + diag(σ 2 )]−1 y , (3.34) σ∗2 = K(x∗ , x∗ ) − K(x∗ , X)[K(X, X) + diag(σ 2 )]−1 K(X, x∗ ) . (3.35) parameter range S0 [25, 45] MeV L [20, 170] MeV m∗s [0.6, 1.0] mN m∗v [0.6, 1.2] mN Table 3.3 Ranges of the parameters used in the Latin hypercube sampling. We have trained the GP emulator to predict the CI n/p ratios for a given input point of (pT /A, x(ImQMD) ). The training data is generated by running the simulation at 75 different x(ImQMD) points, randomly chosen according to the Latin hypercube sampling (LHS) [150] over the ranges listed in Table 3.3. For each x(ImQMD) point, we simulate 2 × 106 collision events. The simulated events are then analyzed in the same way as the experimental data, and the CI n/p ratios are calculated. 73 40 Ca + 58 Ni at 140 MeV/u 48 Ca + 64 Ni at 140 MeV/u 1.6 1.4 1.2 SR(n (CI) /p (CI) ) 1.0 0.8 0.6 0.4200 250 300 350 400 200 250 300 350 400 pT /A [MeV/c] pT /A [MeV/c] Figure 3.17 Coalescence-invariant (CI) n/p single ratios. The data points are the experimental data, replotted from Figure 3.15. The gray lines behind the data points consist of 75 single ratio spectra from ImQMD used to train the GP emulator. In Figure 3.17, we show the CI n/p single ratios, where the data points are the experimental data from Figure 3.15 and the gray lines represent the ImQMD single ratios used to train the 40 GP emulator. For the Ca + 58Ni system, the simulation fails to reproduce the experimental data when pT /A ≥ 350 MeV/c over the entire parameter space specified in Table 3.3. For 48 the Ca + 64Ni system, the experimental data is placed near the edge of the envelope of the simulated spectra. These observations suggest that the ImQMD model is not ideal for predicting the CI single ratios for these two systems. While GP emulator is capable of mathematically extrapolating its predictions beyond the training parameter space, there is no guarantee that the extrapolation is accurate or physically meaningful. Hence, we will only 40 use the GP emulator to predict the CI n/p double ratios for the Ca + 58Ni system, in which most of the systematics present in the model and the analysis procedure are canceled out. The double ratios are presented in Figure 3.18. 74 1.8 1.7 1.6 1.5 DR(n (CI) /p (CI) ) 1.4 1.3 1.2 1.1 1.0200 225 250 275 300 325 350 375 400 pT /A [MeV/c] Figure 3.18 Coalescence-invariant (CI) n/p double ratios. The data points are the experimental data, replotted from Figure 3.16. The gray lines behind the data points consist of 75 double ratio spectra used to train the GP emulator. 3.4.2 Bayesian inference In the context of Bayesian inference, the process of estimating the posterior distribution is given by Bayes’ theorem: P(y | x)P(x) P(x | y) = , (3.36) P(y) where P is probability, x is the parameters of interest (e.g., x(ImQMD) ), and y is the observed data (e.g., the CI double ratio spectrum arranged as a vector). In this equation: • P(x | y) is the posterior, which represents our updated belief about the parameters x after observing the data y. This is the result that we are interested in. • P(y | x) is the likelihood of observing the data y given the parameters x. It measures the compatibility of the observed data with the chosen parameters. This is given by the GP emulator, which was first trained on the simulation data from ImQMD. • P(x) is the prior of the parameters, encapsulating our prior knowledge or beliefs about the parameters before observing the data. This is often the most controversial 75 component of Bayesian inference as the choice of prior is subjective and can have a significant impact on the posterior. • P(y) is the evidence or marginal likelihood. It acts as a normalization constant ensuring that the posterior distribution is a valid probability distribution. However, it is often computationally challenging to calculate directly due to the necessity of integrating over all possible parameter of interest x. The computational challenge of calculating the evidence P (y) had historically limited the application of Bayesian inference. However, the development of Markov Chain Monte Carlo (MCMC) methods and the advancements of modern computers have made Bayesian inference more accessible. MCMC is a suite of algorithms designed to sample from complex, high-dimensional probability distributions, particularly when direct sampling is impractical or impossible. These methods are typically used to sample from the posterior distribution, yielding a representative sample of parameter values for making inferences and predictions. In our current work, we have implemented both the GP emulator and Bayesian inference using the Python package PyMC of version 5.5.0 [151]. Within PyMC, we employ a imple- mentation of the MCMC method known as the No-U-Turn Sampler (NUTS), an extension to the Hamiltonian Monte Carlo algorithm that eliminates the need to manually tune the proposal distribution. Using the CI double ratio spectrum from the Ca + Ni systems at 140 MeV/u, we perform Bayesian inference on the parameters (S0 , L, m∗s , m∗v ). The posterior distributions of these parameters are presented in Figure 3.19. The diagonal plots show the marginal distributions of each parameter, i.e. one-dimensional projections of the posterior distributions, whereas the off-diagonal plots show the pairwise correlations between parameters, i.e. two-dimensional projections of the posterior distributions. Due to the limited number of observables, the posterior distributions are broad, indicating that using CI double ratio solely to constrain the parameters results in large uncertainties. 76 S0 170 140 L [MeV] 110 80 50 L 20 1.0 0.9 ms∗ /mN 0.8 0.7 ms∗ 0.6 1.2 1.0 mv∗ /mN 0.8 mv∗ 0.625 30 35 40 45 20 50 80 110 140 170 0.6 0.7 0.8 0.9 1.0 0.6 0.8 1.0 1.2 S0 [MeV] L [MeV] ms∗ /mN mv∗ /mN Figure 3.19 Posterior distributions of (S0 , L, m∗s , m∗v ) for the Ca + Ni systems at 140 MeV/u. The diagonal plots show the marginal distributions of each parameter, whereas the off-diagonal plots show the pairwise correlations between parameters. The units of S0 and L are MeV, whereas the units of m∗s and m∗v are mN . The red dashed line in the m∗v -m∗s correlation plot indicates the line where m∗v = m∗s . Despite the inherent challenges, our analysis does unveil an intriguing correlation between m∗s and m∗v , as illustrated in the m∗v -m∗s correlation plot in Figure 3.19. A red line is drawn in this plot to indicate the line where m∗v = m∗s . The posterior of (m∗s , m∗v ) is narrowly concentrated above this line, suggesting a preference for mN mN fI ≡ ∗ − ∗ (3.37) ms mv 77 0.14 ± 0.11 −0.17 ± 0.12 0.4 0.2 0.0 0.2 0.4 0.4 0.2 0.0 0.2 0.4 fI ∆mnp∗ /δ Figure 3.20 Left: Posterior distributions of fI . Right: Posterior distributions of ∆m∗np /δ. to be positive, or equivalently, the effective mass splitting, ∆m∗ ≡ m∗n − m∗p , (3.38) is favored to be negative. The posteriors for these two quantities are presented in Figure 3.20, showing that fI = 0.14 ± 0.11 or ∆m∗np /δ = −0.17 ± 0.12. To check the consistency of this result, we perform the same Bayesian inference on the simulation data from ImQMD, i.e. a closure test. Specifically, this involves the application of Bayesian inference to synthetic data generated from the model, a valuable test to assess the efficacy of our approach. Our procedure has produced a coalescence-invariant double ratio spectrum for 75 distinct parameter points x(ImQMD) . This rich set of data allows us to perform leave-one-out (LOO) closure tests. First, we select one parameter point to serve as synthetic experimental data. Next, we train a GP emulator using the remaining 74 parameter points and subsequently perform Bayesian inference on the selected synthetic data. This procedure is repeated for all 75 parameter points, generating us 75 posterior distributions. In Figure 3.21, we display the posterior distributions of u(ImQMD) for one arbitrarily chosen synthetic coalescence-invariant double ratio spectrum. The stars and dash-dotted vertical lines represent the actual parameter values for (S0 , L, m∗s , m∗v ). Similar to the posteriors from the experimental CI double ratio shown earlier in Figure 3.19, the posteriors derived from the synthetic data exhibit substantial spreads for all parameters. Nonetheless, the strong 78 S0 170 140 L [MeV] 110 80 50 L 20 1.0 0.9 ms∗ /mN 0.8 0.7 ms∗ 0.6 1.2 1.0 mv∗ /mN 0.8 mv∗ 0.625 30 35 40 45 20 50 80 110 140 170 0.6 0.7 0.8 0.9 1.0 0.6 0.8 1.0 1.2 S0 [MeV] L [MeV] ms∗ /mN mv∗ /mN Figure 3.21 Posterior distributions of (S0 , L, m∗s , m∗v ) for a synthetic CI double n/p ratio spectrum generated from ImQMD. The star markers in the off-diagonal panels and the dash-dotted vertical lines in the diagonal panels indicate the true parameter values. The red dashed line in the m∗v -m∗s correlation plot indicates where m∗v = m∗s . correlation between m∗s and m∗v can still be observed. Similar patterns are observed for all other 74 LOO closure tests. In Figure 3.22, we present the correlation between the true and predicted values for S0 , L, m∗s , and m∗v . Each point represents the posterior mean of the respective parameter from one of the 75 LOO closure tests, with the error bars indicating the posterior standard deviations. The correlation coefficients between the true and predicted values are 0.73+0.04 +0.03 +0.08 −0.04 , 0.79−0.04 , 0.37−0.09 , and 79 45 .04 170 .03 corr = 0.73 +0 −0.04 corr = 0.79 +0 −0.04 140 40 Predicted S0 [MeV] Predicted L [MeV] 110 35 80 30 50 2525 30 35 40 45 2020 50 80 110 140 170 True S0 [MeV] True L [MeV] 1.0 .08 1.2 .05 corr = 0.37 +0 −0.09 corr = 0.62 +0 −0.05 1.1 0.9 1.0 Predicted ms∗ /mN Predicted mv∗ /mN 0.8 0.9 0.8 0.7 0.7 0.60.6 0.7 0.8 0.9 1.0 0.60.6 0.7 0.8 0.9 1.0 1.1 1.2 True ms∗ /mN True mv∗ /mN Figure 3.22 The correlation between true and predicted values for S0 , L, m∗s , and m∗v from the synthetic data. The black dashed diagonal line indicates the line where the two quantities are equal. Each blue data point represents the result from one of the 75 LOO closure tests. 0.62+0.05 ∗ ∗ −0.05 for S0 , L, ms , and mv , respectively. The asymmetric error bars show the 95% confidence intervals. From the figure, it can be concluded that using only the CI double ratio to constrain S0 , L, m∗s , and m∗v is not particularly effective. If we focus on the correlation between m∗s and m∗v , or equivalently, the quantity fI , a more accurate prediction can be established. In Figure 3.23, we present the correlation between the true and predicted fI values from the synthetic data. Each data point represents the result 80 0.8 .02 corr = 0.86 +0 −0.02 0.4 Predicted fI 0.0 0.4 0.80.8 0.4 0.0 0.4 0.8 True fI Figure 3.23 The correlation between true fI and predicted fI from the synthetic data. The black dashed diagonal line indicates the line where the two quantities are equal. Each blue data point represents the result from one of the 75 LOO closure tests. from one of the 75 LOO closure tests. It can be seen that the data points are concentrated around the diagonal line where the true and predicted values are equal. The correlation coefficient for fI is calculated to be 0.86+0.02 −0.02 , showing not only a higher correlation, but also a smaller spread than of S0 , L, m∗s , and m∗v . This completes the closure test, demonstrating that the inference of fI from the CI double ratio is capable of reproducing the true fI value. If we consider only the signs of the predictions, the predicted fI mean values share the same sign as the true fI values in 71 out of 75 LOO closure tests, giving a 94.7% accuracy. Nonetheless, it is far from clear to understand why the strong correlation between m∗s and m∗v exists. In future investigation, it would be informative to understand if similar correlations between m∗s and m∗v exist in closure tests using different transport models. Each model, with its unique assumptions, approximations, and momentum dependence of the symmetry potential, can greatly affect the results. Using different models can help to assess the robustness of the observed correlation and provide further insights into its underlying physics. 81 Starting from the next chapter, we will shift our focus from charged particles to neutrons. In Chapter 5, we will revisit the Bayesian analysis of the CI double ratio spectrum by replacing the pseudo-neutron in the CI neutron with the actual neutron spectrum from the vLANA neutron detection system. 82 CHAPTER 4 CALIBRATIONS FOR THE NEUTRON DETECTION SYSTEM In this chapter, we discuss the calibration procedures for the neutron detection system, vLANA (LANA with Veto wall). These procedures are performed by treating vLANA as a single detection system, independent of the other two detection systems, the HiRA10 and Microball. The sections are arranged according to the order of the data analysis pipeline. Procedure described in each section may depend on the results from all previous sections. The vLANA is composed of the Forward array (FA), the Veto wall (VW), the Neutron wall A (NWA), and the Neutron wall B (NWB). The two Neutron walls, NWA and NWB, are collectively called the Large Area Neutron Array (LANA) [98, 152, 123]. In this dissertation work, only data from the NWB is used. While adding NWA to the analysis would enhance the overall detection efficiency after accounting for multiple scattering between the two walls, it is not essential for this work as the angular coverage of NWA is overlapping with that of the NWB, hence not providing any additional kinetic phase space for the measurement of neutron spectra. Calibrate procedures are performed for every NWB bar because the bar-to-bar calibration parameters are not expected to be the same. As long as there is sufficient statistics, we also repeat the same procedures for every single run. This allows us to pick up any systematic shifts in the calibration parameters in between runs, which may be caused by changes in the experimental conditions. Nonetheless, due to the stochastic nature of collision events, the calibration parameters would still be fluctuating from run to run despite identical experimental conditions. To distinguish between statistical fluctuations and systematic shifts in the calibration parameters programmatically, we adopt a change point detection algorithm that uses a linearly penalized segmentation (Pelt) [153] implemented in the ruptures Python package [154, 155]. In Figure 4.2, we show a demonstration of the change point detection algorithm, in which p0 and p1 are two synthetic calibration parameters that change as a function of time. The 83 vLANA (neutron detection system) Neutron Neutron wall A Microball wall B (impact parameter Veto detection system) wall LANA 𝜃lab = 39.37∘ beam line Forward array target HiRA10 (charged-particle detection system) Figure 4.1 A top view of the experimental setup, including the vLANA boxed by the blue dashed line. Not to scale. 1.2 p0 p1 Parameter value [arb. units] 1.0 0.8 0.6 Chunk 1 Chunk 2 Chunk 3 0.4 0.2 0.0 0 20 40 60 80 100 Run number Figure 4.2 The change point detection algorithm gives two breakpoints marked by the vertical dashed lines, splitting the data into three chunks labeled as 1, 2, and 3. The horizontal lines represent the average values of p0 and p1 in each chunk. 84 algorithm automatically detects two breakpoints right before run 20 and run 70, marked by the vertical dashed lines. These breakpoints identified by the algorithm are also apparent to a trained eye, indicating a change in the experimental conditions that affects the calibration parameters. These two breakpoints split the data into three “chunks” of runs. The average values of p0 and p1 in each chunk, as shown by the horizontal lines, are then used as the final calibration parameters for the corresponding runs to minimize the statistical fluctuations. 4.1 Hit position calibration When an ionizing radiation hits a neutron wall bar, scintillation occurs in the NE-213 liquid. Light is then produced from the point of interaction and propagates along the NW bar toward both ends until it is collected by the photomultiplier tubes (PMTs). Hit position calibration compares the arrival times of the light pulses at the two PMTs to determine the position of the hit on the NW bar. In the following subsections, we first discuss the coordinate systems used in the analysis of neutron detection system (vLANA). This provides a common basis for any discussion that involves the position of a hit on the NW bar, and it will be used throughout the rest of this dissertation. Following that, we delve into the position calibration procedure, which is performed for every NWB bar and every run. 4.1.1 Coordinate systems To specify the position of a hit on the NW bar, it is necessary to establish a coordinate system. There are two common coordinate systems used in the analysis of neutron detection system (vLANA) — the lab coordinate system and the bar coordinate system. The lab coordinate system sets its origin at the beam-target interaction point. The z-axis is along the beam direction, whereas the y-axis is pointing upward (opposite to gravity). The x-axis is then given by the right-hand rule, which is pointing to the left when looking along the beam direction. When extracting the physics information from the data, it is more convenient to express the lab coordinate system using the spherical coordinates, which are denoted by θ (polar angle), ϕ (azimuthal angle), and r (radial distance). The entire vLANA 85 sits in the northern hemisphere (θ < 90◦ ), with VW and LANA further confined in the −90◦ < ϕ < 90◦ range (+x hemisphere). The bar coordinate system is defined for each NW bar and assumes the NW bar can be modeled as a simple cuboid. Its origin is defined at the geometric center of the NW bar. The x-axis is along the bar’s longest dimension, whereas the y-axis is pointing upward. The z-axis is then given by the right-hand rule, which is roughly pointing toward the target. Figure 4.3 shows a schematic view of the NW bar coordinate system. In this dissertation, we also use the terms “left” and “right” to refer to the x-axis direction. Left is defined as xbar < 0 and right is defined as xbar > 0. The conversion between the two coordinate systems is given by rlab = P rbar + r0 , (4.1) where rlab = (xlab , ylab , zlab )T , rbar = (xbar , ybar , zbar )T , r0 is the NW bar center in the lab Figure 4.3 A schematic view of the NW bar coordinate system. The drawing is not to scale, however the relative sizes of the NW bar dimensions have been preserved. Excluding the Pyrex thickness, each bar has a dimension of 193.04 cm × 7.62 cm × 6.35 cm (x, y, z). 86 coordinate system, and P is the transformation matrix given by   − cos θNW 0 sin θNW    P =  0 1 0  .  (4.2)   − sin θNW 0 − cos θNW Here, θNW ≈ 39.55◦ . This value is slightly different from the 39.37◦ in Figure 4.1 because the NW bars are not oriented perfectly during the experiment to have their front surfaces facing the target normally. This small difference has been taken in account when running various simulations for the analysis. Nevertheless, the effect of this difference is negligible for the physics objectives of this dissertation. 4.1.2 Calibrating hit positions using Veto wall bar shadows Let t0 = 0 be the time when the beam isotope hits the target and tx be the time when the neutron coming from the beam-target nuclear reaction first hits the neutron wall and induces light. Denote the speed of light propagating along the bar as vx . The arrival times at the left and right PMTs are, respectively, ℓ/2 + x tL = tx + (4.3) vx and ℓ/2 − x tR = tx + , (4.4) vx where ℓ = 193.04 cm is the full length of the NW bar and x is the hit position along the bar longest dimension. Given the design of the NW, there is no reliable way to reconstruct the y and z coordinates of the hit positions accurately. In the experimental data, time information may be delayed by some constant amount of time due to cable delays and electronic delays. Denote the delays as τL and τR , respectively. This gives the observed arrival times at the left and right PMTs as, respectively, (obs) ℓ/2 + x tL = tx + + τL (4.5) vx and (obs) ℓ/2 − x tR = tx + + τR . (4.6) vx 87 The position x can then be obtained by canceling out tx , resulting in 1 x = vx · ∆t(obs) + x0 , (4.7) 2 where (obs) (obs) ∆t(obs) ≡ tL − tR (4.8) is the observed arrival time difference and 1 x0 = − vx · (τL − τR ) (4.9) 2 is the calibration offset. To calibrate the position, we want to determine the values of vx and x0 for each NW bar. It is achieved by fitting Equation 4.7 to a given set of (∆t(obs) , x) pairs as calibration points. The calibration points are obtained by using the VW bars to cast “shadows” onto the NW bars. VW is placed right in front of NWB with a separation of 48.3 cm, and it completely blocks the view of the NWB from the target with its overlapping width for every adjacent pair of VW bars. While the VW was originally designed to veto charged particles, we have taken advantage of the fact that the VW bars are installed vertically and the NWB bars are placed horizontally, which gives rise to a grid-like pattern that resembles the design of a double-sided silicon strip detector (DSSD). This allows the VW bar shadows casting onto the NWB to be used for position calibration, providing 18 calibration points for each NW bar. During the experiment, we used a laser tracker [156] to measure the physical locations of both NW and VW. The measured coordinates were then used to refine the original 3D models created in Autodesk Inventor [157] before the experiment, resulting in more accurate models that better reflect the actual setup. Using the refined 3D models, we calculated the expected positions of the VW bar shadows on the NWB bars by performing ray tracing from the target. These calculated positions of the VW bar shadows on the NWB serve as the calibration points for x with high fidelity. Figure 4.4 shows the expected positions of the VW bar shadow centers on the NWB bars. 88 25 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 20 NWB bar number 15 10 5 0 100 50 0 50 100 NWB x position [cm] Figure 4.4 The horizontal rectangles indicate the NWB bars, whereas the vertical lines indicate the expected positions of VW bar shadow centers. VW bar numbers are labeled on the top. Shadows of VW bars 1, 2, and 23 fall outside of the NWB; shadows of VW bars 3 and 22 only cast onto the NWB partially, hence they are discarded for position calibration. Only VW bars 4 to 21, inclusively, are used for position calibration, giving us 18 calibration points per NWB bar. In Figure 4.5, we present a visualization of the VW bar shadows for a single experimental run. The VW bars are designed to overlap with neighboring bars to eliminate inter-bar gaps that may allow charged particles to pass through without triggering the VW. Consequently, the VW bar shadows have to be plotted every two bars. The left panel displays the shadows from the even VW bars, while the right panel displays the shadows from the odd VW bars. In the data, each VW shadow cast onto the NWB follows a normal distribution. To extract the mean and standard deviation of this distribution, we fit a Gaussian function to the shadow. Position calibration is performed for every NWB bar in each experimental run using the VW shadows. The calibrated vx values range from 14.5 cm/ns to 15.5 cm/ns (about half the 89 Figure 4.5 Vertical shadows from the veto wall bars casted onto the neutron wall B. The left panel shows shadows from the even veto wall bars, whereas the right panel shows shadows from the odd veto wall bars. The veto wall bar numbers are labeled on the top, with the respective vertical dashed lines indicating the expected shadow centers. speed of light in vacuum) across different bars. The calibrated ∆τ values vary from −80 cm to +90 cm across different bars. It should be noted that the reported vx is slower than the speed of light in bulk NE-213 volume, vNE213 , which is approximately 19.9 cm/ns for light at a wavelength of 425 nm. The reason for this is because light propagating within the NW bar do not generally travel parallel to the x direction. Instead, the light travels at an angle with respect to the x direction, and may reach the PMT after bouncing off the bar’s inner surface multiple times. This results in a longer path length and a slower effective speed for the light to cover the same displacement along the x direction. Table 4.1 provides the calibrated parameters for run 4300 as an example. NWB bar 1 2 3 4 5 6 7 8 9 10 11 12 vx [cm/ns] 15.1 15.3 14.9 15.3 15.5 15.0 14.9 15.0 15.0 15.5 15.0 15.3 x0 [cm] 15.3 -44.2 -3.1 -5.6 21.3 -37.7 -1.2 71.6 -14.8 84.0 35.7 90.8 NWB bar 13 14 15 16 17 18 19 20 21 22 23 24 vx [cm/ns] 15.3 14.7 14.5 15.2 15.1 15.0 14.6 14.9 15.1 15.0 15.0 15.4 x0 [cm] -78.5 -77.2 88.9 0.5 29.5 58.8 -15.5 49.6 19.0 3.5 9.9 -5.9 Table 4.1 Position calibration parameters for all NWB bars in run 4300. 90 The position resolution, as measured using cosmic muons, is approximately 7 cm [122]. This value is comparable to both the transverse height (7.62 cm) and the longitudinal thickness (6.35 cm) of a NWB bar. However, hit positions in the y and z directions within the NWB bar cannot be accurately reconstructed. As a result, we set both the y and z coordinates to zero, i.e., the center of the transverse cross-sectional area of the NWB bar. Nevertheless, this discrete reconstruction of hit positions may introduce spikes in the distributions of other quantities that depend on the hit positions, causing difficulties in the analysis. To address this issue, we randomize the y and z coordinates within the transverse cross-sectional area of the NWB bar for such use cases. The associated smearing due to this randomization is then treated as an uncertainty. The same procedure is repeated for all NWB bars in all runs. The final calibration parameters are obtained by averaging the calibrated parameters from runs within the same chunk of runs, as suggested by the change point detection algorithm. 4.2 Time-of-flight calibration To determine the energy of neutrons, we rely on time-of-flight (ToF) measurements. This is different from how the energies of charged particles are measured in HiRA10 because neutrons do not interact with scintillators in the same way that charged particles do. Charged particles lose energy as they pass through a scintillator, and the resulting light output monotonically depends on the energy loss. In contrast, neutrons can interact with the scintillator material through various reaction channels, each of which has a different cross section. This leads to a wide range of light output distribution for a given neutron energy, making it difficult to accurately determine the energy of the neutron from the light output. The ToF of a particle hitting the neutron wall B (NWB) is given by ttof ≡ tx − t0 = tx , (4.10) where t0 = 0 is defined to be the nuclear reaction event start time at the target and tx is the arrival time at one of the NWB bars at position x along the bar. tx is not directly measured 91 in the experiment. To rewrite tx , we add Equation 4.5 and Equation 4.6, resulting in  (obs)  ℓ ttof = t − t0 − τ − , (4.11) 2vx where (obs) 1  (obs) (obs)  t ≡ tL + tR (4.12) 2 is the mean PMT time and 1 τ≡ (τL + τR ) (4.13) 2 is the mean delay time. In other words, ToF can be constructed from the difference between (obs) the mean PMT times t and t0 , plus some constant offset due to cable delays, speed of light in NE-213, and the NW bar length. 4.2.1 Using Forward array as a reference for event start time In the experiment, the exact event start time t0 is not directly accessible. Instead, a Forward array (FA) is installed 12 cm downstream after the target, with its central hole aligned with the beam axis, to give a reference time for t0 . A moment after the beam interacts with the target, some charged particles from the collision may reached the FA and trigger a timing signal at tFA . For different events in a run, the values of tFA may be different due to the stochastic nature of the beam-target collision. However, due to the relative proximity of the FA to the target (∼ 12 cm, i.e. 0.4 light-nanosecond) when compared to the distance between the target and the NWB (∼ 450 cm, i.e. 15 light-nanosecond), we may assume the distribution of tFA has a small variance that allows it to be treated as a constant. In other words, it is more useful for us to write ToF as   (obs) ttof = t − tFA + t′ (4.14) instead of Equation 4.11, where t′ is a constant offset that comprises τ and ℓ/(2vx ) as well as the constant delay of tFA with respect to t0 that we have assumed. The FA has 18 scintillators. For each event, the average number of scintillators that 92 trigger is 12.7, with over 90% of the events having at least five scintillators triggering. To compute a single tFA value for each event, we take the average of the triggered scintillator times excluding the extreme values, i.e. scintillator times that are later than 3 ns after the earliest scintillator time in an event. Individual FA scintillators may suffer from the “time walk” effect, which manifests as a shift in the timing of the scintillator signal depending on the amplitude of the signal itself. This effect is particularly significant when dealing with small signals. These timing signals affected by the time walk effect exhibit an undesirable dependence on the signal amplitude for time-of-flight calibration. In our setup, where the flight time of gamma rays detected in NWB is fixed by the distance to the wall and the speed of light, the time walk effect, causing earlier triggers for larger signals, needs to be accounted for and corrected. In Figure 4.6, we plot the prompt gamma events that begin emission at t0 . The x-axis is (obs) the difference between the tNWB and tFA , where tNWB is the mean PMT time t observed at the NWB; the y-axis is the Total ADC values of the scintillator. Before the time walk correction (left panel), the distribution is skewed toward increasing tNWB − tFA as the ADC Figure 4.6 Forward array time walk correction. The left panel and the right panel show the tNWB − tFA distribution before and after the time walk correction, respectively. 93 value increases. This skewed distribution can be fitted with [123] 1 (tNWB − tFA ) − a ADC = ln , (4.15) c b where a, b, c are the fit parameters. After the time walk correction (right panel), the distribution is independent of the ADC value. 4.2.2 Calibrating ToF using prompt gamma events In the analysis, we assume the calibrated ToF for NWB is given by ttof = (tNWB − tFA ) + t′ , (4.16) where t′ is the time offset to be determined from calibration, and tNWB − tFA is known as the “uncalibrated ToF”. To determine the value of t′ , we rely on the prompt gamma events that begin emission at time t0 . The prompt gamma peak is easily identifiable in a ToF spectrum, even if the ToF is not fully calibrated. In Figure 4.7, we present an uncalibrated ToF spectrum for all hits in NWB bar 12 from run 4090, including the gamma rays, neutrons, and charged particles. 600 500 400 Counts 300 200 100 00 10 20 30 40 50 Uncalibrated ToF [ns] Figure 4.7 Time-of-flight spectrum of all hits in NWB bar 12 from run 4090. The first prominent peak at around 25 ns is the prompt gamma peak. The red dashed is the Gaussian fit. 94 Earlier events are on the left and later events are on the right. The first prominent peak around 25 ns, fitted with a red dashed Gaussian curve, corresponds to the prompt gamma peak. Given that the speed of gamma rays in air closely approximates c, the speed of light in a vacuum, owing to the nearly transparent nature of air to these high-energy photons, we proceed with the following constraint: ∥r∥ ttof = , (4.17) c where r is the position vector of the hit position with respect to the target; ∥r∥ is the “distance-of-flight” (DoF). It follows that the time offset t′ can given by ∥r∥ t′ = − (tNWB − tFA ) (4.18) c from the prompt gamma events. The shortest distance between the target and the front surface of NWB is 441.6 cm. Using this distance as the DoF, we get t′ ≈ 14.7 ns. In other words, the prompt gamma peak in Figure 4.7 should appear at ∼ 15 ns instead of ∼ 25 ns if the ToF is fully calibrated. In general, the DoF depends on the hit position on the NWB, e.g. hits near the center of the NWB will have shorter DoF than hits near the edges or corners of the NWB. In fact, a wide range of DoF values are observed in Figure 4.7, ranging from 441.7 cm to 460.0 cm. In Figure 4.8, each panel represents a ToF spectrum gated with a sub-range of DoF values. As DoF increases, the prompt gamma peak shifts to the right according to the constraint in described Equation 4.17. Using the fit results from various DoF sub-ranges, we can fit a straight line to the prompt gamma peak positions as a function of DoF to determine the time offset t′ for every NWB bar in every run, as shown in Figure 4.9. The straight line has a slope fixed to the reciprocal of the speed of light in air, 1/c. Only the intercept is a free parameter, and it is found to be t′ ≈ 10.02 ns for NWB bar 12 in run 4090. 95 DoF = 446.0 ± 4.0 cm 150 DoF = 447.2 ± 4.0 cm DoF = 448.4 ± 4.0 cm 150 125 125 125 100 100 100 Counts 75 Counts 75 Counts 75 50 50 50 25 25 25 0 0 0 120 DoF = 449.6 ± 4.0 cm DoF = 450.8 ± 4.0 cm DoF = 452.0 ± 4.0 cm 80 100 60 80 60 Counts Counts Counts 60 40 40 40 20 20 20 0 0 0 DoF = 453.2 ± 4.0 cm 40 DoF = 454.4 ± 4.0 cm DoF = 455.6 ± 4.0 cm 50 30 40 30 25 20 Counts Counts Counts 30 20 15 20 10 10 10 5 020 22 24 26 28 30 32 020 22 24 26 28 30 32 020 22 24 26 28 30 32 Uncalibrated ToF [ns] Uncalibrated ToF [ns] Uncalibrated ToF [ns] Figure 4.8 Gaussian fits of prompt gamma peaks for NWB bar 12 in run 4090. Each panel represents a different gate on the distance-of-flight (DoF). The red dashed lines are the Gaussian fit results. Notice that the y-axis is scaled differently for each panel. The same procedure is repeated for all NWB bars in all runs. The final time offset values are determined by averaging the t′ values from runs within the same chunk of runs. 96 25.5 25.4 25.3 Uncalibrated ToF [ns] 25.2 25.1 25.0 24.9 24.8 24.7 446 448 450 452 454 456 Distance of flight [cm] Figure 4.9 Linear fit of all Gaussian fits of the prompt gamma peaks for NWB bar 12 in run 4090. The black solid points indicate the prompt gamma peaks. The green solid line is the linear fit with a fixed slope of 1/c and a variable intercept. 4.3 Pre-processing of analog-to-digital converter values In the previous two sections, we focused on hit position and time-of-flight calibration using timing information from the time-to-digital converter (TDC) modules. Starting from this section, we introduce calibration procedures that utilize not only the timing information, but also the light output signals. In particular, we discuss the pre-processing of the analog-to- digital converter (ADC) values that will be used for neutron-gamma pulse shape discrimination (PSD) and light output calibration. These two topics will be discussed in the subsequent sections, section 4.5 and section 4.4, respectively. In order to extract accurate measurements from the NE-213 scintillator bar, it is necessary to pre-process the ADC values obtained from the photomultiplier tubes (PMTs) attached to the scintillator. The ADC values are generated from the photoelectrons emitted by the PMT when light from the scintillator. A simple relation between the light L detected by the PMT cathode and the output voltage V from the PMT as a function of time is given by [158] V (t) = AGηLe−s/λ · f (t) + V0 , (4.19) 97 where A is a proportionality constant that accounts for the geometric factor, optical coupling to the scintillator, and electronics, G is the PMT gain, η is the PMT quantum efficiency, s is the distance between the light source in the scintillator and PMT cathode, λ is the attenuation length of the scintillator, V0 is the baseline voltage, and f (t) represents the time-dependent response of the PMT which gives the pulse shape observed in an oscilloscope. A charge-to-digital converter (QDC) module is used to integrate the pulse over a certain gate width and convert the integrated charge into an ADC value. The Fast gate (Q1 ) integrates the first 30 ns of the pulse, while the Total gate (Q2 ) integrates over a range of 340 ns to capture the entire pulse. Figure 4.10 shows the pulse shapes collected using digital electronics with a sampling rate of 500 MHz, which was installed in parallel with the analog electronics during the experiment. For the neutron analysis presented in this dissertation, we rely on the analog signals, and only use digitized signals for benchmarking and testing of the newly-built digital electronic module [130, 129]. In this experiment, we use CAEN V792 [128] for the Fast ADC and CAEN V862 [127] Fast Q1 Total Q2 1.0 (gate width of 340 ns) 0.8 Normalized ADC 0.6 0.4 0.2 0.0 0 20 40 60 80 Relative time from 20% threshold (ns) Figure 4.10 Fast gate and Total gate on a pulse. The pulse shapes are collected using digital electronics with a sampling rate of 500 MHz. The solid line indicates the median pulse shape interpolated with a cubic spline, and the light shaded area indicates the inter-quartile range. 98 for the Total ADC. Both QDC modules have 12-bit resolution, i.e. the ADC values range from 0 to 4095, inclusively. This finite range of ADC values means that for high enough light output, the ADC value will saturate at 4095, and no further increase in the light output can change the ADC value. Moreover, the ADC values stop responding linearly to the light output beyond ∼ 3800. Both of these effects need to be corrected before proceeding to the next steps of the analysis. 4.3.1 Non-linearity correction To obtain accurate measurements of the light output from the scintillation bars, it is important to correct for non-linearities in the response of the QDC modules. The non-linearity of ADC values is evident in Figure 4.11, which displays the Fast and Total ADC values for NWB bar 13 in run 4224. The left panel shows data points before pedestal subtraction, with values below their respective pedestals removed. The pedestals for the Fast gates range from 20 to 80, while those for the Total gates range from 170 to 300, (L) (L) varying across all PMTs from different NWB bars. The relation between Q2 and Q1 is Figure 4.11 Fast v.s. Total ADC values for NWb bar 13 in run 4224 before and after pedestal subtraction. The red solid line in the right panel represents a piecewise polynomial fit to the data points. 99 (L) (L) linear until Q2 ≈ 3800, after which saturation sets in at Q2 = 4095, leading to a vertical accumulation of data points. The right panel shows the same data points shifted down by their respective pedestals, except for the saturated data points at 4095, which create a visible gap between the non-linear and saturated data points to aid distinction in later analysis since they require different corrections. To correct for the non-linearity, an empirical piecewise polynomial function is used. The data points, excluding the saturated ones, are fitted with this function, which is composed of a linear function, a quadratic function, and a constant function as given by   a0 + a1 Q1    if Q1 < Qnonlinear   P(Q1 ) = b0 + b1 Q1 + b2 Q21 if Qnonlinear ≤ Q1 < Qstationary , (4.20)      c 0  if Q1 ≥ Qstationary where a0 , a1 , b0 , b1 , b2 , and c0 are the coefficients of the polynomial functions to be determined, Qnonlinear is the light output at which the non-linearity starts, and Qstationary is the stationary point of the quadratic function. Not all the parameters are independent. Continuity and smoothness of the polynomial function are enforced at Qnonlinear , resulting in b0 = a0 + b2 Q2nonlinear , (4.21) and b1 = a1 − 2b2 Qnonlinear . (4.22) The stationary point of the quadratic function can be determined by setting the derivative of the quadratic function to zero, resulting in b1 Qstationary = − (4.23) 2b2 and b21 c 0 = b0 − . (4.24) 4b2 To summarize, only four independent parameters are varied in the fit, namely a0 , a1 , b2 , and 100 Qnonlinear . To improve the numerical stability of the fit, a0 and a1 are first determined by fitting the data points below Q1 ≤ 3000. Then, b2 and Qnonlinear are varied to obtain the best fit for all the data points excluding the saturated ones. Finally, to correct for the nonlinear Total ADC values Q2 , we assume that the Fast ADC values Q1 are linear. The Total ADC values Q2 are then corrected by Q2 → Q2 + [(a0 + a1 Q1 ) − P(Q1 )] . (4.25) 4.3.2 Saturation correction For every NW bar, two PMTs are attached to the bar, one on the left side and the other on the right side. To correct for the saturated Total values, we use the Total values collected from the PMT at the opposite end of the same bar. Saturation happens less often for the Fast values, and we do not correct for the Fast values. In order to use the Total values from the opposite PMT for saturation correction, we need to first understand the relation between the Total values from the two PMTs attached Figure 4.12 The  logarithm  of the quotient of the Total values from the right and (R) (L) left PMTs, ln Q2 /Q2 , as a function of the hit position x for NWB bar 5 in run 4300. The black dashed line is a linear fit to the data points. 101 to the same bar. Generally, the performance of the two PMTs is not identical. Within the linear range of the ADC values, the gate values from the left and the right PMTs can be expressed as   (L) (L) ℓ/2 + x (L) Q =g L exp − + Q0 (4.26) λ and   (R) (R) ℓ/2 − x (R) Q =g L exp − + Q0 , (4.27) λ where ℓ is the full length of the NW bar, x is the hit position, g is the overall gain factor, and Q0 is the pedestal value. In the analysis, the pedestal values are subtracted from the gate values, so Q0 may be dropped from the equations. Next, to measure the attenuation length λ, we take the logarithm of the quotient between Equation 4.27 and Equation 4.26, which yields Q(R) g (R)     2 ln = x + ln . (4.28) Q(L) λ g (L) This gives a linear relation between the hit position x and the logarithm of the ratio between the gate values from the left and the right PMTs, as shown in Figure 4.12. In this plot, (L) (R) only the data points with Q2 ≤ 3000 and Q2 ≤ 3000 are used. We derive an attenuation length of λ = 105.7 ± 0.6 cm for NWB bar 5 in run 4300, using the slope 2/λ. This measured “technical attenuation length” (TAL) [159] is considerably shorter than the nominal bulk attenuation length of 250 cm to 300 cm [160]. The primary cause of this discrepancy is the geometry of the NW bars. Due to the light internally reflecting and thus covering more path to travel the same longitudinal distance, the effective attenuation length along the longitudinal direction becomes shorter than the bulk value. Other contributing factors include reflectivity variations and the age-induced conditions of the bars, such as air contamination and deterioration of the scintillator material, given the NW bars have been in commission for over thirty years. Aside from the slope, the intercept is related to the ratio of the gain factor g (R) /g (L) , which is found to be 1.006 ± 0.003. A value close to 1 indicates that the gain matching between the two PMTs is good. (R) (L) With Equation 4.28, the ratio ρ(x) ≡ Q2 /Q2 can be given as a function of position x. 102 We then recover the saturated Total values according to (L) (R) Q2 → Q2 /ρ(x) (4.29) and (R) (L) Q2 → Q2 · ρ(x) , (4.30) assuming that the left PMT or the right PMT is saturated, respectively. In the case when (L) (R) both Q2 and Q2 are saturated, recovery is not performed. Saturation correction is an essential step in the light output calibration process. Due to the presence of attenuation in the scintillation material, the detected Total values from the left and the right PMTs are dependent on the hit position x, as modeled by Equation 4.26 and Equation 4.27. To remove the dependence on x, we take the square root of the product of these two equations with the pedestal subtracted, which yields p p Q(GM) ≡ Q(L) Q(R) = g (L) g (R) · Le−ℓ/λ . (4.31) This quantity is also known as the “geometric mean” (G.M.) of the two values. Notice that g (L) , g (R) , ℓ, and λ are all constants for a given bar and the PMTs attached to it. This means that the geometric mean of Q(L) and Q(R) , Q(GM) , is proportional to the light output L. It is important to note that this complete removal of position-dependence is only valid based on the simple model given by Equation 4.26 and Equation 4.27. In the next section, section 4.4, we will discuss about the light output calibration in more detail and show that Q(GM) removes only the first-order position dependence. In Figure 4.13, we show the Total values from the left and the right PMTs before non- linearity and saturation corrections. The one-side saturated Total values are concentrated (L) (R) into the vertical line at Q2 = 4095 and the horizontal line at Q2 = 4095. The two-side (L) (R) saturated Total values are concentrated into one single point at (Q2 , Q2 ) = (4095, 4095), and it has been enlarged to make it visible in the plot. The vertical and horizontal gaps between the saturated points and the rest of the data points are due to pedestal subtraction. 103 Figure 4.13 The left Total v.s. the right Total for NWB bar 12 from run 4224 to run 4245, inclusively, before non-linearity and saturation corrections. The black dotted curves indicate the contours of Q(GM) . The two-side saturated point is enlarged to make it visible. The black dotted contours indicate the Q(GM) values at 1000, 2000, 3000, 4000, corresponding to calibrated light output of ≈ 14 MeVee, 28 MeVee, 42 MeVee, 56 MeVee, respectively. As can be seen along the contour of Q(GM) = 2000, even for light output as low as ∼ 28 MeVee, both non-linearity and saturation corrections are needed to construct an accurate light output distribution. In Figure 4.14, we show the same plot, but after non-linearity and saturation corrections. Both of these corrections extend the continuum of the Total values beyond the QDC’s finite dynamic range that ends at 4095. The top-right quadrant of the plot remains empty because it corresponds to the two-side saturated Total values that cannot be recovered easily. Nonetheless, we still have partial access to data points of Q(GM) > 4095 using the corrected Total values, giving us access to light output up to ∼ 100 MeVee. 104 Figure 4.14 The left Total v.s. the right Total for NWB bar 12 from run 4224 to run 4245, inclusively, after non-linearity and saturation corrections. The corrected Total values will be used in all subsequent analyses, and any two-side saturated Total values will be discarded unless their inclusion is necessary for a particular analysis. Specifically, in subsection 5.2.3, the two-side saturated Total values will be addressed as we construct light output distribution curves for determining intrinsic efficiency. 4.4 Light output calibration The mechanisms responsible for energy loss in charged particles and neutrons when interacting with scintillator materials are different. Charged particles primarily lose energy through ionization and excitation of atoms. In contrast, neutrons have multiple reaction channels with differing cross sections. These channels can generate charged particles such as protons, deuterons, and alpha particles, which deposit energy and trigger light output. Due to the indirect nature of this energy deposition process, the light output from neutrons is less predictable than that from charged particles. Consequently, no straightforward relationship 105 can be drawn linking the incident energy of the neutron to the light output. A more thorough examination of the light output mechanisms and their simulation will be provided in subsection 5.2.3, where the light output distribution will be utilized to determine the intrinsic efficiency of the NWB. The relationship between energy deposition and light output is dependent on the type of particle involved. This implies that the light output response to a neutron scattering and subsequently producing a proton will differ from the response of a neutron scattering and generating an alpha particle. To standardize the calibration process for a neutron detector, a standard unit is necessary. The convention of denoting light output in MeV electron equivalent, or MeVee, has been widely adopted. MeVee signifies the energy deposited by a particle that yields the same light output as a 1 MeV electron. 4.4.1 Removal of position dependence As mentioned in Equation 4.31, the geometric mean of the left and right Total values, Q(GM) , can be related to the light output L by p p Q(GM) ≡ Q(L) Q(R) = g (L) g (R) · Le−ℓ/λ (4.32) without any position dependence. However, this is only true for a simple model that assumes the NW bar is essentially one-dimensional along the x direction, and the attenuation of light strictly follows the Beer-Lambert law. In reality, the NW bar is not a one-dimensional object, hence position dependence would still present in Q(GM) . Following the calibration technique using Compton kinematics as discussed by E. Siciliano et al. [161], we place an Americium-Beryllium (AmBe) neutron source ∼ 60 cm in front of the center of the NWB during run 3072 and run 3073. The source emits gamma rays of Eγ = 4.438 MeV, which are associated with the neutrons leaving 12 C in the first excited state from 9 Be(α, n)12 C [162]. This gives us a Compton edge at 2Eγ2 ≈ 4.2 MeV , (4.33) me c2 + 2Eγ where me is the electron mass and c is the speed of light. 106 Figure 4.15 The top panel shows the position dependence of Q(GM) for NWB bar 10 from run 3072 and run 3073. The black dotted curve indicates Compton edge of AmBe. The bottom panel shows Q(GM) after the position dependence is removed. In Figure 4.15, we show how the position dependence of Q(GM) is removed for NWB bar 10 from run 3072 and run 3073. The top panel shows a shallow U-shape dependence of Q(GM) on the position x. The curvature of the U-shape is fitted with a quadratic function, (GM) q(x) = a0 + a1 x + a2 x2 , on the Compton edge in the range of 300 < Q2 < 400. This allows us to flatten the U-shape and remove the position dependence of Q(GM) by (GM) (GM) (GM) Q2 → Q2 − [q(x) − q(0)] = Q2 − (a1 x + a2 x2 ) . (4.34) The bottom panel shows the flattened distribution. There are more data points at the center (x ≈ 0 cm) than the edges (x ≈ ±80 cm) because the source is placed at the center of the NWB. 107 4.4.2 Calibration of light output (GM) While it is possible to calibrate Q2 to MeVee using the Compton edge from the AmBe source used in Figure 4.15 by setting the fitted Compton edge line to 4.2 MeVee, we want to obtain some calibration points at higher light outputs at above 10 MeVee, which can be accessed easily with cosmic muons [163, 164]. Cosmic muons are high-energy charged particles that are produced by the interaction of cosmic rays with the Earth’s atmosphere. At sea level, the average energy of cosmic muons is ∼ 2 GeV [165, 166, 167]. These particles can penetrate several meters of material and deposit energy along their path, making them useful for calibrating detectors. The energy deposited by a cosmic muon in a scintillator material is proportional to the path length in the material. Using this property, we simulation the interaction of 2 GeV cosmic muons with the NWB using Geant4 and its default list of physics models [168, 169] at various incident angles. Specifically, we simulation muons entering the NWB at incident angles of 0◦ ± 5◦ , 44.4◦ ± 5◦ , and 56.3◦ ± 5◦ with respect to the gravity direction. According to the simulation, muons with these incident angles are expected to yield light outputs of 11.96 MeVee, 17.05 MeVee, and 22.16 MeVee, respectively. These values are much higher than the light output of 4.2 MeVee from the Compton edge of the AmBe source. More details about the simulation and the construction of cosmic muon tracks are discussed in [123, 122]. In Figure 4.16, we show the calibration curve of flattened Q(GM) as defined in Equation 4.34 into MeVee for NWB bar 12. There are five calibration points for every NWB bar, namely, the pedestal, the Compton edge from the AmBe source, and the three cosmic muon points. A simple linear regression is performed on the data and plotted as the green solid line. This calibration line is then used to convert flattened Q(GM) into MeVee in the analysis. 108 fit muons at ∼ 0 ◦ ◦ 25 pedestal muons at ∼ 44.4 ◦ AmBe muons at ∼ 56.3 Light output L [MeVee] 20 15 10 5 0 0 200 400 600 800 1000 1200 1400 1600 Flattened Q2(GM) [ADC] Figure 4.16 The calibration curve of flattened Q(GM) (position dependence removed) into MeVee for NWB bar 12. 4.5 Pulse shape discrimination Although LANA was primarily designed to detect neutrons, the NE-213 scintillation material it uses is also sensitive to other types of radiation, such as charged particles and gamma rays. In subsection 2.3.3, we discussed how charged particles can be eliminated from the analysis by using the Veto Wall (VW). In this subsection, we discuss how pulse shape discrimination (PSD) is used to remove gamma rays from the neutron analysis. NE-213 is a common organic liquid scintillation material used in neutron detection because of its high performance in n-gamma PSD. The material is based on solvent xylene (C8 H10 ) and emits light when excited by radiation [170, 171]. The light pulses exhibit different fall-times depending on the type of radiation—a neutron or a gamma ray[172, 173]. This radiation-dependent property of NE-213 forms the basis of n-γ PSD. During the experiment, 7 out of 25 bars from NWB were selected for split signals from the photomultiplier tubes (PMTs). These bars are numbered as 01, 02, 03, 10, 11, 12, and 13. This was done because the digital electronic system was new at the time [130], and it was 109 necessary to test it while keeping the more reliable analog electronic system as the primary data acquisition source. Hence, 90% of the signal strength was sent to the analog electronic system, while the remaining 10% was sent to the digital electronic system. Pulse shapes recorded via the digital electronic system allows for high flexibility in the offline analysis. This flexibility is beneficial in this work because it allows us to develop a new algorithm to optimize PSD in a detector with a long geometry [129], such as the Neutron Wall bars used in this experiment, each has a length of 76” or 193.04 cm and a cross section of 2.5” × 3” or 6.35 cm × 7.62 cm. Teh et al. [129] discuss the development and benchmarking of the new PSD algorithm primarily based on the digitized signals and briefly mention its applications on analog signals, whereas this thesis will emphasize the applications on analog signals, which are the primary signals used for experimental data acquisition. Light pulses emitted from NE-213 exhibit different shapes in response to neutrons and gamma rays, as demonstrated in Figure 4.17. The majority of differences in pulse shapes occur Fast Q1 Total Q2 1.0 (till the end of pulse) 0.8 Normalized ADC 0.6 0.4 neutron 0.2 gamma 0.0 0 20 40 60 80 Relative time from 20% threshold (ns) Figure 4.17 The median digitized waveforms collected by the left PMT of NWB bar 01 from AmBe source. Only hits on the closest 50cm to the left PMT are included. The shaded regions indicate the inter-quartile range of the normalized ADC. Replotted from [129]. 110 after reaching the peak. The oscillating waveforms towards the end of the tails are caused by the splitting of signals into analog and digital electronic systems, as well as amplification via the Fast Amps for signal strength recovery. While it is possible to utilize the entire waveform for PSD, it is commonly done by analyzing the difference between two charge-integration gates, namely “Fast” (Q1 ) and “Total” (Q2 ) gates, over different time ranges when using an analog electronic system [98, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183]. (L) (R) Four ADC signals are collected from each NW bar, namely, Q1,2 and Q1,2 , with superscripts indicating the sides of the PMT where the signals are collected. While it is possible to identify neutrons and gamma rays by looking at either side of the PMT, as shown in Figure 4.18, combining information from both sides of the PMTs and the hit position of radiation x would significantly improve the accuracy of n-γ identification. Here, the hit positions have been calibrated according to section 4.1, such that x = 0 cm would correspond to the center of the bar. This is a two-class clustering problem in a five-dimensional space. The goal is to find a Figure 4.18 This figure shows the two-dimensional (Q2 , Q1 ) histogram of NWB bar 13 in run 4224 to run 4245, inclusively. The left panel shows the full range of the histograms. The right panel shows the fitted curves for neutrons and gamma rays, with the vertical axis flattened according to Equation 4.36 for better numerical stability and visualization. Only hits in the range of −95 cm < x < −50 cm are included. 111 (L) (L) (R) (R) function that reduces the dimensionality of (Q1 , Q2 , Q1 , Q2 , x) and maps it to a scalar parameter, also known as the PSD parameter. First, we reduce the dimensionality of the (Q2 , Q1 ) data pair for signals from both PMTs. This is done by simultaneously fitting two curves onto the neutron and gamma ray ridges in the (Q2 , Q1 ) histogram. We denote the fitted curves for neutrons and gamma rays as Pn and Pγ , respectively. Then, for every signal, we assign a scalar value according to Q1 − Pγ (Q2 ) v(Q2 , Q1 ) ≡ . (4.35) Pγ (Q2 ) − Pn (Q2 ) The quantity v, as well as the procedure, is referred to as the “valued-assigned pulse shape discrimination” (VPSD). This construction of v assigns values close to 1 and 0 to (Q2 , Q1 ) points that are close to the neutron and gamma ray ridges, respectively. In practice, simultaneous fitting of two curves is a challenging algorithm to implement directly on the (Q2 , Q1 ) histogram (left panel of Figure 4.18) due to the small separation between the neutron and gamma ray ridges. To ensure the neutron and gamma ray ridges are well separated, we only select hits up to 45 cm away from the PMT; hits further away from the PMT suffer from lower signal-to-noise ratio due to the attenuation, and thus are more difficult to be distinguished. Furthermore, we fit the curves in (Q2 , q1 ) space (right panel of Figure 4.18) instead of the (Q2 , Q1 ) space, where q1 ≡ F(Q1 ) and F is a “flattening” function defined as q1 ≡ F(Q1 ) ≡ Q1 − L(Q2 ) , (4.36) and L is a simple linear fit of the original (Q2 , Q1 ) histogram. This technique [152, 123] facilitates a more robust fitting procedure by subtracting out the correlation between Q1 and Q2 , hence it was given the name “flattening”. The fitted curves in (Q2 , q1 ) space are denoted as pn and pγ for neutrons and gamma rays, respectively, and they are related to Pn and Pγ by pn,γ (Q2 ) = Pn,γ (Q2 ) − L(Q2 ) . (4.37) 112 Using Equation 4.36 and Equation 4.37, we can rewrite Equation 4.35 into q1 − pγ (Q2 ) v(Q2 , q1 ) = . (4.38) pγ (Q2 ) − pn (Q2 ) VPSD parameters calculated from only one single PMT, v (L) or v (R) , can already be used as PSD parameters by themselves. But putting both v (L) and v (R) together in a two- dimensional histogram, as shown in Figure 4.19, reveals that the neutron and gamma ray clusters can be separated even better. The second step is to correct for the hit position dependence of v ≡ (v (L) , v (R) ). This is (L) (R) done by calculating the centroids of neutron clusters, cn ≡ (cn , cn ), for various hit position ranges. In Figure 4.19, the red dashed contours and the blue solid contours indicate the neutrons clusters at x ≈ −80 cm and x ≈ +80 cm, respectively, whereas the black dotted line traces the centroids of neutron clusters over the entire length of a NW bar, from x = −100 cm to x = +100 cm. Similar trace line has been constructed from gamma rays, but it is not (L) (R) shown in the figure due to small displacement of the gamma ray centroids, ⃗cγ ≡ (cγ , cγ ). Position dependence can then be corrected according to both of the trace lines, cn (x) and cγ (x), reducing the spreads of neutron and gamma ray clusters. Mathematically, this is done by transforming (v (L) , v (R) , x) into (νx , νy ), where ν (1) ≡ (v − cγ (x)) · ĉnγ (x) , (4.39) and   0 −1 ν (2) ≡ (v − cγ (x)) ·   ĉnγ (x) . (4.40) 1 0 Here, cnγ (x) ≡ cn (x) − cγ (x) is the relative position vector from cγ (x) to cn (x), and ĉnγ is the corresponding unit vector. By construction, (νx , νy ) forms an orthogonal coordinate system, with νx encoding most of the n-γ separation from the original (v (L) , v (R) ) space, and νy preserving the information in the orthogonal direction. To finalize the position correction, we normalize ν (1) such that the neutron and gamma ray clusters are centered at 1 and 0, respectively. This convention is chosen to be consistent 113 Figure 4.19 The (v (L) , v (R) ) histogram includes all hits in the range of x ∈ [−95, +95] cm. The red dashed contours and the blue solid contours indi- cate the hits in the range of [−90, −70] cm and [+70, +90] cm, respectively. The black dotted line, denoted as ⃗cn (x), traces the centroids of neutron clusters from x = −100 cm to x = +100 cm. with the VPSD parameters. To normalize ν (1) , we identify the neutron and gamma ray peaks (1) (1) in the one-dimensional ν (1) histogram; the peaks are denoted as νn and νγ , respectively. Finally we define “position-corrected VPSD” (PPSD) as (1) ν (1) − νγ PPSD ≡ (1) (1) . (4.41) νn − νγ There is no need to normalize ν (2) ; we shall call it the “perpendicular PPSD”. The final PSD histogram is shown in Figure 4.20, in which only hits with light output above 3.0 MeVee are included. A good PPSD value to separate the neutrons and the gamma rays is around 0.5. PPSD > 0.5 is a neutron hit, and PPSD ≤ 0.5 is a gamma ray hit. 114 Figure 4.20 The final PSD histograms for NWB bar 13 in run 4224 to run 4245, inclusively. The top panel shows the distribution of the PSD parameter, PPSD. The bottom panel shows the distribution in two dimension, including the perpendicular PPSD in the vertical axis. 115 CHAPTER 5 RECONSTRUCTION OF NEUTRON SPECTRA In this chapter, we will describe the construction of neutron spectra using calibrated mea- surements from the neutron detection system, vLANA (LANA with Veto Wall). We will begin by discussing the challenges of background analysis in neutron measurements and the techniques for estimating and subtracting background contributions. We will then explore the importance of efficiency corrections in neutron spectroscopy and the methods used for determining the system response function. Finally, we will present the neutron spectra and discuss the associated uncertainties. 5.1 Background analysis There are two main sources of background in the neutron detection. First, there are scattered neutrons from the surroundings of the experimental setup. Second, there are random neutrons produced at the radio frequency (RF) of the K1200 cyclotron at NSCL. 5.1.1 Scattered neutron background When neutrons scatter off materials from the surroundings of the experimental setup, e.g. concrete, metal, etc., they travel to NWB via indirect paths. These indirect paths would result in time-of-flights (ToF) that are longer than the direct paths, i.e. the straight paths defined by the target’s center and the hit position on the NWB. Because the exact path lengths for indirect paths are not known, this makes it impossible to calculate the associated kinetic energies of the neutrons. These undesirable scattered neutrons make up the scattered background, which must be subtracted from the neutron yield. To measure the scattered background, we insert four brass shadow bars between the walls and the target to stop direct neutrons. Figure 5.1 shows a photograph of one of the shadow bars. Figure 5.2 shows a photograph of the experimental setup with the shadow bars inserted. The shadow bars are placed around 170 cm away from the target, outside the vacuum chamber. They are oriented such that the longitudinal axis is pointing towards 116 Figure 5.1 A photograph of one of the shadow bars, taken with a ruler placed next to it for scale. Adapted from [184]. Figure 5.2 A photograph of the experimental setup with the shadow bars inserted. The four shadows are labeled as A, B, C, and D. the target. The thickness of the shadow bars is 30 cm, which is sufficient to block direct neutrons with energy up to 300 MeV from the target. Neutron shadows would be casted onto the NWB. The scattered background can then be measured by comparing the number of neutrons within and outside the shadows. Figure 5.3 shows the hit pattern of neutrons on the NWB with all four shadow bars inserted. In this plot, charged particles, gamma rays, and neutrons with light output below 117 Figure 5.3 Hit pattern of neutrons on the NWB with the shadow bars inserted. The four shadows are labeled as A, B, C, and D. 3.0 MeVee have been removed. During the experiment campaign, about half of the runs were taken with the shadow bars inserted, and the other half without. The four shadows are label as A, B, C, and D. Shadows A and B affect NWB bars 15, 16, and 17; shadows C and D affect NWB bars 7, 8, 9, and 10. In the analysis, if shadow bars are inserted, we remove neutrons that hit on xbar ∈ [−50 cm, −15 cm] (shadows A and C) and xbar ∈ [15 cm, 50 cm] (shadows B and D) in the affect bars. In order to quantify the scattered background, we will look at the θ spectrum of neutrons for each NWB bar individually. Inside the shadowed region, we assume the θ spectrum exhibit a flat distribution. Outside the shadowed region, the observed neutron yield would be a sum of both the scattered background and the actual kinematics. In practice, the shadow edges are smeared out due to the finite position resolution of the NWB of about 7 cm. This suggests us to construct a Gaussian convoluted flat well function, which is given by      θ − x0 θ − x1 F (θ) = A 1 + erf − √ + erf + √ + (b + cθ) . (5.1) 2σ0 2σ1 118 By varying A, x0 , x1 , σ0 , σ1 , b, c, Equation 5.1 can be used to fit one single shadow on a NWB bar. x0 and x1 are the shadow edges. We then use the error function erf to model the Gaussian smeared edges, with σ0 and σ1 being the widths of the Gaussian. Finally, A indicates the depth of the shadow, and b + cθ term is used to model the spectrum outside the shadowed region. Figure 5.4 An empirical fit to shadow A on NWB bar 16. The data come from runs 4042 to 4078, inclusively, which is 48 Ca +124 Sn at 140 MeV/A. Top panel: The scatter points are the data, and the lines are fits. Bottom panel: The same fits normalized by Ykin (θ). Outside the shadowed region, F (θ) → A + b + cθ. We shall denote this as Ykin (θ) ≡ A + b + cθ . (5.2) Figure 5.4 shows various fits using Equation 5.1 on the data. The data points in the top panel are selected from NWB bar 16 at the vicinity of shadow A in runs 4042 to 4078, inclusively. 119 The upper limit and the lower limit curves are estimated by varying the fit parameters within their uncertainties. In the bottom panel, the same fits are normalized by Ykin (θ). Hence, outside the shadowed region, the normalized curves should be flat at unity. The scattered neutron background is then given by the normalized value at the center of the shadow, which is (x0 + x1 )/2. In this particular plot, we measure the scattered background to be 38.2 ± 0.8%. In Figure 5.4, we have included neutrons from all energies, En . However, the scattered background is expected to be energy dependent. Since the neutron energy is calculated using the time-of-flight (ToF) in this experiment, the shorter the ToF, the higher the neutron energy. Neutrons measured with shorter ToF are less likely to be scattered. This is because scattered neutrons would need to travel a longer distance to reach the NWB. In order to reach the NWB with a short ToF, the scattered neutrons would need to be traveling at a higher speed, hence a higher energy. But there tend to be fewer neutrons at higher energies. Figure 5.5 The scattered neutron background curves as functions of neutron energy En . The left panel shows reaction systems with beam 40 Ca, and the right panel shows reaction systems with beam 48 Ca. Figure 5.5 shows the scattered neutron background curves as functions of neutron energy En . As expected, the scattered background is energy dependent. As the neutron energy En increases, which corresponds to a shorter ToF, the scattered background decreases. There 120 are significant scattered neutrons at low energies. 5.1.2 Random background at radio frequency In this experiment, the beam is delivered by the cyclotron at the National Superconducting Cyclotron Laboratory (NSCL). The cyclotron accelerates charged particles using an alternating electric field, while a constant magnetic field guides them in a spiral path. As the charged particles travel through the cyclotron, they continuously gain energy and velocity until they are eventually extracted and directed toward a target, forming a high-energy beam. The radio frequency (RF) system within the cyclotron plays a pivotal role in the particle 48 acceleration process. The cyclotron operates at 23.11780 MHz (43.26 ns) for Ca beam and 40 23.22390 MHz (43.07 ns) for Ca beam. It generates an oscillating electric field, which is applied across the acceleration gap between two D-shaped electrodes known as “dees”. The RF is carefully tuned to match the time it takes for the particles to complete a half-circle within the cyclotron. This synchronization ensures that the charged particles consistently gain energy as they spiral outwards. One important characteristic of the cyclotron is that it does not deliver a continuous stream of particles. Instead, the particles are organized into distinct groups or “bunches” due to the periodic nature of the alternating electric field. As the field imparts energy to the particles at specific intervals, they become organized into bunches separated by a fixed time interval corresponding to the RF period. This process is known as cyclotron bunching. In this experiment, the Data Acquisition (DAQ) system is used to record the events generated by the interaction between the beam and the target. However, the time gate employed by the DAQ is longer than the RF period, which leads to the observation of periodic noise after the primary event that triggers the system. Due to cyclotron bunching, particles are delivered to the target in discrete groups or bunches. As these bunches interact with the target material, they may produce secondary particles or prompt radiation, which can reach the detector and create background signals. Since the bunches arrive at the target with a periodicity determined by the RF period, these background signals would also exhibit 121 periodicity in the recorded spectrum. However, these events occur independently of the primary beam-target interactions and do not show a clear correlation with the primary signals of the experiment. 0 10 2 Count density [ns−1 ] 10 3 4 5 6 7 8 9 10 4 0 50 100 150 200 250 300 350 400 450 ToF [ns] Figure 5.6 Time-of-flight (ToF) spectrum of the 48Ca + 64Ni reaction system with beam energy at 140 MeV/u. The red triangular markers indicate the prompt gamma peak, with the RF indices labeled on top of each marker. 48 Figure 5.6 displays the time-of-flight (ToF) spectrum of the Ca + 64Ni reaction system with beam energy at 140 MeV/u. This spectrum includes gamma rays and neutrons, but not charged particles. The reason for including gamma rays is their well-defined peak at the beginning of each RF cycle, which allows us to pinpoint the RF period easily. The first gamma peak, RF-0, is the prompt gamma peak that occurs at ToF ∼ 15 ns. As the speed of light is the fastest speed possible, the prompt gamma peak represents the earliest events from the target. From ToF ∼ 20 ns to ∼ 180 ns, the spectrum is dominated by primary neutrons and gamma rays from the beam-target collisions, making it difficult to identify the RF periodicity or the associated gamma peaks. Above ToF ∼ 180 ns, the RF background starts to dominate, and the gamma peaks for RF-4 to RF-9 can be clearly seen, marked by red triangular markers. These identified gamma peaks are used to measure the RF period. For each gamma peak, 122 450 400 350 RF period: 44.23 ± 0.11 ns 300 250 ToF [ns] 200 150 100 50 00 2 4 6 8 10 RF cycle index Figure 5.7 The peak positions of the gamma peaks as functions of RF indices. The green dashed line is a linear fit to the red data points, whose error bars are smaller than the marker size. we perform a Gaussian plus linear background fit, (ttof − t0 )2   fγ (ttof ; A, t0 , σt , B, m) = A exp − + (B + mtToF ) , (5.3) 2σt2 and extract the peak position t0 . These peak positions are then fitted with a linear function to relate the RF indices to the peak positions, as shown in Figure 5.7. The measured RF 48 40 periods are 44.23 ± 0.11 ns for the Ca beam and 43.84 ± 0.13 ns for the Ca beam. These values are consistent with the nominal RF periods of 43.26 ns and 43.07 ns, respectively. Upon determining the RF period PRF , we apply pulse shape discrimination (PSD) to remove the gamma rays, leaving only the neutrons. Consequently, as shown in Figure 5.8, the gamma peaks are no longer present in the spectrum. To establish an RF background profile, we fit a second-order Fourier series, 2      X 2πn 2πn fRF (ttof ) = C0 + An cos ttof + Bn sin ttof , (5.4) n=1 P RF P RF to the neutron ToF spectrum from the 5th RF gamma peak (235.9 ns) to the 9th RF gamma peak (412.8 ns), as indicated by the pink shaded region. This background profile is assumed to hold all the way down to the prompt gamma peak position, represented by the green filled 123 n spectrum RF background 10 2 Count density [ns−1 ] 10 3 10 4 0 50 100 150 200 250 300 350 400 450 ToF [ns] Figure 5.8 ToF spectrum neutrons from the 48 Ca +64 Ni reaction system with beam energy at 140 MeV/u. The green filled histogram is the periodic RF background, which is constructed by fitting a second-order Fourier series to the ToF in the pink shaded region. histogram, and is capped by the neutron ToF spectrum (blue line) to avert negative values after background subtraction. For the construction of a neutron energy spectrum, it is preferable to depict the RF background fraction as a function of neutron energy to facilitate its subtraction from the neutron energy spectrum. Figure 5.9 demonstrates the RF background fraction as functions of both ToF (top panel) and neutron energy (bottom panel). In the top panel, the RF background fraction is computed by dividing the RF background profile by the total neutron spectrum in Figure 5.8. Below ToF ∼ 20 ns and above ToF ∼ 110 ns, the RF background fraction approaches 1, implying the dominance of the RF background in these regions. The bottom panel converts the ToF axis to neutron energy En axis through relativistic kinematics, taking an average distance of 452.2 cm between the target and NWB. To evade excessive RF background contamination, a lower neutron energy threshold of 10 MeV is enforced. Above this threshold, the RF background fraction remains relatively small. There exists a slight step-down in the RF background fraction from 10 MeV to 30 MeV. Between 30 MeV and 200 MeV, the RF background fraction progressively increases from 1.9% to 3.2%, visible 124 1.0 0.8 300 MeV 10 MeV RF background 0.6 0.4 0.2 0.00 100 200 300 400 500 ToF [ns] 1.0 0.1 0.8 RF background RF background 0.05 0.6 0.4 0.0 0 50 100 150 200 250 300 Neutron energy En [MeV] 0.2 0.00 50 100 150 200 250 300 Neutron energy En [MeV] Figure 5.9 Top: RF background fraction as a function of ToF. The RF background fraction is calculated by dividing the RF background profile by the total neutron spectrum in Figure 5.8. Bottom: RF background fraction as a function of neutron energy, calculated by converting the ToF axis in the top panel to neutron energy. more clearly in the inset plot. When we construct a neutron energy spectrum later, it is more convenient to represent the RF background fraction as a function of neutron energy En . In Figure 5.9, we show the RF background fraction both as a function of ToF (top panel) and neutron energy (bottom panel). In the top panel, the RF background fraction is calculated by dividing the RF background profile by the total neutron spectrum in Figure 5.8. Below ToF ∼ 20 ns or above ToF ∼ 105 ns, the background fraction quickly approaches 100%, indicating that the RF background is dominant in these regions. In the bottom panel, we convert the ToF axis into 125 En axis using relativistic kinematics, assuming an average distance of 452.2 cm between the target and NWB. A lower threshold of 10 MeV, as indicated by the red vertical dashed line, must be enforced to avoid excessive RF background contamination. Above 10 MeV, the background fraction is relatively small. From 10 MeV to 30 MeV, there is a slight step-down in the background fraction. From 30 MeV to 200 MeV, the background fraction gradually increases from 1.9% to 3.2%, as can be seen more clearly in the inset plot. In summary, the RF background fraction is small and gradually increases with neutron energy over the energy range of interest (30 MeV to 200 MeV). 5.2 Neutron detection efficiencies The detection process of neutrons in a scintillator is more complex than that of charged particles. Neutrons, being electrically neutral, do not participate in Coulomb interactions with the electrons of the scintillator atoms. Instead, when a neutron enters the scintillator, it primarily interacts with the atomic nuclei through processes such as elastic and inelastic scattering. During these interactions, the neutron transfers some or all of its energy to the recoiling nucleus. The recoiling nucleus subsequently traverses the scintillator material, transferring its energy to the surrounding electrons via ionization and excitation. As the excited electrons revert to their ground state, they emit photons, generating the characteristic scintillation light. This light is collected by a photomultiplier tube (PMT) and converted into an electrical signal, which can then be analyzed to obtain information about the incident neutron. This detection mechanism suggests that neutrons typically exhibit a lower probability of detection compared to charged particles. In this section, we examine methods for estimating neutron detection efficiency. Neutron detection efficiency is a multifaceted concept that can be measured through various techniques [185]. In particular, we will explore three types of efficiencies. We assume the neutron source to be a point source, originating from the target. While the source emission can generally be anisotropic, it shall remain symmetric around the beam axis due to the azimuthal symmetry present in beam-target collisions. In other words, when averaging 126 over many collision events, there is no directional preference in the transverse plane. All neutrons are presumed to travel directly from the target to the neutron wall in straight lines. Scattered neutrons that reach the neutron detector through indirect paths are not considered here, as they have been subtracted as background in the previous subsection, subsection 5.1.1. The absolute efficiency, ϵabs is characterized as the proportion between the number of measured neutrons and the number of neutrons emitted from the source, expressed as Mmeas ϵabs = . (5.5) Memit In this context, “measured” means that the particle generates a detectable signal, and after applying all necessary calibrations and analysis filters, the particle is identified as a neutron. Absolute efficiency serves as the most direct definition of efficiency. However, the use of this efficiency as a sole metric for comparison with other experiments is limited because it is contingent upon experimental configurations that may be arbitrary, including factors such as detector geometry (size and shape) and the relative geometric arrangement of the detector with respect to the source. This leads us to the introduction of intrinsic efficiency, denoted as ϵint . Intrinsic efficiency is defined as the absolute efficiency divided by the solid angle Ω subtended by the detector with respect to the source, expressed as 4π Mmeas ϵint = , (5.6) Ω Memit where a 4π factor is included to ensure that the intrinsic efficiency is normalized to unity for an ideal detector that registers all incident neutrons within its active volume, irrespective of whether the volume subtends the full solid angle of 4π. Intrinsic efficiency represents the probability that a neutron passing through the active volume of the detector will generate a detectable light signal that can be identified as a neutron. Although this efficiency measure reduces the geometric dependence associated with absolute efficiency, it is not entirely independent of geometric factors. For example, the thickness of the active volume may affect the chance a neutron will interact with the scintillator material. Neutrons entering the NW bar 127 at an oblique angle will travel a longer path length through the scintillator material compared to those entering at normal incidence, increasing the probability of interaction. Nevertheless, intrinsic efficiency serves as a valuable parameter for comparing different detector systems or configurations, as it provides a more standardized measure of detection performance after normalizing by the solid angle Ω. This is also the efficiency number that is typically reported in the literature [164, 186, 187]. Lastly, we have the geometric efficiency, ϵgeo , which is defined as the solid angle Ω subtended by the detector from the source relative to 4π, expressed as Z Ω 1 ϵgeo = = sin θ dθ dϕ , (5.7) 4π 4π D where the integral is taken over the entire detector’s projection onto the (θ, ϕ) space, denoted as D. In this dissertation, we construct the absolute efficiency by determining the intrinsic and geometric efficiencies of NWB independently. Then we combine the two to obtain the absolute efficiency according to ϵabs = ϵgeo ϵint . (5.8) 5.2.1 Modifications for calculating the neutron energy spectrum The efficiencies we have discussed thus far have included neutrons of all energies and from all directions. Now, we consider a more specific scenario involving the neutron energy d2 M spectrum, represented as . This context involves neutrons within a specific angular dEn dΩ acceptance and energy range. We adjust the definition of absolute efficiency, ϵabs , transforming it into a function of the neutron energy En , and angles θ and ϕ: , d2 Mmeas d2 Memit ϵabs (En , θ, ϕ) ≡ , (5.9) dEn dΩ dEn dΩ with both differential multiplicities evaluated at the same En , θ, ϕ. The denominator is the differential multiplicity of neutrons emitted from the source that we would like to reconstruct, 128 whereas the numerator is the differential multiplicity of neutrons that are measured by the detector. We calculate and apply all efficiencies within this study in the laboratory frame. Figure 5.10 The effective thickness of NWB as a function of θ and ϕ. The effective thickness is defined as the accumulated path lengths of a neutron through the active volumes of NWB. The black ellipses are selected regions for Monte Carlo simulations to generate Figure 5.11. The color bar is truncated at 6.35 cm, which is the longitudinal thickness of the NWB bars. Effective thicknesses below this value are represented by the same color. Directly using this revised definition of ϵabs requires performing Monte Carlo simulations with a large number of neutrons at different (En , θ, ϕ). Creating a smooth efficiency surface for ϵabs would involve repeating this process for numerous (En , θ, ϕ) points, each necessitating the simulation of a substantial number of neutrons. This approach is both computationally intensive and time-consuming. Furthermore, given the NWB’s limited position resolution, excessive refinement of the efficiency surface in (θ, ϕ) yields diminishing returns. Thus, it is more practical to quantify the absolute efficiency as separate functions of En and (θ, ϕ): ϵabs (En , θ, ϕ) ≈ ϵint (En ) · ϵgeo (θ, ϕ) , (5.10) 129 where ϵint (En ) denotes the intrinsic efficiency, and ϵgeo (θ, ϕ) refers to the geometric efficiency. 0.14 0.12 0.10 Efficiency 0.08 0.06 0.04 0.02 0.00 6.35 6.40 6.45 6.50 6.55 6.60 Effective thickness [cm] Figure 5.11 Neutron efficiencies as a function of thickness from selected regions on Figure 5.10. The black dashed line indicate the average efficiency, at 12.2%. The error bars are statistical uncertainties. We justify this approximation by examining the effective thickness distribution of the NWB as a function of θ and ϕ, as shown in Figure 5.10. The effective thickness is defined as the accumulated path lengths of a neutron through the active volumes of NWB. It is the primary factor that influences the probability of a neutron generating a detectable light signal upon passing through the active volume. However, due to the fact that NWB is placed relatively far from the target, the most substantial incident angle at the corners of the NWB is only ∼ 19◦ , increasing the effective thickness by merely ∼ 5%. To further corroborate this approximation, we select a series of ellipses on Figure 5.10 to cover a range of effective thicknesses. We run a Monte Carlo simulation for each elliptical region using a newly developed neutron simulation code, “neuSIM4” [188], which we will explore in more depth in subsection 5.2.3. As Figure 5.11 illustrates, the absolute efficiency exhibits practically no sensitivity to thickness variation. It is a reasonable approximation to assume that the absolute efficiency depends solely on neutron incident energies and is independent of any θ and ϕ within the NWB’s coverage. The simulations are repeated at 130 neutron energies of En = 20, 40, 60, 80 MeV, all of which exhibit no significant thickness dependence within statistical uncertainty. We thus consider ϵabs (En , θ, ϕ) as a product of ϵint (En ) and ϵgeo (θ, ϕ), which are independent of each other. In the next two subsections, the determination of these two efficiencies will be discussed in detail. 5.2.2 Geometric efficiency The geometric efficiency of a neutron detector is determined by the solid angle it subtends, making it independent of the neutron energy En . Given the azimuthal symmetry assumption d2 Memit of neutron emission, the energy spectrum shall not depend on ϕ. Hence, it is dEn dΩ sufficient to consider a “ring” of solid angle dΩθ at a given θ. This modifies the definition of geometric efficiency to ∆ϕ(θ) ϵgeo (θ) = , (5.11) 2π where Z ∆ϕ(θ) ≡ dϕ . (5.12) dΩθ ∩D Here, the integral is computed over the θ-slice of the detector’s coverage, D. It is important to note that the denominator is 2π radians, not 4π steradians, as we are considering a ring of circumference. To determine ∆ϕ(θ), Monte Carlo simulations are often employed to accommodate detectors with diverse geometries [189, 190, 191]. The hit pattern on NWB, shown in Figure 5.12, is produced by a Monte Carlo simulation. This simulation designates only the NE-213 filled internal cavity of the bar as the active volume, excluding the Pyrex shell. It involves the emission of one billion isotropic rays from the center of the target, used to identify intersection points with the detector. Although some rays can intersect multiple bars, the simulation confirms that a single ray can intersect a maximum of two NWB bars. Among all rays intersecting the active volumes, a mere 0.808% ± 0.001% intersect two bars, with the remainder intersecting exactly one bar. To calculate the geometric efficiency of NWB, we evaluate ∆ϕ at various θ values. The 131 Figure 5.12 Monte Carlo simulation of hits on NWB. Blue and red dashed lines depict the center lines of odd and even bars, respectively. The gray dashed curves that go roughly vertically indicate the bar positions xbar . first method involves Monte Carlo simulations. For a fixed θ, we generate N rays with ϕ ∈ [0, 2π) and tally the number that intersect the detector, Nintersected . This leads to Nintersected ∆ϕ ≈ . (5.13) N 2π Despite the straightforward nature of Monte Carlo simulations in computing ∆ϕ, they can be computationally intensive. A more efficient method is to determine the total ϕ coverage of NWB at a specific θ, termed the “δϕ-method”. Figure 5.13 demonstrates the ϕ coverage of NWB at θ = 30◦ . A comparison with Figure 5.12 reveals a finite number of continuous segments in ϕ along the slice at θ = 30◦ , separated by small gaps. Simple root-finding algorithms like binary search can accurately and efficiently determine these numerical boundaries. Given the boundaries of these segments as, (low) (upp) (low) (upp) [ϕ1 , ϕ1 ], [ϕ2 , ϕ2 ], . . . , [ϕ(low) n , ϕ(upp) n ], (5.14) 132 1.0 Intersected 0.5 0.0 30 20 10 0 10 20 30 Lab φ [deg] Figure 5.13 The ϕ coverage of NWB at θ = 30◦ . The vertical axis represents the intersection state, with 1 indicating intersection and 0 indicating no intersection. (upp) (low) where ϕi < ϕi+1 for all i = 1, 2, . . . , n − 1 and each segment having a width of δϕi ≡ (upp) (low) ϕi − ϕi , the ∆ϕ can then be calculated as Z n X ∆ϕ = dϕ = δϕi . (5.15) dΩθ ∩D i=1 The δϕ-method permits a precision up to the floating-point tolerance, and its computational cost is independent of the number of rays, making it considerably more efficient than Monte Carlo simulations. 0.14 Monte Carlo 0.12 δφ-method 0.10 ∆φ 0.08 2π 0.06 0.04 0.02 0.0025 30 35 40 45 50 55 25 30 35 40 45 50 55 Lab θ [deg] Lab θ [deg] Figure 5.14 Comparison of geometric efficiencies of NWB calculated using Monte Carlo simulations and the δϕ-method. The left panel displays results without shadow bars, and the right panel includes shadow bars. NWB bars 1 to 24, inclusively, are considered in the calculations. Figure 5.14 juxtaposes geometric efficiencies calculated using Monte Carlo simulations 133 and the δϕ-method. These curves are determined from θ = 25◦ to θ = 55◦ , inclusively, in increments of 0.1◦ . For the Monte Carlo method, each increment generates 105 random rays. The efficiency of the δϕ-method facilitates the inclusion of additional geometric cuts into the NWB during the analysis with minimal computational cost. Both panels of Figure 5.14 truncate all NWB bars at xbar = ±90 cm. This edge cut, defined empirically and implemented in the experimental data, excludes hits too close to the NWB bar edges where data quality is compromised. This cut reduces the length by 6.52 cm at each end compared to the full length that sets the two ends at xbar = ±96.52 cm. Moreover, shadow bar cuts can be introduced by removing hits at xbar ∈ [−50 cm, −15 cm] or xbar ∈ [15 cm, 50 cm] in NWB bars 07, 08, 09, 15, 16, and 17. For certain runs where some NWB bars may malfunction, it is necessary to exclude these bars from the analysis and re-calculate the geometric efficiency. Figure 5.15 Neutron energy spectrum of the 48 Ca +64 Ni reaction system from run 4082 to 4116, inclusively, in which shadow bars were present. The left and right panels show the spectrum before and after geometric efficiency correction, respectively. To check the resulting geometric efficiency, we apply it to raw neutron spectrum data. To contrast the effect, we have chosen runs with shadow bars, which are runs 4082 to 4116, inclusively. Sharp position cuts at xbar = ±90 cm as well as shadow bar cuts at 134 xbar ∈ [−50 cm, −15 cm] or xbar ∈ [15 cm, 50 cm] in NWB bars 07, 08, 09, 15, 16, and 17 are applied to the raw neutron spectrum data. The resulting spectra are shown in Figure 5.15. After the geometric efficiency correction, the neutron spectrum is proportionally scaled up, and the strips due to shadow bars at around θ ≈ 35 ± 2◦ and θ ≈ 44 ± 2◦ are no longer visible. The neutron energy spectrum is now a smooth continuum along θ, except statistical fluctuations. 5.2.3 Intrinsic efficiency Intrinsic efficiency, denoted as ϵint , is defined as the probability of a neutron traversing the detector’s active volume, subsequently generating a detectable light signal that is identified as a neutron. As discussed in subsection 5.2.1, we may postulate that ϵint remains independent of the variables θ and ϕ. Consequently, this allows for an expression that solely depends on neutron energy En :  , Z d2 Mmeas d2 Memit Z  ϵint (En ) ≡ dΩ dΩ . (5.16) D dEn dΩ D dEn dΩ Here, both integrals are taken over the entire detector’s projection onto the (θ, ϕ) space, denoted as D. The numerator and denominator are then functions of energy to be evaluated at the same En . In the course of actual computations, this ratio is approximated by conducting a series of Monte Carlo simulations for a substantial number of neutrons at various neutron energies. We emit isotropic neutrons, and only those that intersect the detector are counted toward the denominator. Nonetheless, intersection with the detector does not assure that the neutron will yield a detectable light signal. It is only upon generating a detectable light signal that surpasses a predetermined threshold, referred to as the bias, that it is added to the numerator. To simulate the neutrons, we adopt a simulation code, “neuSIM4”, recently developed by J. Park et al. [188]. This code integrates neutron detection physics from SCINFUL- QMD [192] into Geant4 [168, 169], a versatile simulation toolkit used in nuclear and particle physics. Unlike SCINFUL-QMD, which is limited to a fixed cylindrical detector with a 135 single photomultiplier tube (PMT), neuSIM4 can handle detectors of various geometries and multiple PMTs. The major challenge faced when developing a robust simulation code for neutron detection is to accurately model how neutrons interact with the scintillator molecules and the subsequent scintillation process. Unlike the detection of charged particles, neutron detection involves complex processes. Neutrons, being electrically neutral, interact with atomic nuclei via elastic and inelastic scattering rather than directly with the electrons of the scintillator atoms. 101 1 H(n, n) 1 H 12 C(n, n) 12 C 12 C(n, nγ) 12 C 100 Cross section σ [barn] 12 C(n, 2n) 11 C 12 C(n, p) 12 B 12 C(n, d) 11 B 10 1 12 C(n, np) 11 B 12 C(n, t) 10 B 12 C(n, 3 He) 10 Be 10 2 12 C(n, α) 9 Be 12 C(n, n 0 3α) n + 12 C total 10 3 5 10 50 100 Neutron energy En [MeV] Figure 5.16 Cross sections as functions neutron incident energies En . The gray vertical dashed line indicates the En = 110 MeV, above which the cross section data for many channels are not available and constant interpolation is used. Reproduced from [192, 188]. NW bars are filled with NE-213, a liquid scintillation material based on solvent xylene (C8 H10 ) [170, 171]. When a neutron enters the NW bar, it interacts with the hydrogen or carbon nuclei of the xylene molecules. One of the keys to accurately simulate the neutron detection process is to consider as many reaction channels as possible with accurate cross sections. A total of 11 initial reaction channels originally employed in SCINFUL are considered in neuSIM4, and they are listed in the legend of Figure 5.16. For neutron energies below 136 110 MeV, the cross section data have carefully been compiled during the development of SCINFUL [193]. Despite the most dominant initial channels being the elastic scatterings of 12 1 H(n, n) 1H and C(n, n) 12C, the aggregate cross section of all other inelastic channels is not negligible. We include these cross section data in neuSIM4. Aside from the initial channels, neuSIM4 considers the sequential decays. The code determines the excitation energy of the residual fragments, deduced from the energy difference between the initial and final states of the reaction. If the residual nucleus possesses adequate excitation energy, it either decays into the daughter nuclei or achieves its ground state by emitting gamma rays. During each step, light output data are stored, and the probabilities of all sequential decays are based on the SCINFUL database. Above 110 MeV, the cross section data are not available for all channels. Hence, we switch to Liège intranuclear cascade (INCL) as the simulation model [194, 195]. Switching of simulation models will inevitably introduce discontinuities in the simulation results. To mitigate this effect, we introduce a smooth transition using a continuous weighting function that samples the simulation results from both models. Specifically, we have      1 if En ≤ 80 MeV En − 200  PSCINFUL (En ) = if 80 MeV < En < 200 MeV , (5.17)    80 − 200 if En ≥ 200 MeV   0 and PINCL (En ) = 1 − PSCINFUL (En ) . (5.18) For a given neutron energy En , the probability of using SCINFUL is PSCINFUL (En ), and the probability of using INCL is PINCL (En ). The switching point at En = 80 MeV is chosen because SCINFUL was established to describe the experimental data up to this energy [193]. The switching point at En = 200 MeV is chosen empirically as INCL starts to yield more accurate results than SCINFUL above this energy when compared to the experimental data. This smoothing procedure requires SCINFUL to be run above 80 MeV despite the lack of 137 cross section data for many channels. For that, we simply use constant extrapolation for all channels that do not have cross section data. particle a1 a2 a3 minimum energy [MeV] p 0.81 2.43 0.29 25 d 0.74 3.45 0.20 10 t 0.72 7.19 0.07 10 3 He 0.54 3.97 0.20 5 4 He 0.51 6.42 0.08 5 Table 5.1 Coefficients in Equation 5.19 for various particles. Adopted from [196]. Having obtained the fragments from the nuclear reactions, neuSIM4 then calculates the light output contribution from each fragment. For a fragment with energy above the corresponding minimum energy listed in Table 5.1, the light output is calculated using L(En ) = a1 En − a2 [1 − exp(−a3 En )] . (5.19) For a fragment with energy below the minimum energy, the light output is given by the experimental work by V.V. Verbinski et al. [197]. With the light output from all fragments, neuSIM4 can now calculate the total light output from a neutron. In Figure 5.17, we show the light responses at various energies up to 200 MeV. A light response is defined as the probability distribution of light output from the scintillator material, given a neutron of a specific energy. Unlike charged particles nor gamma rays, neutrons do not directly deposit their energy in the scintillator material. Instead, they transfer their energy to the recoiling nuclei through various channels. As a result, the light response of a neutron do not exhibit a sharp peak, as in the spectroscopy of charged particles or gamma rays, but instead present a broad, smooth continuum of light output. The unit for light response is MeVee−1 . Suppose we denote the light response curve as f (L), then integrating f (L) over L from zero to infinity yields the total probability of a neutron generating non-zero light output, i.e., the intrinsic efficiency. In the actual world, a detector cannot detect arbitrarily small light output due to the presence of noise floor in the electronics. Consequently, researchers need to introduce a bias, Lbias . The bias denotes 138 109 200 MeV (×10 12) 106 190 MeV (×10 11) 180 MeV (×10 10) 170 MeV (×10 9) 103 160 MeV (×10 8) Light response [MeVee−1 ] 150 MeV (×10 7) 140 MeV (×10 6) 100 130 MeV (×10 5) 120 MeV (×10 4) 110 MeV (×10 3) 10 3 100 MeV (×10 2) 90 MeV 80 MeV (×10(0×) 10 ) 1 70 MeV 65 MeV (×(10×−101.5) ) −1 10 6 60 MeV (×10 −2 50 MeV (×10 −3) ) 40 MeV (×10 −4) 30 MeV (×10 −5) 10 9 0 20 40 60 80 100 120 140 Light output [MeVee] Figure 5.17 Neutron light responses at various energies. The solid lines are simulated results using neuSIM4, and the points are experimental data collected by NWB. Statistical uncertainties are smaller than the point size. Colors are alternating between blue and red for easy distinction among neighboring curves. 139 the smallest light output that can be accurately measured by the detector. This results in a revised definition for intrinsic efficiency: Z ∞ ϵint ≡ dL f (L) . (5.20) Lbias The bias is chosen such that it is small enough to capture the majority of the light output, yet large enough to avoid the noise floor. In our study, we set the bias to Lbias = 3.0 MeVee based on empirical observations. Normalizing light response curves from simulation is a straightforward process because the total probability can be calculated directly. The same process, however, becomes challenging with our experimental data. This challenge arises from the fact that the total probability or the intrinsic efficiency is unknown in our experiment. Unlike experiments that were specifically designed to measure the neutron efficiency [198, 199, 200], which were often performed with a monoenergetic neutron source with a well-defined flux, our experiment measures neutrons emitted from heavy-ion collisions with a broad energy spectrum and undetermined cross sections. Thus, in instances where we must compare the light response curves from simulation and experimental data, a normalization process must be applied to the experimental data. We have illustrated this in Figure 5.17, in which we scale the experimental light response curves such that the integral from bias to infinity is equal to the same integral from the simulation. Figure 5.18 consolidates the light responses at various energies into a two-dimensional surface for comparative analysis. For simplicity, the color scale represents histogram counts without any normalization. The left panel is the experimental data gathered by NWB, and the right panel displays results simulated using neuSIM4. We can observe at least two features in these light response surfaces. The first is the high-energy edge of the light response, which signals the maximum light output that a neutron can generate at a given incident energy. The second is a faint peak to the left of the edge, with a similar slope to the edge. For example, in the light response of En = 65 MeV, a small peak around 35 MeVee is visible. 12 The simulation suggests this peak is attributable to the C(n, d) 11B reaction channel. This 140 Figure 5.18 Neutron light response surfaces. The left panel is from experimental data collected by NWB, and the right panel is from neuSIM4 simulations. (n, d) peak is also observed in several experimental light responses [199], including our data collected by NWB, yet it is often missing from simulations [201, 152], except those that utilize SCINFUL’s cross section database [196]. Briefly put, the resemblance between the two surfaces strongly suggests that relevant physics are most likely included in neuSIM4 which can accurately simulate the neutron detection process in NWB. The one-dimensional light response curves in Figure 5.17 are essentially horizontal slices of this two-dimensional surface, appropriately normalized. Note that comparing energy distributions between these two surfaces is not meaningful because the energy distributions in the simulation are set uniformly, while those in the experimental data depend on both the collision kinematics and the detector response. Lastly, we present in Figure 5.19 the intrinsic efficiency curve for the NWB as a function of neutron energy En with a bias of 3.0 MeVee. At lower energies, below 20 MeV, the intrinsic efficiency ϵint is around 10%. Beyond this, it gradually declines and flattens out at around ∼ 3% as En approaches 300 MeV. 141 0.10 0.08 Intrinsic efficiency ²int 0.06 0.04 SCINFUL 0.02 SCINFUL + INCL INCL 0.000 50 100 150 200 250 300 Neutron energy En [MeV] Figure 5.19 The intrinsic efficiency curve for the NWB as a function of neutron energy En with a bias of 3.0 MeVee. 5.3 Neutron spectra We are now ready to apply all the calibrations, corrections, and efficiencies to the experimental data and reconstruct the neutron spectra. To construct the energy spectrum of neutrons in the lab frame at a given θ′ , we first select all events or hits that satisfy the following criteria: • Select hits on xbar ∈ [−90 cm, 90 cm]. This is the edge cut that removes hits too close to the NWB bar edges where data quality is compromised. • Discard shadow bar hits, if applicable by eliminating hits on xbar ∈ [−50 cm, −15 cm] or xbar ∈ [15 cm, 50 cm] in NWB bars 07, 08, 09, 10, 15, 16, and 17. • Veto charged particles by requiring the absence of hits in the Veto wall (VW). • Discriminate neutrons from gamma rays by requiring PPSD > 0.5. • Select hits with light output of L ≥ 3.0 MeVee. 142 • Select central events with Microball multiplicity greater than or equal to the value that corresponds to impact parameter b ≤ 3 fm according to subsection 2.1.2. This allows us to construct the “raw” neutron energy spectrum as a histogram, represented by N (i, j), where i is the bin index of the neutron energy En and j is the bin index of θ. Next, we apply the background subtraction and geometric efficiency correction to the raw spectrum according to d2 M 1 − ηbg (En′ ) N (i′ , j ′ )   1 X 1 1 (En′ ) ≈ · · · · . (5.21) dEn dΩ Nevent ϵint (En′ ) 2π sin θ′ ϵgeo (θ′ ) ∆En · ∆θ |θj′ −θ′ |≤1◦ Here, Nevent is the total number of central events, ηbg is the background fraction, ϵint is the intrinsic efficiency, ϵgeo is the geometric efficiency, ∆En is the width of the energy bin, and ∆θ is the width of the angular bin. The summation is over all bins that are within 1◦ of θ′ . 40 58 In Figure 5.20, we show the neutron energy spectra from Ca + Ni at 140 MeV/u corrected by ϵint and ϵgeo step by step. Starting from the bottom group, the spectra are normalized to the number of central collision events with neutron background subtracted. The middle group are spectra from the bottom group corrected by the geometric efficiency. Because the geometric efficiency is only dependent on the lab θ angle and is independent of neutron energy En , the spectra in the middle group share the same shape as those in the bottom group. The only changes are the vertical scales. For spectrum at θlab = 30◦ ± 1◦ , the average geometric efficiency over the range of θlab ∈ [29◦ , 31◦ ] is ϵgeo ≈ 12.5% (see Figure 5.14), hence the vertical scale of the middle group is about 1/0.125 = 8 times that of the bottom group. For spectrum at θlab = 45◦ ± 1◦ , the average geometric efficiency is about 10%, giving a 1/0.1 = 10 times enhancement. As a result, applying geometric efficiency correction squeezes the spectra from different θ slices closer together. In practice, the geometric efficiency correction is applied to the finer angular bins of 0.1◦ by multiplying the spectra by a factor of 1/ϵgeo (θlab ) before rebinning to the final angular bins of 1◦ . Lastly, we have the top group, which are spectra from the middle group corrected by the intrinsic efficiency. This correction contributes another order of magnitude enhancement to the spectra. Moreover, it alters the 143 10 2 × 1 ²int (En ) 10 3 d 2 M [(MeV · sr) −1 ] × 1 ²geo (θlab ) dEn dΩ 10 4 10 5 θlab = 30 ◦ ± 1 ◦ θlab = 35 ◦ ± 1 ◦ θlab = 40 ◦ ± 1 ◦ θlab = 45 ◦ ± 1 ◦ 10 6 0 50 100 150 200 250 300 Lab En [MeV] Figure 5.20 Lab energy spectra of neutrons from 40Ca + 58Ni at 140 MeV/u corrected step by step to the final spectra. The bottom group are spectra that are normalized to the number of central collision events and subtracted by the background. The middle group are spectra from the bottom group corrected by the geometric efficiency. The top group are spectra from the middle group corrected by the intrinsic efficiency. 144 shapes of the spectra, as the intrinsic efficiency ϵint (En ) is energy-dependent. We see the relative spectra at lower energies (En ≲ 80 MeV) become flatter due to the fact that the intrinsic efficiency is higher at lower energies, as shown in Figure 5.19. 40 Ca + 58 Ni at 140 MeV/u 10 2 d 2M dEn dΩ [(MeV · sr) ] −1 10 3 θlab = 30 ◦ ± 1 ◦ θlab = 35 ◦ ± 1 ◦ θlab = 40 ◦ ± 1 ◦ θlab = 45 ◦ ± 1 ◦ 10 4 0 50 100 150 200 250 300 Lab En [MeV] 48 Ca + 64 Ni at 140 MeV/u 10 2 d 2M dEn dΩ [(MeV · sr) ] −1 10 3 θlab = 30 ◦ ± 1 ◦ θlab = 35 ◦ ± 1 ◦ θlab = 40 ◦ ± 1 ◦ θlab = 45 ◦ ± 1 ◦ 10 4 0 50 100 150 200 250 300 Lab En [MeV] Figure 5.21 Final lab energy spectra of neutrons and protons. The top panel is from the 40Ca + 58Ni reaction system, and the bottom panel is from the 48Ca + 64Ni reaction system. Beam energies are at 140 MeV/u for both systems. The spectra in the final group of Figure 5.20 are the final lab energy spectra of neutrons. 145 40 To better highlight the results, we present them in Figure 5.21, with Ca + 58Ni reaction 48 in the top panel and Ca + 64Ni reaction in the bottom panel. The beam energies for both systems are at 140 MeV/u. In these figures, we have constructed the neutron spectra with neutrons with energies from 20 MeV to 300 MeV. Neutrons below 20 MeV are discarded due to detector noise and excessive background. Unlike charged particles detection using HiRA10 in which the energy is determined from the energy loss in the scintillator, the energy of neutrons is determined from the time-of-flight (ToF). Hence, there is not really an upper bound on the energy of neutrons due to punch-through. However, we are limited by the time resolution of the detectors. For that, we have chosen to present the spectra up to 300 MeV. When constructing quantities such as isoscaling and neutron-to-proton ratios, we will only use neutrons with energies below 200 MeV mainly because protons above that energy have been discarded due to the punch-through cut. At 200 MeV, the neutron energy has a resolution of FWHM ≈ 5 MeV, given by Equation 2.8. 5.3.1 Transverse momentum spectra While the lab energy spectra of neutrons are straightforward to construct, they do not offer an easy way to distinguish nucleons as participants and spectators. The interaction dynamics of nucleons in the participant region are expected to behave differently from those in the spectator region. For participant nucleons, they actively participate in the collision, allowing them to encode information about the collision dynamics for probing nuclear EOS at high density regime. The nucleons would experience squeeze-out motion, resembling the collective flow of a classical fluid [133, 133, 134]. This results in a significant emission in the transverse direction, motivating us to study the transverse momentum spectra. The transverse momentum is defined as pT = p sin θ , (5.22) where p is the magnitude of the momentum vector. To give a complete kinematics phase 146 space like the θ-E phase space, we introduce the rapidity y as   1 E + p cos θ y = ln , (5.23) 2 E − p cos θ p where E = p2 + m2 is the total energy of the particle with mass m. To facilitate the comparison between spectra from reaction systems with different beam energies, we shall normalize the rapidity by the beam rapidity ybeam in the lab frame, which is denoted as y ŷ ≡ . (5.24) ybeam To select neutrons from the participant region, we often apply a “mid-rapidity” cut by selecting neutrons with ŷ ∈ [0.4, 0.6]. Figure 5.22 The two-dimensional pT /A versus ŷ spectra of neutrons. In the left panel, the blue dashed lines the lab energy curves. In the right panel, the lines indicate the phase space boundaries for proton, deuteron, and triton. In Figure 5.22, we present the two-dimensional pT /A versus ŷ neutron spectrum. The shaded vertical region indicates the mid-rapidity cut of ŷ ∈ [0.4, 0.6]. The neutron spectrum is shown as a color map and is identical in both left and right panels. Only the lines and annotations are added differently to highlight different aspects of the spectrum. In the left panel, the blue dashed lines are the energy curves, q En = p2T + m2n · cosh y , (5.25) 147 for En = 10, 30, 65, 140, 200 MeV. As can been see from the intersections of the energy curves and the neutron color map, only neutrons with energies from 30 MeV to 140 MeV are selected by the mid-rapidity cut. In the right panel, we overlay the phase space boundaries for protons, deuterons, and tritons detected by HiRA10. All three charged particles have the same θ coverage from 30◦ to 75◦ , whereas for neutrons, the θ coverage is from 28◦ to 51.5◦ . The red dashed curves, the orange dash-dotted curves, and the brown dotted curves represent the punch-through boundaries and energy thresholds for protons, deuterons, and tritons, respectively. In short, vLANA has a narrower coverage than HiRA10 in θ but a wider coverage in energy. 40 Ca + 58 Ni at 140 MeV/u 48 Ca + 64 Ni at 140 MeV/u 10 1 d 2M d(pT /A) dŷ [(MeV/c) ] −1 10 2 pseudo-neutron neutron 10 100 3 200 300 400 500 100 200 300 400 500 pT /A [MeV/c] pT /A [MeV/c] Figure 5.23 The left and right panels show the pT /A spectra neutrons and pseudo- neutrons from 40Ca + 58Ni and 48Ca + 64Ni reactions at 140 MeV/u, respectively. 40 In Figure 5.23, we present the one-dimensional pT /A neutron spectra from Ca + 58Ni and 48 Ca + 64Ni reactions at 140 MeV/u, along with the pseudo-neutron spectra for comparison. These spectra are constructed by integrating the two-dimensional spectra in Figure 5.22 over ŷ ∈ [0.4, 0.6]. To account for the yield loss due to the incomplete phase space in pT -ŷ, we apply a coverage correction to the transverse momentum spectrum by dividing the spectrum by the coverage fraction within the mid-rapidity slice at each pT bin, as illustrated in Figure 3.3. Given the limited θ coverage of NWB, as shown in Figure 5.22, the pT /A range of neutron spectra is confined up to 400 MeV/c, whereas the pseudo-neutron spectra can extend slightly 148 40 Ca + 58 Ni at 140 MeV/u 48 Ca + 64 Ni at 140 MeV/u 100 d 2M d(pT /A) dŷ [(MeV/c) ] 10 1 −1 10 2 CI n w/ pseudo-neutron CI n w/ actual neutron 10 100 3 150 200 250 300 350 400 100 150 200 250 300 350 400 pT /A [MeV/c] pT /A [MeV/c] Figure 5.24 The left and right panels show the coalescence-invariant (CI) neutron pT /A spectra constructed using pseudo-neutrons and actual neutrons for 40Ca + 58 Ni and 48Ca + 64Ni reactions at 140 MeV/u, respectively. further up to 450 MeV/c due to the broader θ coverage of HiRA10. In Figure 5.24, we also present the coalescence-invariant (CI) neutron spectra constructed using pseudo-neutrons and actual neutrons. 5.3.2 Spectral yield ratios Having constructed the transverse momentum spectra for both Ca + Ni systems, we may update the isoscaling analysis and the Bayesian analysis of charged particles in section 3.2 and section 3.4, respectively, with the new neutron data points from vLANA. This requires us to construct a series of spectral yield ratios involving neutrons. First, we construct the isoscaling ratio spectra for neutrons, Yn ( 48Ca + 64Ni) R21 (n) = . (5.26) Yn ( 40Ca + 58Ni) In the process of calculating R21 (n), the effect of certain corrections will be canceled out due to their independence from the specific experimental run. We shall refer to these corrections as run-independent because they are assumed to share the same exact functional form for all experimental runs. In this work, intrinsic efficiency ϵint (En ) is the only run-independent correction, which depends only on the neutron energy En . Geometric efficiency ϵgeo (θlab ) is subject to variations due to reasons including but not limited to the presence or absence 149 of shadow bars and the removal of bars with poor data quality in the analysis. Similarly, backgrounds are also run-dependent, influenced by a range of factors such as beam-target configurations, beam intensity, trigger conditions, and changes in the experimental setup surroundings. All calibrations are also considered run-dependent due to potential drifts in electronics or detector responses over time. Thus, to ensure the accuracy of the isoscaling ratio determination, we carefully considered and applied these run-dependent corrections and calibrations. 1.8 α = 0.31 β = − 0.24 1.5 C = 1.07 0 Z= R21 (N, Z) 1.0 = 2 1 Z Z= 0.8 0 1 2 N Figure 5.25 The data points are experimental data, and the dashed lines are the isoscaling fits, with fit parameters listed on the top left corner. In Figure 5.25, we present the isoscaling fits on neutrons as well as the previously analyzed charged particles, including p, d, t, 3He, and 4He. The data points are isoscaling ratios from yields between 120 MeV/c and 350 MeV/c in pT /A. The dashed lines are the isoscaling fits, with fit parameters given to be α = 0.31, β = −0.24, and C = 1.07. It can be seen that the fit performs well for all six particles. Similarly to the charged particles, we can also construct the isoscaling ratio including free neutrons as a function of pT /A, In Figure 5.26, we have added the isoscaling ratio spectrum for neutrons as green crosses. The dashed lines are the isoscaling fits to the data points at 150 1.8 t n 1.6 d 4 He p 3 He 1.4 N−Z=1 R21 1.2 N−Z=0 1.0 0.8 N−Z= −1 0.6 0.5 α ∆µ/T 0.0 β 0.5 12 pT /A [MeV/c] 11 TH − He [MeV] 10 9 40 Ca + 58 Ni 8 48 Ca + 64 Ni 70 100 200 300 400 500 pT /A [MeV/c] Figure 5.26 Top: The isoscaling ratio spectra for the Ca + Ni systems at 140 MeV/u. The points are experimental measurements, whereas the dashed lines are the isoscal- ing fits. Middle: The fitted α and β values. Bottom: The TH−He temperatures for both systems. each pT /A bin. The isoscaling ratio of neutrons fall between the isoscaling ratios of triton and helium-4, extending grouping property by the neutron excess number, N − Z, that we observed with the light charged particles: p and 3He belong to the N − Z = −1 group, d and 4He belong to the N − Z = 0 group, t and the newly added n belong to the N − Z = 1 group. The fitted parameters α and β exhibit similar values and trends as those analyzed only with charged particles. Lastly, the temperatures TH−He are identical to previous values because the Albergo thermometer is a measurement that uses only charged particles. 151 48 Ca + 64 Ni 40 Ca + 58 Ni free n free n 1.8 R21 (n) 1.6 10 1 10 1 1.4 1.2 1.0 10 2 10 2 0.8 0.6 free p free p 1.8 R21 (p) 1.6 10 1 10 1 1.4 1.2 1.0 10 2 10 2 0.8 0.6 1.8 free n/p 1.8 free n/p 1.8 free DR 1.6 1.6 1.6 1.4 1.4 1.4 1.2 1.2 1.2 1.0 1.0 1.0 0.8 0.8 0.8 0.6200 250 300 350 400 0.6200 250 300 350 400 0.6200 250 300 350 400 pT /A [MeV/c] pT /A [MeV/c] pT /A [MeV/c] Figure 5.27 The top-left 2 × 2 panels, enclosed by the orange dashed line, are the spectra of free neutrons and protons. The bottom-right panel is the free double ratio. The rest of the panels are free single ratios. 152 Having re-established the isoscaling property with neutrons, we proceed to update the Bayesian analysis in section 3.4 that was performed with charged particles and pseudo- neutrons. We begin by compiling the spectra and ratios of free neutrons and protons in Figure 5.27. The word “free” is emphasized to distinguish them from the coalescence-invariant (CI) versions that will be presented later. The data points in each panel are the experimental measurements, whereas the gray solid lines beneath them are calculations from the 75 Latin hypercube sampling (LHS) of (S0 , L, m∗s , m∗v ) used to train the Gaussian process emulator. This 3 × 3 grid offers a convenient way to inspect all the spectra and ratios at once. The 48 first and second columns represent observables from Ca + 64Ni and 40 Ca + 58Ni reactions, respectively. Each row in the third column can be obtained by dividing the first column by the second column, as hinted by the left-to-right purple arrow at the bottom of the figure; the third column comprises, from top to bottom, the isoscaling ratio of neutrons, the isoscaling ratio of protons, and the double ratio of n/p. The grid can also be interpreted vertically from top to bottom. Each column in the bottom row can be obtained by dividing the top row by the middle row, as indicated by the top-to-bottom purple arrow on the right side of the 48 figure; the bottom row comprises, from left to right, the single n/p ratio of Ca + 64Ni, the 40 single n/p ratio of Ca + 58Ni, and the double ratio of n/p. Notice that the double ratio in bottom-right panel can be obtained in two equivalent ways, R21 (n) SRn/p ( 48Ca + 64Ni) DR = = . (5.27) R21 (p) SRn/p ( 40Ca + 58Ni) As discussed in subsection 3.3.1, one of the major challenges is the inaccuracy of the coalescence model used by many QMD models to form nucleons into clusters (isotopes). In the top-left 2 × 2 panels, the discrepancies between the experimental data and the calculations are evident for the free neutron and proton spectra. Specifically at lower pT /A (≲ 350 MeV/c), the calculations overestimate the spectra up to a factor of 5. The calculations also fail to reproduce the overall shapes of the spectra for both reaction systems. To circumvent this issue, the clusters are decomposed into their constituent nucleons, and these nucleons are used to populate the so-called coalescence-invariant (CI) spectra [90, 89], 153 48 Ca + 64 Ni 40 Ca + 58 Ni CI n CI n 1.8 R21 (CI-n) 1.6 10 1 10 1 1.4 1.2 1.0 10 2 10 2 0.8 0.6 CI p CI p 1.8 R21 (CI-p) 1.6 10 1 10 1 1.4 1.2 1.0 10 2 10 2 0.8 0.6 1.8 CI n/p 1.8 CI n/p 1.8 CI DR 1.6 1.6 1.6 1.4 1.4 1.4 1.2 1.2 1.2 1.0 1.0 1.0 0.8 0.8 0.8 0.6200 250 300 350 400 0.6200 250 300 350 400 0.6200 250 300 350 400 pT /A [MeV/c] pT /A [MeV/c] pT /A [MeV/c] Figure 5.28 The top-left 2 × 2 panels, enclosed by the orange dashed line, are the CI spectra of neutrons and protons. The bottom-right panel is the CI double ratio. The rest of the panels are single ratios. 154 minimizing the systematic errors introduced by the coalescence model. In Figure 5.28, we present the CI n, CI p spectra, as well as all the relevant CI ratios. Starting from the top-left 2 × 2 panels, we have the CI n and CI p spectra now compare better between the experimental data and the calculations than the free spectra in Figure 5.27. Nonetheless, the calculations still tend to overestimate the spectra at lower pT /A (≲ 350 MeV/c). Given that the calculations (gray lines) fail to fully enclose the experimental data points, the Bayesian analysis will give nonphysical out-of-range predictions for the posterior distributions of (S0 , L, m∗s , m∗v ). For the CI isoscaling ratios R21 and CI single ratios of n/p, the calculations are able to compare better with the experimental data. However, the shapes of these ratios are not well reproduced, hence we are not able to extract meaningful Bayesian results from these observables. We note that experimentally, the isoscaling ratios would minimize the systematic errors due to detector responses and efficiencies, whereas the single ratios would minimize the systematic errors due to the merging of different runs. As for the theoretical calculations, the isoscaling ratios would minimize biases in the coalescence model and other inaccuracies in ImQMD, whereas the single ratios would minimize the any reaction-dependent systematics. Together, we construct the CI double ratio of n/p shown in the bottom-right panel, in which the calculations are able to reproduce both the shape and the magnitude of the experimental double ratio, facilitating the use of Bayesian inference. parameter mean std. dev. S0 [MeV] 31.3 1.5 L [MeV] 87 15 m∗s /mN 0.75 0.09 m∗v /mN 0.80 0.07 Table 5.2 Means and standard deviations of the posterior distributions of S0 , L, m∗s , m∗v from the CI-n/p double ratio using free neutron spectra. Figure 5.29 presents posterior distributions of (S0 , L, m∗s , m∗v ) of the Bayesian analysis after updating the CI double ratio by replacing the pseudo-neutrons estimated using HiRA10 with free neutrons from vLANA. The diagonal plots show the marginal distributions of each 155 S0 170 140 L [MeV] 110 80 50 L 20 1.0 0.9 ms∗ /mN 0.8 0.7 ms∗ 0.6 1.2 1.0 mv∗ /mN 0.8 mv∗ 0.625 30 35 40 45 20 50 80 110 140 170 0.6 0.7 0.8 0.9 1.0 0.6 0.8 1.0 1.2 S0 [MeV] L [MeV] ms∗ /mN mv∗ /mN Figure 5.29 Posterior distributions of (S0 , L, m∗s , m∗v ) for the Ca + Ni systems at 140 MeV/u. The diagonal plots show the marginal distributions of each parameter, whereas the off-diagonal plots show the pairwise correlations between parameters. The units of S0 and L are MeV, whereas the units of m∗s and m∗v are mN . The red dashed line in the m∗v -m∗s correlation plot indicates the line where m∗v = m∗s . parameter, whereas the off-diagonal plots show the pairwise correlations between parameters. Uniform priors are used for all parameters; wide truncated normal priors were explored without any significant changes to the results. For reference, we have compiled the means and standard deviations of the posterior distributions for every parameter in Table 5.2. Similarly to Figure 3.19, the inference on CI double ratio using free neutron spectra also exhibits a correlation between m∗s and m∗v . A red line is drawn in the m∗s -m∗v panel of 156 0.11 ± 0.09 −0.10 ± 0.09 0.4 0.2 0.0 0.2 0.4 0.4 0.2 0.0 0.2 0.4 fI ∆mnp∗ /δ Figure 5.30 Left: Posterior distributions of fI . Right: Posterior distributions of ∆m∗np /δ. Figure 5.29 to indicate m∗v = m∗s . To better quantify the preference on the nucleon effective mass splitting, we have plotted Figure 5.30 by populating the histograms from the fI values computed for every posterior sample point in Figure 5.29. The analysis suggests an effective mass splitting of fI = 0.11 ± 0.09 (left panel) or ∆m∗np /δ = −0.10 ± 0.09 (right panel), with the errors representing the standard deviations of the posterior distributions. This result on effective mass splitting is consistent with the earlier result that involves pseudo-neutron spectra, which gives fI = 0.14 ± 0.11 or ∆m∗np /δ = −0.17 ± 0.12. In other words, despite the evident differences between the free neutron spectra and the pseudo-neutron spectra as shown in Figure 5.23, the Bayesian analyses on the CI double ratio of n/p are able to produce comparable results, hinting at a potential way to probe the nucleon effective mass splitting with only charged particles, especially when neutron measurement is often more challenging than charged particles measurement. In Figure 5.31, we have compiled the effective mass splitting ∆m∗np /δ values from this work and two previous works, both of which used the same ImQMD model to generate the theoretical spectra and similar Bayesian analysis to extract the effective mass splitting. The work by C. Y. Tsang et al. under the SπRIT collaboration [202] looks a different set 112 of observables from Sn + 124Sn and 108 Sn + 132Sn at 270 MeV/u, including the directed flow (v1 ), the elliptic flow (v2 ) [203], and the stopping (VarXZ) [204]. This work suggests ∆m∗np /δ = −0.08 ± 0.07, indicated by the open circle. The work by P. Morfouace et al. 112 [147] looks at the n/p ratios of Sn + 112Sn and 124 Sn + 124Sn at 120 MeV/u. The open 157 0.20 Sn + Sn at 270 MeV/u 0.15 Sn + Sn at 120 MeV/u δ 0.10 Ca + Ni at 140 MeV/u SπRIT (2023) 0.05 P. Morfouace, et al. (2019): SR+DR P. Morfouace, et al. (2019): DR This work 0.000.4 0.3 0.2 0.1 0.0 0.1 0.2 0.3 0.4 ∆mnp /δ ∗ Figure 5.31 The effective mass splitting ∆m∗np /δ from this work and previous works [147, 202]. square shows inference result from both the single ratio (SR) and the double ratio (DR), suggesting ∆m∗np /δ = −0.05 ± 0.09; the cross shows result from the DR only, suggesting ∆m∗np /δ = −0.03 ± 0.15. Lastly, we have the result from this work, indicated by a star 40 marker, which extracts ∆m∗np /δ = −0.10 ± 0.09 from the CI double n/p ratio of Ca + 58Ni 48 and Ca + 64Ni at 140 MeV/u. In summary, all three works, including the results from this thesis, suggest a slightly negative effective mass splitting, with m∗n < m∗p . Each of these studies utilized the ImQMD model, which has been thoroughly examined for self-consistency via the closure test described in subsection 3.4.2. However, it is important to note that this model is one of many other, each featuring its unique formulation and implementation of effective masses. Hence, results may exhibit variations when employing different theoretical models. Investigations using optical potentials [205, 206, 207] and those examining giant resonances and the electrical dipole polarizability of Pb [208, 209, 210, 211] yield a positive sign for the effective mass splitting. This contrast underscores how methodological choices, selected observables, and theoretical models significantly influence outcomes, attesting to the complex nature of the nuclear dynamics, rather than reflecting limitations of any particular study 158 or methodology. This, in turn, accentuates the importance of continual collaboration with theorists and the need for diversified investigations using theoretical models. 159 CHAPTER 6 SUMMARY One goal of this dissertation is to extract the nucleon effective mass splitting, a critical aspect in understanding the nuclear equation of state (EOS). The experiment measured the spectral yield ratios of neutrons and protons from heavy-ion collisions (HICs) and to compare them with theoretical predictions. Conducted at the National Superconducting Cyclotron Laboratory (NSCL), the experiment used three detection systems: the Microball impact parameter detector, the HiRA10 charged-particle detector, and the vLANA neutron detector. A total of 16 beam-target configurations involving calcium isotopes as the beam and nickel and tin isotopes as targets at two distinct beam energies, 56 MeV/u and 140MeV/u, were conducted. In this dissertation, beam-target pairs with the lowest and the highest (N − Z)/A 40 values for a given target element, namely, Ca + 58Ni (0.020), 48 Ca + 64Ni (0.143), 40 Ca + 112 48 Sn (0.079), and Ca + 124Sn (0.186), at both beam energies were analyzed. Beam-target pairs with intermediate (N − Z)/A values will be explored in future works. The first half of this work focuses on the results from charged particle detection by HiRA10. Microball was used to determine the impact parameter of the collisions and select the central collisions. In HiRA10, light charged particle spectra of protons, deuterons, tritons, helium-3, and helium-4 were measured and identified. Energy spectra in the lab frame for each of the particles were presented, followed by the transverse momentum spectra gated on the mid-rapidity region to select participant nucleons from the heavy-ion collisions. An isoscaling study was conducted on the transverse momentum spectra. The study revealed the isoscaling property for both beam energies at 56 MeV and 140 MeV. As expected, the magnitude of isoscaling parameters extracted from the lower beam energy is larger. Exploiting the isoscaling properties, “pseudo-neutrons” spectra were obtained and used to construct the coalescence-invariant (CI) neutron-to-proton (n/p) ratios. The CI n/p double ratio from the Ca + Ni systems was analyzed using Bayesian analysis by varying the EOS parameters of 160 the symmetry energy S0 , the slope parameter L, the isoscalar effective mass m∗s , and the isovector effective mass m∗v . Informed by the Improved Quantum Molecular (ImQMD) model, the Bayesian framework showed a high sensitivity to the nucleon effective mass splitting ∆m∗np /δ. To further validate this methodology, a closure test involving only simulated events from the model was performed and showed a strong correlation coefficient between the true and predicted values of isospin effective mass splitting. The second half of this work emphasizes on the data analysis of neutron detection by vLANA, a composite neutron detection system consisting of two neutron walls (NWA and NWB), a charged-particle veto wall (VW), and a Forward Array (FA), as well as the reconstruction of neutron spectra. A major part of the thesis work focused on the calibration of NWB in vLANA, including calibrations for hit position, time-of-flight, and light output, as well as non-linearity and saturation corrections of the analog-to-digital converter (ADC) values and neutron-gamma pulse shape discrimination (PSD). Following the comprehensive analysis and calibration of raw data from vLANA, the neutron spectra were reconstructed. This process began with an assessment of the background contributions to the spectra. Subsequently, a detection efficiency analysis was performed to establish both the geometric and intrinsic efficiencies. Determination of intrinsic efficiency involved sophisticated simulation of the performance of the neutron detector using a new simulation code, neuSIM4, which was developed in a collaboration with Korea University in conjunction with the neutron analysis. The work concluded with the presentation of neutron spectra for the Ca + Ni systems at 140 MeV/u and an updated isoscaling analysis. However, the measured neutron spectra and the pseudo-neutron spectra obtained from charged particles only using isoscaling are quite different. Nonetheless, our results confirmed that systematic errors associated with CI n can be cancelled out in the double ratio involving two reaction systems, but not in the individual single ratios of CI n/p. From the Bayesian analysis for the CI n/p double ratio using actual reconstructed neutrons, the effective mass splitting ∆m∗np /δ was estimated to be −0.10 ± 0.09. This result aligns with findings from the light charged particle analysis using 161 HiRA10 (∆m∗np /δ = −0.05 ± 0.09) and previous experiments (∆m∗np /δ = −0.08 ± 0.07) that used ImQMD, yet contrasts with studies employing optical potentials and those exploring giant resonances and the electrical dipole polarizability of Pb that suggest a positive sign. The implications of these findings permeate a broad spectrum of nuclear physics and astrophysics, influencing aspects such as the nuclear structure of exotic nuclei, the dynamics of heavy-ion collisions, and the description of neutron stars [207]. Evidence has been presented, for instance, showing that the neutrino mean free path, which is vital to the cooling process of neutron stars, is closely connected to the nucleon effective masses [76, 212]. Additionally, J. Estee et al. have presented a correlation between ∆m∗np /δ and the slope parameter L of the symmetry energy [91], a parameter that can be related to the pressure due to asymmetry and is central to the neutron star description. For example, an increase in L has been correlated with enhanced dynamical mass ejection during the post-merger phase of binary neutron star coalescences [213], and a larger radius in neutron stars above 1M⊙ [214]. Moreover, properties of the neutron star crust, such as its thickness, fractional mass, and moment of inertia, have been shown to exhibit significant sensitivity to L through the transition density [215]. In conclusion, this dissertation measured the splitting of the neutron and proton effective masses, important for the understanding of nuclear EOS, by combining observables of neutrons and protons from multiple detection systems. Two primary sources of uncertainty persist: the limited number of observable types, specifically, the CI n/p double ratio, and the use of a single theoretical model, the ImQMD model. Nonetheless, our results are consistent with results obtained from heavy-ion collisions using different observables but the same transport model. Collaboration with theorists to optimize the application of modern Bayesian tools for understanding the experimental data is ongoing. The work presented in this dissertation underscores the enormous effort put into handling over 900 hours of active beam time data to ensure seamless merging of data across different experimental runs. However, with the methodology developed in the thesis, efforts are in progress to finalize results from all reaction 40 systems, not only those with extreme asymmetry values, Ca + 58Ni (0.020) and 48 Ca + 64Ni 162 (0.143) discussed in this thesis, but also systems with intermediate asymmetry values. 163 BIBLIOGRAPHY [1] G. Gamow, “Mass defect curve and nuclear constitution,” Proceedings of the Royal Soci- ety of London. Series A, Containing Papers of a Mathematical and Physical Character, vol. 126, no. 803, pp. 632–644, 1930. doi: 10.1098/rspa.1930.0032. [2] Weizsäcker, CF v, “Zur Theorie der Kernmassen,” Zeitschrift für Physik, vol. 96, pp. 431–458, 1935. doi: 10.1007/BF01337700. [3] Benzaid, Djelloul and Bentridi, Salaheddine and Kerraci, Abdelkader and Amrani, Naima, “Bethe–Weizsäcker semiempirical mass formula coefficients 2019 update based on AME2016,” Nuclear Science and Techniques, vol. 31, no. 1, p. 9, 2020. doi: 10.1007/s41365-019-0718-8. [4] Huang, W J and Audi, G and Wang, Meng and Kondev, F G and Naimi, S and Xu, Xing, “The AME2016 atomic mass evaluation (I). Evaluation of input data; and adjustment procedures,” Chinese Physics C, vol. 41, no. 3, p. 030002, 2017. doi: 10.1088/1674-1137/41/3/030002. [5] Wang, Meng and Audi, G and Kondev, F G and Huang, W J and Naimi, S and Xu, Xing, “The AME2016 atomic mass evaluation (II). Tables, graphs and references,” Chinese Physics C, vol. 41, no. 3, p. 030003, 2021. doi: 10.1088/1674-1137/41/3/030003. [6] A. Zucker, “Nuclear interactions of heavy ions,” Annual Review of Nuclear Science, vol. 10, no. 1, pp. 27–62, 1960. doi: 10.1146/annurev.ns.10.120160.000331. [7] W. Scheid, H. Müller, and W. Greiner, “Nuclear shock waves in heavy-ion collisions,” Physical Review Letters, vol. 32, no. 13, p. 741, 1974. doi: 10.1103/PhysRevLett.32.741. [8] J. P. Bondorf, P. J. Siemens, S. Garpman, and E. Halbert, “Classical microscopic calculation of fast heavy-ion collisions,” Zeitschrift für Physik A Atoms and Nuclei, vol. 279, pp. 385–394, 1976. doi: 10.1007/BF01418134. [9] H. Stöcker, J. A. Maruhn, and W. Greiner, “Collective sideward flow of nuclear matter in violent high-energy heavy-ion collisions,” Physical Review Letters, vol. 44, no. 11, p. 725, 1980. doi: 10.1103/PhysRevLett.44.725. [10] D. T. Khoa, N. Ohtsuka, A. Faessler, M. A. Matin, S. W. Huang, E. Lehmann, and Y. Lofty, “Microscopic study of thermal properties of the nuclear matter formed in heavy-ion collisions,” Nuclear Physics A, vol. 542, no. 4, pp. 671–698, 1992. doi: 10.1016/0375-9474(92)90263-J. [11] P. Russotto, M. D. Cozma, E. De Filippo, A. Le Fèvre, Y. Leifels, and J. Lukasik, “Studies of the equation-of-state of nuclear matter by heavy-ion collisions at intermediate energy in the multi-messenger era: A review focused on GSI results,” La Rivista del Nuovo Cimento, vol. 46, no. 1, pp. 1–70, 2023. doi: 10.1007/s40766-023-00039-4. 164 [12] J. M. Lattimer and D. G. Ravenhall, “Neutron star matter at high temperatures and densities. I-Bulk properties of nuclear matter,” The Astrophysical Journal, vol. 223, pp. 314–323, 1978. doi: 10.1086/156265. [13] J. Boguta, “Remarks on the beta stability in neutron stars,” Physics Letters B, vol. 106, no. 4, pp. 255–258, 1981. doi: 10.1016/0370-2693(81)90529-3. [14] T. Takatsuka and R. Tamagaki, “Superfluidity in neutron star matter and symmetric nuclear matter,” Progress of Theoretical Physics Supplement, vol. 112, pp. 27–65, 1993. doi: 10.1143/PTP.112.27. [15] J. M. Lattimer and M. Prakash, “Nuclear matter and its role in supernovae, neutron stars and compact object binary mergers,” Physics Reports, vol. 333, pp. 121–146, 2000. doi: 10.1016/S0370-1573(00)00019-3. [16] J. M. Lattimer, “Neutron stars and the nuclear matter equation of state,” Annual Review of Nuclear and Particle Science, vol. 71, pp. 433–464, 2021. doi: 10.1146/annurev-nucl- 102419-124827. [17] E. Wigner, “On the Consequences of the Symmetry of the Nuclear Hamiltonian on the Spectroscopy of Nuclei,” Physical Review, vol. 51, no. 2, p. 106, 1937. doi: 10.1103/PhysRev.51.106. [18] D. D. Warner, M. A. Bentley, and P. V. Isacker, “The role of isospin symmetry in collective nuclear structure,” Nature Physics, vol. 2, no. 5, pp. 311–318, 2006. doi: 10.1038/nphys291. [19] E. Bravo and D. Garcı́a-Senz, “Coulomb corrections to the equation of state of nuclear statistical equilibrium matter: implications for SNIa nucleosynthesis and the accretion- induced collapse of white dwarfs,” Monthly Notices of the Royal Astronomical Society, vol. 307, no. 4, pp. 984–992, 1999. doi: 10.1046/j.1365-8711.1999.02694.x. [20] L. Qin, K. Hagel, R. Wada, J. B. Natowitz, S. Shlomo, A. Bonasera, G. Röpke, S. Typel, Z. Chen, M. Huang et al., “Laboratory tests of low density astrophysical nuclear equations of state,” Physical Review Letters, vol. 108, no. 17, p. 172701, 2012. doi: 10.1103/PhysRevLett.108.172701. [21] X. Roca-Maza and N. Paar, “Nuclear equation of state from ground and collective excited state properties of nuclei,” Progress in Particle and Nuclear Physics, vol. 101, pp. 96–176, 2018. doi: 10.1016/j.ppnp.2018.04.001. [22] R. D. Puff, “Ground-state properties of nuclear matter,” Annals of Physics, vol. 13, no. 3, pp. 317–358, 1961. doi: 10.1016/0003-4916(61)90149-X. [23] D. S. Falk and L. Wilets, “Nuclear compressibility and symmetry energy,” Physical Review, vol. 124, no. 6, p. 1887, 1961. doi: 10.1103/PhysRev.124.1887. [24] J. P. Blaizot, D. Gogny, and B. Grammaticos, “Nuclear compressibility and monopole resonances,” Nuclear Physics A, vol. 265, no. 2, pp. 315–336, 1976. doi: 10.1016/0375- 9474(76)90357-2. 165 [25] J.-P. Blaizot, “Nuclear compressibilities,” Physics Reports, vol. 64, no. 4, pp. 171–248, 1980. doi: 10.1016/0370-1573(80)90001-0. [26] D. H. Youngblood, H. L. Clark, and Y.-W. Lui, “Incompressibility of nuclear matter from the giant monopole resonance,” Physical Review Letters, vol. 82, no. 4, p. 691, 1999. doi: 10.1103/PhysRevLett.82.691. [27] T. Li, U. Garg, Y. Liu, R. Marks, B. K. Nayak, P. V. M. Rao, M. Fujiwara, H. Hashimoto, K. Nakanishi, S. Okumura et al., “Isoscalar giant resonances in the Sn nuclei and implications for the asymmetry term in the nuclear-matter incompressibility,” Physical Review C, vol. 81, no. 3, p. 034309, 2010. doi: 10.1103/PhysRevC.81.034309. [28] M. Dutra, O. Lourenço, J. S. S. Martins, A. Delfino, J. R. Stone, and P. D. Stevenson, “Skyrme interaction and nuclear matter constraints,” Physical Review C, vol. 85, no. 3, p. 035201, 2012. doi: 10.1103/PhysRevC.85.035201. [29] J. W. Holt and Y. Lim, “Universal correlations in the nuclear symmetry energy, slope parameter, and curvature,” Physics Letters B, vol. 784, pp. 77–81, 2018. doi: 10.1016/j.physletb.2018.07.038. [30] P. Danielewicz and J. Lee, “Symmetry energy II: Isobaric analog states,” Nuclear Physics A, vol. 922, pp. 1–70, 2014. doi: 10.1016/j.nuclphysa.2013.11.005. [31] A. Klimkiewicz, N. Paar, P. Adrich, M. Fallot, K. Boretzky, T. Aumann, D. Cortina-Gil, U. D. Pramanik, T. W. Elze, H. Emling et al., “Nuclear symmetry energy and neutron skins derived from pygmy dipole resonances,” Physical Review C, vol. 76, no. 5, p. 051603, 2007. doi: 10.1103/PhysRevC.76.051603. [32] A. Carbone, G. Colo, A. Bracco, L.-G. Cao, P. F. Bortignon, F. Camera, and O. Wieland, “Constraints on the symmetry energy and neutron skins from pygmy resonances in 68 Ni and 132 Sn,” Physical Review C, vol. 81, no. 4, p. 041301, 2010. doi: 10.1103/Phys- RevC.81.041301. [33] A. Tamii, I. Poltoratska, P. von Neumann-Cosel, Y. Fujita, T. Adachi, C. A. Bertulani, J. Carter, M. Dozono, H. Fujita, K. Fujita et al., “Complete electric dipole response and the neutron skin in 208 Pb,” Physical review letters, vol. 107, no. 6, p. 062502, 2011. doi: 10.1103/PhysRevLett.107.062502. [34] D. M. Rossi, P. Adrich, F. Aksouh, H. Alvarez-Pol, T. Aumann, J. Benlliure, M. Böhmer, K. Boretzky, E. Casarejos, M. Chartier et al., “Measurement of the Dipole Polarizability of the Unstable Neutron-Rich Nucleus 68 Ni,” Physical review letters, vol. 111, no. 24, p. 242503, 2013. doi: 10.1103/PhysRevLett.111.242503. [35] M. B. Tsang, T. X. Liu, L. Shi, P. Danielewicz, C. K. Gelbke, X. D. Liu, W. G. Lynch, W. P. Tan, G. Verde, A. Wagner et al., “Isospin diffusion and the nuclear symmetry energy in heavy ion reactions,” Physical review letters, vol. 92, no. 6, p. 062701, 2004. doi: 10.1103/PhysRevLett.92.062701. 166 [36] Z. Y. Sun, M. B. Tsang, W. G. Lynch, G. Verde, F. Amorini, L. Andronenko, M. An- dronenko, G. Cardella, M. Chatterje, P. Danielewicz et al., “Isospin diffusion and equilibration for Sn+Sn collisions at E/A = 35 MeV,” Physical Review C, vol. 82, no. 5, p. 051603, 2010. doi: 10.1103/PhysRevC.82.051603. [37] M. Di Toro, V. Baran, M. Colonna, and V. Greco, “Probing the nuclear symmetry energy with heavy-ion collisions,” Journal of Physics G: Nuclear and Particle Physics, vol. 37, no. 8, p. 083101, 2010. doi: 10.1088/0954-3899/37/8/083101. [38] A. W. Steiner, M. Hempel, and T. Fischer, “Core-collapse supernova equations of state based on neutron star observations,” The Astrophysical Journal, vol. 774, no. 1, p. 17, 2013. [39] J. M. Lattimer, “Symmetry energy in nuclei and neutron stars,” Nuclear Physics A, vol. 928, pp. 276–295, 2014. doi: 10.1016/j.nuclphysa.2014.04.008. [40] J. R. Oppenheimer and G. M. Volkoff, “On massive neutron cores,” Physical Review, vol. 55, no. 4, p. 374, 1939. doi: 10.1103/PhysRev.55.374. [41] J. G. Martinez, K. Stovall, P. C. C. Freire, J. S. Deneva, F. A. Jenet, M. A. McLaughlin, M. Bagchi, S. D. Bates, and A. Ridolfi, “Pulsar J0453+1559: A double neutron star system with a large mass asymmetry,” The Astrophysical Journal, vol. 812, no. 2, p. 143, 2015. doi: 10.1088/0004-637X/812/2/143. [42] H. T. Cromartie, E. Fonseca, S. M. Ransom, P. B. Demorest, Z. Arzoumanian, H. Blumer, P. R. Brook, M. E. DeCesar, T. Dolch, J. A. Ellis et al., “Relativis- tic Shapiro delay measurements of an extremely massive millisecond pulsar,” Nature Astronomy, vol. 4, no. 1, pp. 72–76, 2020. doi: 10.1038/s41550-019-0880-2. [43] Y. Suwa, T. Yoshida, M. Shibata, H. Umeda, and K. Takahashi, “On the minimum mass of neutron stars,” Monthly Notices of the Royal Astronomical Society, vol. 481, no. 3, pp. 3305–3312, 2018. doi: 10.1093/mnras/sty2460. [44] F. Özel, D. Psaltis, R. Narayan, and A. S. Villarreal, “On the mass distribution and birth masses of neutron stars,” The Astrophysical Journal, vol. 757, no. 1, p. 55, 2012. doi: 10.1088/0004-637X/757/1/55. [45] J. M. Lattimer and M. Prakash, “Neutron star structure and the equation of state,” The Astrophysical Journal, vol. 550, no. 1, p. 426, 2001. doi: 10.1086/319702. [46] K. Hebeler and A. Schwenk, “Symmetry energy, neutron skin, and neutron star radius from chiral effective field theory interactions,” The European Physical Journal A, vol. 50, pp. 1–7, 2014. doi: 10.1140/epja/i2014-14011-4. [47] A. Worley, P. G. Krastev, and B.-A. Li, “Nuclear constraints on the moments of inertia of neutron stars,” The Astrophysical Journal, vol. 685, no. 1, p. 390, 2008. doi: 10.1086/589823. 167 [48] F. Fattoyev and J. Piekarewicz, “Sensitivity of the moment of inertia of neutron stars to the equation of state of neutron-rich matter,” Physical Review C, vol. 82, no. 2, p. 025810, 2010. doi: 10.1103/PhysRevC.82.025810. [49] A. W. Steiner, “High-density symmetry energy and direct Urca process,” Physical Review C, vol. 74, no. 4, p. 045808, 2006. doi: 10.1103/PhysRevC.74.045808. [50] H. S. Than, D. T. Khoa, and N. Van Giai, “Neutron star cooling: A challenge to the nuclear mean field,” Physical Review C, vol. 80, no. 6, p. 064312, 2009. doi: 10.1103/PhysRevC.80.064312. [51] B. P. Abbott, R. Abbott, T. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya et al., “GW170817: observation of gravitational waves from a binary neutron star inspiral,” Physical review letters, vol. 119, no. 16, p. 161101, 2017. doi: 10.1103/PhysRevX.9.011001. [52] M. B. Tsang, W. G. Lynch, P. Danielewicz, and C. Y. Tsang, “Symmetry energy constraints from GW170817 and laboratory experiments,” Physics Letters B, vol. 795, pp. 533–536, 2019. doi: 10.1016/j.physletb.2019.06.059. [53] C. Y. Tsang, M. B. Tsang, P. Danielewicz, F. J. Fattoyev, and W. G. Lynch, “Insights on Skyrme parameters from GW170817,” Physics Letters B, vol. 796, pp. 1–5, 2019. doi: 10.1016/j.physletb.2019.05.055. [54] B.-A. Li, B.-J. Cai, W.-J. Xie, and N.-B. Zhang, “Progress in constraining nuclear symmetry energy using neutron star observables since GW170817,” Universe, vol. 7, no. 6, p. 182, 2021. doi: 10.3390/universe7060182. [55] B. A. Brown, “Neutron radii in nuclei and the neutron equation of state,” Physical review letters, vol. 85, no. 25, p. 5296, 2000. doi: 10.1103/PhysRevLett.85.5296. [56] M. B. Tsang, J. R. Stone, F. Camera, P. Danielewicz, S. Gandolfi, K. Hebeler, C. J. Horowitz, J. Lee, W. G. Lynch, Z. Kohley et al., “Constraints on the symmetry energy and neutron skins from experiments and theory,” Physical Review C, vol. 86, no. 1, p. 015803, 2012. doi: 10.1103/PhysRevC.86.015803. [57] H. Hergert, “A guided tour of ab initio nuclear many-body theory,” Frontiers in Physics, vol. 8, p. 379, 2020. doi: 10.3389/fphy.2020.00379. [58] D. R. Hartree, “The wave mechanics of an atom with a non-Coulomb central field. Part I. Theory and methods,” in Mathematical Proceedings of the Cambridge Philosophical Soci- ety, vol. 24, no. 1. Cambridge university press, 1928. doi: 10.1017/S0305004100011919. pp. 89–110. [59] D. R. Hartree, “The wave mechanics of an atom with a non-coulomb central field. Part II. Some results and discussion,” in Mathematical Proceedings of the Cam- bridge Philosophical Society, vol. 24, no. 1. Cambridge University Press, 1928. doi: 10.1017/S0305004100011920. pp. 111–132. 168 [60] M. Rayet, “Skyrme parametrization of an effective λ-nucleon interaction,” Nuclear Physics A, vol. 367, no. 3, pp. 381–397, 1981. doi: 10.1016/0003-4916(76)90262-1. [61] L. S. Celenza and C. M. Shakin, “Relativistic many-body theory,” Physical Review C, vol. 24, no. 6, p. 2704, 1981. doi: 10.1103/PhysRevC.24.2704. [62] K. A. Brueckner, “Two-body forces and nuclear saturation. III. Details of the structure of the nucleus,” Physical Review, vol. 97, no. 5, p. 1353, 1955. doi: 10.1103/Phys- Rev.97.1353. [63] J. P. Jeukenne, A. Lejeune, and C. Mahaux, “Many-body theory of nuclear matter,” Physics Reports, vol. 25, no. 2, pp. 83–174, 1976. doi: 10.1016/0370-1573(76)90017-X. [64] C. Mahaux, P. F. Bortignon, R. A. Broglia, and C. H. Dasso, “Dynamics of the shell model,” Physics Reports, vol. 120, no. 1-4, pp. 1–274, 1985. doi: 10.1016/0370- 1573(85)90100-0. [65] M. Jaminon and C. Mahaux, “Effective masses in relativistic approaches to the nucleon- nucleus mean field,” Physical Review C, vol. 40, no. 1, p. 354, 1989. doi: 10.1103/Phys- RevC.40.354. [66] J. M. Pearson and S. Goriely, “Isovector effective mass in the Skyrme-Hartree-Fock method,” Physical Review C, vol. 64, no. 2, p. 027301, 2001. doi: 10.1103/Phys- RevC.64.027301. [67] G. E. Brown, J. H. Gunn, and P. Gould, “Effective mass in nuclei,” Nuclear Physics, vol. 46, pp. 598–606, 1963. doi: 10.1016/0029-5582(63)90631-X. [68] B.-A. Li, B.-J. Cai, L.-W. Chen, W.-J. Xie, J. Xu, and N.-B. Zhang, “A theoretical overview of isospin and EOS effects in heavy-ion reactions at intermediate energies,” Il Nuovo Cimento C, vol. 45, p. 54, 2022. doi: 10.1393/ncc/i2022-22054-3. [69] O. Sjöberg, “On the landau effective mass in asymmetric nuclear matter,” Nuclear Physics A, vol. 265, no. 3, pp. 511–516, 1976. doi: 10.1016/0375-9474(76)90558-3. [70] W. Zuo, A. Lejeune, U. Lombardo, and J.-F. Mathiot, “Microscopic three-body force for asymmetric nuclear matter,” The European Physical Journal A-Hadrons and Nuclei, vol. 14, pp. 469–475, 2002. doi: 10.1140/epja/i2002-10031-y. [71] F. Hofmann, C. Keil, and H. Lenske, “Density dependent hadron field theory for asymmetric nuclear matter and exotic nuclei,” Physical Review C, vol. 64, no. 3, p. 034314, 2001. doi: 10.1103/PhysRevC.64.034314. [72] B. Liu, V. Greco, V. Baran, M. Colonna, and M. Di Toro, “Asymmetric nuclear matter: The role of the isovector scalar channel,” Physical Review C, vol. 65, no. 4, p. 045201, 2002. doi: 10.1103/PhysRevC.65.045201. [73] G. Steigman, “Primordial nucleosynthesis: successes and challenges,” Interna- tional Journal of Modern Physics E, vol. 15, no. 01, pp. 1–35, 2006. doi: 10.1142/S0218301306004028. 169 [74] D. G. Yakovlev, A. D. Kaminker, O. Y. Gnedin, and P. Haensel, “Neutrino emission from neutron stars,” Physics Reports, vol. 354, no. 1-2, pp. 1–155, 2001. doi: 10.1016/S0370- 1573(00)00131-9. [75] D. Page and S. Reddy, “Dense matter in compact stars: theoretical developments and observational constraints,” Annu. Rev. Nucl. Part. Sci., vol. 56, pp. 327–374, 2006. doi: 10.1146/annurev.nucl.56.080805.140600. [76] M. Baldo, G. F. Burgio, H.-J. Schulze, and G. Taranto, “Nucleon effective masses within the Brueckner-Hartree-Fock theory: Impact on stellar neutrino emission,” Physical Review C, vol. 89, no. 4, p. 048801, 2014. doi: 10.1103/PhysRevC.89.048801. [77] J. A. Nolen Jr and J. P. Schiffer, “Coulomb energies,” Annual Review of Nuclear Science, vol. 19, no. 1, pp. 471–526, 1969. doi: 10.1146/annurev.ns.19.120169.002351. [78] P. J. Woods and C. N. Davids, “Nuclei beyond the proton drip-line,” Annual Review of Nuclear and Particle Science, vol. 47, no. 1, pp. 541–590, 1997. doi: 10.1146/an- nurev.nucl.47.1.541. [79] B.-A. Li, C. B. Das, S. D. Gupta, and C. Gale, “Effects of momentum-dependent symmetry potential on heavy-ion collisions induced by neutron-rich nuclei,” Nuclear Physics A, vol. 735, no. 3-4, pp. 563–584, 2004. doi: 10.1016/j.nuclphysa.2004.02.016. [80] J. Rizzo, M. Colonna, and M. Di Toro, “Fast nucleon emission as a probe of the momentum dependence of the symmetry potential,” Physical Review C, vol. 72, no. 6, p. 064609, 2005. doi: 10.1103/PhysRevC.72.064609. [81] V. Giordano, M. Colonna, M. Di Toro, V. Greco, and J. Rizzo, “Isospin emission and flow at high baryon density: A test of the symmetry potential,” Physical Review C, vol. 81, no. 4, p. 044611, 2010. doi: 10.1103/PhysRevC.81.044611. [82] W.-J. Xie and F.-S. Zhang, “Nuclear collective flows as a probe to the neutron– proton effective mass splitting,” Physics Letters B, vol. 735, pp. 250–255, 2014. doi: 10.1016/j.physletb.2014.06.050. [83] Y. Zhang, M. B. Tsang, Z. Li, and H. Liu, “Constraints on nucleon effective mass splitting with heavy ion collisions,” Physics Letters B, vol. 732, pp. 186–190, 2014. doi: 10.1016/j.physletb.2014.03.030. [84] Y. Zhang, M. B. Tsang, and Z. Li, “Covariance analysis of symmetry energy observ- ables from heavy ion collision,” Physics Letters B, vol. 749, pp. 262–266, 2015. doi: 10.1016/j.physletb.2015.07.064. [85] Y. Zhang and Z. Li, “Elliptic flow and system size dependence of transition energies at intermediate energies,” Physical Review C, vol. 74, no. 1, p. 014602, 2006. doi: 10.1103/PhysRevC.74.014602. [86] P. Danielewicz, R. Lacey, and W. G. Lynch, “Determination of the equation of state of dense matter,” Science, vol. 298, no. 5598, pp. 1592–1596, 2002. doi: 10.1126/sci- ence.1078070. 170 [87] D. D. S. Coupland, M. Youngs, Z. Chajecki, W. G. Lynch, M. B. Tsang, Y. X. Zhang, M. A. Famiano, T. K. Ghosh, B. Giacherio, M. A. Kilburn et al., “Probing effective nucleon masses with heavy-ion collisions,” Physical Review C, vol. 94, no. 1, p. 011601, 2016. doi: 10.1103/PhysRevC.94.011601. [88] M. B. Tsang, Y. Zhang, P. Danielewicz, M. Famiano, Z. Li, W. G. Lynch, A. W. Steiner et al., “Constraints on the density dependence of the symmetry energy,” Physical review letters, vol. 102, no. 12, p. 122701, 2009. doi: 10.1103/PhysRevLett.102.122701. [89] Y. Zhang, P. Danielewicz, M. Famiano, Z. Li, W. G. Lynch, and M. B. Tsang, “The influence of cluster emission and the symmetry energy on neutron–proton spec- tral double ratios,” Physics Letters B, vol. 664, no. 1-2, pp. 145–148, 2008. doi: 10.1016/j.physletb.2008.03.075. [90] M. A. Famiano, T. Liu, W. G. Lynch, M. Mocko, A. M. Rogers, M. B. Tsang, M. S. Wallace, R. J. Charity, S. Komarov, D. G. Sarantites et al., “Neutron and proton transverse emission ratio measurements and the density dependence of the asymmetry term of the nuclear equation of state,” Physical review letters, vol. 97, no. 5, p. 052701, 2006. doi: 10.1103/PhysRevLett.97.052701. [91] J. Estee, W. G. Lynch, C. Y. Tsang, J. Barney, G. Jhang, M. B. Tsang, R. Wang, M. Kaneko, J. W. Lee, T. Isobe et al., “Probing the symmetry energy with the spectral pion ratio,” Physical review letters, vol. 126, no. 16, p. 162701, 2021. doi: 10.1103/PhysRevLett.126.162701. [92] Z. Xiao, B.-A. Li, L.-W. Chen, G.-C. Yong, and M. Zhang, “Circumstantial evidence for a soft nuclear symmetry energy at suprasaturation densities,” Physical Review Letters, vol. 102, no. 6, p. 062502, 2009. doi: 10.1103/PhysRevLett.102.062502. [93] M. Zhang, Z.-G. Xiao, B.-A. Li, L.-W. Chen, G.-C. Yong, S.-J. Zhu et al., “System- atic study of the π-/π+ ratio in heavy-ion collisions with the same neutron/proton ratio but different masses,” Physical Review C, vol. 80, no. 3, p. 034616, 2009. doi: 10.1103/PhysRevC.80.034616. [94] G. Ferini, T. Gaitanos, M. Colonna, M. Di Toro, and H. Wolter, “Isospin effects on subthreshold kaon production at intermediate energies,” Physical review letters, vol. 97, no. 20, p. 202301, 2006. [95] C. Hartnack, L. Zhuxia, L. Neise, G. Peilert, A. Rosenhauer, H. Sorge, J. Aichelin, H. Stöcker, and W. Greiner, “Quantum molecular dynamics a microscopic model from UNILAC to CERN energies,” Nuclear Physics A, vol. 495, no. 1-2, pp. 303–319, 1989. doi: 10.1016/0375-9474(89)90328-X. [96] J. Aichelin, ““Quantum” molecular dynamics—a dynamical microscopic n-body ap- proach to investigate fragment formation and the nuclear equation of state in heavy ion collisions,” Physics Reports, vol. 202, no. 5-6, pp. 233–360, 1991. doi: 10.1016/0370- 1573(91)90094-3. 171 [97] Z. Chajecki et al., “Experiment e14030/e15190 - RUN LOG,” 2018, accessed: 2023-05- 01. [Online]. Available: https://github.com/fanurs/data-analysis-e15190-e14030/blob/ 6fff140570940fb1817afdeba5775e191764dfbe/database/runlog/downloads/elog.html [98] Zecher, P D and Galonsky, A and Kruse, J J and Gaff, S J and Ottarson, J and Wang, J and Deak, F and Horvath, A and Kiss, A and Seres, Z and others, “A large-area, position-sensitive neutron detector with neutron/γ-ray discrimination capabilities,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spec- trometers, Detectors and Associated Equipment, vol. 401, no. 2-3, pp. 329–344, 1997. doi: 10.1016/S0168-9002(97)00942-X. [99] D. G. Sarantites, P.-F. Hua, M. Devlin, L. G. Sobotka, J. Elson, J. T. Hood, D. R. LaFosse, J. E. Sarantites, and M. R. Maier, ““The microball” Design, instrumentation and response characteristics of a 4π-multidetector exit channel-selection device for spectroscopic and reaction mechanism studies with Gammasphere,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 381, no. 2-3, pp. 418–432, 1996. doi: 10.1016/S0168- 9002(96)00785-1. [100] P. Chung, N. N. Ajitanand, J. M. Alexander, J. Ames, M. Anderson, D. Best, F. P. Brady, T. Case, W. Caskey, D. Cebra et al., “Centrality and momentum-selected elliptic flow: tighter constraints for the nuclear equation of state,” Physical Review C, vol. 66, no. 2, p. 021901, 2002. doi: 10.1103/physrevc.66.021901. [101] A. Andronic, V. Barret, Z. Basrak, N. Bastid, L. Benabderrahmane, G. Berek, R. Čaplar, E. Cordier, P. Crochet, P. Dupieux et al., “Excitation function of elliptic flow in Au+Au collisions and the nuclear matter equation of state,” Physics Letters B, vol. 612, no. 3-4, pp. 173–180, 2005. doi: 10.1016/j.physletb.2005.02.060. [102] T. Z. Yan and S. Li, “Yield ratios of light particles as a probe of the proton skin of a nucleus and its centrality dependence,” Physical Review C, vol. 101, no. 5, p. 054601, 2020. doi: 10.1103/physrevc.101.054601. [103] S. Schmidt, C. Kelbch, H. Schmidt-Böcking, and G. Kraft, Delta-electron emission in heavy ion collisions. Springer, 1988. [104] F. Benrachi, B. Chambon, B. Cheynis, D. Drain, C. Pastor, D. Seghier, K. Zaid, A. Giorni, D. Heuer, A. Lleres et al., “Investigation of the performance of CsI(Tl) for charged particle identification by pulse-shape analysis,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Asso- ciated Equipment, vol. 281, no. 1, pp. 137–142, 1989. doi: 10.1016/0168-9002(89)91225-4. [105] M. Alderighi, A. Anzalone, R. Basssini, I. Berceanu, J. Blicharska, C. Boiano, B. Bor- derie, R. Bougault, M. Bruno, C. Calı́ et al., “Particle identification method in the CsI(Tl) scintillator used for the CHIMERA 4π detector,” Nuclear Instruments and Meth- ods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 489, no. 1-3, pp. 257–265, 2002. doi: 10.1016/S0168-9002(02)00800-8. 172 [106] S. A. Bass, A. Bischoff, J. A. Maruhn, H. Stoecker, and W. Greiner, “Neural networks for impact parameter determination,” Physical Review C, vol. 53, no. 5, p. 2358, 1996. doi: 10.1103/PhysRevC.53.2358. [107] J. Xu, S. Yu, F. Liu, X. Luo et al., “Cumulants of net-proton, net-kaon, and net-charge √ multiplicity distributions in Au+Au collisions at sN N = 7.7, 11.5, 19.6, 27, 39, 62.4, and 200 GeV within the UrQMD model,” Physical Review C, vol. 94, no. 2, p. 024901, 2016. doi: 10.1103/PhysRevC.94.024901. [108] L. Li, Y. Zhang, Z. Li, N. Wang, Y. Cui, and J. Winkelbauer, “Impact parameter smearing effects on isospin sensitive observables in heavy ion collisions,” Physical Review C, vol. 97, no. 4, p. 044606, 2018. doi: 10.1103/PhysRevC.97.044606. [109] J. D. Frankland, D. Gruyer, E. Bonnet, B. Borderie, R. Bougault, A. Chbihi, J. E. Ducret, D. Durand, Q. Fable, M. Henri et al., “Model independent reconstruction of impact parameter distributions for intermediate energy heavy ion collisions,” Physical Review C, vol. 104, no. 3, p. 034609, 2021. doi: 10.1103/PhysRevC.104.034609. [110] S. Sweany, “Constraining the proton/neutron effective mass splitting through heavy ion collisions,” Ph.D. dissertation, Michigan State University, 2020, dissertation. [111] M. S. Wallace, M. A. Famiano, M.-J. Van Goethem, A. M. Rogers, W. G. Lynch, J. Clifford, F. Delaunay, J. Lee, S. Labostov, M. Mocko et al., “The high resolution array (HiRA) for rare isotope beam experiments,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 583, no. 2-3, pp. 302–312, 2007. doi: 10.1016/j.nima.2007.08.248. [112] S. Sweany, W. G. Lynch, K. Brown, A. Anthony, Z. Chajecki, D. Dell’Aquila, P. Mor- fouace, F. C. E. Teh, C. Y. Tsang, M. B. Tsang et al., “Reaction losses of charged particles in CsI(Tl) crystals,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 1018, p. 165798, 2021. doi: 10.1016/j.nima.2021.165798. [113] M. S. Wallace, “Experimental and Theoretical Challenges in Understanding rp-Process on Accreting Neutron Stars,” Ph.D. dissertation, Michigan State University, 2005, dissertation. [114] P. S. Fisher and D. K. Scott, “A fast computing circuit for particle identification using signals from a ∆E, E counter telescope,” Nuclear Instruments and Methods, vol. 49, no. 2, pp. 301–308, 1967. doi: 10.1016/0029-554X(67)90700-8. [115] A. G. Seamster, R. E. L. Green, and R. G. Korteling, “Silicon detector ∆E, E particle identification: a theoretically based analysis algorithm and remarks on the fundamental limits to the resolution of particle type by ∆E, E measurements,” Nuclear Instruments and Methods, vol. 145, no. 3, pp. 583–591, 1977. doi: 10.1016/0029-554X(77)90590-0. 173 [116] M. Hauschild, “Progress in dE/dx techniques used for particle identification,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrome- ters, Detectors and Associated Equipment, vol. 379, no. 3, pp. 436–441, 1996. doi: 10.1016/0168-9002(96)00607-9. [117] H. A. Bethe, J. Ashkin et al., “Experimental nuclear physics,” Wiley, New York, 1953. [118] O. B. Tarasov and D. Bazin, “LISE++: Radioactive beam production with in-flight separators,” Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, vol. 266, no. 19-20, pp. 4657–4664, 2008. doi: 10.1016/j.nimb.2008.05.110. [119] LISE++ group at NSCL/FRIB, “LISE++,” accessed: 2023-05-01. [Online]. Available: https://lise.nscl.msu.edu/lise.html [120] D. Dell’Aquila, S. Sweany, K. W. Brown, Z. Chajecki, W. G. Lynch, F. C. E. Teh, C.-Y. Tsang, M. B. Tsang, K. Zhu, C. Anderson et al., “Non-linearity effects on the light-output calibration of light charged particles in CsI(Tl) scintillator crystals,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spec- trometers, Detectors and Associated Equipment, vol. 929, pp. 162–172, 2019. doi: 10.1016/j.nima.2019.03.065. [121] J. Manfredi, J. Lee, W. G. Lynch, C. Y. Niu, M. B. Tsang, C. Anderson, J. Barney, K. W. Brown, Z. Chajecki, K. P. Chan et al., “On determining dead layer and detector thicknesses for a position-sensitive silicon detector,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 888, pp. 177–183, 2018. doi: 10.1016/j.nima.2017.12.082. [122] K. Zhu, M. B. Tsang, D. Dell’Aquila, K. W. Brown, Z. Chajecki, W. G. Lynch, S. Sweany, F. C. E. Teh, C. Y. Tsang, C. Anderson et al., “Calibration of large neutron detection arrays using cosmic rays,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 967, p. 163826, 2020. doi: 10.1016/j.nima.2020.163826. [123] K. Zhu, “Measuring neutrons in heavy ion collisions,” Ph.D. dissertation, Michigan State University, 2020, dissertation. [124] M. L. Roush, M. A. Wilson, and W. F. Hornyak, “Pulse shape discrimination,” Nu- clear Instruments and Methods, vol. 31, no. 1, pp. 112–124, 1964. doi: 10.1016/0029- 554X(64)90333-7. [125] Teledyne Technologies, “Teledyne LeCroy,” accessed: 2023-05-01. [Online]. Available: https://teledynelecroy.com/ [126] CAEN, “Technical Information Manual – V1190A/B, VX1190 A/B,” 2006, accessed: 2023-05-01. [Online]. Available: https://web.pa.msu.edu/people/edmunds/HAWC/ Manuals/caen v1190a tdc manual.pdf 174 [127] CAEN, “V862 – 32 Channel Multievent Individual Gate QDC,” accessed: 2023-05-01. [Online]. Available: https://www.caen.it/products/v862/ [128] CAEN, “V792 – 32 Channel Multievent QDC,” accessed: 2023-05-01. [Online]. Available: https://www.caen.it/products/v792/ [129] F. C. E. Teh, J.-W. Lee, K. Zhu, K. W. Brown, Z. Chajecki, W. G. Lynch, M. B. Tsang, A. Anthony, J. Barney, D. Dell’Aquila et al., “Value-assigned pulse shape discrimination for neutron detectors,” IEEE Transactions on Nuclear Science, vol. 68, no. 8, pp. 2294–2300, 2021. doi: 10.1109/TNS.2021.3091126. [130] H. Shim, J.-W. Lee, B. Hong, J. K. Ahn, G. Bak, J. Jo, M. Kim, Y. J. Kim, Y. J. Kim, H. Lee et al., “Performance of prototype neutron detectors for Large Acceptance Multi-Purpose Spectrometer at RAON,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 927, pp. 280–286, 2019. doi: 10.1016/j.nima.2019.02.064. [131] Eljen Technology, “General purpose plastic scintillator EJ-200, EJ-204, EJ-208, EJ-212,” 2021, accessed: 2023-05-01. [Online]. Available: https://eljentechnology.com/ products/plastic-scintillators/ej-200-ej-204-ej-208-ej-212 [132] H. H. Gutbrod, K. H. Kampert, B. W. Kolb, A. M. Poskanzer, H. G. Ritter, and H. R. Schmidt, “A new component of the collective flow in relativistic heavy-ion collisions,” Physics Letters B, vol. 216, no. 3-4, pp. 267–271, 1989. doi: 10.1016/0370-2693(89)91113- 1. [133] H. H. Gutbrod, K. H. Kampert, B. Kolb, A. M. Poskanzer, H. G. Ritter, R. Schicker, and H. R. Schmidt, “Squeeze-out of nuclear matter as a function of projectile energy and mass,” Physical Review C, vol. 42, no. 2, p. 640, 1990. doi: 10.1103/PhysRevC.42.640. [134] M. B. Tsang, P. Danielewicz, W. C. Hsi, M. Huang, W. G. Lynch, D. R. Bowman, C. K. Gelbke, M. A. Lisa, G. F. Peaslee, R. J. Charity et al., “Squeeze-out of nuclear matter in Au+Au collisions,” Physical Review C, vol. 53, no. 4, p. 1959, 1996. doi: 10.1103/PhysRevC.53.1959. [135] P. Danielewicz and G. Odyniec, “Transverse momentum analysis of collective motion in relativistic nuclear collisions,” Physics Letters B, vol. 157, no. 2-3, pp. 146–150, 1985. doi: 10.1016/0370-2693(85)91535-7. [136] H. S. Xu, M. B. Tsang, T. X. Liu, X. D. Liu, W. Lynch, W. P. Tan, A. Vander Molen, G. Verde, A. Wagner, H. F. Xi et al., “Isospin fractionation in nuclear multifrag- mentation,” Physical Review Letters, vol. 85, no. 4, p. 716, 2000. doi: 10.1103/Phys- RevLett.85.716. [137] M. B. Tsang, W. A. Friedman, C. K. Gelbke, W. G. Lynch, G. Verde, and H. S. Xu, “Isotopic scaling in nuclear reactions,” Physical Review Letters, vol. 86, no. 22, p. 5023, 2001. doi: 10.1103/PhysRevLett.86.5023. 175 [138] M. B. Tsang, C. K. Gelbke, X. D. Liu, W. G. Lynch, W. P. Tan, G. Verde, H. S. Xu, W. A. Friedman, R. Donangelo, S. R. Souza et al., “Isoscaling in statistical models,” Physical Review C, vol. 64, no. 5, p. 054615, 2001. doi: 10.1103/PhysRevC.64.054615. [139] J. W. Lee, M. B. Tsang, C. Y. Tsang, R. Wang, J. Barney, J. Estee, T. Isobe, M. Kaneko, M. Kurata-Nishimura, W. G. Lynch et al., “Isoscaling in central Sn+ Sn collisions at 270 MeV/u,” The European Physical Journal A, vol. 58, no. 10, pp. 1–10, 2022. doi: 10.1140/epja/s10050-022-00851-2. [140] J. Randrup and S. E. Koonin, “The disassembly of nuclear matter,” Nuclear Physics A, vol. 356, no. 1, pp. 223–234, 1981. doi: 10.1016/0375-9474(81)90124-X. [141] S. Albergo, S. Costa, E. Costanzo, and A. Rubbino, “Temperature and free-nucleon densities of nuclear matter exploding into light clusters in heavy-ion collisions,” Il Nuovo Cimento A (1965-1970), vol. 89, no. 1, pp. 1–28, 1985. doi: 10.1007/BF02773614. [142] N. Wang, Z. Li, and X. Wu, “Improved quantum molecular dynamics model and its applications to fusion reaction near barrier,” Physical Review C, vol. 65, no. 6, p. 064608, 2002. doi: 10.1103/PhysRevC.65.064608. [143] N. Wang, Z. Li, X. Wu, J. Tian, Y. Zhang, and M. Liu, “Further development of the improved quantum molecular dynamics model and its application to fusion reactions near the barrier,” Physical Review C, vol. 69, no. 3, p. 034608, 2004. doi: 10.1103/PhysRevC.69.034608. [144] Z. Chajecki, M. Youngs, D. D. S. Coupland, W. G. Lynch, M. B. Tsang, D. Brown, A. Chbihi, P. Danielewicz, R. T. Desouza, M. A. Famiano et al., “Scaling properties of light-cluster production,” arXiv preprint arXiv:1402.5216, 2014. doi: 10.48550/arXiv.1402.5216. [145] M. Kaneko, T. Murakami, T. Isobe, M. Kurata-Nishimura, A. Ono, N. Ikeno, J. Bar- ney, G. Cerizza, J. Estee, G. Jhang et al., “Rapidity distributions of Z = 1 iso- topes and the nuclear symmetry energy from Sn+ Sn collisions with radioactive beams at 270 MeV/nucleon,” Physics Letters B, vol. 822, p. 136681, 2021. doi: 10.1016/j.physletb.2021.136681. [146] C. E. Rasmussen, C. K. I. Williams et al., Gaussian processes for machine learning. Springer, 2006, vol. 1. ISBN 026218253X [147] P. Morfouace, C. Y. Tsang, Y. Zhang, W. G. Lynch, M. B. Tsang, D. D. S. Coupland, M. Youngs, Z. Chajecki, M. A. Famiano, T. K. Ghosh et al., “Constraining the symmetry energy with heavy-ion collisions and bayesian analyses,” Physics Letters B, vol. 799, p. 135045, 2019. doi: 10.1016/j.physletb.2019.135045. [148] B. Acharya and S. Bacca, “Gaussian process error modeling for chiral effective-field- theory calculations of np ↔ dγ at low energies,” Physics Letters B, vol. 827, p. 137011, 2022. doi: 10.1016/j.physletb.2022.137011. 176 [149] J. Keller, K. Hebeler, and A. Schwenk, “Nuclear equation of state for arbitrary pro- ton fraction and temperature based on chiral effective field theory and a gaussian process emulator,” Physical Review Letters, vol. 130, no. 7, p. 072701, 2023. doi: 10.1103/PhysRevLett.130.072701. [150] M. D. McKay, R. J. Beckman, and W. J. Conover, “A comparison of three methods for selecting values of input variables in the analysis of output from a computer code,” Technometrics: American Statistical Association, vol. 21, no. 2, pp. 239–245, 1979. doi: 10.2307/1268522. [151] J. Salvatier, T. V. Wiecki, and C. Fonnesbeck, “Probabilistic programming in Python using PyMC3,” PeerJ Computer Science, vol. 2, p. e55, 2016. doi: 10.7717/peerj-cs.55. [152] D. Coupland, “Probing the nuclear symmetry energy with heavy ion collisions,” Ph.D. dissertation, Michigan State University, 2013, dissertation. [153] R. Killick, P. Fearnhead, and I. A. Eckley, “Optimal detection of changepoints with a linear computational cost,” Journal of the American Statistical Association, vol. 107, no. 500, pp. 1590–1598, 2012. doi: 10.1080/01621459.2012.737745. [154] C. Truong, L. Oudre, and N. Vayatis, “Selective review of offline change point detection methods,” Signal Processing, vol. 167, p. 107299, 2020. doi: 10.1016/j.sigpro.2019.107299. [155] C. Truong, L. Oudre, and N. Vayatis, “ruptures: change point detection in Python,” arXiv preprint arXiv:1801.00826, 2018. doi: 10.48550/arXiv.1801.00826. [156] FARO Technologies Inc., “FARO Vantage Laser Trackers,” accessed: 2023-05-01. [Online]. Available: https://www.faro.com/en/Products/Hardware/ Vantage-Laser-Trackers [157] Autodesk, Inc., “Autodesk Inventor,” accessed: 2023-05-01. [Online]. Available: https://www.autodesk.com/products/inventor [158] A. G. Wright, The photomultiplier handbook. Oxford University Press, 2017. [159] A. Artikov, J. Budagov, I. Chirikov-Zorin, D. Chokheli, M. Lyablin, G. Bellettini, A. Menzione, S. Tokar, N. Giokaris, and A. Manousakis-Katsikakis, “Properties of the Ukraine polystyrene-based plastic scintillator UPS 923A,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and As- sociated Equipment, vol. 555, no. 1-2, pp. 125–131, 2005. doi: 10.1016/j.nima.2005.09.021. [160] Eljen Technology, “Neutron / Gamma PSD EJ-301, EJ-309,” 2021, accessed: 2023- 05-01. [Online]. Available: https://eljentechnology.com/products/liquid-scintillators/ ej-301-ej-309 177 [161] E. R. Siciliano, J. H. Ely, R. T. Kouzes, J. E. Schweppe, D. M. Strachan, and S. T. Yokuda, “Energy calibration of gamma spectra in plastic scintillators using Compton kinematics,” Nuclear Instruments and Methods in Physics Research Section A: Acceler- ators, Spectrometers, Detectors and Associated Equipment, vol. 594, no. 2, pp. 232–243, 2008. doi: 10.1016/j.nima.2008.06.031. [162] Z. Liu, J. Chen, P. Zhu, Y. Li, and G. Zhang, “The 4.438 MeV gamma to neutron ratio for the Am–Be neutron source,” Applied Radiation and Isotopes, vol. 65, no. 12, pp. 1318–1321, 2007. doi: 10.1016/j.apradiso.2007.04.007. [163] Z.-H. Yang, Y.-L. Ye, J. Xiao, H.-B. You, H.-N. Liu, Y.-L. Sun, Z.-H. Wang, and J. Chen, “Fast calibration methods using cosmic rays for a neutron detection array,” Chinese Physics C, vol. 36, no. 3, p. 222, 2012. doi: 10.1088/1674-1137/36/3/006. [164] I. Tilquin, Y. El Masri, M. Parlog, P. Collon, M. Hadri, T. Keutgen, J. Lehmann, P. Leleux, P. Lipnik, A. Ninane et al., “Detection efficiency of the neutron modular detector DEMON and related characteristics,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 365, no. 2-3, pp. 446–461, 1995. doi: 10.1016/0168-9002(95)00425-4. [165] “Physical constants,” Physics Letters B, vol. 204, pp. IN3–127, 1988. doi: 10.1016/0370- 2693(88)90508-4. [166] B. C. Rastin, “An accurate measurement of the sea-level muon spectrum within the range 4 to 3000 GeV/c,” Journal of Physics G: Nuclear Physics, vol. 10, no. 11, p. 1609, 1984. doi: 10.1088/0305-4616/10/11/017. [167] C. Hagmann, D. Lange, and D. Wright, “Cosmic-ray shower generator (CRY) for Monte Carlo transport codes,” in 2007 IEEE Nuclear Science Symposium Conference Record, vol. 2. IEEE, 2007. doi: 10.1109/NSSMIC.2007.4437209. pp. 1143–1146. [168] S. Agostinelli, J. Allison, K. a. Amako, J. Apostolakis, H. Araujo, P. Arce, M. Asai, D. Axen, S. Banerjee, G. Barrand et al., “GEANT4—a simulation toolkit,” Nu- clear instruments and methods in physics research section A: Accelerators, Spectrom- eters, Detectors and Associated Equipment, vol. 506, no. 3, pp. 250–303, 2003. doi: 10.1016/S0168-9002(03)01368-8. [169] J. Allison, K. Amako, J. Apostolakis, P. Arce, M. Asai, T. Aso, E. Bagli, A. Bagulya, S. Banerjee, G. Barrand et al., “Recent developments in Geant4,” Nuclear instruments and methods in physics research section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 835, pp. 186–225, 2016. doi: 10.1016/j.nima.2016.06.125. [170] M. Ahmed, “A comparative study of n-γ discrimination properties of scintillators NE 213, C6H6, C6D6 and stilbene,” Nuclear Instruments and Methods, vol. 143, no. 2, pp. 255–257, 1977. doi: 10.1016/0029-554X(77)90603-6. [171] L. J. Perkins and M. C. Scott, “The application of pulse shape discrimination in NE 213 to neutron spectrometry,” Nuclear Instruments and Methods, vol. 166, no. 3, pp. 451–464, 1979. doi: 10.1016/0029-554X(79)90534-2. 178 [172] S. Green, R. Koohi-Fayegh, and M. C. Scott, “A gated single parameter data acquisition system for fast neutron spectroscopy with NE-213,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 317, no. 1-2, pp. 399–401, 1992. doi: 10.1016/0168-9002(92)90635-H. [173] G. F. Knoll, Radiation detection and measurement, 4th ed. John Wiley & Sons, 2010. ISBN 978-0-470-13148-0 [174] T. Kögler, R. Beyer et al., “Light yield and n–γ pulse-shape discrimination of liquid scintillators based on linear alkyl benzene,” Nucl. Instrum. Methods Phys. Res. A, vol. 701, pp. 285–293, 2013. doi: 10.1016/j.nima.2012.10.059. [175] H. Sakai, H. Okamura, S. Ishida, K. Hatanaka, and T. Noro, “Construction and performance of one-and two-dimensional large position-sensitive liquid and plastic scintillation detectors—an application to a neutron polarimeter,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 320, no. 3, pp. 479–499, 1992. doi: 10.1016/0168- 9002(92)90944-Y. [176] U. Bravar, R. S. Woolf, P. J. Bruillard, E. O. Fluckiger, J. S. Legere, A. L. MacKinnon, J. R. Macri, P. C. Mallik, M. L. McConnell, B. Pirard et al., “Calibration of the fast neutron imaging telescope (FNIT) prototype detector,” IEEE Transactions on Nuclear Science, vol. 56, no. 5, pp. 2947–2954, 2009. doi: 10.1109/TNS.2009.2028025. [177] J. Heideman, D. Pérez-Loureiro, R. Grzywacz, C. R. Thornsberry, J. Chan, L. H. Heilbronn, S. K. Neupane, K. Schmitt, M. M. Rajabali, A. R. Engelhardt et al., “Con- ceptual design and first results for a neutron detector with interaction localization capabilities,” Nuclear Instruments and Methods in Physics Research Section A: Accel- erators, Spectrometers, Detectors and Associated Equipment, vol. 946, p. 162528, 2019. doi: 10.1016/j.nima.2019.162528. [178] L. Stuhl, M. Sasano, K. Yako, J. Yasuda, H. Baba, S. Ota, and T. Uesaka, “PAN- DORA, a large volume low-energy neutron detector with real-time neutron–gamma discrimination,” Nuclear Instruments and Methods in Physics Research Section A: Ac- celerators, Spectrometers, Detectors and Associated Equipment, vol. 866, pp. 164–171, 2017. doi: 10.1016/j.nima.2017.06.015. [179] C. L. Morris, J. E. Bolger et al., “A digital technique for neutron-gamma pulse shape discrimination,” Nucl. Instrum. and Methods, vol. 137, no. 2, pp. 397–398, 1976. doi: 10.1016/0029-554X(76)90353-0. [180] E. M. Ellis, C. Hurlbut, C. Allwork, and B. Morris, “Neutron and gamma ray pulse shape discrimination with EJ-270 lithium-loaded plastic scintillator,” in 2017 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC). IEEE, 2017. doi: 10.1109/NSSMIC.2017.8532688. pp. 1–5. [181] A. A. Ivanova, P. V. Zubarev, S. V. Ivanenko, A. D. Khilchenko, A. I. Kotelnikov, S. V. Polosatkin, E. A. Puryga, V. G. Shvyrev, and Y. S. Sulyaev, “Fast neutron flux 179 analyzer with real-time digital pulse shape discrimination,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 827, pp. 13–17, 2016. doi: 10.1016/j.nima.2016.04.088. [182] A. Jančář, Z. Kopecký et al., “Pulse-shape discrimination of the new plastic scintilla- tors in neutron–gamma mixed field using fast digitizer card,” Radiation Physics and Chemistry, vol. 116, pp. 60–64, 2015. doi: 10.1016/j.radphyschem.2015.05.007. [183] S. MacMullin, M. Kidd et al., “Measurement of the elastic scattering cross section of neutrons from argon and neon,” Physical Review C, vol. 87, no. 5, p. 054613, 2013. doi: 10.1103/PhysRevC.87.054613. [184] R. H. Showalter, “Determination of density and momentum dependence of nuclear symmetry potentials with asymmetric heavy ion reactions,” Ph.D. dissertation, Michigan State University, 2015, dissertation. [185] D. S. McGregor and J. K. Shultis, “Reporting detection efficiency for semiconductor neutron detectors: A need for a standard,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 632, no. 1, pp. 167–174, 2011. doi: 10.1016/j.nima.2010.12.084. [186] B. Luther, T. Baumann, M. Thoennessen, J. Brown, P. DeYoung, J. Finck, J. Hin- nefeld, R. Howes, K. Kemper, P. Pancella et al., “Mona—the modular neutron array,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spec- trometers, Detectors and Associated Equipment, vol. 505, no. 1-2, pp. 33–35, 2003. doi: 10.1016/S0168-9002(03)01014-3. [187] G. Jaworski, M. Palacz, J. Nyberg, G. De Angelis, G. De France, A. Di Nitto, J. Egea, M. N. Erduran, S. Ertürk, E. Farnea et al., “Monte Carlo simulation of a single detector unit for the neutron detector array NEDA,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 673, pp. 64–72, 2012. doi: 10.1016/j.nima.2012.01.017. [188] J. Park, F. C. E. Teh, M. B. Tsang, K. W. Brown, Z. Chajecki, B. Hong, T. Lokotko, W. G. Lynch, D. Satoh, J. Wieske, and K. Zhu, “neuSIM4: A GEANT4 based comprehensive neutron simulation code,” 2023, manuscript in preparation. [189] G. Beam, L. Wielopolski, R. Gardner, and K. Verghese, “Monte Carlo calculation of efficiencies of right-circular cylindrical NaI detectors for arbitrarily located point sources,” Nuclear instruments and methods, vol. 154, no. 3, pp. 501–508, 1978. doi: 10.1016/0029-554X(78)90081-2. [190] L. Moens, J. De Donder, L. Xi-Lei, F. De Corte, A. De Wispelaere, A. Simonits, and J. Hoste, “Calculation of the absolute peak efficiency of gamma-ray detectors for different counting geometries,” Nuclear Instruments and Methods in Physics Research, vol. 187, no. 2-3, pp. 451–472, 1981. doi: 10.1016/0029-554X(81)90374-8. 180 [191] P. Van Otten, R. Van de Vyver, E. Van Camp, E. Kerkhove, P. Berkvens, H. Ferdinande, D. Ryckbosch, A. De Graeve, and L. Van Hoorebeke, “A general Monte Carlo calculation for the geometrical efficiency of a detection system,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 267, no. 1, pp. 183–192, 1988. doi: 10.1016/0168-9002(88)90646-8. [192] D. Satoh, T. Sato, N. Shigyo, and K. Ishibashi, “SCINFUL-QMD: Monte Carlo based computer code to calculate response function and detection efficiency of a liquid organic scintillator for neutron energies up to 3 GeV,” JAEA-Data/Code, vol. 23, p. 2006, 2006. doi: 10.11484/jaea-data-code-2006-023. [193] J. K. Dickens, “SCINFUL: A Monte Carlo based computer program to determine a scintillator full energy response to neutron detection for En between 0.1 and 80 MeV: Program development and comparisons of program predictions with experimental data,” Oak Ridge National Lab., TN (USA), Tech. Rep., 1988. [194] A. Boudard, J. Cugnon, J.-C. David, S. Leray, and D. Mancusi, “New potentialities of the Liège intranuclear cascade model for reactions induced by nucleons and light charged particles,” Physical Review C, vol. 87, no. 1, p. 014606, 2013. doi: 10.1103/Phys- RevC.87.014606. [195] D. Mancusi, A. Boudard, J. Cugnon, J.-C. David, P. Kaitaniemi, and S. Leray, “Ex- tension of the Liège intranuclear-cascade model to reactions induced by light nuclei,” Physical Review C, vol. 90, no. 5, p. 054602, 2014. doi: 10.1103/PhysRevC.90.054602. [196] D. Satoh and T. Sato, “Improvements in the particle and heavy-ion transport code system (PHITS) for simulating neutron-response functions and detection efficiencies of a liquid organic scintillator,” Journal of Nuclear Science and Technology, vol. 59, no. 8, pp. 1047–1060, 2022. doi: 10.1080/00223131.2021.2019622. [197] V. V. Verbinski, W. R. Burrus, T. A. Love, W. Zobel, N. W. Hill, and R. Textor, “Calibration of an organic scintillator for neutron spectrometry,” Nuclear Instruments and Methods, vol. 65, no. 1, pp. 8–25, 1968. doi: 10.1016/0029-554X(68)90003-7. [198] S. Meigo, “Measurements of the response function and the detection efficiency of an NE213 scintillator for neutrons between 20 and 65 MeV,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 401, no. 2-3, pp. 365–378, 1997. doi: 10.1016/S0168- 9002(97)01061-9. [199] N. Nakao, T. Kurosawa, T. Nakamura, and Y. Uwamino, “Absolute measurements of the response function of an NE213 organic liquid scintillator for the neutron energy range up to 206 MeV,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 463, no. 1-2, pp. 275–287, 2001. doi: 10.1016/S0168-9002(01)00260-1. [200] M. Sasaki, N. Nakao, T. Nakamura, T. Shibata, and A. Fukumura, “Measurements of the response functions of an NE213 organic liquid scintillator to neutrons up to 181 800MeV,” Nuclear Instruments and Methods in Physics Research Section A: Accelera- tors, Spectrometers, Detectors and Associated Equipment, vol. 480, no. 2-3, pp. 440–447, 2002. doi: 10.1016/S0168-9002(01)00948-2. [201] Z. Kohley, E. Lunderberg, P. A. DeYoung, B. Roeder, T. Baumann, G. Chris- tian, S. Mosby, J. K. Smith, J. Snyder, A. Spyrou et al., “Modeling interactions of intermediate-energy neutrons in a plastic scintillator array with GEANT4,” Nu- clear Instruments and Methods in Physics Research Section A: Accelerators, Spec- trometers, Detectors and Associated Equipment, vol. 682, pp. 59–65, 2012. doi: 10.1016/j.nima.2012.04.060. [202] C. Y. Tsang, M. Kurata-Nishimura, M. B. Tsang, W. G. Lynch, Y. X. Zhang, J. Barney, J. Estee, G. Jhang, R. Wang, M. Kaneko et al., “Constraining nucleon effective masses with flow and stopping observables from the SπRIT experiment,” 2023, manuscript in preparation. [203] S. Voloshin and Y. Zhang, “Flow study in relativistic nuclear collisions by Fourier expansion of azimuthal particle distributions,” Zeitschrift für Physik C Particles and Fields, vol. 70, pp. 665–671, 1996. doi: 10.1007/s002880050141. [204] Reisdorf, W and Andronic, A and Gobbi, A and Hartmann, ON and Herrmann, N and Hildenbrand, KD and Kim, YJ and Kirejczyk, M and Koczoń, P and Kress, T and others, “Nuclear stopping from 0.09 A to 1.93 A GeV and its correlation to flow,” Physical review letters, vol. 92, no. 23, p. 232301, 2004. doi: 10.1103/PhysRevLett.92.232301. [205] C. Xu, B.-A. Li, and L.-W. Chen, “Symmetry energy, its density slope, and neutron- proton effective mass splitting at normal density extracted from global nucleon optical potentials,” Physical Review C, vol. 82, no. 5, p. 054607, 2010. doi: 10.1103/Phys- RevC.82.054607. [206] R. J. Charity, W. H. Dickhoff, L. G. Sobotka, and S. J. Waldecker, “Isospin dependence of nucleon correlations in ground-state nuclei,” The European Physical Journal A, vol. 50, pp. 1–7, 2014. doi: 10.1140/epja/i2014-14023-0. [207] X.-H. Li, W.-J. Guo, B.-A. Li, L.-W. Chen, F. J. Fattoyev, and W. G. Newton, “Neutron–proton effective mass splitting in neutron-rich matter at normal density from analyzing nucleon–nucleus scattering data within an isospin dependent optical model,” Physics Letters B, vol. 743, pp. 408–414, 2015. doi: 10.1016/j.physletb.2015.03.005. [208] Z. Zhang and L.-W. Chen, “Isospin splitting of the nucleon effective mass from gi- ant resonances in 208 Pb,” Physical Review C, vol. 93, no. 3, p. 034335, 2016. doi: 10.1103/PhysRevC.93.034335. [209] H.-Y. Kong, J. Xu, L.-W. Chen, B.-A. Li, and Y.-G. Ma, “Constraining simultaneously nuclear symmetry energy and neutron-proton effective mass splitting with nucleus giant resonances using a dynamical approach,” Physical Review C, vol. 95, no. 3, p. 034324, 2017. doi: 10.1103/PhysRevC.95.034324. 182 [210] J. Xu, W.-J. Xie, and B.-A. Li, “Bayesian inference of nuclear symmetry energy from measured and imagined neutron skin thickness in 116,118,120,122,124,130,132 , 208 Pb, and 48 Ca,” Physical Review C, vol. 102, no. 4, p. 044316, 2020. doi: 10.1103/PhysRevC.102.044316. [211] J. Xu and W.-T. Qin, “Nucleus giant resonances from an improved isospin-dependent Boltzmann-Uehling-Uhlenbeck transport approach,” Physical Review C, vol. 102, no. 2, p. 024306, 2020. doi: 10.1103/PhysRevC.102.024306. [212] P. T. P. Hutauruk, H. Gil, S.-i. Nam, and C. H. Hyun, “Effect of nucleon effective mass and symmetry energy on the neutrino mean free path in a neutron star,” Physical Review C, vol. 106, no. 3, p. 035802, 2022. doi: 10.1103/PhysRevC.106.035802. [213] E. R. Most and C. A. Raithel, “Impact of the nuclear symmetry energy on the post- merger phase of a binary neutron star coalescence,” Physical Review D, vol. 104, no. 12, p. 124012, 2021. doi: 10.1103/PhysRevD.104.124012. [214] R. Cavagnoli, D. P. Menezes, and C. Providencia, “Neutron star properties and the symmetry energy,” Physical Review C, vol. 84, no. 6, p. 065810, 2011. doi: 10.1103/Phys- RevC.84.065810. [215] J. Xu, L.-W. Chen, B.-A. Li, and H.-R. Ma, “Nuclear constraints on properties of neutron star crusts,” The Astrophysical Journal, vol. 697, no. 2, p. 1549, 2009. doi: 10.1088/0004-637X/697/2/1549. 183 APPENDIX A SPECTRAL YIELDS FOR LIGHT CHARGED PARTICLES This appendix is dedicated to the presentation of two-dimensional spectra yield plots, organized as θlab -Elab and pT /A-y/ybeam histograms. The plots encompass results from six types of charged particles, namely, protons (p), deuterons (d), tritons (t), 3He, 4He, and 6He, observed in each examined pair of reaction systems at a particular beam energy. It should be noted that while 6He is included in these histograms for completeness of the helium isotopes set, it was not incorporated in the analysis of neutron to proton (n/p) ratios and isoscaling due to its low statistical presence. On comparing the spectra across different particles within the same reaction system, similarities in spectral patterns are noted for the Ca + Ni and Ca + Sn systems. Beam energy, however, serves as the biggest differentiating factor among the spectra. Finally, we must bring to attention an observable gap at θlab ≈ 64◦ in the histograms of 48 the Ca + 64Ni and 48 Ca + 64Sn systems at 140 MeV/u. This data void, which unfortunately is irretrievable, is due to all the pixels on the HiRA10 detector at those particular angles being inactive during those measurements. 184 Ca + Ni systems at 56 MeV/u: Figure A.1 Two-dimensional θlab -Elab plot for Ca + Ni systems at 56 MeV/u. 185 Figure A.2 Two-dimensional pT /A-y/ylab plot for Ca + Ni systems at 56 MeV/u. 186 Ca + Ni systems at 140 MeV/u: Figure A.3 Two-dimensional θlab -Elab plot for Ca + Ni systems at 140 MeV/u. 187 Figure A.4 Two-dimensional pT /A-y/ylab plot for Ca + Ni systems at 140 MeV/u. 188 Ca + Sn systems at 56 MeV/u: Figure A.5 Two-dimensional θlab -Elab plot for Ca + Sn systems at 56 MeV/u. 189 Figure A.6 Two-dimensional pT /A-y/ylab plot for Ca + Sn systems at 56 MeV/u. 190 Ca + Sn systems at 140 MeV/u: Figure A.7 Two-dimensional θlab -Elab plot for Ca + Sn systems at 140 MeV/u. 191 Figure A.8 Two-dimensional pT /A-y/ylab plot for Ca + Sn systems at 140 MeV/u. 192 APPENDIX B TRANSVERSE MOMENTUM VERSUS RAPIDITY This appendix aims to derive two relations involving transverse momentum and rapidity. These relations outline the phase space boundaries due to energy and θ cuts. Transverse momentum is given as pT = p sin θ , (B.1) where p is the magnitude of the momentum vector and θ is the angle between the momentum vector and the beam direction. On the other hand, the (experimentally modified) rapidity is defined as   1 E + p cos θ y = ln , (B.2) 2 E − p cos θ p where E = p2 + m2 is the total energy of the particle and m is the mass of the particle. To derive the energy-cut boundary, we first rewrite Equation B.2 into p cos θ = E tanh y . (B.3) Then, using Equation B.1 and Equation B.3, we obtain (p cos θ)2 + (p sin θ)2 = (E tanh y)2 + p2T (B.4) q ∴ E = p2T + m2 · cosh y . (B.5) Note that this expression of E is independent of θ, and thus describes the energy-cut boundary for a given value of y. To derive the θ-cut boundary, we take the ratio of Equation B.1 and Equation B.3: p sin θ pT = (B.6) p cos θ E · tanh y ! pT ∴ θ = arctan p . (B.7) p2T + m2 · sinh y Note that this expression of θ is independent of E, and thus describes the θ-cut boundary for a given value of θ. 193