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 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)}_{\hat H_\text{hyb}} + \underbrace{\sum_k \epsilon_k \hat c_k^\dagger \hat c_k}_{\hat H_\text{cond}}.\]
This Hamiltonian represents a located impurity ($\hat H_\text{loc}$), conduction electrons ($\hat H_\text{cond}$) and a hybridization term between the impurity and the conduction electrons ($\hat H_\text{hyb}$). 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 , 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 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{aligned}
&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}},
@@ -69,11 +69,11 @@
\begin{bmatrix}
E_{1,N} \hat b_{1,N}^\dagger \hat b_{1,N} \\ 0 \\0 \\ J_{1,N} \hat b_{1,N}^\dagger \\ J_{1,N} \hat b_{1,N} \\ \hat{\mathbb I}
\end{bmatrix} .
-\end{aligned}\]
The fermionic mapping can be strightforwardly implemented with the methods of MPSDyanmics.jl
. We provide two examples:
examples/anderson_model_double
, simulating the SIAM with the double-chain geometryexamples/anderson_model_interleaved
, simulating the SIAM with the interleaved-chain geometry
They only differ on the way the Hamiltonian MPO is defined. We now briefly review the important bits of the code.
Both of the examples start with the definition of the physical parameters,
N = 40 # number of chain sites
+\end{aligned}\]The fermionic mapping can be strightforwardly implemented with the methods of MPSDynamics.jl
. We provide two examples:
examples/anderson_model_double
, simulating the SIAM with the double-chain geometryexamples/anderson_model_interleaved
, simulating the SIAM with the interleaved-chain geometry
They only differ on the way the Hamiltonian MPO is defined. We now briefly review the important bits of the code.
Both of the examples start with the definition of the physical parameters,
N = 40 # number of chain sites
β = 10.0 # inverse temperature
μ = 0. # chemical potential
Ed = 0.3 # energy of the impurity
-ϵd = Ed - μ # energy of the impurity minus the chemical potential
The chaincoeffs_fermionic
function is needed to compute the chain coefficients. It requires as inputs the number of modes of each chain N
, the inverse temperature β
, a label to specify if the chain modes are empty (label is 1.0
) or filled (label is 2.0
), and both the dispersion relation $\epsilon_k$ and the fermionic spectral density funciton $V_k$.
function ϵ(x)
+ϵd = Ed - μ # energy of the impurity minus the chemical potential
The MPSDynamics.chaincoeffs_fermionic
function is needed to compute the chain coefficients. It requires as inputs the number of modes of each chain N
, the inverse temperature β
, a label to specify if the chain modes are empty (label is 1.0
) or filled (label is 2.0
), and both the dispersion relation $\epsilon_k$ and the fermionic spectral density funciton $V_k$.
function ϵ(x)
return x
end
@@ -87,7 +87,7 @@
method = :DTDVP # time-evolution method
Dmax = 150 # MPS max bond dimension
prec = 0.0001 # precision for the adaptive TDVP
-
and with this we are ready to construct the Hamiltonian MPO and specify the initial state, which will obviously differ depending on the chosen geometry.
The Hamiltonian is defined using the tightbinding_mpo
function, which takes as an input the number of modes of each chain N
, the defect's energy ϵd
, and the chain coefficients of the first chainparams1
and second chainparams2
chain. The MPS for the initial state is a factorized state made of: N filled states, a filled impurity, and N empty states.
H = tightbinding_mpo(N, ϵd, chainparams1, chainparams2)
+
and with this we are ready to construct the Hamiltonian MPO and specify the initial state, which will obviously differ depending on the chosen geometry.
The Hamiltonian is defined using the MPSDynamics.tightbinding_mpo
function, which takes as an input the number of modes of each chain N
, the defect's energy ϵd
, and the chain coefficients of the first chainparams1
and second chainparams2
chain. The MPS for the initial state is a factorized state made of: N filled states, a filled impurity, and N empty states.
H = tightbinding_mpo(N, ϵd, chainparams1, chainparams2)
ψ = unitcol(2,2) # (0,1) filled impurity state
A = productstatemps(physdims(H), state=[fill(unitcol(2,2), N)..., ψ, fill(unitcol(1,2), N)...]) # MPS
To avoid the DTDVP
algorithm from getting stuck in a local minimum, it is better to embed the MPS in a manifold of bond dimension 2 (or more):
mpsembed!(A, 2) # to embed the MPS in a manifold of bond dimension 2
We can now define the observables for the two chains and for the impurity
ob1 = OneSiteObservable("chain1_filled_occup", numb(2), (1,N))
@@ -97,7 +97,7 @@
method = method,
obs = [ob1, ob2, ob3],
convobs = [ob1],
- params = @LogParams(N, ϵd, β, c1, c2),
+ params = @LogParams(N, ϵd, β),
convparams = [prec],
Dlim = Dmax,
savebonddims = true, # we want to save the bond dimension
@@ -162,7 +162,7 @@
# Display the plots
plot(p2, p3, p4, p5, p1, layout = (3, 2), size = (1400, 1200))
-
The Hamiltonian is defined using the interleaved_tightbinding_mpo
function, which takes as an input the number of modes of each chain N
, the defect's energy ϵd
, and the chain coefficients of the first chainparams1
and second chainparams2
chain. The MPS for the initial state is a factorized state (bond dimension 1) made of: a filled impurity, and 2N alternate filled-empty states.
H = interleaved_tightbinding_mpo(N, ϵd, chainparams1, chainparams2)
+
The Hamiltonian is defined using the MPSDynamics.interleaved_tightbinding_mpo
function, which takes as an input the number of modes of each chain N
, the defect's energy ϵd
, and the chain coefficients of the first chainparams1
and second chainparams2
chain. The MPS for the initial state is a factorized state (bond dimension 1) made of: a filled impurity, and 2N alternate filled-empty states.
H = interleaved_tightbinding_mpo(N, ϵd, chainparams1, chainparams2)
ψ = unitcol(2,2) # (0,1) filled impurity state
Tot = Any[]
@@ -178,7 +178,7 @@
method = method,
obs = [ob1, ob2],
convobs = [ob1],
- params = @LogParams(N, ϵd, β, c1, c2),
+ params = @LogParams(N, ϵd, β),
convparams = [prec],
Dlim = Dmax,
savebonddims = true, # we want to save the bond dimension
@@ -266,4 +266,4 @@
end
# Display the plots
-plot(p2, p3, p4, p5, p1, layout = (3, 2), size = (1400, 1200))
________________