# Evaluation of Density Matrix and Helmholtz Free Energy for Harmonic Oscillator Asymmetric Potential via Feynmans Approach

Authors: Piyarut Moonsri; Artit Hutem
Journal: Advances in Physical Chemistry
 (2014)

---

## Abstract

We apply a Feynmans technique for calculation of a canonical density matrix of a single particle under harmonic oscillator asymmetric potential and solving the Bloch equation of the statistical mechanics system. The density matrix(P^u) and kinetic energy per unit length (τ^L) can be directly evaluated from the solving solutions. From the evaluation, it was found that both of the density matrix and kinetic energy per unit length depended on the parameter of the value of asymmetric potential (λ), the value of axes-shift potential (g), and temperature (T). Comparison of the Helmholtz free energy was derived by the Feynmans technique and the path-integral method. The results illustrated are slightly different.


---

## Body

## 1. Introduction

We consider a quantum mechanical ensemble of identical systems in situations where the description is incomplete. When the description of the system is incomplete, the state of the system is described by means of what is called a density matrix. In the statistical point of view, a state of quantum ensemble is described by a density operator (or statistical operator)P^u; the corresponding matrix is called a density matrix. A knowledge of this matrix enables us to calculate the mean value of physical quantity of the system and also the probabilities of various values of such a quantity.March and Murray used the Bloch equation to calculate the canonical density matrix in a single particle framework which may be related directly to the generalized canonical density matrix, containing the Fermi-Dirac function [1]. Howard et al. used the Bloch equation evaluate the electron density P^u and kinetic energy τ^L(x) for the one dimensional sech2(x) potential. Howard et al. used tool employed is the Slater sum, which satisfies a partial differential equation [2].For calculation of the canonical density matrix in terms of the Slater sum and kinetic energy per unit length in terms of the Slater sum see March and Nieto [3]. To study problems of the canonical ensemble, we focus on the calculation the canonical density matrix and the Helmholtz free energy of a single particle under asymmetric harmonic oscillator potential.The scheme of the paper is as follows. In Section2, there is detailing with the density matrix. In Section 3, we show the calculation of the canonical density matrix (or charge bond-order matrix) and the Helmholtz free energy by Feynmans approach. In Section 4, we representation of the evaluation of the density matrix and the Helmholtz free energy by the path-integral method (κ(y)) method. In Section 5 we make a presentation of numerical results for the Helmholtz free energy of a single particle under asymmetric harmonic oscillator potential. The conclusion and discussion are given in Section 6.


## 2. Statistical Density Matrix

A density matrix is a matrix that describes a quantum system in a mixed state, a statistical ensemble of several quantum states. This should be contrasted with a single state vector that describes a quantum system in a pure state. The density matrix is the quantum-mechanical analogue to a phase-space probability measure (probability distribution of position and momentum) in classical statistical mechanics. Explicitly, suppose a quantum system may be found in stateψ1 with probability P1 or it may be found in state |ψ2〉 with probability P2 or it may be found in state |ψ3〉 with probability P3, and so on. The density matrix for this system is
(1)P^=∑n‍Pnψnψn,
where (|ψn〉) need not be orthogonal and ∑nPn‍=1. By choosing an orthonormal basis (|uj〉), one may resolve the density operator into the density matrix, whose elements are
(2)P^i,j=∑n‍Pnuj∣ψnψn∣ui.
For an operator A^ (which describes an observable A^ of the system), the expectation value 〈A^〉 is given by
(3)〈A^〉=∑n‍PnψnA^ψn=∑i‍uiP^A^ui=tr⁡(P^A^).
In the canonical ensemble, the probability that the system is in state |n〉 is given by
(4)Pn=e-βEn∑n‍eβEn,
and the thermal average of an arbitrary operator A^ is
(5)A^=∑n‍nA^ne-βEn∑n‍eβEn=1QNβ∑n‍nA^ne-βEn.
The operator
(6)P^=1QN(β)∑n‍nne-βEn
is commonly referred to as density operator and its matrix representation as the density matrix. Here, QN(β) is the partition function.


## 3. A Linear Harmonic Oscillator Asymmetric Potential by Feynmans Approach

