Matrix Product Toolkit | Main / NanoRing

Hamiltonian

The basic model is the Hubbard Hamiltonian with periodic boundary conditions and a magnetic flux,

{%ham,H = -t \sum_{\sigma,i=1}^{L} \left( e^{i 2\pi \phi / L} c_{i,\sigma}^\dagger c_{i+1,\sigma}^{} + \mathrm{H.c.} + U (n_{i,\uparrow} - 1/2)(n_{i,\downarrow}-1/2) \right) %}

Here {$\phi$} is the total flux in units of the flux quanta. Note that we have used the particle-hole symmetric form of the Coulomb interaction. We have not yet introduced any impurities, they will be added separately.

Splitting into real and imaginary parts, and assuming {$t=1$},

{%H = \cos(2\pi \phi / L) H_t + \sin(2\pi \phi L) H_j + U H_U %}

where we have set {%H_t = -\sum_{\sigma,i=1}^L \left( c_{i,\sigma}^\dagger c_{i+1,\sigma} + c_{i+1,\sigma}^\dagger c_{i,\sigma} \right) %} {%H_j = -i \sum_{\sigma,i=1}^L \left( c_{i,\sigma}^\dagger c_{i+1,\sigma} - c_{i+1,\sigma}^\dagger c_{i,\sigma} \right) %} {%H_U = \sum_{\sigma,i=1}^L \left(n_{i,\uparrow} - 1/2\right)\left(n_{i,\downarrow}-1/2) \right) %}

Note that {$H_j$} is also the current operator. These operators are constructed by the hubbard-periodic-u1su2 program, so the basic Hamiltonian as used by the DMRG code will be of the form

lattice:"cos($phi*2*pi/L)*H_t + sin($phi*2*pi/L)*H_j + $U*H_U"

where $phi is the magnetic flux and $U is the interaction strength. This is enough to do a simple calculation.

Walkthrough simple calculation for a 10-site chain

We firstly construct the lattice file for 10 sites, saved as the file ring-10. This lattice is defined using {$SU(2)$} spin symmetry. This improves the efficiency of the calculation, but we need to be careful when calculating observables as we only have the {$SU(2)$} reduced matrix elements available, not the usual up,down spin operators. Some examples later.

hubbard-periodic-u1su2 10 ring-10

Now construct a random initial state that has the correct quantum numbers. Here we choose half-filling, with zero total spin, so {$N=10, S=0$}, and we save the wavefunction as the file psi.

mp-random -l ring-10 -q 10,0 -o psi

This is probably a good point to introduce the mp-info command. In the simplest mode, it will give some basic information about the state. Running mp-info psi gives as output:

Symmetry list is N:U(1),S:SU(2)
State transforms as 10,0
Number of sites = 10

Mp-info has lots of other uses, for example to see how many states are kept at each bond of the wavefunction, try mp-info -s psi. Since the wavefunction we are using is random, this might be slightly different but a typical output is

#bond    #dimension  #degree
    0             1        1
    1             3        4
    2             8       11
    3             9       16
    4             9       15
    5            10       19
    6             9       16
    7             8       13
    8             6        9
    9             3        4
   10             1        1

Here we can already see the effect of the {$SU(2)$} symmetry. The maximum number of states kept was 10 {$SU(2)$} reps, equivalent to 19 states using {$U(1)$} symmetry.

Now lets obtain the groundstate. Lets choose a flux of 0.25 and U=1, m=300 states kept. Using the simple DMRG program mp-dmrg,

mp-dmrg -w psi -H ring-10:"cos(0.25*2*pi/10)*H_t + sin(0.25*2*pi/10)*H_j + 1.0 * H_U" -m 300

This will do 2 half-sweeps through the system, which won't be enough to get to the groundstate. The useful measure is the variance, {$\langle (H-E)^2 \rangle$}, which will be close to zero when the calculation has converged. At this point, the variance measure will be not much smaller than 1, still way too high. Repeating the same command will run a couple more sweeps. As we approach the groundstate, better convergence is usually achieved by switching to the 2-site DMRG algorithm with no mixing term (see mp-dmrg help page for more info). To do this, we add the options -2 -f 0 to mp-dmrg,

mp-dmrg -w psi -H ring-10:"cos(0.25*2*pi/10)*H_t + sin(0.25*2*pi/10)*H_j + 1.0 * H_U" -m 300 -2 -f 0

