Linear-algebraic bath transformation for simulating complex open quantum systems

Problems associated with open quantum systems are of interest in various research fields [1]. In the theory of open quantum systems, the universe is partitioned into system and bath components. The system of interest is then coupled to the bath degrees of freedom (DOF) by means of an effective Hamiltonian. Solving the full quantum dynamics with currently known exact analytic or numerical methods is not feasible as the system and bath DOF increase. A simple but effective approach to model these vibrations is to treat them as a collection of non-interacting quantum harmonic oscillators bilinearly coupled to the system [1, 2]. The system-bath interaction is then characterized by a spectral density function (SDF) that represents the coupling strength in the frequency domain [1, 2]. The energy transfer in photosynthetic systems is an example of a complex open quantum system, where the pigments involved in the energy transfer interact with a richly structured set of molecular vibrations, and hence a very structured SDF [3].

The spin-boson model [1] is one of the simplest models for studying the dynamics of open quantum systems. In the most common representation, the spin-boson model is mathematically represented as a set of non-interacting oscillators coupled to the system. This can be graphically represented in a star configuration as shown in figure 1(A) [4, 5]. Generalized spin-boson models, such as the Hubbard-Holstein model, have been successfully employed to describe the energy transfer process in photosynthetic antenna complexes [6]. Some numerically-exact methods have been so far developed to solve these models; see, for example, reduced-density-matrix approaches, such as hierarchy equations of motion (HEOM) [7–9], stochastic approaches [10–14], multi-configuration time-dependent Hartree [15], numerical renormalization group [4, 5, 16–19], and path-integral approaches [20], amongst many others. However, the applicability of numerically exact methods is limited by the system size and the bath DOF. For example, the simulation of HEOM with current computers is limited to ≈40 sites of the system with only a single Drude–Lorentzian peak representing the bath [9, 21]. On the other hand, the renormalization group approach [16, 19, 22] could be used for relatively large systems. However, the system and bath size that can be handled is still far from that required for solving problems at biological scales. In this approach, then, the bath transformation from the non-interacting bath model (figure 1(A)) to the 1-D Wilson chain (figure 1(B)) [23–25] is necessary. The collective modes of the bath oscillators in the chain model have been used widely in various fields such as quantum molecular dynamics [26, 27], open quantum dynamics [28–32], quantum information [33] and nuclear physics [34, 35].

Zoom In Zoom Out Reset image size

Quantum simulators defined as controlled quantum devices that can effectively reproduce the dynamics of any other quantum system [36–39] could become an attractive alternative for solving the dynamics of open quantum systems 'directly'. Different platforms can be used for implementing quantum simulators, such as superconducting qubits [21, 33, 40–46], trapped ions [47–55], quantum optics [56–60], nuclear magnetic resonance [12, 61–64] or a system of electrons [65, 66].

The experimental implementation of quantum simulators of open quantum system dynamics, for example, using superconducting circuits [21, 41], poses challenges due to at least two of the main current constraints in the realizable circuits. First, the number of quantum bath oscillators, which are directly coupled to a system operator (qubit), is limited by the physical size of the superconducting loop that embodies the qubit. Hence, a star-model approach with many oscillators coupled to the qubit may pose fabrication challenges. A physical layout that involves fewer oscillators directly coupled to a qubit is more experimentally realizable [41]. In addition, the coupling strength of the qubit to the bath may be limited. In superconducting qubits, the system-bath coupling strength should not exceed a certain percentage of the frequency of the quantum oscillator [21].

In this work, we address the question of how the two afore-mentioned implementation issues can be resolved by a unitary bath transformation which introduces interaction terms among the transformed quantum oscillators. In chain-like bath models (figure 1(B)) [4, 5, 18, 22, 31–33, 67, 68], only one bath oscillator is directly coupled to the system. However, in some cases, one needs to couple more than a single chain to deal with the limitation of the oscillator-qubit coupling mentioned above. Here, we propose a partitioning strategy of the bath modes for multiple parallel chains to reduce primary mode coupling strengths and also the number of the modes directly coupled to the system operator. This is shown in figures 1(C) and (D), respectively. We found that the coupling strength of the primary modes, which are directly coupled to the system, can be reduced as we increase the number of chains; at the same time, we can also shorten the lengths of each chain. In addition to the fabrication and implementation benefits for open quantum simulators using quantum hardware, these methods are also potentially applicable to simulations in classical computers. In this case, perturbative methods may be employed to simulate these chain models with reduced system-bath coupling [31, 32].

