  An accurate description of the nonequilibrium dynamics of systems with
coupled spin and bosonic degrees of freedom remains theoretically challenging,
especially for large system sizes and in higher than one dimension. Phase space
methods such as the Truncated Wigner Approximation (TWA) have the advantage of
being easily scalable and applicable to arbitrary dimensions. In this work we
adapt the TWA to generic spin-boson models by making use of recently developed
algorithms for discrete phase spaces [Schachenmayer, PRX 5, 011022 (2015)].
Furthermore we go beyond the standard TWA approximation by applying a scheme
based on the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy of equations
[Pucci, PRB 93, 174302 (2016)] to our coupled spin-boson model. This allows in
principle to study how systematically adding higher order corrections improves
the convergence of the method. To test various levels of approximation we study
an exactly solvable spin-boson model which is particularly relevant for
trapped-ion arrays. Using TWA and its BBGKY extension we accurately reproduce
the time evolution of a number of one- and two-point correlation functions in
several dimensions and for arbitrary number of bosonic modes.
