Consider a spin-1/2 Ising model with time-dependent transverse field:
$$ H = - \sum_{i<j} J_{i, j} \sigma^z_i \sigma^z_j - \Gamma(t) \sum_i \sigma^x_i$$
Given the initial state $|\psi(0)\rangle$ (usually $|\psi(0)\rangle = \left( \frac{|0\rangle + |1\rangle}{\sqrt{2}} \right)^{\otimes N}$), what are the state-of-the-art methods for numerically calculating the expectation value of the final state with respect to an observable $\sigma^z_i$: $\langle\psi(t')|\sigma^z_i|\psi(t')\rangle$?
My current method
I believe this method is called discrete-time path-integral Monte Carlo method.
Referring to chapter 5.3 of (Binder and Heermann, 2019), by dividing the Hamiltonian into two non-commuting parts and apply Trotterization, we can obtain a corresponding classical Hamiltonian ($z_i^k = \pm 1$) that is one-dimensional higher: $$ H_c = - \sum_{k=1}^P \left( \sum_{i<j} J_{i, j} z_i^k z_j^k + J^\perp(t) \sum_i z_i^k z_i^{k+1} \right)$$ $$ J^\perp(t) = -\frac{PT}{2} \ln \tanh \frac{\Gamma(t)}{PT}$$ and then use a classical Monte Carlo method (e.g. simulated annealing) to sample the classical states at each time step. The expectation value $\langle\psi(t')|\sigma^z_i|\psi(t')\rangle$ is estimated by $\frac{1}{P} \sum_{k=1}^P z_i^k$.
