Time Dependent Variational Principle for Tree Tensor Networks

Daniel Bauernfeind1*, Markus Aichhorn1,

1 Institute of Theoretical and Computational Physics Graz University of
Technology, 8010 Graz, Austria

*<EMAIL_ADDRESS>
## Abstract

We present a generalization of the Time Dependent Variational Principle (TDVP)
to any finite sized loop-free tensor network. The major advantage of TDVP is
that it can be employed as long as a representation of the Hamiltonian in the
same tensor network structure that encodes the state is available. Often, such
a representation can be found also for long-range terms in the Hamiltonian. As
an application we use TDVP for the Fork Tensor Product States tensor network
for multi-orbital Anderson impurity models. We demonstrate that TDVP allows to
account for off-diagonal hybridizations in the bath which are relevant when
spin-orbit coupling effects are important, or when distortions of the crystal
lattice are present.

##  1 Introduction

The development of the Density Matrix Renormalization Group (DMRG) [1, 2] was
an immensely important milestone in our understanding of one-dimensional
quantum systems. The subsequent realizations that DMRG produces Matrix Product
States [3] (MPS) and that it can be formulated as a variational method [4],
ultimately led to the development of numerous approaches using not only MPS
but also general Tensor Networks to handle quantum systems. Notable examples
are the Projected Entangled Pair States (PEPS) [5, 6], the Multi-scale
Entanglement Renormalization Ansatz (MERA) [7] and so-called Tree-Tensor
Networks (TTN) [8, 9, 10, 11, 12, 13, 14, 15, 16, 17] including also the
recently developed Fork Tensor Product States (FTPS) method [18, 19].
Among the most important properties of tensor networks is whether their graph
is loop-free, i.e., whether there exists only a single path from one tensor to
any other. While PEPS and MERA are not loop-free, the TTNs and MPS are.
Cutting any edge of a loop-free network, results in two separated segments and
therefore gives a notion of left and right with respect to this edge. This in
turn allows a controlled truncation scheme based on the Schmidt-decomposition
of quantum states in the spirit of DMRG.
One of the major reasons behind the success of tensor networks are the
celebrated area laws of entanglement [20] stating that the entanglement of
ground states of gapped Hamiltonians with short-range couplings is
proportional to the surface area connecting the two regions. MPS in 1-d and
PEPS in 2-d efficiently encode quantum states obeying these area laws and are
hence efficient parametrizations. In addition, MPS-based time evolution for
one dimensional systems is an important method to calculate dynamical
properties [21, 22, 23, 24]. Approaches to perform the real-time evolution
include, among others, the Time-dependent Density Matrix Renormalization Group
(tDMRG) [25, 26], the closely related Time Evolving Block Decimation (TEBD)
[27, 28] as well as the Time Dependent Variational Principle (TDVP) [29, 30,
31]. An in depth comparison of several time evolution algorithms performed in
Ref. [32] came to the conclusion that while all approaches have strengths and
weaknesses, TDVP is among the most reliable methods to perform the time
evolution.
While time evolution approaches for MPS are well established, much less has
been done for general tensor networks. So far, mostly TEBD (and variations)
have been used, for example for the MERA network [33], for PEPS [34, 35, 36]
and for TTNs [9, 37, 38, 17, 18]. The advantage of TEBD is its relative
simplicity, since it effectively boils down to a repeated application of short
range operators obtained from a Suzuki-Trotter decomposition [39] of the full
time-evolution operator.
However, one of the major disadvantages of TEBD is that it can become
difficult to implement for more complicated Hamiltonians, especially when
long-range couplings are present. One approach to treat such couplings is an
MPO-based approach introduced by Zaletel et al. [44] in which a MPO
approximation of the time-evolution operator is constructed. Alternatively,
TDVP circumvents this problem by only demanding a Hamiltonian represented in
the same tensor network structure as the state which is often easy to find.
Additionally, TDVP in its single-site variant exactly respects conserved
quantities of the Hamiltonian like energy or magnetization [30]. Although some
works applied TDVP to more general tensor networks [45, 46], it is not obvious
how these algorithms work in detail and how it can be generalized. A notable
exception is Ref. [47] which introduces TDVP for binary TTNs. Parallel to
these developments in the tensor network physics community, very similar
approaches to TDVP have been developed in quantum chemistry under the name of
Multi-layer Multi-Configurational Time-Dependent Hartree approach [40, 41, 42,
43]. These methods effectively generate tensor networks by repeatedly grouping
degrees of freedom together and transforming them with (time dependent) basis
transformations into new degrees of freedom.
A more practical motivation for the formulation of TDVP for TTNs are Dynamical
Mean-Field Theory (DMFT) calculations using the FTPS tensor network. So far,
this approach has been used for so-called diagonal hybridizations only. On the
other hand, real materials often exhibit off-diagonal hybridizations, which
can for example come from spin-orbit coupling, or from distortions of the
crystal lattice. For off-diagonal hybridizations, the TEBD approach used so
far [18, 19] is difficult to generalize and we hence choose to use TDVP in
these situations.
Although part of the motivation for this work comes from the FTPS tensor
network, in this paper we formulate TDVP for _general_ loop-free and finite-
size tensor networks. After establishing the relevant concepts of TTNs in Sec.
2, we generalize TDVP to these networks in Sec. 3. Finally in Sec. 4 we show
how this approach can be used for the FTPS tensor network and that it can be
applied to off-diagonal hybridizations.

##  2 Tree Tensor Networks Basics

Figure 1:  Example of a TTN with 7 tensors with different numbers of link-
indices on each site. Each dot represents a tensor and each line an index,
where a connected line implies summation over this index. The open lines are
the physical indices $s_{1}\cdots s_{N}$ ($N=7$), while the connected lines
are the link indices $q_{1}\cdots q_{L}$ ($L=6$). Cutting the link between
sites $4$ and $5$, as indicated by the dashed line, results in two
disconnected tensor network segments and defines a notion of left and right at
each link. In this example, sites $1$, $2$, $3$, $6$ and $7$ are the leaves of
the TTN.

In this section, we discuss concepts of TTNs relevant for the formulation of
TDVP. All these properties are generalizations of the corresponding concepts
for MPS. Although these have been discussed previously in several publications
(see for example Refs. [17, 48, 14]), here, we present them in a format that
will suit us for the subsequent formulation of the TDVP algorithm.

###  2.1 TTNs

x Any state ${\ket{\psi}}$ of a quantum system consisting of $N$ sites with
local basis states $\ket{s_{i}}$ on site $i$ can be expanded in the
corresponding product basis:

${\ket{\psi}}=\sum_{s_{1}\cdots s_{N}}c_{s_{1}\cdots s_{N}}\ket{s_{1}\cdots
s_{N}}.$ (1)

The coefficient $c_{s_{1}\cdots s_{N}}$ is interpreted as a rank-N tensor with
indices $s_{1}\cdots s_{N}$. Tensor networks represent this rank-N tensor as a
product over tensors of much smaller rank:

$\displaystyle c_{s_{1}\cdots s_{N}}=\sum_{q_{1}\cdots
q_{L}}T^{s_{1}}_{Q_{1}}\cdot T^{s_{2}}_{Q_{2}}\cdots T^{s_{N}}_{Q_{N}}$
$\displaystyle\ket{\psi[T]}=\sum_{\stackrel{{\scriptstyle s_{1}\cdots
s_{N}}}{{q_{1}\cdots q_{L}}}}T^{s_{1}}_{Q_{1}}\cdot T^{s_{2}}_{Q_{2}}\cdots
T^{s_{N}}_{Q_{N}}\ket{s_{1}\cdots s_{N}}.$ (2)

Each tensor $T^{s_{i}}_{Q_{i}}\equiv T^{s_{i}}_{q_{1}q_{2}\cdots q_{r_{i}}}$
has a set of auxiliary indices $Q_{i}=\\{q_{k}:q_{k}\text{ is attached to node
$i$ }\\}$ such that each auxiliary index is part of exactly two tensors. We
call $r_{i}=|Q_{i}|$ the number of indices of the tensor on site $i$.
Additionally, we attached to each tensor a physical index as for example in
the FTPS tensor network. While for general TTNs not all tensors have a
physical index, the following results can be straightforwardly generalized by
just removing the physical index from the notation. Alternatively, every
tensor without a physical index could be interpreted as having a dummy index
with just a single entry corresponding to a single state, say $\ket{0}$, onto
which the Hamiltonian acts as an identity $H\ket{0}=\ket{0}$. Note that if all
sites have a physical index, the number of links is $L=N-1$. In the following,
we will often omit sums over auxiliary indices $\sum_{q_{1}\cdots q_{L}}$ and
assume Einstein convention for the summations.

Figure 2:  Gauge degree of freedom in tensor networks. At each link, one can
insert an identity $\mathbb{1}=G\cdot G^{-1}$ without changing the physical
state ${\ket{\psi}}$. By absorbing $G$ into one tensor and $G^{-1}$ into the
other, we obtain a different representation of the same state ${\ket{\psi}}$.
Figure 3: With an SVD, we can orthogonalize a tensor towards one of its
neighbors with which it shares an index. _Top:_ Tensor $T$ is reshaped into a
matrix and the $U$-matrix of its SVD is used as new tensor. $S\cdot
V^{\dagger}$ is absorbed into the neighboring tensor. _Middle:_ Graphical
representation of $U^{\dagger}\cdot U=\mathbb{1}$. _Bottom:_ A tensor that is
normalized towards one of its neighbors is depicted as a triangle pointing in
the direction of this neighbor. The picture shows tensor 4 of the TTN in Fig.
1 orthogonalized towards tensor $2$. Let us call the index connecting tensor
$4$ with tensor $2$, $q_{2}$. In this case, we denote tensor $4$ as
$\left(T^{\mathcal{N}[q_{2}]}\right)^{s_{4}}_{Q_{4}}$.