A recurrence equation derived by Bulla et al [5] has been used in the renormalization group approaches to construct the 1-D Wilson chain (figure 1(B)). This recurrence relation, however, potentially shows numerical instabilities [5, 17, 69]. Recently, Chin et al [22] developed an exact mapping approach for a continuous SDF using orthogonal polynomials without discretizing the SDF. However, this approach may pose challenges applicable to arbitrary structured SDFs. In many applications of chemistry and biology, structured SDFs appear when atomistic details are involved in open quantum dynamics, as already mentioned in the introduction. Therefore, the SDFs may not be well approximated as simple analytic functions such as Ohmic spectral functions.

In this paper, we test a generalized linear algebraic transformation approach for any given discrete SDFs, where a transformation on multiple parallel chains is involved, as shown in figure 1(C). With the multiple-chain–bath transformation described in the following sections, complex open quantum systems, such as photosynthetic antennae, can be studied practically via quantum simulators.

In the next sections we present the model Hamiltonian and the linear algebraic bath transformation. As an example, a two-oscillator bath transformation is presented analytically. This example shares many features of our general scheme of bath partitioning. Numerical stability of the bath transformation methods is discussed in the result section and the results are compared with Bullaʼs transformation approach [5]. Then, we propose a way to partition the bath modes into multiple parallel chains to reduce the system-bath coupling strengths. We apply the proposed leaping partitioning (LP) strategy to a structured spectral density of the chlorosome [70, 71], as an example. The numerical result is compared with a 'standard' sequential partitioning (SP) scheme.

As mentioned in the Introduction, in the theory of open quantum systems, the system-bath Hamiltonian $\hat{H}$ is composed of three parts, namely,

$$\begin{eqnarray}&&\hat{H}={{\hat{H}}_{{\rm S}}}+{{\hat{H}}_{{\rm SB}}}+{{\hat{H}}_{{\rm B}}},\end{eqnarray} \tag{ 1 }$$

where ${{\hat{H}}_{{\rm S}}}$ is the system Hamiltonian. The phonon bath ${{\hat{H}}_{{\rm B}}}$ is approximated as a set of non-interacting harmonic oscillators. The coupling term ${{\hat{H}}_{{\rm SB}}}$ between the system and bath is almost universally treated as a bilinear coupling. More precisely, we write ${{\hat{H}}_{{\rm SB}}}+{{\hat{H}}_{{\rm B}}}$ in a compact form [72, 73] as follows,

$$\begin{eqnarray}{{\hat{H}}_{{\rm SB}}}+{{\hat{H}}_{{\rm B}}}=\mathop{\sum }\limits_{n}\left( \begin{array}{ccccccccccccccc} \hat{L}_{n}^{\dagger } & {\boldsymbol{\hat{a}}} _{n}^{\dagger } \\ \end{array} \right){{{\boldsymbol{ \Gamma }} }_{n}}\left( \begin{array}{ccccccccccccccc} {{{\hat{L}}}_{n}} \\ {{{{\boldsymbol{\hat{a}}} }}_{n}} \\ \end{array} \right),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}{{{\boldsymbol{ \Gamma }} }_{n}}=\left( \begin{array}{ccccccccccccccc} 0 & {\boldsymbol{ \kappa }} _{n}^{\dagger } \\ {{{\boldsymbol{ \kappa }} }_{n}} & {{{\boldsymbol{ \Omega }} }_{n}} \\ \end{array} \right),\end{eqnarray} \tag{ 3 }$$

