To be honest, I'm quite surprised by
arXiv:2102.10188 that it actually seems to work with MPS. They simulate eq. 2, with all-to-all couplings \(H=\sum_{i<j} S_i \cdot S_j \). Naively, I would have expected much more entanglement than only ~log(N) to build up with the time evolution. On the other hand, it kind of makes sense, since it's effectively infinite dimensional and you expect a good description with mean field, so you might hope to not require too much entanglement.
Let me just say as a warning, that a priori, that there's
no guarantee that MPS-based methods will work as well for the more general H of eq. 1.
Are you aiming for the more general eq. 1 with arbitrary J_ij, or "just" trying to reproduce the results for eq. 2 for now?
The method they use is a variant of the ususal TEBD, but to incorporate the long-range interactions, each time step consists of N usual TEBD steps intertwined with Swap gates.
You have two alternatives to try to reproduce their results with TeNPy:
- Directly reproduce the same TEBD-like method with time steps and swap gates as they did.
You'd need to generalize TeNPy's TEBD update for this, following their Fig. 1.
Swap gates are implemented, you can just call permute_sites on the MPS to exchange neighboring sites.
TeNPy's TEBD usually just works for nearest-neighbor models, the M.H_bond
are the corresponing two-site H_{i,i+1} terms (not the M.H_MPO
). If you have the truly uniform eq. 2, you might be able to get away with generating just a nearest-neighbor model, since the couplings are uniform anyways; but I'm not sure about the boundaries i=0,L-1 right now; there's probably special cases to account correctly for the onsite terms there. (I didn't think this through completely...)
- Try another time evolution method based on the MPO, like TDVP.
The *uniform* all-to-all coupling can actually be written *very* efficiently with a bond-dimension 5 Matrix product operator of the form
\( W = \begin{pmatrix}
1 & Sp & Sm & Sz & h_x Sx + h_z Sz \\
0 & 1 & 0 & 0 & J Sm \\
0 & 0 & 1 & 0 & J Sp \\
0 & 0 & 0 & 1 & J Sz \\
0 & 0 & 0 & 0 & 1
\end{pmatrix} \)
This is a special case of the exponentially decaying interactions without decay, so you can use add_exponentially_decaying_interactions with lambda_ = 1.
to generate that MPO, like this:
Code: Select all
class InfiniteRangeHeisenberg(tenpy.models.model.CouplingMPOModel):
def init_sites(self, model_params):
conserve = model_params.get('conserve', 'best')
if conserve == 'best':
# check how much we can conserve
if not model_params.any_nonzero(['hx', 'hy'], "check Sz conservation"):
conserve = 'Sz'
else:
conserve = None
self.logger.info("%s: set conserve to %s", self.name, conserve)
return tenpy.networks.site.SpinHalfSite(conserve=conserve)
def init_terms(self, model_params):
hx = model_params.get('hx', 0.)
hz = model_params.get('hz', 0.)
self.add_onsite(-hx, 0, 'Sx')
self.add_onsite(-hz, 0, 'Sz')
J = model_params.get('J', 1.)
for op_i, op_j in [('Sp', 'Sm'), ('Sm', 'Sp'), ('Sz', 'Sz')]:
self.add_exponentially_decaying_coupling(J, 1., op_i, op_j)
You can then try to evolve the state with methods which directly use the MPO, i.e. TDVP and/or ExpMPOEvolution.
In all cases, you should definitely make sure that you can get the correct results first at small system sizes L=14/16 by comparison with exact diagonalization! In particular, both TDVP and ExpMPOEvolution might fail to correctly incorporate the evolution of the long-range terms. The ExpMPOEvolution since it's not designed for such a use-case - it assumes somewhat limited range interactions. Similarly, the projection of TDVP into the tangent space will probably have strong effects especially during the first few sweeps while you try to grow the bond dimension.
So fter all, the SWAP method might work better - but of course, you'd have to try to really see this.