Entropy production in quantum Yang-Mills mechanics in semi-classical approximation

We discuss thermalization of isolated quantum systems by using the Husimi-Wehrl entropy evaluated in the semiclassical treatment. The Husimi-Wehrl entropy is the Wehrl entropy obtained by using the Husimi function for the phase space distribution. The time evolution of the Husimi function is given by smearing the Wigner function, whose time evolution is obtained in the semiclassical approximation. We show the efficiency and usefullness of this semiclassical treatment in describing entropy production of a couple of quantum mechanical systems, whose classical counter systems are known to be chaotic. We propose two methods to evaluate the time evolution of the Husimi-Wehrl entropy, the test-particle method and the two-step Monte-Carlo method. We demonstrate the characteristics of the two methods by numerical calculations, and show that the simultaneous application of the two methods ensures the reliability of the results of the Husimi-Wehrl entropy at a given time.


Introduction
Thermalization process or entropy creation of isolated quantum systems is a long-standing issue, but not well understood problem. Relevant systems include the early universe where the transition from a vacuum state to a thermalized state occurs at the end of cosmic inflation, and the QCD matter created in the initial stage of relativistic heavy-ion collisions where thermal matter should be formed in a rather short time. It is known that both systems are well described in semiclassical approximation, and moreover a chaotic behavior of the classical limit may play some role in the entropy production. The present paper is concerned with the entropy production of an isolated quantum system for which the semiclassical approximation is valid and the classical counter part may show a chaotic behavior.
To describe entropy in a pure quantum system, one may of course adopt the von Neumann entropy [1] as quantum mechanical entropy given by where ρ is the density matrix. For a pure state, however, ρ is idempotent, ρ 2 = ρ, implying that the eigen value of ρ is 0 or 1, and the von Neumann entropy is zero. Even if we start from a mixed state, the time evolution described by a unitary operator will never lead to entropy growth. On the other hand, the entropy production in a rarefied gas composed of classical or quantum mechanical particles can be well described by an analog of the H function of Boltzmann given in terms of the distribution functionf (q, p): It is noteworthy that a phase-space description is desirable for making classical-quantum correspondence clear, and even natural when the semiclassical approximation is valid. The standard method for such a description is to use the celebrated Wigner function [2], which is defined as a Wigner transform of the density matrix: The Wigner function f W (q, p) can be regarded as a quasi phasespace distribution function. The use of the Wigner function as the phase space distribution function (f (q, p) = f W (q, p) in Eq. (2)), however, has essential drawbacks: First, the Wigner function is, actually, not a genuine distribution function; f W can be negative, which prevents us to calculate the entropy density according to Eq. (2). Second, the entropy defined by Eq. (2) given in terms of the Wigner function f W does not grow in time, because the Wigner transform only gives an equivalent description of the quantum system in terms of, say, the qor p-representation [3,4,5,6,7]. Some coarse graining of the phase space is needed to describe an entropy production. In a classical chaotic system, two adjacent points in the phase space depart from each other exponentially in time. If available phase space volume is limited, the exponentially diffusing classical trajectories have to be folded in a complicated manner in the phase space. After a certain time starting from a localized phase space cell, a given phase space cell (2π ) D consists of the mixture of trajectories stemming from the initially occupied localized cell and vacant regions not yet visited. Since we cannot distinguish the phase space points in a cell due to the uncertainty principle, it is reasonable to define a phase space distribution as a smeared or coarse-grained function over the phase space cell.
We adopt the Husimi function f H (q, p) [8] as such a coarse-grained distribution function, which is defined as the expectation value of the density matrix with respect to a coherent state |z . It is readily shown that f H (q, p) is semi-positive definite, f H ≥ 0, and a coarse-grained function of the Wigner function, as will be shown in a later section. It is shown [9,10] that the Husimi function faithfully describes the characteristic properties of the underlying classical system, and has been utilized to identify the chaotic remnants in quantum systems [9,10,11,12]. Thus a natural candidate of the quantum mechanical entropy is given by (2) with f (q, p) being substituted by the Husimi function f H (q, p). This entropy was introduced by Wehrl [13] and may be called the Wehrl entropy, although he himself called it the classical entropy and failed in identifying the distribution function f H (q, p) with the Husimi function: Such an identification was made later [14]. We refer to the Wehrl entropy obtained by using the Husimi function as the Husimi-Wehrl (HW) entropy [15], It is worth mentioning that the HW entropy can be a good measure for a quantum entanglement of a system including quantum optical systems [16,17]. For a one-dimensional case, there is a minimum of S HW = 1 [18,19], in contrast to the von Neumann entropy, which takes S vN = 0 in the ground state. It is also shown that the HW entropy takes a value close to the von Neumann entropy at high temperature, and its growth rate coincides with the Kolmogorov-Sinaï entropy for the one-dimensional inverted harmonic oscillator [15]. A direct evaluation of the HW entropy for a quantum system is a kind of challenge even for the system with a few degrees of freedom because it involves a large-dimensional integral over the phase space even apart from the cumbersome calculation of the logarithm with precision. Nevertheless the HW entropy and its time evolution have been calculated for some quantum systems [20,21]. The equation of motion (EOM) of the Husimi function is given in [10], which contains a term of the order , and thus has a more complicated form than that of the Wigner function even in the semiclassical approximation; see below. To solve the complicated EOM of the Husimi function, a test-particle method was proposed by Tsai and Muller [21], where the evolution of the test particles are determined to reproduce some of the moments.
As already mentioned, the semiclassical approximation is suitable to reveal the effect of the chaotic nature of the classical counter part. It is noteworthy that the time evolution of the Wigner function in the semiclassical approximation where the O( 2 ) terms are ignored is readily obtained by solving the classical Hamilton equation; quantum mechanical information such as the uncertainty relation is encoded in the initial Wigner function, provided that it is given as the Wigner transform of the quantum density matrix.
The time evolution of the Husimi function is given by smearing the time-evolved Wigner function obtained in the semiclassical approximation. This is the method we adopt in this article. We shall show its efficiency and usefulness in describing entropy production using a couple of quantum mechanical systems whose respective classical counter systems are known to be chaotic. We propose two methods to evaluate the time evolution of the Husimi-Wehrl entropy. One is an adaptation of the usual test-particle method without recourse to the moments of the distribution function. The other is a sequential application of Monte-Carlo integration, which we call the two-step Monte-Carlo method. We shall demonstrate the characteristics of the two methods by numerical calculations, and show that the simultaneous application of the two methods ensures the reliability of the results of the HW-entropy's time evolution. It should be noted that these two methods are, in principle, applicable to systems with large degrees of freedom such as quantum field theories.
The paper is organized as follows. In Sec. 2, we summarize some basic ingredients of the Wigner and Husimi functions together with the HW entropy. In Sec. 3, we introduce the two numerical methods to evaluate the HW entropy in an efficient way. In Sec. 4, the quantum mechanical models are introduced and numerical results of the Husimi-Wehrl entropy are shown. The final section is devoted to a brief summary and concluding remarks.