An example for a TTN with $N=7$ sites and $L=6$ auxiliary indices (links) is
shown in Fig. 1. The property distinguishing a TTN from a general tensor
network is that the graph of a TTN is _loop-free_ , i.e., to move from one
site to any other there is only one unique path along the links. This also
implies that by cutting any link, the tensor networks splits into two
disconnected segments. Therefore, at each link there is a notion of _left_ and
_right_ which is a first hint towards the capability of TTNs to access the
Schmidt decomposition and with it also the reduced density matrix as
demonstrated below. We also define the _leaves_ of the TTN as all tensors with
just a single link index. For convenience, we assume site $N$ to be a leave of
the TTN. Since TTNs are _loop-free_ , one can also define a measure of
distance $d_{ij}$ between two sites $i$ and $j$ given by the number of links
one has to traverse to move from site $i$ to site $j$.

###  2.2 Tensor Gauge and Orthogonality Center

The representation of a quantum state as a tensor network is highly non-
unique. This _gauge degree of freedom_ can be used to obtain useful
representations of the same quantum state as a TTN with certain properties,
which can speed up calculations dramatically. As shown in Fig. 2, at each link
one can insert an identity $\mathbb{1}=G\cdot G^{-1}$ for any invertible
matrix $G$. By absorbing $G$ into one tensor and $G^{-1}$ into the other, a
different representation of the same state is reached. In this part, we make
use of this gauge degree of freedom to define an orthogonality center of the
TTN.
A tensor $T^{s_{i}}_{Q_{i}}$ can be orthogonalized towards one of its
neighbors with which it shares link $q_{k}$ as follows:

  * •

Reshape $T^{s_{i}}_{Q_{i}}$ into a matrix $T_{(s_{i},Q_{i}\setminus
q_{k}),(q_{k})}$ with rows $(s_{i},Q_{i}\setminus q_{k})$ and column
$(q_{k})$.

  * •

Perform an SVD (a QR decomposition is faster): $T_{(s_{i},Q_{i}\setminus
q_{k}),(q_{k})}=\sum_{\alpha}U_{(s_{i},Q_{i}\setminus q_{k}),(\alpha)}\cdot
S_{\alpha}\cdot V^{\dagger}_{(\alpha),(q_{k})}$

  * •

Keep $U_{(s_{i},Q_{i}\setminus q_{k}),(\alpha)}$ as the new local tensor on
site $i$ and absorb $S\cdot V^{\dagger}$ into the corresponding neighbor by
multiplying $S\cdot V^{\dagger}$ onto it (formally also relabel $\alpha\to
q_{k}$).

The SVD as well as the QR decomposition guarantees that the new site tensor
has the property (see Fig. 3)

$\left(U^{\dagger}\cdot
U\right)_{(\alpha),(\alpha^{\prime})}=\sum_{s_{i},Q_{i}\setminus
q_{k}}(U^{\dagger})_{(\alpha),(s_{i},Q_{i}\setminus
q_{k})}U_{(s_{i},Q_{i}\setminus
q_{k}),(\alpha^{\prime})}=\delta_{(\alpha),(\alpha^{\prime})}.$

For tensors orthogonalized towards their neighbor along link $q_{k}$ we
introduce the notation $\left(T^{\mathcal{N}[q_{k}]}\right)^{s_{i}}_{Q_{i}}$
(see Fig. 3 bottom).
As already mentioned, in TTNs there is a unique path between any two tensors.
Therefore, by orthogonalizing a tensor towards one of its neighbors, we also
orthogonalize it towards all other tensors, which can be reached via this
neighbor. For example, to orthogonalize tensors $1$, $2$ and $6$ in Fig. 1
towards tensor $3$, we orthogonalize all of them towards tensor $4$ with the
procedure described above.
Next, let us introduce orthogonality centers. Site $i$ is an orthogonality
center with tensor $C^{s_{i}}_{Q_{i}}$ if all tensors of all other sites are
orthogonalized towards site $i$. To obtain such an orthogonality center, we
can use the following algorithm:

  1. 1.

Find the maximum distance $d_{max}$ between site $i$ and any other site in the
TTN.

  2. 2.

Initialize $d=d_{max}$ and perform the following steps until $d=0$

     * •

Orthogonalize all sites $j$ that are at distance $d$ from site $i$ towards
site $i$, i.e., towards the single neighbor on the path from $j$ to $i$.

     * •

Reduce $d$ by one $d\to d-1$.

For example, to orthogonalize the TTN shown in Fig. 1 towards site $4$, we
first orthogonalize sites $3$ and $7$ towards site $5$ and then sites $1$,
$2$, $6$ and $5$ towards site $4$.

Figure 4: If the orthogonality center of the TTN depicted in Fig. 1 is placed
on site $4$, the center tensor $C$ has four link indices $q_{1}\cdots q_{4}$.
Each of these links corresponds to one of four mutually orthogonal set of
states $\ket{q^{[4]}_{1}}\cdots\ket{q^{[4]}_{4}}$. This orthogonality is a
direct result of the orthogonality property of the $U$-matrices of the SVD
used on all sites except site 4 (see also Fig. 3). The inset is a reminder of
the TTN used in this section.

The wave function of a TTN with orthogonality center $C^{s_{i}}_{Q_{i}}$ can
be written as:

$\displaystyle{\ket{\psi}}$ $\displaystyle=\sum_{\stackrel{{\scriptstyle
q_{1},q_{2},\cdots,q_{r_{i}}\in
Q_{i}}}{{s_{i}}}}C^{s_{i}}_{q_{1},q_{2},\cdots,q_{r_{i}}}\ket{s_{i}}\ket{q^{[i]}_{1}}\ket{q^{[i]}_{2}}\cdots\ket{q^{[i]}_{r_{i}}}$
$\displaystyle\ket{q^{[i]}_{k}}$ $\displaystyle=\sum_{s_{1}\cdots s_{r}\in
S^{i}_{q_{k}}}\left(T^{s_{1}}_{Q_{1}}\cdot T^{s_{2}}_{Q_{2}}\cdots
T^{s_{r}}_{Q_{r}}\right)_{q_{k}}\ket{s_{1}\cdots s_{r}}$ $\displaystyle\langle
q^{[i]}_{k}|q^{[i]^{\prime}}_{k}\rangle$
$\displaystyle=\delta_{q_{k},q_{k}^{\prime}}.$ (3)

Here, $S^{i}_{q_{k}}$ is the segment of the tensor network that is obtained by
cutting index $q_{k}$ and which does not contain site $i$. The states
$\ket{q^{[i]}_{k}}$ form an orthogonal basis and are defined in Fig. 4 for the
TTN of Fig. 1 with orthogonality center on site $i=4$.
Orthogonality centers allow to easily calculate local observables acting on
the orthogonality center. For example, the expectation value of the operator
$\hat{A}=\sum_{s_{i},s_{i}^{\prime}}A^{s_{i}^{\prime},s_{i}}\ket{s_{i}^{\prime}}\bra{s_{i}}$
acting non-trivially only on site $i$, reduces to:

$\langle\psi|\hat{A}|\psi\rangle=\sum_{\begin{subarray}{c}q_{1}q_{2}\cdots
q_{r_{i}}\in Q_{i}\\\
s_{i},s_{i}^{\prime}\end{subarray}}\bar{C}^{s_{i}^{\prime}}_{q_{1}q_{2}\cdots
q_{r_{i}}}\cdot A^{s_{i}^{\prime}s_{i}}\cdot C^{s_{i}}_{q_{1}q_{2},\cdots
q_{r_{i}}},$ (4)

where the bar denotes complex conjugation. Orthogonality centers hence reduce
the costly contraction over the whole tensor network, to a simple contraction
over the center tensor $C$ only.

###  2.3 Truncation of TTNs

A TTN with an orthogonality center allows to calculate any Schmidt
decomposition of the quantum state with respect to the two parts of the system
defined by cutting any of the links of the orthogonality center. To do so, we
reshape the center tensor $C$ into a matrix with physical index $s_{i}$ and
one of the links $q_{k}$, combined into the row index and all other indices
into the column indices, i.e.,
$C^{s_{i}}_{q_{1},q_{2},\cdots,q_{r_{i}}}=C_{(s_{i},q_{k}),(Q_{i}\setminus
q_{k})}$. The Schmidt decomposition then follows from an SVD of this matrix:

$\displaystyle C_{(s_{i}q_{k}),(Q_{i}\setminus
q_{k})}=\sum_{\alpha}U_{(s_{i}q_{k}),(\alpha)}\cdot
S_{\alpha}\cdot(V)^{\dagger}_{(\alpha),(Q_{i}\setminus q_{k})}$
$\displaystyle\Rightarrow{\ket{\psi}}$
$\displaystyle=\sum_{\alpha}S_{\alpha}\underbrace{\sum_{s_{i}q_{k}}U_{(s_{i}q_{k}),(\alpha)}\ket{q^{[i]}_{k}}\ket{s_{i}}}_{\ket{L}_{\alpha}}\cdot\underbrace{\sum_{Q_{i}\setminus
q_{k}}(V)^{\dagger}_{(\alpha),(Q_{i}\setminus q_{k})}\bigotimes_{l\neq
k}\ket{q^{[i]}_{l}}}_{\ket{R}_{\alpha}}$
$\displaystyle=\sum_{\alpha}S_{\alpha}\ket{L}_{\alpha}\ket{R}_{\alpha}.$ (5)

Note that if one is interested solely in the truncation of a given
orthogonality center and keeping it at the same site, an even more efficient
approach would be to perform an SVD on the matrix
$C_{(q_{k}),(s_{i}Q_{i}\setminus q_{k})}$. Again, the orthogonality of the
states $\ket{q_{k}}$ and the orthogonality of the $U$ and $V^{\dagger}$
matrices guarantee that the left and right vectors also form an orthogonal
basis and hence Eq. 5 is a true Schmidt decomposition. Note that this Schmidt
decomposition separates all sites in segment $S^{i}_{q_{k}}$ as well as site
$i$ to the rest of the lattice. From there it is straightforward to calculate
the reduced density matrix for one of these two subsystems and approximate
states by keeping only the largest eigenvalues in the spirit of DMRG.

##  3 TDVP equations for Tree Tensor Networks

