Bipartitioning infinite MPS into 2 semi-infinite chains

[jamarks]

I'm trying to compute the Resta Polarization of a Fermionic SSH chain according to the method described in https://arxiv.org/pdf/1810.07043.pdf equations 13-15. In particular, I need to compute 'particle imbalance' by partitioning the system into 2 semi-infinite chains, and tracing out one of these to get components of the reduced density matrix. When I use the

Code: Select all

get_rho_segment([0], 1)

function, I end up getting a trivial 2x2 matrix, and I'm pretty sure I'm not accurately cutting the system into semi-infinite chains.

Does anyone know how to bipartition an infinite MPS into 2 semi-infinite chains so that I can compute the Resta polarization?

Thank you in advance for the help

[gunheepark0229]

There is a relevant post here!

viewtopic.php?f=7&t=95

[Johannes]

You want to cut the system into two halves, but get_rho_segment extracts the density matrix for the *few* sites you give as argument, so that's the wrong function.

The equations from the cited paper is really what's calculated by average_charge of the MPS (up to the prefactor constants, I believe), take a look at the (two-line) implementation of the function and probablity_per_charge.

Note that for infinite MPS this value is only defined modulo 1, i.e. 0.123 and 1.123 are equivalent: we need to subtract (somewhat) arbitrary integer numbers of charges to get a finite value.
However, adiabatically following the change of the value still makes sence, since then really you are interested in the *differences* of the charge values only - this is what's done in the paper.

The examples/chern_insulators/chiral_pi_flux.py make a similar "experiment" of flux pumping and following the value of the average charge, you can probably base your code on that example.

[jamarks]

Hi Johannes,

Thank you so much for your helpful and detailed response. I am indeed aware that the polarization is only defined modulo 1, but would still like to be able to compute this quantity.

The code you reference, which demonstrates flux pumping, (and the function 'average_charge',) is close to what I am looking for. However, the example in the code is still for finite systems, whereas I am looking to do this for iMPS. I tried to modify the example and played around with 'probability_per_charge' for infinite systems, but was not able to obtain any reasonable results. How does one properly use these functions on iMPS to get the polarization for infinite systems?

[Johannes]

I'm not sure what's the problem here?

The chiral_pi_flux.py example does work directly in the thermodynamic limit with infinite MPS.
Running it generates chiral_pi_flux_charge_pump.pdf, the equivalent of Fig. 4 in the paper you quote, and chiral_pi_flux_ent_spec_flow.pdf, the equivalent of Fig. 8.
The P from the paper is simply what's returned by average_charge.

How did you define the model?
average_charge does only work if you use charge conservation. Is that the issue?

[jamarks]

Hi Johannes,

Thanks again for the reply. The issue is that I'm using GroupedSite, and in addition to the Polarization, I am also interested in obtaining the eigenvalues of the transfer matrix. I've found that when I define the GroupedSite using charges = 'same', the TransferMatrix ends up having no eigenvalues after a few evolution steps, so I am forced to use charges='drop'. On the other hand, average_charge only gives a non-empty array when charges = 'same'.

[Johannes]

Hang on - you're running DMRG to obtain ground states, right? What do you mean with the "evolution" steps?
Why do you need the GroupedSite at all, if you're only running DMRG? Can you specify the hamiltonian?
And the transfer matrix between which states? You might need to explicitly set charge_sector=None, but then the charge conservation should have no effect on the eigenvalues of the transfer matrix.

[jamarks]

I'm using TEBD to get the ground state, and to do the real-time evolution. Here's the model definition:

Code: Select all

from tenpy.models.model import NearestNeighborModel, CouplingMPOModel
from tenpy.models.lattice import Site, Chain
from tenpy.networks.site import FermionSite, GroupedSite

class SSHIntModel(CouplingMPOModel):
	def init_sites(self, model_params):
		fs = FermionSite(conserve = 'N')
		gs = GroupedSite([fs, fs], labels = ['A', 'B'], charges = 'drop')
		gs.add_op('Sz', gs.NA - gs.NB, False)
		return gs

	def init_lattice(self, model_params):
		bc_MPS = model_params.get('bc_MPS', 'infinite')
		bc = 'periodic' if bc_MPS == 'infinite' else 'open'
		nsys = model_params.get('nsys') if bc_MPS == 'finite' else 2
		self.L = nsys
		fs = self.init_sites(model_params)
		lat = Chain(nsys, fs, bc=bc, bc_MPS=bc_MPS)
		return lat

	def init_terms(self, model_params):
		J = model_params.get('J', 0.)
		d = model_params.get('d', 1.)
		tau = model_params.get('tau', 1.)
		V = model_params.get('V', 0.)
		dlt = model_params.get('dlt', 0.)

		#### Staggered Onsite Potential ####
		if dlt != 0.:
			self.add_onsite(dlt, 0, 'Sz')

		#### intracell hopping terms ####
		self.add_onsite(-J, 0, 'CdA CB')
		self.add_onsite(J, 0, 'CA CdB')

		#### intercell hopping terms ####
		hop_plus = (d + tau)/2.
		self.add_coupling(-hop_plus, 0, 'CdB', 0, 'CA', [1])
		self.add_coupling(-hop_plus, 0, 'CdA', 0, 'CB', [-1])

		hop_minus = (d - tau)/2.
		self.add_coupling(-hop_minus, 0, 'CdA', 0, 'CB', [1])
		self.add_coupling(-hop_minus, 0, 'CdB', 0, 'CA', [-1])

		#### interaction terms ####
		if V != 0:
			self.add_coupling(V, 0, 'NA', 0, 'NB', [1])
			self.add_coupling(V, 0, 'NB', 0, 'NA', [1])

