Nuclear liquid-gas phase transition with machine learning

The machine-learning techniques have shown their capability for studying phase transitions in condensed matter physics. Here, we employ the machine-learning techniques to study the nuclear liquid-gas phase transition. We adopt an unsupervised learning and classify the liquid and gas phases of nuclei directly from the final state raw experimental data of heavy-ion reactions. Based on a confusion scheme which combines the supervised and unsupervised learning, we obtain the limiting temperature of the nuclear liquid-gas phase transition. Its value $9.24\pm0.04~\rm MeV$ is consistent with that obtained by the traditional caloric curve method. Our study explores the paradigm of combining the machine-learning techniques with heavy-ion experimental data, and it is also instructive for studying the phase transition of other uncontrollable systems, like QCD matter.


I. INTRODUCTION
The nuclear liquid-gas phase transition is an old and long-last topic [1][2][3][4][5][6]. Since the interaction between nucleons exhibit Van der Waals features similar with that between molecules, i.e., a short-distance repulsive core and a long-distance attractive tail, the nuclei, considered as self-bound Fermi liquid, can experience liquid-gas phase transition as well. Over the past several decades, the nuclear liquid-gas phase transition has been studied based on the heavy-ion collisions at intermediate and relativistic energies and hadron-nucleus collisions at relativistic energies. The information of the reaction products are obtained with powerful multidetectors allowing the detection of a large amount of the fragments and light particles produced in the reaction. A lot of probes by analyzing sophisticatedly the information of the reaction products have been proposed to recognize the liquid-gas phase transition of nuclei [7][8][9][10][11][12][13][14][15][16][17][18][19][20].
The ability of machine-learning techniques [21,22] of recognizing and characterizing complex sets of data stimulates their applications on physics, and brings new possibilities to the study of the nuclear liquid-gas phase transition. Besides the common uses like particle identification and tagging in experiments [23][24][25], machinelearning techniques have various novel applications. Several examples are solving quantum many-body problem [26], analyzing strong gravitational lenses [27], exploring phase properties of quark matter [28][29][30], constraining and studying field theories [31,32], and quantum state tomography [33]. Most notably, in condensed matter physics, machine-learning techniques have been used to classify phases of matter, and identify topological order and phase transitions [34][35][36][37].
The resemblance between condensed matter physics and nuclear physics is remarkable. They both involve large amount of degrees of freedom, and share the same theoretical tools like Hartree-Fock theory and finitetemperature field theory. Experimentally however, unlike the condensed matter physics, the nucleus is an uncontrollable system, and its thermal properties can only be accessed through nuclear reactions. Thus to treat the nuclear liquid-gas phase transition in terms of equilibrium thermodynamics is not practicable. The nuclear liquidgas phase transition is realized through tracing the effect of the spinodal instability, which is intimately related to first-order phase transition, on the reaction dynamics, e.g., by measuring the properties of the intermediate mass fragments (with charge number larger than 3). Heavy-ion reaction experiments may bring the excited nuclei into the spinodal region of the phase diagram in which the spinodal instability may develop exponentially and lead to the break-up of nuclei. This is commonly referred to as the nuclear multifragmentation. In Fig. 1, we arXiv:2010.15043v1 [nucl-th] 28 Oct 2020 draw schematically the phase diagram and typical phase trajectories of an excited nucleus in heavy-ion reactions. At the beginning of the reaction, the vast majority of a ground state nucleus is around the saturation density ρ 0 , which is approximately 0.16 nucleon fm −3 . After hitting the target nucleus, the projectile nucleus is excited and compressed, and is regarded as heated liquid (same for the target nucleus but we will focus on the projectile, which is easier to measure for the detectors). For low excitation energy, the compressed projectile nucleus will expand and then exhibit a damped monopole oscillation accompanied by the emissions of a few light particles, and remains as liquid. As the excitation energy increases, the expansion of the nucleus becomes severer and drives the nucleus into the spinodal region. In the spinodal region, due to the attractive part in the nucleon-nucleon interaction, the high density region will attract nucleons from low density region. This causes the formation of intermediate mass fragments. The spinodal decomposition process thus drives the system towards thermodynamic liquid-gas phase coexistence. If the excitation energy continues to increase, the excited nucleus will expand quickly enough to pass through the spinodal region. The formation of the intermediate mass fragments becomes less important, since dynamically it takes time for the fragments to format. In this case, the entire projectile nucleus ends up with dominant light fragments, which corresponds to a gas phase. Based on the above discussion, the observation of the intermediate mass fragments or nuclear multifragmentation is a sign of the nuclear liquid-gas phase transition.
In the present article, we want to demonstrate that the machine learning techniques are also capable of dealing with the phase transition in realistic nuclear system, other than theoretical models in condensed matter physics, though the way of studying the phase transition of the two systems exhibits essential differences. To that end, we train the neural network instead of, e.g., by the spin configurations from Monte Carlo simulations, but by the final state information of heavy-ion reaction experiment. We then use the trained networks to classify the liquid and gas phases of nuclei, and determine the limiting temperature of the nuclear liquid-gas phase transition.

