Simulation of Matrix Product States For Dissipation and Thermalization Dynamics of Open Quantum Systems

We transform the system/reservoir coupling model into a one-dimensional semi-infinite discrete chain through unitary transformation to simulate the open quantum system numerically with the help of time evolving block decimation (TEBD) algorithm. We apply the method to study the dynamics of dissipative systems. We also generate the thermal state of a multimode bath using minimally entangled typical thermal state (METTS) algorithm, and investigate the impact of the thermal bath on an empty system. For both cases, we give an extensive analysis of the impact of the modeling and simulation parameters, and compare the numerics with the analytics.

Introduction coupling between the system and the environment does not appear as the most general situation; for instance, nanomechanical resonators [22,23], circuit QED setups [16,17], optomechanical systems [24,25] and the impurity affected solid state systems [26]. In these cases, the theory has been implemented after transforming the nonlinear Hamiltonian to a linear one by linearizing the quantum fluctuation over a nonlinear steady state field amplitude. Appreciating its simplicity, this model, however, does not provide a satisfactory platform to obtain exact dynamical behavior. The interesting effects, e.g. non-classical behavior of the nonlinear systems [27], are often overlooked when we cannot handle the interaction between two systems in a perturbative manner. Apart from nonlinear S/B coupling, the analytical model is also limited to providing the exact solution in case of non-Markovian dynamical phenomenon, for instance phase transition in a two level system (TLS) between dynamically localized and delocalized states, at zero temperature for Ohmic and sub-Ohmic reservoirs [3,[28][29][30]. The limitation of the analytics explicitly provokes us to simulate the time evolution numerically. The numerical approach consists of transforming the environmental degrees of freedom to many body systems and simulate it in order to obtain the time evolution. The computational method includes a numerical diagonalization and renormalization process.
The time-adaptive density matrix renormalization group (t-DMRG) is considered as one of the most powerful methods in atomic, optical and condensed matter physics to study strongly correlated many-body quantum system. The method have often been used for some renowned models of quantum mechanics, such as the Ising model [31,32], the Hubbard model [33][34][35] and the Bose-Hubbard model [36][37][38], especially aiming to study the magnetization, quenching dynamics and phase transition properties. In order to implement it in this case, we map the canonical S/B model to a one-dimensional harmonic chain with nearest neighbor interactions.
Here, we present a numerical model for the analysis of the simple coupling between the system and environment, along with the consequences associated with the modeling and the numerical simulation. We start with simulating the dissipative dynamics of an open quantum system, an then we study the thermalization of the system in the presence of a thermal bath.
According to quantum statistical mechanics, the thermal state is a mixed state, and therefore, it is represented by an ensemble of pure states. It is naturally expressed by using density matrix for the inverse temperature β and Hamiltonian H. A few numerical approaches have been employed to study the impact of thermal bath on an interactive system, e. g. quasi adiabatic propagator path integral algorithm (QUAPI) [39,40] or solving hierarchical equations of motion (HEOM) [41], but in all cases, the influence of the bath on the dynamics of the system is taken care analytically using well defined Feynman-Vernon influence functional. The influence functional appears to be different for different types of coupling between system and environment, and in some cases, it becomes extremely difficult to determine, especially when nonlinear coupling appears. However, the method we discuss here has the ability to overcome this problem. This includes generation of thermal bath numerically and evaluation of time dynamics of both the system and environment. Even though the DMRG technique is designed to determine the ground state [42] and the time evolution of many body systems, a different method has been used to study the thermal state. Here, we introduce a complementary approach which includes taking a large number of sample pure states and determine an observable by averaging over them, instead of expressing the state by a density matrix. The states whose ensemble collectively generates the impact of thermal state, are determined through imaginary time evolution and projective measurements, typically known as minimally entangled typical thermal states (METTS) [43]. In this article, we impose the algorithm for the first time to generate the thermal bath, parameterize it and investigate the consequences to study the thermalization of open quantum system.

Theoretical model
In this article, we discuss the dynamics of a simple coupling model between the system and the reservoir. We start with the Hamiltonian of the system coupled to a thermal environment, which is given by is the Hamiltonian of the isolated system, ω c is the frequency, and ( ) † c c are the bosonic creation (annihilation) operators corresponding to the mode of the system. † , is chosen to be symmetric around the system mode (ω c ). The idea behind kind of modeling is the fact that, realistically, the system couples to the few modes of the environment around its resonating mode (ω c ). The model also includes consideration of the linear dispersion relation of the modes of the bath (w µ k k ).
In order to characterize the properties of the bath, we define the spectral density function J(ω k ) [30], implying the hard cutoff limit of the reservoir and wide band limit approximation ( where θ is the Heaviside step function, and g p = Dc 2 0 2 is the decay rate of the system where = = is the input field to the system. If the reservoir remains at zero temperature, no input field has contributed to the system. Therefore, the Heisenberg equation of motion (HEM) of the system mode is ( ) , which determines the free dissipative nature of the system population ( ( ) is the normalization coefficient.
The transformed Hamiltonian of the 1-D lattice chain is . The schematic diagram of the transformation is shown in figure 1(a). Recently, similar mapping was introduced in [45] to simulate open quantum systems aiming to be applied to spin-boson models [46] and biomolecular systems [47]. In all cases, the model had remain successful to overcome the complexity of the deduction of non-Markovian dynamical phenomenon, but the bath was considered to be at zero temperature. But, in the following section, we introduce METTS algorithm for the first time, for the generation of thermal bath at finite temperature and the evolution of system in the presence of the thermal bath.

Real-time evolution
We use time-evolving block decimation algorithm (TEBD) to do the numerical simulation, which requires expressing the state of the full chain as a matrix product state (MPS): a a a a a a a a a = = + - The MPS state is obtained through the Schmidt decomposition of the pure state of N sites where M is the dimension of local Hilbert space and χ is the Schmidt number [48]. The method of numerical simulation for the real-time evolution is shown diagrammatically in figure 1(b), where we choose the second order Suzuki Trotter (ST) expansion [49], which minimizes the error in third order of the time step by evolving the pairs of alternate sites. Using ST expansion, we express the unitary evolution operator as The simulation parameters are estimated by looking at errors which can appear in two ways: while modeling the S/B formalism to a one-dimensional chain and simulating each step during the real-time evolution. The errors are discussed extensively afterwards to estimate the parameters for numerical simulation.

Algorithm for thermal state
We imply the METTS algorithm by sampling over a huge number of pure quantum states [43]. Overall, these samples contain physical properties of the system for a given temperature, which approximates thermal where Z β is the partition function. The CPS | ñ n becomes a matrix product state (MPS) |f ñ n after the imaginary time evolution with the probabilities P n as In the next step, the METTS |f ñ n collapses to a new CPS | ¢ñ n through a projective measurement with an arbitrary measurement basis from which one can subsequently compute a new METTS |f ñ ¢ n , and, this process keeps on going on to generate a large set of MPS which typically represents a thermal state altogether. Thus, the generation of METTS samples undergoes a Markov process which is illustrated in figure 1(c). In this framework, the thermal average is determined from the set of imaginary time evolved normalized MPS states (|f ñ n ) with the probabilities P n /Z βˆ|ˆ| The algorithm has been used widely to simulate the spin chain at finite temperatures [50][51][52]. However, we use the technique to study the thermalization of the open quantum system numerically for the first time, which includes generation of the thermal bath using the METTS algorithm, and afterwards, evolve an empty system in the presence of the thermal bath.
The computation cost of TEBD increases with the entanglement of the quantum state, and hence the CPS is the natural choice to start with for having least entanglement. The entanglement of the obtained MPS states remains relatively low during real time evolution, which makes the simulation efficient.

Free dissipative system
We check the applicability of the TEBD algorithm in the dynamics of open quantum systems by comparing it with the analytics of a simple system/bath coupling model, where we assume that one photon is kept initially in the system and the bath is completely empty at zero temperature. Transforming the modes of the bath to a chain, we see that the first site is populated by a single quantum and all other sites remain empty. In figure 2, the dissipative nature of the system, which is obtained numerically by doing real-time evolution of the full chain, is compared with the analytics determined from the Heisenberg equation of motion (HEM) given by equation (4). We see an increment in the system population obtained in the numerical result after a certain time due to the reflection of the particle from the end of the chain, which is visible explicitly from the plots of the population of the full chain given in the insets of figure 2.

Recurrence time and density of states
The recurrence time decreases with the increment of the group velocity, causing the phonon to travel faster in the lattice. The group velocity is defined by where ω k is the frequency, and k N is the wavenumber determined by the number of lattice sites (k N ∝N). Eventually, the group velocity is inversely proportional to the density of states ( ) µv D g 1 , and therefore, the recurrence time increases linearly with the increment of DOS, which is seen in figure 3(a), where we increased the DOS by increasing the number of sites, keeping the cutoff frequency fixed.

Errors and estimation of parameters
The case of real-time evolution, the ST expansion introduces time error, which tends to concentrate in the overall phase [49]. In figure 3(b), we show how the accuracy of the simulation improved while reducing the time step. Truncation of the Hilbert space is also an issue in TEBD simulation. However, we chose the size of local Hilbert space to be 2 for having a single particle, and therefore the set is complete and we do not expect any error associated with the truncation of local Hilbert space. Because the time evolution started from an initial product state, the entanglement of the sites remains relatively low, and therefore the impact of the truncation in the Schmidt spectrum is also negligible.

Thermalization of an open quantum system
Analytics of the thermalization of a system Accounting for the back action of the environment, we determine the time dynamics of the field operator of the system by integrating QLE of a simple S/B coupling model, given in equation (3) as The initial thermal population distribution of the bath is As the coefficient B goes to zero at the steady state, the population of the system is determined by the function I k given in equation (14). Realistically, the population dynamics of the system are dependent on a few modes of the bath around the resonating mode of the system, and therefore, the cutoff limit (ò) should be considered in such a way that the contribution of the modes of the bath, those which far away from the system mode, can be marginalized. Hence, the function I k is expected to converge when the frequency (ω k ) goes far away from ω c . However, figure 4 shows that even though the function I k exhibits a peak at the resonating frequency of the system, it rises up again when the frequency of the bath mode is much lower than the frequency of the system mode ( . The acceptance up to the limit of the second order perturbation for the theoretical formulation of an open quantum system, is essentially based on the weak coupling between system and environment, ensuring the Lorentzian function acts like a delta function around ω c . It is also seen from figure 4 that I k rises up faster towards the lower cutoff limit in the case of a bath which is at a lower temperature than it does in the case of a high temperature. Such situations can even be bypassed by reducing ò, but that increases the ratio between γ and ò. However, we can solve this issue by reducing the value of γ, but that demands more time for the system to reach the stationary state, and hence it might not be possible at times to reach the steady state before the recurrence of the particle from the boundary. In that case, we increase the recurrence time by increasing the density of states. However, in the case of a zero temperature bath, the thermal population remains zero for all modes. Therefore, the relaxation of the system to the ground state is not affected by the cutoff frequency.
Anticipating the fact that the Lorentzian function becomes a delta function for a given condition g < < , the steady state population of the system is approximated to c which is the thermal population of the bath corresponding to the mode of the system. So the steady state population of the system is approximated to the population of the bath corresponding to the mode of the system.

Generation of the thermal bath
The quality of the thermal state generated by the METTS algorithm is dependent on two crucial parameters: temperature and number of samples. The frequency spectrum of the thermal population of the bath is plotted in figure 5(a), which determines that in the case of lower temperatures, as the thermal population reduces rapidly, the fewer modes are required to be taken into account to express a thermal state. This is also suggested by figure 5(b), which shows that the cumulative probability saturates faster for the low temperature bath, reducing the requirement of number of METTS samples to represent the thermal state. The consequence is observed in figure 5(c), where the plot of population distribution becomes smooth, and therefore defines a significant pattern while reducing the temperature for a fixed number of METTS, which indicates a better quality of the preparation of thermal state. In table 1 we compare the thermal population obtained analytically and numerically by taking average over 50 METTS samples. However, anticipating better performance of the METTS algorithm at low temperatures, it is also seen that the overall thermal population reduces so significantly that after a certain range, the number is not reliable for numerical simulation. Therefore, we prefer to generate a thermal state higher temperature in order to obtain the thermal population to a significant level, which forces us to take a large number of METTS samples into account while doing real time evolution. In figure 5(d), we show how the increment of the number of METTS samples modifies the population distribution over the entire lattice chain. As anticipated, the quality of the preparation of the bath improves while increasing the number of METTS samples, which is indicated by the improvement of the smoothness of the plot. As we increase the density of states for a fixed frequency range, the number of modes and number of lattice sites also increase, which essentially demands more METTS samples to represent a thermal state. Therefore, we see a poor population distribution in figure 5(e) compared to figure 5(d)[iv], when we doubled the DOS and keep the number of METTS fixed. However, the increase of the density of states increases the total population of the bath, which is shown in table 2.
Real-time propagation of systems coupled to thermal bath Hereafter, we study the thermalization of an empty system in the presence of a thermal bath at inverse temperature β=5[1/ω c ]. The time evolution of the system population for different cutoff frequencies and rates of dissipation are shown in figure 6. As anticipated from equation (13), the oscillation in the population of the system is introduced by the left tail of figure 4. The extension of the lower cutoff frequency contributes more oscillation to the dynamics, and more population in the stationary state of the system, which is visible when we compare figures 6(a) and (b). The higher value of γ also contributes more oscillation as an error to the dynamics  of the system population. We see the steady state population of the system is comparable to the thermal population at ω c , which is indicated by equation (16). In both figures 6(a) and (b), the system has not been able to achieve the steady state for the slow dissipation rate (especially γ=0.011 3ω c ). Therefore, we extend the recurrence time by increasing the DOS in figure 6(c), which gives sufficient freedom to the system to relax to the steady state. The numerical technique, therefore proves a promising scheme to study the open quantum dynamics. In order to investigate its applicability in the physics of quantum Brownian motion, we plot real-time dynamics of the quadrature fluctuations in figure 7 with a comparison to its analytics. The arbitrary quadrature is defined as , the quadrature fluctuation becomes phase (θ) independent ( ( ) ( ) ) d = + q X t n t 1 2 , and its time dynamics gives a pattern similar to the population dynamics.

Conclusion
In this article, we intended to investigate the applicability of METTS algorithm in the thermalization dynamics of open-quantum systems, anticipating the fact that the DMRG technique has the ability to extract out exact dynamics without linearizing nonlinear Hamiltonians. The consequences of this approach are demonstrated in terms of the efficiency of the algorithm with a discussion of advantages and disadvantages of this simulation. In this spirit, we also compare the numerical result with analytical result determined using Heisenberg equation of motion. We started with presenting a model that transforms the Hamiltonian of a quantum system coupled linearly to a discrete set of modes of a bosonic reservoir, to a Hamiltonian of a one-dimensional chain with nearest-neighbour interactions. We then used the model to study free dissipation and thermalization of that open quantum system. We found the recurrence time of the real-time evolution increases linearly with the increment of density of states. Our results also show that even though the minimally entangled typical thermal states (METTS) algorithm performs better at lower temperatures, we preferred to work at higher temperature in order to obtain the thermal population at a significant level and avoid unwanted error in the population dynamics of the system contributed by the lower cutoff frequency limit. Therefore, more METTS samples are taken into account, which consume more computation resources. In conclusion, one can say that the numerically generated thermal bath shows promise, but, this requires a compromise between the quality of the result and the computation resources. The numerical scheme presented here was mainly motivated by an attempt to determine the exact solution in the case of nonlinear coupling between the system and the environment [26], non-classical dynamics of non-linear systems [27], and reach out single photon limit in optomechanical systems [24,25]. The combination of real and imaginary time evolution of open quantum system will allow us to investigate quantum Brownian motion of topological quantum matter [53,54]. In addition, the method will be useful to study the non-Markovian dynamics and critical behaviors of the subohmic or ohmic spin-Bosson coupling models [3,4,28,29].