Wigner function, Husimi function, and Husimi-Wehrl entropy
In this section, we briefly review quantum mechanical phase space distribution functions, Wigner [2] and Husimi [8] functions, and the phase space expression of the entropy, Husimi-Wehrl entropy [13]. While we introduce Wigner and Husimi functions in one-dimensional quantum mechanics in Subsec. 2.1 and 2.2, extension to multi-dimensional cases is straightforward.

Wigner and Husimi functions
The Wigner function [2] is defined as a Wigner transform of the density matrix While the Wigner function f W (q, p) can be regarded as a quasi phase space distribution function and provides intuitive picture of the phase space dynamics, it is not semi-positive definite and hence we cannot regard f W (q, p) as the phase space probability density. In order to overcome the above drawbacks of the Wigner function, Husimi introduced a Gaussian smeared Wigner function [8], known as the Husimi function, where ∆ is an arbitrary width parameter that gives the smearing manner in the phase space. The Husimi function is defined also as the expectation value of the density matrix with respect to a coherent state |z : for a one-dimensional case with ∆ being an arbitrary constant. Here the coherent state is given by where |0 is the ground state;â|0 = 0. It is readily shown that f H (q, p) is semi-positive definite, f H ≥ 0 by using Eq. (6); f H = | z|ψ | 2 ≥ 0 for a pure state |ψ , and f H = i w i | z|ψ i | 2 ≥ 0 for a mixed state specified by the density matrix The Husimi function f H (q, p) serves as the probability density to observe the phase space variables (q, p) under a minimum wave packet |z , and is now semi-positive definite, f H ≥ 0. Compared with the Wigner function, the Husimi function is smooth and the peak of the Husimi function often appears around the expectation value of the position and momentum [10,22].

