How can I conceptualize virtual indices used in the time evolving block decimation (TEBD) algorithm?

Question

I am trying to work with the Heisenberg $XX$ model which the Hamiltonien is given by $\hat{H} = -J \sum_i \left(\hat{X}_i\hat{X}_{i+1} + \hat{Y}_i\hat{Y}_{i+1}\right), \quad J > 0.$ Using $s^+ = |{0}\rangle \langle{1}|$ and $s^- = |1\rangle \langle 0|$, I have redefine the Hamiltonien such as $\hat{H} = -2J\sum_i (s_i^+s_{i+1}^- + s_i^-s_{i+1}^+)$. I want to evolve a Matrix Product State (MPS) with the mentionned Hamiltonien using the TEBD algorithm. The block representation of the TEBD uses virtual indices, called bonds. I was wondering how could I associate these indices with the ones from the tensor on each site of my MPS ? More specifically, I seem to have trouble contracting the right (meaning the good one) indices with the axes = () option using np.tensordot().

For example, how to choose the right indices to contract between tensors on adjadcent site with the np.tensordot() only knowing the shape that should have the resulting tensor? It seems to me like there is more than one way to contract indices for two adjacent tensors such that the resulting one have always a specific shape. I hope my question is not too nebulous. Thank you.

EDIT: The figure is taken from: Block representation of the TEBD algorithm (Efficient numerical simulations with Tensor Networks: Tensor Network Python (TeNPy), Johannes Hauschild1 and Frank Pollmann. (2018).

Answer 1

I think this seems to be more a python problem. See an explanation of the numpy tensordot function here: https://stackoverflow.com/questions/41870228/understanding-tensordot

Regarding your question, I think it depends on your convention of ordering virtual and physical bonds. As an example, if you have two adjacent MPS tensors, say $B_n$ with shape $(\alpha_n, j_n, \alpha_{n+1})$ and $B_{n+1}$ with shape $(\alpha_{n+1}, j_{n+1}, \alpha_{n+2})$, where $\alpha_i$ are virtual bonds and $j_i$ physical bonds; sweeping from left to right you now want to sum over the last axes of $B_{n}$ and the first axes of $B_{n+1}$. The first thing you'd need to do according to your figure would be:

theta_n = np.tensordot(B_n, B_n_1, (-1, 0))

Note that $\Theta_n$ will have the shape $(\alpha_n, j_n, j_{n+1},\alpha_{n+2})$. For the next step (i), you need to contract the 2nd and 3rd axes of $\Theta_n$ with the respective axes of the two-side unitary operator.

I hope this helps!