<p>The accurate simulation of wave propagation in real media requires properly taking the attenuation into account, which leads to wave dissipation together with its causal companion, wave dispersion. In this study, to obtain a weak formulation of heterogenous viscoacoustic wave propagation in an infinite domain, the viscoacoustic medium is first characterized by its frequency-dependent complex bulk compliance instead of the classically used complex bulk modulus. Then, a mechanical model using serially connected standard linear solids (SSLS) is built to obtain the rational approximation of the complex bulk compliance whose parameters are calculated via an adapted nonlinear optimization method. Utilizing the obtained bulk compliance-based constitutive relation, a novel second-order viscoacoustic wave equation in the frequency domain is derived, of which the weak formulation can be physically explained as the virtual work equation and can thus be discretized using a continuous spectral element method in space. Additionally, a new method is introduced to address the convolution terms involved in the inverse Fourier transform, whose accurate time integration can then be achieved using an explicit time scheme, which avoids the transient growth that exists in the classical method. The resulting full time-space decoupling scheme can handle wave propagation in arbitrary heterogeneous media. Moreover, to treat the wave propagation in an infinite domain, a perfectly matched layer in weak formulation is derived for the truncation of the infinite domain via complex coordinate stretching of the virtual work equation. With only minor modification, the resulting perfectly matched layer can be implemented using the same time scheme as for the wave equation inside the truncated domain. The accuracy, numerical stability, and versatility of the new proposed scheme are demonstrated with numerical examples.</p>