Time evolution in semiclassical approximation
The equation of motion (EOM) for the Wigner function f W is obtained from the Wigner transform of the von Neumann equation for the density matrix, ∂ρ/∂t = [H, ρ]/i . By applying the Wigner transform of the operator product, The Wigner transform H W of a Hamiltonian with the form of H = p 2 /2m + U (q) does not change its form. We note that the O( 2 ) term in (8) is proportional to the third derivative of H W or U . Thus the EOM (8) without the O( 2 ) term turns out to be exact for some simple models such as a (an inverted) harmonic oscillator.
The semiclassical EOM for f W is given by retaining the terms up to O( ) in Eq. (8), which reads We remark that the semiclassical EOM is exact for the linear systems mentioned above. Equation (9) asserts that f W is constant along the classical trajectory: Let us see this. Let (q(t;q), p(t;p)) is a solution of the classical EOM, i.e., Hamilton's equation; with an initial condition (q(0) =q, p(0) =p). Then we have for f W (q(t;q), p(t;p), t), Thus we can obtain the semiclassical time evolution of the Wigner function by solving the classical equation of motion. Note that the quantum mechanical effects are taken into account through the distribution of the initial value in the phase space encoded in the Wigner function f W (q, p, 0) constructed from the initial density matrix. It is worth mentioning that the exact analytical solution of the time evolution of f W for some linear systems including a (stable) harmonic oscillator potential [3,15], an inverted (unstable) harmonic oscillator potential [15,23] and an external potential [15] can be obtained. Then even the analytic form of the Husimi function f H (q, p, t) for these systems are readily obtained [15] by the Gaussian smearing of f W (q, p, t), which is easy to perform analytically.
We note here that one may obtain the time evolution of the Husimi function f H (q, p, t) by solving the EOM for f H (q, p, t), which involves terms proportional to , and thus has a more complicated structure than that for f W (q, p, t) even in the semiclassical approximation [10]. If one sticks to solve the EOM for f H directly, some numerical method would be necessary. A test-particle method is adopted as such a numerical method by Tsai and Muller [21], where the time evolution of test particles are determined so as to reproduce some of moments. We remark that there are some ambiguities in such an approach inherent in the moment method.
In this work, we do not adopt this direct method for obtaining the time evolution of the Husimi function f H (q, p, t). We take advantage of the fact that the EOM of the Wigner function f W (q, p, t) in this regime is given simply by solving the classical EOM, and obtain f H (q, p, t) by the Gaussian smearing of thus obtained f W (q, p, t). This strategy should be workable and natural when the semiclassical approximation is meaningful. The remaining task that we have to do for obtaining the Husimi function is just the multi-dimensional integrations over the phase space with the Gaussian kernel for the smearing, which should be feasible by standard methods such as the Monte-Carlo integration.

Husimi-Wehrl entropy
Since the Wigner function f W is merely the Weyl transform of the density matrix, any observable is calculable in terms of f W in principle, and it is also the case with the Husimi function f H . A drawback of the f W is that it can have negative values, and hence is not suitable for the calculation of entropy. As is mentioned in Introduction and the previous subsection, the Husimi function is, in contrast, a semi-positive definite coarse-grained phase space distribution function smeared by a minimum wave packet, and hence a good candidate for the phase space distribution f (q, p) to evaluate the entropy of a quantum system, as the H function of Boltzmann in the classical system, Eq. (2), or equivalently the Husimi-Wehrl entropy given in Eq. (3) [13].
An explicit form of the HW entropy in terms of the Wigner function is given by substituting the D-dimensional extension of Eq. (5) into Eq. (3), One may now recognize some difficulty of the numerical evaluation of the HW entropy: It involves repeated numerical integrations over the multi-dimensional phase space, and in particular one of them appears as an argument of logarithm, which turns out to be quite problematic in the Monte-Carlo integration.

