Time-dependent approaches to open quantum systems

Couplings of a system to other degrees of freedom (that is, environmental degrees of freedom) lead to energy dissipation when the number of environmental degrees of freedom is large enough. Here we discuss quantal treatments for such energy dissipation. To this end, we discuss two different time-dependent methods. One is to introduce an effective time-dependent Hamiltonian, which leads to a classical equation of motion as a relationship among expectation values of quantum operators. We apply this method to a heavy-ion fusion reaction and discuss the role of dissipation on the penetrability of the Coulomb barrier. The other method is to start with a Hamiltonian with environmental degrees of freedom and derive an equation which the system degree of freedom obeys. For this, we present a new efficient method to solve coupled-channels equations, which can be easily applied even when the dimension of the coupled-channels equations is huge.


INTRODUCTION
Open quantum systems are ubiquitous in many branches of science. In general, a system is never isolated but couples to other degrees of freedom, which are often referred to as the environment. The couplings to the environmental degrees of freedom can strongly affect the dynamics of the system. When the number of environmental degrees of freedom is huge, the couplings lead to energy dissipation. It has been demonstrated by Caldeira and Leggett that such couplings suppress the tunneling rate [1,2], going into a transition from quantum to classical regimes. In nuclear physics, it has been well known that a large amount of the relative energy and angular momentum is dissipated during collisions of heavy nuclei at energies close to the Coulomb barrier, known as deep inelastic collisions [3]. In this case, the dissipation occurs due to the couplings between the relative motion of two colliding nuclei and nucleonic degrees of freedom in those nuclei. A classical Langevin equation [4] has been successfully applied to describe such collisions [3]. The Langevin approach has also been employed in order to discuss fusion reactions for syntheses of superheavy elements [5][6][7][8][9][10][11].
The classical Langevin approach, by definition, is not applicable at energies around the Coulomb barrier, at which quantum effects play an important role [12,13]. One can then ask: what is a quantum model that, in the classical Limit, is equivalent to a classical Langevin equation? There are two approaches to address this question. One is to use a phenomenological quantum friction model in which the expectation values of operators obey the classical equation of motion with friction [14][15][16][17]. Recently, we solved such quantum friction Hamiltonians with a time-dependent wave packet approach in order to discuss the effect of friction on quantum tunneling [18]. The other approach to a quantum Langevin equation is to start from a system-plus-bath Hamiltonian, that is, a Hamiltonian that consists of the system and the environmental degrees of freedom and eliminates the environmental degrees of freedom. For instance, one can employ the Caldeira-Leggett Hamiltonian [1] since the classical Langevin equation can be derived from it [3,4]. This approach is more microscopic, and a computation would thus be more involved than the quantum friction model. It has been known that, in the Markovian limit, the time evolution for the reduced density matrix for the system degree of freedom in general takes the so-called Lindblad form [19,20].
In this paper, we discuss both of these approaches for open quantum systems from the point of view of the timedependent method. In the next section, we first discuss the phenomenological quantum friction models using a timedependent wave packet approach. We apply them to heavy-ion fusion reactions around the Coulomb barrier and discuss a role of friction in fusion dynamics. In section III. we solve the Calderira-Leggett Hamiltonian using a time-dependent coupled-channels approach. Using a quantum damped harmonic oscillator, we discuss how one can deal with a large number of degrees of freedom. A summary of the paper is then given in section IV.

