Proposal of a computational approach for simulating thermal bosonic fields in phase space

When a quantum field is in contact with a thermal bath, the vacuum state of the field may be generalized to a thermal vacuum state, which takes into account the thermal noise. In thermo field dynamics, this is realized by doubling the dimensionality of the Fock space of the system. Interestingly, the representation of thermal noise by means of an augmented space is also found in a distinctly different approach based on the Wigner transform of both the field operators and density matrix, which we pursue here. Specifically, the thermal noise is introduced by augmenting the classical-like Wigner phase space by means of Nos\'e-Hoover chain thermostats, which can be readily simulated on a computer. In this paper, we illustrate how this may be achieved and discuss how non-equilibrium quantum thermal distributions of the field modes can be numerically simulated.

fields [27-34] on a computer, which exploits the Wigner formulation of quantum field theory [35][36][37][38][39][40]. When considering a free field, it is useful to first express the field operators and density matrices in terms of creation and annihilation operators. The creation and annihilation operators are then mapped onto canonical position and momentum operators, which are finally Wigner transformed. When the initial state of the field is pure, the Wigner distribution function of the field is expressed as an infinite product of the Wigner distribution functions [41][42][43] of the single field modes. This is also the case for a quantum thermal state of an infinite set of free field modes, i.e., the thermal Wigner function of the field is an infinite product of the thermal Wigner functions of the modes. In practice, it is possible to simulate quantum thermal distributions with a different temperature for each mode and to define a time-dependent temperature. This is achieved by means of a technique called massive Nosé-Hoover chain (NHC) thermostatting [44,45], which, in our case, involves coupling a separate NHC thermostat [44,45] to each mode. Since one requires a method that can be practically implemented on a computer, the Fock space is truncated to become a many-dimensional (but finite) Hilbert space.
It should be noted that, in spite of their theoretical differences, the Wigner representation of the thermal vacuum that is obtained in our approach is analogous to the one given in thermo field dynamics [24][25][26]. As in thermo field dynamics, where the additional degrees of freedom represent the thermal bath, the Wigner phase space of the thermal system is augmented by the NHC degrees of freedom [44,45]. In practice, because each NHC can often be realized with just two Nosé-Hoover thermostats chained to each other, the dimension of the thermal Wigner phase space is tripled.
A Wigner representation of the modes of the field [46] has certain technical and conceptual advantages. First, since each mode is independently thermostatted, it is easy to generate a non-equilibrium quantum thermal distribution of the field by assigning each mode a different thermodynamic temperature. Secondly, there is no technical difficulty in making the mode temperatures time dependent. As will be explained below, this feature can be exploited to simulate the passage from a quantum thermal distribution of the modes to a classical one. The generalization of our formalism to situations where the field is coupled to a spin system [47][48][49][50][51][52][53][54][55][56][57] is left for future work. Such an endeavour could lead to new computational experiments on spin dynamics, which to the best of our knowledge, would be novel. We propose one such computational experiment in the Conclusions.
The paper is organized as follows. In Section II, we show how to represent the dynamics and averages of bosonic field modes in Wigner phase space. In Section III, we discuss how a bosonic thermal state is described within the Wigner representation of quantum field theory.
We also illustrate how to control the temperature of each field mode on a computer. Finally, our conclusions and perspectives are given in Section IV.

