The state vector | ψ ( t ) 〉 of an isolated system evolves according to the Schr o ̈ dinger equation, or the corresponding density operator ρ ^ ( t ) = | ψ ( t ) 〉 〈 ψ ( t ) | obeys the quantum Liouville equation [12]. i ℏ d ρ ^ t d t = H ^ , ρ ^ t , (1) where H ^ is the Hamiltonian of the system. Equation 1 holds a fundamental position in quantum mechanics and serves as the starting point for studying quantum master equations of open systems. The solution of the abovementioned equation is ρ ^ t = e − i ℏ H ^ t ρ ^ 0 e i ℏ H ^ t , (2) where ρ ^ ( 0 ) = | ψ ( 0 ) 〉 〈 ψ ( 0 ) | is the initial density matrix. In the eigenstate basis of the Hamiltonian { | ϕ n 〉 } , Equation 2 can be expressed as follows: ρ ^ t = ∑ n , m c n c m * e − i ℏ E n − E m t | ϕ n 〉 〈 ϕ m | , where the coefficient c n = ⟨ ϕ n | ψ ( 0 ) ⟩ , the asterisk denotes complex conjugation, and E n is the eigenvalue corresponding to eigenstate | ϕ n 〉 . The average value of any observable O ^ evolves as follows: ⟨ O ^ ⟩ = T r ρ ^ t O ^ = ∑ n , m c n c m * e − i ℏ E n − E m t O m n , where O m n = ⟨ ϕ m | O ^ | ϕ n ⟩ is a matrix element of the observable O in the energy basis.

It is evident that a generic pure state cannot exhibit the thermal equilibrium behavior. The time-independent diagonal terms depend on the initial state’s coefficients, while the off-diagonal terms exhibit harmonic time dependence that generally does not decay to 0. On the other hand, the thermodynamic entropy of a system (commonly referred to as Boltzmann or Gibbs entropy in statistical physics) vanishes for a pure state but is positive for a mixed state. From the perspective of quantum mechanics, all thermodynamic systems—including so-called isolated systems—are effectively open and must be described by mixed states, thus possessing positive thermodynamic entropy [13]. However, an isolated quantum system undergoing unitary evolution remains in a time-dependent pure state, implying that its global thermodynamic entropy is rigorously 0. This observation reinforces the conclusion that a genuinely isolated system cannot thermalize globally. Nevertheless, as discussed in Section 3.1, a subsystem of such an isolated system can undergo thermalization by treating the remainder as an effective heat reservoir at fixed energy [14].

When an isolated system interacts with an external environment, it becomes an open system. Due to the entanglement between the system and its environment, the system’s evolution becomes non-unitary, and its dynamics are governed by the quantum master equation. The combined system–environment pair S + E forms a composite isolated system, with S itself acting as a subsystem of S + E .

The total Hamiltonian H ^ of S + E can be expressed as follows: H ^ = H ^ S ⊗ I E + I S ⊗ H ^ E + H ^ I , (3) where “ ⊗ ” represents the direct product; H ^ S and H ^ E are the Hamiltonians of S and E , respectively; H ^ I represents the interaction between S and E ; and I S and I E are the identity operators in the Hilbert spaces of S and E , respectively.

Let O ^ S be an observable referring to S ; it should be expressed as O ^ = O ^ S ⊗ I E when considering S as a subsystem of S + E . Let ρ ^ ( t ) and ρ ^ S ( t ) be the density matrices of S + E and S at time t , respectively. Then, physical consistency requires that measurement averages be the same whether computed via ρ ^ ( t ) or ρ ^ S ( t ) , ⟨ O ^ ⟩ = T r ρ ^ t O ^ S ⊗ I E = T r S ρ ^ S t O ^ S , (4) where ρ ^ S t = T r E ρ ^ t

is also known as the reduced density matrix of S , T r S denotes the trace over the Hilbert space of S , and T r E denotes the partial trace over the degrees of freedom of E . The time evolution of the reduced density matrix ρ ^ ( t ) is derived by taking the partial trace over the degrees of freedom of E on both sides of the quantum Liouville equation as follows: i ℏ d ρ ^ S t d t = T r E H ^ , ρ ^ t , where H ^ is given by Equation 3. Starting from the abovementioned equation, the quantum master equation in Lindblad form—a first-order linear differential equation for the reduced density matrix—can be derived under the weak-coupling limit, Born approximation, Markov approximation, and rotating wave approximation [15]. d ρ ^ S t d t = − i ℏ H ^ S , ρ ^ S t + D ρ ^ S t , (5) where D ( ρ ^ S ( t ) ) is the dissipator, which describes the possible transitions, dissipations, and decoherence of the system due to the interaction between the system and the environment.

