Excitation Dynamics in Chain-Mapped Environments

The chain mapping of structured environments is a most powerful tool for the simulation of open quantum system dynamics. Once the environmental bosonic or fermionic degrees of freedom are unitarily rearranged into a one dimensional structure, the full power of Density Matrix Renormalization Group (DMRG) can be exploited. Beside resulting in efficient and numerically exact simulations of open quantum systems dynamics, chain mapping provides an unique perspective on the environment: the interaction between the system and the environment creates perturbations that travel along the one dimensional environment at a finite speed, thus providing a natural notion of light-, or causal-, cone. In this work we investigate the transport of excitations in a chain-mapped bosonic environment. In particular, we explore the relation between the environmental spectral density shape, parameters and temperature, and the dynamics of excitations along the corresponding linear chains of quantum harmonic oscillators. Our analysis unveils fundamental features of the environment evolution, such as localization, percolation and the onset of stationary currents.


Introduction
The thorough understanding of transport of energy, heat, particle, or mass in complex quantum systems is of utmost importance both from a fundamental and technological point of view. Such a relevance is witnessed by the enormous efforts invested by the scientific community over the last decades on the theoretical and experimental investigation of the unique features of transport at the quantum regime.
Open quantum systems (OQS) formalism [18,19] has been widely employed for the description of quantum transport in the, often unavoidable, presence of additional and uncontrollable degrees of freedom interacting with the system under study. The tools provided by open quantum system theory led to the derivation of fundamental results allowing to understand and control, or at least mitigate, environmental effects. Such control, for example, is of utmost importance to preserve the quantum resources, as entanglement and coherence, that could enable the development of quantum devices possibly outperforming their classical counterparts. On the other side, the analysis of certain open quantum systems has unveiled the delicate interplay between coherence and sources of decoherence, as in the paradigmatic case of energy transport in disordered lattices [2,3,16,20].
The simulation of open quantum systems, on the other hand, represents a formidable task. Even when a microscopic description of the environment surrounding a quantum system is available,

