# Simulation of Nuclear Quantum Effects in Condensed Matter Systems via Quantum Baths
## 3 Enforcing Quantum Statistics via Generalized Bath
### The Quantum Thermostat  
Independently from the work by [PERSON] and coworkers, [PERSON], [PERSON] and [PERSON] also proposed to approximate NQEs via a non-Markovian GLE. In two papers in 2009 [13; 14], they introduced a general framework to implement and propagate GLEs of the form of (6) using auxiliary variables. The authors aim at thermalizing different degrees of freedom with distinct characteristic frequencies at different effective temperatures. This tool can be used for different purposes: in particular, the authors exploit it to build a quantum thermostat (QT) and thermalize vibrational modes with an effective temperature that includes a zero-point motion contribution [14], similarly to the QTB. They also propose to apply this method in the context of Car-Parrinello molecular dynamics [33], where the ionic degrees of freedom should evolve at a given (physical) temperature while the faster electronic degrees of freedom should oscillate around the ground state that corresponds  
Figure 2: \(\theta(\omega,T)\) (left panel) and its temperature derivative \(d\theta(\omega,T)/dT\) (right panel) for three representative frequencies \(\
u=\omega/2\pi\). \(v_{1}=4\) THz is typical of acoustic modes in crystals; \(v_{2}=20\) THz is in the range of optical modes in crystals formed by light mass elements, as well as of OH and CH bending modes; \(v_{3}=100\) THz could correspond to C-H and O-H stretching modes. For the latter frequency, one can see that the classical behavior is never reached for \(T\leq 1000\) K.  
to the instantaneous ionic configuration, at much smaller temperatures [13]. The authors successfully applied the colored bath thermostat on a model polarizable system and verified that the fast degrees of freedom (i.e., the motion of the centers of mass of the shells with respect to their nuclei) are correctly thermalized. More generally, they showed that the colored bath thermostat can be adapted to specific cases once the density of vibrational states is known and that it allows a faster thermalization than the Nose-Hoover chains [13; 22]. Starting from the Orstein-Uhlenbeck theory of stochastic processes, [PERSON] and coworkers show how a non-Markovian GLE can be implemented via a Markovian Langevin dynamics for an extended set of variables of which two are the physical position and momentum (for a one-dimensional physical system) and \(N\) are additional auxiliary momenta [34]. In this way, the colored thermostat GLE can be efficiently applied by integrating the Markovian dynamics in the \((N+2)\)-dimensional space via usual molecular dynamics techniques. This dynamics is governed by the following equations:  
\[ \begin{bmatrix}\dot{\pi}\\ \dot{p}\\ \dot{\mathbf{s}}\end{bmatrix}  = - \begin{bmatrix}0&-1/m&\mathbf{0}\\ V^{\prime}(x)&a_{pp}&\mathbf{a_{p}^{T}}\\ \mathbf{0}&\mathbf{\bar{a}_{p}}&\mathbf{A}\end{bmatrix}\begin{bmatrix}x\\ p\\ \mathbf{s}\end{bmatrix} + \begin{bmatrix}0&0&\mathbf{0}\\ \mathbf{0}&b_{pp}&\mathbf{b_{p}^{T}}\\ \mathbf{0}&\mathbf{\bar{b}_{p}}&\mathbf{B}\end{bmatrix}\begin{bmatrix}0\\ r_{p}\\ \mathbf{r}\end{bmatrix}  \tag{16}\]  
where \(\mathbf{s}\) and \(\mathbf{\hat{s}}\) are the \(N\)-component vectors of the auxiliary momenta and their time derivatives, \(\mathbf{A}\) and \(\mathbf{B}\) are \((N\times N)\) matrices of adjustable parameters that characterize respectively the friction and the random force memory kernels, together with the \(N\)-component parameter vectors \(\mathbf{a_{p}}\), \(\mathbf{\bar{a}_{p}}\) and \(\mathbf{b_{p}}\), \(\mathbf{\bar{b}_{p}}\). \(V^{\prime}(x)\) is the spatial derivative of the potential energy and the \((N+1)\)-component vector \((r_{p},\mathbf{r})\) contains normalized white noise variables: \(\langle r_{i}(t)r_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})\). The dynamics in the extended \((x,p,\mathbf{s})\)-space is therefore Markovian and it can be solved analytically for a harmonic potential \(V(x)=m\omega^{2}x^{2}/2\).