where each N-dimensional creation (annihilation) operator vector ${\boldsymbol{\hat{a}}} _{n}^{\dagger }\;\;({{{\boldsymbol{\hat{a}}} }_{n}})$ of oscillators for the site n are coupled to the operator ${{\hat{L}}_{n=n\gt \lt n}}$ that acts on the system. We note here that the bath transformation is independent of the system-bath coupling when the coupling is bilinear. Therefore the bath transformation can be applied to the spatially correlated noise as well, where the system coupling operator is ${{\hat{L}}_{nm}}=|n\rangle \langle m|$, $n\ne m$. Lower case bold and capital bold fonts are used for a column vector and a matrix, respectively. The bosonic operators satisfy a commutation relation, $[{{\hat{a}}_{n;i}},\hat{a}_{m;j}^{\dagger }]={{\delta }_{nm}}{{\delta }_{ij}}$ where ${{{\boldsymbol{\hat{a}}} }_{n}}={{({{\hat{a}}_{n;1}},\ldots ,{{\hat{a}}_{n;N}})}^{{\rm t}}}$. Here ${{{\boldsymbol{ \Omega }} }_{n}}$ is a diagonal matrix, which has the harmonic frequencies as the elements, i.e. ${{{\boldsymbol{ \Omega }} }_{n}}={\rm diag}({{\omega }_{n;1}},\ldots ,{{\omega }_{n;N}})$. ${{{\boldsymbol{ \kappa }} }_{n}}$ is the system-bath coupling strength vector. Accordingly, the SDF ${{J}_{n}}(\omega )$ is defined [74] as,

$$\begin{eqnarray}&&{{J}_{n}}(\omega )=\pi \mathop{\sum }\limits_{j}|{{\kappa }_{n;j}}{{|}^{2}}\;\delta (\omega -{{\omega }_{n;j}}).\end{eqnarray} \tag{ 4 }$$

In this work, the coupling strength vector ${{{\boldsymbol{ \kappa }} }_{n}}$, which is in general a complex vector, is chosen to be a positive real vector. That is, the vector elements are given by ${{\kappa }_{n;j}}=\sqrt{|{{\kappa }_{n;j}}{{|}^{2}}}\geqslant 0$. The non-interacting bath in equation (2) is the star-bath model (figure 1(A)), where the independent harmonic oscillators are all coupled directly to the system.

2.1. Linear algebraic bath transformation

With a suitable choice of unitary transformation on the bath oscillators, one can turn a star-bath into a multiple-chain-bath. The multiple-chain-bath has a few primary bath oscillators and the remaining oscillators (secondary bath modes) are coupled to the primary bath modes in a chain as depicted in figure 1(C). A mixture of star and chain models is also possible as shown in figure 1(D). Burghardt et al [31, 32, 75] exploited the latter model to develop a perturbative truncated bath model, which approximates the terminal star-coupled yellow oscillators in figure 1(D) as Markovian baths.

The bath transformation from the star model (equation (2)) to the 1-D Wilson chain (figure 1(B)) can be simply obtained by a unitary transformation of the matrix ${{{\boldsymbol{ \Gamma }} }_{n}}$ that keeps the system operators unchanged. We introduce, here, an arbitrary unitary transformation (${{{\bf U}}_{n}}{\bf U}_{n}^{\dagger }={\bf I}$) satisfying the following conditions:

$$\begin{eqnarray}{{\hat{H}}_{{\rm SB}}}+{{\hat{H}}_{{\rm B}}}=\mathop{\sum }\limits_{n}\left( \begin{array}{ccccccccccccccc} \hat{L}_{n}^{\dagger } & {\boldsymbol{\hat{b}}} _{n}^{\dagger } \\ \end{array} \right){{\tilde{{\boldsymbol{ \Gamma }} }}_{n}}\left( \begin{array}{ccccccccccccccc} {{{\hat{L}}}_{n}} \\ {{{{\boldsymbol{\hat{b}}} }}_{n}} \\ \end{array} \right),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}{{\tilde{{\boldsymbol{ \Gamma }} }}_{n}}=\left( \begin{array}{ccccccccccccccc} 0 & \tilde{{\boldsymbol{ \kappa }} }_{n}^{{\rm t}} \\ {{\tilde{{\boldsymbol{ \kappa }} }}_{n}} & {{\tilde{{\boldsymbol{ \Omega }} }}_{n}} \\ \end{array} \right)=\left( \begin{array}{ccccccccccccccc} 1 & {{0}^{{\rm t}}} \\ 0 & {\bf U}_{n}^{\dagger } \\ \end{array} \right){{{\boldsymbol{ \Gamma }} }_{n}}\left( \begin{array}{ccccccccccccccc} 1 & {{0}^{{\rm t}}} \\ 0 & {{{\bf U}}_{n}} \\ \end{array} \right),\end{eqnarray} \tag{ 6 }$$