Tedopa
Here and in what follows we consider a system interacting with a bosonic environment. The complete Hamiltonian reads (h = 1): where H S is the (arbitrary) free system Hamiltonian, H E describes the free evolution of the bosonic environmental degrees of freedom, and H I is the bilinear system -environment interaction Hamiltonian [50], and A S , O ω are self-adjoint operators. This last assumption is necessary for the Thermalized-TEDOPA (T-TEDOPA) mapping [51], that we will introduce later in this section. We assume that h(ω) has finite support [ω min , ω max ], with ω min < ω max , and define the spectral density (SD), namely the positive valued function J : [ω min , ω max ] → R + , as J(ω) = h 2 (ω). (2) As shown in References [41,42,46] the Hamiltonian (1) can be unitarily mapped into an equivalent one through by defining a countably infinite set of new operators where U n (ω) = h(ω)p n (ω).
The operators b n and b † n satisfy the bosonic commutation relations [b n , b † m ] = δ nm ; moreover, the polynomials p n (ω) are orthogonal with respect to the measure dµ = J(ω)dω and satisfy three-term recursion relations [42,46]. Thanks to these properties, the Hamiltonian (1) is mapped into the one dimensional Hamiltonian where, for the sake of definiteness, we have specialized the operator O ω in (1) to X ω = a ω + a † ω . After the mapping, the system interacts with the first site of a linear (infinite) chain of bosonic modes; the system-chain interaction strength is given by [42,46] whereas the frequency of the first TEDOPA chain oscillator is namely the first moment of the normalized measure J(ω)/κ 2 0 dω on [ω min , ω max ]. The remaining coefficients ω n and κ n are defined by the above mentioned three-terms recursion relations; while in certain cases it is possible to analytically determine their value [42], a numerically stable procedure is in general used [52,53].
For the following analysis, it is important to stress that the chain Hamiltonian H C E is made up of exchange terms b n−1 b † n + H.c. and therefore conserves the "number" operator, that is, Excitations can therefore be added or subtracted from the chain because of the interaction with the system. The initial joint system-environment state is assumed factorized ρ SE (0) = ρ S (0) ⊗ ρ E,β (0), with ρ E,β (0) a thermal state at inverse temperature β = 1/k B T, namely with Z ω = Tr[exp(−βωa † ω a ω )] the partition function. The initial state after the chain mapping is a factorized state ρ C If the environment is initially at zero temperature, its initial state is the vacuum state, and the initial state of the chain is also a factorized vacuum state |0 C E (i.e., b k |0 C E = 0, k = 1, 2, . . .): the chain contains therefore no excitations. This case provides us with the simplest setting where to analyze the transport properties of the chain corresponding to some representative spectral densities. As recently shown in Reference [51], however, by the spectral density transformation with J (ω) = sign(ω)J(|ω|), it is always possible to replace the thermal state of the original environment with the vacuum state of an extended environment, comprising negative frequencies.
As the spectral density (13) is now temperature dependent, the TEDOPA chain coefficients ω n,β , κ n,β will be temperature dependent as well. In the following we will drop the β dependence wherever clear from the context. From now on will therefore always consider the factorized vacuum state as the initial chain state without loss of generality. In our analysis we will consider, in particular, the Lorentzian spectral density and Ohmic spectral densities defining a very important class of environments entering in the study of many systems, such as microscopic models leading to a Lindblad master equation for an harmonic oscillator in a weakly coupled high temperature environment, or a particle undergoing quantum Brownian motion [18,54,55]. From now on frequencies will be in cm −1 and temperatures in Kelvin. We remark that, because of the relation (8), if two spectral densities differ only for the overall coupling constant λ, their mappings (i.e., all of the chain coefficients ω n , κ n ) will be identical, with the exception of κ 0 , namely the coupling strength between the system ant the first TEDOPA chain mode. We also observe that the chain coefficients ω n , k n , n ≥ 1 are independent of the specific system-environment interaction term, that is, they are independent of the choice of A S , O ω of Equation (1). As customary for chain mappings, in what follows, we will moreover impose a hard cutoff ω hc to the considered spectral densities, thus limiting their support to the interval [0, ω hc ] for T = 0 and to the interval [−ω hc , ω hc ] for T > 0. The value of ω hc is suitably chosen as to keep the neglected relative reorganization energy ∞ ω hc dω J(ω)/ω ∞ in the order of 10 −4 for all the considered instances. If the considered spectral density belongs to the Szegö class, the asymptotic relations hold (see Theorem 47 of Woods et al. [46]). Clearly enough, in our setting ω max (and, at finite temperature, ω min ) depends on the imposed hard-cutoff ω hc so that both ω ∞ and κ ∞ are functions of ω hc . For any suitably fixed ω hc , however, the relations (16) allow for the simple heuristic estimation L = 2κ ∞ t max of the maximal distance travelled within the time t max by an excitation initially located at the first TEDOPA chain site. For fixed time t max , therefore, the effective environment is made up of L oscillators within the "light-cone". Interestingly enough, the width of such light cone depends only on the "artificially" imposed hard cutoff and, as long as the choice ω hc is sensibly chosen, different choices of the hard-cutoffs do not impact on the reduced dynamics of the system. On the other side, different spectral densities with the same support will have the same asymptotic coefficients, and the differences in the reduced dynamics of the system will be due to a (typically quite small) finite number of modes, as we will see in the following sections.

Chain Dynamics
We start by analyzing the dynamics of a single excitation moving along the chain-mapped environment produced by the (T-)TEDOPA mapping. To this end, we can disregard the system and the interaction term H I , or equivalently set κ 0 = 0, and restrict our attention to the single excitation sector of the TEDOPA-chain Hilbert space. The set {|k , k = 1, 2, ...}, where |k indicates the Fock state |n 1 = 0, . . . , n k−1 = 0, n k = 1, n k+1 = 0, . . . with the single excitation located at the k-th chain site, is a basis for the considered single excitation subspace. In what follows we will assume that the excitation is initially located at site 1, namely the initial state is |1 .