II. WIGNER FORMULATION OF BOSONIC FIELD THEORY
Let us consider a quantum bosonic free field, confined to a finite region of space. This field possesses discrete mode frequencies, ω J , and bosonic creation and annihilation operators,â J andâ † J , where J is an integer. The Hamiltonian operator of the free field reads [38,46]: Let us now define the following canonical transformation: whereQ J andP J are Hermitian operators obeying the canonical commutation relation The scaling factor λ J is equal to √ m J ω J when the field consists of oscillators of mass m J or is equal to λ J = √ (ω J /c) in the case of an electromagnetic field [46]. Upon substituting Equations (2) and (3) in Equation (1), we obtain another representation for the discretized quantum free field Hamiltonian, Since any initial state of the field can be expanded in the basis of Fock states, we consider a multimode Fock state. This can be written as: where n J is the occupation number of state J. To each state |Φ J , one can associate a single mode density matrix, i.e.,ρ J = |Φ J Φ J |. Therefore, the total density matrix of the field is an infinite tensor product of the single mode density matrices: We can now define the Wigner transform [41][42][43] of each mode's density matrix as: where d is the dimension of the spatial region in which the field is confined. The Wigner distribution function of the field is therefore given by: where the field-mode coordinates, (Q, P ) = (Q 1 , Q 2 , . . . , P 1 , P 2 , . . .), are c-numbers in Wigner phase space [42,43].
We now consider the quantum averages of generic bosonic field operators,Ô, in the second quantized form: IfÔ is symmetrically ordered with respect to theâ J andâ † J operators [46], that is, where b nm are expansion coefficients, then we can introduce the Wigner transform ofÔ as: Under the restriction for the field observables mentioned above, averages can then be calculated as [46]: where, in practice, in order to obtain the Wigner transform of operatorsÔ obeying the condition mentioned above, one first uses the canonical transformations in Equations (2) and (3) and then replaces the operators (Q J ,P J ) with the c-numbers (Q J , P J ). For example, where µ J = λ 2 J /ω J . Equation (13) is just the c-number version of Equation (4). In the Wigner representation, the quantum vacuum field equations in the Heisenberg picture: become: where X = (X 1 , X 2 , . . . , X J , . . .) (with X J = (Q J , P J )), ∇ = ∂/∂X, and the direction of the arrow indicates the direction in which the operator acts. The matrix B is the constant symplectic matrix [58]: whose specific form gives rise to the Poisson bracket [58] on the right hand side of Equation (16) (where 1, for example, denotes an infinite-dimensional block diagonal matrix of ones). Equation (16) is equivalent to Equations (14) and (15), showing that the quantum dynamics of a harmonic system may be exactly mapped onto a classical-like time evolution. In this case, the quantum character of the system enters the description through the initial conditions.
In the next sections, we will see how the use of the Wigner representation allows one to exploit the many numerical algorithms originating from molecular dynamics simulation in order to propagate the degrees of freedom represented in phase space.