Figure 5: Definition of the vertical subspace for site $4$ in Eq. 9 with end
point site $7$. The space of all tensors defining the kernel of the map from
the tangent space to the physical hilbert space, is spanned by a matrix
$X_{q_{l}^{\prime}q_{l}}$ (yellow squares) for each link of the TTN. The
physical index $s_{4}$ is labeled separately to distinguish it from the bond
indices. The inset is a reminder of the TTN used in this section.

In this section, we generalize the derivation of the tangent space projector
presented for MPS in Ref. [29] to general TTNs. While the overall approach is
very similar to the derivation for MPS, the lack of a clear start and end
point in the tensor network geometry will make the derivation and the
subsequent integration of the equations quite different from standard MPS.
TDVP amounts to the solution of the Schrödinger equation in the space spanned
by the tensor network without ever leaving this manifold (at least in its
single-site variant). In TDVP, one solves a modified Schrödinger equation by
projecting its right-hand side onto the so-called tangent space:

$\frac{d\ket{\psi[T]}}{dt}=-i\mathcal{P}_{T_{\ket{\psi[T]}}}H\ket{\psi[T]}.$
(6)

In the following, we want to find a representation of the tangent space
projection operator $\mathcal{P}_{T_{\ket{\psi[T]}}}$, which not only depends
on the current state $\ket{\psi[T]}$ but importantly also on the structure of
the TTN.

###  3.1 Tangent Space Projector

Any element of the tangent space $\ket{\Theta[B]}$ is parametrized by a set of
tensors $B^{s_{i}}_{Q_{i}}$:

$\ket{\Theta[B]}=\sum_{i=1}^{N}B^{s_{i}}_{Q_{i}}\frac{d\ket{\psi[T]}}{dT^{s_{i}}_{Q_{i}}}.$
(7)

Importantly, for each summand we use the representation of the state
$\ket{\psi[T]}$ in which site $i$ is the orthogonality center, i.e., all
tensors $T^{s_{j}}_{Q_{j}}$ are orthogonalized towards site $i$ such that

$\ket{\Theta[B]}=\sum_{i=1}^{N}\sum_{\stackrel{{\scriptstyle
Q_{i}}}{{s_{i}}}}B^{s_{i}}_{Q_{i}}\ket{s_{i}}\ket{q^{[i]}_{1}\cdots
q^{[i]}_{r_{i}}}$ (8)

The gauge degree of freedom in the TTN, reflects itself in the tangent space
that not all linearly independent choices of $B^{s_{i}}_{Q_{i}}$ result in
different tangent vectors. Ref. [29] solves this problem by first defining the
so-called vertical subspace, i.e., all tensors $B^{s_{i}}_{Q_{i}}$ that give
the zero-state $\ket{\Theta[B]}=0$ and hence define the kernel of the map from
the tensors to the physical Hilbert space. Then, imposing a gauge
prescription, they fix this kernel to a single element which guarantees that
the resulting parametrization is unique.
In order to arrive at a result that resembles the MPS algorithm, we first need
to define a fixed end point of the TTN with the restriction that it should be
a leave. Note however that any site of the tensor network can be used as end
point. Without loss of generality, we choose site $N$ as end point. The
vertical subspace, i.e., all tensors $B^{s_{i}}_{Q_{i}}$ for which
$\ket{\Theta[B]}=0$ can then be parametrized by matrices
$X_{q_{k}^{\prime}q_{k}}$ such that:

$\displaystyle
B^{s_{i}}_{Q_{i}}=\sum_{l=1}^{r_{i}}\sum_{q_{l}^{\prime}}\left(T^{\mathcal{N}[q_{l}]}\right)^{s_{i}}_{q_{1}\cdots
q_{l}^{\prime}\cdots q_{r_{i}}}X_{q_{l}^{\prime}q_{l}}\cdot\text{sgn}(q_{l}\to
N)$ $\displaystyle\text{sgn}(q_{l}\to N)=\begin{cases}\mbox{1,}&\mbox{if
}q_{l}\text{ points towards $N$}\\\ \mbox{-1,}&\mbox{otherwise}\end{cases}.$
(9)

$(T^{\mathcal{N}[q_{l}]})^{s_{i}}_{Q_{i}}$ is the unique tensor of the state
$\ket{\psi[T]}$ with site $i$ orthogonalized towards the neighbor on the other
end of the link $q_{l}$. This definition of the vertical subspace is depicted
in Fig. 5 for the tensor $B^{s_{4}}_{Q_{4}}$.
The factor $\text{sgn}(q_{l}\to N)$ is $1$ if link $q_{l}$ points towards the
end point and $-1$ otherwise. This construction guarantees that for any choice
of $B^{s_{i}}_{Q_{i}}$ in the vertical subspace, $\ket{\Theta[B]}=0$, because
the single term with positive sign $(q_{l}\to N)$ is exactly canceled by one
negative term of its neighbor (since there $(q_{l}\nrightarrow N)$). Note that
this definition of the vertical subspace reduces in the case of MPS to the
definition used in Ref. [29] if the right-most site of the MPS is chosen as
the end point.
To uniquely specify the kernel, we impose the following matrix-valued (with
indices $q_{k}$ and $q_{k}^{\prime}$) gauge fixing condition for the
$B$-tensors of the tangent space:

$\displaystyle\sum_{\stackrel{{\scriptstyle Q_{i}\setminus
q_{k}}}{{s_{i}}}}\bar{B}^{s_{i}}_{q_{1}\cdots q_{k}^{\prime}\cdots
q_{r_{i}}}\cdot\left(T^{\mathcal{N}[q_{k}]}\right)^{s_{i}}_{q_{1}\cdots
q_{k}\cdots q_{r_{i}}}=0\quad\forall\quad i\neq N.$ (10)

Again, the bar denotes complex conjugation. Above, $q_{k}$ is the single index
pointing towards the end point $N$. These are $N-1$ matrix-valued constraints,
for the $X$-matrices living on $L=N-1$ indices. This implies that no ambiguity
is left in the definition of the kernel, if we choose $B$-tensors according to
Eq. 10.
It also guarantees that the overlap between two tangent vectors reduces to a
contraction over local tensors only:

$\langle\Theta[B^{\prime}]|\Theta[B]\rangle=\sum_{i=1}^{N}\sum_{\stackrel{{\scriptstyle
Q_{i}}}{{s_{i}}}}\bar{B}^{s_{i}}_{Q_{i}}\cdot B^{s_{i}}_{Q_{i}}.$ (11)

Similar to MPS, we can now reformulate the projection problem of an arbitrary
state $\ket{\Xi}$ onto the tangent space
$\ket{\Theta[B]}=\mathcal{P}_{T_{\ket{\psi[T]}}}\ket{\Xi}$ as a minimization
problem:

$\min_{B}||\ket{\Theta[B]}-\ket{\Xi}||^{2},$ (12)

under the constraints given by Eq. 10. With Eq. 8 and using a Lagrange
multipliers $\lambda^{[i]}_{q_{k}q_{k}^{\prime}}$ to account for the
constraints, the minimization can be reformulated as:

$\displaystyle\min_{B}\Bigg{[}\sum_{i=1}^{N}\sum_{\stackrel{{\scriptstyle
Q_{i}}}{{s_{i}}}}\left(\bar{B}^{s_{i}}_{Q_{i}}\cdot
B^{s_{i}}_{Q_{i}}-\bar{B}^{s_{i}}_{Q_{i}}\cdot
F^{s_{i}}_{Q_{i}}-\bar{F}^{s_{i}}_{Q_{i}}\cdot B^{s_{i}}_{Q_{i}}\right)$
$\displaystyle-\sum_{i=1}^{N-1}\sum_{q_{k}q_{k}^{\prime}}\lambda^{[i]}_{q_{k}q_{k}^{\prime}}\sum_{\stackrel{{\scriptstyle
Q_{i}\setminus q_{k}}}{{s_{i}}}}\bar{B}^{s_{i}}_{q_{1}\cdots
q_{k}^{\prime}\cdots
q_{r_{i}}}\cdot\left(T^{\mathcal{N}[q_{k}]}\right)^{s_{i}}_{q_{1}\cdots
q_{k}\cdots q_{r_{i}}}\Bigg{]}$ (13)

with $F^{s_{i}}_{Q_{i}}=\langle s_{i}q_{1}\cdots q_{r_{i}}\ket{\Xi}$. The
solution to this minimization problem can be found by setting the derivative
with respect to $\bar{B}^{s_{i}}_{Q_{i}}$ as well as
$\lambda^{[i]}_{q_{k}q_{k}^{\prime}}$ to zero. Using some algebra we find the
minimum for all sites $i\neq N$:

$\displaystyle B^{s_{i}}_{Q_{i}}=$ $\displaystyle
F^{s_{i}}_{Q_{i}}-\sum_{\stackrel{{\scriptstyle Q_{i}^{\prime\prime}\setminus
q_{k}^{\prime\prime},q_{k}^{\prime}}}{{t}}}\left(T^{\mathcal{N}[q_{k}]}\right)^{s_{i}}_{q_{1}\cdots
q_{k}^{\prime}\cdots q_{r_{i}}}$
$\displaystyle\cdot\left(T^{\mathcal{N}[q_{k}]}\right)^{t}_{q_{1}^{\prime\prime}\cdots
q_{k}^{\prime}\cdots q_{r_{i}}^{\prime\prime}}\cdot
F^{t}_{q_{1}^{\prime\prime}\cdots q_{k}\cdots q_{r_{i}}^{\prime\prime}},$ (14)

while for $i=N$ it is just $B^{s_{N}}_{Q_{N}}=F^{s_{N}}_{Q_{N}}$. This allows
us to obtain a representation of the tangent space projector
$\ket{\Theta[B]}=\mathcal{P}_{T_{\ket{\psi[T]}}}\ket{\Xi}$ as:

