By using a small number of Gaussian basis functions, one can synthesize the wave fields radiated from planar and focused piston transducers in the form of a superposition of Gaussian beams. Since Gaussian beams can be transmitted through complex geometries and media, such multi-Gaussian beam models have become powerful simulation tools. In previous studies the basis function expansion coefficients of multi-Gaussian beam models have been obtained by both spatial domain and $k$ -space domain methods. Here, we will give an overview of these two methods and relate their expansion coefficients. We will demonstrate that the expansion coefficients that have been optimized for circular piston transducers can also be used to generate improved field simulations for rectangular probes. It will also be shown that because Gaussian beams are only approximate (paraxial) solutions to the wave equation, a multi-Gaussian beam model is ultimately limited in the accuracy it can obtain in the very near field.