Next, we consider the case of a linear harmonics oscillator of asymmetric potential, whose Hamiltonian is given by [4]
(7)H^S=P^22m+12mω2x^2+mωℏℏω2λx^-gℏmω;
we note that this is a time-independent Hamiltonian of harmonic oscillator asymmetric axes-shift potential system. Substituting in the time-independent Hamiltonian equation (7) into the Bloch equation [3], we can rewrite the differential equation completely in terms of P^u(ξ,ξ´,f) as
(8)-∂P^u∂β=-ℏ22m∂2P^u∂x2+12mω2x^2P^u+mωℏℏω2λx^-gℏmωP^u.
We define some new dimensionless variable. The position variable x^=x is replaced with the dimensionless variable ξ:
(9)ξ=mωℏx,x2=ℏmωξ2,β=2fℏω.
Substituting in for x in terms of ξ, we can rewrite the differential equation completely in terms of ξ and P^u as
(10)-∂P^u∂f=-∂2P^u∂ξ2+ξ2P^u+λξP^u-λgP^u.
Knowing the Gaussian dependence of the density matrix P^u in case of the free particle, we try the following form:
(11)P^uξ,ξ´,f≈e-(a(f)ξ2+b(f)ξ+c(f)),λgP^u=λge-(a(f)ξ2+b(f)ξ+c(f)),-∂P^u∂f=ξ2dafdf+ξdbfdf+dcfdf×e-(a(f)ξ2+b(f)ξ+c(f)),-∂2P^u∂ξ2=2af-4a2fξ2-4afbfξ-b2f×e-(a(f)ξ2+b(f)ξ+c(f)),ξ2P^u=ξ2e-(a(f)ξ2+b(f)ξ+c(f)),λξP^u=λξe-(a(f)ξ2+b(f)ξ+c(f)).
Substituting expansion equation (11) into (10) leads us to the equation
(12)da(f)dfξ2+db(f)dfξ+dc(f)df=1-4a2fξ2+λ-4afbfξ+2af-λg-b2f.
Equating the coefficients for ξ2, ξ, and ξ0=1 on both sides of the last equation (12) gives
(13)da(f)df=(1-4a2(f)),(14)db(f)df=λ-4afbf,(15)dc(f)df=(2a(f)-λg-b2(f)),
respectively. Computing (13) is integrated to give
(16)a(f)=12coth⁡(2f).Equation (14) is the linear first-order ordinary differential equation to give [5, 6]
(17)λ=db(f)df+2coth⁡(2f)b(f),bf=e-∫‍2coth⁡⁡2fdf∫‍λe∫‍2coth⁡⁡2fdfdf+A,b(f)=λ2coth⁡(2f)+Acsch2f.
Equation (15) is integrated with respect to f as
(18)dcfdf=coth⁡2f-λg-λ2coth⁡2f+Acsch2f2,cf=∫coth⁡2fdf-∫λgdfcf=-∫λ2coth⁡2f+Acsch2f2df,cf=18-2fλ4g+λ+4A2+λ2coth⁡2f+4ln⁡⁡sinh⁡2f+2Aλcschfsech⁡f-ln⁡⁡B4A2+λ2,
where B, λ, and g are constants. Substituting a(f), b(f), and c(f) of (16), (17), and (18), respectively, into (11), we may rewrite density matrix as
(19)P^u=Bsinh⁡⁡2f=×exp⁡⁡-ξ22+ξλ2+184A2+λ2coth⁡⁡2f+ξA+Aλ2csch2f-fλ4(4g+λ)ξ22+ξλ2+184A2+λ2-ξ22+ξλ2+184A2+λ2coth⁡⁡2f.
Thus, the Taylor series expansion of sinh⁡(2f), coth⁡(2f) about f→0 is given by
(20)lim⁡f→0sinh⁡2f≈2f,lim⁡f→0coth⁡2f≈12f;
we obtain the density matrix final form:
(21)P^u≈B2f×exp⁡⁡-ξ22+ξλ2+184A2+λ212f+ξA+Aλ212f-λf4(4g+λ).
The initial condition is
(22)P^u=δx-x´forf=0
or
(23)P^u=mωℏδξ-ξ´forf=0.
For low f=ℏω/2β=ℏω/2kBT (high temperature) the particle should act almost like a free particle, as its probable kinetic energy is so high. Therefore, we expect that, for low f, the density matrix for a harmonic oscillator will be given approximately by [7]
(24)P^ufre(ξ,ξ´,f)≈mω4πℏfe-ξ-ξ´2/4f.Comparing the value ofA, B in (21) with (24) identifies the solution as
(25)B2f=mω4πℏf⟶B=mω2πℏ,ξ2-2ξξ´+ξ´2=ξ2+λ+2Aξ+A2+λ24+Aλ-4gλf2-λ2f2,A=-ξ´+λ2.Substituting (25) into (19) and setting ξ=ξ´ produce
(26)P^uξ,f≈mω2πℏsinh⁡⁡2f×e-[ξ+λ/22cosh⁡⁡2f-1/sinh⁡⁡2f-λf/4(λ+4g)].
Using the identity of trigonometric tanh(x/2)=cosh⁡(x)-1/sinh⁡(x) and setting ξ=mω/ℏx into (26), we can rewrite the canonical density matrix or charge bond-order matrix [8] as
(27)P^ux,f≈mω2πℏsinh⁡⁡2f×e-mω/4ℏ[(2x+ℏ/mωλ)2tanh(f)-λfℏ/mω(λ+4g)].The obtained results are neatly summarized in Figure1. It is the probability density for finding the system at x; we note as a Gaussian distribution in x. Let us define four new the variables γ, η, μ, and ϕ as function of f by
(28)γ(f)=mω2πℏsinh⁡(2f),η(f)=mω4ℏtanh(f),μf=λf4λ+4g,ϕ=ℏmωλ;
we can rewrite the density matrix completely in terms of the new variable (γ,η,μ,ϕ) as
(29)P^u(x,f)≈γe-(η(2x+ϕ)2-μ).
The expectation value of x2 in the density matrix P^u(x,f) can be determined from the definition of the expectation values as
(30)x2=∫-∞∞‍x^2P^u(x,f)dx∫-∞∞‍P^u(x,f)dx=γπeμ16η3/22ηγπeμ(1+2ηϕ2),〈x2〉=ℏ2mωcoth⁡(ℏωβ2)+ℏ4mωλ2.
The average of the potential energy (Es) is thus
(31)Es=12mω2x2,〈Es〉=ℏω4coth⁡ℏωβ2+ℏω8λ2.Plot of the canonical density matrix of a single particle for harmonic oscillator asymmetric potential system.
(a)
Density Matrix,λ=1, g=1
(b)
Density Matrix,g=1, T=15
(c)
Density Matrix,λ=1, g=+, T=15
(d)
Density Matrix,λ=1, g=+, T=15
(e)
Density Matrix,g=1, T=15The kinetic energy per unit lengthτ^L(x,f) is defined from the canonical density matrix P^u(x,f) as [3]
(32)τ^L(x,f)=-ℏ22m∂2∂x2P^u(x,f).
By substituting (29) into (32), we can write the kinetic energy per unit length completely as follows
(33)τ^Lx,T=ℏω2tanh⁡ℏω2kBTmω2πℏsinh⁡⁡ℏω/kBT×eλℏω/8kBTλ+4g×2-mωℏtanh⁡ℏω2kBT2x+ℏmωλ2×e-mω/4ℏtanh⁡(ℏω/2kBT)(2x+ℏ/mωλ)2.
The obtained results from (33) are summarized in Figure 2.Illustration the kinetic energy per unit length under harmonic oscillator asymmetric potential system.
(a)
Kinetic Energy,g=1, T=5K
(b)
Kinetic Energy,g=1, T=5K
(c)
Kinetic Energy,λ=1, T=5K
(d)
Kinetic Energy,λ=1, T=5K
(e)
Kinetic Energy,λ=1, g=1Now, the Helmholtz free energy of the harmonic oscillator asymmetric potential system(FNew) is given by
(34)FNew=-1βln⁡QN(β)
or
(35)QN(β)=e-βFNew.The density matrix of the system in the integration representation will be given by(36)e-βFnew=∫-∞∞‍P^u(x,β)dx.
Substituting (29) into (36), we obtain for the partition function of the particle in harmonic oscillator asymmetric potential system
(37)e-βFnew=∫-∞∞‍γe-(η(2x+ϕ)2-μ)dx=γ2πη,=eλf/4(λ+4g)2sinh⁡(f).
Thus, we have the following results:
(38)Fnew=1βln⁡2sinh⁡f-λf4β(λ+4g).
This is the Helmholtz free energy for a one-dimensional harmonic oscillator asymmetric potential (Fnew), as already derived. When we write the Helmholtz free energy in new variable is Gnew=2Fnew/ℏω,  f=ℏωβ/2 into (38), we obtain that
(39)Gnew=1fln⁡2sinh⁡f-λ4(λ+4g).
One may summarize the numerical results of Helmholtz free energy for harmonic oscillator asymmetric potential by Feynmans approach in Section 5. In Section 5, comparison of Helmholtz free energy for harmonic oscillator asymmetric potential obtain from Feynmans approach and the path-integral method κ(y). This work aims at computing the entropy (S) of the harmonic oscillator asymmetric potential using the Helmholtz free energy in (38). Using the relation
(40)Snew(T)=-∂FnewT∂T=kB∂∂TTln⁡QNβ
the entropy of a system can be write as
(41)SnewT=kBℏω4kBTcosh⁡⁡ℏω/2kBTsinh⁡⁡ℏω/2kBT-kBln⁡(2sinh⁡(ℏω/2kBT))-kBλωℏ8kBT2(λ+4g).
Equation (41) is plotted in Figure 3(a) assuming the value of axes-shift potential g equal -2 and vary the value of asymmetric potential λ. From (41) are plotted in Figure 3 assuming the value of asymmetric potential λ equal 0.05 and vary the value of axes-shift potential g. Figure 3 shows that, as λ increases, the depth of the entropy decreases and its minimum moves farther away the origin.Illustration of the entropy(S) of a single particle under harmonic oscillator asymmetric potential.
(a)
Entropy(Snew(T)), g=-2
(b)
Entropy(Snew(T)), g=-4
(c)
Entropy(Snew(T)), λ=0.05
(d)
Entropy(Snew(T)), λ=0.05