Numerical methods to analyze the semiclassical time evolution of Husimi-Wehrl entropy
Here, two numerical methods are introduced to calculate the time dependence of the HW entropy as given by the Gaussian smearing of the Wigner function obtained in the semiclassical approximation. Both methods are based on an adaptation of the Monte-Carlo integration over the phase-space. We call the two methods the test-particle (TP) and two-step Monte-Carlo (tsMC) methods, respectively. In this section, we deal with the D-dimensional system described by the Hamiltonian H = H(q, p), where q and p denote the D-dimensional vector, respectively, i.e., q = (q 1 , q 2 , . . . , q D ) and p = (p 1 , p 2 , . . . , p D ).

Test-particle method
In the test-particle method [24,25,26,27], the Wigner function is represented as a sum of the delta functions, with the initial function where N TP is the total number of the test particles, and their coordinates are given by (q i (t), p i (t)).
The initial distribution of the test particles (q i (0), p i (0)) (i = 1, 2, . . . , D) is chosen so as to well sample that of f W (q, p, 0): Hence N TP is called the sampling number. The time evolution of the coordinates (q i (t), p i (t)) is determined by the EOM for f W (q, p, t), which is reduced to the canonical equation of motion, in the semiclassical approximation. For the test-particle representation of the Wigner function Eq. (14), the Husimi function is readily expressed as It is noteworthy that the Husimi function here is a smooth function in contrast to the corresponding Wigner function in Eq. (14). Inserting the Wigner function (14) into Eq. (13), the HW entropy in the test-particle method is given as, Now note that the integral over (q, p) i for each i has a support only around the positions of the test particles (q i (t), p i (t)) due to the Gaussian function, and then we can effectively perform the Monte-Carlo integration as follows; By generating a set of random numbers (Q, P ) i with standard deviations of /2∆ and ∆/2, Monte-Carlo sampling point (q, p) i for each i is obtained as (q, p) i = (Q, P ) i + (q i , p i ). Thus we reach the formula to be used in the actual evaluation of the HW entropy in the test-particle method: where the amount of the sample number of (Q, P ) i is denoted by N MC .

Two-step Monte-Carlo method
The second method is a direct Monte-Carlo evaluation of the multi-dimensional integrals. We rewrite Eq. (13) as where (Q k , P k ) and (Q ′ l , P ′ l ) are Gaussian random numbers for the Monte-Carlo (MC) integration to compute Husimi function f H (q, p). For the (q, p)-integration, we generate MC samples (q ′ , p ′ ) at t = 0 according to the initial distribution, and obtain the corresponding phase space sample points (q(q ′ , p ′ , t), p(q ′ , p ′ , t)) at t by solving the canonical equation of motion. Under the semiclassical approximation, f W is constant and the Jacobian is unity along the classical trajectory, J(q(t), p(t)/q ′ (0), p ′ (0)) = 1. Then we can replace the integral over (q, p) in the first line of Eq. (19) with the integral at t = 0 by using the initial distribution and the Liouville theorem as, where (q ′ , p ′ ) are the phase space coordinates at t = 0, and (q(q ′ , p ′ , t), p(q ′ , p ′ , t)) are those at t evolved from (q ′ , p ′ ). The Wigner function at t in the log in Eq. (19) can be obtained by the trace back of the trajectory from t to t = 0 as shown in Eq. (12). Equation (19) contains an MC integral of a function obtained by an MC integral; we first generate (q ′ , p ′ ) at t = 0 according to the distribution f W (q, p, 0) and (Q, P ) as Gaussian random numbers, and then perform the MC integral in the log by generating MC samples (Q ′ , P ′ ). We call this procedure two-step Monte-Carlo (tsMC).
In the following sections, we show the characteristic properties of the two methods and demonstrate numerically how they work using two-dimensional quantum-mechanical systems.

Numerical calculation of Husimi-Wehrl entropy in quantum Yang-Mills model
In this section, we show the numerical results of the HW entropy in "quantum Yang-Mills system" [28], obtained by the two distinct methods, TP and tsMC methods.

