# Computational Quantum Physics

## Matrix-Product-States-based Methods

### DMRG

**Task:** Given Hamiltonian $H$, find its ground state, i.e.

$E_0 = min {\langle\psi|H|\psi\rangle \over \langle \psi|\psi \rangle}$

**Basis:** Renormalziation group. In other words, the microscopical characteristics of a system are of minor importance near a critical point or in a bulk phase.

#### General Procedure

* Instead of miminzing the equation $E_0 = min {\langle\psi|H|\psi\rangle \over \langle \psi|\psi \rangle}$ with repect to all the $\Gamma$ and $\Lambda$ matrices, we can minimize the expectation value $\langle\phi|H|\phi\rangle$ under the condition that $\langle|\phi|\phi\rangle = 1$. For this purpose, we introduced a Lagrangian multiplier to extremize the functional

$f[|\psi\rangle]=\langle\psi|H|\psi\rangle - \lambda (1-\langle\psi|\psi\rangle)$

It is computationally most stable to split up the optimization process and optimize one site $l$ at a time, while keeping the others fixed.

Here we contracted all tensors around the site $l$ into $\mathcal{H}_l$ and $\mathcal{N}_l$. These operators act on the combiend space of left and right virtual and physical indices, so that their matrix dimensions are $\chi_l\chi_{l+1}d\times\chi_l\chi_{l+1}d$. If we write all the parameters of site $l$ as a vector $\vec{a}^{[l]}$, then the functional becomes

$f_l[\vec{a}^{[l]}]=(\vec{a}^{[l]})^\dagger \mathcal{H}_l \vec{a}^{[l]} + \lambda [1-(\vec{a})^{[l]}\mathcal{N}_l\vec{a}^{[l]}]$.

* To proceed, we groupd the $\Gamma$ and $\Lambda$ martrices in such a way that all the matrirces on the left of $l$ are written as $A^{[i<l]}=\Lambda^{[i]}\Gamma^{[i]}$ and all the matrices on the right of $l$ are $B^{[i>l]}=\Gamma$. This makes $\mathcal{N}_i=I$, minimizing the quadratic functional now becomes
$\mathcal{H}_i \vec{a}^{[i]}=\lambda \vec{a}^{[i]}$
which can be solved using standard numerical libraries.

* To optimzie the entire MPS, we start from the left edge of the system working our way to the right edge and then back. This righ-left moving defines one *sweep*. By sweeping until the energy $\lambda$ converges to a predefined accuracy, we obtain an excellent approximations of the ground state in case of a one-dimensional gapped system.

### DMRG Algorithm

#### Single Site