where ${{{\boldsymbol{\hat{b}}} }_{n}}={\bf U}_{n}^{\dagger }{{{\boldsymbol{\hat{a}}} }_{n}}$. The first column of ${{{\bf U}}_{n}}$ is $||{{{\boldsymbol{ \kappa }} }_{n}}||_{2}^{-1}{{{\boldsymbol{ \kappa }} }_{n}}$ and the remaining columns are constructed using the Gram–Schmidt process with random vectors (or unit vectors) [76]. As a result, ${{\tilde{{\boldsymbol{ \Gamma }} }}_{n}}$ is a dense symmetric matrix and ${{\tilde{{\boldsymbol{ \kappa }} }}_{n}}={{({{\tilde{\kappa }}_{n;1}},0,\ldots ,0)}^{{\rm t}}}$ is the new system-bath coupling strength vector. (See appendix A for the details, and down panels of figure 1 for the structures of ${{\tilde{{\boldsymbol{ \Gamma }} }}_{n}}$.)

Now we have new sets of interacting harmonic oscillators while the system operators remain unchanged. The tridiagonalization of ${{{\boldsymbol{ \Gamma }} }_{n}}$ in equation (2) for the Wilson chain (figure 1(B)) can be performed numerically by Householder or Lanczos procedures [76]. Alternatively, we use here tridiagonalization of ${{\tilde{{\boldsymbol{ \Omega }} }}_{n}}$ with a Hessenberg reduction of a symmetric matrix [76]. We refer to the later transformation method as Gram–Schmidt–Hessenberg (GSH). The Hessenberg reduction of a symmetric matrix produces a tridiagonal matrix and then the numerical procedures for the reduction, such as Householder, Lanczos and Gauss transform, can be applied. These numerical algorithms are standard numerical linear algebraic techniques, see e.g. [76].

2.2. Multiple chain transformation

In this subsection we explain the multiple parallel chain transformation that is depicted in figure 1(C). We introduce a unitary transformation,

$$\begin{eqnarray}&&{{\tilde{{\bf U}}}_{n}}={{{\bf P}}_{n}}{{{\bf U}}_{n}},\end{eqnarray} \tag{ 7 }$$

that additionally rearranges (using a permutation matrix ${{{\bf P}}_{n}}$) the non-interacting bath oscillators as multiple groups of several interacting oscillators, i.e. ${{{\boldsymbol{\hat{b}}} }_{n}}=\tilde{{\bf U}}_{n}^{\dagger }{{{\boldsymbol{\hat{a}}} }_{n}}$. The unitary transformation matrix ${{{\bf U}}_{n}}$ is block diagonal and does not allow the interaction between oscillators from different groups (an example of the rearrangement is given in appendix B). We also define the following relations for the unitary transformation:

$$\begin{eqnarray}&&{{\tilde{{\boldsymbol{ \Omega }} }}_{n}}=\tilde{{\bf U}}_{n}^{\dagger }{{{\boldsymbol{ \Omega }} }_{n}}{{\tilde{{\bf U}}}_{n}}\quad {\rm and}\quad {{\tilde{{\boldsymbol{ \kappa }} }}_{n}}=\tilde{{\bf U}}_{n}^{\dagger }{{{\boldsymbol{ \kappa }} }_{n}}\ .\end{eqnarray} \tag{ 8 }$$