This should bring the variance down to close to 1E-8, perhaps two additional sweeps will get there. (you can also add more sweeps by using the -s option. See mp-dmrg for more details).

Now we can calculate the persistent current. This is defined as the derivative of the ground state energy with respect to the flux,

{%I = - \frac{dE}{d\phi} %}

It would be possible to numerically evaluate this derivative by calculating the energy for two nearby values of the flux. A better approach is to utilize the Hellmann-Feynman theorem {$ dE / dx = \langle dH / dx \rangle $}, and calculate the expectation value

{%I = \frac{2\pi}{L} \left(-\sin(2\pi \phi / L) H_t + \cos(2\pi \phi L) H_j \right)%}

which translates into

mp-expectation -r psi ring-10:"(2*pi/10) * (-sin(0.25*2*pi/10)*H_t + cos(0.25*2*pi/10)*H_j)"

This should give a number very close to 1.197129, which appears to be consistent with the curve given in figure 1 of Czakja et al, in Phys. Rev. B 72, 035320 (2005).

A figure calculated in this fashion with more data points:

There is some smoothing of the curve, especially for U=0 that should be a perfect triangle waveform. Perhaps this can be improved by keeping more states (around 550 will turn it into an exact diagonalization). But this is probably unavoidable on longer chains.

Adding impurities

Any impurities can be added by adding impurity terms to the Hamiltonian. For example, a site-dependent chemical potential (proportional to the particle number N is simple to add. The sites are numbered from 1 to L, so we just add terms like w * N(1). To modify the Coulomb term at a particular site, the symmetric version of the coulomb operator is denoted by P/4 (this arises because P is the single-site parity operator {$(-1)^N$}). So, to modify U at site 3 we could add a term like U*P(3)/4 to the Hamiltonian. Modifying the hopping term is simple, but we need to remember that in the {$SU(2)$} formulation we don't have access to the {$c^\dagger_{\uparrow,\downarrow}$} operators, but only the reduced matrix elements. In effect, this means we have the {$c,c^\dagger$} operators represented as two-component spinors, and we can only ever manipulate as a whole (eg, by taking the dot product). This means that the expression, for example for hopping from site 2 to site 3 takes the form

dot(CH(2), C(3))

and the Hermitian conjugate is dot(C(2), CH(3)). This might look strange - since the C operators should anticommute, one might expect the Hermitian conjugate of dot(CH(2), C(3)) to be dot(C(3), CH(2)) = -dot(C(2), CH(3)). In the {$SU(2)$} formulation, it is still true that {$c,c^\dagger$} anticommute (so dot(C(3), CH(2)) = -dot(C(2), CH(3)) always). But for a spinor operator, taking the hermitian conjugation twice does not give you the original operator, but gives an additional minus sign. So {$(c^\dagger)^\dagger = -c$} and this cancels out with the fermion sign.

In summary, a hopping term at phase angle 0.25 between sites 2 and 3 appears as the expression

-(exp(0.25*2i*pi/10) * dot(CH(2), C(3)) + exp(-0.25*2i*pi/10) * dot(C(2), CH(3)))

If we calculate the expectation value of this operator on the state psi calculated above, we should get -1.2602653, which is exactly 1/10 of the kinetic energy (calculated as cos(0.25*2*pi/10)*H_t + sin(0.25*2*pi/10)*H_j).

where we make use of complex arithmetic to save expanding out the sin and cos terms, the expression evaluator understands complex numbers like 2i. Remember also that the coefficient of the hopping term is {$-t$}, hence the minus sign out the front. Modifying the magnitude of the hopping at a particular bond should be easy. Modifying the phase is a bit more attention to get the arithmetic right. Note that the expression evaluator will optimally factorize the expression, so there is no harm in repeating operators multiple times. There is no reason to simplify by hand an expression just for the sake of making it efficient: the DMRG code does this anyway.

Other commands that will be useful

For 'serious' DMRG calculations (say, in a batch queue), use Tools.MpDmrgInitResume. To see some other info on a wavefunction (including the Hamiltonian that was used), use Tools.MpAttr. It would be very interesting if there is a way of using the fidelity, ie the overlap between groundstates calculated at different interactions {$\langle \phi_1 \vert \phi_2 \rangle $}. This can be done using Tools.MpOverlap. Of course, Tools.MpExpectation can also calculate expectation values between unrelated wavefunctions as well. I don't know if this is useful for this project, but I don't think anyone ever considered this before.