{'fraction_non_alphanumeric': 0.03334583518019708, 'fraction_numerical': 0.011047086501414982, 'mean_word_length': 4.676407973679117, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computational approaches to many-body dynamics of unstable nuclear systems

The goal of this presentation is to highlight various computational techniques used to study dynamics of quantum many-body systems. We examine the projection and variable phase methods being applied to multi-channel problems of scattering and tunneling; here the virtual, energy-forbidden channels and their treatment are of particular importance. The direct time-dependent solutions using Trotter-Suzuki propagator expansion provide yet another approach to exploring the complex dynamics of unstable systems. While presenting computational tools, we briefly revisit the general theory of the quantum decay of unstable states. The list of questions here includes those of the internal dynamics in decaying systems, formation and evolution of the radiating state, and low-energy background that dominates at remote times. Mathematical formulations and numerical approaches to time-dependent problems are discussed using the quasi-stationary methods involving effective Non-Hermitian Hamiltonian formulation.


Introduction
There is no physical system that is truly isolated from the rest of the world, the closed system idealization may be convenient but becomes poor or completely invalid for many questions of modern-day science. In nuclear physics, as interests shift towards weakly bound, unbound or even dynamically evolving reaction states, the theoretical approaches for dealing with unstable dynamics of open quantum systems with multiple degrees of freedom must be revisited. The availability of advanced computational technologies calls forth innovative thinking and new philosophies in addressing these types of quantum many-body problems. In this presentation, using different models and realistic examples from the world of nuclear physics, we discuss computational strategies and techniques for dealing with dynamically unstable many-body systems. The Nuclear Theory in the Supercomputing Era venue is especially timely and allows us to put emphasis on some of the techniques, that due to their computational nature, remained behind the curtains in a number of recent investigations [1][2][3]. see also Refs. [2,4]. Consider a model of scattering illustrated in Fig. 1. In this one-dimensional problem two particles with masses µ 1 and µ 2 comprise a composite system of unit mass µ 1 +µ 2 = 1. The system can be described with the center-of-mass and relative coordinates, X = µ 1 x 1 + µ 2 x 2 and x = x 1 − x 2 , respectively. The two particles are confined by a potential v(x). The intrinsic Hamiltonian is assumed to have a complete set of discrete eigenstates ψ n (x) with corresponding intrinsic energies n : hψ n (x) = n ψ n (x), n = 0, 1, 2, ...
Here the reduced mass is µ = µ 1 µ 2 and we select our units so that 2 /2 = 1. We assume that this system scatters off an infinite wall and the wall interacts only with the second particle. Therefore the full Hamiltonian is As illustrated in Fig. 1, we assume that the incident beam is traveling from the left and contains the projectiles in an intrinsic state (channel) n. A complete set of reflected waves is characterized by the amplitudes R nm defined here so that |R mn | 2 represents the probability for the initial beam in channel n to reflect in channel m; R nm = R mn due to time-reversal invariance. The scattering wave function is where is the center-of-mass momentum of the two-particle system while in the n th intrinsic state, and E is the total energy. Figure 1: Schematic picture of scattering.
A composite system of two particles bound by a harmonic oscillator potential scatters off an infinite wall. One of the particles does not interact with the wall, at the same time the wall is impenetrable for the second particle.
A channel n is considered to be open if E ≥ n and the corresponding momentum K n is real. The conservation of particle-number in all the open channels necessitates m∈open |R mn | 2 = 1. The channel is closed if E < n , in which case K n is purely imaginary. We stress that the principal value of the square root is implied in Eq. (4).
The boundary condition set by an impenetrable wall is to be used for determining the set of coefficients R mn . Since at x 2 = 0 the center-ofmass coordinate X = µ 1 x, the boundary condition can be expressed in the intrinsic coordinate x only Φ(µ 1 x, x) = 0. Therefore we can project the reaction problem onto a complete set of intrinsic basis states, which leads to the following linear equation where the matrix D is defined as Eq. (6) represents a typical mathematical challenge associated with the formulation of reaction problems where reaction states are projected onto the intrinsic states; see also Sec. 3. It is a linear algebra problem where the construction of the scattering matrix amounts to matrix inversion in the projected space. The scattering energy E is a running parameter here, and studies of scattering at different energies is therefore time consuming. And, finally, the underlying matrix is highly singular and there are issues with convergence. The latter difficulty is the one that we would like to illustrate using this example. If the two particles forming a composite system are bound by a harmonic oscillator confinement, v(x) = µω 2 x 2 /2 in Eq. (1), the D-matrix is then known analytically [2]. Then to solve the problem we truncate the channel space at some large number N of oscillator quanta, and solve Eq. (6) using standard numerical techniques. This turns out to be a difficult task; the matrix element D mn (κ) for virtual channels, where κ is real, are exponentially large, making the process of matrix inversion difficult and numerically unstable [2,5,6]. As shown in Fig. 2, left panel, the absolute values of the reflection amplitudes, R n ≡ R 0n exponentially diverge for increasingly remote virtual channels.
While it is possible to overcome the numerical issues, further examination shows that the approach has fundamental flaws. In Fig. 2, right panel, the phase shift, defined as e 2iδ = −R 00 , is shown as a function of N . While satisfactory and seemingly convergent results can be easily found for the cases where the mass of the noninteracting particle is small, in general, as N increases, the results start oscillating; situations where the non-interacting particle is heavy and therefore deeply penetrates the wall are particularly difficult to handle. It was emphasized in Refs. [2,4] that there is no numerical convergence with increasing N .

Variable Phase Method
The above example shows that reaction problems call for new techniques. One approach, based on the Variable Phase Method (VPM), see Ref. [7], is proposed in Ref. [2]. The VPM is an effective technique for solving the coupled-channel problem of the form where scattering observables are to be expressed relative to free-space solutions normalized to unit current the ± sign corresponds to a wave moving in the right/left direction. In the VPM approach the coupled-channel Schödinger's equation (8) is reformulated as a set of first order differential equations for dynamic reflection and transmission amplitude matrices R nn (X ) and T nn (X ). These amplitudes correspond to a potential that is cut at X , namely to V nn (X)θ(X − X ) :  , is plotted as a function of truncation N. The problems with the approach are highlighted by an unstable and oscillatory behavior of the phase shifts. The problem is particularly severe when the non-interacting particle of mass µ 1 is heavy. Different curves show phase-shifts for different mass ratios µ 1 /µ 2 = 1, 2, 3, 5, 10, as labeled; the exact values obtained with Variable Phase Method (see Sec. 2.2) are shown by the horizontal grid lines with the tic-marks on the right. Inset shows the case when µ 1 /µ 2 = 3 extending the study to considerably large values of N and emphasizing that for any choice of parameters the approach fails at some point.
These equations, being solved from X = +∞ towards X → −∞, recover the reflection and transmission amplitudes R nn (−∞) = R nn and T nn (−∞) = T nn . Using factorization of the form the Schrödinger's equation for the scattering problem described in Fig. 1 can be transformed into a coupled-channel equation (8) for the center-of-mass wave-functions Ψ n (X) where the folded potentials are Some representative results for the scattering problem where an oscillator-bound system interacts with an infinite wall are shown in Figs Probability density X µ 1 , E _____ 0.5, 0.5 0.9, 0.5 0.5, 7.5 0.9, 7.5 Figure 4: The density of probability for the center of mass of the projectile to be at a location X when it is reflected from an infinite wall at X = 0.
energy. The kinetic energy is expressed in units of oscillator's ω and therefore for each integer value thereof a new channel opens. One can notice typical cusps at thresholds associated with the loss of flux into newly opened channels. In Fig. 4 the probability distribution for the center of mass is shown. The four curves show four of the most representative situations; low and high incident kinetic energies E = 0.5 ω and E = 7.5 ω, respectively, and two different mass-ratios µ 1 = 0.5 and 0.9.

Time-dependent approach
Turning to a time-dependent approach is a natural strategy for dealing with nonstationary systems. There are various computational techniques; see Ref. [8] for some recent tests and comparisons of methods being applied to one dimensional Schrödinger's equation. In time-dependent techniques preservation of unitarity is often at the core of computational difficulties: lack of unitarity could lead to exponential amplification of numerical noise even for single channel, while in multi-channel problems discontinuities near thresholds are particularly challenging. Here we propose and demonstrate another approach that is computationally efficient, even in multi-variable cases, and preserves unitarity exactly. The time propagation can be performed by considering, separately, the potential and kinetic parts of the can be readily applied. Similarly, in the conjugate momentum space p = {p 1 , p 2 . . . } the propagation with kinetic energy operator, which is diagonal, is also easy to perform. While the operators K and V do not commute, the time evolution (13) with the combined Hamiltonian can be done efficiently with the Trotter-Suzuki approach [9,10].
In this approach the propagation is done in small time steps ∆T ; for each of these steps the evolution operator is approximated as The Fast Fourier Transform allows for an efficient transition between coordinate and momentum representations so that exponentials of operators are always applied in the diagonal form. Even with the finite time steps the unitarity is fully retained; the method is applicable to time-dependent Hamiltonians. The computational cost of two back and forth Fourier transforms involved in each step is N log(N ) assuming the coordinate space is discretized into N points. While this at first appears to be higher than the typical O(N ) scaling of the traditional methods, in practice the cost cN of any high quality method involves a constant factor c that often exceeds log(N ). Moreover, modern computer hardware often comes with signal processing tools which are optimized at hardware and software level to perform Fast Fourier Transform with incredible efficiency.
Let us return to the problem of scattering illustrated in Fig. 1. The time dependent picture of the scattering process is shown in Fig. 5 with a series of four plots showing the two-dimensional wave function using a density plot at four different times. The plot of the density projection onto the center of mass coordinate X, which is the time dependent analog of Fig. 4, is shown below each of the four snapshots. The initial wave function at t = 0, shown on the first panel, is selected as the ground state wave function for the intrinsic potential, and as a moving Gaussian wave packet for the center of mass coordinate, In this example σ 0 = 2, X 0 = −5, and initial momentum K 0 = 1, all quantities are being expressed here in dimensionless units of distance as defined earlier. While this time dependent consideration is different from the stationary state formulation studied above, the series of snapshots for different times shown in Fig. 5 highlights some similar features. At high energies the dynamics of virtual excitations is complex; this is illustrated in Fig. 6, where the initial wave packet is selected to have K 0 = 5. Some semiclassical interpretation can be given to the stages of the process. Initial compression at t = 1 is followed by two particles bouncing apart at t = 2. Having equal masses, their center of mass remains at the origin but the relative separation x becomes large so that the particles are positioned roughly symmetrically on the opposite sides of the wall. Next at t = 3 the center of mass moves into X < 0 region pressing the interacting particle A.Volya against the wall. Finally, the system is reflected at t = 4 with the initial wave packet being considerably distorted.
In comparison to the projection and VPM techniques discussed earlier, the timedependent approach is substantially faster numerically; moreover, any potential U (x 1 , x 2 ) can be considered with ease in this approach. One has to keep in mind, however, that it is not always easy to provide quantitative answers to stationary state questions, such as determination of scattering phase shifts in this example, using time-dependent techniques. The exact choice of the initial state as well as the energy uncertainty of the initial state can be important for some stages of time evolution.
The physics of decay of unstable states represents a particularly important class of time-dependent process. The familiar exponential decay law is only an incomplete picture, requiring some subtle approximations, and being valid only within certain time limits. The complex intrinsic dynamics that can occur in the decaying manybody system further complicates the time evolution. The non-exponential decay laws in quantum mechanics have been studied and revisited by many authors ( see Ref. [3] and references therein). The presence of three regimes, namely, initial, exponential, and long-time power law, appears to be a universal feature of the decay processes. The transitions from one regime to another are accompanied by the interference of corresponding quantum amplitudes that is seen as oscillations on the decay curve.
As another demonstration of the time-dependent technique based on the Trotter-Suzuki expansion and as an introduction to the section that follows, we demonstrate in Fig. 7 the decay process in Winter's model [11], which has been a very popular tool for exploring non-exponential features in decays. In this model a particle is confined to a region x ≥ 0 by an impenetrable wall at x = 0 and is held by a delta barrier at x = 1. The initial state at t = 0 is taken as Ψ(x, 0) = √ 2 sin(πx). The survival probability shown in the solid red line on the left panel of Fig. 7 illustrates the three general regimes: pre-exponential, exponential, and post-exponential. Oscillations can be seen in transitional regions. The snapshots of the wave function at different times are shown on the right.
The pre-exponential behavior at very early times is influenced by the memory of how the state was created and, in particular, by the high energy components in the state. Later in time the internal structure and transitions between the intrinsic states become relevant. Short times correspond to remote energy components where the presence of other resonant states is to be considered. The high energy components have much shorter lifetimes and decay quickly leading to an exponential decay phase. This phase is dominated by a single resonant component, the radiating state, so that the wave function retains its shape while decreasing in amplitude. This can be seen on the right panel of Fig. 7. In the same figure one can also trace a moving away background component. The background contains very low energy particles; being far off-resonance, they essentially do not interact but move slowly away from the interaction region. Near the decay threshold the number of such particles with a certain energy is determined by the available phase space, which for neutral particles scales with energy following a power-law E l+1/2 where l is the angular momentum quantum number. This type of scaling leads to non-resonant components that follow a power-law decay S(t) ∼ 1/t 2 +3 . While the non-resonant component can be very small in the initial state, eventually it becomes dominant due to its slower-thanexponential decay. Further discussion of decay processes in quantum mechanics and other examples can be found in Ref. [3]. The near-threshold phase space scaling with energy which leads to power-law decay at remote times is an important consideration in the Time Dependent Continuum Shell Model approach that is discussed in the following section, see also Refs. [1,12,13], as well as in more complicated sequential decay processes [14].
3 Time dependent continuum shell model 3

.1 Continuum Shell Model
A seamless transition between structure physics and reactions is one of the central present-day theoretical problems. The computational aspect associated with transitions from discrete levels to a continuum of reaction states is especially challenging. The Continuum Shell Model approach [12,13] and its time-dependent version, in particular, is one among several theoretical tools confronting these issues. In the Continuum Shell Model the Feshbach projection formalism [15,16] is used to express the exact dynamics in the full Hilbert space using an effective Hamiltonian in the projected intrinsic subspace of interest, Q : Here the effective Hamiltonian contains H QQ which is the part of the original Hamiltonian that acts in the space Q, and the energy-dependent non-Hermitian termH(E), that emerges from the coupling of the space Q to an external space containing a continuum of reaction states, P.
In practical applications the intrinsic space Q is assumed to represent the configuration space of the traditional shell model, built from states |1 that are Slater determinants constructed from bound-state single-particle wave functions. The space P contains continua of reaction states |c, E characterized by the channel index, c, and the continuous energy parameter E. There is a certain threshold energy E (c) thr for each channel c. The energy-dependent non-Hermitian effective Hamiltonian (16) is then represented by a matrix H 12 (E) ≡ 1|H(E)|2 , and the channel amplitudes are the matrix elements A c 1 (E) = 1|H|c, E . The traditional shell model Hamiltonian is recovered when the internal space Q is isolated and thus is decoupled, A c 1 (E) = 0. The computational challenges of the traditional shell model approach are well known, they are mainly associated with the need to find some selected eigenvalues and eigenvectors of the Hamiltonian matrix H 12 . The matrix is generally sparse, thanks to few-body nature of the underlying nucleon-nucleon interactions which inhibits mixing of very remote configurations, thus iterative techniques such as Lanczos approach are commonly used.
The physics of weakly-bound and unstable nuclear systems is much more rich as questions of interest span from properties of bound states to features in scattering cross sections. Narrow resonances are well characterized by the usual properties of bound states with the decay width being an additional characteristic. This requires the non-Hermitian eigenvalue problem H(E)|I = E|I to be solved. The resulting complex energies E represent positions of resonances, E = Re(E), and their widths, Γ = −2 Im(E). The most practical technique here is to start with the perturbative treatment and evaluate the termH(E) associated with continuum, using the wave functions of the traditional shell model Hamiltonian H QQ . As coupling to the continuum increases the states become broad and one is forced to treat the non-Hermitian energy-dependent eigenvalue problem as an iterative non-Hermitian diagonalization process. In this limit a major problem is associated with the physical interpretation of the resonances and their widths.
Formally, the energy-dependent non-Herminitan Hamiltonian provides an exact propagator for the intrinsic space and therefore the scattering matrix is Here ξ c is a potential (direct-reaction) phase. The matrix is unitary (see Ref [1]) and the unitarity is related to a factorized form of the imaginary W 12 (E) in Eq. (17). The eigenvalues of the non-Hermitian Hamiltonian are therefore poles of the scattering matrix. In the limit of broad resonances one has to address the reaction problem where obtaining a reaction cross section is the main goal. There are several numerical challenges associated with Eq. (18), many of these challenges being similar to the ones discussed in Sec. 2.1. First, the size of the Hamiltonian matrix and the complex arithmetic involved are not making this problem simpler as compared to matrix diagonalization. Second, the scattering energy E represents a running parameter so that the procedure should be repeated for all energies of interest. Finally, the problem is numerically unstable: bound states, as well as resonances with widths ranging by many orders of magnitude, may be encountered and should be treated consistently. All of these technical issues are resolved by the Time-dependent Continuum Shell Model approach which we discuss next.

Time-dependent many-body evolution operator
The many-body wave function follows the time evolution which is a Fourier image of the retarded propagator involved in the scattering matrix (18): Here H is an arbitrary Hamiltonian, but as discussed below, it is advantageous to include a factorized imaginary part W using a different procedure described in Sec. 3.3. Thus, we view H as being a Hermitian Hamiltonian of the traditional shell model in which case it is set to have an infinitesimal negative-definite imaginary part. The time-dependent evolution operator can be factorized using a Chebyshev polynomial expansion method, see Ref. [1,17,18].
where J n is the Bessel function of the first kind and T n represents Chebyshev polynomials. The Chebyshev polynomials, defined as T n [cos(θ)] = cos(nθ) or, in explicit form, provide a complete set of orthogonal functions covering uniformly the interval [-1, 1]. In contrast, Taylor expansion relies on power functions which favor the edges of the interval and thus are more sensitive to extreme eigenvalues. The "angular addition" identity which follows from the definition, allows one to obtain these polynomials using the recurrence relation Therefore, the process of evaluation of Chebyshev polynomials of the Hamiltonian operator is an iterative procedure, similar to the one in Lanczos approach. For a given initial state |λ ≡ |λ 0 , a sequence |λ n = T n (H)|λ can be constructed as |λ 0 = |λ , |λ 1 = H|λ , and |λ n+1 = 2H|λ n − |λ n−1 .
For overlap functions, assuming Hermitian H, one can also use the following identity  (20) needed for convergence is denoted as n max (τ ). The asymptotic of Bessel functions J n (x) ≈ 1/(2πn)[ex/(2n)] n suggests n max (τ ) ≈ eπτ /2 ≈ 4τ. At fixed values of n but for large times the convergence remains stable due to J n (t) ≈ 2/(πt) cos(t − πn/2 − π/4) in this limit. For the desired energy resolution ∆E/N the propagation in time has to be extended up to ≈ τ N which requires n max ≈ 4N ; therefore 2N matrix-vector multiplications are required if one also uses Eq. (25).
The time-dependent approach provides the Green's function for all energies at once; it is also exceptionally stable numerically when dealing with very narrow resonances or with stable states. Indeed, the time-dependent behavior of stationary states is regular and the corresponding delta function in energy is well handled by Fourier transform, which at the desired energy resolution properly conservers the integrated strength.
In order to illustrate the approach, let us consider strength and integrated strength functions defined for a given state |λ as In Fig. 8 both strength (left) and integrated strength (right) functions are shown for 15 N for neutron channels where |λ corresponds to different angular momentum channels constructed from 1 + ground state in 14 N coupled to a single nucleon on either d 5/2 (top panels) or d 3/2 (bottom panels) single-particle states. This theoretical study follows recent experimental work in Ref. [19]. The full p-sd valence space is used with the Hamiltonian from Ref. [20]. With about 10 7 m-scheme basis states, obtaining and computing strength functions in energy regions around 20 MeV of excitation is impractical; the time-dependent method provides an excellent alternative.

Sherman-Morrison-Woodbury relations
It is certainly possible to implement the Chebyshev polynomial expansion procedure for a full non-Hermitian Hamiltonian using Eq. (20); however the factorized structure ofH offers a different alternative which is much more computationally advantageous. The two propagators corresponding to Eq. (16) can be related through Dyson's equation Since the contribution from the continuum emerges in the factorized form the expression for the full propagator can be found in a closed form in the space spanned by the channel states The operators here are represented by matrices in the channel subspace G ab = a|G(E)|b and G ab = a|G(E)|b . In computer science these relations are known as Sherman-Morrison-Woodbury matrix inversion equations [21]. The unitarity of the scattering matrix immediately follows from these relations, see [1]. We illustrate the TDCSM approach in its complete form in Figs. 9 and 10 where the resonances in 24 O are considered. The system is treated in the sd valence space using the USD shell model Hamiltonian [22]. In Fig. 9 the norm of the survival amplitude is shown as a function of time for the following set of most representative states 2 + 1 (4180, 2.7), 1 + 1 (5291, 195.1), 4 + 1 (6947, 0.0), 2 + 3 (8107, 92.5), and 2 + 4 (9673, 17.5). The states are listed here with their excitation energies followed by the decay widths, both in keV. The initial wave functions at t = 0 are taken as eigenstates of the traditional shell model. For the states such as 4 + 1 , which cannot decay in this model due to high angular momentum, the norm of the survival amplitude remains constant. Narrow states, exhibit a nearly exponential decay, for the state 2 + 4 the survival amplitude expected in exponential decay is shown. The decay is non-exponential for broad states such as 1 + 1 and 2 + 3 . In Fig. 10   . They are eigenstates of the traditional USD SM but are non-stationary resonances in the TDCSM, except for the 4 + 1 state which due to its high spin does not decay within the sd valence space. To emphasize the non-exponentiality in the decay law the unmarked solid line shows the exp(−Γ α t/2) function with parameters for the 2 + 4 state.
The time-dependent approach provides an effective computational strategy for treating many-body systems that feature both bound and unbound states. In contrast to the stationary state formalism, the time dependent approach addresses the evolution of states in a natural way, thus providing a computationally robust and stable strategy, where experimental observables are easily recovered and fundamental principles of quantum mechanics, such as linearity and unitarity, are followed. From the computational perspective, the matrix-vector multiplication, the most efficient operation available, is utilized in building the time evolution operator with full con-Computational approaches to many-body dynamics of unstable nuclear systems 15 trol of the desired energy and time resolution. The specifics of the terms that emerge due to coupling to continuum in Feshbach projection formalism can be used to build the full evolution operator using Sherman-Morrison-Woodbury relations. TDCSM found broad practical applications, see Refs. [23][24][25] for example.

Conclusions
As our interests shift towards open, reacting, decaying, and otherwise evolving quantum many-body systems, new theoretical and computational techniques must be developed to address multiple new challenges that emerge. The goal of this presentation is to highlight some of the methods used in the recent scientific projects. We use a simple model to demonstrate three distinctly different techniques. The most straightforward method involves projecting the dynamics onto a set of basis states, allowing subsequently for the well-developed methods of linear algebra to be used; in certain reaction problems this method appears to have significant drawbacks associated with numerical instabilities and poor convergence. We demonstrate the Variable Phase Method that can treat reaction problems efficiently in a discretized coordinate space. Finally, we consider explicitly time-dependent techniques that are perhaps most adequate for the time-dependent dynamics associated with decay. We put forward the Time Dependent Continuum Shell Model approach, as a practical tool and demonstrate its application to realistic problems in nuclear physics.