Model Hamiltonian and setup of initial condition
The Hamiltonian of the system is given by We have restricted ourselves to the two-dimensional case here. The name, "quantum Yang-Mills (qYM)", is originated from the fact that the spatially uniform Yang-Mills system is reduced to a (0 + 1)-dimensional system, i.e., a quantum mechanical system, and its Hamiltonian is just given by Eq. (21). We adopt the initial condition given by a minimal wave packet centered at (q 1 , q 2 , p 1 , p 2 ) = (0, 0, 10, 10), This initial condition is also adopted in Ref. [21].
In the following, we show numerical results calculated by using the TP and tsMC methods. We show the results in the unit with m = 1 and = 1, and take ∆ = 1 for the wave packet width. In the case of ∆ = 1, the smearing Gaussian is not symmetric in p and q directions. But the results do not change qualitatively. We have confirmed that the results with ∆ = 0.1 and 10 are qualitatively the same as those with ∆ = 1. Time dependence of the HW entropy by using TP method in qYM, with N TP = 100, 1000, 5000 and 15000, and N MC = 500. The arrow shows how the calculated HW entropy changes as N TP increases.

Numerical results with TP method
First, we show the numerical results of the HW entropy in the qYM system calculated in the TP method using Eq. (18). Figure 1 shows the time evolution of the HW entropy calculated in the TP method with the following test-particle numbers, N TP = 100, 1000, 5000 and 15000. The MC sample number is taken to be N MC = 500. The statistical errors are estimated for N MC samples from a standard deviation. We note that the calculated HW entropy at each t tends to increase along with increasing N TP , which is an artifact due to the small number of the test particles N TP and discussed later. Apart from tiny fluctuations, all the calculation show that the HW entropy first increases in time with a small oscillatory behavior being accompanied; its local maxima are seen around t ≃ 0.5 and 1.7. We note that a similar behavior is also seen in Ref. [21].
Entropy evaluated by the TP method has a (unphysical) maximum depending on N TP , which causes apparent saturation at large t in Fig. 1. In fact, when the system is chaotic and the phase space volume is very large, all the test particles will be so separated from each other in the phase space at later time that only the i = j terms in Eq. (18) will remain. In this limiting case, the HW entropy as given in (18) is evaluated as follows; which gives the inevitable upper limit of S (TP) HW . In Appendix A, we examine the HW entropy of an inverted harmonic oscillator, for which S HW can be calculated analytically and is found to increase permanently. At later times, S HW is underestimated with small N TP values because of the upper limit discussed above. By comparison, S HW at early times is calculated precisely in the TP method, as long as N TP is large enough for S HW to converge.
From the above argument, S HW (t) would be obtained reliably as an extrapolated value in the limit of N TP → ∞. The extrapolation should be made in the N TP range, where the limiting value is larger than the HW entropy to be obtained. The limiting values are S for N TP = 100, 1000, 5000 and 15000, respectively. The large-t values found in Fig. 1 are close to these limiting values for smaller N TP , i.e., N TP = 100 and 1000. Thus we see that the saturation behavior seen for smaller values of N TP may be an artifact of the TP method. In contrast , the larget values for N TP = 5000 and 15000 in Fig. 1 are well below the limiting values (9.1 and 10.2), found free from the above mentioned artifact, and can be used to obtain the extrapolated value at N TP → ∞, as discussed later in Subsec. 4.4. Thus we conclude that the entropy production of the "quantum Yang-Mills" system can be well described with the use of HW entropy as calculated with the TP method with sufficiently large number of the test particles.

Numerical results with tsMC method
Next, we show the numerical results of the HW entropy in qYM in the tsMC method using the formula Eq. (19). Figure 2 shows the time evolution of the HW entropy calculated in the tsMC method with the sample numbers N in = 1200, 2400, 4800 and 12000. N out is taken to be the same as N in . The errors attached to S HW in the present figure is estimated only for the Monte-Carlo integrals outside of log in Eq. (19), and those from the integral inside the log is not taken into account, which causes an additional systematic error.
We see that the larger the value of N in , the smaller the HW entropy, which is an opposite dependence on the sample number to that in the TP method. Nevertheless the gross behavior in the time evolution of the HW entropy is quite similar in the two methods apart from the tiny fluctuations; After showing an oscillatory behavior in a first short period, it increases in a monotonous way and its growth rate decreases gradually. More quantitative comparison of the two methods will be presented in the next subsection. Figure 3 shows the HW entropy at t = 10 as a function of N TP (N in ) in the TP (tsMC) method. We fit a linear function f (t) = at + b to the calculated S HW (t) data in the range 10 − ∆t ≤ t ≤ 10 + ∆t (∆t = 1), and adopt f (t = 10) as the HW entropy value at t = 10. This procedure provides a smoother curve and reduces the errors coming from fluctuations compared to directly using the raw data. The HW entropy in the TP method becomes larger with increasing N TP as already mentioned; At t = 10, S HW ≃ 5.1 for N TP = 100 and S HW ≃ 8.7 for N TP = 15000. We also show the fit results to the data for larger samples, say N T P ≥ 5000, with a fit function,