* Start from an intial guess for $|\phi\rangle$ which is right-normalized, i.e., consists of $B$ matrices only. This is often simply a random state. Note that ther should be a non-zero overlap with the ground state.
* Calculate the right environments, or $R$-expression, for all site positions $L$ through 2 iteratively as discussed before.
* Right sweep through sites $l=1,...,L-1$:
* Solve the standard eigenvalue problem by an iterative eigensolver for $M^{[l]}$. Explicitly, we need to solve
$\sum_{\sigma_1'a_l'a_{l+1}'b_lb_{l+1}}(L_{a_l a_{l}'}^{b_l}H_{b_lb_{l+1}}^{\sigma_l'\sigma_{l}}R_{a_{l+1}a_{l+1}'}^{b_{l+1}})M_{a_l'a_{l+1}'}^{\sigma_l'}-\lambda M_{a_la_{l+1}}^{\sigma_l}=0$
We only need to compute the **lowest eigenvalue**, it provides an estimate for the ground state energy.
* We perform SVD on the new $M^{[l]}$ to left normalize the site tensor into $A^{[l]}$. We then multiply the remaining matrices from SVD to $B^{[l+1]}$, which yields the new $M^{[l+1]}$ as a *strating guess* for the next site.
* Build $L^{[l]}$ iteratively.
* Move by one site $l\to l+1$ and repeat.
* Left sweep through sites $l=L,...,2$
* Solve eigenvalue probelm for $M^{[l]}$
* Right normalzie $M^{[l]}$
* Build $R^{[l]}$
* Repeart left and right sweeps until the energy is converged, or beteer, until $var(H)<\epsilon$

##### Advantages

Single-site tends to achieve slightly better result than the two-site version fir a fixed $\chi$.

##### Drawbacks
It tends to get stuck in local minima. A more stable implementation, two -site DMRG, groups two sites together to update them simulataneously.

#### Two-Site DMRG

##### Advantage also Drawbacks

It allows for dynamically adjusting the bond dimension $\chi$ on the fly as one optimizes the wavefunction. This flexibility comes about because in splitting the optimized two-site tensor into single-site tensors using SVD, we can tune the number of states $\chi_{max}$ kept. Note, however, that keeping a fixed $\chi_{max}$ in this approach requires finding the MPS with bond dimension $\chi_{max}$ that best approximates the MPS after SVD, which in general has bond dimension $d\chi_{max}$.

### Time-Evoling block decimation

We can use imaginary time evolution to find the ground state:

$|\phi_0\rangle = \lim_{\tau\to \infty} {e^{-\tau H}|\psi_i\rangle \over ||e^{-H\tau}|\psi_i\rangle||}$

#### Basis
* The split-operators
* Trotter-Suzuki decomposition

$H = \sum_{k=1}^K h_k$

For the XXZ model, we have seen that we can split the Hamiltonian into an even and odd part. Each of these two then consistes of commuting two-site terms $H_i$, such that the time-evolution operator for the odd parts, as an example, takes the form

$U_{odd}(\Delta t)=\Pi_{i\in \mathcal{O}}U_i(\Delta t)$。

We are applying two-site gates to the MPS.

We could not in principle apply all the individual gates to the MPS and decouple the individual gates to the MPSs and decouple the resulting tensor into its MPS form including *truncation* of the bond dimension.

The individual gate operations can be performed very efficiently within the MPS formalism and importantly, the unitary evolution operations preserve the canonical form.

#### Algorithm

* Calculate the two-site tensor $\theta^{[l]}$ around bond $l$ by contracting the diagram.
In particular, we conrtact here the tensors $\Lambda^{[l]}$,$\Gamma^{[l]}$,$\Lambda^{[l+1]}$,$\Gamma^{[l+1]}$,$\Lambda^{[l+2]}$.
* Apply the gate $U^{[l]}$ at bound $l$ to find the new two-site tensor $\tilde{\theta}^{[l]}$
* Reshape the two-site tensor $\tilde{\theta}^{[l]}$ and perform an SVD to split it into single-site tensors.
In the last step, we have reshaped the matrices into tensors and separated out the matrices $\Lambda$ matrices by multiplying $A(B)$ from the left (right) with $\Gamma^{-1}$ to recover the canoncial form.
* Finally, to avoid a growth of the bond dimension to $d·\chi_{max}$, we discard the smallest Schmidt values in $\Lambda^{[l+1]}$ keeping only the largest $\chi_{max}$ ones. In order for our state to stay normalzied, we renomralzie $\Lambda^{[l+1]}$ such that $\sum_{i}\Lambda_{ii}=1$ after truncation. We have thus recovered an MPS states in canonical form with desired bond dimension.

#### Computation Errors in TEBD

* Truncation Error: this error is usually the main error when performing real-time evolution as the required bond dimension for exact time evolution grows expoentially.
* Trotter Error: relatively harmless and can be reduced by reducing $\Delta t$ or performing a higher-order expansion.
* To recover the canoncial form, we multiply both A and B in step 3 with $\Lambda^{-1}$. For small singular values, this operation is not stable. However, this sort of instability can be avoided with little extra computational cost.
* If we perform imaginary-time evolution, we need to keep in mind that the canonical form is only retained when $\Delta \tau \to 0$.

#### Implementation of TEBD

* Write a function that generates a matrix product state in the canonical form.

* Implement the two-site Hamiltoian operator for the Heisenberg model
* We also need a cuntion which computes the unitary evolution operator
* Now we want to evolve the wave function by applying the operator $U$ to our initial MPS. To this end, we first calculate the action on even bonds, followed by all odd bonds. For each bond update, we perform the following procedure:
* Construct the two-site tensor at the sites $i,i+1$ by contracting the corresponding $\Gamma$ and $\Lambda$ matrices over their shared virtual bonds.
* The resulting $\theta$ should be rank 4 with two virtual and two physical indices.
* Apply the bond operator $U$
* To ensure that the wavefunction remains compressed with our chosen bond dimension $\chi_{max}$, we have to truncate it by performing a singular value decomposition on the tensor $\theta$ and keeping only the largest $\chi_max$ singular values.
* The SVD results in three new tensor, the new $\tilde{\Gamma}$ matrix of Schmdit values at bonds $i+1$ and two tensors $A,B$ left and right of it.
* Bring the MPS back to canonical form, we multiply both tensors $A,B$ with the inverse of $\Gamma$ at the corresponding sites.
* The final update of the MPS changes the $\Gamma$ tensors at site $i,i+1$ and updates the $\Lambda$ values as obtaiend from the truncation.

#### Mixed states and open-quantum system

Markovian and Lindblad equation

$\dot{\rho}=-i[\rho,H]+\gamma \sum_i [L_i^\rho L_i^\dagger -{1\over2}\{L^\dagger L, \rho\}] $

To treat this evolution, we rearrange the legs of the $d^N\times d^N$ operator matrix. We arrive at a superstate.

The Lindbladian $\mathcal{L}$ is then a linear superoperator acting on $\rho$, such that we can formally integrate the time evolution $|\rho(t)\rangle_{\#} = exp(-i\mathcal{L}_{\#}t)|\rho_0\rangle_{\#}$.

### TDVP

TDVP: Time-dependent variational principle (TDVP) in quantum mechanics. Introducing the action functional

$S=\int_{t_1}^{t_2}dt L[\psi(t),\bar{\psi}(t)]$

with

$L[\psi(t),\bar{\psi}(t)] = \langle \psi(t)|i\partial_t-H |\psi(t)\rangle$

yiedling upon extremizing the time-dependent Schrodinger equation. This variational principle can be used to simulate the time evolution of a state also in a subspaceof the full Hilbert space, for our purposes, the space of fixed-bond-dimension amtrix product states.

**TDVP is useful in simulating long-range Hamiltonians given in MPO form**, for which **TEBD is not avialable**.

Another useful property of the TDVP is that by construction, it conserves both the norm and the energy of the state under real-time evolution.

### Symmerty

For a state with definite total magneticzation, we can force Schmidt values accordingly. We can then perform a full state decomposition, always grouping Schmidt values that corresponds to a fixed magnetization to the left and to the right of the bond.

### Two-Dimensional System

MPS allows us to calculate the ground state, even in **Fermion system**, where sign-problem-free calculation can be performed.

To implement it in higher dimension, a comon tactic is to consider strip geometries or cylinders and unwrap the transverse direction.

The entanglement expected from the area law noe scales as $L_y$, and w expect the bond dimension to $L_y~exp(L_y)$. Nevertheless, system of width up to $12$ can be simulated that way. The alternative is to study the natural extension of the MPS, named projected entangled pair states (PEPS). But they are hard to compute.