II. EXPERIMENTAL DATA
The reactions of 40 Ar on 27 Al and 48 Ti at 47 MeV/nucleon were performed using beams from the TAMU K500 super-conducting cyclotron. The charge and momentum of charged fragments are probed with the 4π detector, NIMROD [38] (Neutron Ion Multidetector for Reaction Oriented Dynamics). Generally speaking, phase transitions of small system (about 10 constituents) are still well defined, distinguishable [39], and even detectable [40]. The spectator matter in the nuclear reaction is argued to be ideally suited to investigate the nuclear liquid-gas phase transition since it can largely avoid the effect of dynamical evolution. Because of the essential binary nature of such reactions, we applied a reconstruction method for the quasi-projectile (QP) source developed by Ma et al. [15], which perform a three-source (i.e., a QP source, an intermediate velocity source and a quasi-target source) fit for light particles with Z 3, and then use the probability of QP particles to identify the QP light fragments in an event-by-event basis. For heavier fragments, they can be assigned to the QP source directly through a rapidity cut. The QP fragments are supposed to come from the excited projectile nucleus. By the above QP-labeled light and heavier fragments, the mass and charge numbers of QP source could be reconstructed in each event and its excitation energy and other bulk properties can be obtained. Details of the Ma's QP reconstruction method can be found in Ref. [15]. Due to the existence of an intermediate velocity source, the total charge of QP fragments Z QP is generally less than the charge number of the projectile nucleus. We choose events with Z QP 12 as good QP events, and focus on those Z QP = 12 events since they have the largest statistics (using events with other Z QP does not change our conclusion qualitatively). The universality of the QP fragment distributions [15] indicates the memory of the entrance channel dynamics is lost prior to the decay of the excited spectators. Thus the final state information of the QP fragments of the above two reactions is combined to form a single event-by-event data set. This leads to 40081 valid QP events with Z QP = 12. Physically, the excitation energy per nucleon E ex /A and apparent temperature T ap of the QP nucleus are used to characterize each QP event. E ex is deduced event-byevent through [8,15]  The three terms represent the kinetic energy of the charged QP fragments and neutrons, and the mass excess, respectively. The apparent temperature of the QP nucleus can be obtained by measuring the light particles evaporated from its surface. This can be achieved via different thermometers, e.g., particle kinetic energy or isotope yield ratios [41]. In the present work, the quadruple momentum fluctuation [42] is used as the thermometer.
The quadruple momentum is defined as Q xy = p 2 x − p 2 y , where p x and p y are the transverse components of the emitted particle momentum in lab frame. Since the apparent temperature is derived from the fluctuation of Q xy , determining the reaction plane is not necessary. When the momenta distribute in a Maxwellian form, the average temperature of the events in a given E ex /A bin is related to the variance of Q xy , i.e., where m represents the mass of the probe particle (deuteron in the present work). The event-by-event T ap is obtained through a Monte Carlo method based on the standard deviation of T ap (details can be found in the appendix). We draw the scatter plot of T ap verse E ex /A, or the caloric curve, of the events with Z QP = 12 in Fig. 2. The red dashed curve in the figure represents the T ap as a function of E ex /A.