$\displaystyle\mathcal{P}_{T_{\ket{\psi[T]}}}$
$\displaystyle=\sum_{i=1}^{N}\mathbb{1}_{s_{i}}\otimes\sum_{Q_{i}}\ket{q^{[i]}_{1}\cdots
q^{[i]}_{r_{i}}}\bra{q^{[i]}_{1}\cdots q^{[i]}_{r_{i}}}$
$\displaystyle-\sum_{<i,j>_{q_{k}}}\sum_{q_{k}q_{k}^{\prime}}\ket{q^{[j]\prime}_{k}}\bra{q^{[j]\prime}_{k}}\otimes\ket{q^{[i]}_{k}}\bra{q^{[i]}_{k}},$
(15)

Figure 6: Definition of the states used in the projection operator onto the
link $q_{k}$ defined in the second line of Eq. 15 for link $q_{1}$ connecting
sites $i=4$ and $j=5$.

Where $\sum_{<i,j>_{q_{k}}}$ denotes a sum over all nearest neighbors $i$ and
$j$ with the corresponding index $q_{k}$ connecting these two sites. The
graphical representation of the states in the second line of the tangent space
projector for the bond connecting sites $i=4$ and $j=5$ is shown in Fig. 6.
Formally, this result resembles the projection operator obtained for MPS [29].
The first term with positive sign corresponds to the forward time propagation
of the site tensor. The second term on the other hand is the evolution
backwards in time of the bonds between two site tensors and is a direct
consequence of the gauge fixing of the tangent vectors used in Eq. 10.

###  3.2 Single-Site TDVP

With the representation of the projection operator in Eq. 15, we can go back
to the projected time dependent Schrödinger equation (Eq. 6) and integrate
each term one by one using Trotter breakups [39]. First, let us discuss a
first-order update, which can later easily be modified to perform a second
order integration. Since each term in the projection operator keeps all but
one tensor fixed, the integration can be performed locally. Therefore, we
define effective Hamiltonians for the sites $i$ and for the links $q_{k}$:

$\displaystyle H_{(s_{i}Q_{i}),(s_{i}^{\prime}Q_{i}^{\prime})}=\langle
s_{i}q^{[i]}_{1}\cdots q^{[i]}_{r_{i}}|H|s_{i}^{\prime}q^{[i]\prime}_{1}\cdots
q^{[i]\prime}_{r_{i}}\rangle$ (16a) $\displaystyle
K_{(q^{[i]}_{k}q^{[j]}_{k})(q^{[i]\prime}_{k}q^{[j]\prime}_{k})}=\langle
q^{[i]}_{k}q^{[j]}_{k}|H|q^{[i]\prime}_{k}q^{[j]\prime}_{k}\rangle$ (16b)

and solve equations of the form:

$\displaystyle\mathbf{\dot{A}}$ $\displaystyle=\pm
iH^{\text{eff}}\cdot\mathbf{A}$ $\displaystyle\mathbf{A}(t+{\Delta t})$
$\displaystyle=e^{\pm iH^{\text{eff}}{\Delta t}}\mathbf{A}(t),$ (17)

where $\mathbf{A}$ is either a site-tensor or a link tensor and
$H^{\text{eff}}$ either $H_{(s_{i}Q_{i}),(s_{i}^{\prime}Q_{i}^{\prime})}$
(negative sign) or
$K_{(q^{[i]}_{k}q^{[j]}_{k})(q^{[i]\prime}_{k}q^{[j]\prime}_{k})}$ (positive
sign). In matrix form, the exponential of these effective Hamiltonians can be
efficiently calculated using Krylov exponentiation.
A full TDVP step is then given by a series of $N-1$ local updates of a site
tensor and the corresponding link tensor connecting the site to the end point
as shown below. The single local update on site $i$ and link $q_{k}$ is

  * •

Orthogonalize the TTN such that site $i$ is the orthogonality center.

  * •

Calculate the one-site effective Hamiltonian
$H^{\text{eff}}=H_{(s_{i}Q_{i}),(s_{i}^{\prime}Q_{i}^{\prime})}$ (Eq. 16a) and
forward time evolve (negative sign) according to Eq. 17 with
$\mathbf{A}=C^{s_{i}}_{Q_{i}}$. If site $i$ is the chosen end point, stop
here; otherwise continue.

  * •

Reshape the time evolved tensor into a matrix
$C^{s_{i}}_{Q_{i}}=C_{(s_{i}Q_{i}\setminus q_{k}),(q_{k})}$ and perform an SVD
(QR-decomposition suffices) $C_{(s_{i}Q_{i}\setminus
q_{k}),(q_{k})}=\sum_{q^{[i]}_{k}}U_{(s_{i}Q_{i}\setminus
q_{k}),(q^{[i]}_{k})}\cdot\underbrace{S_{q^{[i]}_{k}}\cdot(V^{\dagger})_{(q^{[i]}_{k}),(q_{k})}}_{L_{q^{[i]}_{k}q_{k}}}$.
As usual, take the $U$-tensor as new tensor on site $i$.

  * •

Calculate the link effective Hamiltonian
$H^{\text{eff}}=K_{(q^{[i]}_{k}q^{[j]}_{k})(q^{[i]\prime}_{k}q^{[j]\prime}_{k})}$
(Eq. 16b) for link $q_{k}$. To do so, use the time evolved tensor obtained in
the previous step for site $i$. Then evolve tensor
$\mathbf{A}=L_{q^{[i]}_{k}q_{k}}\equiv L_{q^{[i]}_{k}q^{[j]}_{k}}$ from the
previous step backwards in time (positive sign) according to Eq. 17. Finally,
absorb the $C$-tensor onto the neighbor of site $i$ along $q_{k}$ by
multiplying it onto its site tensor.

A full TDVP time step can then achieved by the following sweeping procedure:

  1. 1.

Choose a start and an end point; initialize site $i$ as the chosen start
point.

  2. 2.

Perform the following steps until $i$ is the chosen end point:

     * •

Find the link $q_{k}\in Q_{i}$ that connects site $i$ to the end point.

     * •

If any tensor attached to the other links $Q_{i}\setminus q_{k}$ has not been
updated, choose one of these links and choose one of the leaves attached to
the corresponding segment of the TTN as new site $i$.

     * •

Otherwise, perform a local update on site $i$ as described above and choose
the neighbor of site $i$ along link $q_{k}$ as new site $i$.

  3. 3.

Perform one last local update for the endpoint $i=N$ as described above.

A depiction of the sweeping order for the TTN in Fig. 1 is shown in Fig. 7.
The procedure described above defines a first-order time step. A second-order
method can easily be obtained by performing the first order time step with
$\frac{{\Delta t}}{2}$ and then simply perform the exact same steps in reverse
order corresponding to repeated second order Trotter breakups
$e^{\tau(A+B+C)}=e^{\frac{\tau}{2}C}e^{\frac{\tau}{2}B}e^{\tau
A}e^{\frac{\tau}{2}B}e^{\frac{\tau}{2}C}$ used on Eq. 15. Importantly, this
means that during the local update, the link update has to be performed
_before_ the site update (see also caption of Fig. 7).
Note that for a given TTN, there can be several versions of this algorithm
depending on the sequence of chosen indices in step 2. Very often though, the
TTN structure itself defines some natural order when to time evolve which
sites, as we will see in the next section for the FTPS tensor network.

Figure 7: Update sequence to perform a first-order single-site TDVP time step
from time $t$ to $t+{\Delta t}$ for the TTN shown in Fig. 1. Start point is
site $1$ and end point site $7$. Yellow denotes tensor that are updated in the
current step. Red and blue tensors indicate whether this tensor is taken at
time $t+{\Delta t}$ (red) or $t$ (blue). Triangles indicate the
orthogonalization of each tensor. Updates on site-tensors are in forward
direction (negative sign), while updates on bond-tensors are backwards time
evolutions (positive sign in Eq 17) For a second order update, first perform
all steps $(1)\to(13)$ with time step $\frac{dt}{2}$ in the order shown and
then reapply them in the reverse order $(13)\to(1)$, again with time step
$\frac{dt}{2}$.

###  3.3 Two-Site TDVP

Figure 8: Update sequence to perform a first-order two-site TDVP time step
from time $t$ to $t+{\Delta t}$ for the TTN shown in Fig. 1. Start point is
site $1$ and end point site $7$. Yellow denotes tensors that are updated in
the current step. Red and blue tensors indicate whether this tensor is taken
at time $t+{\Delta t}$ (red) or $t$ (blue). Triangles indicate the
orthogonalization of each tensor. Updates on two-site tensors are in forward
direction (negative sign), while updates on a single site are backwards time
evolutions (positive sign in Eq. 17) For a second-order update, first perform
all steps $(1)\to(11)$ with time step $\frac{dt}{2}$ in the order shown and
then reapply them in the reverse order $(11)\to(1)$, again with time step
$\frac{dt}{2}$.

It is also straightforward to generalize the single-site TDVP approach
presented above to a two-site TDVP integration scheme which allows to
dynamically adapt the necessary bond dimensions. To do so, we need to define
the two-site effective Hamiltonian $H^{\text{2-site}}$ for two sites $i$ and
$j$ connected by the index $q_{l}$:

$\displaystyle H^{\text{2-site}}$ $\displaystyle=\langle
s_{i}s_{j}Q^{\text{red}}_{i}Q^{\text{red}}_{j}|H|s_{i}^{\prime}s_{j}^{\prime}Q^{\text{red}}_{i}Q^{\text{red}}_{j}\rangle$
$\displaystyle\ket{Q^{\text{red}}_{i}}$ $\displaystyle=\bigotimes_{q_{n}\in
Q_{i}\setminus q_{l}}\ket{q^{[i]}_{n}}$
$\displaystyle\ket{Q^{\text{red}}_{j}}$ $\displaystyle=\bigotimes_{q_{n}\in
Q_{j}\setminus q_{l}}\ket{q^{[j]}_{n}}.$ (18)

With this, only small modifications to the algorithm presented above are
necessary. The single update for sites $i$ and $j$ sharing link $q_{l}$
becomes:

  * •

Orthogonalize the TTN such that site $i$ is the orthogonality center.

  * •

Calculate the two-site effective Hamiltonian $H^{\text{2-site}}$ according to
Eq. 18. and forward time evolve (negative sign) with
$\mathbf{A}=\sum_{q_{l}}C^{s_{i}}_{Q_{i}}T^{s_{j}}_{Q_{j}}$.

  * •

