Nuclear ground states in a consistent implementation of the time-dependent density matrix approach
NNuclear ground states in a consistent implementation of the time-dependent densitymatrix approach
Matthew Barton,
1, 2
Paul Stevenson, ∗ and Arnau Rios
1, 3, † Department of Physics, Faculty of Engineering and Physical Sciences,University of Surrey, Guildford, Surrey GU2 7XH, United Kingdom Nuclear Theory Group, Faculty of Physics, Warsaw University of Technology, 00-662 Warsaw, Poland Departament de F´ısica Qu`antica i Astrof´ısica, Institut de Ci`encies del Cosmos (ICCUB),Universitat de Barcelona, Mart´ı i Franqu`es 1, E08028 Barcelona, Spain (Dated: January 20, 2021)
Background:
Time-dependent techniques in nuclear theory often rely on mean-field or Hartree-Fock descriptions.Beyond mean-field dynamical calculations within the time-dependent density matrix (TDDM) theory have ofteninvoked symmetry restrictions and ignored the connection between the mean-field and the induced interaction.
Purpose:
We study the ground states obtained in a TDDM approach for nuclei from A = 12 to A = 24, includingexamples of even and odd-even nuclei with and without intrinsic deformation. We overcome previous limitationsusing three-dimensional simulations and employ density-independent Skyrme interactions self-consistently. Methods:
The correlated ground states are found starting from the Hartree-Fock solution, by adiabaticallyincluding the beyond-mean-field terms in real time.
Results:
We find that, within this approach, correlations are responsible for ≈ −
5% of the total energy. Radiiare generally unaffected by the introduction of beyond mean-field correlations. Large nuclear correlation entropiesare associated to large correlation energies.By all measures, C is the most correlated isotope in the mass region considered.
Conclusions:
Our work is the starting point of a consistent implementation of the TDDM technique for appli-cations into nuclear reactions. Our results indicate that correlation effects in structure are small, butbeyond-mean-field dynamical simulations could provide new insight into several issues of interest.
I. INTRODUCTION
The study of the time evolution of nuclei provides keyinsights into their structure, their excitation and their as-sociate reactions. Several techniques have been devisedto tackle numerically the dynamics of nuclear many-bodysystems [1–4]. Traditionally, non-stationay simulationsin nuclear physics have been implemented in the time-dependent Hartree-Fock approximation (TDHF) [5–7] or,in more modern terms, the time-dependent density func-tional approach [8–10]. This approach assumes thatnucleons move only under an (instantaneous) averagepotential generated by the other nucleons and consis-tently takes the Pauli exclusion principle into account[11]. Using Skyrme density functionals, simulations arenowadays routinely implemented in unrestricted three-dimensional (3D) geometries and have been used to de-scribe a plethora of different nuclear phenomena [12–21].In the past, there have been attempts to move beyondthis mean-field approximation. There are a handful ofmethods that introduce genuine two-body correlationsin the dynamics [22]. These include, amongst others,the Balian-V´eneroni approach to incorporate particle-number fluctuations [23–25] and the Kadanoff-Baym ap-proach for infinite [26, 27] and finite [28, 29] systems.However, the Time-Dependent Density Matrix (TDDM) ∗ [email protected] † [email protected] approach is probably the most widely applied beyond-mean-field method in nuclear physics and, in this sense,it is the most successful to date.TDDM was first introduced by Cassing and Wang in1980 [30], and has been applied extensively in nuclearphysics by Tohyama [31–43] and others [44–46]. In thiscontext, the TDDM equations are often projected into amoving basis dictated by a TDHF-like equation, plus atime-dependent term that depends on correlations [44].Successful implementations with different levels of con-sistency have also been used to study breakup [47, 48]and, recently, fusion reactions in an energy-conservingapproach [49, 50].TDDM allows one to go beyond TDHF by trun-cating the Bogoliubov–Born–Green–Kirkwood–Yvon(BBGKY) hierarchy of quantum mechanical many-bodydensity matrices order by order [22, 51]. Here, and inthe following, we define uncorrelated systems as thosewhere the probability distributions of two particlesare independent. Two-body correlations are thereforea measure of the “lack” of independence of the twoprobability distributions. This approach has also seenwidespread use within other areas of physics, such ascondensed matter and quantum optics, where TDDMoften goes by the name of reduced density matrix theory[52–54].In addition to the many-body truncation, otherapproximations are typically implemented in nuclearTDDM simulations. Due to computational intensivenessof the calculations, previous implementations of TDDM a r X i v : . [ nu c l - t h ] J a n have worked in spherical or, more recently, axial sym-metry [37, 41]. Moreover, one often assumes that a de-coupling applies to nuclear systems, so that the inter-action acting at the mean-field level is different to thatacting at the beyond mean-field level [26]. The latter isoften dubbed the “residual interaction” and often takesthe form of a simplified δ − function [38, 47, 49]. On theone hand, the geometric restrictions preclude the appli-cation of TDDM methods to triaxially deformed nucleiand to general dynamical settings between multiple nu-clei. On the other, the inconsistent use of mean-field andresidual interactions hampers the possibility of discussingsystematics in the TDDM expansion. For instance, if onewere to find an improvement in the dynamical descriptionwhen going from TDHF to TDDM, it would be difficultto unambiguously ascribe the improvement to the use ofTDDM when the employed Hamiltonian is not the sameat all levels.In view of these limitations, we have implemented afully unrestricted 3D implementation of TDDM that usesthe Skyrme interaction for both the mean-field and theresidual channels [55]. Our truncation of the BBGKY hi-erarchy includes two-body correlations only [42]. In thispaper, we provide details and results for our implementa-tion of this method. Additional information can be foundin the PhD thesis of Ref [55].Our physics focus is the generation of correlatedground states within the TDDM approach. We obtainthese from the dynamical equations by means of an adi-abatic switching-on technique, as explained below. Ourcalculations extend from light systems, like C, up toexotic nuclei, like O. We do not expect the calcula-tions to provide a good match to experimental data, be-cause the interactions have been fitted at the mean-fieldlevel. However, the simulations provided here are infor-mative in terms of the structure, size and mass evolutionof two-body correlations within the TDDM formalismwith Skyrme forces. Ultimately, our aim is to use theground states described here to study nuclear dynamicswithin a fully consistent TDDM approach.This paper is laid out as follows. Section II gives a briefoutline of the theoretical background and the numericalimplementation of our TDDM approach. Some furtherdetails are provided in the Appendix.In Sec. III, we discuss nuclear ground states obtainedwithin HF calculations, which are necessary for compari-son to the correlated TDDM results provided in Sec. IV.Section V concludes this paper with a summary andshort discussion on areas of future research.
II. TIME-DEPENDENT DENSITY MATRIXMETHODA. Formalism
The BBGKY hierarchy relates the time evolution ofan A − body density matrix, ρ A , to the Hamiltonian, ˆ H , and the ( A + 1) density matrix [30]. A truncation of thehierarchy is necessary to make the dynamical equationsof the density matrix numerically tractable for practicalimplementations. Depending on the truncation, one findsdifferent coupled differential equations for the evolutionof the density matrix that obey conservation laws [52, 53].By assuming that A = 3 body correlations can be castin terms of A = 2 and A = 1 density matrices only,one recovers the most popular implementation of TDDM[42, 45, 56].If we denote generally the space coordinates of a nu-cleon by x i , the two-body density matrix (with no ref-erence to spin or isospin) is a tensor in four positions, ρ ( x (cid:48) , x (cid:48) ; x , x ). In a 3D mesh of N x points in each di-rection, this quantity scales like N x , quickly overcomingpresent computational capabilities. To avoid this limita-tion, we solve the TDDM equations in a moving TDHF-like basis [44]. This has several advantages. First, be-cause part of the dynamics is incorporated in the ba-sis, the TDDM equations are simplified with respect tostatic basis approaches [45]. More importantly, the sizeof the correlation tensor is dictated by the total num-ber of single-particle orbitals, N max , and scales with thefourth power of this variable, N . In addition, we canuse already existing computational capabilities at theTDHF level to evolve the basis states in a fully unre-stricted 3D geometry [17]. We note, however, that thereare instances, particularly in fusion reactions in the merg-ing phase, where a finite basis set may be insufficient toguarantee energy conservation [49, 50].In this approach, the one-body density matrix is ex-panded into a finite set of HF-like single-particle orbitalsthat depend on time, ρ ( x , x (cid:48) ; t ) = N max (cid:88) αα (cid:48) n αα (cid:48) ( t ) ψ ∗ α (cid:48) ( x (cid:48) , t ) ψ α ( x , t ) . (1)In this subsection, we denote by N max the total numberof such states, including neutrons and protons . Fromnow on we omit the time variable t as it is clear thatall quantities depend on it and all our summations runfrom the lowest single-particle index up to N max . Thesingle-particle orbitals follow the dynamics dictated by aTDHF-like equation, i (cid:126) ddt ψ λ = (cid:88) α (cid:15) αλ ψ α , (2)where (cid:15) αβ = t αβ + (cid:88) γδ ν αγβδ n δγ . (3) In the numerical implementations discussed below, we also de-note the total number of neutron and of proton states, indepen-dently, by N max . The factor of 2 between the two definitionshould not cause any confusion. is the so-called energy matrix. This includes a kineticcontribution, t αβ , and an interaction term. We give moredetails on the calculation of the matrix elements ν αγβδγ below. If the energy matrix is diagonal, the elements (cid:15) αα ,are the single-particle energies associated with a givenstate α .The matrix n αα (cid:48) is known as the occupation matrix.When the occupation matrix is diagonal, the diagonal el- ements correspond to the occupations of the associatedsingle-particle orbitals. The time evolution of n αα (cid:48) isdictated by the correlation tensor, C . The latter cor-responds to the correlated part of the two-body densitymatrix, C = A ( ρ ρ ) − ρ , where the operator A an-tisymmetrizes with respect to exchanges between singleindices x i and x j . The correlation tensor can also bedecomposed into a time-dependent single-particle basis, C ( x , x ; x (cid:48) , x (cid:48) ) = N max (cid:88) αβα (cid:48) β (cid:48) C αβα (cid:48) β (cid:48) ψ ∗ α (cid:48) ( x ) ψ ∗ β (cid:48) ( x ) ψ α ( x (cid:48) ) ψ β ( x (cid:48) ) . (4)Upon making this decomposition, one finds that the evo-lution of n αα (cid:48) becomes: i (cid:126) dn αα (cid:48) dt = N max (cid:88) γδσ (cid:104) ν ασγδ C γδα (cid:48) σ − C αδγσ ν γσα (cid:48) δ (cid:105) . (5)As opposed to a static basis projection, the right-handside of this equations does not have any Hartree-Fock(HF) term [44, 45]. As one can clearly see, when corre-lations are not active ( C = 0), the occupation probabili-ties become static. Further, in a pure mean-field picturewithout pairing correlations, one can prove that theseoccupation probabilities are either 1 or 0 depending onwhether the orbital is below or above the Fermi surface,respectively.In contrast, when correlations are active, one expectsthat the occupation numbers take values between 0 and1, to abide with the Pauli principle and their probabilisticnature. The truncation in the TDDM equations does not always mathematically guarantee that this is the case, asnumerically corroborated in early nuclear physics appli-cations [57, 58] and more recent strongly correlated elec-tronic simulations [52]. We have observed this anomalousbehaviour in a handful of simulations, but it is difficultto discriminate their origin that could partially be due tonumerical issues.When the hierarchy is truncated at some level, theevolution of the correlation tensor C is dictated by anequation which depends on occupation numbers; inter-action matrix elements; and the correlation tensor itself.We work under the assumption that genuine three-bodycorrelations are negligible [43, 44, 52]. In other words, thethree-body density matrix is a properly antisymmetrizedproduct of one-body and two-body density matrices only,but does not include any genuine C terms. Under thisapproximation, the equation of motion for the correlationtensor is [55] : i (cid:126) dC αβα (cid:48) β (cid:48) dt = 12 (cid:88) λµ ν βαλµ ( n λβ (cid:48) n µα (cid:48) − n λα (cid:48) n µβ (cid:48) + C µλα (cid:48) β (cid:48) ) + 12 (cid:88) λµ ν λµα (cid:48) β (cid:48) ( n βλ n αµ − n βµ n αλ + C αβµλ ) − (cid:88) δλµ ν βδλµ (cid:34) n λβ (cid:48) ( n µα (cid:48) n αδ − n αα (cid:48) n µδ + C µαα (cid:48) δ ) − n λα (cid:48) C µαβ (cid:48) δ − n µβ (cid:48) C λαα (cid:48) δ + n µα (cid:48) C λαβ (cid:48) δ + n αδ C λµβ (cid:48) α (cid:48) (cid:35) + 12 (cid:88) δλµ ν δλβ (cid:48) µ (cid:34) n βλ ( n µδ n αα (cid:48) − n αδ n µα (cid:48) − C αµδα (cid:48) ) + n βδ C αµλα (cid:48) − n αδ C βµλα (cid:48) + n αλ C βµδα (cid:48) + n µα (cid:48) C βαδλ (cid:35) + 12 (cid:88) δλµ ν αδλµ (cid:34) n λβ (cid:48) ( n µα (cid:48) n βδ − n βα (cid:48) n µδ + C µβα (cid:48) δ ) − n λα (cid:48) C µββ (cid:48) δ − n µβ (cid:48) C λβα (cid:48) δ + n µα (cid:48) C λββ (cid:48) n + n βn C λµβ (cid:48) α (cid:48) (cid:35) − (cid:88) δλλ ν δλα (cid:48) λ (cid:34) n βλ ( n λδ n αβ (cid:48) − n αδ n λβ (cid:48) − C αλδβ (cid:48) ) + n βδ C αλλβ (cid:48) − n αδ C βλλβ (cid:48) + n αλ C βλδβ (cid:48) + n λβ (cid:48) C βαδλ (cid:35) . (6)We note that this equation uses antisymmetrized matrix elements [see Eq. (10) below] unlike other implementa-tions [44, 45].A brute force implementation of the previous equationswould scale as N . We exploited the symmetries of thecorrelation tensor elements to reduce this computationalcost [55].Also, certain matrix elements are zero based on isospinconservation arguments.The two-body density matrix (or, equivalently, the cor-relation tensor C ) provides direct access to the total en-ergy of the system, which is customarily split into twocontributions, E = E MF + E c . The mean-field term isexpressed in terms of occupation matrix elements onlyand is already active at the HF level, E MF = (cid:88) αβ t αβ n βα + 12 (cid:88) αβγδ ν αβγδ n γα n δβ . (7)The correlation energy term, in contrast, is directly pro-portional to C , E c = 14 (cid:88) αβγδ ν αβγδ C δγβα , (8)and is only active in beyond mean-field calculations. Thecorrelation energy can therefore be used as a metricto quantify correlations in the TDDM approach. TheTDDM approach based on Eqs. (2), (5) and (6) con-serves the number of particles, the total momentum andthe total energy [44]. B. Interaction matrix elements
In the past, the implementation of the TDDM ap-proach has often relied on approximations. An often-usedassumption is that the beyond mean field interaction (the“residual interaction”) is a delta function multiplied by aconstant. This reduces substantially the computationalcost required to calculate the matrix elements, ν αβα (cid:48) β (cid:48) .This approximation however ignores the self-consistencybetween the mean-field interaction and residual interac-tion which, from a first-principles perspective, should bebased on the same Hamiltonian. In this work, we insteaduse the Skyrme interaction, V ( r ) = t (1 + x P s ) δ ( r ) + 16 t (1 + x P s ) ρ α ( R ) δ ( r )+ 12 t (1 + x P s ) (cid:104) k (cid:48) δ ( r ) + δ ( r ) k (cid:105) + t (1 + x P s ) [ k · δ ( r ) k ]+ iW ( σ + σ ) · [ k (cid:48) × δ ( r ) k ] , (9)to model the nucleon-nucleon interaction, both at themean-field and the residual interaction level. In thisequation, r = r − r represents the relative distancebetween two nucleons at positions r and r ; R =( r + r ) / k = ( ∇ − ∇ ) / i the relative momentumacting on the right and k (cid:48) its conjugate acting on the left. P s = (1 + σ · σ ) / is the spin exchange operator.The last term, proportional to W , corresponds to thezero-range spin-orbit term.The use of interactions with density-dependent termsin beyond mean-field implementations may be problem-atic. These do not constitute true forces and hencemust be treated with care in many-body approaches toavoid pathologies [55]. We therefore preclude from us-ing standard Skyrme parametrizations, but employ twoparametrizations of this force, SV [59] and SHZ2 [60],that do not have a density dependent term. In otherwords, t = 0 in the notation of the original Skyrmeforce [61, 62]. We note that SHZ2 is in fact a slight refitof SV, with very similat t i parameters and a very small x term [60]. These two different fits therefore allow usto minimally explore the parameterization-dependence ofour results.The matrix elements of the interaction need to be pro-jected into the single-particle orbitals so they can be usedin Eqs. (5) and (6). This is achieved by means of a double3D integral ν αβα (cid:48) β (cid:48) = (cid:90) d x d x ψ ∗ α ( x ) ψ ∗ β ( x ) V ( x − x ) × [ ψ α (cid:48) ( x ) ψ β (cid:48) ( x ) − ψ β (cid:48) ( x ) ψ α (cid:48) ( x )] . (10)However, because of the zero-range nature of the Skyrmeforce, these integrals simplify substantially [55]. The cal-culation of all elements of ν αβα (cid:48) β (cid:48) scales as N x N max .This quantity is calculated 4 times at each time-step,which becomes a numerical bottleneck for very fine gridsor large boxes, and for heavy systems. We note that thesematrix elements are antisymmetrized from the outset. C. Adiabatic switching and asymptoticconvergence of correlation energies
We obtain nuclear ground states employing an adi-abatic real time switching technique that continuouslytransitions from the mean-field to the correlated groundstate [55]. This is achieved by multiplying the matrixelements of the interaction that appear in Eqs. (5) and(6) by afactor, γ ( t ), that slowly goes from 0 to 1.We use a Gaussian factor, γ ( t ) = 1 − e − t τ , (11)which has finite derivatives at t = 0 and t (cid:29)
1. Wework under the assumption that the switch-on procedureallows the Gell–Mann Low Theorem [63] to be applied.In other words, if the residual interaction is switched onslowly enough, the final state should become an eigen-state of the TDDM approach. Numerical tests indicatethat the value τ = 32000 fm c − is sufficient to guaran-tee a converged correlated ground state. This correspondto physical changes on a timescale of t = √ τ ≈
180 fmc − .We indicate that a test run with A = 4 and τ =64000 fm c − provided no significant differences in termsof asymptotic energies. The asymptotic energy valuesprovided below are obtained either from a converged finalresult, or from fits of the different energy componentsassuming a time dependence proportional to γ ( t ). Moredetails of this procedure can be found in Ref. [55]. D. Numerical details and bottlenecks
In unrestricted 3D TDHF simulations, one typicallyworks with as many single-particle orbitals as nucleonsin the system, N max = A [17]. As nucleons are al-lowed to scatter off each other in TDDM, the numberof single-particle orbitals must necessarily be larger thanthe number of nucleons, N max > A . Our TDDM simula-tions are projected into a TDHF-like basis with an equalmaximum number of neutron and proton states, N max .In the following, we provide results for N max = 14 and20 to explore what in ab initio terms is typically calledthe “model-space” dependence of our results. In a shellmodel language, N max = 20 corresponds to a model spacespanning the full sd shell.There are two major numerical bottlenecks in our ap-proach, that affect the ability to propagate over time andrestrict the size of nuclei that can be tackled. First,simulations are expensive in terms of memory require-ments, since the correlation tensor C and the interactionmatrix elements both scale like N in number of ele-ments. Large amounts of memory are required to storethese tensors. Second, the calculations of both C α (cid:48) β (cid:48) αβ and ν α (cid:48) β (cid:48) αβ are time-consuming. As reported before,these scale as N and N N x , respectively. Gener-ally speaking, for a small model space ( N max < N x ), thecalculation of the interaction matrix elements takes mostof the computational time. For larger model spaces, itis the calculation of C α (cid:48) β (cid:48) αβ that dominates the com-putational cost. We note that parallelization helps incomputing some of these tensors at each time step.All calculations presented here were performed on aCartesian 3D grid with spacings ∆ x = ∆ y = ∆ z = 1 fmfrom − . . N x = 20. For the relativelylight nuclei in consideration here, we operate in a regimewhere N x ≈ N max . We note that for He, a smaller gridspacing of 0 . C and ν are computed at eachtime-step, which makes dynamical simulations slow. Thethree differential equations for the evolution of the sin-gle particle orbitals [Eq. (2)], occupations [Eq. (5)] andcorrelation tensor elements [Eq. (6)] are solved using a4 point Runge-Kutta method. We note that the tradi-tional midpoint method commonly used in TDHF [17, 64] SV Charge radius[fm] B.E. [MeV] Point radius[fm]Nucleus HF Exp HF Exp Proton Neutron C 2.724 2.4702 69.432 92.160 2.604 2.587 O 2.765 2.6991 113.536 127.616 2.647 2.629 Ne 3.058 3.0055 138.860 160.64 2.951 2.927 Ne 3.046 2.9695 146.79 167.391 2.939 2.991 Na 3.129 3.0136 142.926 163.065 3.025 2.919 Na 3.109 2.9852 153.626 174.130 3.004 2.978 O 2.800 N/A 144.552 168.96 2.683 3.433SHZ2 Charge radius[fm] B.E. [MeV] Point radius[fm]Nucleus HF Exp HF Exp Proton Neutron O 2.762 2.6991 113.648 127.616 2.644 2.624 Ne 3.054 3.0055 139.140 160.64 2.947 2.919TABLE I. Charge radii (columns 2-3) and binding energies(columns 4-5) obtained from HF calculations using the SV(top values) and SHZ2 (bottom values) Skyrme forces, along-side experimental values. Proton and neutron point radii arereported in columns 6-7. Different nuclei are listed in eachrow. Experimental data taken from Refs. [65, 66]. provided unstable results in the TDDM implementation(unless a very small time step was used). In all calcula-tions performed in this work, a value of dt = 0 . − was used.Further details about the time stepping procedure areprovided in the Appendix. III. MEAN-FIELD GROUND STATES
We start the discussion of results by providing some ofthe bulk properties of the HF ground states with the SV[59] and SHZ2 [60] Skyrme interactions. These results actas a baseline and allow us to quantify the importance ofthe correlations induced by the TDDM approach. Theground states are obtained using Sky3D [17] and theirproperties are summarised in Table I. We investigate awide range of nuclei from A = 12 to A = 24. The toppart of Table I shows results for the SV force, whereasthe bottom part shows SHZ2 results for O, Ne and O. We provide experimental data where known.Our final aim is not so much to produce a set of reliableground states to compare to experiments, but rather todevelop an understanding of the size and structure of thecorrelations induced by the TDDM approach.The charge radii provided in column 2 are generallyoverestimated with respect to the experimental results(column 3) by about 3% on average. The HF bindingenergies per nucleon are provided in column 4 of TableI. The theoretical results underestimate the experimentalones (column 5) by about 12% on average. We can put C SV N max
14 20 E c [MeV] − . − . E MF [MeV] − . − . E [MeV] − . − . E c /E [%] 7 .
72 11 . .
644 2 . ± . .
627 2 . ± . Cground states for different values of N max (columns 2 and 3).Results are provided for the SV interaction forward some explanations for this discrepancy. First,we note that our results do not include a center of masscorrection, which will be relevant for the energetics of thelightest isotopes. In fact, the binding energies are some-what closer to experimental results as A increases, sug-gesting this is the case. Second, these effective interac-tions were fitted to the ground-state properties of spher-ical system from A = 16 to 208, with very many moreheavy systems than light isotopes in the fitting protocol.This naturally biases the parametrizations towards heav-ier nuclei. Finally, the HF approximation is expected towork better for heavier than for light systems on generalgrounds.Overall, however, the HF simulations produce rea-sonable values of the energy. The mass dependence ofthe simulations follows reasonably the energy and radiusdata. We also stress that there are relatively small differ-ences between the results obtained with the two Skyrmeinteractions. The charge radii (energies) obtained withSHZ2 are negligibly smaller (larger) than those of SV,in agreement with the fact that this force has a slightlylarger saturation density. The relative differences are ofthe order ≈ . − . IV. CORRELATED GROUND STATES
We now discuss the results obtained for the correlatedTDDM eigenstates. We aim at providing as much of asystematic discussion as possible by focusing on bind-ing energies and radii. We discuss the results isotope byisotope, in order to provide a more detailed explanationand a clearer characterisation of the role of correlationsin each of these systems. -73-72-71-70-69-68-67-66-65-640 50 100 150 200 250 300 350 400 450 50012 C E n e r gy [ M e V ] Time [fm/c]
Total N max =14Total N max =20Mean-feld N max =14Mean-feld N max =20 FIG. 1. Time evolution of the mean-field (open symbols)and the total energy (filled symbols) of C from a HF to acorrelated ground state. Data are provided for N max = 14(squares) and N max = 20 (circles). A. C The ground state structure of C is of considerableinterest for a variety of reasons. In particular, C is rele-vant because of its possible cluster structure, in which in-dividual nucleons may be correlated with others in a waythat cannot be easily captured in a mean-field descrip-tion [67]. It is conceivable that the correlations inducedby TDDM can capture some of the clusterization mech-anisms and provide significantly different ground states.We summarise the TDDM results for the structure of C in Table II. These results are obtained with the SVparametrization. The uncorrelated, HF ground state hasa total energy of E = − . N max = 14 (20). As correlations areintroduced in the system, the energy changes. The totalenergy becomes more attractive, whereas the mean-fieldcontribution yields more repulsive results. The total en-ergy drops to ≈ −
71 MeV. In contrast, the mean-fieldcomponent increases by about 4 to 5 MeV. Importantly,the final energy is not completely stationary after theevolution finishes at t = 500 fm/c.In relative terms, the correlation energy shown in II asa percentage of the total energy is between 7 and 11%.We anticipate that this is over twice that of any of theother nuclei discussed in the following, which we take asan indication of the importance of correlations for thisspecific isotope. It would be interesting to find quanti-tative measures of clustering in these simulations, in linewith what has been achieved at the mean-field level [68].Columns 2 and 3 of Table II show results for both N max = 14 and 20, respectively. In a traditional shellmodel picture, the former would include the ground state1 s / and 1 p / configurations of C as well as the 1 p / and 1 d / levels. The larger model space N max = 20completely fills in the sd shell. We note that, in goingfrom the N max = 14 to N max = 20 configuration, thesystem gains about 2 . N max corresponds to a largercorrelation energy, since the interaction does not have anatural cut-off, and those levels nearest the Fermi energycan be scattered into most available levels.The oscillations in energies found at long times inFig. 1 are evidence of the fact that the system is evolv-ing into the correlated eigenstate too quickly. Turningon the residual interaction more slowly by increasing τ in Eq. (11) may remedy the oscillations at the end ofthe calculation, at increased computational cost. Note,however, that the oscillations in the mean field energyare compensated by anti-phase oscillations in the cor-relations energy (not shown), giving an overall smoothtotal energy as a function of time.These oscillation are also reflected in the rms radii,which oscillates with a typical size of order 0 .
01 fm forthe N max = 20 simulation. This uncertainty for radiiis reported in the bottom rows of Table II. Comparingthe rms radii to the HF values reported in I and the two N max values with one another,we find that the collisions allow nucleons to scatter fur-ther from the nucleus. We note that the proton rms ra-dius increases by about 0 . n αα , for C is shown in Fig. 2. The results areshown for the N max = 20 simulation. Left (right) panelscorrespond to neutron (proton) states. Top panels dis-play the 6 hole states for both species. Within TDHF,the protons and neutrons completely fill the 1s and 1p sub-shells. Bottom panels display particle states instead.We find that the mean-field picture is still mostly rele-vant for the correlated eigenstate in C. Hole state oc-cupations reach a value of about 95%. Here, there areclear differences between the more bound 1 s / states,which remain populated to a 99 .
5% level, and the 1 p / substates, that are substantially more depleted. The oc-cupations of particle states are of order 10 − or lower.Particle states closer to the Fermi surface, with smallervalues of α , are more occupied than states further away.We note that some states, like the hole α = 0 (1 s / )and particle α = 6 (1 p / ) are clearly well converged, inthe sense that they reach a constant occupation num-ber at large times in the adiabatic switching. Others,in contrast, are still evolving at the end of the simula- -6 -5 -4 -3 -2 O cc u p a t i o n n u m b e r s , n αα Time [fm/c] 0 100 200 300 400 500(d) Time [fm/c] α = α = α = α = α = α = α = O cc u p a t i o n n u m b e r s , n αα α = α = α = (b) α = α = α = FIG. 2. Occupation numbers of (a) neutron and (b) protonhole states as a function of time for the adiabatic switchingof C with N max = 20. Panels (c) and (d) give the corre-sponding occupations of particle states in a logarithmic scale.States that are degenerate in spin are not shown for simplicity.The index α denotes energy levels in increasing order. tion. This is the case of the hole α = 2 state (one of thetwo 1 p / states shown in short-dashed lines), but alsoof the α = 18 neutron state (double-dashed-dotted linein panel c) which has a very low occupation that turnsnegative just before the end of the evolution. We notethat states with large values of α are unbound (eg suchthat (cid:15) αα >
0, see next paragraph) and hence may besubstantially affected by box discretization issues.As discussed above, C has the largest relative correla-tion energy of the nuclei we discuss in the following. Thismay be surprising in the context of Fig. 2, which suggestsa relatively small redistribution of single-particle strengththat could be interpreted as having little impact in thenuclear structure (although, as we shall see later, thechanges are not insignificant). In addition, the diagonalelements of the single-particle energies themselves do notchange much. For C, the time evolution of the (cid:15) αα ele-ments are shown in panel (a) of Fig. 3 for the N max = 14simulations. The changes in these single-particle ener-gies are imperceptible in the scale of the graph, in linewith previous TDDM implementations [48]. To be quan-titative, the maximum change across the 500 fm/c of thesimulation, for the most bound α = 0 (1 s / ) state, isless than 0 . B. O O is a benchmark nucleus, as a light doubly-magicsystem which is open to calculation by many beyondmean-field methods. As with C, it is also an nα sys- -60-50-40-30-20-100100 100 200 300 400 500(b) O S i n g l e - p a r t i c l e e n e r g i e s , ε αα Time [fm/c] α = α = α = α = α = α = α = -50-40-30-20-10010 (a) C S i n g l e - p a r t i c l e e n e r g i e s , ε αα α = α = α = α = α = α = α = FIG. 3. Evolution of the single-particle energies (cid:15) αα with timefor (a) C and (b) O. Particle and hole states are shownin different colors and linestyles. These results are obtainedwith the SV interaction with N max = 14. tem, where clustering correlations may play a substantialrole. A graph of the mean field and total energy of Owith the SV Skyrme force parameterisation as it evolvesfrom the HF eigenstate to the correlated eigenstate, forvarious N max , is shown in Fig. 4. Simulations start in theHF ground-state at around ≈ − . − . − . N max = 14 (20)simulation, whereas the mean-field energy becomes about ≈ . . − .
5% to thetotal energy - far less than in the case of C. We alsonote that no oscillations appear in the total energy or itscomponents in the large-time limit. This may indicatethat the transition to the TDDM is somehow “easier” inthis less correlated nucleus.We simulate the correlated eigenstate of O using thetwo chosen Skyrme forces, SV and SHZ2. The summaryof results shown in Table III indicates an insignificantdifference between the two interactions, for both valuesof N max . For N max = 14, the correlation energy is − . N max = 20,the correlation energy increases to − . − . -116-115-114-113-112-1110 50 100 150 200 250 300 350 400 450 50016 O E n e r gy [ M e V ] Time [fm/c]
Total N max =14Total N max =20Mean-feld N max =14Mean-feld N max =20 FIG. 4. Same as Fig. 1 for O. ascribe these minutes differences to the fact that the twoSkyrme forces, themselves, are very similar to each other.As we have already seen, the occupation probabilities n αα provide a way of characterising correlations. Theirtime evolution within the adiabatic switch-on processfor O is shown in Fig. 5. These have been obtainedwith the SV interaction in the N max = 14 model space.For O, HF simulations include 16 single-particle or-bitals corresponding to the 1 s / , 1 p / and 1 p / neu-tron and proton hole subshells. Like the C case, thedeeply bound 1 s state (solid line in panels (a) and (b)of Fig. 5) remains practically fully occupied. In contrastto the previous case, the occupation of the hole 1 p states[dashed and dotted lines in panels (a) and (b)] is at levelof ≈ C , where the same orbital was depleted by almost5%.The particle states of panels (c) and (d) tell a similarstory. Whereas levels close to the Fermi surface for Creached relatively large occupations of order 10 − , all O SV O SHZ2 N max
14 20 14 20 E c [MeV] − . − . − . − . E MF [MeV] − . − . − . − . E [MeV] − . − . − . − . E c /E [%] 3 .
55 4 . .
55 4 . .
660 N/A 2 .
657 N/Arms [fm]Neutron 2 .
640 N/A 2 .
636 N/Arms [fm]TABLE III. The same as Table II for O. Results for theSHZ2 parametrization are also provided in columns 4 and 5. -6 -5 -4 -3 O cc u p a t i o n n u m b e r s , n αα Time [fm/c] 0 100 200 300 400 500(d) Time [fm/c] α = α = α = O cc u p a t i o n n u m b e r s , n αα α = α = α = α = (b) α = α = α = α = FIG. 5. Same as Fig. 2 for O with N max = 14. (degenerate 1 d / ) particle states of O show a muchsmaller final occupation, close to 0 . C, the diagonal elements of thesingle-particle energy matrix shown in panel (b) of Fig. 3are relatively constant across the adiabatic evolution forthe N max = 14 model space. Unlike the high-lying Cparticle states, all particle states in O are bound so wedo not anticipate any continuum discretization issues inour simulations. O has been used as a benchmark nuclear sys-tem in the past, including different implementations ofTDDM [35, 39, 40, 47, 48]. These studies have generallyrelied on different mean-field and residual interactions;have neglected spin-orbit couplings in the residual chan-nel and/or have restricted the relevant correlation modelspace to p and d subshell orbitals. The results typicallyobtained in these models are much more correlated thanthose discussed here. Typical p − shell ( d -shell) orbitaloccupations are closer to 90% (10%), and correlation en-ergies are (cid:39) −
10 MeV. Without a more in-depth analysisand lacking unbiased benchmarks, it is difficult to find aclear explanation for the origin of these discrepancies. C. Ne We discuss the isotope Ne as the first of a seriesof examples centered around A = 20. This region ofthe chart has received significant experimental attention[69, 70] due, among other things, to its relevance for as-trophysics [71]. In theoretical studies, this region is typ-ically accessed theoretically using the shell model and isof particular interest in the context of isospin symmetrybreaking [72]. In a standard shell model picture, Ne isbuilt from an O core by adding two neutrons and two Ne SV Ne SHZ2 N max
14 20 14 20 E c [MeV] − . − . − . − . E MF [MeV] − . − . − . − . E [MeV] − . − . − . − . E c /E [%] 1 . .
15 1 .
43 4 . .
955 2 .
980 2 .
950 2 . .
930 2 .
946 2 .
922 2 . Ne. Results for SHZ2are also provided. -142-141-140-139-138-137-136-135-1340 50 100 150 200 250 300 350 400 450 50020 Ne E n e r gy [ M e V ] Time [fm/c]
Total N max =14Total N max =20Mean-feld N max =14Mean-feld N max =20 FIG. 6. Same as Fig. 1 for Ne. protons. It is, of course, yet another nα system.We provide a figure for the time evolution of the mean-field (dashed lines) and total (solid lines) energies of Nein Fig. 6. As with previous cases, the results are shownfor two values of N max for the SV force. The N max = 14results (squares) converge well with time, and show nosigns of oscillations at late times. The N max = 20 simula-tion stopped some time before 500 fm/c, but the resultsappear to be relatively well converged at this level. Akey difference between the results shown in this figureand those of previous isotopes is the relatively large dif-ference in energies between the results obtained with thetwo model spaces. When going from N max = 14 to 20,the total energy decreases by almost 2 MeV - more thandouble the result observed in other isotopes. Anotherstriking feature is the large increase in the ratio E C /E ,which more than doubles when going from one modelspace to the other. We interpret these results as a signthat the 1 d / subshell closure obtained in the N max = 14model space is relatively weak in this nucleus. The full0 Ne SV N max
14 20 E c [MeV] − . − . E MF [MeV] − . − . E [MeV] − . − . E c /E [%] 0 .
88 4 . Ne ground states for differentvalues of N max obtained with the SV force. sd shell of the N max = 20 simulation provides a muchmore complete picture that substantially enhances thecorrelation of the system.Table IV provides a breakdown of the energy contri-butions for the two Skyrme forces, SV and SHZ2. For N max = 14, the results obtained with these two param-eterisations are almost indistinguishable in terms of cor-relation energy. As one increases the model space to N max = 20, both parametrisations predict the aforemen-tioned substantial increase in correlation energies, from − . − . − .
9) MeV for the SV (SHZ2) force. D. Ne We turn our attention to an open-shell, odd-even anddeformed system, Ne, to confirm that such systems canbe tackled with our approach.For Ne, Table V summarises the energy contributionsfor the two values of N max . Just like in the previous case,we find a substantial increase of the correlation energy(more than a factor of 4) when moving from N max = 14to 20. In turn, the relative contribution to the energy in-creases from below 1% to over 4%. This clearly indicatesthe importance of sd shell contributions for this regionof the nuclear chart.We also find interesting systematics when comparingthe Ne results of Table V to the Ne simulations pre-sented in Table IV. With the addition of one neutronon top of Ne, for instance, we observe that the cor-relation energy drops by 0 . .
2) MeV, in the case of N max = 14 (20). This drop in correlation energy can beunderstood naively, in terms of a reduction in the num-ber of neutron levels available to scatter into. As for thetotal energy, the Hartree-Fock prediction for Ne is ≈ Ne. The TDDM ground statesenergies of the two isotopes differ by 7 . Ne is such that, forboth N max = 14 and 20, the proton and neutron radiusdid not converge to a static result. This indicates that thetransition from the mean-field to the correlated state is Na SV N max
14 20 E c [MeV] − . − . E MF [MeV] − . − . E [MeV] − . − . E c /E [%] 0 .
97 4 . Na. more difficult than in some of the previous examples, pos-sibly because of the odd-even nature of the isotope. Wealso performed an analysis of some of the mean-field en-ergy components, not provided here for brevity. The datafor the t component of the mean-field (which is propor-tional to the overall density of the system) for N max = 14shows an increase of ≈ Ne, in contrast, went up by over 3 MeV. We take thisas an indication that the single addition of a neutroncan change significantly how different components of theSkyrme force change beyond the mean field limit. E. Na We continue our analysis by considering Na, the mir-ror nucleus to Ne with an odd proton number. Thisprovides an interesting insight into the nature of isospinsymmetry not only at the mean-field, but also at theTDDM level. The different energy contributions for Naare provided in Table VI for two values of N max . We findresults that bode well with those observed for the isospinpartner nucleus. First, as observed for the two previousisotopes, we find that the correlation energy increasessubstantially with the model space size: form − . N max = 14 to − . N max = 20. This correspondsto almost a factor of 4 in the relative contribution ofthe correlation energy, which increases from about 1% to4%. Second, comparing the correlation energy obtainedfor Ne with that of Na for N max = 14 (20) one seesthat the addition of one proton reduces the magnitude ofthe correlation energy by 0 . .
4) MeV. This mirrors thereduction we found for Ne compared to Ne. Again,this is presumably due to the reduction in the number oflevels available for the nucleons to scatter into.The results in tables V and VI allow us to analyse thelevel of isospin symmetry in our TDDM simulations. Atthe mean-field level, the results of Table I indicate a bind-ing energy difference between the two isotopes of ≈ . ≈ . Na SV N max E c [MeV] − . E MF [MeV] − . E [MeV] − . E c /E [%] 4 . Na, but for a singlevalue of N max . mean-field contribution including the density-dependentexchange term. As a result, we find that the total energydifference between the two isotopes remains very close to4 MeV, regardless of the model space, the same value asobtained in the mean-field simulation within our methoduncertainties.Looking at the specifics, we find that, with N max =20, Ne produces -6 . Nayields -5 . F. Na To finish the discussion in the A = 20 −
22 mass re-gion, we discuss the proof-of-principle case of Na, a nu-cleus with an odd number of neutrons and protons. Weencountered some technical issues in attempting to runsimulations of this isotope with N max = 14. We found,for instance, energy level crossings, which precluded usfrom identifying the final state of the evolution with theground state of the system. In addition, as the single-particle energies crossed, the occupations of both levelsapproached ≈ .
5, indicating a strong departure fromthe single-particle picture. Finally, the time evolutionof the energy departed significantly from the expected γ ( t ) dependence associated to the Born term. All inall, the N max = 14 results indicate that correlations aresignificantly changing the structure of this nucleus. It ispossible that an insufficiently large model space cannotcapture this significant changes in the switching proce-dure.In contrast, the numerics for the N max = 20 case wereremarkably more stable. Table VII provides a summaryof the energetics obtained for this isotope, focusing onlyon the N max = 20 results. We find a relatively largecorrelation energy of 7 MeV, which is about ≈ A = 21 isotopes and in goodagreement with the Ne result. In relative terms, this isabout ≈ .
5% of the total energy, close to the value thatwe have observed across this mass region. In other words,it appears that the instability in the N max = 14 resultsdoes not significantly reflect in the converged results. O SV N max E c [MeV] − . E MF [MeV] − . E [MeV] − . E c /E [%] 3 . O. G. O We finish our discussion with an exotic, neutron-richisotope: O. This provides a test case for a nucleusrelatively far from stability with a significant asymme-try between proton and neutrons. This isotope is in-deed close to the neutron drip line and is at the cen-ter of a series of contemporary experimental develop-ments [73–76]. Importantly, some results for this isotopehave been previously reported in other TDDM imple-mentations [38, 47, 48]. Our calculations were performedwith the SV parameterisation and a model space with N max = 20. The results are summarised in Table VIII.We predict a correlation energy in O which is ≈ . ≈
250 fm/c. Tohyama andUmar report a correlation energy of − . O inRef. [38], whereas Assi´e and Lacroix find − . P approach. Both values bode relativelywell with our finding, even though they have been ob-tained with different mean-field (and residual) interac-tions. More importantly, these predictions rely on usingonly a handful of active orbitals and, typically, an Oinert core.Compared to the equivalent results for the symmetricisotope O in Table III, the correlation energy has de-creased by about 0 . N = 8 to 16. This follows the qualita-tive trend discussed in previous isotopes, which indicatesa reduction of correlation energy as neutron number in-creases. These findings bode well with the idea that,in increasing neutron number, there are less orbitals toscatter into and, hence, less of a correlation energy. Thisdecrease is also consistent with the isotopic evolution re-ported in oxygen both in Refs. [38] and [48]. The resultsshown in Ref. [38] when going from O to O show adecrease of almost 1 MeV. The no-core simulations inRef. [48] also show a decrease of E c with neutron num-ber, although the order of magnitude of the correlationis different.Our simulations are also influenced by the closenessto the drip line. In the HF simulation, all 16 neutronstates are bound. Upon switching correlations on withTDDM, however, some of the unbound HF states be-come occupied through beyond mean-field scattering. Inparticular, there are 4 levels that are very close to beingbound with energies ≈ .
25 MeV. These almost-bound2 O C o rr e l a t i o n e n t r o p y [] Time [fm/c]N max =14N max =20
FIG. 7. Correlation entropy of O as a function of time fortwo different model spaces with N max = 14 (dashed line) and20 (solid line). levels also have a relative occupation which is two ordersof magnitude larger than any of the remaining unboundstates. In a naive shell model interpretation, one wouldinterpret these 4 additional as those filling the neutron sd shell. The treatment of a box-discretized continuum maybe somewhat inadequate here, but it does not precludethe convergence of our simulations. While we have notexplored these effects further, it is possible that addingmore neutrons to oxygen isotopes may shift occupationsand single-particle energies in a way that the mean-fieldsimulation cannot capture, providing perhaps a differentneutron drip line in TDDM than in HF. H. Correlation entropy
Up to this point, we have looked at the effect of cor-relations on specific, measurable single-particle and bulknuclear properties. There are however other theoreticalmeasures of correlations that provide independent char-acterisations. One of such measures is the so-called cor-relation entropy [77, 78], which is computed from thediagonal occupation numbers n αα as S cor = − A (cid:88) α n αα ln n αα , (12)where A is the number of nucleons. This quantity is,strictly speaking, not an entropy from a thermodynamicpoint of view, but rather an approximation of one as,among other things, one neglects the off-diagonal occu-pation matrix elements [78]. The correlation entropy isexactly 0 for the HF state, and necessarily increases asone goes towards a correlated eigenstate. The absolute numerical quantity does not have a direct physical mean-ing, but comparison of values between different calcula-tions can be instructive. Naively, one expects more “cor-related” ground states to yield larger values for S cor , inthe sense that they depart more from the HF eigenstatesof 0 entropy.A graph of this quantity as a single O nucleus goesfrom the HF eigenstate to the correlated eigenstate, forboth N max = 14 (dashed line) and 20 (solid), is shown inFig. 7. As expected, in the initial HF state both simula-tions yield 0 correlation entropy. As the system evolvestowards the correlated ground state, the entropy steadilyincreases until it levels off around 300 fm/c into the evo-lution. We note that the entropy does not completelyconverge at the end of the calculation, particularly for N max = 20, where the simulations suggest the occur-rence of a maximum of S cor at intermediate times. Wetake this as an indication that the occupation numbers n αα are not entirely converged, which in turn suggeststhat the adiabatic switching time is relatively small. Nei-ther the energy nor the occupation numbers reported inFigs. 4 and 5, however, showed a clear non-stationarityat the end of the simulation.Figure 7 also suggests that the correlation entropy in-creases with N max . We find that this is a generic fea-ture at least in the the two model spaces explored here.Table IX shows the correlation entropy obtained numer-ically at the end of the adiabatic evolution for the var-ious nuclei studied in this work. In all cases, the en-tropy computed with N max = 14 states is smaller thanthat computed with N max = 20. In a sense, this can beunderstood naively, in that additional levels necessarilyprovide more contributions to the correlation entropy. Inthis sense, the correlation energy is not a good measureof the model-space convergence of the results.Two more features stand out from the results on Ta-ble IX. First, we find that the nuclei with large correlationenergies, like C, also have large correlation entropies.Second, we find that the relative increase in correlationentropy when going from N max = 14 to 20 is very sim-ilar to the corresponding relative increase in correlationenergies. In other words, we find that both the correla-tion energy and the correlation entropy provide relativelysimilar measures of correlations in the systems that wehave studied. We note that this is not necessarily triv-ial a priori . The calculation of the correlation entropyrelies entirely on one-body occupation numbers, whereasthe correlation energy is the result of the contraction of2 seemingly different two-body objects - the interactionand the correlation tensors, see Eq. (8). V. CONCLUSIONS AND FUTURE OUTLOOK
In this paper, we implemented the TDDM approachincluding up to two-body correlations to study nuclearground states. Unlike some of the previous work in thefield, our simulations are performed without any spatial3 SV N max
14 20 C 0.088 0.132 O 0.041 0.053 Ne 0.020 0.073 Ne 0.010 0.041 Na 0.010 0.038SHZ2 N max
14 20 O 0.044 0.055 Ne 0.020 0.065TABLE IX. Correlation entropies in the TDDM ground stateof all the isotopes considered in our work for two values of N max and two Skyrme parametrizations. symmetry restrictions, following well-established precur-sors in TDHF [17]. We also use a self-consistent interac-tion, both at the mean-field level and at the residual one.To this end, we work with density-independent Skyrmeinteractions to avoid any issues related to rearrangementterms. We compute the ground states using a dynamicalTDDM code. Our starting point is the corresponding HFground state. We then switch-on beyond mean-field cor-relations adiabatically, by slowly incorporating beyond-mean-field terms over time.With this approach, we investigate light nuclei withmasses ranging from A = 12 to A = 24. Our approachcan tackle closed-shell nuclei, like C and O, but alsoopen-shell isotopes, like − Ne, − Na or O. We findcorrelated ground states for all these isotopes. We studythe effect of TDDM correlations using a variety of quan-tities, including single-particle energies and occupations,but the main focus of our analysis is on binding ener-gies. For the majority of nuclei, the correlation energyaccounts for ≈ C, where two-body correlationsare significantly stronger and account for ≈
11% of thetotal energy. A quantitative metric based on the corre-lation entropy provides similar results. Where the com-parisons are possible, our results provide qualitativelysimilar predictions to previous TDDM implementations.We also confirm a trend of diminishing correlations whenthe neutron number increases.We find two key limitations in this initial study, thatcould be improved in the future. On the one hand, theSkyrme parametrizations that we have used are relativelypoor. Among other things, they have not been fitted tothis mass region or to account for beyond mean-field cor-relations and, in this sense, our predictions can only indi-cate qualitatively the size of correlations in these nuclei.On the other hand, we find that our adiabatic switch-ing produces final states that may appear static whenit comes to one observable, like the energy, but are notstationary in others, like the correlation entropy. This in-dicates that longer evolution times are required, although this comes at a significant larger numerical cost.This work opens several potential avenues for imme-diate future work. When it comes to computing groundstate properties, we have demonstrated that TDDM pro-vides a stable description of relatively light nuclear sys-tems. The extension to higher mass numbers is straight-forward, if numerically challenging. One could fur-ther characterise these TDDM ground-state by analysingtheir cluster structure or by exploiting the connectionsbetween TDDM and different pairing approximations.Furthermore, time-dependent techniques are particularlysuitable for the study of excitations on top of theseground states. It may be interesting to excite and time-evolve different modes using TDDM, to test the validityof mean-field approaches but also to identify correlationeffects on resonances. Finally, dynamical simulations canalso tackle nuclear collisions of interest for a variety ofapplication, including heavy-ion fusion reactions [49].We can also envisage some additional formal devel-opments that could be useful in the context of nuclearphysics. One could attempt to overcome the limitationsassociated to density-independent forces by extendingthe TDDM formalism to include genuine three-nucleoninteractions. This is presumably challenging, since theBBGKY hierarchy would likely have to be modified. Thetreatment of genuine three-body correlations may be rel-evant in nuclear systems too [43]. To break away fromthe adiabatic evolution picture, one could also attempt todevise an energy minimisation process that included two-body density matrices [79–82]. By implementing beyondmean-field simulations, like those presented here, to-gether with the aforementioned developments, one wouldopen the door to a truly first-principles understanding ofthese applications.
ACKNOWLEDGMENTS
This work was supported by the Polish NationalScience Centre (NCN) under Contract No. UMO-2016/23/B/ST2/01789; by STFC, through Grants NosST/M503824/1, ST/L005743/1 and ST/P005314/1; andby the Spanish State Agency for Research of the SpanishMinistry of Science and Innovation through the “Ram´ony Cajal” program with grant RYC2018-026072 and the“Unit of Excellence Mar´ıa de Maeztu 2020-2023” awardto the Institute of Cosmos Sciences (CEX2019-000918-M). This work used the DiRAC Complexity system, op-erated by the University of Leicester IT Services, whichforms part of the STFC DiRAC HPC Facility. Thisequipment is funded by BIS National E-Infrastructurecapital grant ST/K000373/1 and STFC DiRAC Opera-tions grant ST/K0003259/1. DiRAC is part of the Na-tional E-Infrastructure.4
Appendix A: th order Runge-Kutta timestepimplementation In typical implementations of TDHF, the time-stepping procedure involves an integration via the mid-point method. We have found that the solution of theTDDM equations necessarily requires time stepping algo-rithms that provide more accurate results for values of dt ,the time-step size, that are not prohibitively small. Wehave therefore implemented an explicit 4 th order Runge-Kutta (RK4) algorithm to solve the set of coupled differ-ential equations of relevance. We note that RK4 carriesa cumulative error of order dt [83].In the case of TDDM, one has a set of three differentialequations for the evolution of the single particle orbitals,Eq. (2); occupations, Eq. (5); and correlation tensor el-ements, [Eq. (6). We can schematically write this set ofequations as follows: dψdt = P ( t, ψ, n, C ) , (A1) dndt = N ( t, ψ, n, C ) , (A2) dCdt = C ( t, ψ, n, C ) . (A3)Using a RK4 algorithm, given the initial conditions( t p , ψ p , n p , C p ), the estimates for the functions at t p +1 read: ψ p +1 = ψ p + dt (cid:34) d + 2 d + 2 d + d (cid:35) , (A4) n p +1 = n p + dt (cid:34) e + 2 e + 2 e + e (cid:35) , (A5) C p +1 = C p + dt (cid:34) f + 2 f + 2 f + f (cid:35) . (A6) The coefficients d · · · d are given by the following 4equations evaluated either at the initial step, at the mid-points or at the final one: d = P ( t p , ψ p , n p , C p ) (A7) d = P (cid:18) t p + dt , ψ p + dt d , n p + dt e , C p + dt f (cid:19) (A8) d = P (cid:18) t p + dt , ψ p + dt d , n p + dt e , C p + dt f (cid:19) (A9) d = P ( t p + dt, ψ p + dt d , n p + dt e , C p + dt f ) . (A10)The remaining coefficients e · · · e and f · · · f are foundanalogously using the replacements P → N and
P → C ,respectively.We note that, in addition to the single particle orbitals,occupation matrices and correlation tensors, other auxil-iary quantities, such as densities and mean-fields, need tobe recalculated in the 4 steps involving d i , e i and f i . Thisguarantees the stability of the RK4 algorithm within theTDDM method. In all the simulations performed in thiswork, we found that dt ≤ . − provided acceptableand numerically stable results. [1] P. Ring and P. Schuck, The Nuclear Many-Body Problem ,Physics and astronomy online library (Springer-Verlag,Berlin, 1980).[2] J.-P. Blaizot and G. Ripka,
Quantum Theory of FiniteSystems (The MIT Press, Cambridge, 1986).[3] C. Simenel, B. Avez, and D. Lacroix, “Micro-scopic approaches for nuclear Many-Body dynam-ics: Applications to nuclear reactions,” (2008),https://arxiv.org/abs/0806.2714, arXiv:0806.2714 [nucl-th].[4] C. Simenel, Eur. Phys. J. A , 1 (2012).[5] Y. Engel, D. Brink, K. Goeke, S. Krieger, and D. Vau-therin, Nucl. Phys. A , 215 (1975).[6] P. Bonche, S. E. Koonin, and J. Negele, Phys. Rev. C , 1226 (1976).[7] K.-H. Kim, T. Otsuka, and P. Bonche, J. Phys. G: Nucl.Part. Phys. , 1267 (1997). [8] T. Nakatsukasa, K. Matsuyanagi, M. Matsuo, and K. Ya-bana, Rev. Mod. Phys. , 045004 (2016).[9] A. Bulgac, Ann. Rev. Nucl. Part. Sci. , 97 (2013).[10] A. Bulgac, P. Magierski, K. J. Roche, and I. Stetcu,Phys. Rev. Lett. , 122504 (2016).[11] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965).[12] A. S. Umar and V. E. Oberacker, Phys. Rev. C ,054607 (2006).[13] I. Stetcu, A. Bulgac, P. Magierski, and K. J. Roche,Phys. Rev. C , 051309 (2011).[14] G. Scamps and D. Lacroix, Phys. Rev. C , 034314(2014).[15] G. Scamps, C. Simenel, and D. Lacroix, Phys. Rev. C , 011602 (2015).[16] P. D. Stevenson, E. B. Suckling, S. Fracasso, M. C. Bar-ton, and A. S. Umar, Phys. Rev. C , 054617 (2016).[17] J. Maruhn, P.-G. Reinhard, P. Stevenson, and A. Umar, Comp. Phys. Comm. , 2195 (2014).[18] Y. Iwata and P. Stevenson, New J. Phys. , 043010(2019).[19] T. Nishikawa, Y. Iwata, and S. Chiba, Proceedings ofthe 15th International Conference on Nuclear ReactionMechanisms, NRM 2018 , 165 (2018).[20] P. D. Stevenson and J. L. Willerton, SciPost Physics Pro-ceedings , 047 (2020).[21] K. Godbey and A. S. Umar, Front. Phys. , 40 (2020).[22] G. Stefanucci and R. van Leeuwen, NonequilibriumMany-Body Theory of Quantum Systems: A Modern In-troduction (Cambridge University Press, 2013).[23] R. Balian and M. V´en´eroni, Phys. Rev. Lett. , 1353(1981).[24] J. M. A. Broomfield and P. D. Stevenson, J. Phys. G:Nucl. Part. Phys. , 095102 (2008).[25] C. Simenel, Phys. Rev. Lett. , 112502 (2011).[26] P. Danielewicz, Ann. Phys. (N. Y). , 305 (1984).[27] H. S. K¨ohler, N. Kwong, and H. a. Yousif, Comput.Phys. Commun. , 123 (1999).[28] M. H. Mahzoon, P. Danielewicz, and A. Rios, in AIPConf. Proc. , Vol. 1912 (2017) p. 020012.[29] H. Lin, H. Mahzoon, A. Rios, and P. Danielewicz, Ann.Phys. , 168272 (2020).[30] W. Shun-jin and W. Cassing, Ann. Phys. , 328(1985).[31] M. Tohyama, Phys. Rev. C , 187 (1987).[32] M. Gong and M. Tohyama, Z. Phys. A , 153 (1990).[33] M. Tohyama, Prog. Theor. Phys. , 905 (1994).[34] M. Tohyama, Prog. Theo. Phys. , 1293 (1998).[35] M. Tohyama, Phys. Rev. C , 2603 (1998).[36] M. Tohyama, Nucl. Phys. A , 343 (1999).[37] M. Tohyama and A. S. Umar, Phys. Rev. C , 037601(2002).[38] M. Tohyama and A. Umar, Phys. Lett. B , 72 (2002).[39] M. Tohyama, Phys. Rev. C , 054330 (2013).[40] M. Tohyama, Phys. Rev. C , 017301 (2015).[41] M. Tohyama and A. S. Umar, Phys. Rev. C , 034607(2016).[42] M. Tohyama and P. Schuck, Eur. Phys. J. A , 74(2019).[43] M. Tohyama, Front. Phys. , 10 (2020).[44] M. Gong, Time-Dependent Density-Matrix Theory ,Ph.D. thesis, Michigan State University. (1990).[45] W. Cassing and A. Pfitzner, Z. Phys. A , 175 (1990).[46] W. Cassing and A. Pfitzner, Z. Phys. A , 161 (1992).[47] M. Assi´e and D. Lacroix, Phys. Rev. Lett. , 202501(2009).[48] M. Assi´e,
Influence des corr´elations entre les nucl´eonssur les r´eactions de cassure nucl´eaire: aspects th´eoriqueset exp´erimentaux , Ph.D. thesis, Universit´e de Paris-Sud,Paris 11 (2008).[49] K. Wen, M. C. Barton, A. Rios, and P. D. Stevenson,Phys. Rev. C , 014603 (2018).[50] K. Wen, M. C. Barton, Barton, A. Rios, and P. D.Stevenson, Act. Phys. Pol. B , 567 (2019).[51] E. Gross, E. Runge, and O. Heinonen, Many-ParticleTheory, (Taylor & Francis, 1991).[52] A. Akbari, M. J. Hashemi, A. Rubio, R. M. Nieminen,and R. van Leeuwen, Phys. Rev. B , 235121 (2012).[53] A. Akbari, Development And Applications of Time-Dependent Density Matrix Functional Theory , Ph.D. the- sis, University of the Basque Country (2012).[54] H. P. Breuer and F. Petruccione,
The theory ofopen quantum systems (Oxford University Press, GreatClarendon Street, 2002).[55] M. C. Barton,
Symmetry Unrestricted Self Consis-tent Time Dependent Density Matrix Theory with aSkyrme force , Ph.D. thesis, University of Surrey (2018),http://epubs.surrey.ac.uk/id/eprint/845812.[56] M. Tohyama and P. Schuck, Eur. Phys. J. A , 77(2014); Eur. Phys. J. A , 186 (2017).[57] K. J. Schmitt, P. G. Reinhard, and C. Toepffer, Z. Phys.A , 123 (1990).[58] T. Gherega, R. Krieg, P.-G. Reinhard, and C. Toepffer,Nucl. Phys. A , 166 (1993).[59] M. Beiner, H. Flocard, N. V. Giai, and P. Quentin, Nucl.Phys. A , 29 (1975).[60] W. Satu(cid:32)la, J. Dobaczewski, W. Nazarewicz, and T. R.Werner, Phys. Rev. C , 054316 (2012).[61] T. Skyrme, Nucl. Phys. , 615 (1958).[62] D. Vautherin and D. M. Brink, Phys. Rev. C , 626(1972).[63] M. Gell-Mann and F. Low, Phys. Rev. , 350 (1951).[64] B. Schuetrumpf, P.-G. Reinhard, P. D. Stevenson, A. S.Umar, and J. A. Maruhn, Comput. Phys. Commun. , 69 (2013).[67] M. Freer, H. Horiuchi, Y. Kanada-En’Yo, D. Lee, andU. G. Meißner, Rev. Mod. Phys. , 35004 (2018).[68] J.-P. Ebran, E. Khan, T. Niksic, and D. Vretenar, Nature , 341 (2012).[69] M. A. Schumaker et al. , Phys. Rev. C , 044321 (2008).[70] M. A. Schumaker et al. , Phys. Rev. C , 044325 (2009).[71] J. M. D’Auria et al. , Phys. Rev. C , 065803 (2004).[72] M. Bentley and S. Lenzi, Prog. Part. Nuc. Phys. , 497(2007).[73] R. Kanungo et al. , Phys. Rev. Lett. , 152501 (2009).[74] C. R. Hoffman et al. , Phys. Rev. Lett. , 152502(2008).[75] M. D. Jones et al. , Phys. Rev. C , 051306 (2015).[76] K. Tshoo et al. , Phys. Rev. Lett. , 022501 (2012).[77] H. Appel and E. K. U. Gross, EPL (Europhysics Letters) , 23001 (2010).[78] P. Ziesche, J. Mol. Struc.-THEOCHEM , 35 (2000).[79] C. Garrod and J. K. Percus, J. Math. Phys. , 1756(1964).[80] C. Garrod, M. V. Mihailovi´o, and M. Rosina, J. Math.Phys. , 868 (1975).[81] B. Verstichel, PhD Thesis: Variational determination ofthe two-particle density matrix as a quantum many-bodytechnique , Ph.D. thesis (2012).[82] W. Poelmans,
Variational determination of the two-particle density matrix: the case of doubly-occupied space ,Ph.D. thesis, Ghent University (2015).[83] E. S¨uli and D. F. Mayers,