The primary modes, which are directly coupled to the system operators, are defined as collective oscillator modes by choosing the first column of the lth subblock ${\bf U}_{n}^{(l)}$ to be $||{\boldsymbol{g}} _{n}^{(l)}||_{2}^{-1}{\boldsymbol{g}} _{n}^{(l)}$. The normalized vector corresponds to the rearranged coupling strength vector of

$$\begin{eqnarray}&&{{{\boldsymbol{g}} }_{n}}={\bf P}_{n}^{\dagger }{{{\boldsymbol{ \kappa }} }_{n}}=\left( \begin{array}{ccccccccccccccc} {\boldsymbol{g}} _{n}^{(1)} \\ \vdots \\ {\boldsymbol{g}} _{n}^{({{N}_{{\rm eff}}})} \\ \end{array} \right),\end{eqnarray} \tag{ 9 }$$

with N_eff being the number of subblocks (or the number of group of oscillators). The chain-bath model can be obtained via tridiagonalization of the lth subblock

$$\begin{eqnarray}&&\tilde{{\boldsymbol{ \Omega }} }_{n}^{(l)}={{{\bf T}}^{(l)}}{{{\boldsymbol{ \Xi }} }^{(l)}}{{{\bf T}}^{(l)\dagger }}\ \end{eqnarray} \tag{ 10 }$$

with the application of the Hessenberg transform [76] via the Householder procedure. ${{{\boldsymbol{ \Xi }} }^{(l)}}$ is a tridiagonal matrix that defines the frequencies (diagonal elements) of the transformed bath modes and coupling strengths (off-diagonal elements) between the oscillators in the chain model. ${{{\bf T}}^{(l)}}$ is a Hessenberg unitary transform matrix that makes no transformation to the primary bath mode such that the first column of the matrix is ${{(1,0,\ldots ,0)}^{{\rm t}}}$. The resulting transformed bath coupling vector ${{{\bf T}}^{(l)}}\tilde{{\boldsymbol{ \kappa }} }_{n}^{(l)}$ of the lth subblock has only a single non-zero first element, which corresponds to the primary mode coupling strength. In appendix C, we provide a MATLAB [77] code for the GSH with the LP scheme.

The alternative numerical transformation from the star-bath model to the chain-bath model can be obtained by Bullaʼs recursion method [4, 5]. The two methods will be compared numerically later.

3.1. Numerical stability of the transformations

To test the numerical stability of the 1-D Wilson chain transformation methods, we perform back transformations from the chain-bath Hamiltonian to the star-bath Hamiltonian (equation (2)) by a straightforward diagonalization of ${{\tilde{{\boldsymbol{ \Omega }} }}_{n}}$. The structured SDFs are reconstructed by the back transformation and then compared with the original SDF of the chlorosome, as an example system. The chlorosome is a giant light-harvesting antenna complex of green sulfur bacteria [78]. The excitation energy is transferred within the antenna via the fluctuating environment. Various models were developed to study the system from the open quantum dynamics perspective [6, 70, 71]. Here we use the SDF of the chlorosome that some of us [70] obtained via quantum mechanics/molecular mechanics calculations that contain 253 peaks corresponding to the quantum bath oscillators. In figure 2(A), only the peaks in 1300–1700 cm⁻¹ are shown for clarity. The original SDF is plotted as a black line and filled black circles. Bullaʼs method (blue line) suffers from numerical instability as the iteration increases. Therefore, we also test extended precision (EP; 100 digits) with Bullaʼs method (green line). The GSH curve is shown as a red line and unfilled red circles. Householder and Lanczos transformations of ${{{\boldsymbol{ \Gamma }} }_{n}}$ in equation (3) are plotted in brown and purple lines, respectively. As expected, Bullaʼs original method generates a curve that deviates significantly from the original one, especially around 1550–1600 cm⁻¹. The EP improves the result but is still in disagreement with the original. The Lanczos curve seems to agree well with the original but the discrete data points do not match with the original data points in the frequency domain and it produces negative frequencies that are not shown in the figure. The GSH and Householder, on the other hand, can reproduce the original SDF with high accuracy. Both methods are based on the Householder procedure, which has an unconditional stability [79].

