Hardcore two-species bosonic system implementable within tJ sitetype?

I wish to implement a hardcore boson model with two species. In particular, I am interested in dynamics of the Bose-Hubbard Hamiltonian in the large U/t limit with additional nearest neighbour Ising terms:

\hat{H} = -\sum_{\braket{{\bf i}, {\bf j}},\sigma}t_{\bf ij}(\hat{a}^\dagger_{{\bf i}\sigma}\hat{a}_{{\bf j}\sigma} + {\rm h.c.}) - \sum_{\braket{{\bf i}, {\bf j}}}J_{\bf ij}\left( \Delta_{\bf ij} \hat{S}^z_{\bf i} \hat{S}^z_{\bf j} + \frac{1}{2}\left(\hat{S}^+_{\bf i} \hat{S}^-_{\bf j} + \hat{S}^+_{\bf j}\hat{S}^-_{\bf i}\right) + \frac{3\hat{n}_{\bf i}\hat{n}_{\bf j}}{4}\right ),

where \Delta_{\bf ij} = (J_{\bf ij} + V)/J_{\bf ij}, with J_{\bf ij} = 4t_{\bf ij}^2 / U and V the additional nearest-neighbour Ising term. Here, the S-operators are the usual pseudo-spin operators. E.g. S^z_{\bf j} = (\hat{a}^\dagger_{{\bf j}\uparrow}\hat{a}_{{\bf j}\uparrow} - \hat{a}^\dagger_{{\bf j}\downarrow}\hat{a}_{{\bf j}\downarrow})/2.

Now, I wish to look at dynamics of this model in a two-leg ladder of size N = N_x \times 2 using TEBD with ITensor in Julia. I am aware that one can in principle implement a two-component Bose-Hubbard model by considered a system of 2N sites where every even (odd) site refers to bosons in the \uparrow (\downarrow) state. However, since the system is on a two-leg ladder, I already need to perform SWAP gates on every rung of the ladder as e.g. done here for a fermionic t-J model: https://doi.org/10.1103/PhysRevLett.115.056401. If I were also to double the system size to accommodate the two species of the boson, I would, therefore, need to perform even more SWAP gates to implement the dynamics (as far as I can tell). I would like to avoid this.

However, I was wondering if one could not just use the tJ sitetype and then instead of using fermionic “Cdagup”, “Cdagdn” etc. which incorporates the appropriate Jordan-Wigner string for fermions, I simply use the “Adagup”, “Adagdn” operators instead? These should, naively be hardcore bosons. The SWAP gates (used to swap the legs of the ladder) should then be implemented by the following operator as far as I can tell:

ITensors.op(::OpName"boson_SWAP", ::SiteType"tJ") = [
+1 0 0 0 0 0 0 0 0
0 0 0 +1 0 0 0 0 0
0 0 0 0 0 0 +1 0 0
0 +1 0 0 0 0 0 0 0
0 0 0 0 +1 0 0 0 0
0 0 0 0 0 0 0 +1 0
0 0 +1 0 0 0 0 0 0
0 0 0 0 0 +1 0 0 0
0 0 0 0 0 0 0 0 +1
];

Here, the basis is ordered as \{\ket{00},\ket{0\uparrow},\ket{0\downarrow},\ket{\uparrow 0},\ket{\uparrow\uparrow},\ket{\uparrow\downarrow},\ket{\downarrow 0},\ket{\downarrow\uparrow},\ket{\downarrow 0} \}.

So my basic question is this: I am missing something here, or should this “just work”? Moreover, can I use the built in spin operators straight away, or should I write them explicitly in terms of the “A” operators [like S^z_{\bf j} = (\hat{a}^\dagger_{{\bf j}\uparrow}\hat{a}_{{\bf j}\uparrow} - \hat{a}^\dagger_{{\bf j}\downarrow}\hat{a}_{{\bf j}\downarrow})/2]?

I made some quick tests. When the hopping, t_{\bf ij}, across the ladder is 0, the system behaves just like the fermionic case, as far as I can tell. This should also be the case, I believe. However, when the hopping across the ladder is nonzero, the hardcore boson case seems to deviate very strongly from the fermionic case.