Lorentzian Spectrum
The Lorentzian spectral density (14) provides a paradigmatic example. For γ/Ω 1, such spectrum well approximates that of an environment made up of a single harmonic oscillator with frequency Ω and dissipating into the vacuum at rate γ [18,56]. In all the following examples a hard cutoff frequency ω hc = 10Ω has been enforced. Frame (a) of Figure 1 shows the frequencies ω n and couplings coefficients κ n at T = 0 for Ω = 100 and γ = 0.001, 1, 10 (see (14)). We first observe that, for all values of γ, the first and the second TEDOPA chain modes are equally far detuned. The main difference between the three selected cases lies in the coupling strength κ 1 between the same two modes, which is directly proportional to γ. The effect on the system dynamics is remarkable. As shown in Figure 1b the population of the first TEDOPA chain is well approximated by p 1 (t) = exp(−2γt), namely the decay rate of an harmonic oscillator damped into the vacuum at a rate γ. As frame (c) of the same figure shows, the portion of excitation that propagates beyond the first site propagates on the TEDOPA chain at a speed which is independent of γ: the chain coefficients are essentially equal to each other in the three cases for n ≥ 3, and their value is determined by the hard cutoff frequency ω hc through (16).
We turn now our attention to the finite temperature case. As exemplified in Figure 2a, after the thermalization procedure [51] the thermalized spectral density (13) presents two peaks at ±Ω. The system will be thus effectively coupled to two damped modes, with temperature dependent coupling strength proportional to 1 + n β (Ω) resp. n β (Ω).    (13) and (14) It is thus not surprising that the chain dynamics for the case γ = 0.001 is essentially confined to the first two chain modes, as frames (b) and (c) of Figure 2 show. Indeed, the same plots suggest a clear relation between the temperature and the relative occupation of the modes: as T increase, the difference between the maxima of the populations of the first and the second TEDOPA chain sites decreases, and is expected to vanish as T → ∞, that is, when the thermalized spectral density becomes symmetric with respect to the origin.
It is interesting to see that a mechanism very similar to the one discussed for the zero temperature case is at play also at finite temperature. Frame (d) of Figure 2 shows the chain coefficients for γ = 0.001, 1 and 10 at T = 300. This time the detuning between the first and the second TEDOPA chain sites is relatively small and the coupling between the two sites is independent of γ. This time it is the detuning between the second and the third chain site that is considerable, and the coupling κ 2 is monotone with γ. As shown in frame (e) of the same figure, the result is that the population of the first TEDOPA chain site presents damped beatings: the excitation moves forth and back between the first two chain sites, and percolates toward the right part of the chain at a rate exp(−2γt). In the zero temperature case, the "escaped" population travels toward the right part of the chain at a speed which is independent of γ, and keeps trace of such beatings, as shown in Figure 2f, but this time the propagation speed is twice that of the zero temperature case because of the enlarged support [−ω hc , ω hc ] (see (16)).
In order to provide an insight on how the chain dynamics depends on the temperature, in Figure 3 we show the population of the first TEDOPA chain site for different values of T. As already observed, the decay rate and the frequency of the population oscillations are independent of T, which determines instead the amplitude of such oscillations. This leads us to the conclusion that the oscillation frequency must be determined by the parameter Ω, as expected.

