Theoretical models of spins coupled to bosons provide a simple setting for
studying a broad range of important phenomena in many-body physics, from
virtually mediated interactions to decoherence and thermalization. In many
atomic, molecular, and optical systems, such models also underlie the most
successful attempts to engineer strong, long-ranged interactions for the
purpose of entanglement generation. Especially when the coupling between the
spins and bosons is strong---such that it cannot be treated
perturbatively---the properties of such models are extremely challenging to
calculate theoretically. Here, exact analytical expressions for nonequilibrium
spin-spin correlation functions are derived for a specific model of spins
coupled to bosons. The spatial structure of the coupling between spins and
bosons is completely arbitrary, and thus the solution can be applied to systems
in any number of dimensions. The explicit and nonperturbative inclusion of the
bosons enables the study of entanglement generation (in the form of spin
squeezing) even when the bosons are driven strongly and near-resonantly, and
thus provides a quantitative view of the breakdown of adiabatic elimination
that inevitably occurs as one pushes towards the fastest entanglement
generation possible. The solution also helps elucidate the effect of finite
temperature on spin squeezing. The model considered is relevant to a variety of
atomic, molecular, and optical systems, such as atoms in cavities or trapped
ions. As an explicit example, the results are used to quantify phonon effects
in trapped ion quantum simulators, which are expected to become increasingly
important as these experiments push towards larger numbers of ions.