It is always possible to find a purification of a Gaussian state in terms of a Gaussian state.

This is most easily seen in terms of covariance matrices. A convenient way of doing so is to bring your covariance matrix into normal form by sympectic transformations. This decouples the system into a set of decoupled normal modes at some temperature, which can be easily purified by entangling each of them with another bosonic mode (i.e., a two-mode squeezed states). Undoing the sympectic transformation yields a covariance matrix for a Gaussian purification of the original state.

It is always possible to find a purification of a Gaussian state in terms of a Gaussian state.

This is most easily seen in terms of covariance matrices. A convenient way of doing so is to bring your covariance matrix into normal form by sympectic transformations. This decouples the system into a set of decoupled normal modes at some temperature, which can be easily purified by entangling each of them with another bosonic mode (i.e., a two-mode squeezed states). Undoing the sympectic transformation yields a covariance matrix for a Gaussian purification of the original state.

More details, including an explicit form for the purification, can be found e.g. in A. Holevo and R. Werner, Evaluating capacities of bosonic Gaussian channels, Phys. Rev. A 63, 032312 (2001).