Zoom In Zoom Out Reset image size

Figure 2(B) indicates the Huang–Rhys (HR) factors of secondary bath oscillators with nearest-neighbor couplings in a chain, that are obtained from different methods. The HR factor ${{\chi }_{j}}$ of a harmonic oscillator with frequency ${{\omega }_{j}}$ is a normalized coupling strength given by ${{\kappa }_{j}}={{\omega }_{j}}\sqrt{{{\chi }_{j}}}$. For the secondary modes, HR factors are defined with the frequencies of the oscillators and the coupling strengths with the nearest neighbors in the chain. As one compares the results from Bullaʼs method (green cross) and Bullaʼs method with EP (blue triangle), they are significantly different from each other. This shows that Bullaʼs recursion relation is numerically unstable. Comparing the Bullaʼs and GSH methods, one can see the Bullaʼs method (green cross) produces larger HR factors than the GSH method (red circle). The GSH bath transformation can produce oscillators with frequencies distributed on a narrower band. The reason for this difference is that the unitary transformation for the chain model is not unique. Bullaʼs method does not allow the secondary bath modes to have negative interactions, while the GSH has no such constraint. Since the unitary matrix for the GSH transform is constructed via orthogonalization of the random vectors, the signs of the interactions in a single chain can vary in each transformation, but their magnitudes are invariant.

3.2. Multiple-chain-bath model

In this subsection, we apply the GSH transformation method to a couple of illustrative examples. First, to get some insights into the weakly coupled multiple chain model, we employ this method analytically to transform a bath with two oscillator modes. Then, we continue with a multidimensional example of a chlorosome.

In the chain-bath model for two oscillators, the mixing of the two modes with frequencies ${{\omega }_{1}}\leqslant {{\omega }_{2}}$ leads to a HR factor for the primary mode ${{\chi }_{{\rm c}}}$:

$$\begin{eqnarray}&&{{\chi }_{{\rm c}}}={{R}_{{\rm HR}}}{{\chi }_{1}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{R}_{{\rm HR}}}=\frac{{{(1+{{f}^{2}})}^{3}}}{{{(1+\frac{{{\omega }_{2}}}{{{\omega }_{1}}}{{f}^{2}})}^{2}}},\end{eqnarray} \tag{ 12 }$$

where $f={{\kappa }_{2}}/{{\kappa }_{1}}$ and ${{\chi }_{1}}=\kappa _{1}^{2}/\omega _{1}^{2}$. In figure 3, the normalized HR factor ${{R}_{{\rm HR}}}$ of equation (12) is plotted while varying the frequency ratio (${{\omega }_{2}}/{{\omega }_{1}}$) at fixed coupling strength ratios f, 0.5, 1 and 2. For fixed f, the values of ${{R}_{{\rm HR}}}$ decrease as the frequency ratio increases. f determines the slopes of the curves. Larger f makes the curve decrease faster as ${{\omega }_{2}}/{{\omega }_{1}}$ increases. When the oscillators have similar frequencies ${{\omega }_{2}}/{{\omega }_{1}}\simeq 1$, ${{R}_{{\rm HR}}}$ is larger than 1, which makes the chain-bath couplings stronger than the star-bath model. However, as the frequency ratio increases ${{R}_{{\rm HR}}}$ drops down and it can become arbitrarily small as the frequency difference increases. This gives a hint of how to mix the oscillator modes and reduce the coupling strength of the primary modes by forming weakly coupled multiple chains. The oscillators, which have large frequency differences, should be mixed to make the weakly coupled multiple-chain-bath models.

Zoom In Zoom Out Reset image size