Comparison of the two methods
The extrapolated value to N TP → ∞ is 9.19 ± 0.
With increasing N in , the HW entropy calculated in the tsMC method decreases, which is an opposite behavior to that in the TP method as noted before. At t = 10, S HW ≃ 13.2 for N in = 1200 and S HW ≃ 9.5 for N in = 12000. We also show the fit results to the data. We adopt Eq. (24) for the fit function. From the fit results, the HW entropy in the tsMC method is found to be S (tsMC) HW (t = 10) = 9.01 ± 0.21 (stat.) ± 0.06 (syst.) , where the central value and the statistical error are obtained from the fit using Eq. (24), and the systematic error is evaluated from the fits using several fit functions as done in the TP method.

Discussions
The time evolution of the HW entropies obtained in the TP and tsMC methods shows a similar behavior with each other: The HW entropy increases with an oscillatory behavior in the early stage, then shows a monotonous increase with a decreasing rate. The HW entropy at each t in the TP method increases along with N TP , while it decreases with increasing N in in the tsMC method. Thus we can guess that the real value of the HW entropy lies between the results in the TP and  (t = 10) = 9.01 ± 0.21 ± 0.06 at N in → ∞ in the TP and tsMC methods respectively, are consistent with each other within the error. These results are also in agreement with that in Ref. [21].
These two methods, TP and tsMC methods, give consistent results after N → ∞ extrapolation. On the other hand, with finite number of N TP and N in , they could give seemingly inconsistent results depending on the dynamics. We here have a deeper look at this issue. In the tsMC method, the entropy seems to keep increasing even for the later time, in contrast to the results in the TP method with finite N TP and in Ref. [21]. The discrepancy may come from the special shape of the potential: there are two flat directions in the potential for the qYM system, although the width of them tends to shrink at large distances. Then, the classical trajectory can keep growing along the flat direction, which would cause an unlimited spreading of the Husimi function and a permanent increase of the HW entropy calculated in the semiclassical approximation. (In the case of the TP methods, there exists limiting value of the HW entropy depending on N TP , which gives rise to the apparent saturation of S at large t. ) By comparison, it is shown that the exact energy spectra of the qYM are all discrete ones, because of the shrinking width leading to an increase of the kinetic energy due to the uncertainty relation, although the volume of {(p, q)|H(p, q) ≤ E} is infinite [29]. Note that the discrete spectra implies that the wave functions of the energy eigen states are all bound. Thus the corresponding Husimi function would not have a support at the infinite distance due to the quantum effect, and the HW entropy may not show the ever increasing behavior but have a saturated value. This plausible conjecture can only be confirmed by a full quantum calculation beyond the semiclassical approximation. Such a calculation is beyond the scope of the present work and will be left as a future work. Instead, we shall take another model, which is a modified version of the qYM one free from flat directions in its potential.