III. COMPUTER SIMULATION OF THERMAL FIELD STATES
The field Hamiltonian in Equation (1) is isomorphic to a Hamiltonian of a collection of non-interacting harmonic oscillators [59][60][61][62]. The thermal Wigner function, W T J , of a harmonic mode has the following analytical form [46]: where:β is a frequency-dependent inverse temperature with β = 1/k B T , and the Hamiltonian H J W (Q J , P J ) of the J th mode is defined in Equation (13). The thermal Wigner function of the field is: Equations (18)- (20) show that, in the Wigner representation, the thermal state of a free field is given in terms of an infinite product of single mode Wigner functions. Each of these single mode Wigner functions is a Boltzmann factor with a frequency-dependent temperature, as defined in Equation (19). Thus, the thermal vacuum state requires each oscillator to have a minimum average energy consistent with the temperature constraints fixed by Equation (19). In a numerical simulation, the number of oscillators in Equation (20) for the thermal state of the free field must be truncated to a finite number N, i.e., The dynamics of this field can be simulated by first sampling the initial (Q J , P J ) from the Boltzmann factors W T J (Q J , P J ) for each of the N modes and then propagating the modes in time by numerically solving Equation (16). Since Equation (16) generates constant energy trajectories, the thermal equilibrium distribution in Equation (21) is an invariant of the energy conserving phase space flow of the free field. It should be noted, however, that the invariance is broken when the field is coupled to another system. The lack of invariance of the thermal distribution arises ultimately from the non-conservation of the free field Hamiltonian in the presence of the coupling.
In addition to the case when the field is coupled to another system [47][48][49][50][51][52][53][54][55][56][57] there are a number of cases when it is desirable to simulate the dynamics of the thermal field state [27-34] on a computer. To this end, we employ a Nosé-Hoover chain (NHC) thermostat [44,45], which may be theoretically defined in terms of a quasi-Lie bracket [63][64][65][66]. In general, NHC thermostats are used to increase the ergodicity of the dynamics of non-ergodic systems, such as a system of harmonic oscillators. In NHC thermostatted dynamics, the phase space point of oscillator J is extended as follows: where ξ J K and χ J K denote the position and momentum, respectively, of the K th thermostat in the chain attached to mode J and n defines the length of the chain. The energy of the NHC thermostat is given by: The temperature control of the physical coordinates, X, of the field by the fictitious NHC phase space coordinates, denoted collectively by (χ, ξ), is realized by solving a set of quasi-Hamiltonian equations of motion (shown below in Equation (28)). The inertial parameters M K (where K = 1, . . . , n) control the speed of the response of the thermostat variables to the imbalance between the instantaneous kinetic energy of each mode, P 2 J /2µ J , and the kinetic energy corresponding to its desired temperature T J , k B T J , where k B denotes the Boltzmann constant. Finally, the Hamiltonian of the field together with the NHC thermostat is: In order to illustrate the theory, for convenience, we set the number of thermostats in the NHC to n = 2. The quasi-Hamiltonian equations of motion are given by: where ∇ = (∂/∂Q 1 , . . . , ∂/∂Q J , ∂/∂ξ J 1 , ∂/∂ξ J 2 , ∂/∂P J , ∂/∂χ J 1 , ∂/∂χ J 2 , . . . , ∂/∂χ N 2 ) and: is the antisymmetric matrix that generalizes the symplectic matrix, defining the quasi-Hamiltonian NHC phase space flow (where each element of the matrix in Equation (26) is an N ×N block diagonal matrix, e.g., P is an N ×N block diagonal matrix containing the P J 's along the diagonal). The right hand side of Equation (25) defines the quasi-Lie bracket [63], which generalizes the Poisson (Lie) bracket to the case of NHC dynamics [44,45]. It should be noted that the coupling between the field and the NHC thermostat does not arise from the Hamiltonian in Equation (24), but from the the matrix B NHC in the quasi-Hamiltonian equations of motion in Equation (26). Finally, the initial thermal Wigner function in the extended space is given by: Its equation of motion reads [67]: where: is the compressibility of the phase space. The emergence of the compressibility of the phase space in Equation (28) Upon integrating by parts, one obtains the corresponding expression in the Schrödinger picture: where the adjoint Wigner-Liouville operator is given by: From a very different perspective, Equations (27)  (which corresponds to the classical limit) in the NHC dynamics generate a nonequilibrium situation, from which the field would ultimately reach its classical thermal state. In a situation where the field is no longer isolated, it would be interesting to investigate how the time dependent classical limit of the field would modify the dynamics of the coupled system.
Moreover, the classical limit may be applied to groups of modes so that one can study how each group affects the dynamics of the coupled system.

IV. CONCLUSIONS
In this paper, we discussed how pure and thermal states of free bosonic fields [27-34] can be represented theoretically and simulated computationally. The simulation protocol, which involves an NHC thermostat [44,45] coupled to each field mode, propagates the dynamics of a thermal state living in an extended Wigner space [67]. This is conceptually similar to what is done in thermo field dynamics [24][25][26].
For simplicity, the theory was only applied to the free field case [38], in order to demonstrate the technicalities associated with coupling a different NHC thermostat [44,45] to each field mode. However, this theory is intended for situations where the field is coupled to a spin system [47][48][49][50][51][52][53][54][55][56][57]. In this case, the use of NHC thermostats in the dynamics of the thermal state of the field would make it possible to simulate processes that, to our knowledge, have not been investigated so far. For example, NHC dynamics could allow one to simulate the transition from a vacuum field state to a thermal vacuum state (or to thermally excited states of the field). Then, one could study changes in the spin system's transport properties in response to various field thermal processes. Such applications are left for future work.