Next, we apply the GSH method to the example of chlorosome. The SP approach divides the bath oscillators into a sequence of groups of oscillators, as illustrated in figure 4(A) using a different color for each group. Figure 4(B) represents the LP scheme, where only the peaks below 500 cm⁻¹ are shown for clarity. As indicated here, the oscillators are partitioned into six groups to use the LP scheme. In the LP scheme, the elements of the kth group are composed of $k+{{N}_{{\rm eff}}}j$-th ($\leqslant N$) modes, where ${{N}_{{\rm eff}}}$ is the number of groups, j is an integer ($\geqslant 0$) and N is the total number of oscillators. For example, when one has ten modes and three groups to partition, the SP approach groups modes of (1, 2, 3), (4, 5, 6) and (7, 8, 9, 10) as group 1, 2 and 3, respectively. On the other hand, the LP strategy makes the partition (1, 4, 7, 10), (2, 5, 8) and (3, 6, 9) for group 1, 2 and 3, respectively.

Zoom In Zoom Out Reset image size

In figure 5, the maximum HR factor and the corresponding coupling strength of the primary modes are plotted by varying the number of chains from one to six. The LP and SP schemes are used for this calculation. The maximum HR factor of the star-bath model is 0.0315. The single chain model has an even larger HR factor of the primary mode than the star-bath model for both partitioning schemes. The maximum HR factors from the SP scheme (blue cross) do not decrease as the number of chains increases. The figure shows, however, that the maximum HR factor (black circle) and the corresponding coupling strength (red triangle) decrease as the number of chains increases for the LP scheme. The maximum HR factor from the LP scheme with six chains is below 0.01, which corresponds to 10% of the corresponding harmonic frequency. Therefore six chains make the chain model suitable for implementation on quantum simulators, since a quantum simulator can have only a few parallel chains and, in addition, the primary mode HR factor is limited. In principle, the values can be reduced further as long as mixing modes is still possible.

Zoom In Zoom Out Reset image size

The LP and SP schemes are further compared in figure 6. This plot shows the HR factor of the primary modes of each chain. The results are indicated in blue crosses and black circles for the SP and LP schemes, respectively. As evidenced by the figure, the LP scheme gives all values below 0.01 while the SP strategy produces bigger values and the maximum value is more than four times larger than the LP values. Another important aspect of the LP scheme is that all of the primary HR factors have nearly similar values while the results of the SP scheme deviate greatly from each other.

Zoom In Zoom Out Reset image size

In this paper we show that a multiple-chain-bath model, in combination with the leaping bath partitioning scheme, may lead to a practical implementation of quantum simulators for complex open quantum dynamics. We have shown that the multiple-chain-bath model can be employed for the realization of quantum simulators for open quantum systems or for numerical studies in classical computers. Furthermore, the LP scheme can reduce the primary mode coupling strength almost homogeneously for all parallel chains. The reason is that the mixing of oscillators with large frequency differences can result in small HR factors. The two-oscillator model was presented with an analytic expression for the chain transformation and provides a hint for the bath partitioning scheme, i.e. LP. We also tested the unitary transformation algorithm that exploits the GSH transformation, and compared the results with the values from Bullaʼs recursion method [5]. The GSH transformation can produce smaller HR factors of the secondary bath modes with oscillators being in a narrower band. The numerical stability of Bullaʼs method was discussed and the GSH method was shown to be numerically stable.

Our bath transformation method could be useful for the perturbative approaches as well, because of the resulting weak system-bath couplings. The effective spectral density can be obtained based on the chain-model transformation. It can also be used for the reduced density matrix methods [31, 32] for the simulation of non-Markovian dynamics. The effective bath approach with the parallel chain-bath model can also be useful for the HEOM method. The effective bath spectral density can reduce the number of Drude–Lorentzian peaks, then the HEOM method can handle larger systems [21]. Further work in this direction will be conducted.

We acknowledge Jarrod McClean, Florian Schroeder, Dr. Christoph Kreisbeck and Dr. Semion K Saikin for discussions. J H, T F and A A-G acknowledge support from the Center for Excitonics, an Energy Frontier Research Center funded by the US Department of Energy, Office of Science and Office of Basic Energy Sciences under award DE-SC0001088. J H and S M and A A-G also acknowledge Defense Threat Reduction Agency grant HDTRA1-10-1-0046 and the Air Force Office of Scientific Research grant FA9550-12-1-0046. Further, A A-G is grateful for the support from Defense Advanced Research Projects Agency grant N66001-10-1-4063, and the Corning Foundation for their generous support. M-H Y acknowledges the support from the National Basic Research Program of China Grant 2011CBA00300, 2011CBA00301, the National Natural Science Foundation of China Grant 61033001, 61361136003, and the Youth 1000-talent program.