Modified quantum Yang-Mills model
Let us consider the model in which quartic potential terms are added to the qYM Hamiltonian; We call the system "modified quantum Yang-Mills (mqYM)". The system is studied in Ref. [11,12] with g 2 < 0 in the context of chaos. It is apparent that there is no flat direction in the potential due to the quartic terms. We take g 2 = 1 and ǫ = 0.1 in the Hamiltonian, Eq. (27). The mqYM system is found to be integrable with ǫ/g 2 = 1, 1/3 and ∞ [12,30]. Our choice of ǫ/g 2 = 0.1 is well apart from the integrable region. Since ǫ is not very large, the HW entropy shows a similar behavior to that in qYM at early times, as shown later.
In this section, we shall calculate the HW entropy of the mqYM system in the TP and tsMC methods . The analyses are carried out in a similar way to those for the qYM system.
In Figs. 4 and 5, we show the time evolution of the HW entropy in mqYM calculated using the TP (N TP = 500, 1000, 5000 and 15000 with N MC = 500) and tsMC (N in = 600, 1200, 2400 and 12000) methods, respectively. N out is taken to be the same as N in for tsMC.
The distribution function in Eq. (22) is used as the initial condition, and the statistical errors are estimated for N MC (N in ) samples from a standard deviation in the TP (tsMC) method, as in the qYM cases.
Both of the calculated results show that the HW entropy first increases with an oscillatory behavior and tends to saturate at later times, t 6. The later-time S HW values depend on the sample number, N TP and N in ; With increasing N TP (N in ), the HW entropy increases (decreases) in the TP (tsMC) method. These are the features also found in qYM. By comparison, it should be noted that there seems to be saturation of S HW both in the TP and tsMC methods in mqYM, in contrast to qYM. This may be originated from the finite phase space volume where the Husimi function has a support.
In Fig. 6, we show the HW entropy at t = 10 as a function of N TP or N in . We fit a linear function to calculated S HW (t) results in the range 9 < t < 11, and adopt f (t = 10) as the HW entropy value at t = 10. In the TP method, S in the TP and tsMC methods, respectively. The central values and the statistical errors are obtained from the fit using Eq. (24), and the systematic error is evaluated from the fits using several fit functions. These two values are consistent with each other within the error. The observation shows that the two methods, tsMC and TP, are especially effective for such a potential which bounds Husimi function in finite region. Thus, we are confident of the validity of the two methods in the mqYM system.

Summary
We have discussed entropy creation in isolated quantum systems by using the Husimi-Wehrl entropy evaluated in a semiclassical treatment. The semiclassical treatment is known to be useful in some of the systems such as the inflation in early universe and the early stage of relativistic heavy ion collisions. These systems are expected to bear instabilities and/or chaoticities in their classical counter systems, then the smearing of the phase space distribution by the minimal wave packet causes the entropy production in terms of the Wehrl entropy or the H function of Boltzmann even in isolated quantum systems. This is nothing but the Husimi-Wehrl entropy, the Wehrl entropy obtained by using the Gaussian smeared Wigner function (Husimi function) for the phase space distribution.
The semiclassical time evolution of the Husimi function is given by solving a classical equation of motion and smearing with a Gaussian packet. Combining this semiclassical treatment with the Monte-Carlo numerical integral technique, we have developed two methods, the test-particle (TP) method and the two-step Monte Carlo (tsMC) method. We have applied these two methods to quantum mechanical systems in two dimensions, the quantum Yang-Mills (qYM) and the modified quantum Yang-Mills (mqYM) systems. The classical counter systems of these are known to be chaotic. We have demonstrated that the Husimi-Wehrl entropy obtained in the TP (tsMC) method approaches the converged value from below (from above) with an increasing sample number, then we can guess the true value of HW entropy. We have further found that the results of the TP and tsMC methods in the infinite sampling number limit are consistent within the error. Therefore, the simultaneous application of the two methods ensures the reliability of the results of the Husimi-Wehrl entropy at a given time.
The extension of our methods to a multidimensional system is straightforward. We expect that these methods are useful in systems with many degrees of freedom such as the quantum field theory. These methods are, in principle, applicable to higher-dimensional problems, and we have confirmed that they actually work in three and four dimensional systems. In higher dimensions, we need much more Monte-Carlo samples to obtain statistically reliable results, and it would be necessary to make some approximations for practical purposes. Work in this direction is in progress.  Figure A1 shows the time evolution of the HW entropy of IHO calculated in the TP method with N TP = 50 − 800. We find that the TP method can well describe the time evolution of the HW entropy at early times, and that numerical results show saturated behavior in later times. Since there exists a limiting value of S HW in the TP method as discussed in Subsec. 4.2, we need to take a large number of N TP to describe a large amount of entropy production. It should be noted that numerical results converge in the limit of N TP → ∞, and the converged result well describe the analytic result.  Figure A2 shows the time evolution of the HW entropy of IHO in the tsMC method with N in = N out = 100 and 1000. We find that numerical results are consistent with the analytic solution at early times t ≤ 3, but that the numerical results tend to overestimate the analytic results and numerical errors become very large at later times. The large error would come from the poor overlap between the Wigner function and the coarse-graining Gaussian function at later time, which makes