Ohmic Spectrum
We now consider spectral densities belonging to the Ohmic family, defined as in Equation 15. More in particular, we will study the chain dynamics on TEDOPA chains corresponding to the choice s = 0.5, 1 and 2, representative, respectively, of sub-Ohmic, Ohmic, and super-Ohmic spectral densities. In all the following examples we will set ω c = 100, and enforce a hard cutoff ω hc = 10ω c .
We start by the T = 0 case. Frames (a) and (b) of Figure 4 show Ohmic spectral densities for the selected values of s and the corresponding chain coefficients. The chain dynamics shows that an excitation leaves its initial location faster in the super-Ohmic case than in the Ohmic and sub-Ohmic case (see Figure 4c). This can be justified by the higher coupling coefficient and smaller detuning between the first sites of the TEDOPA chain in the s = 2 case with respect to the s = 0.5, 1 cases. Moreover, even if the front of the excitation wavepacket travels at the same speed in the three cases, the delocalization degree of the wavepacket is higher in the sub-Ohmic case, while it remains more "compact" in the super-Ohmic case, as examplified by the inset of Figure 4c, showing the TEDOPA chain site populations p x (t) at t = 0.1. Considered that the chain coefficients for s = 0.5, 1, 2 are very close to each other for n ≥ 4, this difference is explained by the first chain coefficients. Roughly speaking, the higher coupling and smaller detuning between the first chain sites in the s = 2 case allows for more compact evolution of the wavepacket in the momentum space.
In the high-temperature regime T = 300, the main features of the chain dynamics are preserved, though with some differences. The decrease of population the first TEDOPA chain oscillator is still slower in the sub-Ohmic case; for the Ohmic SD, the first site population decay is similar to the T = 0 case, whereas for the super-Ohmic SD such decay is faster than in the zero temperature scenario (compare frames (c) and (f) of Figure 4). As already discussed before, this behaviour is mainly due to the detuning |ω 1 − ω 2 | and the coupling strength κ 1 between the first and second TEDOPA chain oscillators. Interestingly enough, for s = 0.5 part of the wavepacket remains localized at the first chain site, as shown in the inset of Figure 4f and, as in the T=0 case discussed above, the wavepacket is more delocalized in the sub-Ohmic case than in the super-Ohmic case, with the Ohmic case lying in between.
As a last remark, we observe that, similarly to the finite temperature Lorentzian case, the propagation speed of the wavepacket is about twice as large as in the zero temperature case; as already discussed, this is due to the asymptotic relations (16). 20 40 60 80 100    Figure 5 provides more details. As we did for the Lorentzian SD case, we now inspect the dynamics of the first TEDOPA chain population for the three considered spectral densities at different temperatures. It clearly shows that, while for the Ohmic spectral density such population is only slightly affected by the value of T, the temperature has opposite effects on super-and sub-Ohmic SDs. As a matter of fact, whereas for the sub-Ohmic case, an increasing temperature leads to a slower decrease of the first site population, in the for s = 2 the first site empties at a rate which is directly proportional to the temperature. The snapshots on the populations p x (t) for t = 0.02 in the insets of frames (a)-(c) of the same figure, allows us to better appreciate the partial trapping at finite temperature of the wavepacket at the first TEDOPA chain site and the more pronounced spreading of the wavepacket in the s = 0.5 case.