We explain here in detail how to construct the unitary transform matrix ${{{\bf U}}_{n}}$ of GSH transformation in equation (6) for ${{\tilde{{\boldsymbol{ \Omega }} }}_{n}}$. A MATLAB [77] script is given in appendix C.

• 1.  
  Set the first column of ${{{\bf U}}_{n}}$ to be the normalized coupling strength vector $||{{{\boldsymbol{ \kappa }} }_{n}}||_{2}^{-1}{{{\boldsymbol{ \kappa }} }_{n}}$.
• 2.  
  Assign random vectors to the remaining columns (from 2 to N). We use a normal distribution with the first column to be a mean vector for the random vectors.
• 3.  
  Perform the Gram–Schmidt orthogonalization to ${{{\bf U}}_{n}}$.
• 4.  
  Compute ${\bf U}_{n}^{\dagger }{{{\boldsymbol{ \Omega }} }_{n}}{{{\bf U}}_{n}}$ for ${{\tilde{{\boldsymbol{ \Omega }} }}_{n}}$.

Here we represent a permutation matrix ${{{\bf P}}^{\dagger }}$, as an example. It permutes a vector of length 4 (${\boldsymbol{v}}$) to arrange odd and even elements sequentially,

$$\begin{eqnarray}{{{\bf P}}^{\dagger }}{\boldsymbol{v}} =\left( \begin{array}{ccccccccccccccc} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right)\left( \begin{array}{ccccccccccccccc} {{v}_{1}} \\ {{v}_{2}} \\ {{v}_{3}} \\ {{v}_{4}} \\ \end{array} \right)=\left( \begin{array}{ccccccccccccccc} {{v}_{1}} \\ {{v}_{3}} \\ {{v}_{2}} \\ {{v}_{4}} \\ \end{array} \right).\end{eqnarray} \tag{ B.1 }$$

We provide here a MATLAB [77] script for the LP scheme with the GSH procedure. The script produces a single chain for a subblock, which is assigned by the user. Multiple chains can be obtained by using this script for each subblock. The details of the script can be found in the comments of the script, which are indicated with %.

• function [HRF,Dn,Xi]=LP(SDF,Sub,Neff)
• %-Inputs-
• %SDF: Spectral density function = [frequency spectral-density]
• %Sub: Subblock of interest (starting index in the full SDF)
• %Neff: Number of subblocks (index step size for LP scheme)

• %-Outputs-
• %HRF: HR factor of the primary mode of the subblock
• %Dn: norm of the coupling vector (equivalent to the coupling strength of the primary mode)
• %Xi: Tridiagonal coupling bath coupling matrix
• % equation (10), OmegaTilde=T*Xi*T', T*T'=I

• %Scale SDF with 1/pi
• SDF(:,2)=SDF(:,2)/pi;
• %dimension of the SDF
• dim=max(size(SDF));

• %LP scheme for the subgroup
• SDF=spd(Sub:Neff:dim,:);
• %new dimension
• dim=max(size(SDF));

• %t1: first column of U, U*U'=I
• t1=sqrt(SDF(:,2));
• Dn=norm(t1);
• t1=t1/Dn;

• %Random vectors for the Gram-Schmidt orthogonalization
• %are assigned to columns from 2 to N.
• U0=randn(dim,dim-1);
• U0=[t1 U0];
• for k=2:dim
• U0(:,k)=U0(:,k)+t1;
• end
• %Gram-Schmidt orthogonalization, user should provide the GS routine.
• U=GS(U0);
• %OmegaTilde in equation 8
• OmegaTilde=U'*diag(SDF(:,1))*U;
• %Hessenberg transformation
• [T,Xi]=hess(OmegaTilde);
• HRF=Dn2ˆ/Xi(1,1)2ˆ;

• end