The flow chart will go here
Contributions are welcome! Don't hesitate to contact us if you
You can contact us by raising an issue on Github, or by writing to one of the developers.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
The flow chart will go here
Contributions are welcome! Don't hesitate to contact us if you
You can contact us by raising an issue on Github, or by writing to one of the developers.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Here we give some context on the examples provided in MPSDynamics/examples/anderson_model_double.jl and MPSDynamics/examples/anderson_model_interleaved.jl. In these two examples, we use the fermionic chain mapping proposed in [khon_efficient_2021] to show how to perform tensor network simulations of the Single Impurity Anderson Model (SIAM) with the MPSDynamics.jl library.
Before giving the code implementation, a short recap on the problem statement and on the fermionic mapping used to solved it. The SIAM Hamiltonian is defined as:
\[ \hat H^\text{SIAM} = \hat H_\text{loc} + \hat H_\text{hyb} + \hat H_\text{cond} = \overbrace{\epsilon_d \hat d^\dagger \hat d}^{\hat H_\text{loc}} + \underbrace{\sum_{k} V_k \Big( \hat d^\dagger \hat c_k + \hat c_k^\dagger \hat d \Big)}_{H_\text{hyb}} + \underbrace{\sum_k \epsilon_k \hat c_k^\dagger \hat c_k}_{H_I^\text{chain}}.\]
All of the operators obey to the usual fermionic anti-commutation relations: $\{\hat c_i, \hat c_j^\dagger \} = \delta_{ij}$, $\{\hat c_i, \hat c_j \} =\{\hat c_i^\dagger, \hat c_j^\dagger \} =0$ $\forall i,j$. The chain mapping is based on a thermofield-like transformation [devega_thermo_2015], performed with fermions: ancillary fermionic operators $\hat c_{2k}$ are defined, one for each of the original fermionic modes $\hat c_{1k}$. A Bogoliubov transformation is then applied, so that two new fermionic modes $\hat f_{1k}$ and $\hat f_{2k}$ are defined as a linear combination of $\hat c_{1k}$ and $\hat c_{2k}$. Two chains are defined: the chain labelled $1$ for the empty modes:
\[\hat f_{1k}=e^{-iG} \hat c_k e^{iG}= \cosh(\theta_k) \hat c_{1k} -\sinh(\theta_k) \hat c_{2k}^\dagger \\ +
Here we give some context on the examples provided in MPSDynamics/examples/anderson_model_double.jl and MPSDynamics/examples/anderson_model_interleaved.jl. In these two examples, we use the fermionic chain mapping proposed in [khon_efficient_2021] to show how to perform tensor network simulations of the Single Impurity Anderson Model (SIAM) with the MPSDynamics.jl library.
Before giving the code implementation, a short recap on the problem statement and on the fermionic mapping used to solved it. The SIAM Hamiltonian is defined as:
\[ \hat H^\text{SIAM} = \hat H_\text{loc} + \hat H_\text{hyb} + \hat H_\text{cond} = \overbrace{\epsilon_d \hat d^\dagger \hat d}^{\hat H_\text{loc}} + \underbrace{\sum_{k} V_k \Big( \hat d^\dagger \hat c_k + \hat c_k^\dagger \hat d \Big)}_{H_\text{hyb}} + \underbrace{\sum_k \epsilon_k \hat c_k^\dagger \hat c_k}_{H_I^\text{chain}}.\]
All of the operators obey to the usual fermionic anti-commutation relations: $\{\hat c_i, \hat c_j^\dagger \} = \delta_{ij}$, $\{\hat c_i, \hat c_j \} =\{\hat c_i^\dagger, \hat c_j^\dagger \} =0$ $\forall i,j$. The chain mapping is based on a thermofield-like transformation [devega_thermo_2015], performed with fermions: ancillary fermionic operators $\hat c_{2k}$ are defined, one for each of the original fermionic modes $\hat c_{1k}$. A Bogoliubov transformation is then applied, so that two new fermionic modes $\hat f_{1k}$ and $\hat f_{2k}$ are defined as a linear combination of $\hat c_{1k}$ and $\hat c_{2k}$. Two chains are defined: the chain labelled $1$ for the empty modes:
\[\hat f_{1k}=e^{-iG} \hat c_k e^{iG}= \cosh(\theta_k) \hat c_{1k} -\sinh(\theta_k) \hat c_{2k}^\dagger \\ \hat f_{2k}=e^{-iG} \hat c_k e^{iG}= \cosh(\theta_k) \hat c_{1k} +\sinh(\theta_k) \hat c_{2k}^\dagger.\]
We remark that this is the same Bogoliubov transformation used in the thermofield[devega_thermo_2015] for the bosonic case: the only thing that changes is a plus sign, that takes into account the fermionic anti-commutation relations. With these new fermionic operators we obtain the thermofield-transformed Hamiltonian, where the system interacts with two environments at zero temperature, allowing for pure state simulations (and thus the employement of MPS).
The thermofield-transformed Hamiltonian is then mapped on two chains, defined and constructed using the TEDOPA chain mapping: the chain labelled $1$ is for the empty modes, the chain labelled $2$ for the filled modes. The following relations are used to define the functions equivalent to the spectral density of the bosonic case, one for each chain:
\[\begin{split} &V_{1k} = V_{k} \sin \theta_k = \sqrt{\frac{1}{e^{\beta \epsilon_k}+1}} \\ &V_{2k} = V_{k} \cos \theta_k = \sqrt{\frac{1}{e^{-\beta \epsilon_k}+1}}, @@ -266,4 +266,4 @@ end # Display the plots -plot(p2, p3, p4, p5, p1, layout = (3, 2), size = (1400, 1200))
________________
Kohn, L.; Santoro, G. E. Efficient mapping for anderson impurity problems with matrix product states. Phys. Rev. B 2021, 104 (1), 014303. https://doi.org/10.1103/PhysRevB.104.014303.
de Vega, I.; Banuls, M-.C. Thermofield-based chain-mapping approach for open quantum systems. Phys. Rev. A 2015, 92 (5), 052116. https://doi.org/10.1103/PhysRevA.92.052116.
Kohn L.; Santoro G. E. J. Stat. Mech. (2022) 063102 DOI 10.1088/1742-5468/ac729b.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
________________
Kohn, L.; Santoro, G. E. Efficient mapping for anderson impurity problems with matrix product states. Phys. Rev. B 2021, 104 (1), 014303. https://doi.org/10.1103/PhysRevB.104.014303.
de Vega, I.; Banuls, M-.C. Thermofield-based chain-mapping approach for open quantum systems. Phys. Rev. A 2015, 92 (5), 052116. https://doi.org/10.1103/PhysRevA.92.052116.
Kohn L.; Santoro G. E. J. Stat. Mech. (2022) 063102 DOI 10.1088/1742-5468/ac729b.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Here we give some context on the example script provided in MPSDynamics/example/bath-observables.jl.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Here we give some context on the example script provided in MPSDynamics/example/bath-observables.jl.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
The Pure-Dephasing Model describes a two-level system interacting linearly with an environment characterised by a spectral density (SD) $J(\omega)$. The coupling only acts on diagonal terms through the $\sigma_z$ operator. The Hamiltonian reads
\[ \hat{H} = \frac{\omega_0}{2}\hat{\sigma}_z + \int_0^{+\infty}\omega \hat{a}^\dagger_\omega\hat{a}_\omega \mathrm{d}\omega + \frac{\hat{\sigma}_z}{2}\int_0^{+\infty}\sqrt{J(\omega)}(\hat{a}_\omega + \hat{a}^\dagger_\omega)\mathrm{d}\omega\]
Although this interaction will not change the population of the two-level system, the coherences between the two states will vary due to the environment. Introducing the two-level system reduced density matrix $\rho_{ij}(t)$ with $i,j \in (0,1)$, the diagonal terms $\rho_{ii}(t)$ are the populations of the states and the anti-diagonal terms $\rho_{ij}(t)$ with $i \neq j$ are the coherences between the two states. The effect of the $\sigma_z$ bath interaction is to decouple the two states $|0\rangle$ and $|1\rangle$.
The density matrix can be calculated within the MPS formalism with the MPO $|\psi\rangle \langle \psi|$. Tracing out the environment leads to the reduced density matrix. This can be done with the function MPSDynamics.rhoreduced_1site.
An analytical formula can be found for the decoherence function $\Gamma(t)$, taking into account the SD as well as the temperature of the environment [breuer]
\[ \Gamma(t) = - \int_0^{\omega_c} \mathrm{d} \omega J(\omega)\frac{(1 - \cos(\omega t))}{\omega^2} \coth(β\omega/2) ,\]
with $\beta = (k_B T)^{-1}$. For the case where $\beta \longrightarrow \infty$, the integral reads
\[ \Gamma(t) = - \int_0^{\omega_c} \mathrm{d} \omega J(\omega)\frac{(1 - \cos(\omega t))}{\omega^2} ,\]
The time-dependent anti-diagonal terms of the reduced density matrix are then expressed as
\[ \rho_{12}(t) = \rho_{21}(t)^* =\rho_{12}(0) \exp(\Gamma(t)) \]
Setting up the initial two-level system as a cat state $|\psi\rangle_S(0) = \frac{|0\rangle \pm |1\rangle}{\sqrt{2}}$, this leads to $\rho_{12}(0)=\frac{1}{2}$.
Here we break out and comment the script in MPSDynamics/examples/puredephasing.jl and MPSDynamics/examples/puredephasing_temperature.jl to show how to simulate this model with an Ohmic SD (hard cut-off) using the T-TEDOPA method as implemented in MPSDynamics.jl.
The T-TEDOPA method relies on a truncated chain mapping that transform the initial Hamiltonian into
\[ \hat{H} = \frac{\omega_0}{2} \hat{\sigma}_z + c_0 \frac{\hat{\sigma}_z}{2}(\hat{b}_0^\dagger + \hat{b}_0) + \sum_{i=0}^{N-1} t_i (\hat{b}_{i+1}^\dagger \hat{b}_i + \mathrm{h.c.}) + \sum_{i=0}^{N-1} \epsilon_i \hat{b}_i^\dagger \hat{b}_i\]
First, we load the MPSdynamics.jl package to be able to perform the simulation, the Plots.jl one to plot the results, the LaTeXStrings.jl one to be able to use $\LaTeX$ in the plots and eventually QuadGK.jl to perform the analytical integral calculations. A notation is introduced for +∞.
using MPSDynamics, Plots, LaTeXStrings, QuadGK
+Pure-Dephasing · MPSDynamics.jl Pure-Dephasing
Context
The Pure-Dephasing Model describes a two-level system interacting linearly with an environment characterised by a spectral density (SD) $J(\omega)$. The coupling only acts on diagonal terms through the $\sigma_z$ operator. The Hamiltonian reads
\[ \hat{H} = \frac{\omega_0}{2}\hat{\sigma}_z + \int_0^{+\infty}\omega \hat{a}^\dagger_\omega\hat{a}_\omega \mathrm{d}\omega + \frac{\hat{\sigma}_z}{2}\int_0^{+\infty}\sqrt{J(\omega)}(\hat{a}_\omega + \hat{a}^\dagger_\omega)\mathrm{d}\omega\]
Although this interaction will not change the population of the two-level system, the coherences between the two states will vary due to the environment. Introducing the two-level system reduced density matrix $\rho_{ij}(t)$ with $i,j \in (0,1)$, the diagonal terms $\rho_{ii}(t)$ are the populations of the states and the anti-diagonal terms $\rho_{ij}(t)$ with $i \neq j$ are the coherences between the two states. The effect of the $\sigma_z$ bath interaction is to decouple the two states $|0\rangle$ and $|1\rangle$.
The density matrix can be calculated within the MPS formalism with the MPO $|\psi\rangle \langle \psi|$. Tracing out the environment leads to the reduced density matrix. This can be done with the function MPSDynamics.rhoreduced_1site.
An analytical formula can be found for the decoherence function $\Gamma(t)$, taking into account the SD as well as the temperature of the environment [breuer]
\[ \Gamma(t) = - \int_0^{\omega_c} \mathrm{d} \omega J(\omega)\frac{(1 - \cos(\omega t))}{\omega^2} \coth(β\omega/2) ,\]
with $\beta = (k_B T)^{-1}$. For the case where $\beta \longrightarrow \infty$, the integral reads
\[ \Gamma(t) = - \int_0^{\omega_c} \mathrm{d} \omega J(\omega)\frac{(1 - \cos(\omega t))}{\omega^2} ,\]
The time-dependent anti-diagonal terms of the reduced density matrix are then expressed as
\[ \rho_{12}(t) = \rho_{21}(t)^* =\rho_{12}(0) \exp(\Gamma(t)) \]
Setting up the initial two-level system as a cat state $|\psi\rangle_S(0) = \frac{|0\rangle \pm |1\rangle}{\sqrt{2}}$, this leads to $\rho_{12}(0)=\frac{1}{2}$.
Here we break out and comment the script in MPSDynamics/examples/puredephasing.jl and MPSDynamics/examples/puredephasing_temperature.jl to show how to simulate this model with an Ohmic SD (hard cut-off) using the T-TEDOPA method as implemented in MPSDynamics.jl.
The T-TEDOPA method relies on a truncated chain mapping that transform the initial Hamiltonian into
\[ \hat{H} = \frac{\omega_0}{2} \hat{\sigma}_z + c_0 \frac{\hat{\sigma}_z}{2}(\hat{b}_0^\dagger + \hat{b}_0) + \sum_{i=0}^{N-1} t_i (\hat{b}_{i+1}^\dagger \hat{b}_i + \mathrm{h.c.}) + \sum_{i=0}^{N-1} \epsilon_i \hat{b}_i^\dagger \hat{b}_i\]
The code
First, we load the MPSdynamics.jl package to be able to perform the simulation, the Plots.jl one to plot the results, the LaTeXStrings.jl one to be able to use $\LaTeX$ in the plots and eventually QuadGK.jl to perform the analytical integral calculations. A notation is introduced for +∞.
using MPSDynamics, Plots, LaTeXStrings, QuadGK
const ∞ = Inf
We then define variables for the physical parameters of the simulation. Among these, two are convergence parameters:
d is the number of states we retain for the truncated harmonic oscillators representation of environmental modesN is the number of chain (environmental) modes we keep. This parameters determines the maximum simulation time of the simulation: indeed excitations that arrive at the end of the chain are reflected towards the system and can lead to unphysical results
The value of β determines whether the environment is thermalized or not. The example as it is is the zero-temperature case. For the finite-temperature case, 'β = ∞' has be commented instead of the line above that precises a β value.
#----------------------------
# Physical parameters
@@ -93,4 +93,4 @@
plot(time_analytical, t->Decoherence_ohmic(t), label="Analytics", title=L"Pure Dephasing, Ohmic $s=%$s$, $\beta = %$β ~\mathrm{K}$", linecolor=:black, xlabel="Time (arb. units)", ylabel=L"Coherence $|\rho_{12}(t)|$", linewidth=4, titlefontsize=16, legend=:best, legendfontsize=16, xguidefontsize=16, yguidefontsize=16, tickfontsize=10)
-plot!(dat["data/times"], ρ12, lw=4, ls=:dash, label="Numerics")
Bibliography
- breuer
Breuer, H.;Petruccione, F. The Theory of Open Quantum Systems; Oxford University Press, 2002.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Breuer, H.;Petruccione, F. The Theory of Open Quantum Systems; Oxford University Press, 2002.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
The Spin-Boson Model (SBM) is a prototypical model in the theory of open quantum systems where a two level system interacts linearly with a bosonic bath
\[ \hat{H} = \frac{\omega_0}{2}\hat{\sigma}_z + \Delta\hat{\sigma}_x + \int_0^{+\infty}\omega \hat{a}^\dagger_\omega\hat{a}_\omega \mathrm{d}\omega + \hat{\sigma}_x\int_0^{+\infty}\sqrt{J(\omega)}(\hat{a}_\omega + \hat{a}^\dagger_\omega)\mathrm{d}\omega\]
Even though this model is fairly simple it is physically very rich and it is not analytically solvable. For these reason it has become a test-bed for numerical methods simulating open quantum systems dynamics in the non-perturbative non-Markovian regime.
For instance when the SD is Ohmic, this model presents a phase transition between a so called localised and a delocalised phase for $\alpha \approx 1.2$.
Here we break out and comment the script in MPSDynamics/examples/sbm_zero_temperature.jl to show how to simulate this model with an Ohmic SD (hard cut-off) using the T-TEDOPA method as implemented in MPSDynamics.jl.
The T-TEDOPA method relies on a truncated chain mapping that transform the initial Hamiltonian into
\[ \hat{H} = \frac{\omega_0}{2} \hat{\sigma}_z + \Delta \hat{\sigma}_x + c_0 \hat{\sigma}_x(\hat{b}_0^\dagger + \hat{b}_0) + \sum_{i=0}^{N-1} t_i (\hat{b}_{i+1}^\dagger \hat{b}_i + \mathrm{h.c.}) + \sum_{i=0}^{N-1} \epsilon_i \hat{b}_i^\dagger \hat{b}_i \]
First a multi-line comment introduces the model and the aim of the script.
#=
+The Spin-Boson Model · MPSDynamics.jl The Spin-Boson Model
Context
The Spin-Boson Model (SBM) is a prototypical model in the theory of open quantum systems where a two level system interacts linearly with a bosonic bath
\[ \hat{H} = \frac{\omega_0}{2}\hat{\sigma}_z + \Delta\hat{\sigma}_x + \int_0^{+\infty}\omega \hat{a}^\dagger_\omega\hat{a}_\omega \mathrm{d}\omega + \hat{\sigma}_x\int_0^{+\infty}\sqrt{J(\omega)}(\hat{a}_\omega + \hat{a}^\dagger_\omega)\mathrm{d}\omega\]
Even though this model is fairly simple it is physically very rich and it is not analytically solvable. For these reason it has become a test-bed for numerical methods simulating open quantum systems dynamics in the non-perturbative non-Markovian regime.
For instance when the SD is Ohmic, this model presents a phase transition between a so called localised and a delocalised phase for $\alpha \approx 1.2$.
Here we break out and comment the script in MPSDynamics/examples/sbm_zero_temperature.jl to show how to simulate this model with an Ohmic SD (hard cut-off) using the T-TEDOPA method as implemented in MPSDynamics.jl.
The T-TEDOPA method relies on a truncated chain mapping that transform the initial Hamiltonian into
\[ \hat{H} = \frac{\omega_0}{2} \hat{\sigma}_z + \Delta \hat{\sigma}_x + c_0 \hat{\sigma}_x(\hat{b}_0^\dagger + \hat{b}_0) + \sum_{i=0}^{N-1} t_i (\hat{b}_{i+1}^\dagger \hat{b}_i + \mathrm{h.c.}) + \sum_{i=0}^{N-1} \epsilon_i \hat{b}_i^\dagger \hat{b}_i \]
The code
First a multi-line comment introduces the model and the aim of the script.
#=
Example of a zero-temperature Spin-Boson Model with an hard cut-off Ohmic spectral density J(ω) = 2αω when ω < ωc and 0 otherwise
The dynamics is simulated using the T-TEDOPA method that maps the normal modes environment into a non-uniform tight-binding chain.
@@ -78,4 +78,4 @@
method == :DTDVP && heatmap(dat["data/times"], collect(0:N+1), dat["data/bonddims"], xlabel=L"t",ylabel="bond index")
-heatmap(dat["data/times"], collect(1:N), abs.(dat["data/SXdisp"][1,:,:]), xlabel=L"t",ylabel="chain mode")
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
To simulate a drive or a laser pulse, a time-dependent Hamiltonian is often needed. Here we explain thanks to the script in MPSDynamics/examples/sbm_Htimedependent.jl how this type of Hamiltonian can be integrated within the method and applied on a wavefunction. We will take a Spin-Boson Model (SBM) Hamiltonian with an Ohmic spectral density. However, the time-dependency is not model specific and can be adapated for other models. For more information about the SBM, see the dedicated example. For this example, we simulate a drive of the form
\[ \hat{H}_\text{drive}(t) = \epsilon \hat{\sigma}_x \sin(\omega_\text{drive} t)\]
with $\epsilon = \frac{2 \pi}{T_\text{Rabi}}$ the amplitude of the drive, $T_\text{Rabi}$ the Rabi period, and $\omega_\text{drive}$ the frequency of the drive. The drive is set up to be on resonance with the two-level system.
First we load the MPSdynamics.jl package to be able to perform the simulation, the Plots.jl one to plot the results, and the LaTeXStrings.jl one to be able to use $\LaTeX$ in the plots. The function MPSDynamics.disp is also imported.
using MPSDynamics, Plots, LaTeXStrings
+Time-dependent Hamiltonian · MPSDynamics.jl Time-dependent Hamiltonian
Context
To simulate a drive or a laser pulse, a time-dependent Hamiltonian is often needed. Here we explain thanks to the script in MPSDynamics/examples/sbm_Htimedependent.jl how this type of Hamiltonian can be integrated within the method and applied on a wavefunction. We will take a Spin-Boson Model (SBM) Hamiltonian with an Ohmic spectral density. However, the time-dependency is not model specific and can be adapated for other models. For more information about the SBM, see the dedicated example. For this example, we simulate a drive of the form
\[ \hat{H}_\text{drive}(t) = \epsilon \hat{\sigma}_x \sin(\omega_\text{drive} t)\]
with $\epsilon = \frac{2 \pi}{T_\text{Rabi}}$ the amplitude of the drive, $T_\text{Rabi}$ the Rabi period, and $\omega_\text{drive}$ the frequency of the drive. The drive is set up to be on resonance with the two-level system.
The code
First we load the MPSdynamics.jl package to be able to perform the simulation, the Plots.jl one to plot the results, and the LaTeXStrings.jl one to be able to use $\LaTeX$ in the plots. The function MPSDynamics.disp is also imported.
using MPSDynamics, Plots, LaTeXStrings
import MPSDynamics: disp
We then define variables for the physical parameters of the simulation. Among these, three are convergence parameters:
d is the number of states we retain for the truncated harmonic oscillators representation of environmental modesN is the number of chain (environmental) modes we keep. This parameters determines the maximum simulation time of the simulation: indeed excitations that arrive at the end of the chain are reflected towards the system and can lead to unphysical results
The variable Ndrive represents the site of the MPO where the operator of the time-dependent part acts. For this example, the two-level system is at the first site of the MPS.
#----------------------------
# Physical parameters
@@ -79,4 +79,4 @@
# Plots
#----------
-plot(dat["data/times"], dat["data/sz"], label="Dmax = $(D...)", xlabel=L"t",ylabel=L"\sigma_z", title="")
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
The MPSDynamics.jl package provides an easy to use interface for performing tensor network simulations on matrix product states (MPS) and tree tensor network (TTN) states. Written in the Julia programming language, MPSDynamics.jl is a versatile open-source package providing a choice of several variants of the Time-Dependent Variational Principle (TDVP) method for time evolution. The package also provides strong support for the measurement of observables, as well as the storing and logging of data, which makes it a useful tool for the study of many-body physics. The package has been developed with the aim of studying non-Markovian open system dynamics at finite temperature using the state-of-the-art numerically exact Thermalized-Time Evolving Density operator with Orthonormal Polynomials Algorithm (T-TEDOPA) based on environment chain mapping. However the methods implemented can equally be applied to other areas of physics.
The documentation is currently undergoing massive restructurations/improvement. It's a work in progress until the next release scheduled for May, 2024.
The package may be installed by typing the following into a Julia REPL
] add https://github.com/shareloqs/MPSDynamics.gitIf you use the package in your research, please consider citing it. You can add the Zenodo record to your BibTex file:
@misc{mpsdynamics_zenodo2021,
+Introduction · MPSDynamics.jl Introduction
The MPSDynamics.jl package provides an easy to use interface for performing tensor network simulations on matrix product states (MPS) and tree tensor network (TTN) states. Written in the Julia programming language, MPSDynamics.jl is a versatile open-source package providing a choice of several variants of the Time-Dependent Variational Principle (TDVP) method for time evolution. The package also provides strong support for the measurement of observables, as well as the storing and logging of data, which makes it a useful tool for the study of many-body physics. The package has been developed with the aim of studying non-Markovian open system dynamics at finite temperature using the state-of-the-art numerically exact Thermalized-Time Evolving Density operator with Orthonormal Polynomials Algorithm (T-TEDOPA) based on environment chain mapping. However the methods implemented can equally be applied to other areas of physics.
Warning The documentation is currently undergoing massive restructurations/improvement. It's a work in progress until the next release scheduled for May, 2024.
Installation
The package may be installed by typing the following into a Julia REPL
] add https://github.com/shareloqs/MPSDynamics.git
Table of Contents
- Introduction
- User Guide
- Theoretical Background
- Chain-Mapping of bosonic environments
- Zero Temperature
- Finite Temperature with T-TEDOPA
- Finite temperature with the thermofield transformation
- Computation of the chain coefficients
- Tensor Networks
- Time-Dependent Variational Principle
- Bibliography
- List of all methods
- Developers
Citation
If you use the package in your research, please consider citing it. You can add the Zenodo record to your BibTex file:
@misc{mpsdynamics_zenodo2021,
title = {shareloqs/{MPSDynamics}},
shorttitle = {{MPSDynamics.jl}},
url = {https://zenodo.org/record/5106435},
@@ -8,4 +8,4 @@
author = {Dunnett, Angus and Lacroix, Thibaut and Le Dé, Brieuc and Riva, Angela},
year = {2021},
doi = {10.5281/zenodo.5106435},
-}
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
MPSDynamics.OneSiteObservable — MethodOneSiteObservable(name,op,sites)Computes the local expectation value of the one-site operator op on the specified sites. Used to define one-site observables that are obs and convobs parameters for the runsim function.
MPSDynamics.MPOtoVector — MethodMPOtoVector(mpo::MPO)Convert an ITensors chain MPO into a form compatible with MPSDynamics
MPSDynamics.addchild! — Methodaddchild!(tree::Tree, id::Int)Add child to node id of tree.
MPSDynamics.addchildren! — Methodaddchildren!(tree::Tree, id::Int, n::Int)Add n children to node id of tree.
MPSDynamics.chaincoeffs_ohmic — Methodchaincoeffs_ohmic(N, α, s; ωc=1, soft=false)Generate chain coefficients $[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$ for an Harmonic bath at zero temperature with a power law spectral density given by:
soft cutoff: $J(ω) = 2παω_c (\frac{ω}{ω_c})^s \exp(-ω/ω_c)$
hard cutoff: $J(ω) = 2παω_c (\frac{ω}{ω_c})^s θ(ω-ω_c)$
The coefficients parameterise the chain Hamiltonian
$H = H_S + c_0 A_S⊗B_0+\sum_{i=0}^{N-1}t_i (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N} ϵ_ib_i^\dagger b_i$
which is unitarily equivalent (before the truncation to N sites) to
$H = H_S + A_S⊗\int_0^∞dω\sqrt{\frac{J(ω)}{π}}B_ω + \int_0^∞dωωb_ω^\dagger b_ω$
MPSDynamics.chainmps — Methodchainmps(N::Int, site::Int, numex::Int)Generate an MPS with numex excitations on site
The returned MPS will have bond-dimensions and physical dimensions numex+1
MPSDynamics.chainmps — Methodchainmps(N::Int, sites::Vector{Int}, numex::Int)Generate an MPS with numex excitations of an equal super-position over sites
The returned MPS will have bond-dimensions and physical dimensions numex+1
MPSDynamics.chainprop — Methodchainprop(t, cparams...)Propagate an excitation placed initially on the first site of a tight-binding chain with parameters given by cparams for a time t and return occupation expectation for each site.
MPSDynamics.dynamap — Methoddynamap(ps1,ps2,ps3,ps4)Calulate complete dynamical map to time step at which ps1, ps2, ps3 and ps4 are specified.
MPSDynamics.electron2kmps — Functionelectronkmps(N::Int, k::Vector{Int}, spin=:Up, chainparams=[fill(1.0,N), fill(1.0,N-1)])Generate an MPS with 2 electrons in k-states k1 and k2.
MPSDynamics.electronkmps — Functionelectronkmps(N::Int, k::Int, spin=:Up, chainparams=[fill(1.0,N), fill(1.0,N-1)])Generate an MPS for an electron with momentum k.
MPSDynamics.elementmpo — Methodelementmpo(M, el...)Return the element of the MPO M for the set of physical states el...
MPSDynamics.elementmps — Methodelementmps(A, el...)Return the element of the MPS A for the set of physical states el...
Examples
julia> A = chainmps(6, [2,4], 1);
+Methods · MPSDynamics.jl List of all methods
MPSDynamics.OneSiteObservable — MethodOneSiteObservable(name,op,sites)
Computes the local expectation value of the one-site operator op on the specified sites. Used to define one-site observables that are obs and convobs parameters for the runsim function.
sourceMPSDynamics.MPOtoVector — MethodMPOtoVector(mpo::MPO)
Convert an ITensors chain MPO into a form compatible with MPSDynamics
sourceMPSDynamics.addchild! — Methodaddchild!(tree::Tree, id::Int)
Add child to node id of tree.
sourceMPSDynamics.addchildren! — Methodaddchildren!(tree::Tree, id::Int, n::Int)
Add n children to node id of tree.
sourceMPSDynamics.chaincoeffs_ohmic — Methodchaincoeffs_ohmic(N, α, s; ωc=1, soft=false)
Generate chain coefficients $[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$ for an Harmonic bath at zero temperature with a power law spectral density given by:
soft cutoff: $J(ω) = 2παω_c (\frac{ω}{ω_c})^s \exp(-ω/ω_c)$
hard cutoff: $J(ω) = 2παω_c (\frac{ω}{ω_c})^s θ(ω-ω_c)$
The coefficients parameterise the chain Hamiltonian
$H = H_S + c_0 A_S⊗B_0+\sum_{i=0}^{N-1}t_i (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N} ϵ_ib_i^\dagger b_i$
which is unitarily equivalent (before the truncation to N sites) to
$H = H_S + A_S⊗\int_0^∞dω\sqrt{\frac{J(ω)}{π}}B_ω + \int_0^∞dωωb_ω^\dagger b_ω$
sourceMPSDynamics.chainmps — Methodchainmps(N::Int, site::Int, numex::Int)
Generate an MPS with numex excitations on site
The returned MPS will have bond-dimensions and physical dimensions numex+1
sourceMPSDynamics.chainmps — Methodchainmps(N::Int, sites::Vector{Int}, numex::Int)
Generate an MPS with numex excitations of an equal super-position over sites
The returned MPS will have bond-dimensions and physical dimensions numex+1
sourceMPSDynamics.chainprop — Methodchainprop(t, cparams...)
Propagate an excitation placed initially on the first site of a tight-binding chain with parameters given by cparams for a time t and return occupation expectation for each site.
sourceMPSDynamics.dynamap — Methoddynamap(ps1,ps2,ps3,ps4)
Calulate complete dynamical map to time step at which ps1, ps2, ps3 and ps4 are specified.
sourceMPSDynamics.electron2kmps — Functionelectronkmps(N::Int, k::Vector{Int}, spin=:Up, chainparams=[fill(1.0,N), fill(1.0,N-1)])
Generate an MPS with 2 electrons in k-states k1 and k2.
sourceMPSDynamics.electronkmps — Functionelectronkmps(N::Int, k::Int, spin=:Up, chainparams=[fill(1.0,N), fill(1.0,N-1)])
Generate an MPS for an electron with momentum k.
sourceMPSDynamics.elementmpo — Methodelementmpo(M, el...)
Return the element of the MPO M for the set of physical states el...
sourceMPSDynamics.elementmps — Methodelementmps(A, el...)
Return the element of the MPS A for the set of physical states el...
Examples
julia> A = chainmps(6, [2,4], 1);
julia> elementmps(A, 1, 2, 1, 1, 1, 1)
0.7071067811865475
@@ -11,17 +11,17 @@
0.0
julia> elementmps(A, 1, 1, 1, 1, 1, 1)
-0.0
sourceMPSDynamics.entanglemententropy — Methodentanglemententropy(A)
For a list of tensors A representing a right orthonormalized MPS, compute the entanglement entropy for a bipartite cut for every bond.
sourceMPSDynamics.findchainlength — Methodfindchainlength(T, cparams...; eps=10^-6)
Estimate length of chain required for a particular set of chain parameters by calculating how long an excitation on the first site takes to reach the end. The chain length is given as the length required for the excitation to have just reached the last site after time T. The initial number of sites in cparams has to be larger than the findchainlength result.
sourceMPSDynamics.findchild — Methodfindchild(node::TreeNode, id::Int)
Return integer corresponding to the which number child site id is of node.
sourceMPSDynamics.hbathchain — Methodhbathchain(N::Int, d::Int, chainparams, longrangecc...; tree=false, reverse=false, coupletox=false)
Generate MPO representing a tight-binding chain of N oscillators with d Fock states each. Chain parameters are supplied in the standard form: chainparams $=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$. The output does not itself represent a complete MPO but will possess an end which is open and should be attached to another tensor site, usually representing the system.
Arguments
reverse: If reverse=true create a chain were the last (i.e. Nth) site is the site which couples to the systemcoupletox: Used to choose the form of the system coupling. coupletox=true gives a non-number conserving coupling of the form $H_{\text{I}}= A_{\text{S}}(b_{0}^\dagger + b_0)$ where $A_{\text{S}}$ is a system operator, while coupletox=false gives the number-converving coupling $H_{\text{I}}=(A_{\text{S}} b_{0}^\dagger + A_{\text{S}}^\dagger b_0)$tree: If true the resulting chain will be of type TreeNetwork; useful for construcing tree-MPOs
Example
One can constuct a system site tensor to couple to a chain by using the function up to populate the tensor. For example, to construct a system site with Hamiltonian Hs and coupling operator As, the system tensor M is constructed as follows for a non-number conserving interaction:
u = one(Hs) # system identity
+0.0
sourceMPSDynamics.entanglemententropy — Methodentanglemententropy(A)
For a list of tensors A representing a right orthonormalized MPS, compute the entanglement entropy for a bipartite cut for every bond.
sourceMPSDynamics.findchainlength — Methodfindchainlength(T, cparams...; eps=10^-6)
Estimate length of chain required for a particular set of chain parameters by calculating how long an excitation on the first site takes to reach the end. The chain length is given as the length required for the excitation to have just reached the last site after time T. The initial number of sites in cparams has to be larger than the findchainlength result.
sourceMPSDynamics.findchild — Methodfindchild(node::TreeNode, id::Int)
Return integer corresponding to the which number child site id is of node.
sourceMPSDynamics.hbathchain — Methodhbathchain(N::Int, d::Int, chainparams, longrangecc...; tree=false, reverse=false, coupletox=false)
Generate MPO representing a tight-binding chain of N oscillators with d Fock states each. Chain parameters are supplied in the standard form: chainparams $=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$. The output does not itself represent a complete MPO but will possess an end which is open and should be attached to another tensor site, usually representing the system.
Arguments
reverse: If reverse=true create a chain were the last (i.e. Nth) site is the site which couples to the systemcoupletox: Used to choose the form of the system coupling. coupletox=true gives a non-number conserving coupling of the form $H_{\text{I}}= A_{\text{S}}(b_{0}^\dagger + b_0)$ where $A_{\text{S}}$ is a system operator, while coupletox=false gives the number-converving coupling $H_{\text{I}}=(A_{\text{S}} b_{0}^\dagger + A_{\text{S}}^\dagger b_0)$tree: If true the resulting chain will be of type TreeNetwork; useful for construcing tree-MPOs
Example
One can constuct a system site tensor to couple to a chain by using the function up to populate the tensor. For example, to construct a system site with Hamiltonian Hs and coupling operator As, the system tensor M is constructed as follows for a non-number conserving interaction:
u = one(Hs) # system identity
M = zeros(1,3,2,2)
M[1, :, :, :] = up(Hs, As, u)
The full MPO can then be constructed with:
Hmpo = [M, hbathchain(N, d, chainparams, coupletox=true)...]
Similarly for a number conserving interaction the site tensor would look like:
u = one(Hs) # system identity
M = zeros(1,4,2,2)
-M[1, :, :, :] = up(Hs, As, As', u)
And the full MPO would be
Hmpo = [M, hbathchain(N, d, chainparams; coupletox=false)...]
sourceMPSDynamics.heisenbergmpo — Functionheisenbergmpo(N::Int, J=1.0) = xyzmpo(N; Jx=J)
Generate MPO for the N-spin Heisenberg XXX model, defined by the Hamiltonian
$H = \sum_{n=1}^{N-1} -J σ_x^{n} σ_x^{n+1} - J σ_y^{n} σ_y^{n+1} - J σ_z^{n} σ_z^{n+1}$
with $σ_x^{n}, σ_y^{n}, σ_z^{n}$ the Pauli spin-1/2 matrices of the $n^\text{th}$ site.
sourceMPSDynamics.ibmmpo — Methodibmmpo(ω0, d, N, chainparams; tree=false)
Generate MPO for a spin-1/2 coupled to a chain of harmonic oscillators with the interacting boson model (IBM), defined by the Hamiltonian
$H = \frac{ω_0}{2}σ_z + c_0σ_z(b_0^\dagger+b_0) + \sum_{i=0}^{N-1} t_i (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N} ϵ_ib_i^\dagger b_i$.
The spin is on site 1 of the MPS and the bath modes are to the right.
This Hamiltonain is unitarily equivalent (before the truncation to N sites) to the spin-boson Hamiltonian defined by
$H = \frac{ω_0}{2}σ_z + σ_z\int_0^∞ dω\sqrt{\frac{J(ω)}{π}}(b_ω^\dagger+b_ω) + \int_0^∞ dω ωb_ω^\dagger b_ω$.
The chain parameters, supplied by chainparams=$[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$, can be chosen to represent any arbitrary spectral density $J(ω)$ at any temperature.
sourceMPSDynamics.isingmpo — Methodisingmpo(N; J=1.0, h=1.0)
Generate MPO for the N-spin 1D Ising model with external field $\vec{h} = (0,0,h)$, defined by the Hamiltonian
$H = \sum_{n=1}^{N-1} -J_x σ_x^{n} σ_x^{n+1} + \sum_{n=1}^{N}(- h_z σ_z^{n})$
with $σ_x^{n}, σ_y^{n}, σ_z^{n}$ the Pauli spin-1/2 matrices of the $n^\text{th}$ site.
sourceMPSDynamics.longrange_isingmpo — Functionlongrange_isingmpo(N::Int, α::Float64=0.; J=1.0, h=1.0) = longrange_xyzmpo(N, α; Jx=J, Jy=0., Jz=0., hz=h, hx=0.)
sourceMPSDynamics.longrange_xyzmpo — Functionlongrange_xyzmpo(N::Int, α::Float64=0.; Jx=1.0, Jy=Jx, Jz=Jx, hx=0., hz=0.)
Gennerate MPO for the N-spin long-range XYZ model with external field $\vec{h}=(h_x, 0, h_z)$, , defined by the Hamiltonian
sourceMPSDynamics.measure — Methodmeasure(A, O; kwargs...)
Measure observable O on mps state A
sourceMPSDynamics.measure1siteoperator — Methodmeasure1siteoperator(A::Vector, O, sites::Vector{Int})
For a list of tensors A representing a right orthonormalized MPS, compute the local expectation value of a one-site operator O for every site or just one if it is specified.
For calculating operators on single sites this will be more efficient if the site is on the left of the mps.
sourceMPSDynamics.measure1siteoperator — Methodmeasure1siteoperator(A::Vector, O)
For a list of tensors A representing a right orthonormalized MPS, compute the local expectation value of a one-site operator O for every site.
sourceMPSDynamics.measure2siteoperator — Method measure2siteoperator(A::Vector, M1, M2, j1, j2)
Caculate expectation of M1*M2 where M1 acts on site j1 and M2 acts on site j2, assumes A is right normalised.
sourceMPSDynamics.modemps — Functionmodemps(N::Int, k::Vector{Int}, numex::Int, chainparams=[fill(1.0,N), fill(1.0,N-1)])
Generate an MPS with numex excitations of an equal superposition of modes k of a bosonic tight-binding chain.
sourceMPSDynamics.modemps — Functionmodemps(N::Int, k::Int, numex::Int, chainparams=[fill(1.0,N), fill(1.0,N-1)])
Generate an MPS with numex excitations of mode k of a bosonic tight-binding chain.
chainparams takes the form [e::Vector, t::Vector] where e are the on-site energies and t are the hoppping parameters.
The returned MPS will have bond-dimensions and physical dimensions numex+1
sourceMPSDynamics.mpsembed! — Methodmpsembed(A::Vector, Dmax::Int)
Embed MPS A in manifold of max bond-dimension Dmax
sourceMPSDynamics.mpsleftnorm! — Functionmpsleftnorm!(A::Vector, jq::Int=length(A))
Left orthoganalise MPS A up to site jq.
sourceMPSDynamics.mpsmixednorm! — Methodmpsmixednorm!(A::TreeNetwork, id::Int)
Normalise tree-MPS A such that orthogonality centre is on site id.
sourceMPSDynamics.mpsmixednorm! — Methodmpsmixednorm!(A::Vector, OC::Int)
Put MPS A into mixed canonical form with orthogonality centre on site OC.
sourceMPSDynamics.mpsmoveoc! — Methodmpsmoveoc!(A::TreeNetwork, id::Int)
Move the orthogonality centre of right normalised tree-MPS A to site id.
This function will be more efficient than using mpsmixednorm! if the tree-MPS is already right-normalised.
sourceMPSDynamics.mpsrightnorm! — Functionmpsrightnorm!(A::Vector, jq::Int=1)
Right orthoganalise MPS A up to site jq.
sourceMPSDynamics.mpsrightnorm! — Methodmpsrightnorm!(A::TreeNetwork)
When applied to a tree-MPS, right normalise towards head-node.
sourceMPSDynamics.mpsshiftoc! — Methodmpsshiftoc!(A::TreeNetwork, newhd::Int)
Shift the orthogonality centre by one site, setting new head-node newhd.
sourceMPSDynamics.multiply — Methodmultiply(M1::Vector, M2::Vector)
Calculates M1*M2 where M1 and M2 are MPOs
sourceMPSDynamics.nearestneighbourmpo — Functionnearestneighbourmpo(N::Int, h0, A, Ad = A')
sourceMPSDynamics.nearestneighbourmpo — Functionnearestneighbourmpo(tree_::Tree, h0, A, Ad = A')
sourceMPSDynamics.normmps — Methodnormmps(net::TreeNetwork; mpsorthog=:None)
When applied to a tree-MPS mpsorthog=:Left is not defined.
sourceMPSDynamics.normmps — Methodnormmps(A::Vector; mpsorthog=:None)
Calculate norm of MPS A.
Setting mpsorthog=:Right/:Left will calculate the norm assuming right/left canonical form. Setting mpsorthog=OC::Int will cause the norm to be calculated assuming the orthoganility center is on site OC. If mpsorthog is :None the norm will be calculated as an MPS-MPS product.
sourceMPSDynamics.orthcentersmps — Methodorthcentersmps(A)
Compute the orthoganality centres of MPS A.
Return value is a list in which each element is the corresponding site tensor of A with the orthoganility centre on that site. Assumes A is right normalised.
sourceMPSDynamics.physdims — Methodphysdims(M)
Return the physical dimensions of an MPS or MPO M.
sourceMPSDynamics.productstatemps — Functionproductstatemps(physdims::Dims, Dmax=1; state=:Vacuum, mpsorthog=:Right)
Return an MPS representing a product state with local Hilbert space dimensions given by physdims.
By default all bond-dimensions will be 1 since the state is a product state. However, to embed the product state in a manifold of greater bond-dimension, Dmax can be set accordingly.
The indvidual states of the MPS sites can be provided by setting state to a list of column vectors. Setting state=:Vacuum will produce an MPS in the vacuum state (where the state of each site is represented by a column vector with a 1 in the first row and zeros elsewhere). Setting state=:FullOccupy will produce an MPS in which each site is fully occupied (ie. a column vector with a 1 in the last row and zeros elsewhere).
The argument mpsorthog can be used to set the gauge of the resulting MPS.
Example
julia> ψ = unitcol(1,2); d = 6; N = 30; α = 0.1; Δ = 0.0; ω0 = 0.2; s = 1
+M[1, :, :, :] = up(Hs, As, As', u)
And the full MPO would be
Hmpo = [M, hbathchain(N, d, chainparams; coupletox=false)...]
sourceMPSDynamics.heisenbergmpo — Functionheisenbergmpo(N::Int, J=1.0) = xyzmpo(N; Jx=J)
Generate MPO for the N-spin Heisenberg XXX model, defined by the Hamiltonian
$H = \sum_{n=1}^{N-1} -J σ_x^{n} σ_x^{n+1} - J σ_y^{n} σ_y^{n+1} - J σ_z^{n} σ_z^{n+1}$
with $σ_x^{n}, σ_y^{n}, σ_z^{n}$ the Pauli spin-1/2 matrices of the $n^\text{th}$ site.
sourceMPSDynamics.ibmmpo — Methodibmmpo(ω0, d, N, chainparams; tree=false)
Generate MPO for a spin-1/2 coupled to a chain of harmonic oscillators with the interacting boson model (IBM), defined by the Hamiltonian
$H = \frac{ω_0}{2}σ_z + c_0σ_z(b_0^\dagger+b_0) + \sum_{i=0}^{N-1} t_i (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N} ϵ_ib_i^\dagger b_i$.
The spin is on site 1 of the MPS and the bath modes are to the right.
This Hamiltonain is unitarily equivalent (before the truncation to N sites) to the spin-boson Hamiltonian defined by
$H = \frac{ω_0}{2}σ_z + σ_z\int_0^∞ dω\sqrt{\frac{J(ω)}{π}}(b_ω^\dagger+b_ω) + \int_0^∞ dω ωb_ω^\dagger b_ω$.
The chain parameters, supplied by chainparams=$[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$, can be chosen to represent any arbitrary spectral density $J(ω)$ at any temperature.
sourceMPSDynamics.isingmpo — Methodisingmpo(N; J=1.0, h=1.0)
Generate MPO for the N-spin 1D Ising model with external field $\vec{h} = (0,0,h)$, defined by the Hamiltonian
$H = \sum_{n=1}^{N-1} -J_x σ_x^{n} σ_x^{n+1} + \sum_{n=1}^{N}(- h_z σ_z^{n})$
with $σ_x^{n}, σ_y^{n}, σ_z^{n}$ the Pauli spin-1/2 matrices of the $n^\text{th}$ site.
sourceMPSDynamics.longrange_isingmpo — Functionlongrange_isingmpo(N::Int, α::Float64=0.; J=1.0, h=1.0) = longrange_xyzmpo(N, α; Jx=J, Jy=0., Jz=0., hz=h, hx=0.)
sourceMPSDynamics.longrange_xyzmpo — Functionlongrange_xyzmpo(N::Int, α::Float64=0.; Jx=1.0, Jy=Jx, Jz=Jx, hx=0., hz=0.)
Gennerate MPO for the N-spin long-range XYZ model with external field $\vec{h}=(h_x, 0, h_z)$, , defined by the Hamiltonian
sourceMPSDynamics.measure — Methodmeasure(A, O; kwargs...)
Measure observable O on mps state A
sourceMPSDynamics.measure1siteoperator — Methodmeasure1siteoperator(A::Vector, O, sites::Vector{Int})
For a list of tensors A representing a right orthonormalized MPS, compute the local expectation value of a one-site operator O for every site or just one if it is specified.
For calculating operators on single sites this will be more efficient if the site is on the left of the mps.
sourceMPSDynamics.measure1siteoperator — Methodmeasure1siteoperator(A::Vector, O)
For a list of tensors A representing a right orthonormalized MPS, compute the local expectation value of a one-site operator O for every site.
sourceMPSDynamics.measure2siteoperator — Method measure2siteoperator(A::Vector, M1, M2, j1, j2)
Caculate expectation of M1*M2 where M1 acts on site j1 and M2 acts on site j2, assumes A is right normalised.
sourceMPSDynamics.modemps — Functionmodemps(N::Int, k::Vector{Int}, numex::Int, chainparams=[fill(1.0,N), fill(1.0,N-1)])
Generate an MPS with numex excitations of an equal superposition of modes k of a bosonic tight-binding chain.
sourceMPSDynamics.modemps — Functionmodemps(N::Int, k::Int, numex::Int, chainparams=[fill(1.0,N), fill(1.0,N-1)])
Generate an MPS with numex excitations of mode k of a bosonic tight-binding chain.
chainparams takes the form [e::Vector, t::Vector] where e are the on-site energies and t are the hoppping parameters.
The returned MPS will have bond-dimensions and physical dimensions numex+1
sourceMPSDynamics.mpsembed! — Methodmpsembed(A::Vector, Dmax::Int)
Embed MPS A in manifold of max bond-dimension Dmax
sourceMPSDynamics.mpsleftnorm! — Functionmpsleftnorm!(A::Vector, jq::Int=length(A))
Left orthoganalise MPS A up to site jq.
sourceMPSDynamics.mpsmixednorm! — Methodmpsmixednorm!(A::TreeNetwork, id::Int)
Normalise tree-MPS A such that orthogonality centre is on site id.
sourceMPSDynamics.mpsmixednorm! — Methodmpsmixednorm!(A::Vector, OC::Int)
Put MPS A into mixed canonical form with orthogonality centre on site OC.
sourceMPSDynamics.mpsmoveoc! — Methodmpsmoveoc!(A::TreeNetwork, id::Int)
Move the orthogonality centre of right normalised tree-MPS A to site id.
This function will be more efficient than using mpsmixednorm! if the tree-MPS is already right-normalised.
sourceMPSDynamics.mpsrightnorm! — Functionmpsrightnorm!(A::Vector, jq::Int=1)
Right orthoganalise MPS A up to site jq.
sourceMPSDynamics.mpsrightnorm! — Methodmpsrightnorm!(A::TreeNetwork)
When applied to a tree-MPS, right normalise towards head-node.
sourceMPSDynamics.mpsshiftoc! — Methodmpsshiftoc!(A::TreeNetwork, newhd::Int)
Shift the orthogonality centre by one site, setting new head-node newhd.
sourceMPSDynamics.multiply — Methodmultiply(M1::Vector, M2::Vector)
Calculates M1*M2 where M1 and M2 are MPOs
sourceMPSDynamics.nearestneighbourmpo — Functionnearestneighbourmpo(N::Int, h0, A, Ad = A')
sourceMPSDynamics.nearestneighbourmpo — Functionnearestneighbourmpo(tree_::Tree, h0, A, Ad = A')
sourceMPSDynamics.normmps — Methodnormmps(net::TreeNetwork; mpsorthog=:None)
When applied to a tree-MPS mpsorthog=:Left is not defined.
sourceMPSDynamics.normmps — Methodnormmps(A::Vector; mpsorthog=:None)
Calculate norm of MPS A.
Setting mpsorthog=:Right/:Left will calculate the norm assuming right/left canonical form. Setting mpsorthog=OC::Int will cause the norm to be calculated assuming the orthoganility center is on site OC. If mpsorthog is :None the norm will be calculated as an MPS-MPS product.
sourceMPSDynamics.orthcentersmps — Methodorthcentersmps(A)
Compute the orthoganality centres of MPS A.
Return value is a list in which each element is the corresponding site tensor of A with the orthoganility centre on that site. Assumes A is right normalised.
sourceMPSDynamics.physdims — Methodphysdims(M)
Return the physical dimensions of an MPS or MPO M.
sourceMPSDynamics.productstatemps — Functionproductstatemps(physdims::Dims, Dmax=1; state=:Vacuum, mpsorthog=:Right)
Return an MPS representing a product state with local Hilbert space dimensions given by physdims.
By default all bond-dimensions will be 1 since the state is a product state. However, to embed the product state in a manifold of greater bond-dimension, Dmax can be set accordingly.
The indvidual states of the MPS sites can be provided by setting state to a list of column vectors. Setting state=:Vacuum will produce an MPS in the vacuum state (where the state of each site is represented by a column vector with a 1 in the first row and zeros elsewhere). Setting state=:FullOccupy will produce an MPS in which each site is fully occupied (ie. a column vector with a 1 in the last row and zeros elsewhere).
The argument mpsorthog can be used to set the gauge of the resulting MPS.
Example
julia> ψ = unitcol(1,2); d = 6; N = 30; α = 0.1; Δ = 0.0; ω0 = 0.2; s = 1
julia> cpars = chaincoeffs_ohmic(N, α, s)
julia> H = spinbosonmpo(ω0, Δ, d, N, cpars)
-julia> A = productstatemps(physdims(H), state=[ψ, fill(unitcol(1,d), N)...]) # MPS representation of |ψ>|Vacuum>
sourceMPSDynamics.productstatemps — Functionproductstatemps(N::Int, d::Int, Dmax=1; state=:Vacuum, mpsorthog=:Right)
Return an N-site MPS with all local Hilbert space dimensions given by d.
sourceMPSDynamics.randisometry — Methodrandisometry([T=Float64], dims...)
Construct a random isometry
sourceMPSDynamics.randmps — Functionrandmps(N::Int, d::Int, Dmax::Int, T=Float64)
Construct a random, N-site, right-normalised MPS with all local Hilbert space dimensions given by d.
sourceMPSDynamics.randmps — Functionrandmps(tree::Tree, physdims, Dmax::Int, T::Type{<:Number} = Float64)
Construct a random, right-normalised, tree-MPS, with structure given by tree and max bond-dimension given by Dmax.
The local Hilbert space dimensions are specified by physdims which can either be of type Dims{length(tree)}, specifying the dimension of each site, or of type Int, in which case the same local dimension is used for every site.
sourceMPSDynamics.randmps — Methodrandmps(physdims::Dims{N}, Dmax::Int, T::Type{<:Number} = Float64) where {N}
Construct a random, right-normalised MPS with local Hilbert space dimensions given by physdims and max bond-dimension given by Dmax.
T specifies the element type, eg. use T=ComplexF64 for a complex valued MPS.
sourceMPSDynamics.randtree — Methodrandtree(numnodes::Int, maxdegree::Int)
Construct a random tree with nummodes modes and max degree maxdegree.
sourceMPSDynamics.readchaincoeffs — Methodreadchaincoeffs(fdir, params...)
sourceMPSDynamics.rmsd — Methodrmsd(dat1::Vector{Float64}, dat2::Vector{Float64})
Calculate the root mean squared difference between two measurements of an observable over the same time period.
sourceMPSDynamics.spinbosonmpo — Methodspinbosonmpo(ω0, Δ, d, N, chainparams; rwa=false, tree=false)
Generate MPO for a spin-1/2 coupled to a chain of harmonic oscillators, defined by the Hamiltonian
$H = \frac{ω_0}{2}σ_z + Δσ_x + c_0σ_x(b_0^\dagger+b_0) + \sum_{i=0}^{N-1} t_i (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N} ϵ_ib_i^\dagger b_i$.
The spin is on site 1 of the MPS and the bath modes are to the right.
This Hamiltonain is unitarily equivalent (before the truncation to N sites) to the spin-boson Hamiltonian defined by
$H = \frac{ω_0}{2}σ_z + Δσ_x + σ_x\int_0^∞ dω\sqrt{\frac{J(ω)}{π}}(b_ω^\dagger+b_ω) + \int_0^∞ dω ωb_ω^\dagger b_ω$.
The chain parameters, supplied by chainparams=$[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$, can be chosen to represent any arbitrary spectral density $J(ω)$ at any temperature.
The rotating wave approximation can be made by setting rwa=true.
sourceMPSDynamics.spinchainmpo — Methodspinchainmpo(N::Int; J=1.0, hz=1.0, hx=0.0, i=div(N,2))
sourceMPSDynamics.svdmps — Methodsvdmps(A)
For a right normalised mps A compute the full svd spectrum for a bipartition at every bond.
sourceMPSDynamics.svdtrunc — MethodU, S, Vd = svdtrunc(A; truncdim = max(size(A)...), truncerr = 0.)
Perform a truncated SVD, with maximum number of singular values to keep equal to truncdim or truncating any singular values smaller than truncerr. If both options are provided, the smallest number of singular values will be kept. Unlike the SVD in Julia, this returns matrix U, a diagonal matrix (not a vector) S, and Vt such that A ≈ U * S * Vt
sourceMPSDynamics.tightbinding_mpo — Methodtightbinding_mpo(N, ϵd, chainparams1, chainparams2)
+julia> A = productstatemps(physdims(H), state=[ψ, fill(unitcol(1,d), N)...]) # MPS representation of |ψ>|Vacuum>
sourceMPSDynamics.productstatemps — Functionproductstatemps(N::Int, d::Int, Dmax=1; state=:Vacuum, mpsorthog=:Right)
Return an N-site MPS with all local Hilbert space dimensions given by d.
sourceMPSDynamics.randisometry — Methodrandisometry([T=Float64], dims...)
Construct a random isometry
sourceMPSDynamics.randmps — Functionrandmps(N::Int, d::Int, Dmax::Int, T=Float64)
Construct a random, N-site, right-normalised MPS with all local Hilbert space dimensions given by d.
sourceMPSDynamics.randmps — Functionrandmps(tree::Tree, physdims, Dmax::Int, T::Type{<:Number} = Float64)
Construct a random, right-normalised, tree-MPS, with structure given by tree and max bond-dimension given by Dmax.
The local Hilbert space dimensions are specified by physdims which can either be of type Dims{length(tree)}, specifying the dimension of each site, or of type Int, in which case the same local dimension is used for every site.
sourceMPSDynamics.randmps — Methodrandmps(physdims::Dims{N}, Dmax::Int, T::Type{<:Number} = Float64) where {N}
Construct a random, right-normalised MPS with local Hilbert space dimensions given by physdims and max bond-dimension given by Dmax.
T specifies the element type, eg. use T=ComplexF64 for a complex valued MPS.
sourceMPSDynamics.randtree — Methodrandtree(numnodes::Int, maxdegree::Int)
Construct a random tree with nummodes modes and max degree maxdegree.
sourceMPSDynamics.readchaincoeffs — Methodreadchaincoeffs(fdir, params...)
sourceMPSDynamics.rmsd — Methodrmsd(dat1::Vector{Float64}, dat2::Vector{Float64})
Calculate the root mean squared difference between two measurements of an observable over the same time period.
sourceMPSDynamics.spinbosonmpo — Methodspinbosonmpo(ω0, Δ, d, N, chainparams; rwa=false, tree=false)
Generate MPO for a spin-1/2 coupled to a chain of harmonic oscillators, defined by the Hamiltonian
$H = \frac{ω_0}{2}σ_z + Δσ_x + c_0σ_x(b_0^\dagger+b_0) + \sum_{i=0}^{N-1} t_i (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N} ϵ_ib_i^\dagger b_i$.
The spin is on site 1 of the MPS and the bath modes are to the right.
This Hamiltonain is unitarily equivalent (before the truncation to N sites) to the spin-boson Hamiltonian defined by
$H = \frac{ω_0}{2}σ_z + Δσ_x + σ_x\int_0^∞ dω\sqrt{\frac{J(ω)}{π}}(b_ω^\dagger+b_ω) + \int_0^∞ dω ωb_ω^\dagger b_ω$.
The chain parameters, supplied by chainparams=$[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$, can be chosen to represent any arbitrary spectral density $J(ω)$ at any temperature.
The rotating wave approximation can be made by setting rwa=true.
sourceMPSDynamics.spinchainmpo — Methodspinchainmpo(N::Int; J=1.0, hz=1.0, hx=0.0, i=div(N,2))
sourceMPSDynamics.svdmps — Methodsvdmps(A)
For a right normalised mps A compute the full svd spectrum for a bipartition at every bond.
sourceMPSDynamics.svdtrunc — MethodU, S, Vd = svdtrunc(A; truncdim = max(size(A)...), truncerr = 0.)
Perform a truncated SVD, with maximum number of singular values to keep equal to truncdim or truncating any singular values smaller than truncerr. If both options are provided, the smallest number of singular values will be kept. Unlike the SVD in Julia, this returns matrix U, a diagonal matrix (not a vector) S, and Vt such that A ≈ U * S * Vt
sourceMPSDynamics.tightbinding_mpo — Methodtightbinding_mpo(N, ϵd, chainparams1, chainparams2)
Generate MPO for a tight-binding chain of N fermionic sites with a single impurity site (fermionic as well)
of energy ϵd at the center. The impurity is coupled to two leads, each described by a set of chain parameters.
@@ -36,4 +36,4 @@
* chainparams1::Array{Real,1}: chain parameters for the first lead
* chainparams2::Array{Real,1}: chain parameters for the second lead
-The chain parameters are given in the standard form: `chainparams` ``=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]``.
sourceMPSDynamics.tightbindingmpo — Methodtightbindingmpo(N::Int, d::Int; J=1.0, e=1.0)
sourceMPSDynamics.tunnelingmpo — Methodtunnelingmpo(ϵ, delta, α, s, β, d::Int, nummodes::Int; tree=false, ωc=1)
sourceMPSDynamics.twobathspinmpo — Functiontwobathspinmpo(ω0, Δ, Nl, Nr, dl, dr, chainparamsl=[fill(1.0,N),fill(1.0,N-1), 1.0], chainparamsr=chainparamsl; tree=false)
Generate MPO for a spin-1/2 coupled to two chains of harmonic oscillators, defined by the Hamiltonian
$H = \frac{ω_0}{2}σ_z + Δσ_x + c_0^rσ_x(b_0^\dagger+b_0) + \sum_{i=0}^{N_r-1} t_i^r (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N_r} ϵ_i^rb_i^\dagger b_i + c_0^lσ_x(d_0^\dagger+d_0) + \sum_{i=0}^{N_l-1} t_i^l (d_{i+1}^\dagger d_i +h.c.) + \sum_{i=0}^{N_l} ϵ_i^l d_i^\dagger d_i$.
The spin is on site $N_l + 1$ of the MPS, surrounded by the left chain modes and the right chain modes.
This Hamiltonain is unitarily equivalent (before the truncation to N sites) to the spin-boson Hamiltonian defined by
$H = \frac{ω_0}{2}σ_z + Δσ_x + σ_x\int_0^∞ dω\sqrt{\frac{J(ω)}{π}}(b_ω^\dagger+b_ω) + \int_0^∞ dω ωb_ω^\dagger b_ωi + σ_x\int_0^∞ dω\sqrt{\frac{J^l(ω)}{π}}(d_ω^\dagger+d_ω) + \int_0^∞ dω ωd_ω^\dagger d_ω$.
The chain parameters, supplied by chainparams=$[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$, can be chosen to represent any arbitrary spectral density $J(ω)$ at any temperature. The two chains can have a different spectral density.
sourceMPSDynamics.xxzmpo — Functionxxzmpo(N::Int, Δ = 1.0, J=1.0) = xyzmpo(N; Jx=J, Jy=J, Jz=J*Δ)
Generate MPO for the N-spin XXZ model, defined by the Hamiltonian
$H = \sum_{n=1}^{N-1} -J σ_x^{n} σ_x^{n+1} - J σ_y^{n} σ_y^{n+1} - \Delta J σ_z^{n} σ_z^{n+1}$
with $σ_x^{n}, σ_y^{n}, σ_z^{n}$ the Pauli spin-1/2 matrices of the $n^\text{th}$ site.
sourceMPSDynamics.xyzmpo — Methodxyzmpo(N::Int; Jx=1.0, Jy=Jx, Jz=Jx, hx=0., hz=0.)
Generate MPO for the N-spin XYZ model with external field $\vec{h}=(h_x, 0, h_z)$, , defined by the Hamiltonian
$H = \sum_{n=1}^{N-1} -J_x σ_x^{n} σ_x^{n+1} - J_y σ_y^{n} σ_y^{n+1} - J_z σ_z^{n} σ_z^{n+1} + \sum_{n=1}^{N}(- h_x σ_x^{n} - h_z σ_z^{n})$
with $σ_x^{n}, σ_y^{n}, σ_z^{n}$ the Pauli spin-1/2 matrices of the $n^\text{th}$ site.
sourceSettings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
MPSDynamics.tightbindingmpo — Methodtightbindingmpo(N::Int, d::Int; J=1.0, e=1.0)MPSDynamics.tunnelingmpo — Methodtunnelingmpo(ϵ, delta, α, s, β, d::Int, nummodes::Int; tree=false, ωc=1)MPSDynamics.twobathspinmpo — Functiontwobathspinmpo(ω0, Δ, Nl, Nr, dl, dr, chainparamsl=[fill(1.0,N),fill(1.0,N-1), 1.0], chainparamsr=chainparamsl; tree=false)Generate MPO for a spin-1/2 coupled to two chains of harmonic oscillators, defined by the Hamiltonian
$H = \frac{ω_0}{2}σ_z + Δσ_x + c_0^rσ_x(b_0^\dagger+b_0) + \sum_{i=0}^{N_r-1} t_i^r (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N_r} ϵ_i^rb_i^\dagger b_i + c_0^lσ_x(d_0^\dagger+d_0) + \sum_{i=0}^{N_l-1} t_i^l (d_{i+1}^\dagger d_i +h.c.) + \sum_{i=0}^{N_l} ϵ_i^l d_i^\dagger d_i$.
The spin is on site $N_l + 1$ of the MPS, surrounded by the left chain modes and the right chain modes.
This Hamiltonain is unitarily equivalent (before the truncation to N sites) to the spin-boson Hamiltonian defined by
$H = \frac{ω_0}{2}σ_z + Δσ_x + σ_x\int_0^∞ dω\sqrt{\frac{J(ω)}{π}}(b_ω^\dagger+b_ω) + \int_0^∞ dω ωb_ω^\dagger b_ωi + σ_x\int_0^∞ dω\sqrt{\frac{J^l(ω)}{π}}(d_ω^\dagger+d_ω) + \int_0^∞ dω ωd_ω^\dagger d_ω$.
The chain parameters, supplied by chainparams=$[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]$, can be chosen to represent any arbitrary spectral density $J(ω)$ at any temperature. The two chains can have a different spectral density.
MPSDynamics.xxzmpo — Functionxxzmpo(N::Int, Δ = 1.0, J=1.0) = xyzmpo(N; Jx=J, Jy=J, Jz=J*Δ)Generate MPO for the N-spin XXZ model, defined by the Hamiltonian
$H = \sum_{n=1}^{N-1} -J σ_x^{n} σ_x^{n+1} - J σ_y^{n} σ_y^{n+1} - \Delta J σ_z^{n} σ_z^{n+1}$
with $σ_x^{n}, σ_y^{n}, σ_z^{n}$ the Pauli spin-1/2 matrices of the $n^\text{th}$ site.
MPSDynamics.xyzmpo — Methodxyzmpo(N::Int; Jx=1.0, Jy=Jx, Jz=Jx, hx=0., hz=0.)Generate MPO for the N-spin XYZ model with external field $\vec{h}=(h_x, 0, h_z)$, , defined by the Hamiltonian
$H = \sum_{n=1}^{N-1} -J_x σ_x^{n} σ_x^{n+1} - J_y σ_y^{n} σ_y^{n+1} - J_z σ_z^{n} σ_z^{n+1} + \sum_{n=1}^{N}(- h_x σ_x^{n} - h_z σ_z^{n})$
with $σ_x^{n}, σ_y^{n}, σ_z^{n}$ the Pauli spin-1/2 matrices of the $n^\text{th}$ site.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
We consider, in the Schrödinger picture, a general Hamiltonian where a non-specified system interacts linearly with a bosonic environments
\[\begin{aligned} +
We consider, in the Schrödinger picture, a general Hamiltonian where a non-specified system interacts linearly with a bosonic environments
\[\begin{aligned} \hat{H} =& \hat{H}_S + \int_0^{+\infty} \hbar\omega\hat{a}^\dagger_\omega\hat{a}_\omega\ \mathrm{d}\omega + \hat{A}_S\int_0^{+\infty}\sqrt{J(\omega)}\left(\hat{a}_\omega + \hat{a}^\dagger_\omega\right)\mathrm{d}\omega \end{aligned}\]
where $\hat a_\omega$ ($\hat a^\dagger_\omega$) is a bosonic annihilation (creation) operator for a normal mode of the environment of energy $\hbar\omega$, $\hat{A}_S$ is a system operator, and $J(\omega) = \sum_k |g_k|^2\delta(\omega - \omega_k)$ is the bath spectral density (SD), defined with the microscopic system-environment coupling strength $g_k$. The SD quantifies the coupling strengths of the different normal modes of the environment with the system. Any SD that is not flat corresponds to a non-Markovian environment.
Let us consider the Hamiltonian presented in Eq.(1). To investigate its rich phenomenology, using the powerful tools provided by tensor network techniques, it is convenient to map the system and its environment from a star-like to a chain-like configuration, with local interactions. To do so, we introduce a unitary transformation of the continuous normal modes $\hat{a}_\omega$ to an infinite discrete set of interacting modes $\hat{b}_n$[chin_exact_2010].
\[ \hat{a}_\omega = \sum_{n=0}^{+\infty} U_n(\omega)\hat{b}_n = \sum_{n=0}^{+\infty} \sqrt{J(\omega)}P_n(\omega)\hat{b}_n\ ,\]
where $P_n(\omega)$ are orthonormal polynomials such that
\[ \int_{0}^{+\infty}P_n(\omega)P_m(\omega)J(\omega)\mathrm{d}\omega = \delta_{n,m}\ ;\]
and the inverse transformation is
\[ \hat{b}_n = \int_0^{+\infty} U_n(\omega)\hat{a}_\omega\mathrm{d}\omega\ .\]
Note that the orthonormality of the polynomials ensures the unitarity of the transformation defined in Eq.(2). The mapping from a continuous set of modes to a (still infinite) discrete set might seem counter-intuitive, however it is a direct consequence of the separability of the underlying Hilbert space.
Under this transformation, the Hamiltonian in Eq.(1) becomes
\[ \hat{H}= \hat{H}_S + \sum_{n=0}^{+\infty}\varepsilon_n\hat{b}_n^\dagger\hat{b}_n + t_n(\hat{b}_{n+1}^\dagger\hat{b}_n + \mathrm{h.c.}) + \kappa\hat{A}_S(\hat{b}_0 + \hat{b}_0^\dagger)\ .\]
Hence, this mapping transforms the normal bath Hamiltonian into a tight-binding Hamiltonian with on-site energies $\varepsilon_n$ and hopping energies $t_n$. Another important consequence of this mapping is that now the system only interacts with the first mode $n = 0$ of the chain-mapped environment. The chain coefficients $\varepsilon_n$, $t_n$, and the coupling $\kappa$ depend solely on the SD.
This makes chain mapping a tool of choice for describing systems coupled to environment with highly structured SD (e.g. experimentally measured or calculated ab initio)[chin_role_2013][alvertis_nonequilibrium_2019][dunnett_influence_2021][caycedosoler_exact_2022]. In this new representation, the Hamiltonian in Eq.(5) has naturally a 1D chain topology. This makes its representation as a Matrix Product Operator (MPO) and the representation of the joint {System + Environment} wave-function as a Matrix Product State (MPS) suited [orus_practical_2014][paeckel_timeevolution_2019].
The orthogonal polynomial-based chain mapping and the subsequent representation of the joint wave-function as a MPS (and the operators as MPO) are the building blocks of the Time-dependent Density operator with Orthonormal Polynomials Algorithm (TEDOPA) one of the state-of-the-art numerically exact method to simulate the dynamics of open quantum systems especially in the non-Markovian, non-perturbative regimes both at zero and finite temperatures [prior_efficient_2010][woods_simulating_2015][tamascelli_efficient_2019][dunnett_simulating_2021][lacroix_unveiling_2021].
Assuming a unitary evolution for both the system and environment, the system's dynamics can be isolated by tracing out the environmental degrees of freedom. The density operator for the system at time $t$ is described as
\[\hat{\rho}_S(t) = \text{T}r_E\left\{\hat{U}(t) \hat{\rho}_S(0) \otimes \hat{\rho}_E(0) \hat{U}^\dagger(t)\right\}.\]
Initially, the system state, $\hat{\rho}_S(0)$, can be pure or mixed, and the environment is in a thermal state defined by the inverse temperature $\beta = (k_B T)^{-1}$. This state is represented by a product of Gaussian states
\[\hat{\rho}_E(0) = \bigotimes_\omega \frac{e^{-\beta \omega \hat{a}_\omega^\dagger \hat{a}_\omega}}{Z_\omega(\beta)} .\]
The system's evolution is dictated by the environment's two-time correlation function $S(t)$, which in turn is determined by the spectral density function $J$ and the temperature $\beta$ via the Bose-Einstein distribution $n_\omega(\beta)$
\[S(t) = \int_0^\infty d\omega J(\omega)\left[e^{-i\omega t}(1 + n_\omega(\beta)) + e^{i\omega t} n_\omega(\beta)\right] .\]
To simulate finite temperature effects using a zero-temperature model with the T-TEDOPA method [tamascelli_efficient_2019], we extend the spectral density function to cover both positive and negative frequencies, allowing us to use a pure state description for the environment. This is achieved by defining a new spectral density function $J(\omega, \beta)$ that incorporates the Boltzmann factors, supporting the entire real axis
\[J(\omega, \beta) = \frac{\text{sign}(\omega)J(\left|\omega\right|)}{2} \Big(1 + \coth\Big(\frac{\beta \omega}{2}\Big)\Big).\]
This modified bath has the same correlation function $S(t)$ and thus allows us to maintain a pure state description of the environment, represented as a vacuum state, and avoid the computational complexities of density matrices
\[\left|\text{vac}\right\rangle = \bigotimes_\omega \left|0_\omega\right\rangle .\]
The Hamiltonian of the system interacting with this extended bath now includes temperature-dependent interactions:
\[\hat{H} = \hat{H}_S + \int_{-\infty}^{+\infty} d\omega \omega \hat{a}_\omega^\dagger \hat{a}_\omega + \hat{A}_S \otimes \int_{-\infty}^{+\infty} d\omega \sqrt{J(\omega,\beta)}\left(\hat{a}_\omega^\dagger+\hat{a}_\omega\right),\]
This method simplifies the simulation of finite temperature effects by treating them within an effective zero-temperature framework, thereby keeping the computational advantages of using pure states. In conclusion: the dynamics of the system resulting from the interaction with the original bath, starting in a thermal state at finite temperature, is exactly the same as the one resulting from the interaction with the extended environment, starting in the vacuum state at zero temperature. Once computed the chain coefficients at a given inverse temperature $\beta$, the time evolution of the vacuum state interacting with the extended environment can be efficiently simulated using MPS time evolution methods.
Finite temperature initial states of the system are challenging: in the previous section we have seen that a very efficient solution is to reabsorb the effect of the temperature inside of the spectral density function that characterizes the environment. De Vega and Bañuls present[devega_thermo_2015] another suggestive approach based on the idea of the thermofield transformation. The idea is to double the environmental degrees of freedom and then to apply a Bogoliubov transformation: the real environment in a thermal state is transformed into two virtual environments in the vacuum state, defined as the thermofield vacuum. For any operator of the real environment, the expectation values in the thermal state are equivalent to those calculated in the thermofield vacuum. Behind the thermofield approach there is the concept of purification: the initial mixed state of the environment can be represented as the partial trace of a pure state (the thermofield vacuum) defined on the larger Hilbert space of the modes of the two virtual environments at zero temperature. The interaction of the system and each one of the two environments is then mapped on two separate chains, using the TEDOPA chain mapping to define a unitary transformation that can be applied to the bath modes, in order to define new modes for the environment, mapping it from a star-like to a one dimensional chain-like configuration, which is very well suited for the application of tensor network techniques.
In the thermofield approach, the first step is to introduce the auxiliary environment $E^{\text{aux}}$ of non-interacting bosonic modes of negative frequencies:
\[ \hat H^\text{aux} = \hat H - \sum_k \omega_k \hat c_k^\dagger \hat c_k,\]
where the Hamiltonian $\hat H$ is the bosonic Hamiltonian defined in Eq. \ref{eq:bosonic_ham}. The two environments, of positive and negative frequencies, are assumed to be in a thermal environment at inverse temperature $\beta$; the second step is to apply a thermal Bogoliubov transformation to change the basis. The applied transformation produces two-modes squeezed states:
\[ \hat a_{1k}=e^{-iG} \hat b_k e^{iG}= \cosh(\theta_k) \hat b_k -\sinh(\theta_k) \hat c_k^\dagger \\ \hat a_{2k}=e^{-iG} \hat c_k e^{iG}= \cosh(\theta_k) \hat c_k -\sinh(\theta_k) \hat b_k^\dagger,\]
where the exponent of the squeeze operator is $G = i \sum_k \theta_k(\hat b_k^\dagger \hat c_k^\dagger-\hat c_k \hat b_k)$, and $\theta_k$ is dependent on the temperature as in the following relations, where the number of excitations in the $k$-th mode is $n_k = 1/(e^{\beta \omega_k}-1)$:
\[ \cosh(\theta_k) = \sqrt{1+n_k} = \sqrt{\frac{1}{1-e^{-\beta \omega_k}}} \\ \sinh(\theta_k) =\quad\sqrt{n_k}\quad= \sqrt{\frac{1}{e^{\beta \omega_k}-1}}.\]
The Bogoliubov transformation defines a new squeezed vacuum state, which we write in terms of the vacuum state $\ket{\Omega_0}$ of the operators $\hat b_k$, $\hat c_k$:
\[ |\Omega\rangle = e^{iG} |\Omega_0\rangle, \quad \text{such that: }\quad \hat a_{1k} |\Omega\rangle = 0, \hat a_{2k} |\Omega\rangle = 0.\]
From the vacuum state, we can obtain the thermal state of the original environment:
\[ \hat \rho_E = \text{Tr}_{\text{aux}}\{ |\Omega\rangle\langle\Omega| \},\]
and it can be now used as pure an initial state for both of the environments. Moreover, the expectation value on $\ket{\Omega}$ of the number of physical modes $\hat n_k$ does not vanish:
\[ n_k= \langle\Omega| \hat b_k^\dagger \hat b_k |\Omega\rangle = \sinh^2(\theta_k).\]
Therefore, solving the dynamics given by the original Hamiltonian $\hat H$, starting from the initial condition $\hat \rho_S(0) \otimes \hat \rho_E(\beta)$, is equivalent to solving the dynamics given by the following Hamiltonian:
\[ \hat H = \hat H_S +\hat H_E +\hat H_I = \\ - = \overbrace{\hat A_S}^{\hat H_S} + \overbrace{\sum_k \omega_k \big(\hat a_{1k}^\dagger \hat a_{1k} - \hat a_{2k}^\dagger \hat a_{2k} \big)}^{H_E} + \overbrace{\hat L_S \otimes \sum_k g_{1k}(\hat a_{1k}^\dagger + \hat a_{1k})+\hat L_S \otimes \sum_k g_{2k}(\hat a_{2k}^\dagger + \hat a_{2k})}^{H_I}.\]
where $\hat L_S = \hat L_S^\dagger$, considering $\hat \rho_S(0) \otimes \ket{\Omega}\bra{\Omega}$ as the initial state of system and environment. It is this Hamiltonian that is mapped on two chains with TEDOPA, to be able to perform the time evolution using tensor network techniques.
A useful property of the orthonormal polynomials is that they obey a recurrence relation
\[ P_n(\omega) = (C_{n-1}\omega - A_{n-1})P_{n-1}(\omega) + B_{n-1}P_{n-2}(\omega)\ ,\]
where $A_n$ is related to the first moment of $P_n$, $B_n$ and $C_n$ to the norms of $P_n$ and $P_{n-1}$[appel_mathematics_2007]. This recurrence relation can be used to construct the polynomials with the conditions that $P_0(\omega) = ||p_0||^{-1} = \left(\int_{\mathbb{R}^{+}} J(\omega)\mathrm{d}\omega \right)^{-\frac{1}{2}}$ and $P_{-1}(\omega) = 0$, with $||\bullet||$ the norm of $\bullet$ with respect to the measure $J(\omega)$, and $P_n(\omega) = p_n(\omega)||p_n||^{-1}$ ; where the polynomials $\{p_n\}_{n\in\mathbb{N}}$ are the so called monic polynomials where the factor $a_n$ in front of $\omega^{n}$ is equal to 1.
The energy of the chain mode $n$ is given by $\varepsilon_n = A_n C_n^{-1}$ and $t_n=C_n^{-1}$ is the coupling between mode $n$ and $n+1$[chin_exact_2010].
The system couples only to the first mode with the coupling strength $\kappa = ||p_0||$.
Explain that for some weight function/SD they are known analytically and that for others we can use the build-in routines inspired by Gautschi or the PolyChaos.jl package.
A multipartite quantum state $|\psi\rangle$, e.g. a $N$-site system where the sites can each be in a state $|\phi_i\rangle$ belonging to a $d$-dimensional Hilbert space, can be written as follows
\[ |\psi\rangle = \sum_{\{i_k\}}c_{i_1\ldots i_N}|\phi_{i_1}\rangle\otimes\ldots\otimes|\phi_{i_N}\rangle\ ,\]
where the complex numbers $c_{i_1\ldots i_N}$ are the amplitudes of each state $|\phi_{i_1}\rangle\otimes\ldots\otimes|\phi_{i_N}\rangle$ whose superpositions form in full generality the state $|\psi\rangle$. Thus the state $|\psi\rangle$ can be completely represented by a rank-$N$ tensor $c$ that is the collection of all possible amplitudes $c_{i_1\ldots i_N}$. Here by the rank of a tensor, we simply mean the number of indices it has.
The tensor $c$ of a quantum state $|\psi\rangle$ corresponding to a one-dimensional system can be decomposed into a product of $N$ smaller rank-3 tensors $T_{k}$ (except for the first and last sites where the tensors will have a rank-2)
\[ c_{i_1\ldots i_N} = \sum_{\{\alpha\}} T^{\alpha_1}_{i_1}T^{\alpha_1\alpha_2\ }_{i_2}T^{\alpha_2\alpha_3\ }_{i_3}\ldots T^{\alpha_{N-1}}_{i_N} \ .\]
In this form, the local tensor $T_k$ contains the information on the quantum state on site $k$ and its relation (especially the entanglement) with the neighbouring sites.
The decomposition of the tensor of the amplitudes of a quantum state into a product of smaller rank tensors is called a Matrix Product State decomposition.
The contracted indices $\alpha_k$ between the tensors are called virtual indices and carry information about the correlations between bi-partitions of the state at bond $k$. The number of different values a virtual index can take is called the bond dimension and is denoted $D$. The free indices $i_k$ associated with local quantum states are called physical indices. Thus, they can take $d$ values (with $d$ the dimension of the local Hilbert space).
Any state in the Hilbert space of a one-dimensional many-body system can in principle be represented by a MPS by choosing a sufficiently large value for the bond dimension $D$ [Orus]. On top of this intellectually satisfying property of MPSs being a dense set of states for a 1d-system, they can also be used as a practical Ansätze for a many-body quantum states by setting a maximal allowed value $\chi$ for the bond dimension $D$. In doing so, we restrict ourselves to a corner of the total Hilbert space. The rationale behind this Ansatz is the following: if the initial quantum state of a many-body system has a low bond dimension (typically if the initial state is a product state with $D = 1$), then in a finite time it will only be able to explore a region of the Hilbert space that is not to far away from its starting point. Thus, the bond dimension will not have the time to diverge exponentially [poulin_quantum_2011]. However, depending on the physical system at hand, this sub-manifold of the Hilbert space could still be "too large". There is an additional reason that explains why MPSs are good Ansätze for 1d physical systems. Most many-body Hamiltonians we (physicists) are interested in are local, meaning that the interactions they describe involve objects that are "neighbours". For such Hamiltonians, the ground states (outside of potential critical phases) follow the so called area law for the entanglement entropy [srednicki_entropy_1993][vidal_entanglement_2003][wolf_area_2008]. This law states that the entanglement entropy $S_{vN}$ of a bi-partition of the system is proportional, not to the volume of the partition as one might expect, but to the hyper-surface of the partition's boundary; hence the name "area law". For a 3d system this corresponds to an actual surface area $A$, $S_{vN} \sim A$; for a 2d system it corresponds to the length $L$ of the partition's boundary, $S_{vN} \sim L$; and in 1d the boundary reduces to a point, thus the entropy will be independent of the size of the system $S_{vN} \sim \text{constant}$. The MPSs are states that satisfy this area law.
An application of the Singular Value Decomposition is to create efficient approximations of quantum states to perform computations. The main idea is to reduce the content of the MPS to keep only the parts that contain the physics of interest. One method to realise this approximation is to do a SVD on each of the tensors of the MPS after each time step of the state time-evolution and to trim the smallest singular values in order to decrease the bond dimension of the MPS down to a chosen maximal value $\chi$. The corresponding columns and rows of the unitary matrices $U$ and $V^\dagger$ are also removed. Then, the trimmed matrices $\tilde{U}$, $\tilde{S}$ and $\tilde{V}^\dagger$ are contracted back to give an approximated tensor $T$ with a smaller bond dimension. Another way to apply the restricted rank approximation is to restrict oneself into working in a manifold of fixed bond dimension $D$ and to use methods that can enforce this constraint.
In order to compute expectation values of observables or apply unitary transformations to a quantum state, we need a TN representation of operators. In the same fashion as a one-dimensional quantum state can be represented as a MPS, operators acting on those states can be represented as Matrix Product Operators (MPO). For an operator $\hat{O}$, its MPO can be defined as follows
\[ \hat{O} = \sum_{\{i_k\}\{i_k^{'}\} \{w\}} W^{i_1\ i^{'}_1}_{1\ w_0w_1}\ldots W^{i_N\ i^{'}_N}_{N\ w_{N-1}w_N} |\phi_{i_1^{'}}\ldots \phi_{i_N^{'}}\rangle\langle\phi_{i_1}\ldots \phi_{i_N}| \]
The contracted indices between the tensors are called virtual indices. The free indices are called physical indices and correspond to the different input and output local quantum states. They can take $d$ values (with $d$ the dimension of the local Hilbert space).
A natural extension to the MPS is the (loop-free) tree tensor network. A TTN is a generalisation of the MPS wherein each site, instead of being connected to only one other site to its right, may be connected to any arbitrary number of child sites. Provided the tree does not contain any loops, everything that one can do to an MPS/MPO can be extended straight-forwardly to TTN states and TTN operators. The generalisation to trees introduces no new conceptual complexity (only implementational complexity). The sites of a TTN are usually referred to as nodes. For our purposes, every node of a TTN state and operator has one parent leg, and any number (including zero) of child legs. The first node is known as the head-node and has a dummy parent leg with dimension 1.
The original idea behind TDVP goes back to Dirac [dirac_note_1930] and Frenkel [frenkel_wave_1934]. The main point, in the modern tensor networks formulation, is that instead of solving the Schrödinger equation and then truncating the MPS representation of the quantum state, one can solve the equations of motion projected into a space of restricted bond dimension [haegeman_timedependent_2011][haegeman_unifying_2016].
The general formulation of the Dirac-Frenkel Variational Principle [raab_diracfrenkelmclachlan_2000] is that one looks for a solution $|\varphi\rangle \in \mathcal{M}$ of the Schrödinger equation where $\mathcal{M} \subset \mathcal{H}$ is a manifold of the total Hilbert space $\mathcal{H}$ in which we think that the relevant physical states `live'.
We define $T_{|\varphi\rangle}\mathcal{M}$ the tangent space of $\mathcal{M}$ around the state $|\varphi\rangle$. The criterion to find $|\varphi\rangle$ is that for every state $|\chi\rangle \in T_{|\varphi\rangle}\mathcal{M}$
\[ \langle\chi|\left(\frac{\mathrm{d}}{\mathrm{d}t} - \frac{1}{\mathrm{i}\hbar}\hat{H}\right)|\varphi\rangle =0\ ,\]
which can be interpreted as saying that the time evolution procedure should keep $|\varphi\rangle$ inside of the manifold $\mathcal{M}$.
The term variational in the name of the method comes from the fact that in practice one aims at minimising the right-hand side of the hereinabove Eq. to find $|\varphi\rangle$.
Introducing $\hat{P}_{T_{|\varphi\rangle}\mathcal{M}}$ the projector onto the tangent space $T_{|\varphi\rangle}\mathcal{M}$, we can write the state $|\chi\rangle = \hat{P}_{T_{|\varphi\rangle}\mathcal{M}}|\phi\rangle$ with $|\phi\rangle$ a state in $\mathcal{H}$. Leading to
\[ \forall |\phi\rangle \in \mathcal{H}, \ \langle\phi|\hat{P}_{T_{|\varphi\rangle}\mathcal{M}}\left(\frac{\mathrm{d}}{\mathrm{d}t} - \frac{1}{\mathrm{i}\hbar}\hat{H}\right)|\varphi\rangle =0\ .\]
Because the time derivation and the projector commute, we have
\[ \forall |\phi\rangle \in \mathcal{H}, \ \langle\phi|\left(\frac{\mathrm{d}}{\mathrm{d}t} - \frac{1}{\mathrm{i}\hbar}\hat{P}_{T_{|\varphi\rangle}\mathcal{M}}\hat{H}\right)|\varphi\rangle =0\ .\]
This equation must be true for any $|\phi\rangle \in \mathcal{H}$, Eq.~(\ref{eq:DiracFrenkel1}) can thus be written
\[ \left(\frac{\mathrm{d}}{\mathrm{d}t} - \frac{1}{\mathrm{i}\hbar}\hat{P}_{T_{|\varphi\rangle}\mathcal{M}}\hat{H}\right)|\varphi\rangle =0\ .\]
In the context of MPS, the manifold $\mathcal{M}$ will correspond to the space of full-ranked MPS of a given bond dimension $D$, and the tangent space will be the space spanned by variations of single MPS tensors. The projector defines an effective Hamiltonian under which the dynamics are constrained on $\mathcal{M}$. Constraining the dynamics on a manifold, introduces a projection error: the time evolution will obey to an effective Hamiltonian different from the starting one. After the introduction of TDVP as a time evolution method for MPS, Haegeman et al. pointed out [haegeman_unifying_2016] that there exist an analytical decomposition for the projector operator $\hat{\mathcal{P}}$ that simplifies the resolution of the equation, turning the problem into one where each matrix $A_i$ can be updated with an effective on site Hamiltonian $\hat H_\text{eff}$ via a Schroedinger like equation. The effective Hamiltonian $\hat H_\text{eff}$ is a contraction of the Hamiltonian MPO and the current state of the other matrices composing the MPS. This allows to do a sequential update.
There exist different versions of the TDVP algorithm. In MPSDynamics.jl three methods have been so far implemented:
The main advantage of the one-site 1TDVP algorithm is that it preserves the unitarity of the MPS during the time evolution. Its main problem, conversely, is that the time evolution is constrained to happen on a manifold constituted by tensors of fixed bond dimension, a quantity closely related to the amount of entanglement in the MPS, and such a bond dimension has therefore to be fixed before the beginning of the time evolution. This strategy will necessarily be non optimal: the growth of the bond dimensions required to describe the quantum state should ideally mirror the entanglement growth induced by the time evolution. 2TDVP allows for such a dynamical growth of the bond dimensions, and therefore better describes the entanglement in the MPS. It suffers however of other drawbacks: first of all, a truncation error is introduced (by the means of an SVD decomposition), which entails a loss of unitarity of the time-evolved MPS. Moreover, 2TDVP has bad scaling properties with the size of the local dimensions of the MPS: this is a major issue when dealing with bosons. The DTDVP algorithm combines the best features of 1TDVP and 2TDVP: it preserves unitarity, it has the same scaling properties of 1TDVP, and it adapts the bond dimensions to the entanglement evolution at each site and at each time-step. DTDVP does not suffer from a truncation error, but introduces only a projection error.
Chin, A. W.; Rivas, Á.; Huelga, S. F.; Plenio, M. B. Exact Mapping between System-Reservoir Quantum Models and Semi-Infinite Discrete Chains Using Orthogonal Polynomials. Journal of Mathematical Physics 2010, 51 (9), 092109. https://doi.org/10.1063/1.3490188.
Chin, A. W.; Prior, J.; Rosenbach, R.; Caycedo-Soler, F.; Huelga, S. F.; Plenio, M. B. The Role of Non-Equilibrium Vibrational Structures in Electronic Coherence and Recoherence in Pigment–Protein Complexes. Nature Phys 2013, 9 (2), 113–118. https://doi.org/10.1038/nphys2515.
Alvertis, A. M.; Schröder, F. A. Y. N.; Chin, A. W. Non-Equilibrium Relaxation of Hot States in Organic Semiconductors: Impact of Mode-Selective Excitation on Charge Transfer. J. Chem. Phys. 2019, 151 (8), 084104. https://doi.org/10.1063/1.5115239.
Dunnett, A. J.; Gowland, D.; Isborn, C. M.; Chin, A. W.; Zuehlsdorff, T. J. Influence of Non-Adiabatic Effects on Linear Absorption Spectra in the Condensed Phase: Methylene Blue. J. Chem. Phys. 2021, 155 (14), 144112. https://doi.org/10.1063/5.0062950.
Caycedo-Soler, F.; Mattioni, A.; Lim, J.; Renger, T.; Huelga, S. F.; Plenio, M. B. Exact Simulation of Pigment-Protein Complexes Unveils Vibronic Renormalization of Electronic Parameters in Ultrafast Spectroscopy. Nat Commun 2022, 13 (1), 2912. https://doi.org/10.1038/s41467-022-30565-4.
Orus, R. A Practical Introduction to Tensor Networks: Matrix Product States and Projected Entangled Pair States. Annals of Physics 2014, 349, 117–158. https://doi.org/10.1016/j.aop.2014.06.013.
Paeckel, S.; Köhler, T.; Swoboda, A.; Manmana, S. R.; Schollwöck, U.; Hubig, C. Time-Evolution Methods for Matrix-Product States. Annals of Physics 2019, 411, 167998. https://doi.org/10.1016/j.aop.2019.167998.
Prior, J.; Chin, A. W.; Huelga, S. F.; Plenio, M. B. Efficient Simulation of Strong System-Environment Interactions. Phys. Rev. Lett. 2010, 105 (5), 050404. https://doi.org/10.1103/PhysRevLett.105.050404.
Woods, M. P.; Cramer, M.; Plenio, M. B. Simulating Bosonic Baths with Error Bars. Phys. Rev. Lett. 2015, 115 (13), 130401. https://doi.org/10.1103/PhysRevLett.115.130401.
Tamascelli, D.; Smirne, A.; Lim, J.; Huelga, S. F.; Plenio, M. B. Efficient Simulation of Finite-Temperature Open Quantum Systems. Phys. Rev. Lett. 2019, 123 (9), 090402. https://doi.org/10.1103/PhysRevLett.123.090402.
Dunnett, A. J.; Chin, A. W. Simulating Quantum Vibronic Dynamics at Finite Temperatures With Many Body Wave Functions at 0 K. Front. Chem. 2021, 8. https://doi.org/10.3389/fchem.2020.600731.
Lacroix, T.; Dunnett, A.; Gribben, D.; Lovett, B. W.; Chin, A. Unveiling Non-Markovian Spacetime Signaling in Open Quantum Systems with Long-Range Tensor Network Dynamics. Phys. Rev. A 2021, 104 (5), 052204. https://doi.org/10.1103/PhysRevA.104.052204.
Appel, W. Mathematics for Physics and Physicists; Princeton University Press, 2007.
de Vega, I.; Banuls, M-.C. Thermofield-based chain-mapping approach for open quantum systems. Phys. Rev. A 2015, 92 (5), 052116. https://doi.org/10.1103/PhysRevA.92.052116.
Orus, R. A Practical Introduction to Tensor Networks: Matrix Product States and Projected Entangled Pair States. Annals of Physics 2014, 349, 117–158. https://doi.org/10.1016/j.aop.2014.06.013.
Poulin, D.; Qarry, A.; Somma, R.; Verstraete, F. Quantum Simulation of Time-Dependent Hamiltonians and the Convenient Illusion of Hilbert Space. Phys. Rev. Lett. 2011, 106 (17), 170501. https://doi.org/10.1103/PhysRevLett.106.170501.
Srednicki, M. Entropy and Area. Phys. Rev. Lett. 1993, 71 (5), 666–669. https://doi.org/10.1103/PhysRevLett.71.666.
Vidal, G.; Latorre, J. I.; Rico, E.; Kitaev, A. Entanglement in Quantum Critical Phenomena. Phys. Rev. Lett. 2003, 90 (22), 227902. https://doi.org/10.1103/PhysRevLett.90.227902.
Wolf, M. M.; Verstraete, F.; Hastings, M. B.; Cirac, J. I. Area Laws in Quantum Systems: Mutual Information and Correlations. Phys. Rev. Lett. 2008, 100 (7), 070502. https://doi.org/10.1103/PhysRevLett.100.070502.
Dirac, P. A. M. Note on Exchange Phenomena in the Thomas Atom. Mathematical Proceedings of the Cambridge Philosophical Society 1930, 26 (3), 376–385. https://doi.org/10.1017/S0305004100016108.
Frenkel, Y. Wave Mechanics, Advanced General Theory; Oxford, 1934; Vol. 1.
Haegeman, J.; Cirac, J. I.; Osborne, T. J.; Pižorn, I.; Verschelde, H.; Verstraete, F. Time-Dependent Variational Principle for Quantum Lattices. Phys. Rev. Lett. 2011, 107 (7), 070601. https://doi.org/10.1103/PhysRevLett.107.070601.
Haegeman, J.; Lubich, C.; Oseledets, I.; Vandereycken, B.; Verstraete, F. Unifying Time Evolution and Optimization with Matrix Product States. Phys. Rev. B 2016, 94 (16), 165116. https://doi.org/10.1103/PhysRevB.94.165116.
Raab, A. On the Dirac–Frenkel/McLachlan Variational Principle. Chemical Physics Letters 2000, 319 (5), 674–678. https://doi.org/10.1016/S0009-2614(00)00200-1.
Dunnett, A. J.; Chin A. W. Efficient bond-adaptive approach for finite-temperature open quantum dynamics using the one-site time-dependent variational principle for matrix product states. Physical Review B, 104(21):214302, Dec 2021. doi:10.1103/PhysRevB.104.214302.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
where $\hat L_S = \hat L_S^\dagger$, considering $\hat \rho_S(0) \otimes \ket{\Omega}\bra{\Omega}$ as the initial state of system and environment. It is this Hamiltonian that is mapped on two chains with TEDOPA, to be able to perform the time evolution using tensor network techniques.
A useful property of the orthonormal polynomials is that they obey a recurrence relation
\[ P_n(\omega) = (C_{n-1}\omega - A_{n-1})P_{n-1}(\omega) + B_{n-1}P_{n-2}(\omega)\ ,\]
where $A_n$ is related to the first moment of $P_n$, $B_n$ and $C_n$ to the norms of $P_n$ and $P_{n-1}$[appel_mathematics_2007]. This recurrence relation can be used to construct the polynomials with the conditions that $P_0(\omega) = ||p_0||^{-1} = \left(\int_{\mathbb{R}^{+}} J(\omega)\mathrm{d}\omega \right)^{-\frac{1}{2}}$ and $P_{-1}(\omega) = 0$, with $||\bullet||$ the norm of $\bullet$ with respect to the measure $J(\omega)$, and $P_n(\omega) = p_n(\omega)||p_n||^{-1}$ ; where the polynomials $\{p_n\}_{n\in\mathbb{N}}$ are the so called monic polynomials where the factor $a_n$ in front of $\omega^{n}$ is equal to 1.
The energy of the chain mode $n$ is given by $\varepsilon_n = A_n C_n^{-1}$ and $t_n=C_n^{-1}$ is the coupling between mode $n$ and $n+1$[chin_exact_2010].
The system couples only to the first mode with the coupling strength $\kappa = ||p_0||$.
Explain that for some weight function/SD they are known analytically and that for others we can use the build-in routines inspired by Gautschi or the PolyChaos.jl package.
A multipartite quantum state $|\psi\rangle$, e.g. a $N$-site system where the sites can each be in a state $|\phi_i\rangle$ belonging to a $d$-dimensional Hilbert space, can be written as follows
\[ |\psi\rangle = \sum_{\{i_k\}}c_{i_1\ldots i_N}|\phi_{i_1}\rangle\otimes\ldots\otimes|\phi_{i_N}\rangle\ ,\]
where the complex numbers $c_{i_1\ldots i_N}$ are the amplitudes of each state $|\phi_{i_1}\rangle\otimes\ldots\otimes|\phi_{i_N}\rangle$ whose superpositions form in full generality the state $|\psi\rangle$. Thus the state $|\psi\rangle$ can be completely represented by a rank-$N$ tensor $c$ that is the collection of all possible amplitudes $c_{i_1\ldots i_N}$. Here by the rank of a tensor, we simply mean the number of indices it has.
The tensor $c$ of a quantum state $|\psi\rangle$ corresponding to a one-dimensional system can be decomposed into a product of $N$ smaller rank-3 tensors $T_{k}$ (except for the first and last sites where the tensors will have a rank-2)
\[ c_{i_1\ldots i_N} = \sum_{\{\alpha\}} T^{\alpha_1}_{i_1}T^{\alpha_1\alpha_2\ }_{i_2}T^{\alpha_2\alpha_3\ }_{i_3}\ldots T^{\alpha_{N-1}}_{i_N} \ .\]
In this form, the local tensor $T_k$ contains the information on the quantum state on site $k$ and its relation (especially the entanglement) with the neighbouring sites.
The decomposition of the tensor of the amplitudes of a quantum state into a product of smaller rank tensors is called a Matrix Product State decomposition.
The contracted indices $\alpha_k$ between the tensors are called virtual indices and carry information about the correlations between bi-partitions of the state at bond $k$. The number of different values a virtual index can take is called the bond dimension and is denoted $D$. The free indices $i_k$ associated with local quantum states are called physical indices. Thus, they can take $d$ values (with $d$ the dimension of the local Hilbert space).
Any state in the Hilbert space of a one-dimensional many-body system can in principle be represented by a MPS by choosing a sufficiently large value for the bond dimension $D$ [Orus]. On top of this intellectually satisfying property of MPSs being a dense set of states for a 1d-system, they can also be used as a practical Ansätze for a many-body quantum states by setting a maximal allowed value $\chi$ for the bond dimension $D$. In doing so, we restrict ourselves to a corner of the total Hilbert space. The rationale behind this Ansatz is the following: if the initial quantum state of a many-body system has a low bond dimension (typically if the initial state is a product state with $D = 1$), then in a finite time it will only be able to explore a region of the Hilbert space that is not to far away from its starting point. Thus, the bond dimension will not have the time to diverge exponentially [poulin_quantum_2011]. However, depending on the physical system at hand, this sub-manifold of the Hilbert space could still be "too large". There is an additional reason that explains why MPSs are good Ansätze for 1d physical systems. Most many-body Hamiltonians we (physicists) are interested in are local, meaning that the interactions they describe involve objects that are "neighbours". For such Hamiltonians, the ground states (outside of potential critical phases) follow the so called area law for the entanglement entropy [srednicki_entropy_1993][vidal_entanglement_2003][wolf_area_2008]. This law states that the entanglement entropy $S_{vN}$ of a bi-partition of the system is proportional, not to the volume of the partition as one might expect, but to the hyper-surface of the partition's boundary; hence the name "area law". For a 3d system this corresponds to an actual surface area $A$, $S_{vN} \sim A$; for a 2d system it corresponds to the length $L$ of the partition's boundary, $S_{vN} \sim L$; and in 1d the boundary reduces to a point, thus the entropy will be independent of the size of the system $S_{vN} \sim \text{constant}$. The MPSs are states that satisfy this area law.
An application of the Singular Value Decomposition is to create efficient approximations of quantum states to perform computations. The main idea is to reduce the content of the MPS to keep only the parts that contain the physics of interest. One method to realise this approximation is to do a SVD on each of the tensors of the MPS after each time step of the state time-evolution and to trim the smallest singular values in order to decrease the bond dimension of the MPS down to a chosen maximal value $\chi$. The corresponding columns and rows of the unitary matrices $U$ and $V^\dagger$ are also removed. Then, the trimmed matrices $\tilde{U}$, $\tilde{S}$ and $\tilde{V}^\dagger$ are contracted back to give an approximated tensor $T$ with a smaller bond dimension. Another way to apply the restricted rank approximation is to restrict oneself into working in a manifold of fixed bond dimension $D$ and to use methods that can enforce this constraint.
In order to compute expectation values of observables or apply unitary transformations to a quantum state, we need a TN representation of operators. In the same fashion as a one-dimensional quantum state can be represented as a MPS, operators acting on those states can be represented as Matrix Product Operators (MPO). For an operator $\hat{O}$, its MPO can be defined as follows
\[ \hat{O} = \sum_{\{i_k\}\{i_k^{'}\} \{w\}} W^{i_1\ i^{'}_1}_{1\ w_0w_1}\ldots W^{i_N\ i^{'}_N}_{N\ w_{N-1}w_N} |\phi_{i_1^{'}}\ldots \phi_{i_N^{'}}\rangle\langle\phi_{i_1}\ldots \phi_{i_N}| \]
The contracted indices between the tensors are called virtual indices. The free indices are called physical indices and correspond to the different input and output local quantum states. They can take $d$ values (with $d$ the dimension of the local Hilbert space).
A natural extension to the MPS is the (loop-free) tree tensor network. A TTN is a generalisation of the MPS wherein each site, instead of being connected to only one other site to its right, may be connected to any arbitrary number of child sites. Provided the tree does not contain any loops, everything that one can do to an MPS/MPO can be extended straight-forwardly to TTN states and TTN operators. The generalisation to trees introduces no new conceptual complexity (only implementational complexity). The sites of a TTN are usually referred to as nodes. For our purposes, every node of a TTN state and operator has one parent leg, and any number (including zero) of child legs. The first node is known as the head-node and has a dummy parent leg with dimension 1.
The original idea behind TDVP goes back to Dirac [dirac_note_1930] and Frenkel [frenkel_wave_1934]. The main point, in the modern tensor networks formulation, is that instead of solving the Schrödinger equation and then truncating the MPS representation of the quantum state, one can solve the equations of motion projected into a space of restricted bond dimension [haegeman_timedependent_2011][haegeman_unifying_2016].
The general formulation of the Dirac-Frenkel Variational Principle [raab_diracfrenkelmclachlan_2000] is that one looks for a solution $|\varphi\rangle \in \mathcal{M}$ of the Schrödinger equation where $\mathcal{M} \subset \mathcal{H}$ is a manifold of the total Hilbert space $\mathcal{H}$ in which we think that the relevant physical states `live'.
We define $T_{|\varphi\rangle}\mathcal{M}$ the tangent space of $\mathcal{M}$ around the state $|\varphi\rangle$. The criterion to find $|\varphi\rangle$ is that for every state $|\chi\rangle \in T_{|\varphi\rangle}\mathcal{M}$
\[ \langle\chi|\left(\frac{\mathrm{d}}{\mathrm{d}t} - \frac{1}{\mathrm{i}\hbar}\hat{H}\right)|\varphi\rangle =0\ ,\]
which can be interpreted as saying that the time evolution procedure should keep $|\varphi\rangle$ inside of the manifold $\mathcal{M}$.
The term variational in the name of the method comes from the fact that in practice one aims at minimising the right-hand side of the hereinabove Eq. to find $|\varphi\rangle$.
Introducing $\hat{P}_{T_{|\varphi\rangle}\mathcal{M}}$ the projector onto the tangent space $T_{|\varphi\rangle}\mathcal{M}$, we can write the state $|\chi\rangle = \hat{P}_{T_{|\varphi\rangle}\mathcal{M}}|\phi\rangle$ with $|\phi\rangle$ a state in $\mathcal{H}$. Leading to
\[ \forall |\phi\rangle \in \mathcal{H}, \ \langle\phi|\hat{P}_{T_{|\varphi\rangle}\mathcal{M}}\left(\frac{\mathrm{d}}{\mathrm{d}t} - \frac{1}{\mathrm{i}\hbar}\hat{H}\right)|\varphi\rangle =0\ .\]
Because the time derivation and the projector commute, we have
\[ \forall |\phi\rangle \in \mathcal{H}, \ \langle\phi|\left(\frac{\mathrm{d}}{\mathrm{d}t} - \frac{1}{\mathrm{i}\hbar}\hat{P}_{T_{|\varphi\rangle}\mathcal{M}}\hat{H}\right)|\varphi\rangle =0\ .\]
This equation must be true for any $|\phi\rangle \in \mathcal{H}$, Eq.~(\ref{eq:DiracFrenkel1}) can thus be written
\[ \left(\frac{\mathrm{d}}{\mathrm{d}t} - \frac{1}{\mathrm{i}\hbar}\hat{P}_{T_{|\varphi\rangle}\mathcal{M}}\hat{H}\right)|\varphi\rangle =0\ .\]
In the context of MPS, the manifold $\mathcal{M}$ will correspond to the space of full-ranked MPS of a given bond dimension $D$, and the tangent space will be the space spanned by variations of single MPS tensors. The projector defines an effective Hamiltonian under which the dynamics are constrained on $\mathcal{M}$. Constraining the dynamics on a manifold, introduces a projection error: the time evolution will obey to an effective Hamiltonian different from the starting one. After the introduction of TDVP as a time evolution method for MPS, Haegeman et al. pointed out [haegeman_unifying_2016] that there exist an analytical decomposition for the projector operator $\hat{\mathcal{P}}$ that simplifies the resolution of the equation, turning the problem into one where each matrix $A_i$ can be updated with an effective on site Hamiltonian $\hat H_\text{eff}$ via a Schroedinger like equation. The effective Hamiltonian $\hat H_\text{eff}$ is a contraction of the Hamiltonian MPO and the current state of the other matrices composing the MPS. This allows to do a sequential update.
There exist different versions of the TDVP algorithm. In MPSDynamics.jl three methods have been so far implemented:
The main advantage of the one-site 1TDVP algorithm is that it preserves the unitarity of the MPS during the time evolution. Its main problem, conversely, is that the time evolution is constrained to happen on a manifold constituted by tensors of fixed bond dimension, a quantity closely related to the amount of entanglement in the MPS, and such a bond dimension has therefore to be fixed before the beginning of the time evolution. This strategy will necessarily be non optimal: the growth of the bond dimensions required to describe the quantum state should ideally mirror the entanglement growth induced by the time evolution. 2TDVP allows for such a dynamical growth of the bond dimensions, and therefore better describes the entanglement in the MPS. It suffers however of other drawbacks: first of all, a truncation error is introduced (by the means of an SVD decomposition), which entails a loss of unitarity of the time-evolved MPS. Moreover, 2TDVP has bad scaling properties with the size of the local dimensions of the MPS: this is a major issue when dealing with bosons. The DTDVP algorithm combines the best features of 1TDVP and 2TDVP: it preserves unitarity, it has the same scaling properties of 1TDVP, and it adapts the bond dimensions to the entanglement evolution at each site and at each time-step. DTDVP does not suffer from a truncation error, but introduces only a projection error.
Chin, A. W.; Rivas, Á.; Huelga, S. F.; Plenio, M. B. Exact Mapping between System-Reservoir Quantum Models and Semi-Infinite Discrete Chains Using Orthogonal Polynomials. Journal of Mathematical Physics 2010, 51 (9), 092109. https://doi.org/10.1063/1.3490188.
Chin, A. W.; Prior, J.; Rosenbach, R.; Caycedo-Soler, F.; Huelga, S. F.; Plenio, M. B. The Role of Non-Equilibrium Vibrational Structures in Electronic Coherence and Recoherence in Pigment–Protein Complexes. Nature Phys 2013, 9 (2), 113–118. https://doi.org/10.1038/nphys2515.
Alvertis, A. M.; Schröder, F. A. Y. N.; Chin, A. W. Non-Equilibrium Relaxation of Hot States in Organic Semiconductors: Impact of Mode-Selective Excitation on Charge Transfer. J. Chem. Phys. 2019, 151 (8), 084104. https://doi.org/10.1063/1.5115239.
Dunnett, A. J.; Gowland, D.; Isborn, C. M.; Chin, A. W.; Zuehlsdorff, T. J. Influence of Non-Adiabatic Effects on Linear Absorption Spectra in the Condensed Phase: Methylene Blue. J. Chem. Phys. 2021, 155 (14), 144112. https://doi.org/10.1063/5.0062950.
Caycedo-Soler, F.; Mattioni, A.; Lim, J.; Renger, T.; Huelga, S. F.; Plenio, M. B. Exact Simulation of Pigment-Protein Complexes Unveils Vibronic Renormalization of Electronic Parameters in Ultrafast Spectroscopy. Nat Commun 2022, 13 (1), 2912. https://doi.org/10.1038/s41467-022-30565-4.
Orus, R. A Practical Introduction to Tensor Networks: Matrix Product States and Projected Entangled Pair States. Annals of Physics 2014, 349, 117–158. https://doi.org/10.1016/j.aop.2014.06.013.
Paeckel, S.; Köhler, T.; Swoboda, A.; Manmana, S. R.; Schollwöck, U.; Hubig, C. Time-Evolution Methods for Matrix-Product States. Annals of Physics 2019, 411, 167998. https://doi.org/10.1016/j.aop.2019.167998.
Prior, J.; Chin, A. W.; Huelga, S. F.; Plenio, M. B. Efficient Simulation of Strong System-Environment Interactions. Phys. Rev. Lett. 2010, 105 (5), 050404. https://doi.org/10.1103/PhysRevLett.105.050404.
Woods, M. P.; Cramer, M.; Plenio, M. B. Simulating Bosonic Baths with Error Bars. Phys. Rev. Lett. 2015, 115 (13), 130401. https://doi.org/10.1103/PhysRevLett.115.130401.
Tamascelli, D.; Smirne, A.; Lim, J.; Huelga, S. F.; Plenio, M. B. Efficient Simulation of Finite-Temperature Open Quantum Systems. Phys. Rev. Lett. 2019, 123 (9), 090402. https://doi.org/10.1103/PhysRevLett.123.090402.
Dunnett, A. J.; Chin, A. W. Simulating Quantum Vibronic Dynamics at Finite Temperatures With Many Body Wave Functions at 0 K. Front. Chem. 2021, 8. https://doi.org/10.3389/fchem.2020.600731.
Lacroix, T.; Dunnett, A.; Gribben, D.; Lovett, B. W.; Chin, A. Unveiling Non-Markovian Spacetime Signaling in Open Quantum Systems with Long-Range Tensor Network Dynamics. Phys. Rev. A 2021, 104 (5), 052204. https://doi.org/10.1103/PhysRevA.104.052204.
Appel, W. Mathematics for Physics and Physicists; Princeton University Press, 2007.
de Vega, I.; Banuls, M-.C. Thermofield-based chain-mapping approach for open quantum systems. Phys. Rev. A 2015, 92 (5), 052116. https://doi.org/10.1103/PhysRevA.92.052116.
Orus, R. A Practical Introduction to Tensor Networks: Matrix Product States and Projected Entangled Pair States. Annals of Physics 2014, 349, 117–158. https://doi.org/10.1016/j.aop.2014.06.013.
Poulin, D.; Qarry, A.; Somma, R.; Verstraete, F. Quantum Simulation of Time-Dependent Hamiltonians and the Convenient Illusion of Hilbert Space. Phys. Rev. Lett. 2011, 106 (17), 170501. https://doi.org/10.1103/PhysRevLett.106.170501.
Srednicki, M. Entropy and Area. Phys. Rev. Lett. 1993, 71 (5), 666–669. https://doi.org/10.1103/PhysRevLett.71.666.
Vidal, G.; Latorre, J. I.; Rico, E.; Kitaev, A. Entanglement in Quantum Critical Phenomena. Phys. Rev. Lett. 2003, 90 (22), 227902. https://doi.org/10.1103/PhysRevLett.90.227902.
Wolf, M. M.; Verstraete, F.; Hastings, M. B.; Cirac, J. I. Area Laws in Quantum Systems: Mutual Information and Correlations. Phys. Rev. Lett. 2008, 100 (7), 070502. https://doi.org/10.1103/PhysRevLett.100.070502.
Dirac, P. A. M. Note on Exchange Phenomena in the Thomas Atom. Mathematical Proceedings of the Cambridge Philosophical Society 1930, 26 (3), 376–385. https://doi.org/10.1017/S0305004100016108.
Frenkel, Y. Wave Mechanics, Advanced General Theory; Oxford, 1934; Vol. 1.
Haegeman, J.; Cirac, J. I.; Osborne, T. J.; Pižorn, I.; Verschelde, H.; Verstraete, F. Time-Dependent Variational Principle for Quantum Lattices. Phys. Rev. Lett. 2011, 107 (7), 070601. https://doi.org/10.1103/PhysRevLett.107.070601.
Haegeman, J.; Lubich, C.; Oseledets, I.; Vandereycken, B.; Verstraete, F. Unifying Time Evolution and Optimization with Matrix Product States. Phys. Rev. B 2016, 94 (16), 165116. https://doi.org/10.1103/PhysRevB.94.165116.
Raab, A. On the Dirac–Frenkel/McLachlan Variational Principle. Chemical Physics Letters 2000, 319 (5), 674–678. https://doi.org/10.1016/S0009-2614(00)00200-1.
Dunnett, A. J.; Chin A. W. Efficient bond-adaptive approach for finite-temperature open quantum dynamics using the one-site time-dependent variational principle for matrix product states. Physical Review B, 104(21):214302, Dec 2021. doi:10.1103/PhysRevB.104.214302.
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Here we explain the different steps to perform a simulation.
Examples with detailed explanations can be found in Examples.
The initial many-body state of the {System + Environment} can be described easily as Matrix Product States (MPS) or Tree Tensor Networks (TTN) state (e.g. when a system is coupled to several environments).
A MPS can be initialized with several methods.
The productstatemps method enables to instantiate arbitrary MPS of fixed uniform bond dimension with non-uniform physical dimensions. The individual states of the MPS sites can be provided by setting state to a list of column vectors. Setting state=:Vacuum will produce an MPS in the vacuum state. Setting state=:FullOccupy will produce an MPS in which each site is fully occupied. The gauge of the MPS can also be set using a keyword argument.
julia> ψ = unitcol(1,2) # system initial state
+User Guide · MPSDynamics.jl User Guide
Here we explain the different steps to perform a simulation.
Examples with detailed explanations can be found in Examples.
Initial State
The initial many-body state of the {System + Environment} can be described easily as Matrix Product States (MPS) or Tree Tensor Networks (TTN) state (e.g. when a system is coupled to several environments).
Matrix Product States
A MPS can be initialized with several methods.
The productstatemps method enables to instantiate arbitrary MPS of fixed uniform bond dimension with non-uniform physical dimensions. The individual states of the MPS sites can be provided by setting state to a list of column vectors. Setting state=:Vacuum will produce an MPS in the vacuum state. Setting state=:FullOccupy will produce an MPS in which each site is fully occupied. The gauge of the MPS can also be set using a keyword argument.
julia> ψ = unitcol(1,2) # system initial state
julia> d = 6; N = 30; α = 0.1; Δ = 0.0; ω0 = 0.2; s = 1
@@ -34,4 +34,4 @@
"nchain" => [0.0 0.23466 … 1.84319 1.76098; 0.0 0.00231507 … 0.83105 0.9033…
"sy" => [0.0, -0.0133489, -0.0588887, -0.0858181, -0.0759996, -0.048539…
"times" => [0.0, 0.0005, 0.001, 0.0015, 0.002, 0.0025, 0.003, 0.0035, 0.00…
-
- Gautschi
Gautschi, W. Algorithm 726: ORTHPOL–a package of routines for generating orthogonal polynomials and Gauss-type quadrature rules. ACM Trans. Math. Softw. 20, 21–62 (1994).
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.
Gautschi, W. Algorithm 726: ORTHPOL–a package of routines for generating orthogonal polynomials and Gauss-type quadrature rules. ACM Trans. Math. Softw. 20, 21–62 (1994).
Settings
This document was generated with Documenter.jl on Tuesday 14 May 2024. Using Julia version 1.7.3.