## 4. Evaluate Helmholtz Free Energy via Path-Integral Method (κ(y))

In the Helmholtz free energy the quantity of the thermodynamics of a given harmonic oscillator asymmetric potential system is derived from its the path-integral method:(42)κ(y)=6mkBTπℏ2∫-∞∞‍e-(6mkBT(y-z)2/ℏ2)V(z)dz.
In this case of the harmonic oscillator asymmetric potential is
(43)V(z)=12kz2+τλ(z-φg),
where setting k=mω2, τ=mω/ℏℏω/2, φ=ℏ/mω and substituting (43) into (42), we can write classical κ(y) to simply produce
(44)κ(y)=Φ∫-∞∞‍12kz2+τλz-φge-α(z-y)2dz,
where α=6m/βℏ2, Φ=6m/βπℏ2. The final the path-integral method is
(45)κy=Φπk+2ky2α+4αλτy-φg4α3/2.
In such case it is more useful to write the partition function of κ(y) as
(46)e-βFcl-K=Ω∫-∞∞‍e-βκ(y)dy,
where Ω=m/2πℏU. Substituting (45) into (46) can obtain the partition function of path-integral method (κ(y)) in case of harmonic oscillator asymmetric potential
(47)e-βFcl-K=Ω∫-∞∞‍e-βΦπk+2ky2α+4αλτy-φg/4α3/2dy=Ωe-βπΦk/4α3/2-4gβατλφ/4α3/2×∫-∞∞‍e-βπΦ(2kαy2+4ατλy)/4α3/2dy,e-βFcl-K=2π4ΩβkΦ/α[e-βΦπ/4kα3/2(k2-2αλ2τ2-4gkαλτφ)].
Substituting α=6m/βℏ2, k=mω2, τ=mω/ℏℏω/2, φ=ℏ/mω, Φ=6m/βπℏ2, and Ω=m/2πℏU into (47) leads to
(48)e-βFcl-K=1ℏβω×e-β2ℏ2/24m2ω2(m2ω4-3λ2m2ω3/βℏ-12gλm2ω2/βℏ),-βFcl-K=-ln⁡(2f)-f26-fλ24-gfλ.
Accordingly, the Helmholtz free energy the path-integral method (κ(y)) is given by
(49)Fcl-K(f)=1βln⁡2f+f26-fλ24-gfλ.
After introducing new variable of the Helmholtz free energy path-integral method (κ(y)) as Gcl-K=2Fcl-K/ℏω, f=ℏωβ/2, the Helmholtz free energy of the path-integral method (κ(y)) is then transformed to
(50)Gcl-K(f)=1fln⁡(2f)+f6-λ24-gλ.
This is known as the Helmholtz free energy of the path-integral method (κ(y)). We can calculate numerical of the Helmholtz free energy Feynmans (Gnew(f) approach in Section 3) and the Helmholtz free energy of the path-integral method (κ(y))(Gcl-K(f) in Section 4) for the harmonic oscillator asymmetric potential by mathematica program to give Tables 1 and 2 as shown in Section 5.Table 1
Comparison of the Helmholtz free energy of the harmonic oscillator asymmetric potential obtained from Feynmans approach(Gnew), path-integral method κ(y)(Gcl-K(T)), where λ is the value of asymmetric potential, g is the value of axes-shift potential system, and T is the temperature, setting g=λ=1.
T
(
K
)
G
new
(
T
)
G
c
l
-
K
(
T
)
%Difference
1
−1.167350290774
−1.166666666667
0.058579189
2
−3.981008518269
−3.980922055573
0.002171902
3
−7.813921629128
−7.813895954231
0.000328579
4
−12.319532395568
−12.319521555626
0.000087989
5
−17.327718009705
−17.327712457674
0.000032041
6
−22.737227955451
−22.737224741848
0.000014133
7
−28.480839348833
−28.480837324869
7.10640
×
1
0
-
6
8
−34.510649356211
−34.510648000211
3.92922
×
1
0
-
6
9
−40.790784085205
−40.790783132793
2.33487
×
1
0
-
6
10
−47.293369220882
−47.293368526548
1.46814
×
1
0
-
6
11
−53.996120765666
−53.996120243988
9.66139
×
1
0
-
7
12
−60.880815552301
−60.880815150468
6.60032
×
1
0
-
7
13
−67.932273353648
−67.932273053847
4.65256
×
1
0
-
7
14
−75.137653101331
−75.137652848275
3.36789
×
1
0
-
7
15
−82.485950683258
−82.485950477511
2.49433
×
1
0
-
7Table 2
Comparison of the Helmholtz free energy of the harmonic oscillator asymmetric potential obtained from Feynmans approach(Gnew), path-integral method κ(y)(Gcl-K(T)), where we may vary the value of g,λ,T.
λ
g
T
(
K
)
G
new
(
T
)
G
c
l
-
K
(
T
)
%Difference
−1
1
1
0.832649709225
0.833333333333
0.082068555
−2
2
2
0.268991481731
0.269077944427
0.032138119
−3
3
3
0.186078370872
0.186104045769
0.013796943
−4
4
4
0.930467604432
0.930478444374
0.001164993
−5
5
5
2.672281990295
2.672287542326
0.000207763
−6
6
6
5.512772044548
5.512775258152
0.000058294
−7
7
7
9.519160651167
9.519162675130
0.000021262
−8
8
8
14.739350643789
14.739351999789
9.199869
×
1
0
-
6
−9
9
9
21.209215914795
21.209216867207
4.490557
×
1
0
-
6
−10
10
10
28.956630779118
28.956631473452
2.397841
×
1
0
-
6
−11
11
11
38.003879234334
38.003879756012
1.372696
×
1
0
-
6
−12
12
12
48.369184447699
48.369184849532
8.307624
×
1
0
-
7
−13
13
13
60.067726646352
60.067726962410
5.261694
×
1
0
-
7
−14
14
14
73.112346898668
73.112347151725
3.461208
×
1
0
-
7
−15
15
15
87.514049316743
87.514049522489
2.351017
×
1
0
-
7