Reshape the time evolved tensor into a matrix $A_{(s_{i}Q_{i}\setminus
q_{l}),(s_{j}Q_{j}\setminus q_{l})}$ and perform an SVD
$A_{(s_{i}Q_{i}\setminus q_{l}),(s_{j}Q_{j}\setminus
q_{l})}=\sum_{q_{l}}U_{(s_{i}Q_{i}\setminus
q_{l}),(q_{l})}\cdot\underbrace{S_{q_{l}}\cdot(V^{\dagger})_{(q_{l}),(s_{j}Q_{j}\setminus
q_{l})}}_{C^{s_{j}}_{Q_{j}}}$. In this step one can also truncate the smallest
Schmidt values. As usual, keep the $U$-tensor to update site $i$ and
$C^{s_{j}}_{Q_{j}}$ as site tensor on site $j$, shifting the orthogonality
center to site $j$. If site $j$ is the chosen end point, stop here; otherwise
continue.

  * •

Calculate the one-site effective Hamiltonian
$H^{\text{eff}}=H_{(s_{j}Q_{j}),(s_{j}^{\prime}Q_{j}^{\prime})}$ (Eq. 16a) for
site $j$. To do so, use the time-evolved tensor obtained in the previous step
for site $i$. Then evolve tensor $\mathbf{A}=C^{s_{j}}_{Q_{j}}$ backwards in
time (positive sign in Eq. 17).

A full two-site TDVP step can then be performed by:

  1. 1.

Choose a start and end point. Initialize site $i$ as the chosen start point.

  2. 2.

Perform the following steps until $i$ is the chosen end point:

     * •

Find the link $q_{k}$ and the corresponding neighbor $j$ that connects site
$i$ to the end point.

     * •

If any tensor attached to the other links $Q_{i}\setminus q_{k}$ has not been
updated, choose one of these links and choose one of the leaves attached to
the corresponding segment of the TTN as new site $i$.

     * •

Otherwise, perform a local update on site $i$ and $j$ as described above and
go to site $j$, i.e., $i\to j$.

A depiction of the necessary sweeping order of this two-site scheme is shown
in Fig. 8

##  4 TDVP for FTPS

An FTPS is a special TTN designed to efficiently encode states of multi-
orbital Anderson Impurity Models (AIMs). An AIM consists of an interacting
impurity coupled to a bath of free fermions with Hamiltonian

$\displaystyle H$
$\displaystyle={H_{\text{loc}}}+{H_{\text{bath}}}+{H_{\text{hyb}}}$
$\displaystyle{H_{\text{loc}}}~{}~{}$
$\displaystyle=\sum_{m\sigma}\epsilon_{m\sigma 0}n_{m\sigma
0}+H_{\text{int}}.$ $\displaystyle{H_{\text{bath}}}$
$\displaystyle=\sum_{m\sigma}\sum_{k}\epsilon_{m\sigma k}n_{m\sigma k}$
$\displaystyle{H_{\text{hyb}}}~{}$
$\displaystyle=\sum_{m\sigma}\sum_{k}V^{[k]}_{m\sigma}\left(c_{m\sigma
0}^{\dagger}c_{m\sigma k}+\text{h.c.}\right)$ (19)

$c_{m\sigma k}^{\dagger}$ ($c_{m\sigma k}$) creates (annihilates) an electron
in chain $m$ with spin $\sigma$ on site $k$, where $k=0$ denotes the impurity
site (see Fig. 9 (b)). $n_{m\sigma k}$ are the corresponding particle number
operators. $H_{\text{int}}$ is the interaction Hamiltonian that only couples
impurity degrees of freedom and for which we choose the Kanamori Interaction
[49, 18] without the spin-flip and pair-hopping terms parametrized by two
interaction strengths $U$ and $J$

$H_{\text{int}}=U\sum_{m}n_{m\uparrow 0}n_{m\downarrow
0}+(U-2J)\sum_{m^{\prime}>m\sigma}n_{m\sigma
0}n_{m^{\prime}\bar{\sigma}0}+(U-3J)\sum_{m^{\prime}>m\sigma}n_{m\sigma
0}n_{m^{\prime}\sigma 0}$ (20)

where $\bar{\sigma}$ is the opposite spin direction of $\sigma$. In the
following, we will use a combined index $l=(m\sigma)$ to denote the orbital
and spin-degrees of freedom.
For a single orbital, an FTPS reduces to a MPS, while for multiple orbitals it
has tensors with three link indices as depicted in Fig. 9 for a two-orbital
model. It consists of a single MPS-like chain for the bath tensors of each
orbital/spin and impurity tensors connecting the different chains. An FTPS for
a $N_{\text{orb}}$-orbital AIM has a total of $N_{C}=2N_{\text{orb}}$ chains.
For simplicity we assume that each chain has the same number of bath sites
$N_{b}$.

Figure 9: Graphical representation of a FTPS tensor network for a two orbital
model. For each orbital, we use two chains, one for each spin-species. _(a)_
one way to label the sites is just to numerate them in ascending order. _(b)_
a different way to label sites is to specify the chain (orbital $m$ and spin
$\sigma$) as well as an index (bath index $k$). This way to label sites
resembles the labels used for the operators of the Hamiltonian in Eq. 4.

According to the algorithm presented in the previous section, we first need to
choose a start and end point. We choose to start at the outermost bath site of
the first chain (site $4$ in Fig. 9 (a)) and the outermost bath site of the
last chain as end point (site $16$ in Fig. 9 (a)). To actually perform the
time evolution, we choose to employ a hybrid TDVP scheme using 2-site TDVP for
the bath tensors as well as for the bath-impurity link, and 1-site TDVP for
the impurity tensors itself and the corresponding impurity-impurity links. We
choose to use 1-site TDVP for the impurity links, since 2-site TDVP becomes
computationally expensive, since one would have to deal with tensors with four
link indices. This leads to the following algorithm for a single time step:

  1. 1.

For $l=1:N_{C}-1$ perform the following steps:

     * •

For $k=N_{b}:1$:

       * –

Perform a two-site step on sites $i=(l,k)$ and $j=(l,k-1)$ (see Fig. 9 for the
definition of the site-labeling).

     * •

Perform a one-site step on the impurity tensor $i=(l,0)$; $q_{k}$ connects
site $(l+1,0)$

  2. 2.

$l=N_{C}$, for $k=0:N_{b}-1$ perform the following steps:

     * •

Perform a two-site step on sites $i=(l,k)$ and $j=(l,k+1)$.

For the actual calculations, we apply the second order version of this
algorithm by using only the half time step and reapplying each step in reverse
order. Again, this also means that the order in the local updates changes.
Note that the backwards propagation during the two-site involving an impurity
site cancels with the subsequent forwards time evolution of the one-site step
on the same impurity tensor. Therefore, these two steps can be omitted.

$\displaystyle 0$$\displaystyle 5$$\displaystyle 10$$\displaystyle
15$$\displaystyle t$$\displaystyle 0.0$$\displaystyle
0.5$$\displaystyle\mathbf{\Re G^{>}(t)}$TEBDTDVP Figure 10: Comparison of the
impurity greater Greens function $G^{>}(t)$ for a two orbital model between
the TEBD time-evolution approach used in Ref. [18] and TDVP. The calculation
was performed for a spin- and orbital degenerate model using $9$ bath sites
per orbital and spin with parameters $\epsilon_{k}=-0.8,-0.6,\cdots,0.6,0.8$
and $V^{[k]}_{l}=0.1~{}\forall~{}k,l$. Therefore, only diagonal entries of the
Green’s function are non-zero and for the diagonals there is only one
independent function, i.e.,
$G^{>}_{l^{\prime},l}(t)=\delta_{l^{\prime}l}G^{>}(t)$. Interaction parameters
were $U=1$, $J=0.1$. The impurity on-site energy was chosen to obtain particle
hole symmetry, i.e., $\epsilon_{m\sigma 0}=-\frac{3U-5J}{2}$. The time step
for TEBD was ${\Delta t}=0.01$ and for TDVP ${\Delta t}=0.1$, since TDVP
generally allows to use larger time steps [50]. Truncated weight (sum of all
truncated Schmidt values) for TEBD was $10^{-12}$ without restricting the bond
dimension, and for TDVP it was $10^{-9}$ for all links except the impurity-
impurity links which were not truncated, but restricted to a maximal dimension
of 50.

As a first demonstration of this algorithm, let us compare the TDVP time
evolution to the TEBD-like approach used in Refs. [18, 19]. Therefore, we look
at the greater Greens function of the impurity defined by:

$G^{>}_{l^{\prime},l}(t)=\langle\psi_{0}|c_{l^{\prime}0}e^{-iHt}c_{l0}^{\dagger}|\psi_{0}\rangle
e^{iE_{0}t}.$ (21)

$\ket{\psi_{0}}$ is the ground state of Hamiltonian $H$ with ground state
energy $E_{0}$. For a degenerate two orbital model, Fig. 10 shows that the
TDVP time evolution indeed produces the correct result. In a recent
publication, the authors have shown, that for diagonal hybridizations, TDVP
has larger errors than TEBD for the bath geometry chosen here [50]. This means
that for such systems, the TEBD approach is most likely preferable over TDVP.
For more involved baths on the other hand, TEBD can become difficult to
formulate as discussed next.

One of the major advantages of TDVP is that it allows to perform the time
evolution for arbitrary couplings in the Hamiltonian between the sites, as
long as an _MPO_ with the same tensor network structure as the state can be
found. Eq. 4 is in fact not the most general AIM, since the bath only couples
_diagonally_ to its impurity. Often, one is also interested in so-called off-
diagonal hybridizations which can be encoded as hoppings from impurity $l$ to
a different bath $l^{\prime}$. Therefore, we can account for off-diagonal
hybridizations by replacing the hybridization terms in Eq. 4 with:

$\sum_{ll^{\prime}k}V^{[k]}_{ll^{\prime}}(c_{l0}^{\dagger}c_{l^{\prime}k}+\text{h.c}).$
(22)