def get_sshint_model(model_params):
	M = SSHIntModel(model_params)
	return NearestNeighborModel.from_MPOModel(M)			

[Johannes]

Why are you trying to do the time evolution? The paper you cited does *not* use TEBD, it uses DMRG!

We define \(\theta = 2\pi t/T\), assume adiabaticity and work in the instantaneous eigenbasis of \(H(\theta)\). Hence, the charge transport becomes independent of the time scale \(T\).

"Adiabatically" implies that the time evolution is slow enough to always stay in the ground state of the system.
So you really just run DMRG (which is way more efficient than imaginary TEBD!) to find the ground state for increasing \(\theta = 0, ..., 2\pi\) in sufficiently small steps; probably something like 20-30 steps in total is enough to capture the physics. To ensure the "adiabaticity" for the DMRG, you should initialize each run of DMRG with the initial state and environment from the previous \(\theta\).

[jamarks]

Thanks again for the reply. I fully recognize the different between the 'adiabatic' evolution that is entailed by charge pumping, and the real time evolution that I am doing. The basic idea I am looking at is quenches in order-obstructed SPT phases. As such, I am interested in looking at the evolution of the polarization during the course of a quench - so I do actually need to do the real-time evolution. And I know that the particular values taken by the polarization during the quench are not meaningful, but their changes over time are.

Also, the reason I am using TEBD rather than DMRG is that I am doing both real and imaginary evolution, working with infinite systems. And from what I can tell, TeNPy's dmrg/tdvp do not allow for real-time evolution of iMPS. Am I somehow mistaken?

Thanks,
Jacob

[Johannes]

I see. Sorry, if I appeared rude; I just wasn't sure what you wanted to do.

I'd still recommend running DMRG to find the ground state, which gives you an MPS, and then just take that as initial state for the TEBD, as outlined below.

On a side note, you can derive your model directly from the NearestNeighborModel, making the get_sshint_model funciton obsolete:

Python: Select all

class SSHIntModel(CouplingMPOModel, NearestNeighborModel):
    ...

Regarding the charges, your model should work perfectly fine with charges='same' for the GroupedSite.
After these changes, I can run

Python: Select all

from tenpy.networks.mps import MPS
from tenpy.algorithms import dmrg, tebd
# SSHIntModel defined above, with charges='same'
M = SSHIntModel({'bc_MPS': 'infinite', 'tau': 1., 'nsys': 4})
pstate = ['empty_A full_B'] * M.lat.N_sites
psi = MPS.from_product_state(M.lat.mps_sites(), pstate, bc=M.lat.bc_MPS)
print("before dmrg:", psi.average_charge())
dmrg.run(psi, M, {})
print("after dmrg:", psi.average_charge())
eng = tebd.Engine(psi, M, {'dt': 0.001, 'N_steps':1, 'order': 2})
for i in range(5):
    eng.run()
    print("t =", eng.evolved_time, "average charge =", psi.average_charge())

Doing that, I noticed that there was indeed an issue with the gauge fixing of the charge values, basically increasing the charge numbers on each bond with each TEBD update.
(It only affected the infinite case with non-zero qtotal in the unit cell.)
Anyways, I've fixed in in f93f27b63d89a9669e95f2553afbde8d6ac140d7.
After this commit, the average charge stays constant -0.5, as expected.

This was probably also the reason why you got only zero eigenvalues for the transfer matrix, if you have initialized it with charge_sector=0, which is the default and assumes that the dominant eigenvalue is in charge sector 0.
If you would have put charge_sector=None (which is the default if you use psi.overlap(other_psi), it checks *all* possible charge sectors.

While charge_sector=0 probably works after the fix, I think you should check for large eingenvalues in the other charge sectors as well; I'd expect that if you do a flux pump such that psi.average_charge() changes by one, the charge sector of the dominant eigenvector of the mixed transfer matrix will also change.

I hope this settles the problems

[jamarks]

Hi Johannes,

Thank you so much! Using charge_sector = None and pulling the most recent commit from TeNPy, I am now able to extract both the Polarization and the TM eigs when using charges='same'. Also switched to using DMRG for the ground state as you suggested. Everything works perfectly

Thanks for bearing with me and helping solve this issue. I really appreciate it!

Best,
Jacob