For the time evolution of Heisenberg model (J_{2}=0), the trotter gates can be applied in the order of bonds (1, 2), (2, 3), …, (N-1, N), (N-1, N), …, (2, 3), (1, 2) for numerical stability. Is there any good arrangement of bond operators for the J_{1}-J_{2} Heisenberg chain? I know the “apply” function in ITensors can automatically handle the swap gates involved in applying the nonlocal gates.

It’s a good question and not totally obvious. But I believe a good way to think about this is to imagine grouping your sites together until all the interactions become at most nearest-neighbor.

Let’s say you group site 1 with site 2, then group site 3 with site 4, etc. As a result, some of the J_1 interactions will become completely local (on-site) and other J_1 interactions will become nearest-neighbor. All of the J_2 interactions will become nearest-neighbor for these new “grouped sites”.

Then you would apply the usual, nearest-neighbor Trotter decomposition for these grouped sites. That suggests the following pattern (where I hope the notation is clear to you):