PHENOMENOLOGICAL QUANTUM FRICTION MODELS
We first considered a phenomenological approach to quantum friction. In this approach, one treats the environmental degrees of freedom implicitly and introduce a phenomenological Hamiltonian with which the classical equation of motion with a frictional force is reproduced as expectation values. For this purpose, several model Hamiltonians have been proposed so far [14][15][16][17]. Among these, we focused in this paper on the one introduced by Kostin [16].
Consider a particle of mass m moving in a one-dimensional space q under the influence of a potential V(q). With a friction coefficient γ , the phenomenological Schrödinger equation in the Kostin model is given by [16]: where S(q, t) is the phase of the wave function, ψ(q, t) = |ψ(q, t)| exp(iS(q, t)/h). From this equation, it is easy to confirm that one can derive the equation of motion with a frictional force: as is desired. Here, the expectation value of an operator O is denoted as O = dq ψ * (q, t)Oψ(q, t), and p is the momentum operator.
When one simulates an energy dissipation in heavy-ion collisions by means of friction, a frictional force should be active only when the colliding nuclei are close to each other. In other words, one needs to deal with a coordinate-dependent friction coefficient, γ = γ (q). An extension of the Kostin model along this line has been proposed in Immele et al. [21], Bargueno and Miret-Artes [22], with which the modified Schrödinger equation is given by: To apply the phenomenological model to realistic collision problems, one further needs an extension to a three-dimensional space, q = q(r, θ , φ). To this end, we first must expand the wave function with the Legendre polynomials P l (x) as ψ( q, t) = ∞ l=0 u l (r, t)P l (cos θ )/r. One can then modify the Schrödinger equation for u l (r, t) in the same way as Equation (3) to incorporate a frictional force, where S l (r, t) is the phase of the radial wave function u l (r, t) = |u l (r, t)| exp(iS l (r, t)/h). We have here assumed a spherically symmetric potential, V( q) = V(r). Notice that only the radial friction is taken into account here, while the angular momentum dissipation is neglected.
In applying Equation (4) to scattering problems, one needs to use the time-dependent approach since the Hamiltonian depends explicitly on time. To this end, we propagate a wave packet and observe how it bifurcates after it crosses the potential region. Since a wave packet is a superposition of various energy waves, one has to choose the initial condition carefully to get scattering quantities at certain initial energy. Notice that the energy projection approach [23] is inapplicable for our purpose since the energy is not conserved.
In the initial condition, the width of the energy distribution must be small enough to get reasonable results. In this context, the energy refers to the expectation value of the asymptotic Hamiltonian, H 0 . If we restrict to the s-wave scattering, l = 0, and if V(r) rapidly vanishes as r → ∞, one can simply take the kinetic energy operator as H 0 , that is, H 0 = −h 2 /2m ∂ 2 /∂r 2 . The minimum uncertainty wave packet in this case has been discussed in Bracher [24], which reads: where r 0 and σ r are related to the mean position and the width of the wave function in the coordinate space, respectively, and k 0 is related to the mean initial energy.
In heavy-ion collisions, on the other hand, the potential is a sum of the nuclear potential V N and the Coulomb potential V C (r) = Z P Z T e 2 /r with the projectile charge Z P and the target charge Z T . Since the Coulomb potential is a long range potential, the asymptotic Hamiltonian H 0 has to include it: H 0 = −h 2 /2m ∂ 2 /∂r 2 + V C . Thus, the minimum uncertainty wave packet in the form of Equation (5) would not be efficient in this case. Instead, one needs to construct a wave packet from the energy distribution, f C (E), of H 0 in the presence of the Coulomb potential. In analogy to the spherical Bessel functions, we find that this can be achieved as: with E =h 2 k 2 /2m, where η = mZ P Z T e 2 /h 2 k is the Sommerfeld parameter and F 0 (η, kr) is the regular Coulomb wave function.
With the initial condition given by Equation (6), we compute the penetrability of the Coulomb barrier for the 16 O + 208 Pb system in the presence of friction. For the nuclear potential, we employ the optical potential in Evers et al. [25]: .76 fm, and a w = 0.4 fm. With this potential, the Coulomb barrier height V B is 74.5 MeV. For a friction coefficient γ (r), we employ the surface friction model [3], with the Woods-Saxon from ). This is a general form of the friction coefficient obtained perturbatively [26], and it has successfully been applied to the above barrier fusion reactions and to deep inelastic scatterings [3]. We arbitrarily set γ 0 = 4.7 × 10 −23 s/MeV, and this is used in the classical calculations. We compute the phase of the wave function in the same way as in Tokieda and Hagino [18]. For the initial energy distribution in Equation (6), we assume the Gaussian form, where E 0 and σ E are the mean and the width of the initial energy distribution. We have confirmed that σ E = 0.5 MeV is sufficient in the present parameter set. Figure 1 compares the penetrability obtained with and without friction. One finds that the penetrability with friction is shifted to higher energies around the barrier. That is, in the presence of friction, a particle needs additional energy to penetrate the barrier, which originates from the energy dissipation. One can also see that the penetrability does not reach unity at high energies, but it is almost saturated at around 0.9. This means that the exit channel is in a quantum superposition state of absorption and reflection even at sufficiently abovebarrier energies. Notice that, in classical mechanics without a random force, the penetrability can be only 0 or 1. In this sense, this is peculiar to the quantum friction model.
For a practical application to fusion reactions, one needs to take into account explicitly low-lying collective excitations, as they play a crucial role [27]. This can be achieved by extending the above method to the coupled-channels formalism. In that treatment, low-lying collective excitations are taken into account explicitly, while other degrees of freedom, such as non-collective excitations and nucleon transfers, are treated by means of a frictional force.
In Tokieda and Hagino [28], we have applied the method to the 16 O+ 208 Pb system. The nuclear potential is the Woods-Saxon form as Equation (7). For low-lying collective states, the first excited state of both 16 O and 208 Pb were taken into account. The channel-coupling effect was treated in the same way as the CCFULL code, and the iso-centrifugal approximation was adopted [29]. The surface friction model was employed for a friction coefficient [with a replacement of V B with V 0 in Equation (8)], treating γ 0 as an adjustable parameter. See Tokieda and Hagino [28] and reference therein for details of the parameters.
To compute the penetrability at certain energy, we used the same initial condition as Equation (9). It has turned out that σ E = 0.5 MeV was sufficient for the employed parameter set. By taking a sum of the penetrability at each angular momentum, one could calculate fusion cross sections, which could then be compared with experimental data.
In Figure 2, we compare the energy dependence of fusion cross sections. Notice that the present method is reduced to the conventional coupled-channels method when friction is turned off. Thus, the no-friction result is nothing but the result of the conventional coupled-channels approach. To reproduce the experimental data, we set γ 0 = 0.6 × 10 −23 s/MeV. In the left panel of Figure 2, on one hand, one finds that the fusion cross sections at above barrier energies are suppressed in the friction model (the solid line) compared to that in the no-friction model (the dashed line). In the right panel of Figure 2, on the other hand, one sees that the sub-barrier fusion cross sections in both the models give almost the same result. Considering overall energies, the present friction model provided a more consistent description of fusion reactions around the Coulomb barrier as compared to the conventional coupled-channels approach. In this calculation, the same behavior as in Figure 1 has been found, and this may be a key to achieve a consistent description [28].
In comparison with the classical Langevin approach, the present method did not contain a random force originating from thermal fluctuation. For the 16 O+ 208 Pb system, the compound nucleus was formed once the projectile and the target nuclei touched each other, and the fluctuation played a much more minor role compared to the massive systems. In that situation, fusion cross sections were given as an averaged quantity, and it was expected that the presence of a random force does not change the result much. However, when one deals with phenomena in which thermal fluctuation plays a crucial role, such as deep inelastic scattering or a synthesis of superheavy elements, the thermal fluctuation should be explicitly taken into account. One could in fact directly add a random force to the Schrödinger equation Equation (1) as was done in the original paper [16]. Alternatively, friction and fluctuation will naturally emerge by explicitly treating environmental degrees of freedom, which we have discussed in the next section.