## 5. The Numerical Result

In this section, we have evaluated numerical of the Helmholtz free energy Feynmans (Gnew(f)) approach and the Helmholtz free energy of the path-integral method (κ(y))(Gcl-K(f)) for the harmonic oscillator asymmetric potential is (39) and (50). By substituting f=ℏω/2kBT into (39) and (50), we obtain the Helmholtz free energy Feynmans (Gnew(T)) approach as follows:
(51)Gnew(T)=2kBTℏωln⁡2sinh⁡ℏω2kBT-λ4(λ+4g)
and the Helmholtz free energy of the path-integral method (κ(y))(Gcl-K(T))(52)Gcl-KT=2kBTℏωln⁡ℏωkBT+ℏω12kBT-λ24-gλ.
For the beginning of the numerical shooting method, we need to input the parameters λ, g, and T and (39) and (50) into the mathematical program [9]. Next, we obtain the numerical values of Gnew(T) and Gcl-K(T) for the harmonic oscillator asymmetric potential as displayed in Tables 1 and 2.From Table1 we show that magnitude of the Helmholtz free energy decrease from -1.1673 to -82.4859 with increasing the temperature of system. From Table 2 we illustrate magnitude of the Helmholtz free energy increase from 0.8326 to 87.5140 with increasing the temperature of system.Nevertheless,Gcl-K(T) is not as useful as Tables 1 and 2 might indicate. The first, it can not be used when quantum mechanics exchange effect exist. The second, it fails in its present form when the potential V(y) has a very large derivative as in the case of hard-sphere interatomic potential (see Lennard-Jones potential [10]).


## 6. Conclusion

This work is interested in finding the canonical density matrix by Feynman’s technique and kinetic energy per unit length and evaluation of the Helmholtz free energy. The results show that the magnitude of the density matrix(P^u) and the kinetic energy per unit length (τ^L) vary according to the parameters λ, g, and T, respectively. The magnitude of (P^u) and (τ^L) for these cases λ is positive with the graph of (P^u) and (τ^L) moving from right hand side to left hand side. Moreover, we found that as the values of g have increased and the amplitude of (P^u) and (τ^L) has decreased.



---

*Source: 912054-2014-11-02.xml*