Assume that the environment is in a stationary state; then, ρ ^ E is invariant under the reservoir dynamics. Furthermore, if the environment is in thermal equilibrium with inverse temperature β , the Lindblad equation predicts that the density of the matrix of the canonical ensemble, ρ ^ t h = e − β H ^ S T r e − β H ^ S , (6)

is a stationary solution of the quantum master (Equation 5) [15]. This means that regardless of the initial state, S will asymptotically approach thermodynamic equilibrium after interacting with R . Equation 6 can also be directly derived from traditional statistical physics: any smaller system coupled to a large thermal reservoir will eventually thermalize, and its density matrix adopts the canonical ensemble form. In this study, we present the density matrix in this way to demonstrate how the temperature-dependent pure state for a composite isolated system consisting of a thermodynamic system and its heat reservoir can be constructed from first principles.

The abovementioned discussion indicates that the thermal equilibrium density matrix of S commutes with the Hamiltonian, indicating that they share common eigenstates. Consequently, the equilibrium density matrix of S is diagonal in the energy eigenbasis, ρ ^ S = ∑ n ρ n | ϕ n 〉 〈 ϕ n | , (7) where ρ n = e − β E n Z (8)

is the probability that the system is in the energy eigenstate | ϕ n 〉 of S , where Z = ∑ m e − β E m is referred to as the partition function [16].

According to the quantum Liouville equation, the composite isolated system S + R is always in a pure state. Let { | φ n 〉 } and { | Φ μ 〉 } be any fixed orthonormal bases for S and R , respectively; then, a pure state of S + R at any time can be expanded as follows [17]: | ψ 〉 = ∑ n , μ a n , μ | φ n 〉 ⊗ | Φ μ 〉 , (9) where the coefficient a n , μ = 〈 φ n | ⊗ ⟨ Φ μ | ψ ⟩ .

As established earlier, in the long-time limit, S as a subsystem will thermalize such that its density matrix and Hamiltonian commute with each other and, hence, share common eigenstates. Consequently, the eigenstates of the density matrix appearing in the Schmidt decomposition are also the eigenstates of the Hamiltonian. As we shall demonstrate, this property extends to R . Let | φ n 〉 correspond to the Hamiltonian eigenstate | ϕ n 〉 of S ; this allows us to express | ψ 〉 in Equation 9 in the form of Schmidt decomposition [15]: | ψ 〉 = ∑ n ρ n 1 / 2 | ϕ n 〉 ⊗ | Φ ̃ n 〉 , (10) where ρ n is calculated using Equation 8, and | Φ ̃ n 〉 is defined as follows: | Φ ̃ n 〉 = ∑ μ ρ n − 1 / 2 a n , μ | Φ μ 〉 . (11)

It can be proven that the set | Φ ̃ n 〉 of vectors with index n is orthonormal, ⟨ Φ ̃ m | Φ ̃ n ⟩ = δ m n ,

and it is the set of eigenstates of the following reduced density matrix for R , ρ ^ R β = T r S ρ ^ S R = T r S | ψ 〉 〈 ψ | = ∑ n ρ n | Φ ̃ n 〉 〈 Φ ̃ n | .

It should be emphasized that if { | ϕ n 〉 } is not the eigenstate basis of the Hamiltonian and density matrix of S , the basis introduced in Equation 11 may not be necessarily orthonormal or even not orthogonal in general [18].

Since the temperature parameter β appears on the right side of Equation 10, | ψ 〉 is now a pure state with temperature. Denoting | ϕ n 〉 by | n 〉 and | Φ ̃ n 〉 by | n ̃ 〉 , respectively, Equation 10 can be rewritten as follows: | ψ β 〉 = ∑ n ρ n 1 / 2 | n 〉 ⊗ | n ̃ 〉 . (12)

The corresponding density matrix of the temperature-dependent pure state ( | ψ ( β ) 〉 ) is ρ ^ S R β = | ψ β 〉 ⟨ ψ β | = ∑ n , m ρ n 1 / 2 ρ m 1 / 2 | n ⟩ ⊗ | n ̃ 〉 〈 m | ⊗ 〈 m ̃ | .

It is easy to check that the reduced density of S is ρ ^ S β = T r R ρ ^ S R β = ∑ n ρ n | n 〉 〈 n | .

