  We perform classical simulations of the 127-qubit kicked Ising model, which
was recently emulated using a quantum circuit with error mitigation [Nature
618, 500 (2023)]. Our approach is based on the projected entangled pair
operator (PEPO) in the Heisenberg picture. Its main feature is the ability to
automatically identify the underlying low-rank and low-entanglement structures
in the quantum circuit involving Clifford and near-Clifford gates.
  We assess our approach using the quantum circuit with 5+1 trotter steps which
was previously considered beyond classical verification. We develop a Clifford
expansion theory to compute exact expectation values and use them to evaluate
algorithms. The results indicate that PEPO significantly outperforms existing
methods, including the tensor network with belief propagation, the matrix
product operator, and the Clifford perturbation theory, in both efficiency and
accuracy. In particular, PEPO with bond dimension $\chi=2$ already gives
similar accuracy to the CPT with $K=10$ and MPO with bond dimension
$\chi=1024$. And PEPO with $\chi=184$ provides exact results in $3$ seconds
using a single CPU.
  Furthermore, we apply our method to the circuit with 20 Trotter steps. We
observe the monotonic and consistent convergence of the results with $\chi$,
allowing us to estimate the outcome with $\chi\to\infty$ through
extrapolations. We then compare the extrapolated results to those achieved in
quantum hardware and with existing tensor network methods. Additionally, we
discuss the potential usefulness of our approach in simulating quantum
circuits, especially in scenarios involving near-Clifford circuits and quantum
approximate optimization algorithms. Our approach is the first use of PEPO in
solving the time evolution problem, and our results suggest it could be a
powerful tool for exploring the dynamical properties of quantum many-body
systems.