It turns out that for each $k$, the matrix $V^{[k]}_{ll^{\prime}}$ can be
chosen as a lower-triangular matrix. This means, that for a spin-symmetric,
two-orbital model there are three free parameters for each value of $k$
(instead of two for the diagonal hybridization).

Figure 11: Comparison of impurity greater Greens function
$G_{ll^{\prime}}^{>}(t)$ for a two-orbital model between TDVP and the exact
solution. In each segment, the top panel shows the Green’s function itself,
while the bottom panel shows the absolute value of the difference
$\epsilon(t)=|G^{>}_{\text{exact}}(t)-G^{>}_{\text{TDVP}}(t)|$, i.e., the
numerical error. We obtained the exact solution from diagonalization of the
hopping matrix at $U=J=0$. The calculation was performed for a spin-degenerate
model using $9$ bath sites per orbital and spin. We allowed off-diagonal
hopping terms only between the orbital degrees of freedom, i.e.,
$V^{[k]}_{(m\sigma)(m^{\prime}\sigma^{\prime})}=\delta_{\sigma\sigma^{\prime}}V^{[k]}_{mm^{\prime}}$
and therefore
$G_{(m\sigma)(m^{\prime}\sigma^{\prime})}=\delta_{\sigma\sigma^{\prime}}G_{mm^{\prime}}$.
The parameters were $\epsilon_{k}=-0.8,-0.6,\cdots,0.6,0.8$ for all orbitals,
diagonal hybridizations $V^{[k]}_{mm}=0.1~{}\forall~{}k,l$ and off-diagonal
hybridizations $V^{[k]}_{mm^{\prime}}=0.05~{}\forall~{}k$ for
$m=2,m^{\prime}=1$. These off-diagonal terms correspond to a hopping processes
from the impurity of orbital $2$ to the bath of orbital $1$. The TDVP time
step was chosen $\Delta t=0.1$ and on-site energies were $\epsilon_{m\sigma
0}=0$. Note that the off-diagonal hybridizations break the orbital degeneracy,
albeit for the parameters chosen only slightly and the differences between the
two orbitals are barely visible. Truncated weight (sum of all truncated
Schmidt values) during DMRG and the time evolution was $10^{-9}$, except
during the time evolution of the impurity impurity links where no truncation
was performed. During DMRG as well as time evolution, the impurity-impurity
links were restricted to 140.

$\displaystyle 10^{-1}$$\displaystyle 10^{0}$$\displaystyle\Delta
t$$\displaystyle 10^{-5}$$\displaystyle 10^{-4}$$\displaystyle
10^{-3}$$\displaystyle 10^{-2}$max($\displaystyle\epsilon(t)$)TDVPFit
$\displaystyle\sim(\Delta t)^{2}$

2070140200$\displaystyle m_{Imp-Imp}$$\displaystyle 10^{-4}$$\displaystyle
10^{-3}$max($\displaystyle\epsilon(t)$)TDVP

Figure 12:  Error as a function of control parameters for the same AIM used in
Fig. 11. We plot the maximum value of the error obtained in all four Green’s
functions $max(\epsilon(t))$.
_Left:_ As function of step size $\Delta t$, the error shows the expected
scaling $\sim(\Delta t)^{2}$ for larger values of $\Delta t$. The deviations
for smaller values can be explained by the other sources of error, like the
truncation during time evolution and a not perfect representation of the
ground state. Additionally, we found TDVP in the star geometry to be quite
sensitive to a too small time step in combination with a too large truncation.
The parameters used in the truncation of the tensor network were exactly the
same as discussed above.
_Right:_ Error as a function of impurity-impurity bond dimension. All other
parameters were the same as above.

As a second demonstration of the TDVP approach for FTPS we calculate the
$2\times 2$ matrix of the greater Green’s function of such a spin symmetric
two-orbital model. Since the TEBD approach we compared with in Fig. 10 is
difficult to generalize to such off-diagonal hybridizations, we perform the
calculation in the non-interacting case $U=J=0$ and note that for tensor
network based approaches this is a highly non-trivial situation. This is
because the bipartitions defined by the links of the FTPS structure have non-
trivial entanglement also for non-interacting systems, and the off-diagonal
hoppings $V^{[k]}_{ll^{\prime}}$ for $l\neq l^{\prime}$ introduce entanglement
_between_ the orbitals, i.e., non-trivial links between the impurities. The
results of such a comparison can be seen in Fig. 11. Having access to the
exact solution, we also plot the difference between the exact and numerical
Green’s functions in the bottom panels. Again we find very good agreement
between TDVP and the reference calculations.
Finally let us demonstrate that the results indeed converge with respect to
the control parameters. The left plot of Fig. 12 shows the scaling of the
error as a function of $\Delta t$ and we indeed observe the expected
$~{}\sim{\Delta t}^{2}$ behavior at larger values of ${\Delta t}$. The
deviation of this behavior at smaller ${\Delta t}$ can be understood from the
additional errors due to the truncations of the tensor network in the ground
state as well as during time evolution. Additionally, we frequently observed
that when TDVP is used in the star-geometry representation of the bath (with
long-range couplings $V_{m\sigma}$), a good balance between truncation and
${\Delta t}$ is necessary. Surprisingly we found that it is often advantageous
to use rather large time steps compared to what one would use in TEBD
calculations. In the right plot of Fig. 12 we show the convergence of the
error as a function of dimension of the impurity-impurity links which is
usually the bottle-neck of FTPS calculations as these links need to transport
the entanglement between the different orbitals $l$. Also here, we observe
convergence with the control parameter, showing that TDVP indeed can be
efficiently used to account for off-diagonal hybridizations.

##  5 Conclusion

We presented a generalization of the Time Dependent Variational Principle
(TDVP) to general loop-free tensor networks (TTNs). The major advantage of
TDVP over the commonly used TEBD approach is that the latter is often
difficult to implement if long-range couplings are present in the Hamiltonian.
TDVP on the other hand allows to perform the time evolution (either in
imaginary- or real-time) for any Hamiltonian for which a representation in the
same TTN structure can be found, which is often possible for long-range
couplings. Using a similar derivation as in Ref. [29], we were able to find
the projection operator onto the tangent space for any TTN - the central
object in TDVP. Integrating the terms in the tangent space projector one after
the other, equivalent to a Suzuki Trotter breakup, we were able to formulate
TDVP in its single-site as well as two-site variant. We then applied TDVP to
the FTPS tensor network which is a TTN especially suited for multi-orbital
Anderson impurity models. For FTPS, TDVP is particularly appealing if the
hybridizations with the bath are off-diagonal. In DMFT calculations, off-
diagonal hybridizations are of significance to account for spin-orbit coupling
effects as well as distortions of the crystal lattice. We verified the TDVP
approach by comparing first to TEBD using a diagonal bath including
interactions, and second to the exact solution in the non-interacting case for
an off-diagonal bath.
When finalizing this manuscript we became aware of an independent publication
by Kohn et al. [51], describing the TDVP applied to a TTN for periodic
boundary conditions in a one-dimensional system.

The authors would like to thank Florian Maislinger, Hans Gerd Evertz and Jutho
Haegeman for fruitful discussions. This work was supported by the Austrian
Science Fund (FWF) through the START program Y746, as well as by NAWI Graz.