SYSTEM-PLUS-BATH MODEL
We next considered a more microscopic model for quantum friction, employing a system-plus-bath model. To be more specific, we considered the Caldeira-Leggett model [1,2], whose Hamiltonian is given by, where H S and H B are the Hamiltonians for the system and the bath degrees of freedom, respectively, while V coup is the coupling Hamiltonian between the system and the bath. Here, the bath degree of freedom is assumed to be a set of harmonic oscillators, whose creation and annihilation operators are denoted by an a † i and a i . The coupling Hamiltonian is assumed to be separable between the system and the bath degrees of freedom. In there, d i is the coupling strength, and h(q) is the coupling form factor, where q is the coordinate of the system. There are several ways to solve the Calderira-Leggett Hamiltonian. In Caldeira and Leggett [1,2], the bath degrees of freedom were integrated out using the path integral in order to obtain an effective action for the system degree of freedom [see also Takigawa and Bertsch [31]]. One can also introduce the influence functional [32]. Here we discussed the coupledchannels approach [33].
In the conventional coupled-channels approach [27], one expands the total wave function in terms of the eigen-wave functions of H B : where the basis states |{n i } are given by: Here, |0 is the vacuum state defined as a i |0 = 0. One can derive the coupled equations for ψ {n i } (q, t) by evaluating the equation: that is, The coupled-channels equations, Equation (15), can be numerically solved when the number of the oscillator modes is not large [27,29]. However, in general, the number of the oscillator modes can be huge, or the distribution of the frequency of the oscillator may even be given as a continuous function. In that situation, it is almost hopeless to solve the coupled-channels equations directly. In order to overcome this problem, we introduced a more efficient basis to expand the total wave function [33]. To this end, we first expanded the function exp(−iωt) with a finite basis set as: where u k (t) is a known function such as a Bessel function, and η k (ω) is the expansion coefficient. We then introduced a new phonon creation operator as: Notice that the number of k is finite, k running from 1 to K, even though the number of i may be infinite. We then constructed the basis states using the operators b † k and expanded the total wave function with them. That is, instead of Equation (12), we expanded the total wave function as: with, One can then obtain the coupled-channels equations similar to Equation (15): We once again emphasize that the dimension of the coupledchannels equations, Equation (20), is much smaller than that of the original equations, Equation (15). The structure of the coupled-channels equations, Equation (20), becomes simple when the basis functions u k (t) satisfy the following two conditions. 1. The matrix D defined as: is diagonal with respect k and k ′ . That is, D kk ′ = λ k δ k,k ′ . Notice that the matrix D can be expressed also as: with the spectral density given by: 2. The basis function u k (t) is closed under differentiation: Notice that Bessel functions satisfy this condition since the following relation holds, (with J −k (x) = (−1) k J k (x) for an integer value of k).
See Equation (31) in Tokieda and Hagino [33] for the explicit form of the coupled-channels equations, Equation (20). Tokieda and Hagino [33] also provides an alternative derivation of the coupled-channels equations, which used the influence functional of the path integral method. This allows one to extend the present formalism to finite temperatures. Before the present method is applied to heavy-ion collisions, we should make sure that it works in principle. To this end, we considered a damped harmonic oscillator in which the exact solution can be obtained easily [33]. The Hamiltonian for the system, H S in Equation (10), is now given by: where M and ω S are the mass and the frequency of the system, respectively, and the last term represents the so-called counter term. In the following, we measured the length of the system in units of the oscillator length q S defined by q S ≡ h /Mω S , and took the coupling form factor, h(q), to be h(q) = q/q S . We assumed that the bath oscillators are distributed according to the spectral density (see Equation 23) as: In the numerical calculations shown below, we tookhω S = 2 eV, V I = 1 eV andh = 4 eV. At t = 0, we assumed thatψ {ñ k } (q, t = 0) = 0 for N ≡ K k=1ñk = 0. For N = 0, that is, forñ k = 0 for all k, we assumed that the wave function is given by: with q 0 /q S = −1, σ 0 /q S = 1/ √ 2, and p 0 q S /h = 0. Figures 3, 4 show the time evolution. The upper panel of Figure 3 shows the norm for each phonon state N as a function of ω S t. Here, the norm is defined as: To draw this figure, we took Bessel functions, J k ( t), for u k (t) in Equation (16) with K = 20. A new basis was then constructed by diagonalizing the matrix D in Equation (21). With this basis, we solved the coupled-channels equations by including the phonon states with N ≤ 5. The expectation value of the norm was also shown in the lower panel. As is expected, the number of phonons in the bath gradually increases as a function of time. Notice that the contribution of the 5-phonon states is small in the whole time range shown in the figure. This justifies the truncation at N max = 5 for the present parameter set. One can also see that the contribution of each phonon reaches its equilibrium at around ω S t = 6. Figure 4 compares the results of the present method with the exact solution for the quantum damped harmonic oscillator. To this end, we evaluated the expectation values for the following four quantities: ξ q ≡ q /q S , ξ p ≡ p q S /h, ξ qq ≡ q − q 2 /q 2 S , and ξ pp ≡ p − p 2 q 2 S /h 2 . We carried out the calculations with three different values of N max , that is, N max = 3, 4, and 5, and compare them with the exact results shown by the solid lines. One can see that all of the calculations with N max = 3, 4, and 5 reproduce the exact results up to ω S t ∼ 5, for which J 20 ( t) is negligibly small, and the expansion in Equation (16) with Bessel functions up to K = 20 [that is, up to J 19 ( t)] is therefore reasonable. The deviation from the exact results becomes significant for larger values of ω S t, especially for the second order moments, ξ qq and ξ pp . This is a natural consequence of the fact that the larger number of phonon states are required to describe the finer structures.
We would like to make a few comments on an application of the present method to heavy-ion collisions. Unlike the phenomenological quantum friction models discussed in the previous section, the total energy conservation was assured with the Hamiltonian given by Equation (10). Therefore, one can utilize the energy projection method [23] to calculate the penetrability for a given energy. In the present method, FIGURE 4 | Comparison between the present method and the exact results for the expectation values of several quantities, that is, (A) ξ q = q /q S , (B) ξ p = p q S /h, (C) ξ qq = (q− q ) 2 /q 2 S , and (D) ξ pp = (p− p ) 2 q 2 S /h 2 . The solid lines show the exact results, while the solid lines with squares, triangles, and circles are the results of the present method with N max = 3, 4, and 5, respectively. the number of N max required to achieve convergence was an important parameter that controlled the numerical cost. We expected that N max was small for collision problems since the system-bath coupling was active only while a wave packet overlapped with the potential region. Another issue in an application to heavy-ion collisions was that one needs to model properly the environmental degrees of freedom. Note that they are described solely by the spectral density (see Equation 23) in the Caldeira-Leggett model. That is, one needs to model a suitable spectral density for heavy-ion collisions. A similar problem has already been discussed in the linear response approach to heavy-ion collisions [34] as well as to fission reactions [35]. We anticipate that these approaches provide a useful means for our future works of an application of the present method to heavy-ion reactions.

SUMMARY
We have discussed two time-dependent methods for quantum friction. The first method was based on an effective Hamiltonian, which was constructed so that expectation values of operators obey a classical equation of motion with friction. Such Hamiltonian is in general time-dependent, and we have solved it with a time-dependent wave packet method. The other method is to start with a total Hamiltonian with both the system and the environmental degrees of freedom and then eliminate the environmental degree of freedom to derive an equation that the system degree of freedom obeys. For this approach, we have presented a new efficient basis for coupled-channels equations. These two methods were complimentary to each other. In the first method, whereas several parameters had to be determined phenomenologically, a required computational time was much shorter than the second method. On the other hand, the second method was based on a more microscopic Hamiltonian, and fewer empirical inputs were thus required even though the computational time may have been large. By combining these two approaches appropriately, one may be able to achieve a quantum description of heavy-ion deep inelastic collisions as well as fusion reactions to synthesize superheavy elements.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.