Full Dynamics
So far we focused our analysis on the dynamics of a single excitation moving along TEDOPA chains. This allowed us to isolate the main features of such dynamics for representative spectral densities and to investigate the dependence of the kinematic properties of TEDOPA chains on the specific form of the SD and on the temperature. Clearly enough, the single excitation subspace we restricted ourselves to is not suited to describe the chain dynamics in the presence of a system interacting with the environment. As a matter of fact, the interaction with the system will dynamically inject in (and subtract from) the chain excitation, at a rate that depends, among other things, on the system-environment coupling strength.
In this section, therefore, we extend our analysis by considering a two-level system interacting with a bosonic environment described by either Lorentzian or Ohmic spectral densities. Given the spectral density, the spin-boson model is fully specified once the system and the system-environment interaction Hamiltonian are fixed. In what follows, we specialize the Hamiltonian (1) to with σ x , σ z Pauli matrices, describing, for example, an homo-dimer interacting with a vibronic environment [57]. The resulting dynamics is therefore not a pure dephasing dynamics, and is representative of the class of physical systems for which numerically exact approaches are required. Considered that the interaction term does not change the system's populations but affects only its coherences, we will initialize the system to the state |+ = 1/ √ 2(1, 1) T , namely the eigenstate of σ x belonging to the eigenvalue +1, representative of the maximally coherent states in the σ z basis. The initial state of the environment will be instead a thermal state (11) at temperature T. In the following examples we will set ∆ = 70cm −1 , and tune the parameter λ of Equations (14) and (15) so that the system-TEDOPA chain coupling κ 0 (see (8)) is the same at T = 0 for all the considered spectral densities. More precisely, by definition, the k 0 coefficient of the Ohmic spectral density is independent of s so that, in the Ohmic cases, we set λ = 1; for Lorentzian spectral densities we set to λ = 60. Before presenting our results it is important to remark that we are not so much interested in the reduced dynamics of the system, but rather on the TEDOPA chain dynamics in the presence of an interaction with the open system. In particular, we will try to understand which of the features discussed in the preceding section persist in the presence of an interaction with the system. To this end we will use the average occupation number of the k-th chain oscillator where ρ C (t) is the system+chain state at time t determined via TEDOPA simulation. We first discuss the chain dynamics for Lorentzian spectral densities. The γ = 0.001 case is still paradigmatic. At T = 0 only the first TEDOPA chain oscillator is essentially involved in the dynamics. By comparing the purple lines in frames (a) and (b) of Figure 6, we can clearly see the beatings between the system and the first TEDOPA chain site. For T > 0 a the second TEDOPA chain mode enters into play. The average occupation number n 1,2 (t) of the first two chain sites depend on the temperature. Interestingly enough, in the high (T = 300) temperature regime the both n 1 (t) and n 2 (t) present small and fast out of phase oscillations, imprinting on the system dynamics a much more erratic dynamics than the T = 77 environment, for which such oscillations are slower and almost in phase. Figure 7 shows instead the system and chain dynamics for γ = 10. Analogously to the γ = 0.001 the average occupation of the first two TEDOPA chain sites is temperature dependent. The larger value of γ implies that, loosely speaking, more environmental modes are interacting with the system. While the first two sites are still the highest occupied ones, some excitations can percolate to the right part of the chain, as we already observed in the chain dynamics analysis of the previous section (see Figure 2f). Since the system-TEDOPA chain coupling is about the same for the two considered values of γ, it is such percolation responsible for the faster relaxation of the system. Now we turn our attention to Ohmic spectral densities. As it happens for the Lorentian case discussed above, the main features of the excitations dynamics presented in Section 3 provides a key to understanding the results. We observed (see Figure 4c,f) that an excitation located at the first chain site will leave its initial location more slowly in the sub-Ohmic case than in the Ohmic and super-Ohmic case. Moreover, the excitation wavepacket tends for s = 0.5 to be more spread over the chain than for s = 2, with the case s = 1 showing an intermediate behaviour. This features translate to the chain dynamics in the presence of an interaction with the system, as comparison between Figures 8-10 shows.
In more detail, we observe that at T = 0 the excitations leave the first chain sites almost ballisticaly for s = 2 (Figure 10b), whereas for s = 0.5 there is an accumulation of excitations in the very first part of the chain (Figure 8b). The diagonal fringes appearing in the sub-Ohmic (and less pronounced in the Ohmic) case at zero temperature are easily explained in therms of the (moving in time) population profile shown in the inset of Figure 4c. The inclination of the fringes, is instead related the the coupling coefficients beteween the TEDOPA chain oscillators that, as already pointed out, do not depend on s but only on the spectral density support.
At finite T the situation changes quite drastically. First of all we observe that for all the chosen values of s vertical fringes appear in frames (c) of Figures 8-10. Such vertical fringes can be associated to a the alternation of higher and lower average occupation number in nearest-neighbor sites, and allow to appreciate the onset of a stationary current when the state of the system gets close to its stationary state. A comparison between frames (a) of the same figures shows that in the sub-Ohmic the average occupation of the first TEDOPA chain sites is much higher than in the Ohmic and super-Ohmic cases. It must be noticed that, while the system-TEDOPA chain coupling κ 0 is equal for T = 0 for all values of s, at finite temperature such coupling is inversely proportional to s. The sub-Ohmic TEDOPA chain is therefore more strongly coupled to the system, and this justifies the faster system dynamics at short times.