## References

  * [1] S. R. White,  _Density matrix formulation for quantum renormalization groups_ ,  Phys. Rev. Lett. 69, 2863 (1992),  10.1103/PhysRevLett.69.2863.
  * [2] U. Schollwöck,  _The density-matrix renormalization group in the age of matrix product states_ ,  Ann. Phys. 326, 96 (2011),  10.1016/j.aop.2010.09.012.
  * [3] S. Östlund and S. Rommer,  _Thermodynamic limit of density matrix renormalization_ ,  Phys. Rev. Lett. 75, 3537 (1995),  10.1103/PhysRevLett.75.3537.
  * [4] J. Dukelsky, M. A. Martin-Delgado, T. Nishino and G. Sierra,  _Equivalence of the variational matrix product method and the density matrix renormalization group applied to spin chains_ ,  EPL (Europhysics Letters) 43(4), 457 (1998).
  * [5] F. Verstraete and J. I. Cirac,  _Renormalization algorithms for Quantum-Many Body Systems in two and higher dimensions_ ,  arXiv e-prints cond-mat/0407066 (2004),  cond-mat/0407066.
  * [6] V. Murg, F. Verstraete and J. I. Cirac,  _Variational study of hard-core bosons in a two-dimensional optical lattice using projected entangled pair states_ ,  Phys. Rev. A 75, 033605 (2007),  10.1103/PhysRevA.75.033605.
  * [7] G. Vidal,  _Entanglement renormalization_ ,  Phys. Rev. Lett. 99, 220405 (2007),  10.1103/PhysRevLett.99.220405.
  * [8] M. A. Martín-Delgado, J. Rodriguez-Laguna and G. Sierra,  _Density-matrix renormalization-group study of excitons in dendrimers_ ,  Phys. Rev. B 65, 155116 (2002),  10.1103/PhysRevB.65.155116.
  * [9] S. Depenbrock and F. Pollmann,  _Phase diagram of the isotropic spin- $\frac{3}{2}$ model on the $z=3$ bethe lattice_,  Phys. Rev. B 88, 035138 (2013),  10.1103/PhysRevB.88.035138.
  * [10] H. Otsuka,  _Density-matrix renormalization-group study of the spin- $\frac{1}{2}$ $\mathrm{XXZ}$ antiferromagnet on the bethe lattice_,  Phys. Rev. B 53, 14004 (1996),  10.1103/PhysRevB.53.14004.
  * [11] B. Friedman,  _A density matrix renormalization group approach to interacting quantum systems on cayley trees_ ,  Journal of Physics: Condensed Matter 9(42), 9021 (1997),  10.1088/0953-8984/9/42/016.
  * [12] M. Gerster, P. Silvi, M. Rizzi, R. Fazio, T. Calarco and S. Montangero,  _Unconstrained tree tensor network: An adaptive gauge picture for enhanced performance_ ,  Phys. Rev. B 90, 125154 (2014),  10.1103/PhysRevB.90.125154.
  * [13] P. Silvi, V. Giovannetti, S. Montangero, M. Rizzi, J. I. Cirac and R. Fazio,  _Homogeneous binary trees as ground states of quantum critical hamiltonians_ ,  Phys. Rev. A 81, 062335 (2010),  10.1103/PhysRevA.81.062335.
  * [14] V. Murg, F. Verstraete, O. Legeza and R. M. Noack,  _Simulating strongly correlated quantum systems with tree tensor networks_ ,  Phys. Rev. B 82, 205105 (2010),  10.1103/PhysRevB.82.205105.
  * [15] K. Gunst, F. Verstraete, S. Wouters, Ö. Legeza and N. D. Van,  _T3ns: Three-legged tree tensor network states_ ,  Journal of Chemical Theory and Computation 14(4), 2026 (2018),  https://doi.org/10.1021/acs.jctc.8b00098,  Doi: 10.1021/acs.jctc.8b00098.
  * [16] L. Tagliacozzo, G. Evenbly and G. Vidal,  _Simulation of two-dimensional quantum systems using a tree tensor network that exploits the entropic area law_ ,  Phys. Rev. B 80, 235127 (2009),  10.1103/PhysRevB.80.235127.
  * [17] Y.-Y. Shi, L.-M. Duan and G. Vidal,  _Classical simulation of quantum many-body systems with a tree tensor network_ ,  Phys. Rev. A 74, 022320 (2006),  10.1103/PhysRevA.74.022320.
  * [18] D. Bauernfeind, M. Zingl, R. Triebl, M. Aichhorn and H. G. Evertz,  _Fork tensor-product states: Efficient multiorbital real-time dmft solver_ ,  Phys. Rev. X 7, 031013 (2017),  10.1103/PhysRevX.7.031013.
  * [19] D. Bauernfeind,  _Fork Tensor Product States: Efficient Multi-Orbital Impurity Solver for Dynamical Mean Field Theory_ ,  dissertation, Graz University of Technology (2018).
  * [20] J. Eisert, M. Cramer and M. B. Plenio,  _Colloquium_ ,  Rev. Mod. Phys. 82, 277 (2010),  10.1103/RevModPhys.82.277.
  * [21] Eisert J., Friesdorf M. and Gogolin C.,  _Quantum many-body systems out of equilibrium_ ,  Nature Physics 11, 124 (2015),  https://doi.org/10.1038/nphys3215 10.1038/nphys3215.
  * [22] M. Gruber and V. Eisler,  _Magnetization and entanglement after a geometric quench in the XXZ chain_ ,  arXiv e-prints arXiv:1902.05834 (2019),  1902.05834.
  * [23] V. Eisler and D. Bauernfeind,  _Front dynamics and entanglement in the xxz chain with a gradient_ ,  Phys. Rev. B 96, 174301 (2017),  10.1103/PhysRevB.96.174301.
  * [24] M. Collura, A. De Luca and J. Viti,  _Analytic solution of the domain-wall nonequilibrium stationary state_ ,  Phys. Rev. B 97, 081111 (2018),  10.1103/PhysRevB.97.081111.
  * [25] A. J. Daley, C. Kollath, U. Schollwöck and G. Vidal,  _Time-dependent density-matrix renormalization-group using adaptive effective hilbert spaces_ ,  Journal of Statistical Mechanics: Theory and Experiment 2004(04), P04005 (2004).
  * [26] S. R. White and A. E. Feiguin,  _Real-time evolution using the density matrix renormalization group_ ,  Phys. Rev. Lett. 93, 076401 (2004),  10.1103/PhysRevLett.93.076401.
  * [27] G. Vidal,  _Efficient Classical Simulation of Slightly Entangled Quantum Computations_ ,  Phys. Rev. Lett. 91, 147902 (2003),  10.1103/PhysRevLett.91.147902.
  * [28] G. Vidal,  _Efficient simulation of one-dimensional quantum many-body systems_ ,  Phys. Rev. Lett. 93, 040502 (2004),  10.1103/PhysRevLett.93.040502.
  * [29] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken and F. Verstraete,  _Unifying time evolution and optimization with matrix product states_ ,  Phys. Rev. B 94, 165116 (2016),  10.1103/PhysRevB.94.165116.
  * [30] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde and F. Verstraete,  _Time-Dependent Variational Principle for Quantum Lattices_ ,  Phys. Rev. Lett. 107(7), 070601 (2011),  10.1103/PhysRevLett.107.070601.
  * [31] C. Lubich, I. V. Oseledets and B. Vandereycken,  _Time integration of tensor trains_ ,  SIAM Journal on Numerical Analysis 53(2), 917 (2015),  10.1137/140976546,  https://doi.org/10.1137/140976546.
  * [32] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck and C. Hubig,  _Time-evolution methods for matrix-product states_ ,  arXiv preprint arXiv:1901.05824 (2019).
  * [33] M. Rizzi, S. Montangero and G. Vidal,  _Simulation of time evolution with multiscale entanglement renormalization ansatz_ ,  Phys. Rev. A 77, 052328 (2008),  10.1103/PhysRevA.77.052328.
  * [34] H. N. Phien, I. P. McCulloch and G. Vidal,  _Fast convergence of imaginary time evolution tensor network algorithms by recycling the environment_ ,  Phys. Rev. B 91, 115137 (2015),  10.1103/PhysRevB.91.115137.
  * [35] C. Hubig and J. I. Cirac,  _Time-dependent study of disordered models with infinite projected entangled pair states_ ,  SciPost Phys. 6, 31 (2019),  10.21468/SciPostPhys.6.3.031.
  * [36] P. Czarnik, J. Dziarmaga and P. Corboz,  _Time evolution of an infinite projected entangled pair state: An efficient algorithm_ ,  Phys. Rev. B 99, 035115 (2019),  10.1103/PhysRevB.99.035115.
  * [37] W. Li, J. von Delft and T. Xiang,  _Efficient simulation of infinite tree tensor network states on the bethe lattice_ ,  Phys. Rev. B 86, 195137 (2012),  10.1103/PhysRevB.86.195137.
  * [38] D. Nagaj, E. Farhi, J. Goldstone, P. Shor and I. Sylvester,  _Quantum transverse-field ising model on an infinite tree from matrix product states_ ,  Phys. Rev. B 77, 214431 (2008),  10.1103/PhysRevB.77.214431.
  * [39] M. Suzuki,  _Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations_ ,  Physics Letters A 146(6), 319 (1990).
  * [40] H.-D. Meyer, U. Manthe and L. Cederbaum,  _The multi-configurational time-dependent hartree approach_ ,  Chemical Physics Letters 165(1), 73 (1990),  https://doi.org/10.1016/0009-2614(90)87014-I.
  * [41] U. Manthe, H. Meyer and L. S. Cederbaum,  _Wave‐packet dynamics within the multiconfiguration hartree framework: General aspects and application to nocl_ ,  The Journal of Chemical Physics 97(5), 3199 (1992),  10.1063/1.463007,  https://doi.org/10.1063/1.463007.
  * [42] U. Manthe,  _A multilayer multiconfigurational time-dependent hartree approach for quantum dynamics on general potential energy surfaces_ ,  The Journal of Chemical Physics 128(16), 164116 (2008),  10.1063/1.2902982,  https://doi.org/10.1063/1.2902982.
  * [43] U. Manthe,  _Wavepacket dynamics and the multi-configurational time-dependent hartree approach_ ,  Journal of Physics: Condensed Matter 29(25), 253001 (2017),  10.1088/1361-648x/aa6e96.
  * [44] M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore and F. Pollmann,  _Time-evolving a matrix product state with long-ranged interactions_ ,  Phys. Rev. B 91, 165112 (2015),  10.1103/PhysRevB.91.165112.
  * [45] M. M. Rams and M. Zwolak,  _Breaking the entanglement barrier: Tensor network simulation of quantum transport_ ,  arXiv e-prints arXiv:1904.12793 (2019),  1904.12793.
  * [46] F. A. Y. N. Schröder and A. W. Chin,  _Simulating open quantum dynamics with time-dependent variational matrix product states: Towards microscopic correlation of environment dynamics and reduced system evolution_ ,  Phys. Rev. B 93, 075105 (2016),  10.1103/PhysRevB.93.075105.
  * [47] C. Lubich, T. Rohwedder, R. Schneider and B. Vandereycken,  _Dynamical approximation by hierarchical tucker and tensor-train tensors_ ,  SIAM Journal on Matrix Analysis and Applications 34(2), 470 (2013),  10.1137/120885723,  https://doi.org/10.1137/120885723.
  * [48] P. Silvi, F. Tschirsich, M. Gerster, J. Jünemann, D. Jaschke, M. Rizzi and S. Montangero,  _The Tensor Networks Anthology: Simulation techniques for many-body quantum lattice systems_ ,  SciPost Phys. Lect. Notes p. 8 (2019),  10.21468/SciPostPhysLectNotes.8.
  * [49] J. Kanamori,  _Electron correlation and ferromagnetism of transition metals_ ,  Progress of Theoretical Physics 30(3), 275 (1963),  10.1143/PTP.30.275.
  * [50] D. Bauernfeind, M. Aichhorn and H. G. Evertz,  _Comparison of MPS based real time evolution algorithms for Anderson Impurity Models_ ,  arXiv e-prints arXiv:1906.09077 (2019),  1906.09077.
  * [51] L. Kohn, P. Silvi, M. Gerster, M. Keck, R. Fazio, G. E. Santoro and S. Montangero,  _Superfluid to Mott transition in a Bose-Hubbard ring: Persistent currents and defect formation_ ,  arXiv e-prints arXiv:1907.00009 (2019),  1907.00009.


Suppose $f$ is partially hyperbolic, dynamically coherent, volume preserving,
and all the center exponents are negative at almost every point. Then the
center stable foliation of $f$ is upper leafwise absolutely continuous. This
is a fairly direct consequence of Pesin theory. Indeed, if all the Lyapunov
exponents are negative then the Pesin local stable manifold of almost every
point is a neighborhood of the point inside its center stable leaf. Then the
absolute continuity of Pesin laminations [35] implies that the center stable
foliation is upper leafwise absolutely continuous.

