Hi @fretchen! Welcome to the forum. For the hamiltonian you and @nlwach mention in the PhD thesis above you can follow two strategies.

A first strategy is to decompose the collective angular momentum operators into qubits as it is done in Eq. 2.42 of the same thesis. In this case for a spin j you will need 2j qubits. This strategy is somewhat inefficient in that you would use a Hilbert space of size 2^{2j} where you only ever visit the symmetric subspace of dimension 2j+1.

A second strategy is to use a dual rail encoding, i.e., two harmonic oscillators with cutoff 2j+1 to encode a spin j system (a nice treatment of this can be found in Sec. 3.8 of the book by Barnett and Radmore). In this case the collective operators are mapped according to

J_x \to \tfrac12 \left(a b^\dagger + a^\dagger b \right)

J_z \to \tfrac12 \left(a^\dagger a - b^\dagger b \right)

As mentioned above the gate generated by J_z^2 can be decomposed into Kerr and Cross Kerr gates while the gates generated by J_x and J_z correspond to a symmetric beamsplitter and a product of two rotation gates respectively.

If you use two harmonic oscillator with cutoff 2j+1 you will have a Hilbert space of dimensions (2j+1)^2 which is significantly better than using bare qubits.

Finally, the Hamiltonian you are trying to simulate H = H_1+H_2+H_3 = \chi J_z^2-\Omega J_x+\delta J_z is made of non commuting pieces [J_z,J_x] \neq 0 thus one strategy you can use is to approximate

\exp(i H) \approx \left[\exp(i H_1/N) \exp(i H_2/N) \exp(i H_3/N) \right]^N

via Trotterization.

Hope this helps,

Nicolas