As shown in Section 2.2, for an observable O ^ S of S , it should be expressed as O ^ S ⊗ I R when considering S as a subsystem of S + R . According to Equation 4, the expectation value of the pure state of S + R coincides with the statistical average of a mixed state of S , T r ρ ^ S R β O ^ S ⊗ I R = T r S ρ ^ S β O ^ S .

In the energy representation, ρ ^ S ( β ) is obtained using Equation 7. Omitting the identity operator I R and the subscript “ S ,” we finally obtain 〈 ψ β | O ^ | ψ β 〉 = ∑ n ρ n 〈 n | O ^ | n 〉 . (13)

In this way, a question in quantum field theory realized through a temperature-dependent pure state and a question in statistical mechanics concerning a subsystem in thermal equilibrium can be converted into each other.

As demonstrated in Section 2.3, a system interacting with a heat reservoir will undergo thermalization in the long-time limit. However, this framework requires the heat reservoir to be predefined as part of a composite isolated system. It seems that this prior assumption introduces limitations in constructing | ψ ( β ) 〉 from the first principles.

Recent studies, however, demonstrate that for an isolated system satisfying the ETH, its subsystem can indeed thermalize [19–22]. The ETH serves as the primary conceptual framework for understanding how quantum mechanics leads to thermalization. According to the ETH, each eigenstate of the Hamiltonian of the composite system S + E always implicitly contains a thermal state; consequently, the reduced density matrix of such an eigenstate becomes indistinguishable from a canonical ensemble density matrix [23]. This yields the following relation: T r E | ϕ k 〉 〈 ϕ k | = e − β k H ^ S T r e − β k H ^ S ,

and the temperature parameter β k is determined using the following constraint: T r e − β k H ^ S T r e − β k H ^ S H ^ S = E k , where E k is the eigenvalue corresponding to the eigenstate | ϕ k 〉 of the S + E Hamiltonian. This implies that under unitary evolution of the full system, a subsystem can exhibit thermal behavior when its complement acts as an effective heat reservoir—a key assumption in Section 2.3. Notably, the ETH has not only been theoretically predicted but also experimentally verified [24].

The microscopic origin of temperature might still be considered an open question. Although the ETH has been widely validated as a sufficient condition for thermalization in many finite-sized systems through numerical simulations [14], it has not been universally proven to be a necessary condition. Due to the presence of local conserved quantities that prevent thermalization, there exist systems that do not thermalize at all, such as integrable systems or many-body localized states [25].

As shown in Equation 12, the eigenstates | n 〉 and | n ̃ 〉 appear as a pair in | ψ ( β ) 〉 . Consequently, when calculating thermodynamic quantities using the expectation value method, the degrees of freedom must be doubled compared to those in the statistical ensemble approach. Takahashi et al. termed | ψ ( β ) 〉 the “thermal vacuum state” and designated the tilde system as the “fictitious system” [7]. They further interpreted particles in the fictitious system as “holes.” Other scholars have proposed alternative interpretations, such as the “mirror of the real world” [8], “complement to the dissipating system” [10], and “thermal ghost” [11].

Equation 13 serves as the foundational starting point of TFD. Notably, in TFD, the density matrix ∑ n ρ n | n 〉 〈 n | is assumed a priori, whereas the pure state | ψ ( β ) 〉 must be explicitly constructed. This construction corresponds to the purification of the thermodynamic system’s mixed state. However, purification is often viewed as a mathematical procedure with the auxiliary reference system (the fictitious system) lacking direct physical significance [17]. It can be shown that there exist different pure states satisfying the given condition, all are related to one another by unitary transformations. Thus, for the sole purpose of computing thermodynamic quantities, the “fictitious system” in TFD requires no physical interpretation. However, diverse interpretations remain possible precisely because none affect computational outcomes, leading to the current plurality of explanations.

If | ψ ( β ) 〉 is determined as discussed in Section 2.4, the reduced density matrix ∑ n ρ n | n 〉 〈 n | becomes uniquely defined. The tilde system representing the heat reservoir in this framework directly corresponds to the “fictitious system” in TFD. Therefore, assigning a physical interpretation to the “fictitious system” is most logically consistent when identifying it as the system’s heat reservoir. This perspective elevates the purification of the mixed state beyond a mere mathematical tool, endowing it with physical meaning.