We close with a couple of conjectures on the issue of absolute continuity. The
first one deals with dissipative systems.

######  Conjecture 7.6.

Let $k>1$ and ${\mathcal{C}}_{k}$ be the space of partially hyperbolic,
dynamically coherent $C^{k}$ diffeomorphisms with mostly contracting center
direction. Then, for an open and dense subset,

  * •

if there is a unique physical measure then the center stable foliation is
absolutely continuous;

  * •

if there is more than one physical measure then the center stable foliation is
not upper leafwise absolutely continuous.

Examples of the second situation will appear in a forthcoming paper [49].

Figure 4.

######  Conjecture 7.7.

Let $k>1$ and ${\mathcal{V}}_{k}$ be the space of partially hyperbolic,
dynamically coherent, volume preserving $C^{k}$ diffeomorphisms whose center
Lyapunov exponents are negative at almost every point. Then, for an open and
dense subset, the center stable foliation is absolutely continuous.

Figure 4 outlines a scenario for these issues in a relevant special case,
namely near the map $f_{0}=g_{0}\times\operatorname{id}$ as in Theorem A.
Accessibility is assumed throughout (but is not needed for the negative
results in $\lambda^{c}\neq 0$). Generically means for open and dense in
$C^{k}$ topology, $k\geq 1$. Upper leafwise absolute continuity of the center
unstable is known for $\lambda^{c}>0$, as we have seen, and we have also found
an open subset with (full) absolute continuity of the center unstable.

## References

  * [1] F. Abdenur and M. Viana.  Flavors of partial hyperbolicity.  Preprint www.preprint.impa.br 2009.
  * [2] J. F. Alves.  SRB measures for non-hyperbolic systems with multidimensional expansion.  Ann. Sci. École Norm. Sup., 33:1–32, 2000.
  * [3] J. F. Alves.  Statistical analysis of non-uniformly expanding dynamical systems.  Lecture Notes 24th Braz. Math. Colloq. IMPA, Rio de Janeiro, 2003.
  * [4] J. F. Alves and V. Araújo.  Hyperbolic times: frequency versus integrability.  Ergodic Theory Dynam. Systems, 24:329–346, 2004.
  * [5] J. F. Alves, C. Bonatti, and M. Viana.  SRB measures for partially hyperbolic systems whose central direction is mostly expanding.  Invent. Math., 140:351–398, 2000.
  * [6] J. F. Alves, S. Luzzatto, and V. Pinheiro.  Markov structures for non-uniformly expanding maps on compact manifolds in arbitrary dimension.  Electron. Res. Announc. Amer. Math. Soc., 9:26–31 (electronic), 2003.
  * [7] J. F. Alves, S. Luzzatto, and V. Pinheiro.  Markov structures and decay of correlations for non-uniformly expanding dynamical systems.  Ann. Inst. H. Poincaré Anal. Non Linéaire, 22:817–839, 2005\.
  * [8] J. F. Alves and M. Viana.  Statistical stability for robust classes of maps with non-uniform expansion.  Ergodic Theory Dynam. Systems, 22:1–32, 2002.
  * [9] M. Andersson.  Robust ergodic properties in partially hyperbolic dynamics.  Trans. Amer. Math. Soc., 362:1831–1867, 2010.
  * [10] D. V. Anosov.  Geodesic flows on closed Riemannian manifolds of negative curvature.  Proc. Steklov Math. Inst., 90:1–235, 1967.
  * [11] A. Avila.  Density of positive Lyapunov exponents for quasiperiodic ${\rm SL}(2,\mathbb{R})$-cocycles in arbitrary dimension.  J. Mod. Dyn., 3:631–636, 2009.
  * [12] A. Avila, J. Santamaria, and M. Viana.  Cocycles over partially hyperbolic maps.  Preprint www.preprint.impa.br 2008.
  * [13] A. Avila and M. Viana.  Extremal Lyapunov exponents: an invariance principle and applications.  Inventiones Math., 181:115–178, 2010.
  * [14] A. Avila, M. Viana, and A. Wilkinson.  Absolute continuity, Lyapunov exponents, and rigidity.  Preprint www.preprint.impa.br 2009.
  * [15] A. Baraviera and C. Bonatti.  Removing zero central Lyapunov exponents.  Ergod. Th. & Dynam. Sys., 23:1655–1670, 2003.
  * [16] C. Bonatti, L. J. Díaz, and R. Ures.  Minimality of strong stable and unstable foliations for partially hyperbolic diffeomorphisms.  J. Inst. Math. Jussieu, 1:513–541, 2002.
  * [17] C. Bonatti, L. J. Díaz, and M. Viana.  Dynamics beyond uniform hyperbolicity, volume 102 of Encyclopaedia of Mathematical Sciences.  Springer-Verlag, 2005.
  * [18] C. Bonatti, X. Gómez-Mont, and M. Viana.  Généricité d’exposants de Lyapunov non-nuls pour des produits déterministes de matrices.  Ann. Inst. H. Poincaré Anal. Non Linéaire, 20:579–624, 2003\.
  * [19] C. Bonatti and M. Viana.  SRB measures for partially hyperbolic systems whose central direction is mostly contracting.  Israel J. Math., 115:157–193, 2000.
  * [20] M. Brin and Ya. Pesin.  Partially hyperbolic dynamical systems.  Izv. Acad. Nauk. SSSR, 1:177–212, 1974.
  * [21] M. Brin and G. Stuck.  Introduction to dynamical systems.  Cambridge University Press, 2002.
  * [22] K. Burns and A. Wilkinson.  On the ergodicity of partially hyperbolic systems.  Annals of Math.
  * [23] J. Buzzi, O. Sester, and M. Tsujii.  Weakly expanding skew-products of quadratic maps.  Ergod. Th. & Dynam. Sys., 23:1401–1414, 2003.
  * [24] D. Dolgopyat and A. Wilkinson.  Stable accessibility is $C^{1}$ dense.  Astérisque, 287:33–60, 2003.
  * [25] S. Gouezel.  Decay of correlations for nonuniformly expanding systems.  Bull. Soc. Mat. France, 134:1–31, 2006.
  * [26] F. Rodriguez Hertz, M. A. Rodriguez Hertz, A. Tahzibi, and R. Ures.  Creation of blenders in the conservative setting.  Nonlinearity, 23:211–223, 2010.
  * [27] F. Rodriguez Hertz, M. A. Rodriguez Hertz, and R. Ures.  Some results on the integrability of the center bundle for partially hyperbolic diffeomorphisms.  In Partially hyperbolic dynamics, laminations, and Teichmüller flow, volume 51 of Fields Inst. Commun., pages 103–109. Amer. Math. Soc., 2007.
  * [28] M. Hirsch, C. Pugh, and M. Shub.  Invariant manifolds, volume 583 of Lect. Notes in Math. Springer Verlag, 1977.
  * [29] A. Katok.  Lyapunov exponents, entropy and periodic points of diffeomorphisms.  Publ. Math. IHES, 51:137–173, 1980.
  * [30] R. Mañé.  Contributions to the stability conjecture.  Topology, 17:383–396, 1978.
  * [31] R. Mañé.  A proof of the ${C}^{1}$ stability conjecture.  Publ. Math. I.H.E.S., 66:161–210, 1988.
  * [32] V. Niţică and A. Török.  An open dense set of stably ergodic diffeomorphisms in a neighborhood of a non-ergodic one.  Topology, 40:259–278, 2001.
  * [33] D. Ornstein and B. Weiss.  On the Bernoulli nature of systems with some hyperbolic structure.  Ergodic Theory Dynam. Systems, 18:441–456, 1998.
  * [34] Ya. Pesin and Ya. Sinai.  Gibbs measures for partially hyperbolic attractors.  Ergod. Th. & Dynam. Sys., 2:417–438, 1982.
  * [35] Ya. B. Pesin.  Families of invariant manifolds corresponding to non-zero characteristic exponents.  Math. USSR. Izv., 10:1261–1302, 1976.
  * [36] V. Pinheiro.  Sinai-Ruelle-Bowen measures for weakly expanding maps.  Nonlinearity, 19:1185–1200, 2006.
  * [37] C. Pugh and M. Shub.  Ergodic attractors.  Trans. Amer. Math. Soc., 312:1–54, 1989.
  * [38] C. Pugh and M. Shub.  Stably ergodic dynamical systems and partial hyperbolicity.  J. Complexity, 13:125–179, 1997.
  * [39] C. Pugh, M. Viana, and A. Wilkinson.  Absolute continuity of foliations.  In preparation.
  * [40] V. A. Rokhlin.  On the fundamental ideas of measure theory.  A. M. S. Transl., 10:1–52, 1952.  Transl. from Mat. Sbornik 25 (1949), 107–150.
  * [41] D. Ruelle and A. Wilkinson.  Absolutely singular dynamical foliations.  Comm. Math. Phys., 219:481–487, 2001.
  * [42] M. Shub.  Global stability of dynamical systems.  Springer Verlag, 1987.
  * [43] M. Shub and D. Sullivan.  Expanding endomorphisms of the circle revisited.  Ergodic Theory Dynam. Systems, 5:285–289, 1985.
  * [44] M. Shub and A. Wilkinson.  Pathological foliations and removable zero exponents.  Invent. Math., 139:495–508, 2000.
  * [45] S. Smale.  Differentiable dynamical systems.  Bull. Am. Math. Soc., 73:747–817, 1967.
  * [46] M. Tsujii.  Physical measures for partially hyperbolic surface endomorphisms.  Acta Math., 194:37–132, 2005.
  * [47] C. H. Vásquez.  Statistical stability for diffeomorphisms with dominated splitting.  Ergodic Theory Dynam. Systems, 27:253–283, 2007.
  * [48] M. Viana.  Almost all cocycles over any hyperbolic system have nonvanishing Lyapunov exponents.  Ann. of Math., 167:643–680, 2008.
  * [49] M. Viana and J. Yang.  Contributions to the theory of maps with mostly contracting center direction.  In preparation.