III. RESULTS
The power of the machine-learning techniques lies in their ability to classify the phases of matter prior to the knowledge of the characteristic quantities, namely, E ex /A and T ap . In that sense, the event-by-event chargeweighted charge multiplicity distribution of QP fragments ZM c (Z) from experiment is used as the input to train the neural network (the charge weighting is for normalization). Among the 40081 valid QP events with Z QP = 12, 2/3 are used for training and others for testing. To obtain an intuitive impression of M c (Z), we average for testing events their M c (Z) in three typical E ex /A bins, namely, low excitation (0.9 -2. We first adopt an unsupervised learning, the autoen-coder method [43], to study the nuclear liquid-gas phase transition. We show the construction of the autoencoder network used in the present work in Fig. 4. The neural network consists of two main parts, the encoder part encodes the inputted event-by-event ZM c (Z) to a latent variable (or code [44]), and the decoder part decodes the latent variable to ZM c (Z), and tries to restore the original ZM c (Z). The neural network is trained to best restore the encoded information, which means the network is trained to minimize the difference between ZM c (Z) and ZM c (Z). There are two layers in the encoder and decoder parts respectively, and all layer are fully connected, more details can be found in the appendix. For the testing QP events, the reconstructed M c (Z) are averaged and compared with the original M c (Z) in Fig. 3. We also show the mean and standard deviation of the event-by-event reconstruction loss in the appendix. We notice from Fig. 3 that the autoencoder network succeeds in capturing essential information of the inputted event-by-event ZM c (Z). Once we finish training the autoencoder network, through the charge multiplicity distribution, each QP event is mapped to a number (latent variable). We plot in Fig. 5 the latent variable as a function of T ap and E ex /A, with each data point averaged over 500 testing events. The vertical error bars in the figure represent the standard deviations of the latent variable of these events, while the horizontal error bars the standard deviations of characteristic parameter (T ap or E ex /A). Although there are large errors due to the event-by-event fluctuations of the experimental charge multiplicity distribution, the averaged latent variable as a function of T ap or E ex /A exhibits a sigmoid pattern, which indicates the trained autoencoder network treats the low and high temperature (or low and high excitation energy) regions differently. Considering that the autoencoder network is trained prior to any knowledge of the characteristic quantities, i.e., E ex /A and T ap , the autoencoder network is capable of classifying different phases of nuclei directly from the final state information of the heavy-ion experi-ment. The area in the midst of the two phases represents those liquid-gas coexistence events that enter the spinodal region and affected by the spinodal instability. It is interesting for further studies to find out how the latent variable is related to physical quantities.

B. Limiting temperature from confusion scheme
Traditionally, the nuclear liquid-gas phase transition can be recognized from the relation between E ex /A and the apparent temperature T ap , namely, the caloric curve [7]. As the excitation energy increases, more energy is transferred to the internal energy of the liquid phase projectile nucleus, and the apparent temperature increases. Dramatic change happens if the excited projectile nucleus goes into the spinodal region. Part of the excitation energy is consumed to form the fragments, in other words, transfer to latent heat. As a consequence, the increase of the apparent temperature slows down significantly, and even a plateau in the caloric curve is observed [7,13]. The specific heat capacity of the collision processc is defined to describe quantitatively the effect of the spinodal instability, i.e., which can be obtained from the caloric curve shown in Fig. 2. Note its difference with c P and c V since the external conditions on pressure and volume is not practicable due to the complex nature of such finite, self-bound object like nucleus. The specific heat capacityc will reach a peak when the spinodal instability affects the excited projectile nucleus the most severely. The corresponding apparent temperature at the maximum is called limiting temperature, which can be used to deduce the critical temperature of infinite nuclear matter [12]. Although the limiting temperature is different from the first-order critical temperature of isobaric process, the non-monotonic structure ofc(T ap ) indicates the existence of the spinodal region in nuclear matter phase diagram, which is a convincing evidence of the existence of nuclear liquidgas (first-order) phase transition. The above picture of the nuclear liquid-gas phase transition is consistent with the feature of the latent variable shown in Fig. 5, and we can further obtain the limiting temperature by a confusion scheme [35]. In the confusion scheme, the neural network is trained with data that are deliberately labelled incorrectly according to a proposed critical point, and the phase transition properties can be deduced from the performance curve, i.e., the total testing accuracy as a function of the proposed critical point, of the neural network [35]. In the original confusion scheme [35], a W-shape performance curve of the Ising model is observed, and the proposed critical point corresponding to the local maximum in the middle of the curve is recognized as the realistic critical point. As mentioned above, the nucleus is an uncontrollable system, thus the events used to train the neural network can not come from separate phases like those in the Ising model. The change from liquid to gas phase is gradual through a liquid-gas phase coexistence. Therefore as we will demonstrate below, a W-shape performance curve in the case of the Ising model, is replaced by a V-shape, when adopting the confusion scheme in the nuclear liquid-gas phase transition. Taking the example of T ap , the picture is as follows. In order to properly include the uncertainty of the obtained limiting temperature, we construct a Bayesian neural network (BNN) which contains two fully connected layers, as shown in Fig. 6, to perform this supervised classification. The QP events are divided into two categories and labelled as liquid-like or gas-like according to a proposed transition temperature T ap . When T ap increases from low temperature (vice versa when T ap decreases from high temperature), the feature of the liquid phase (a large fragment accompanied by several small fragments) emerges on both the liquid-like and gas-like category. As a consequence, the testing accuracy of low temperature region below T ap is relatively low, and the total performance of the network P (T ap ) starts to decrease, as shown in the left panel of Fig. 7. When T ap continues to increase, the features of the intermediate mass fragments begin to show up in the liquid-like category and the total performance continues to decrease. As T ap approaches the realistic limiting temperature, the features of the intermediate mass fragments become evident in both the low and high temperature categories, thus significantly reduces the total efficiency of the neural network. The total testing accuracy then reaches its minimum at T ap ≈ T lim . The error bars in the figure are the standard deviation of the testing accuracy, and are obtained based on the trained BNN by performing 10 test runs, with each consists 1000 random selected testing events. The limiting temperature is obtained by a parabolic fit of the lowest five data points with errors. The limiting temperature through the confusion scheme is 9.24 ± 0.04 MeV, which is consistent with the 9.0 ± 0.4 MeV obtained from the traditional analysis of caloric curve [45]. Besides that, it also locates in the intermediate temperature region in the left panel of Fig. 5 (cyan vertical dashed line), which indicates both the autoencoder network and confusion scheme learn the basic feature of liquid and gas phase from the raw event-by-event charge multiplicity distribution.
In the right panel of Fig. 7, similar analysis is carried out on E ex /A, and we obtain an analogical characteristic value of E ex /A = 5.79 ± 0.02 MeV. Since the network classifies each event purely based on ZM c (Z), the liquidlike events and gas-like events divided by the characteristic value of T ap (horizontal dotted cyan line in Fig. 2) correspond to almost same events that divided by the characteristic value of E ex /A (vertical dotted cyan line in Fig. 2). Considering that the latent variable obtained in Sec. III A is correlated with T ap or E ex /A, performing the above analysis with the latent variable will leads to similar result.
The performance curve P (T ap ) and P (E ex /A), i.e., the testing accuracy as a function of the proposed temperature T ap and transition excitation energy E ex , respectively. The yellow dotted lines represent a parabolic fit of the lowest five data points with errors.

IV. SUMMARY AND OUTLOOK
We have shown that the machine-learning techniques can be employed to a traditional nuclear physics topic, the nuclear liquid-gas phase transition. Based on the experiment event-by-event charge multiplicity distribution, the neural networks are capable of classifying the liquid and gas phases, and determining the limiting temperature of the nuclear liquid-gas phase transition. The hidden parameters identified by machine learning may acquire physical significance. The latent variable can be related to the order parameters for certain systems [46,47]. A new field called 'softness', which characterizes the local structure, is used to study the correlations between structure and dynamics in glassy liquids [48]. To relate the latent variable in Fig. 5 to certain physical quantities might be meaningful for future studies. The analysis performed here can also be applied to study the first-order phase transition of the QCD matter by choosing proper final state observables, since the picture of the first-order phase transition of the QCD matter is quite similar with that of the nuclear matter [49]. We anticipate more sophisticated observables like the kinetic energy spectra of the final state particles, along with more advanced machine-learning techniques will provide us new features in the nuclear liquid-gas phase transition, even in other fields of nuclear physics, that beyond the present knowledge.

ACKNOWLEDGMENTS
We thank Meisen Gao, Jie Pu and Ying Zhou for the maintenance of the GPU severs, Xiaopeng Zhang for technical support, and Feng Li for useful discussion. This work is partially supported by the National Natural Science Foundation of China under Contracts No. 11890714, No. 11421505  The event-by-event excitation energy is obtained through Eq. (1). In Eq. (1), E kin i represents the kinetic energy of the i-th charged QP fragment, which is obtained through the measured data directly. The contribution of the neutron kinetic energy is taken as 3/2M n T . Since the NIMROD does not provide the momentum of neutrons, the associated QP neutrons can not be separated through the three source fit. Therefore the neutron multiplicity M n is approximated as the difference between the assumed total QP mass and the sum of the detected masses of the QP fragments, i.e., M n = A QP − A QP i . The total mass of the QP spectator A QP is determined from Z QP (the sum of Z QP i , with Z QP i being the charge of the i-th measured fragments), by assuming the QP spectator has the same neutron-proton ratio as the initial projectile. The temperature T is assumed to equal to the temperature of QP protons which is taken from the three source fit parameters. The non-observed light charged particles from the QP spectator are calculated from the extracted three source fit parameters and added to the QP in mass and energy.
Appendix B: Event-by-event apparent temperature The apparent temperature of the QP nucleus is obtained via the quadrupole momentum fluctuation. For the QP events in a given excitation energy bin, the experimental deuteron quadrupole momentum fluctuation temperature T ap is obtained through Eq. (2), and its standard deviation values ∆ T ap are evaluated by the TProfile class of ROOT in the CERN data analysis library [50]. Both T ap and ∆ T ap are fitted by poly-nomial functions to obtain their E ex /A dependence, i.e., T pol (E ex /A) and ∆T pol (E ex /A). Since the standard deviation ∆ T ap is small, the event-by-event apparent temperature of the reconstructed QP nucleus is evaluated through where N (0, 1) is a zero-mean Gaussian random number with unit variance.
Appendix C: Details of the neural network The neural network contains successively one input layer, several hidden layers, and one output layer. Each layer provides its output z through a matrix multiplication of its input x, i.e., z = W · x + b. The elements in the matrix W are known as weights and in the vector b as biases. In a normal full-connected neural network these parameters are single values, while in a Bayesian neural network (BNN) they become distributions, usually assumed to be Gaussian form. The layer is then followed by an activation function, f (z), which turns a linear transform to a non-linear one. Commonly used activation functions are sigmoid, tanh, and ReLU (rectified linear unit). f (z) is then used as the input of the next layer. The neural network can be treated as a functioñ y = g(x; W, b), which transfers non-linearly a given input x to an output predictionsỹ. The cost function of the network C(ỹ, y) is used to measure the difference between the network predictionsỹ and their true values y. The neural network is trained to minimize C(ỹ, y), by adjusting its parameters W and b. For the autoencoder network used in the present work shown in Fig. 4, the optimization is fulfilled by the Adam [51] package in Tensorflow for this full-connected network. When training the network, we use an exponential decreasing learning rate α = 10 −3 + (10 −3 − 10 −6 ) exp(−i/10000), with i the training epoch. The cost function is defined as C(ỹ, y) = (ỹ − y) 2 . To prevent the network from over fitting the data, we adopt a standard l 2 regularization term in the cost function of this neural network. The cost function is adjusted to include the norm of the weight W and the bias b, i.e.,C(ỹ, y) = C(ỹ, y) + l 2 ( W 2 /2 + b 2 /2), with l 2 a positive number. The l 2 regularization prevents the weights and biases from increasing to arbitrary large values during the optimization. We list the information of the autoencoder network in Table I. The encoder part and decoder part are mirror symmetric and each consists of two layers. Since the output layer represents the positive defined reconstructed charge-weighted charge multiplicity distribution ZM c (Z), a ReLU is used to connect the decoder and the output layer. Figs. 3 and 5 in the main text show the result with 8 and 4 neurons in the first and second encoder layer, respectively, and ReLU activation between encoder part and latent variable (with regularization l 2 = 0.1). Note that not all the network constructions shown in the second column of Table I restore properly the inputted charge multiplicity distributions. For those who restore properly the charge multiplicity distributions, their feature of the latent variable shown in Fig. 5 are quite similar. As a supplement to Fig. 3, we show the mean and standard deviation of the event-by-event reconstruction loss in Fig. 8. The eventby-event reconstruction loss is defined as the sum of absolute deviation of the charge multiplicity distribution, 12 Z=1 |ZM c (Z) − ZM c (Z)|. We note that the average reconstruction loss for the event with low temperature is larger than that with higher temperature. Considering ZMc(Z) is normalized to 12 and the event-by-event nature, we think this is still acceptable. For the supervised learning BNN used in the present work shown in Fig. 6, we use two hidden layers, each consists of 100 neurons. The input layer and the hidden layers are followed by ReLU. Weights and biases of the neurons are represented as Gaussian distributions. The maximization of the evidence lower bound (ELBO) is employed to solve the Bayes' formula. Here we do not employ the l 2 regularization since the problem of over-fit is not severe in BNN.