This reinterpretation of the “fictitious system” may prompt a re-examination of previously controversial or even meaningless questions. For instance, Takahashi et al. argued that a “one particle state” is built up from the thermal equilibrium state either by adding a particle to the physical system or eliminating a particle from the “fictitious system,” thereby interpreting a particle in the “fictitious system” as a “hole” of a physical particle. In our framework, however, this corresponds to the physical equivalence between annihilating a particle in the heat reservoir and creating a particle in the system.

As a further example, we can now analyze the quantum entanglement between the physical system and the so-called “fictitious system”—a question previously devoid of physical significance. The Schmidt number for the state | ψ ( β ) 〉 is greater than one, which indicates that | ψ ( β ) 〉 is necessarily an entangled state. The quantity that measures the quantum entanglement between two subsystems in a composite isolated system is the entanglement entropy. Usually, the entanglement entropy is defined as the von Neumann entropy corresponding to any reduced density matrix of the composite isolated system [26]. For instance, for a free particle of mass m in a cube with a side length of L at temperature T , the entanglement entropy between the particle and its heat reservoir is given by S = − T r ρ ^ S l n ρ ^ S = 3 2 ln 2 π m k T L 2 h 2 + 1 , where ρ ^ S is given by Equation 7 with the partition function Z = L 3 ( 2 π m k T h 2 ) 3 / 2 .

It is essential to clarify that for a composite isolated system consisting of two subsystems in thermal equilibrium, the vanishing thermodynamic entropy of the global pure state does not contradict the non-zero thermodynamic entropy of its subsystems. This arises because while the quantum information is preserved, it becomes hidden locally throughout the system and is only accessible through measurements of global observables that do not thermalize [16].

As an entangled state, | ψ ( β ) 〉 is represented as a linear superposition of the direct product states | n 〉 ⊗ | n ̃ 〉 ( n = 0,1,2 , ⋯   ) . The common eigenstate | n 〉 of the density matrix and Hamiltonian of S appear as one component state in every product state, while the common eigenstate | n ̃ 〉 of the density matrix and Hamiltonian of R appear as the other component state. Crucially, the component state is also the eigenstate of any mechanical operator that commutes with both the density operator and the Hamiltonian. Specifically, if the particle number operator satisfies this commutation property, the component state within the direct product state is also the eigenstate of the particle number operator. This structure allows | ψ ( β ) 〉 to be dealt with directly using the methods and techniques of quantum field theory since the particle number operator and its eigenstates play a fundamental role in the usual quantum field theory.

To illustrate the construction of such a state, consider a system of free bosons. Since the Hamiltonian commutes with the particle number operator, the density operator, Hamiltonian, and particle number operator share common eigenstates. Let a k † and a k denote the creation and annihilation operators for a free boson with energy ϵ k , respectively; the orthonormal basis vectors of the Fock space are then expressed as follows: | n 〉 = 1 n ! a k † n | 0 〉 , n = 0,1,2 , … ,

and the corresponding eigenvalues of the particle number operator, density operator, and Hamiltonian are n , ρ n , and n ϵ k ( n = 0,1,2 , ⋯   ) , respectively.

On introducing the eigenstates of the particle number operator of R , characterized by a tilde, | ψ ( β ) 〉 can be constructed as follows: | ψ β 〉 = ∑ n = 0 + ∞ ρ n 1 / 2 | n 〉 ⊗ | n ̃ 〉 = ∑ n = 0 + ∞ 1 − e − β ϵ k e − β n ϵ k 1 / 2 | n 〉 ⊗ | n ̃ 〉 .

Note that the identity operator in the Hilbert spaces of R is I R = ∑ n = 0 + ∞ | n ̃ 〉 〈 n ̃ | , and the average particle number is obtained as follows: f k − B = T r ρ ^ a k † a k ⊗ I R , = ∑ l , n , m = 0 + ∞ 〈 l | ⊗ 〈 l ̃ | ρ n 1 / 2 ρ m 1 / 2 | n 〉 ⊗ | n ̃ 〉 〈 m | ⊗ 〈 m ̃ | a k † a k ⊗ I R | l 〉 ⊗ | l ̃ 〉 , = ∑ l = 0 + ∞ 1 − e − β ϵ k l ⋅ e − β l ϵ k = 1 e β ϵ k − 1 ,

which is the well-known Bose distribution.

It is important to note that for analytically solvable problems, the pure-state expectation value method does not inherently simplify thermodynamic calculations. However, for systems requiring approximate treatments, this formulation provides novel perspectives and techniques for studying finite-temperature phenomena through quantum field theory [1].