# Simulation of Nuclear Quantum Effects in Condensed Matter Systems via Quantum Baths
## 4 Critical Analysis of Quantum Baths: Model Systems
### Harmonic Systems  
A pedagogical introduction to the QTB has been given by [PERSON] and [PERSON], with several useful remarks, in particular, concerning the application of the QTB to a harmonic system, \(V(x)=\frac{1}{2}m\omega_{0}^{2}x^{2}\)[39]. In that case, as the force is linear with position, an explicit expression can be obtained for the Fourier transform of the position and the velocity and hence for their respective correlation spectra:  
\[C_{xx}(\omega)=\frac{2\gamma\ \theta(\omega,T)/m}{(\omega^{2}-\omega_{0}^{2})^{ 2}+\omega^{2}\gamma^{2}}\ \ \text{and}\ \ C_{vv}(\omega)=\frac{2\gamma\ \omega^{2}\theta(\omega,T)/m}{(\omega^{2}-\omega_{0}^{2})^{2}+\omega^{2} \gamma^{2}} \tag{19}\]  
These spectra correspond to Lorentzian peaks centered at \(\omega_{0}\), and with a width \(\gamma\). They are good approximations to the quantum symmetrized correlation spectra \(C_{xx}^{s}(\omega)\) and \(C_{vv}^{s}(\omega)\) apart from a spectral broadening with respect to the Dirac \(\delta\)--function expected for the 1D quantum system, which is only recovered in the \(\gamma\to 0\) limit. The integrals of these two spectra provide, respectively, the average potential and kinetic energies as:  
\[\langle E_{pot}\rangle = \int\frac{d\omega}{2\pi}\ \frac{\gamma\ \omega_{0}^{2}}{(\omega^{2}-\omega_{0}^{2})^{2}+\omega^{2}\gamma^{2}}\ \theta(\omega,T) \tag{20}\] \[\langle E_{kin}\rangle = \int\frac{d\omega}{2\pi}\ \frac{\gamma\ \omega^{2}}{(\omega^{2}- \omega_{0}^{2})^{2}+\omega^{2}\gamma^{2}}\ \theta(\omega,T) \tag{21}\]  
As the zero-point energy term in \(\theta(\omega,T)\) is linear for large \(\omega\), the integral in (21) diverges logarithmically. [PERSON] and [PERSON] therefore suggested the introduction of a cutoff frequency \(\omega_{max}\) in the random force spectrum, that should be larger than the highest frequency of the physical system. (Note that a finite integration time step \(\delta t\) implicitly imposes a maximal frequency \(\delta t^{-1}/2\). However, the use of a noise power spectrum without an explicit high-frequency cutoff \(\omega_{max}\) poses various numerical problems). In Appendix A, the practical noise generation as well as the choice of the algorithms to integrate the Langevin equation is briefly discussed. With this precaution, in the limit of small values of \(\gamma\), the Lorentzian factors in Equations (20) and (21) tend to Dirac \(\delta\)--functions and both the kinetic and the potential energy equal \(\theta(\omega_{0},T)/2\) as expected for the quantum harmonic oscillator (for nonzero \(\gamma\), the average energies are slightly different from \(\theta(\omega_{0},T)/2\), although this difference can be corrected using spectral deconvolution techniques [18; 40]). The position and the momentum probability distributions obtained in QTB simulations of the harmonic oscillator are Gaussian, with widths fixed by Equations (20) and (21), therefore the method provides exact estimates (at least in the \(\gamma\to 0\) limit) for the quantum average of any static observable depending only on position or momentum, including zero-point energy effects. As mentioned above, it also yields a good approximation to some dynamical properties such as position and velocity autocorrelation spectra. This distinction is important in the context of nuclear quantum effects simulations as path-integral techniques can provide exact references for static properties (though at an elevated computational cost, see Section 6), but only approximate methods are available to simulate the quantum dynamical properties of complex systems with large numbers of atoms. The paragraphs below examine the effect of anharmonicity on the accuracy of QTB simulations.