Conclusions and Outlook
While chain mapping has been recognized as a powerful tool for the efficient simulation of open quantum system dynamics, the subtle role of excitation dynamics on the determination of such reduced dynamics has never been investigated in detail. This work represents a first step in this direction. While the single excitation dynamics is unable to capture the full complexity of the evolution of TEDOPA chains put in interaction with the system, it provides a most useful key to understand such evolutions, as in the case of the Lorentzian spectral density we considered. It moreover provides a mean to sensibly set DMRG parameters, such as the chain truncation length and the local dimension of the chain oscillators: for super-Ohmic SDs, for example, the local dimension of the first TEDOPA chain oscillators must be set large enough as to host all the excitations that will accumulate in proximity of the system because of localization, while in the super-Ohmic case the local dimension of the first chain oscillators can be kept much smaller, since there is no signature of localization. While an analysis along the same lines for a specific spectral density was already presented [51], in this work we systematically compared and contrasted the features of the chain and full dynamics for a larger and very representative class of spectral densities. This allowed for example to shed light on the mechanisms allowing oscillators chain obtained by the T-TEDOPA procedure, and therefore starting from the vacuum state, to mimic an environment in the thermal state. For the Lorentzian case study such mechanism emerged quite clearly, and provided an key for the interpretation of the chain dynamics for SDs belonging to the Ohmic family.
We moreover observed that, while the asymptotic values of the TEDOPA coefficients determine the maximum distance reachable within a given time by an excitation initially located at the beginning of the TEDOPA chain, or light-cone, the features of a specific spectral density are typically determined by a very small number of coefficients. Indeed, as it happens in the γ = 0.001 Lorentzian SD case, the propagation of excitations in the light-cone can be hindered by an "effective" decoupling of the first sites of the chain from the remaining one. The analysis of the Ohmic SD instances, on the other side, showed that different (s-dependent) chain coefficients in the very fist part of the chain lead to quite different occupation probability profiles of the sites within the light-cone.
One of the, so far unexploited, advantages of chain mapping is the possibility of acquiring information on the state of the environment, something not meaningful when effective dynamics of Lindblad or Bloch-Redfield type are employed. While the number of chain modes perturbed by the interaction with the system is, in general, increasing with time, at any finite time it is in line of principle possible to make measurements on the oscillators in the light-cone. This could allow to understand, for example, which environmental modes are more involved in the dynamics and properly select the environmental reaction coordinates [58]. Moreover, in the presence of a fast convergence of the chain coefficients toward the asymptotic values, one expects a very small number of such coordinates. This represents a possible line of future research.
There are features of the TEDOPA chain evolution that remained quite obscure. For example, the fringes that appear in the Ohmic scenario at finite temperature are not present in the Lorentzian case. Considered that, as already observed, for γ/Ω 1 a Lorentzian environment can be assimilated to a damped harmonic oscillator undergoing a Lindblad-type dynamics, an therefore incoherently dissipating into an memoryless environment, one could read the lack of fringes as a signature of incoherent dynamics. A further analysis is therefore needed to better qualify the coherence dynamics in structured environments, and will be the focus of future work.
Funding: This research received no external funding.