text:
# 1 Introduction

## 1 Introduction

The study of quantum phase transitions represents a wide area of research in condensed matter physics . These transitions occur at $`T=0`$ when a parameter of the hamiltonian is varied across a critical value and are driven by quantum fluctuations. Among them, those that are due to the effect of disorder are more difficult to analyze both analytically and numerically and have been the focus of much attention recently.

In this paper we consider the spin $`1/2`$ XXZ chain in the presence of a non-symmetry breaking random magnetic field. This system is particularly interesting since a renormalization group analysis  suggests the existence of a critical region in the space of parameters where the disorder is not relevant. This idea is also supported by numerical calculations . In  the system was studied using an exact diagonalization method up to $`16`$ sites. Although this size is enough to study the system far from the critical region, near the transition, as it is explained there, the calculations were less reliable due to the small size of the chains considered. Larger chains can be handled using a density matrix renormalization group (DMRG) algorithm  at the expense of truncating the Hilbert space of states. This method was used in  to obtain further information about the critical region by computing the phase sensitivity, which is the difference in energy between periodic and antiperiodic boundary conditions.

In this paper we use the DMRG method to compute the correlation length and the superfluid density which give direct information about the existence of a quasi-long-range order region since that region is defined as the one where the correlation length diverges and the superfluid density is non-zero. In particular, to compute the superfluid density we used a superblock configuration that allowed us to introduce an arbitrary twist at the boundary.

Since we perform DMRG calculations, our results contain two approximations. The first is that we consider a finite chain and the second that we truncate the Hilbert space (keeping a finite number of density matrix eigenstates). These approximations are related since with a truncated Hilbert space the correlation length is always finite and, in fact, the value corresponding to an infinite chain (what we are going to call the ’truncated system correlation length’) becomes smaller as we reduce the Hilbert space. To minimize the finite size effects in the calculation of the spin correlations, it is necessary to consider chains larger than the truncated system correlation length. However, in the critical region this correlation length increases rapidly as we keep more DMRG states, therefore, reducing the truncation error by keeping a large number of states soon becomes prohibitive since we need to consider very long chains 111For a more quantitative statement see Fig. 16. Thus, most of our calculations of the spin correlations are done keeping up to 32 states and give only an estimation to the shape of the critical region. As we will discuss later, the main source of error in this estimation is the truncation of the Hilbert space. To further support these qualitative results, for the special case of anisotropy $`\mathrm{\Delta }=0.75`$ and selected values of the random magnetic field, we keep up to 96 states and extrapolate the computed correlation length to infinite number of sites and infinite number of DMRG states. The extrapolations allow us to conclude the existence of a quasi-long-range order region and put bounds on the maximum critical field.

In the calculations of the superfluid density we also keep up to 32 states in most of the cases and estimate the shape of the critical region using a criterion in which, again, the main approximation is the truncation of the Hilbert space. As we did for the correlations, for the special case of $`\mathrm{\Delta }=0.75`$ and selected values of the random magnetic field we keep up to 128 states and extrapolate the superfluid density to infinite length and infinite number of DMRG states. We find an agreement between the correlation length and the superfluid density results.

The remainder of this paper is organized as follows. In section 2 we summarize the known results for this system and in section 3 we describe the numerical method used. We present and discuss the results in sections 4 and 5 and finally, in section 6 we give our conclusions.

## 2 The spin $`1/2`$ $`XXZ`$ chain in a random $`z`$ magnetic field

In this section we summarize previous results obtained for this system by C. Doty and D. Fisher .

In the absence of disorder the spin $`1/2`$ XXZ chain is defined by the Hamiltonian

$$H_0=\underset{i=1}{\overset{L}{}}(S_i^xS_{i+1}^x+S_i^yS_{i+1}^y+\mathrm{\Delta }S_i^zS_{i+1}^z).$$
(1)
When the anisotropy parameter, $`\mathrm{\Delta }`$, satisfies $`1<\mathrm{\Delta }1`$, the system is in a non-magnetic gapless phase where the spin-spin correlation functions in the ground state decay as a power law (quasi-long-range order (quasi-LRO) phase) . On the other hand, the spectrum has a gap for $`\mathrm{\Delta }>1`$ ($`\mathrm{\Delta }1`$) and possesses a long-range Ising like antiferromagnetic (ferromagnetic) order.

In this work we study the effects produced by the addition of the following term:

$$H_1=\underset{i=1}{\overset{L}{}}h_i^zS_i^z,$$
(2)
where $`h_i^z`$ is a random magnetic field with an independent and uniform distribution in each site ($`h_i^z=0`$, $`h_i^zh_j^z=\delta _{ij}h_0^2/3`$). Note that the addition of this term does not break the symmetry in the $`xy`$ plane.

The main conclusion of  is that for $`1<\mathrm{\Delta }<1/2`$ the quasi-long-range order that the system possesses at $`h_z^i=0`$ persists for finite values of the random magnetic field.

To reach this conclusion, the long-distance behavior of the system was analyzed by studying the renormalization group flow. In the region $`|\mathrm{\Delta }|<1`$ one can perform a Jordan-Wigner transformation to spinless fermionic variables. The component $`S^z`$ is converted into the fermion occupation number and $`S^\pm `$ are creation and annihilation operators:

$`S_n^z`$ $`=`$ $`\mathrm{\Psi }_n^{}\mathrm{\Psi }_n{\displaystyle \frac{1}{2}},`$ (3)
$`S_n^{}`$ $`=`$ $`\mathrm{exp}\left(i\pi {\displaystyle \underset{m<n}{}}(S_m^z+{\displaystyle \frac{1}{2}})\right)\mathrm{\Psi }_n.`$ (4)
For the pure XXZ chain when $`\mathrm{\Delta }=0`$, the ground state is a half-filled band with Fermi momenta $`k_F=\pm \pi /2a`$ where $`a`$ is the lattice spacing. Low energy excitations around the Fermi surface can be described in terms of bosonic fields using a standard procedure :

$$\mathrm{\Psi }_n=\sqrt{\frac{a}{L}}\left[e^{ik_Fxi\mathrm{\Phi }_R(x)}+e^{ik_Fx+i\mathrm{\Phi }_L(x)}\right].$$
(5)
Here $`L`$ is the length of the system. The bosonic fields $`\mathrm{\Phi }_{L,R}`$ describe left and right moving excitations.

The random perturbation (2), can also be expressed in terms of the bosonic fields and, at low energies, the continuum limit can be taken which results in an action :

$$S=𝑑x𝑑\tau \left(\frac{1}{2}\kappa (_\tau \mathrm{\Phi })^2+\frac{1}{2}\kappa (_x\mathrm{\Phi })^2+g\mathrm{cos}(2\mathrm{\Phi })+\eta (x)_x\mathrm{\Phi }+\rho (x)e^{i\mathrm{\Phi }}+\rho ^{}(x)e^{i\mathrm{\Phi }}\right),$$
(6)
for the field $`\mathrm{\Phi }(x)=\mathrm{\Phi }_R+\mathrm{\Phi }_L`$. Besides, $`g`$ is a function of the anisotropy $`\mathrm{\Delta }`$, and $`\eta (x)`$ and $`\rho (x)`$ are composed of the Fourier components of the random field near zero and $`2k_F`$, respectively. Also, $`\kappa `$ is the spin stiffness of the ordered system given by the Bethe-ansatz solution :

$$\kappa =\frac{1}{2\pi }\left(1\frac{1}{\pi }\mathrm{cos}^1\mathrm{\Delta }\right).$$
(7)
Afterwards, the replica formalism can be used to compute disorder averages of ground-state expectation values giving the effective action:

$`\overline{S_{eff}}`$ $`=`$ $`{\displaystyle \underset{\alpha =1}{\overset{n}{}}}{\displaystyle 𝑑x𝑑\tau \left[\frac{1}{2}\kappa (_\tau \mathrm{\Phi }_\alpha )^2+\frac{1}{2}\kappa (_x\mathrm{\Phi }_\alpha )^2\right]}`$ (8)
$`{\displaystyle \underset{\alpha ,\beta }{}}D{\displaystyle 𝑑x𝑑\tau 𝑑\tau ^{}\mathrm{cos}[\mathrm{\Phi }_\alpha (x,\tau )\mathrm{\Phi }_\beta (x,\tau ^{})]},`$
where $`D=<(h_i^z)^2>=h_0^2/3`$. Thus, the disorder induces interactions between pairs of replicas. The interaction comes from averaging over $`\rho `$, that is the term due to the staggered component of the magnetic field. The term due to $`\eta `$ is irrelevant at large distances and is not included in (8).

Using the fact that the dimension of the operator $`\mathrm{cos}(\mathrm{\Phi })`$ is $`1/4\pi \kappa `$, the RG eigenvalue $`\lambda `$, which determines the rescaling of $`D`$ is given by :

$$\lambda =1+2z2\zeta ,$$
(9)
where the dynamical exponent $`z=1`$ and $`\zeta =1/(4\pi \kappa )`$.

The sign of $`\lambda `$ determines the relevance of the randomness $`D`$. In the interval $`1/2<\mathrm{\Delta }1`$, $`\lambda `$ is positive, indicating that the disorder is relevant and that the quasi-LRO is unstable under any small amount of randomness. On the other hand, $`\lambda `$ is negative in the region $`1<\mathrm{\Delta }<1/2`$, so there the disorder is irrelevant and the quasi-LRO is expected to persist for small values of $`D`$. Therefore, in the latter case, a quantum phase transition between a quasi-LRO phase and a localized disordered phase has to occur for finite $`D`$ . By means of a RG treatment for weak disorder  the RG equations governing the critical properties can be obtained . The flow lines in the $`\kappa D`$ plane, near the point $`\kappa =1/6\pi `$ are parabolas, $`(\kappa 1/6\pi )^2D=0`$ being the critical one. The phase transition is of the Kosterlitz-Thouless type  and the correlation length corresponding to the spin correlation functions diverges like:

$$\xi \mathrm{exp}(A/\sqrt{\mathrm{\Delta }\mathrm{\Delta }_c}),$$
(10)
when the quasi-LRO phase is approached along a path in the phase diagram at fixed disorder. In the region where the pure system is an Ising like ferromagnet (antiferromagnet) the Imry-Ma argument  can be used to conclude that an arbitrarily weak random $`z`$ field will destroy the long-range order.

The phase diagram $`h_0\mathrm{\Delta }`$ proposed by Doty and Fisher is shown in Fig. 1. Note that the top of the region where the quasi-LRO persists i.e the maximum critical field is not established since the perturbative renormalization group result is valid only near $`\mathrm{\Delta }=0.5`$.

A numerical work analyzing this phase diagram was performed by Runge and Zimanyi . They use an exact diagonalization procedure up to $`16`$ sites and compute, among other things, the superfluid density $`\rho _s`$ which measures the sensitivity of the system to changes in the boundary conditions:

$$\rho _s=\frac{1}{L}\frac{^2E}{(\theta )^2}.$$
(11)
Here, $`\theta `$ is a twist introduced at the boundary by adding to the Hamiltonian (1) the term

$$\frac{1}{2}\left(S_1^+S_L^{}e^{i\theta }+S_1^{}S_L^+e^{i\theta }\right)+\mathrm{\Delta }S_1^zS_L^z.$$
(12)
Doing this, the chain is closed in such a way that a boson acquires a phase $`e^{i\theta }`$ ($`e^{i\theta }`$) if it is hopping to the right (left) through the boundary.

This superfluid density is the order parameter of the transition ($`\rho _s>0`$ in the quasi-LRO phase and $`\rho _s=0`$ in the localized phase).

In , the authors find a qualitative agreement with the RG predictions . However they conclude that it is necessary to study larger chains to obtain more accurate results near the critical region.

In this paper we compute the superfluid density and the correlation length using the DMRG approach which allows us to study longer chains. In that way we expect to obtain further evidence and complementary results.

## 3 Numerical method

We compute the properties of the system using an infinite-size DMRG algorithm . The disorder is included by doing the calculations for a single configuration of the random field and then, averaging over many realizations. Since the ground-state is non-magnetic in the quasi-LRO regime for the pure system and is expected to be also non-magnetic in the whole range of $`\mathrm{\Delta }`$ when the disorder is turned on, we perform our calculations constrained to the Hilbert subspace where the $`z`$-component of the total spin is zero. We have checked that this approximation is justified after sufficient averaging is performed.

We calculate the spin correlations

$$C_j(r)=S_0^xS_r^x+S_0^yS_r^y,$$
(13)
using open boundary conditions and the same superblock configuration $`BB`$ as in  (where $`B`$ represents a block and $``$ a single site). In (13), $`\mathrm{}`$ denotes quantum expectation value in the ground state and $`j`$ is a given realization of the random field. To avoid the boundary effects, we first grow the system until the length of each block is around 30 sites. After that we start to compute the correlations taking the spin operator $`S_0^+`$ corresponding to the left site added at that DMRG step. In the next iteration this operator is rotated and truncated using the eigenvectors of the density matrix and so it will correspond to the right-most site of the left block. Iterating this procedure for the same operator we obtain after $`r`$ DMRG steps an operator $`S_r^+`$ associated to a site which is at a distance $`r`$ from the center. Simultaneously, the correlation $`C_j(r)`$ is computed at each of those DMRG steps as $`S_0^{}S_r^++S_0^+S_r^{}`$. Using this standard procedure it is only necessary to store one spin matrix ($`S_r^+`$). Afterwards, we average the logarithm of the correlation functions $`C_j(r)`$ over disorder to obtain what is called the typical correlation . Usually, more than 150 configurations of the random field are required.

At large distances $`C_j(r)`$ should be exponential outside the critical region where the disorder is relevant (finite correlation length) and power law inside the critical region, where it is irrelevant (infinite correlation length). However, due to the truncation of the DMRG and the finite size effects, we also find an exponential behavior inside the critical region. This fact makes possible to define a typical correlation length from the slopes of the curves $`\mathrm{ln}C_j(r)`$ vs $`r`$ and analyze its behavior in the phase diagram $`\mathrm{\Delta }h_0`$. The correlation length computed in this way is of course approximated and consequently gives qualitative information about the phase diagram.

If chains of various length are considered and, if for each of them, various number of states are kept, an extrapolation to infinite length and infinite number of DMRG states can be performed to decrease the error. This is discussed in section 5 where we do the extrapolations for a few special points in the phase diagram.

To compute the superfluid density which is the order parameter of the transition between the quasi-LRO and the localized states, we need to consider a closed system with twisted boundary conditions. Closed systems can also be treated with the density matrix algorithm as shown in . In that work two superblock configurations were analyzed: $`BB`$ and $`BB`$. It was found that the first one is better since it has a smaller truncation error (because the two blocks are not connected directly), although larger than the error obtained with open boundary conditions. In our case, we found the configuration $`BB`$ more convenient to introduce twisted boundary conditions since the twist can be introduced between the two single sites that we add at each step of the procedure (in that case the operators $`S^+`$ and $`S^{}`$ which appear in (12) are not truncated). On the other hand since we are also interested in minimizing the truncation error, for the actual calculation we chose the configuration in Fig. 2 which combines the advantages of both. Note that the left and right blocks have different number of sites although the same number of states are kept on both of them. This asymmetry implies that we have to compute two density matrices, one for the left block and one for the right one. However this is necessary anyway since the reflection symmetry is already broken by the random magnetic field.

To approximate eq. (11) numerically, given a configuration of the random magnetic field, we perform the DMRG procedure twice (with and without twist) and compute for each length $`L`$ the finite difference

$$\rho _s(L)\frac{1}{L}\frac{E(L\theta )E(0)}{(\theta )^2},$$
(14)
where we used that the first derivative $`E/\theta `$ vanishes at $`\theta =0`$ due to the symmetry $`E(\theta )=E(\theta )`$. Afterwards, we average this quantity over many configurations of the field. The h larger the strength of the field the more configurations needed (in most of the cases more than 200 configurations were considered).

To test the convergence of the DMRG method we computed $`\rho _s`$ keeping different number of eigenstates of the density matrix. It was found that the DMRG procedure converges if we keep $`\theta `$ constant i.e. the twist $`\theta `$ increases with the length of the chain. This is shown in Fig. 3, where $`\rho _s`$ vs. $`L`$ for $`\mathrm{\Delta }=0.75`$ is plotted for a random field $`h_0=0.1`$ keeping up to 128 eigenstates of the density matrix. Since we need $`L\theta `$ to be small, in the calculations we took $`\theta =0.01`$ and grew the system up to $`50`$ sites.

We also checked (see Fig. 4) that without disorder our numerical results for $`\rho _s`$ agree with the exact solution

$$\rho _s(\mathrm{\Delta },h_0=0)=\frac{\pi \mathrm{sin}\mu }{2\mu (\pi \mu )},$$
(15)
found in . Here, $`\mu =\mathrm{cos}^1\mathrm{\Delta }`$.

Finally, since we are truncating the Hilbert space it is important to estimate how good this approximation is. A standard measure of that is the truncation error defined as the sum of the eigenvalues of the density matrix that are discarded. Keeping up to $`32`$ eigenstates of the density matrix we obtained a truncation error less than $`10^7`$ for open boundary conditions (spin correlations) and less than $`10^5`$ for twisted boundary conditions (superfluid density).

## 4 Results

### 4.1 $`S_i^xS_j^x+S_i^yS_j^y`$ correlation

Before summarizing the correlation functions results, let us discuss what should be the behavior of the correlation length in a phase diagram like the one proposed in Fig. 1. As we said before, in the localized phase it is possible to define a correlation length, i.e. the correlations functions display an exponential behavior at large distances. Instead, in the phase with quasi-LRO this is not possible because, since the region is critical, the correlation length is infinite. This means that, if we were able to solve the problem exactly, we should observe in a 3-dimensional plot of correlation length vs ($`\mathrm{\Delta }`$,$`h_0`$) that this quantity diverges when we approach the line that divides the two regimes. This is not the case when we consider a finite chain and moreover when we truncate the Hilbert space. What we should obtain instead is that the correlation length has a peak located in a region that coincides approximately with the critical region. This peak should become sharper, higher and closer to the transition as we consider larger chains and keep more DMRG states.

Having this in mind, we computed the spin-spin correlations

$$\overline{C_t(L)}=\overline{\mathrm{ln}|S_0^xS_L^x+S_0^yS_L^y|},$$
(16)
transverse to the random field for several points of the phase diagram $`\mathrm{\Delta }h_0`$. Here, the bar indicates average over disorder.

In Fig. 5, $`\overline{C_t(L)}`$ is plotted for $`h_0=0.3`$ and different values of $`\mathrm{\Delta }`$. The linear behavior at large distances means the correlations decay exponentially.

From the slope of these lines we estimate the typical correlation length $`\xi `$ and analyzing its behavior as a function of $`\mathrm{\Delta }`$ and $`h_0`$ we obtain the 3-dimensional plot shown in Fig. 6 for $`16`$ and $`32`$ DMRG states. There, we can see that the correlation length not only increases when we approach the axis $`h_0=0`$ (pure system) keeping $`\mathrm{\Delta }`$ constant but also in a region with finite disorder centralized in a negative value of $`\mathrm{\Delta }`$ between $`1`$ and $`0.5`$. This is an indication of the presence of a critical region for finite disorder where the quasi-LRO found for $`h_0=0`$ persists.

To obtain an approximation to the shape of the critical region, we propose to analyze the contour plot of the typical correlation length. Following the previous discussion, in the limit of infinite length and infinite number of DMRG states, the contour lines will concentrate on the border of the critical region. If, for a given truncation, we use the input that the critical region starts at $`\mathrm{\Delta }=0.5`$, the choice of a contour line such that $`h_00`$ at $`\mathrm{\Delta }=0.5`$ will give us qualitative information about the shape of the critical region. This choice will be justified in the next section when we discuss the superfluid density results. Fig. 7 shows the contour plot of the typical correlation length for $`16`$ and $`32`$ DMRG states. From the chosen contour lines (bold lines in the figure) we estimate $`h_{\mathrm{crit}.}(\mathrm{\Delta }=0.75)0.17`$ for $`16`$ states kept and $`h_{\mathrm{crit}.}(\mathrm{\Delta }=0.75)0.20`$ keeping $`32`$ states. We want to remark here that in Fig. 7, the contour lines used to estimate the boundary of the critical region correspond to values of $`\xi `$ smaller than 100 (see Fig 6) while the chains considered were larger than 200 sites. Thus, the finite size effects are small and the main approximation is in the truncation of the Hilbert space.

### 4.2 Superfluid density

We computed the superfluid density using (14) and performed afterwards the average over disorder. In Fig. 8, we show $`\rho _s`$ vs $`L`$ for different values of the strength of the random field $`h_0`$ and $`\mathrm{\Delta }=0.75`$. As expected, for a fixed length of the chain, $`\rho _s`$ decreases with $`h_0`$. Note that we obtain non zero values of $`\rho _s`$ because we study finite systems. For the infinite system the superfluid density is different from zero only for ($`h_0=0`$,$`|\mathrm{\Delta }|<1`$) and inside the critical region, that is, where the correlation length diverges. For that reason $`\rho _s`$ can be used to study the phase diagram of Fig. 1.

First we are going to check that the quasi-LRO region extends up to $`\mathrm{\Delta }=0.5`$. This point is characterized by the fact that there, the correlation length critical exponent diverges. In the region $`1/2<\mathrm{\Delta }<1`$ where the disorder is relevant, the dependence of the correlation length with the disorder $`D=h_0^2/3`$ near the critical line corresponding to $`h_0=0`$, should be $`\xi D^{\varphi _s}`$ and according to the perturbative renormalization group results , the critical exponent $`\varphi _s`$ diverges at $`\mathrm{\Delta }=0.5`$ like

$$\varphi _s=\frac{1}{3(1\frac{1}{\pi }\mathrm{arccos}\mathrm{\Delta })^1}.$$
(17)
To extract this exponent from our numerical data, we analyzed the derivative $`\rho _s/D`$ at $`D=0`$. Assuming that $`\rho _s=\rho _s(L/\xi )`$ for $`L1`$ , the behavior of $`\mathrm{ln}(\rho _s/D|_{D=0})`$ with $`\mathrm{ln}L`$ should be linear with the slope corresponding to $`\varphi _s^1`$. In the inset of Fig. 9 we show this linear behavior for different values of $`\mathrm{\Delta }`$. The critical exponent extracted from our data and the analytical prediction (eq. 17) are shown in Fig. 9 where a good agreement can be observed. This result is in line with the existence of a critical region for $`\mathrm{\Delta }<0.5`$ where the disorder is irrelevant and justify our choice made in the $`\xi `$ contour plots.

For a given truncation, the border of this region can be estimated studying the dependence of $`\rho _s`$ with the strength of the random field for different lengths of the chain, $`L`$. Figure 10 is a plot of $`\rho _s`$ vs $`h_0`$ for $`\mathrm{\Delta }=0.75`$ and different lengths. The exact solution (infinite length and infinite DMRG states) should show a discontinuity in a transition of the Kosterlitz-Thouless type ($`\rho _s(h_0<h_{\mathrm{crit}.})>0`$ and $`\rho _s(h_0>h_{\mathrm{crit}.})=0`$).

Using this fact, an estimation of the critical field can be obtained analyzing the behavior of the function $`\delta \rho _s(\mathrm{\Delta },h_0)=\rho _s(\mathrm{\Delta },h_0,L)\rho _s(\mathrm{\Delta },h_0,L^{})`$ where $`L`$ and $`L^{}`$ are fixed and $`L,L^{}1`$. This corresponds, in Fig. 10, to taking two of the curves and computing the difference. If the two curves correspond to lengths that are $`L,L^{}1`$, the difference between them will start to increase considerably only when $`h_0>h_{\mathrm{crit}.}`$. This criterion assumes that for $`L1`$ $`\rho _s(L,\xi )=\rho _s(L/\xi )`$ and then, if the truncated system correlation length is much larger than $`L`$, $`\rho _s`$ becomes almost independent on the length $`L`$. In this way, this estimation is affected mainly by the truncation of the DMRG method and then, we expect to have $`\delta \rho _s`$ almost constant in the region where the truncated system correlation length is much larger than $`L=50`$ (which is the maximum length considered in the superfluid density calculations). In Fig. 11 we show $`\delta \rho _s(\mathrm{\Delta },h_0)`$ vs. $`h_0`$ (taking $`L=30`$ and $`L^{}=50`$) for $`\mathrm{\Delta }=0.5`$ and $`\mathrm{\Delta }=0.75`$. The points in these curves are obtained by averaging over disorder and have a statistical error which is shown in the figure. As it is possible to observe, for $`\mathrm{\Delta }=0.75`$ the function stays almost constant until a value of $`h_00.25`$. We have checked that this behavior is the same taking different values of $`L`$ and $`L^{}`$.

We also analyzed $`\delta \rho _s(\mathrm{\Delta },h_0)`$ for other anisotropy values. The results are shown in Fig. 12. Remarkably, all the curves corresponding to values of $`\mathrm{\Delta }`$ between $`1`$ and $`0.5`$ fall between those corresponding to $`\mathrm{\Delta }=0.5`$ and $`\mathrm{\Delta }=0.75`$ supporting the statement that the maximum critical field occurs for $`\mathrm{\Delta }=0.75`$. An estimation of the critical field is then given by the intersection of the curve with the $`h_0`$-axis. To obtain that intersection we use the first two points which differ from zero by more than the statistical error. Also from the error in these points we estimate the error in the critical field so computed.

In Fig. 13, we plot these values of the critical field in the phase diagram giving an approximation to the boundary of the critical region.

## 5 Extrapolations for $`\mathrm{\Delta }=0.75`$

In this section we study the special case of $`\mathrm{\Delta }=0.75`$ and perform the extrapolations to infinite length and infinite number of DMRG states of the computed properties for selected values of the random magnetic field.

For the superfluid density $`\rho _s`$, we considered two points on the phase diagram: ($`\mathrm{\Delta }=0.75`$, $`h_0=0.1`$) and ($`\mathrm{\Delta }=0.75`$, $`h_0=0.5`$). From the approximated phase diagram of Fig. 13 we expect $`h_0=0.1`$ to be inside the critical region and $`h_0=0.5`$ to be outside. In Fig. 14 we plot $`\rho _s`$ as a function of $`1/L`$. At large $`L`$ the two cases display opposite behavior. For $`h_0=0.1`$, $`\rho _s`$ increases as we keep more DMRG states whereas for $`h_0=0.5`$, it decreases. This is already an indication that the point ($`\mathrm{\Delta }=0.75`$, $`h_0=0.1`$) is inside the quasi-LRO region whereas the point ($`\mathrm{\Delta }=0.75`$, $`h_0=0.5`$) is outside. To find the asymptotic values $`\rho _0(s)=\rho _s(s,L\mathrm{})`$, we fitted each curve in Fig. 14 using the expression $`\rho _s(s,L)=\rho _0(s)+\alpha (s)/L^{\beta (s)}`$. The solid lines in Fig. 14 are the results of these fittings. Finally, the values of $`\rho _0(s)`$ vs. $`1/s`$ are plotted in Fig. 15. Again, we considered the expression $`\rho _0(s)=\rho _{\mathrm{}}+a/s^b`$ to fit the curves. For $`h_0=0.1`$ the extrapolated value is $`\rho _{\mathrm{}}=0.586\pm 0.003`$, confirming that the point ($`\mathrm{\Delta }=0.75`$, $`h_0=0.1`$) is inside the critical region whereas for $`h_0=0.5`$ the best fit was obtained for $`\rho _{\mathrm{}}=0`$, indicating the point ($`\mathrm{\Delta }=0.75`$, $`h_0=0.5`$) is in the localized phase as we expected from the previous analysis.

We proceeded in the same way with the correlation length. In Fig. 16 we plot the typical correlation length $`\xi `$ as a function of $`1/L`$ for $`\mathrm{\Delta }=0.75`$ and $`h_0=0.1`$, $`h_0=0.3`$ and $`h_0=0.5`$. In the $`h_0=0.1`$ case, the correlation length increases drastically as we keep more and more DMRG states. As more states are kept, the larger has to be the chain to reach the asymptotic value corresponding to the infinite system. Instead, for $`h_0=0.3`$ and $`h_0=0.5`$, the behavior is different. Not only the asymptotic values for the correlation length are one order of magnitude smaller but also the curves become closer as we increase the number of DMRG states indicating that the correlation length is finite. As we did for the superfluid density, we found the asymptotic values $`\xi _0(s)=\xi (s,L\mathrm{})`$, fitting each curve in Fig. 16 with the expression $`\xi (s,L)=\xi _0(s)+\alpha (s)/L^{\beta (s)}`$. The lines in Fig. 16 are the results of these fittings. The asymptotic values $`\xi _0(s)=\xi (s,L\mathrm{})`$ found in this way are plotted as a function of $`1/s`$ in Fig. 17. In all the cases, the errors in the data are smaller than $`1\%`$ and then, smaller than the size of the symbols.

To fit the curves in Fig. 17 we used again the expression $`\xi _0(s)=\xi _{\mathrm{}}+a/s^b`$. For $`h_0=0.1`$ the best fit was obtained for a negative exponent ($`b=0.155`$) implying an infinite correlation length for infinite number of DMRG states. On the other hand, for $`h_0=0.3`$ and $`h_0=0.5`$ we find $`b`$ possitive and the extrapolated values for the correlation length are $`\xi _{\mathrm{}}(h_0=0.3)=108\pm 18`$, and $`\xi _{\mathrm{}}(h_0=0.5)=26.45\pm 2.95`$.

In summary, from the extrapolations performed for $`\mathrm{\Delta }=0.75`$, we conclude the existence of a quasi-LRO and put bounds on the exact value of the critical field. This value should be larger than $`h_0=0.1`$ and smaller than $`h_0=0.3`$.

## 6 Conclusions

In this work we analyzed numerically the quantum phase diagram of the spin $`1/2`$ XXZ chain in a random magnetic field pointing in the $`z`$ direction.

The two-point spin-spin correlation functions in the XY plane and the superfluid density were computed using an infinite-size density matrix renormalization group (DMRG) algorithm.

The asymptotic exponential behavior of the spin correlations allowed us to extract the typical correlation length and study its dependence with the anisotropy $`\mathrm{\Delta }`$ and the strength of the random magnetic field $`h_0`$. The contour plots of the correlation length as a function of $`\mathrm{\Delta }`$ and $`h_0`$ give numerical evidence of the existence of a critical region for finite disorder and negative values of the anisotropy. In this region the correlation length increases as we keep more states in the DMRG method. Using our data for $`\rho _s`$ we extracted the behavior of the correlation length critical exponent. A good agreement with the analytical prediction was found  supporting the idea that the critical region starts in $`\mathrm{\Delta }=0.5`$. Then, an approximation to the boundary of the critical region can be obtained from the contour plots of the typical correlation length taking a contour line such that for $`\mathrm{\Delta }=0.5`$ the value of the field is $`h_00`$. In this way, we estimated $`h_{\mathrm{crit}.}(\mathrm{\Delta }=0.75)0.20`$ for $`32`$ states kept. From the analysis of the function $`\delta \rho _s(\mathrm{\Delta },h_0)=\rho _s(\mathrm{\Delta },h_0,L)\rho _s(\mathrm{\Delta },h_0,L^{})`$ for $`L,L^{}1`$, we obtained an independent estimation of the critical field $`h_{\mathrm{crit}.}(\mathrm{\Delta })`$ for several values of $`\mathrm{\Delta }`$ and drew the corresponding points of the line that separates the quasi-LRO regime from the localized one. We found that the largest value of the critical field occurs at $`\mathrm{\Delta }=0.75`$ and is approximately $`h_{\mathrm{crit}.}(\mathrm{\Delta }=0.75)=0.25`$ for $`32`$ states kept, in close agreement with the results obtained from the correlation length analysis. Since these results are approximated and are mainly affected by the truncation of the Hilbert space, to further support them, we extrapolated the typical correlation length and the superfluid density to infinite length and infinite number of DMRG states for the special case of anisotropy, $`\mathrm{\Delta }=0.75`$, and selected values of the random magnetic field. The extrapolations allow us to conclude the existence of a quasi-long-range order region and put bounds on the critical field. The range of possible values for the critical field ($`0.1<h_{\mathrm{crit}.}<0.3`$) contains the approximated results before obtained.

## Acknowledgments

We are grateful to P. Henelius and M. Gulacsi for useful discussions. This work was supported by The Swedish Natural Science Research Council, and with computing resources by The Swedish Council for Planning and Coordination of Research (FRN) and Parallelldatorcentrum (PDC), Royal Institute of Technology, Sweden.


format: text

text:
# 1 Introduction

## 1 Introduction

In this paper we analyze a class of *dissipative quantum systems*, modeling the motion of a particle ensemble, say electrons, interacting with a heat bath of oscillators. The resulting irreversible dynamics for the electrons, is a typical example of a so-called *open quantum system* \[Da\], i.e. a system in which the interaction with the environment is taken into account. The evolution equation, sometimes called *Master equation*, for the density (matrix) operator $`R(t)`$ of the particles reads

$`{\displaystyle \frac{d}{dt}}R={\displaystyle \frac{i}{\overline{}h}}[H,R]+A(R),`$ (1.1)
$`R(t=0)=R_0,`$ (1.2)
where $`H`$ is the self-adjoint *Hamiltonian operator* of the free system and $`A`$ models the effects introduced by the heat bath. The right hand side of (1.1) constitutes the formal generator $`L`$ of the *quantum dynamical semigroup* acting on $`R`$.

Assuming that the time evolution of the system is *Markovian*, G. Lindblad \[Li\] gave the most general form of a bounded operator $`L`$, such that the semigroup preserves the positivity, hermiticity and the normalization (unit trace) of the density operator $`R`$. However, for unbounded operators $`L`$, which is the case in our work, the so-called *Lindblad condition* is necessary but *not sufficient* to guarantee the conservation of these properties (see, e.g., \[CF\] and the references therein).

Using the *Wigner transform* \[Wi\], dissipative quantum models can be equivalently represented in *phase space*, resulting in a *kinetic transport equation* with interaction terms for the *quasiprobability distribution* of the particles. In this paper we assume that the mechanism coupling particles and environment can be described by a linear scattering operator $`_q`$ of *Fokker-Planck type*. The Lindblad condition therefore reduces to the assumption of a *positive definite diffusion matrix* $`D`$.

Up to now, a mathematically rigorous derivation of such a *Quantum Fokker-Planck equation* (QFP) from many-body quantum mechanics is still missing. To the authors knowledge, the only result in this direction is given in \[CEFM\], which however justifies only a particular case of the class of models considered in this work. Nevertheless there exists a huge amount of a somewhat phenomenological physical literature on this type of equations, which play a relevant role within the areas of quantum optics (laser physics) \[Da\], \[De\], \[Ri\], microelectronics \[St\], quantum brownian motion \[CaLe\], \[Di\], \[HuMa\], and the description of decoherence and diffusion of quantum states \[AnHa\], \[DGHP\].

Rigorous well-posedness and existence of local in time solutions of the frictionless QFP equation, with self-consistent Coulomb interaction, have been studied in a precedent paper of one of the authors \[ALMS\]. The present work investigates the long time behavior of the linear QFP equation in the presence of friction and an exterior time-independent potential $`V`$. The phase space formalism provided by the Wigner transform proves to be particularly helpful for this task, since it allows the use of certain *entropy techniques* established for classical dissipative equations (for an overview on these techniques, see, e.g., \[MaVi\], \[AMTU\]). The word ”entropy” is used here in a mathematically technical sense and can be seen as a generalization of the classical entropy concept of L. Boltzmann. It should *not* be confused with the physical correct *von Neumann entropy* of quantum states.

By comparison with the classical *Fokker-Planck equation* (FP), we expect the solution of the QFP equation to approach a thermal *equilibrium state* in the long time limit, provided the friction term appearing in $`_q`$ is positive and the exterior potential $`V`$ is *confining*, i.e. $`V(x)\mathrm{}`$ as $`|x|\mathrm{}`$, fast enough. In our work we shall specify the potential $`V`$ to be *harmonic*. This particular choice allows explicit calculations and is maybe the most fundamental one, from a physical point of view \[CEFM\]. Using the entropy approach, we will prove the convergence of the solution towards the steady state with a precise exponential rate, under the assumption that the initial data $`w_0`$ has bounded ”entropy,” relatively to the equilibrium state.

This paper is organized as follows. In section 2 we set up the model and collect some preliminaries. In section 3 we specify the potential $`V`$ to be of harmonic oscillator type and explicitly calculate the corresponding equilibrium state. Exponential decay towards it will be proved in section 4, where we also give precise dispersion estimates in the unconfined cases.

## 2 The model: preliminaries

We consider a linear dissipative equation modeling the motion of particles, say electrons, under the influence of an electric scalar potential $`V`$ and a thermal bath of harmonic oscillators in thermal equilibrium. In the following we denote by $`\rho (t,)L^2(^{2d})`$ the N-particle *density matrix function* of the electrons, which is the kernel of the self-adjoint trace-class *density (matrix) operator* $`R(t):L^2(^d)L^2(^d)`$, i.e.

$$(R(t)f)(x):=_^d\rho (x,y,t)f(y)𝑑y.$$
(2.1)
The evolution equation of the electrons is given as a PDE for the density matrix

$`_t\rho =`$ $`{\displaystyle \frac{i}{\overline{}h}}(H_xH_y)\rho \gamma (xy)(_x_y)\rho `$
$`+\left(D_{qq}|_x+_y|^2{\displaystyle \frac{D_{pp}}{\overline{}h^2}}|xy|^2+{\displaystyle \frac{2iD_{pq}}{\overline{}h}}(xy)(_x+_y)\right)\rho ,`$ (2.2)
subject to the initial condition

$`\rho (t=0,x,y)=\rho _0(x,y),x,y^d.`$ (2.3)
Here $`H_x`$ (resp. $`H_y`$) denotes the electron *Hamiltonian*

$$H_x:=\frac{\overline{}h^2}{2m}\mathrm{\Delta }_x+V(x),$$
(2.4)
acting on the $`x`$ (resp. $`y`$) variable. The constant $`m`$ stands for the mass of the individual particles. Equation (2) is a generalization of the Caldeira-Leggett master equation for medium temperatures, see \[Di\], \[De\].

On a kinetic level this model reads (QFP equation)

$`_tw+\xi _xw+\mathrm{\Theta }[V]w=`$ $`_qwx,\xi ^d,t^+`$ (2.5)
$`w(t=0,x,\xi )=`$ $`w_0(x,\xi ),`$ (2.6)
where the scattering operator $`_q`$ is defined by

$$_qw:=\frac{D_{pp}}{m^2}\mathrm{\Delta }_\xi w+2\gamma \mathrm{div}_\xi (\xi w)+D_{qq}\mathrm{\Delta }_xw+2\frac{D_{pq}}{m}\mathrm{div}_x(_\xi w).$$
(2.7)
Here $`w(t,)L^2(_x^d\times _\xi ^d)`$ is the *Wigner transform* \[Wi\] of the corresponding density matrix $`\rho (t,)L^2(_x^d\times _y^d)`$, i.e.

$$w(x,\xi ,t):=\frac{1}{(2\pi )^d\overline{}h}_^d\rho (x+\frac{\overline{}h}{2m}y,x\frac{\overline{}h}{2m}y,t)e^{i\xi y}𝑑y.$$
(2.8)
In the literature $`w(t,)`$ is often referred to as a quasiprobability distribution because it generally assumes negative values too \[Fo\], \[Hu\].

The (non-local) pseudo-differential operator $`\mathrm{\Theta }[V]`$ is defined by

$`\mathrm{\Theta }[V]f(x,\xi ):={\displaystyle \frac{i}{(2\pi )^d\overline{}h}}{\displaystyle _{^d\times ^d}\left[V\left(x+\frac{\overline{}h}{2m}y\right)V\left(x\frac{\overline{}h}{2m}y\right)\right]}`$
$`f(x,\xi ^{})e^{iy(\xi \xi ^{})}d\xi ^{}dy.`$ (2.9)
In order to be consistent with the usual density matrix formulation of *open quantum systems* \[Da\] in the class of *Lindblad operators* \[Li\], we assume for the *diffusion constants* $`D_{pp}>0`$ and $`D_{qq},D_{pq}0`$. On the other hand, the *friction parameter* $`\gamma `$ has to be nonnegative. Additionally we impose the following relation (for more details see \[ALMS\], \[Li\])

$$D_{pp}D_{qq}D_{pq}^2\frac{\overline{}h^2\gamma ^2}{4}\text{and}D_{pp}>0\text{if}\gamma =0.$$
(A1)
Note that if $`\gamma >0`$, condition (A1), the so-called *Lindblad condition*, implies that (2.7) is a uniformly elliptic operator.

Using the Wigner transform the *charge* and *flux densities* associated to the density matrix $`\rho `$ can be defined (formally, since $`wL^1`$ in general) in the same way as in classical statistical mechanics. Namely they are given as moments of the Wigner transfom, (for details see, e.g., \[GaMa\], \[LiPa\]),

$`n(x,t):={\displaystyle _^d}w(x,\xi ,t)𝑑\xi ,`$ (2.10)
$`j(x,t):={\displaystyle _^d}\xi w(x,\xi ,t)𝑑\xi .`$ (2.11)
Although not obvious from the above definition, the necessary positivity of $`n`$ is guaranteed for *physical quasiprobabilities* $`w`$, i.e. for $`w`$’s, which are indeed the Wigner transformed kernels $`\rho `$ of density operators $`R`$ (see \[Ar\], \[LiPa\], \[Ta\] for a more complete account on this).

With the above definitions we obtain, after formally integrating the QFP equation (2.5) w.r.t. $`\xi `$, the associated ”*continuity equation*

$$_tn+\mathrm{div}j=D_{qq}\mathrm{\Delta }_xn,$$
(2.12)
which obviously implies (again on a formal level) the *conservation of mass*, i.e.

$$M(t):=_{^d\times ^d}w(x,\xi ,t)𝑑x𝑑\xi =_{^d\times ^d}w_0(x,\xi )𝑑x𝑑\xi .$$
(2.13)
In view of this property, we assume for simplicity that

$$M_0=_{^d\times ^d}w_0(x,\xi )𝑑x𝑑\xi =1.$$
(A2)
This is of course not a restriction as long as the equation is linear. The continuity equation (2.12) suggests that in this dissipative model the usual definition of the flux density (2.11) needs to be replaced by

$`J(x,t):`$ $`=j(x,t)D_{qq}_xn(x,t)`$ (2.14)
$`={\displaystyle _^d}(\xi D_{qq}_x)w(x,\xi ,t)𝑑\xi .`$ (2.15)
Thus, instead of (2.12), we obtain the following conservation law

$$_tn+\mathrm{div}J=0$$
(2.16)
associated to the QFP equation.

Remark. In physical units, the friction and diffusion constants are usually given by (see, e.g., \[Di\], \[De\])

$$\gamma =\frac{\lambda }{2m},D_{pp}=\lambda k_BT,D_{qq}=\frac{\lambda \overline{}h^2}{12m^2k_BT},D_{pq}=\frac{\lambda \mathrm{\Omega }\overline{}h^2}{12\pi mk_BT}.$$
(2.17)
Here $`\lambda >0`$ is the coupling constant of the heat bath, $`k_B`$ the Boltzmann constant, $`T`$ the temperature of the bath and $`\mathrm{\Omega }`$ the cut-off frequency of the reservoir oscillators. In terms of these constants the Lindblad condition (A1) reads

$$\mathrm{\Omega }\frac{k_BT}{\overline{}h}$$
(2.18)
which implies the validity of our model at medium/high temperatures.

In the classical limit $`\overline{}h0`$: $`D_{qq}=D_{pq}=0`$ and (at least formally) the pseudo-differential operator simplifies to

$$\mathrm{\Theta }[V]f\frac{1}{m}_xV_\xi f.$$
(2.19)
Thus we recover the well known kinetic Fokker-Planck equation for the limiting classical phase space probability distribution $`w^{cl}(t,)^+(_x^d\times _\xi ^d)`$ (the cone of positive bounded Borel measures)

$$_tw^{cl}+\xi _xw^{cl}\frac{1}{m}_xV_\xi w^{cl}=\frac{D_{pp}}{m^2}\mathrm{\Delta }_\xi w^{cl}+2\gamma \mathrm{div}_\xi (\xi w^{cl}).$$
(2.20)
Also, the flux density $`J`$ formally simplifies to the classical one, i.e. $`Jj`$, as $`\overline{}h0`$. For a rigorous theory of such *homogenization limits* see, e.g., \[LiPa\], \[GMMP\] and for details on the extensively studied FP equation, see for example \[Ri\].

Furthermore note that the Lindblad condition (A1) disqualifies the classical FP scattering operator \[CEFM\] as a relevant quantum mechanical model for the environment interaction.

In this work we are concerned with the following solution concept for our IVP.

###### Definition 2.1.

A function $`wC(_0^+;L^p(_x^d\times _\xi ^d))`$, with $`1p<\mathrm{}`$, is a mild solution of the IVP (2.5), (2.6) with $`w_0L^p(_x^d\times _\xi ^d)`$, if and only if

$$w(t,x,\xi )=_{^d\times ^d}w_0(x_0,\xi _0)G(t,x,\xi ,x_0,\xi _0)𝑑x_0𝑑\xi _0,$$
(2.21)
where the Green’s function $`G`$ satisfies equation (2.5) for all fixed $`(x_0,\xi _0)^{2d}`$, all $`t>0`$, with an initial condition

$`\underset{t0}{lim}G(t,x,\xi ,x_0,\xi _0)=\delta (xx_0,\xi \xi _0).`$ (2.22)
that has to be understood in a weak sense. If additionally $`wC^1(^+;C^2(_x^d\times _\xi ^d))`$, then $`w`$ is called a classical solution.

Remark. In classical physics the theory of kinetic equations focuses on $`L^1`$-solutions. This implies mass conservation, since classical phase-space distributions are pointwise positive functions. Having in mind the definition of the Wigner transform, the $`L^2`$-norm is more convenient in our quantum mechanical context.

## 3 Harmonic oscillator potentials

In the next two sections we choose a normalization such that $`\overline{}h=m=1`$, for simplicity. We moreover assume the confining potential to be of the following class

$$V(x)=\frac{\omega _0^2}{2}|x|^2+ax+b,a,b,\omega _00.$$
(A3)
An easy calculation shows that, maybe after an appropriate shift in the $`x`$-variable, the pseudo-differential operator $`\mathrm{\Theta }[V]`$ is given by

$$\mathrm{\Theta }[V]w=\omega _0^2x_\xi w.$$
(3.1)
The QFP equation thus simplifies to

$`_tw+\xi _xw\omega _0^2x_\xi w=`$ $`_qwx,\xi ^d,t^+`$ (3.2)
$`w(t=0,x,\xi )=`$ $`w_0(x,\xi ),`$ (3.3)
which can be equivalently written in the more compact form

$`_tw=\mathrm{div}_{(x,\xi )}\left(D_{(x,\xi )}w+P(x,\xi )w\right),`$ (3.4)
$`w(0,x,\xi )=w_0(x,\xi ),`$ (3.5)
where the *diffusion matrix* $`D`$ and the vector-valued *drift* $`P`$ are given by

$$D:=\left(\begin{array}{cc}D_{qq}\mathrm{II}_\mathrm{d}& D_{pq}\mathrm{II}_\mathrm{d}\\ D_{pq}\mathrm{II}_\mathrm{d}& D_{pp}\mathrm{II}_\mathrm{d}\end{array}\right),P(x,\xi ):=\left(\begin{array}{c}\xi \\ \omega _0^2x+2\gamma \xi \end{array}\right).$$
(3.6)
Here $`\mathrm{II}_\mathrm{d}`$ denotes the idendity matrix in $`^d`$. Note that the Lindblad condition (A1) guarantees that the diffusion matrix is positive definite if $`\gamma >0`$ and then (3.2) is parabolic in the phase-space coordinates $`(x,\xi )^{2d}`$.

### 3.1 Fundamental solution

Now consider only the first order part of the operator (3.2), resp. (3.4). The associated characteristic ODE’s are given by

$`\dot{X}=\xi ,`$ $`X(t=0)=x_0,`$ (3.7)
$`\dot{\xi }=(\omega _0^2X+2\gamma \xi ),`$ $`\xi (t=0)=\xi _0.`$ (3.8)
This system defines the characteristic flow $`\mathrm{\Phi }_t(x_0,\xi _0)=[X_t(x_0,\xi _0),\dot{X}_t(x_0,\xi _0)]`$ in phase space $`_x^d\times _\xi ^d`$. This flow can be explicitly calculated, depending on the size of the friction constant $`\gamma `$.

###### Lemma 3.1.

Consider the system (3.7), (3.8) with $`\gamma 0`$, then the dissipative flow $`\mathrm{\Phi }_t:^{2d}^{2d}`$ is given by:

1. If $`0\gamma <\omega _0`$, then, with $`\omega :=\sqrt{\omega _0^2\gamma ^2}`$,

$`\mathrm{\Phi }_t(x_0,\xi _0)={\displaystyle \frac{e^{\gamma t}}{\omega }}[`$ $`x_0(\omega \mathrm{cos}(\omega t)+\gamma \mathrm{sin}(\omega t))+\xi _0\mathrm{sin}(\omega t),`$
$`\xi _0(\omega \mathrm{cos}(\omega t)\gamma \mathrm{sin}(\omega t))x_0\omega _0^2\mathrm{sin}(\omega t)].`$ (3.9)
2. If $`\gamma >\omega _0`$, then, with $`\omega :=\sqrt{\gamma ^2\omega _0^2}`$,

$`\mathrm{\Phi }_t(x_0,\xi _0)={\displaystyle \frac{e^{\gamma t}}{\omega }}[`$ $`x_0(\omega \mathrm{cosh}(\omega t)+\gamma \mathrm{sinh}(\omega t))+\xi _0\mathrm{sinh}(\omega t),`$
$`\xi _0(\omega \mathrm{cosh}(\omega t)\gamma \mathrm{sinh}(\omega t))x_0\omega _0^2\mathrm{sinh}(\omega t)].`$ (3.10)
3. If $`\gamma =\omega _0`$ then

$$\mathrm{\Phi }_t(x_0,\xi _0)=e^{\gamma t}[(\gamma t+1)x_0+t\xi _0,(1\gamma t)\xi _0\gamma ^2tx_0].$$
(3.11)

###### Proof.

The proof follows from straightforward calculations. ∎

Now, using lemma 2.1, we obtain an explicit representation of Green’s function $`G`$ of the QFP equation with harmonic oscillator potential, depending on the size of $`\gamma `$.

###### Proposition 3.1.

Let $`\gamma 0`$ and let conditions (A1), (A3) hold, then the Green’s function $`G`$ associated to (3.2) is, for every fixed $`t>0`$, a pointwise positive function, given by

$$G(t,x,\xi ,x_0,\xi _0)=e^{2d\gamma t}F(t,X_t(x,\xi )x_0,\dot{X}_t(x,\xi )\xi _0),$$
(3.12)
with

$$F(t,x,\xi )=\frac{\mathrm{exp}\left(\frac{\nu (t)|x|^2+\lambda (t)|\xi |^2+\mu (t)(x\xi )}{4\lambda (t)\nu (t)\mu ^2(t)}\right)}{(2\pi )^d(4\lambda (t)\nu (t)\mu ^2(t))^{d/2}}C^1(^+;𝒮(_x^d\times _\xi ^d)).$$
(3.13)
In (3.12) we denote by $`X_t,\dot{X}_t`$ the components of the inverse characteristic flow $`\mathrm{\Phi }_t`$, which satisfies $`\mathrm{\Phi }_t\mathrm{\Phi }_t=\mathrm{id}`$:

$$\mathrm{\Phi }_t(x,\xi )=[X_t(x,\xi ),\dot{X}_t(x,\xi )].$$

The associated functions $`\lambda ,\nu ,\mu :^+`$ are defined by the following expressions:

$`\lambda (t)`$ $`:={\displaystyle _0^t}\left(D_{qq}\alpha ^2(s)+D_{pp}\beta ^2(s)+2D_{pq}\alpha (s)\beta (s)\right)𝑑s,`$ (3.14)
$`\nu (t)`$ $`:={\displaystyle _0^t}\left(D_{qq}\dot{\alpha }^2(s)+D_{pp}\dot{\beta }^2(s)+2D_{pq}\dot{\alpha }(s)\dot{\beta }(s)\right)𝑑s,`$ (3.15)
$`\mu (t)`$ $`:=2{\displaystyle _0^t}\left(D_{qq}\alpha (s)\dot{\alpha }(s)+D_{pp}\beta (s)\dot{\beta }(s)+D_{pq}{\displaystyle \frac{d}{ds}}(\alpha (s)\beta (s))\right)𝑑s,`$ (3.16)
where the functions $`\alpha ,\beta :^+`$ are given by:

$`\alpha (t):=d^1\mathrm{div}_x\left(X_t(x,\xi )\right),\beta (t):=d^1\mathrm{div}_\xi \left(X_t(x,\xi )\right).`$ (3.17)
###### Proof.

The proof is similar to the calculations given in \[Bo\], \[Ho\]. First note that by definition of $`\alpha `$ and $`\beta `$, we can write $`X_t(x,\xi )=\alpha (t)x+\beta (t)\xi `$. Thus if $`G`$ is the fundamental solution of (3.2), the linear transformation (3.12) guarantees that $`F`$ has to be a fundamental solution of the following PDE with time dependent coefficients

$$_tF(t,x,\xi )=\left(\frac{d\lambda }{dt}(t)\mathrm{\Delta }_x+\frac{d\nu }{dt}(t)\mathrm{\Delta }_\xi \frac{d\mu }{dt}(t)_x_\xi \right)F(t,x,\xi ).$$

where the functions $`\lambda ,\nu ,\mu `$ are calculated depending on the choice of $`\gamma `$. A Fourier transform now shows that

$$(F)(t,k,\eta )\widehat{F}(t,k,\eta ):=_{^d\times ^d}F(t,x,\xi )e^{i(xk+\xi \eta )}𝑑x𝑑\xi $$

is a solution of

$$_t\widehat{F}(t,k,\eta )=\left(\frac{d\lambda }{dt}(t)|k|^2+\frac{d\nu }{dt}(t)|\eta |^2\frac{d\mu }{dt}(t)(k\eta )\right)\widehat{F}(t,k,\eta ).$$

This equation can easily be integrated and thus, using that $`\widehat{F}(t,0,0)=1`$ for all $`t0`$, we obtain

$`\mathrm{ln}\left(\widehat{F}(t,k,\eta )\right)`$ $`=\left(\lambda (t)|k|^2+\nu (t)|\eta |^2\mu (t)(k\eta )\right).`$
After some lengthy and tedious calculations (where one checks that $`4\lambda \nu \mu ^2`$), an inverse Fourier transform gives

$`F(t,x,\xi )`$ $`=(2\pi )^{2d}{\displaystyle _{^d\times ^d}e^{i(xk+\xi \eta )}\widehat{F}(t,k,\eta )𝑑k𝑑\eta }`$
$`=(2\pi )^d{\displaystyle \frac{e^{\frac{|x|^2}{4\lambda }}}{(4\pi \lambda )^{d/2}}}{\displaystyle _^d}e^{i\frac{\mu }{2\lambda }(x\eta )\frac{1}{4\lambda }(4\lambda \nu \mu ^2)|\eta |^2+i\xi \eta }𝑑\eta `$
$`={\displaystyle \frac{\mathrm{exp}\left(\frac{\nu (t)|x|^2+\lambda (t)|\xi |^2+\mu (t)(x\xi )}{4\lambda (t)\nu (t)\mu ^2(t)}\right)}{(2\pi )^d(4\lambda (t)\nu (t)\mu ^2(t))^{d/2}}},`$
which is the desired result. ∎

Note that the quantum mechanical effects in $`F`$ and consequently in $`G`$ only enter in form of the constants $`D_{qq},D_{pq}\overline{}h^2`$, which appear in the auxiliary functions $`\lambda ,\nu ,\mu `$. In other words we obtain Green’s function for the classical Fokker-Planck equation in a square-well potential by setting $`D_{qq}=D_{pq}=0`$ in the above expressions. The pointwise positivity of $`G`$ is a consequence of the minimum principle for parabolic equations of Fokker-Planck type \[Ev\].

From the above proposition we draw the following consequences (among which we obtain the conservativity of the quantum dynamical semigroup, associated to the harmonically confined QFP equation).

###### Corollary 3.1.

Let $`\gamma 0`$ and assume (A1)-(A3). Then for every initial condition $`w_0L^p(_x^d\times _\xi ^d)`$ with $`1p<\mathrm{}`$, there exists a unique classical solution $`wC(_0^+;L^p(_x^d\times _\xi ^d))`$ $``$ $`C^1(^+;C^{\mathrm{}}(_x^d\times _\xi ^d))`$ with $`M(t)=M_01`$. Moreover if $`w_0`$, $`n_0`$ are non-negative a.e., so are $`w(t,)`$, $`n(t,)`$, for all $`t^+`$.

###### Proof.

For convenience we use the notation: $`y:=(x,\xi )`$ as well as $`y_0:=(x_0,\xi _0)`$. With the following linear change of variables

$$G(t,\mathrm{\Phi }_t(y),y_0)=e^{2d\gamma t}F(t,yy_0)C^1(^+,𝒮(^{2d})),$$

we can express our solution in the form

$$w(t,\mathrm{\Phi }_t(y))=e^{2d\gamma t}(w_0F(t,y)).$$

A straightforward computation now shows that the Jacobian determinant of the mapping $`\mathrm{\Phi }_t()`$ is given by

$$det\left(\frac{\mathrm{\Phi }_t(y)}{y}\right)=\mathrm{exp}(2d\gamma t).$$

With these preparations and using Young’s inequality \[LiLo\], we obtain

$`w(t)_p=w_0F(t)_pw_0_pF(t)_1<\mathrm{},`$
since for each fixed $`t^+`$: $`F(t,)𝒮(^{2d})L^1(^{2d})`$. More precisely we have

$$F(t)_p^p=_{^{2d}}|F(t,y)|^p𝑑y=p^1,$$

for all $`1p<\mathrm{}`$. For $`p=1`$ this implies, after a simple calculation, that $`M(t)=M_01`$, by (A2). Since $`G(t,)`$ is pointwise positive, we clearly obtain that if $`w_0`$, resp. $`n_0`$ are a.e. non-negative, so are $`w(t,)`$, $`n(t,)`$. ∎

Remark. If the initial condition $`w_0`$ is the Wigner transform of a pure quantum state $`\psi _0`$, it is well known \[Hu\], \[LiPa\] that $`w_00`$ pointwise, if and only if $`\psi _0`$ is a Gaussian. A similar characterization for mixed states has not been found yet.

From the above result it is easy to deduce the *conservativity* (i.e. conservation of hermiticity, positivity and normalization of the density matrix $`\rho `$) of the quantum dynamical semigroup corresponding to (3.2), using the inverse Wigner transform

$$\rho (x,y,t)=(2\pi )^d_^dw(\frac{x+y}{2},\xi ,t)e^{i\xi (xy)}𝑑\xi ,$$
(3.18)
which is defined in the sense of the usual $`L^2`$-Fourier transform \[Fo\], \[LiLo\].

### 3.2 Stationary states.

As in the classical case, we expect that the competing effects of the confining potential plus the positive friction and the dissipating behavior of the operator result in a *thermal equilibrium state* as $`t\mathrm{}`$.

###### Definition 3.1.

A (thermal) equilibrium state is a steady state, i.e. stationary solution $`wL^2(_x^d\times _\xi ^d)`$ of the QFP equation (2.5), such that

$$J(x)=_^d(\xi D_{qq}_x)w(x,\xi )𝑑\xi =0.$$
(3.19)
In principle it could be possible that there exist stationary solutions of the QFP equation, which are not equilibrium states, however the following proposition and corollary show that this is not the case.

Furthermore, we shall see that the physical intuitive assumption, that the mass of the initial state is equal to the mass of the steady state, guarantees its uniqueness.

###### Proposition 3.2.

Let $`\gamma >0`$, $`\omega _0>0`$ and assume (A1), (A3), then the unique solution $`w_{\mathrm{}}`$ of the stationary QFP equation, satisfying

$$M_{\mathrm{}}:=_^d_^dw_{\mathrm{}}(x,\xi )𝑑x𝑑\xi =1,$$
(3.20)
is given by the following non-isotropic Gaussian function

$$w_{\mathrm{}}(x,\xi )=\frac{\gamma \omega _0}{(2\pi )^d\sqrt{Q}}\mathrm{exp}\left(\frac{\gamma }{Q}\left[Q_{11}\omega _0^2|x|^2+2Q_{12}\omega _0x\xi +Q_{22}|\xi |^2\right]\right),$$
(3.21)
such that $`w_{\mathrm{}}𝒮(_x^d\times _\xi ^d)`$, with $`Q:=Q_{11}Q_{22}Q_{12}^2`$ and

$`Q_{11}`$ $`:=D_{pp}+\omega _0^2D_{qq},`$ (3.22)
$`Q_{12}`$ $`:=2\omega _0\gamma D_{qq},`$ (3.23)
$`Q_{22}`$ $`:=D_{pp}+\omega _0^2D_{qq}+4\gamma (D_{pq}+\gamma D_{qq}).`$ (3.24)
###### Proof.

The proof is much simpler in Fourier-space where the steady state $`\widehat{w}_{\mathrm{}}`$ is explicitly given by

$`\widehat{w}_{\mathrm{}}(k,\eta )=\mathrm{exp}`$ $`({\displaystyle \frac{D_{pp}}{4\gamma \omega _0^2}}(|k|^2+\omega _0^2|\eta |^2){\displaystyle \frac{D_{pq}}{\omega _0^2}}|k|^2`$
$`D_{qq}[({\displaystyle \frac{\gamma }{\omega _0^2}}+{\displaystyle \frac{1}{4\gamma }})|k|^2+{\displaystyle \frac{\omega _0^2}{4\gamma }}|\eta |^2+k\eta ]).`$
First it is straightforward to check that this is indeed a solution of the Fourier transformed stationary QFP equation

$$\omega _0^2\eta _k\widehat{w}_{\mathrm{}}+(2\gamma \eta k)_\eta \widehat{w}_{\mathrm{}}=\left(D_{pp}|\eta |^2+D_{qq}|k|^2+2D_{pq}k\eta \right)\widehat{w}_{\mathrm{}}$$

and then that it satisfies the Fourier transformed mass normalization condition $`w_{\mathrm{}}𝑑x𝑑\xi =1`$, i.e.

$$\widehat{w}_{\mathrm{}}(0,0)=1.$$
(3.25)
Thus, it remains to prove the uniqueness of the steady state. To do so we consider the corresponding characteristic system

$`\dot{k}=\omega _0^2\eta ,`$ $`k(s=0)=k_0`$ (3.26)
$`\dot{\eta }=2\gamma \eta k,`$ $`\eta (s=0)=\eta _0.`$ (3.27)
Having in mind that, by assumption, $`\gamma >0`$ we check that the real parts of the eigenvalues $`\lambda _{1,2}:=\gamma \pm \sqrt{\gamma ^2\omega _0^2}`$ of this system are positive, which, by the standard theory of ODE’s, implies that $`(k,\eta )=(0,0)`$ is a source of the characteristic flow (a similar argument holds if $`\gamma =\omega _0`$). Thus the condition (3.25) is necessary and sufficient to guarantee the uniqueness of the solution $`\widehat{w}_{\mathrm{}}`$. ∎

From the above proposition we draw the following consequences:

###### Corollary 3.2.

Under the same assumptions as above we have:

1. The state $`w_{\mathrm{}}`$ is also an equilibrium state, i.e. $`J_{\mathrm{}}(x)=0`$ for all $`x^d`$.
2. If either $`\omega _0=0`$ (free motion case) or $`\gamma =0`$ (frictionless case) or both are equal to zero, no nontrivial $`L^1`$-steady state exists.

###### Proof.

It is a lengthy but straightforward calculation to show that the current-density $`J_{\mathrm{}}`$ associated to $`w_{\mathrm{}}`$ vanishes identically.

In the limiting cases $`\gamma =0`$ or $`\omega _0=0`$, we consider the characteristic curves given by the ODE system (3.26), (3.27). Along these curves the steady state varies according to

$$\frac{d}{ds}\widehat{w}_{\mathrm{}}(k_s,\eta _s)=\vartheta _s(k_0,\eta _0)\widehat{w}_{\mathrm{}}(k_s,\eta _s),$$
(3.28)
where $`\eta _s=\eta _s(k_0,\eta _0)`$, $`k_s=k_s(k_0,\eta _0)`$ denote the (vector valued) characteristic curves starting for $`s=0`$ at the point $`(k_0,\eta _0)`$. By $`\vartheta _s(k_0,\eta _0)`$, we mean

$`\vartheta _s(k_0,\eta _0):=D_{pp}|\eta _s(k_0,\eta _0)|^2+D_{qq}|k_s(k_0,\eta _0)|^2+2D_{pq}(k_s\eta _s)(k_0,\eta _0).`$
We integrate equation (3.28) and obtain the Fourier transformed steady state parametrized by $`s`$

$`\widehat{w}_{\mathrm{}}(k_s(k_0,\eta _0),\eta _s(k_0,\eta _0))=\widehat{w}_{\mathrm{}}(k_0,\eta _0)\mathrm{exp}\left({\displaystyle _0^s}\vartheta _\tau (k_0,\eta _0)𝑑\tau \right).`$ (3.29)
In case the potential vanishes, i.e. $`\omega _0=0`$ (and $`\gamma 0`$), we have for all $`s`$ that $`k_s(k_0,\eta _0)=k_0`$ (free motion) and

$$\eta _\tau (k_0,\eta _0)=\{\begin{array}{cc}\eta _0k_0\tau \hfill & \text{if }\gamma =0,\hfill \\ \eta _0\mathrm{exp}(2\gamma \tau )\frac{k_0}{2\gamma }\hfill & \text{if }\gamma >0\text{.}\hfill \end{array}$$

Now let $`\omega _0=\gamma =0`$ and set $`k_0=0`$ in equation (3.29). We obtain

$$\widehat{w}_{\mathrm{}}(0,\eta _0)=\widehat{w}_{\mathrm{}}(0,\eta _0)\mathrm{exp}\left(D_{pp}|\eta _0|^2s\right),\eta _0^d,s,$$

which implies

$$\widehat{w}_{\mathrm{}}(0,\eta _0)=0,\eta _0^d,$$

in contradiction to (3.25). The same type of argument holds in the case $`\omega _0=0`$, $`\gamma >0`$, if we set $`\eta _0=k_0/2\gamma `$ in equation (3.29).

In the frictionless case, i.e. $`\gamma =0`$, $`\omega _0>0`$, the characteristics are circles in the $`(k,\eta )`$ plane. Thus for $`s=2\pi /\omega _0`$ we have $`(k_s(k_0,\eta _0),\eta _s(k_0,\eta _0))=(k_0,\eta _0)`$. However if we integrate over one period, i.e. setting $`s=2\pi /\omega _0`$ in (3.29), we obtain for all $`k_0,\eta _0^d`$

$$\widehat{w}_{\mathrm{}}(k_0,\eta _0)=\widehat{w}_{\mathrm{}}(k_0,\eta _0)\mathrm{exp}\left(\frac{\pi }{\omega _0^3}(D_{pp}+\omega _0^2D_{qq})(|k_0|^2+\omega _0^2|\eta _0|^2)\right)$$

which clearly implies, since $`D_{pp},D_{qq}>0`$, that $`\widehat{w}_{\mathrm{}}`$ identically vanishes. ∎

The above corollary in particular shows that, as expected, the presence of a positive friction together with a confining potential are crucial to guarantee the existence of an admissible, i.e. with finite mass, equilibrium state.

Note that although the solution $`w`$ of (2.5) in general will not be nonnegative, the steady state $`w_{\mathrm{}}`$ is nevertheless a pointwise *positive* function, because of Assumption (A2). The associated density matrix $`\rho _{\mathrm{}}`$ is a particular example of a *mixed* quantum state with positive Wigner transform. It is explicitely given by

$$\rho _{\mathrm{}}(x,y)=\frac{\gamma \omega _0}{(16\pi ^3)^{d/2}\sqrt{\gamma Q_{22}}}e^{\frac{1}{4\gamma Q_{22}}\left[\gamma ^2\omega _0^2(x+y)^2+Q(xy)^2\right]+i\omega _0\frac{Q_{12}}{Q_{22}}\left(\frac{x^2y^2}{2}\right)}.$$
(3.30)
Observe that the equilibrium density $`n_{\mathrm{}}`$, which is obtained from the diagonal (i.e. $`x=y`$) of the steady state density matrix $`\rho _{\mathrm{}}`$, is a real valued Gaussian function.

Remark. As a special limit case, we observe that in the classical limit, i.e. $`D_{qq}=D_{pq}=0`$, the equilibrium state simplifies to the well known stationary solution $`w_{\mathrm{}}^{cl}`$ of the classical kinetic Fokker-Planck equation (2.20) with a harmonic oscillator potential (c.f. \[Ri\]), i.e.

$$w_{\mathrm{}}^{cl}(x,\xi )=\frac{\omega _0\gamma }{(2\pi )^dD_{pp}}e^{\frac{\gamma }{D_{pp}}(|\xi |^2+\omega _0^2|x|^2)}.$$
(3.31)
Note that, in contrast to the classical case, the quantum steady state is *not* a function of the classical energy

$$H(x,\xi ):=\frac{1}{2}|\xi |^2+\frac{1}{2}\omega _0^2|x|^2.$$
(3.32)
For the classical kinetic FP equation, the fact that $`w_{\mathrm{}}^{cl}=w_{\mathrm{}}^{cl}(H)`$, implies that the transport and the classical FP scattering operator vanish independently when applied to $`w_{\mathrm{}}^{cl}`$. In our quantum mechanical framework this is no longer true, since there the steady state $`w_{\mathrm{}}`$ results from a cancellation of the transport operator and the scattering term $`_q`$.

## 4 Long time behavior

We are now in the position to describe the long time behavior of the linear QFP equation in the cases were both friction and an external potential of harmonic oscillator type are present, and in cases where one is missing. We will start with the latter situation.

### 4.1 Dispersion estimates in the unconfined case

We want to address the unconfined or dispersive cases, i.e. either $`\gamma =0`$ or $`\omega _0=0`$ (or both). By comparison with the classical FP equation, we expect that the particles escape to infinity and thus that the macroscopic density $`n`$ decays to $`0`$ as $`t\mathrm{}`$. More precisely we have the following theorem.

###### Theorem 4.1.

Let be either $`\gamma =0`$ or $`\omega _0=0`$, or both, and assume (A1), (A3):

1. If $`w_0L^1L^2`$, then the solution $`w`$ of the QFP equation satisfies

$$w(t)_pC_pR_w(t)^{\frac{d}{2q}}w_0_1,1p\mathrm{},$$
(4.1)
where $`C_p`$ is a constant independent of $`w_0`$, $`p^1+q^1=1`$ and

$$R_w(t):=e^{4\gamma t}(4\lambda (t)\nu (t)\mu ^2(t)),$$
(4.2)
which is a pointwise positive function with $`R_w(t)\mathrm{}`$ as $`t\mathrm{}`$.
2. Consequently we have for the corresponding density, that

$$n(t)_pC_pR_n(t)^{\frac{d}{2q}}w_0_1,1p\mathrm{},$$
(4.3)
where again the rate $`R_n(t)\mathrm{}`$ as $`t\mathrm{}`$. Explicitly

$$R_n(t):=2\left(\lambda (t)\stackrel{~}{\alpha }^2(t)+\mu (t)\stackrel{~}{\alpha }(t)\stackrel{~}{\beta }(t)+\nu (t)\stackrel{~}{\beta }^2(t)\right),$$
(4.4)
using

$$\stackrel{~}{\alpha }(t):=d^1\mathrm{div}_{x_0}\left(X_t(x_0,\xi _0)\right),\stackrel{~}{\beta }(t):=d^1\mathrm{div}_{\xi _0}\left(X_t(x_0,\xi _0)\right).$$

###### Proof.

The estimate in claim no. 1 is trivial in the case $`p=\mathrm{}`$ from the formula of the fundamental solution in (3.12) and (2.21). The estimate for $`p=1`$ is due to the following property of the Green’s function

$$_{^d\times ^d}G(t,x,\xi ,x_0,\xi _0)𝑑x_0𝑑\xi _0=1$$

and a trivial estimate over (2.21). The estimates for $`1<p<\mathrm{}`$ are then obtained by interpolation. The function $`R_w(t)`$ can be computed explicitly in each of the three cases and is given by:

1. If $`\omega _0>\gamma =0`$, then with $`\phi =\phi (t)=2t\omega _0`$,

$`R_w(t)`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(D_{qq}^2+{\displaystyle \frac{1}{\omega _0^4}}D_{pp}^2\right)\left(\phi ^2+2\mathrm{cos}\phi 2\right)`$
$`+{\displaystyle \frac{2}{\omega _0^2}}D_{pq}^2\left(\mathrm{cos}\phi 1\right)+{\displaystyle \frac{1}{2\omega _0^2}}D_{pp}D_{qq}\left(\phi ^22\mathrm{cos}\phi +2\right).`$
2. If $`\gamma >\omega _0=0`$, then with $`\chi =\chi (t)=e^{2t\gamma }`$,

$`R_w(t)`$ $`=`$ $`{\displaystyle \frac{1}{4\gamma ^{\mathrm{\hspace{0.17em}4}}}}(D_{pp}^2+\mathrm{\hspace{0.17em}4}\gamma D_{pp}D_{pq})\left(t\gamma (1\chi ^2)(1\chi )^2\right)`$
$`{\displaystyle \frac{1}{\gamma ^2}}D_{pq}^2(1\chi )^2+{\displaystyle \frac{t}{\gamma }}D_{qq}(1\chi ^2)`$
3. If $`\gamma =\omega _0=0`$, then

$$R_w(t)=4D_{pq}^2t^2+4D_{pp}D_{qq}t^2+\frac{1}{3}D_{pp}^2t^4.$$

In all cases, $`R_w`$ diverges as $`t\mathrm{}`$ and thus claim no. 1 is proved.

To prove claim no. 2 we first compute the integral of the fundamental solution w.r.t. $`\xi `$, which gives

$$_^dG(t,x,\xi ,x_0,\xi _0)𝑑\xi =R_n^{d/2}(t)𝒩\left(\frac{xX_t(x_0,\xi _0)}{\sqrt{R_n(t)}}\right),$$
(4.5)
where $`X_t`$ is defined as in lemma 3.1 and

$$𝒩(\sigma ):=(2\pi )^{d/2}\mathrm{exp}\left(\frac{\sigma ^2}{2}\right).$$

We can then deduce a formula for the evolution of the macroscopic density, using (3.12), (2.21) and the above computation (4.5), to obtain

$$n(t,x)=R_n^{d/2}(t)_{^d\times ^d}𝒩\left(\frac{xX_t(x_0,\xi _0))}{\sqrt{R_n(t)}}\right)w_0(x_0,\xi _0)𝑑x_0𝑑\xi _0.$$

The dispersion estimates are then straightforward from this expression. On the other hand, the function $`R_n(t)`$ again can be computed explicitly:

1. If $`\omega _0>\gamma =0`$, then with $`\phi =\phi (t)=2t\omega _0`$,

$$2\omega _0^3R_n(t)=2\omega _0D_{pq}\left(1\mathrm{cos}\phi \right)+D_{pp}\left(\phi \mathrm{sin}\phi \right)+\omega _0^2D_{qq}\left(\phi +\mathrm{sin}\phi \right).$$
2. If $`\gamma >\omega _0=0`$, then with $`\chi =\chi (t)=e^{2t\gamma }`$,

$$8\gamma ^3R_n(t)=D_{pp}\left(4\chi +4\gamma t3\chi ^2\right)+16\gamma ^3tD_{qq}+8\gamma D_{pq}\left(\chi +2\gamma t1\right).$$
3. If $`\gamma =\omega _0=0`$, then

$$R_n(t)=2D_{qq}t+2D_{pq}t^2+\frac{2}{3}D_{pp}t^3.$$

In all three cases $`R_n`$ diverges as $`t\mathrm{}`$ and thus claim no. 2 is also proved. ∎

The behavior at $`t=+\mathrm{}`$ of the rate functions can be obtained directly from the explicit computations in the previous result, depending on the diffusion constants. In the following corollary we will consider only the physical most important case, namely $`D_{pp}>0`$.

###### Corollary 4.1.

Assume $`D_{pp}>0`$ and (A1), (A3), then we have, as $`t+\mathrm{}`$:

1. If $`\omega _0>\gamma =0`$, then $`R_w(t)=O(t^2)`$ and $`R_n(t)=O(t)`$.
2. If $`\gamma >\omega _0=0`$, then $`R_w(t)=O(t)`$ and $`R_n(t)=O(t)`$.
3. If $`\gamma =\omega _0=0`$, then $`R_w(t)=O(t^4)`$ and $`R_n(t)=O(t^3)`$.

Remark. The computation of $`n`$ can also be used to measure the convergence rates towards the steady state macroscopic density $`n_{\mathrm{}}`$ in the confined case. This computation however is quite involved and we leave the details to the reader, since we will take a more elegant approach in the next subsection.

Nevertheless the advantage of such an explicit calculation would be that one can prove exponential decay of the solution, even for the classical FP equation (this can also be seen by comparison with spectral theoretical approaches \[Ri\]), whereas the entropy method used in section 4.2 only gives a suboptimal rate in this case \[DeVi\] .

### 4.2 Exponential decay towards equilibrium

We now assume the presence of friction and of the confining potential. Like in the classical case, we expect exponential decay of the solution towards the equilibrium state, which is usually proved using spectral theory.

However in this work we shall follow a different approach, the so-called *entropy-entropy-dissipation method* for classical FP type equations (see \[AMTU\] and references therein). As we shall see it can be successfully applied in our quantum mechanical context too (for an overview of the classical applications, see \[MaVi\]).

First let us rewrite the drift-vector $`P`$ given by (3.6) in the following way

$$P(x,\xi )=D(A(x,\xi )+F(x,\xi )),$$
(4.6)
where the gradient is taken w.r.t. $`(x,\xi )^{2d}`$ and $`A`$ is defined as the normalized potential appearing in the expression of the equilibrium state, more precisly

$$w_{\mathrm{}}(x,\xi )\mathrm{exp}\left(A(x,\xi )\right).$$
(4.7)
Consequently the QFP equation (3.2) takes the standard form of a *non-sym**metric drift-diffusion equation*

$`_tw=\mathrm{div}\left[D(w+w(A+F))\right],t^+`$ (4.8)
$`w(t=0,x,\xi )=w_0(x,\xi ).`$ (4.9)
The vector field $`F`$ is explicitly given by

$`F(x,\xi )`$ $`:=D^1P(x,\xi )A(x,\xi ),`$ (4.10)
where $`D^1`$ is the inverse of the diffusion matrix, $`A`$ is defined by (4.7) and $`P`$ is the drift given in (3.6).

The rewritten QFP equation (for harmonic oscillator potentials) with equilibrium state given by (4.7) can now be identified as a special case of FP type equations, see \[AMTU\]. Therefore we introduce the concept of *relative entropies* in the same sense as in the quoted work:

###### Definition 4.1.

Let $`𝒜`$ be either $``$ or $`^+`$. Let $`\phi C(\overline{𝒜})C^4(𝒜)`$ satisfy

$`\phi (1)=0,\phi ^{\prime \prime }0\text{ with }\phi ^{\prime \prime }0\text{ and }(\phi ^{\prime \prime \prime })^2{\displaystyle \frac{1}{2}}\phi ^{\prime \prime }\phi ^{IV}\text{ on }𝒜\text{.}`$ (4.11)
Assume moreover that $`fL^1(^n)`$, $`gL_+^1(^n)`$ with $`_^nf𝑑x=_^ng𝑑x=1`$ and $`f/g\overline{𝒜}`$ $`g(dx)`$ a.e. Then

$$e_\phi (fg):=_^n\phi \left(\frac{f}{g}\right)g(dx)$$
(4.12)
is an admissible relative entropy of $`f`$ w.r.t. $`g`$ and generating function $`\phi `$.

Two relative entropies which are frequently used are the *logarithmic relative entropy*, associated to the generating function

$$\phi _1(\alpha ):=\alpha \mathrm{ln}\alpha \alpha +1,\alpha ^+.$$
(4.13)
and the *quadratic relative entropy* with generator

$$\phi _2(\alpha ):=k(\alpha 1)^2,\alpha ,k>0.$$
(4.14)
Indeed it is known, see \[AMTU\], that (up to a positive multiplicative constant) for every admissible entropy generator $`\phi `$, there exist constants $`K_1`$, $`K_2>0`$ such that for all $`\alpha ^+`$, it holds that

$$K_1\phi _1(\alpha )\phi (\alpha )K_2\phi _2(\alpha ).$$
(4.15)
Remark. The above defined entropies should not be confused with the quantum mechanical *von Neumann entropy* $`S:=\mathrm{Tr}(R\mathrm{ln}R)`$, where $`R`$ is the density operator of the particle ensemble \[Th\] (also note that in contrast to the physical convention the minus sign in front of (4.12) is dropped).

There is however a notion of entropy for quantum states, called *Wehrl entropy* \[We\], which is closely related to the logarithmic entropy defined above (using the pointwise nonnegative *Husimi transform* of $`w`$, see, e.g., \[LiPa\]) and which can be interpreted as a measure of *coherence and localization of quantum states*. For details, see \[AnHa\], \[GnZy\], \[SlZy\] and the references given therein.

Since the solution of the QFP equation in general is not pointwise positive, it seems that we need to look for an admissible entropy on all of $``$. This would imply, see \[AMTU\], that the the only admissible entropy useful for our purpose is the quadratic one (4.14). We circumvent this shortcoming by decomposing the initial condition $`w_0`$ into its negative and positive parts. More precisely we write

$$w_0(x,\xi )=w_0^+(x,\xi )w_0^{}(x,\xi )\text{ a.e.}$$
(4.16)
where now $`w_0^+`$, $`w_0^{}`$ are both non-negative functions with mass $`M_0^+`$, $`M_0^{}`$ respectively. Clearly condition (A2) implies

$$1=M_0^+M_0^{}.$$
(4.17)
Having in mind that the QFP equation is linear, we denote by

$$w^\pm (t,x,\xi ):=_{^d\times ^d}w_0^\pm (x_0,\xi _0)G(t,x,\xi ,x_0,\xi _0)𝑑x_0𝑑\xi _0$$
(4.18)
the mild solution of (4.8) corresponding to $`w_0^\pm `$. Corollary 2.1 implies that $`w^\pm (t,)0`$ a.e. and we shall now apply the entropy-entropy-dissipation method to each of the two functions $`w^+`$, $`w^{}`$. Note that the steady state associated to $`w^\pm (t,)`$ is given by $`w_{\mathrm{}}^\pm =M_0^\pm w_{\mathrm{}}`$.

In order to use the results established in \[AMTU\], we first need to check that the following property is fulfilled.

###### Lemma 4.1.

Assume $`\gamma >0`$, $`\omega _0>0`$ as well as (A1), (A3), $`y:=(x,\xi )^{2d}`$. Then it holds that

$$\mathrm{div}_y\left(DFw_{\mathrm{}}\right)=0on^{2d}$$
(4.19)
where $`F`$ is defined by (4.10).

###### Proof.

To prove the claim, we note, using equation (4.7), that condition (4.19) is equivalent to

$$\mathrm{div}_y\left(Dw_{\mathrm{}}+Pw_{\mathrm{}}\right)=0,$$

which is the stationary QFP equation for potentials of the form (A3). Thus the claim is true by definition of $`A,F`$ and $`w_{\mathrm{}}`$. ∎

With these preparations we can now state the main result of this section.

###### Theorem 4.2.

Assume $`\gamma >0`$, $`\omega _0>0`$ and assume (A1)-(A3). Let the initial data $`w_0`$ be such that $`w_0(x,\xi )=w_0^+(x,\xi )w_0^{}(x,\xi )`$ a.e., $`w_0^\pm L^1L^2`$ and assume

$$e_\phi (w_0^\pm M^\pm w_{\mathrm{}})<\mathrm{}.$$
(A4)
Then there exists a $`\kappa >0`$ such that

$$e_\phi (w^\pm (t,)M_0^\pm w_{\mathrm{}})e^{2\kappa t}e_\phi (w_0^\pm M_0^\pm w_{\mathrm{}}),t>0.$$
(4.20)
As a consequence, the solution of (3.2) converges exponentially towards the equilibrium state. More precisely it holds that:

1. If $`e_{\phi _1}(w_0^\pm w_{\mathrm{}})<\mathrm{}`$, with $`\phi _1`$ defined in (4.13), then

$$w(t,)w_{\mathrm{}}_1Ce^{\kappa t},C^+,t>0.$$
(4.21)
2. More generally, if $`e_{\phi _2}(w_0^\pm w_{\mathrm{}})<\mathrm{}`$, with $`\phi _2`$ defined in (4.14), then

$$w(t,)w_{\mathrm{}}_pCe^{\kappa t}C^+,t>0,$$
(4.22)
with $`1p2`$.

###### Proof.

The proof of the first claim (4.20) is a consequence of the following *convex Sobolev inequality* \[AMTU\] (in which we use the notation $`d\rho _{\mathrm{}}=\rho _{\mathrm{}}dx`$)

$$_{^{2d}}\phi \left(\frac{\rho }{\rho _{\mathrm{}}}\right)𝑑\rho _{\mathrm{}}\frac{1}{2\kappa }_{^{2d}}\phi ^{\prime \prime }\left(\frac{\rho }{\rho _{\mathrm{}}}\right)\left|\left(\frac{\rho }{\rho _{\mathrm{}}}\right)\right|^2𝑑\rho _{\mathrm{}},$$

which holds for every admissible entropy generator $`\phi `$ (c.f. definition 4.1), every function $`\rho L^1(^{2d})`$, with $`M=M_{\mathrm{}}`$ and $`\rho _{\mathrm{}}:=\mathrm{exp}(A(x,\xi ))L_+^1(^{2d})`$, with uniformly convex potential $`A`$, i.e.

$$\frac{^2A}{(x,\xi )^2}\kappa _1\mathrm{II}_{2\mathrm{d}},$$
(4.23)
where $`\kappa _1>0`$ and $`\mathrm{II}_{2\mathrm{d}}`$ denotes the identity matrix in $`^{2d}`$. Using this inequality, we can now estimate (recall that $`w_{\mathrm{}}^\pm =M^\pm w_{\mathrm{}}`$)

$$_{^{2d}}\phi \left(\frac{w^\pm }{w_{\mathrm{}}^\pm }\right)𝑑w_{\mathrm{}}^\pm \frac{1}{2\delta \kappa _1}_{^{2d}}\phi ^{\prime \prime }\left(\frac{w^\pm }{w_{\mathrm{}}^\pm }\right)^{}\left(\frac{w^\pm }{w_{\mathrm{}}^\pm }\right)D\left(\frac{w^\pm }{w_{\mathrm{}}^\pm }\right)dw_{\mathrm{}}^\pm ,$$

where $`\delta >0`$ is the smallest eigenvalue of $`D`$ given by

$$\delta :=\frac{1}{2}\left(D_{pp}+D_{qq}\sqrt{(D_{pp}D_{qq})^2+4D_{pq}^2}\right),$$
(4.24)
since by definition of $`\delta `$, it holds that $`\delta \mathrm{II}_{2\mathrm{d}}D`$. Provided lemma 4.1, the results of \[AMTU\] imply the exponential decay of the relative entropy with a rate $`\kappa \delta \kappa _1`$ and thus (4.20) is proved.

Using the well known *Csiszár-Kullback inequality* \[AMTU\], \[Cs\], \[Ku\], we further obtain

$$w^\pm (t,)M_0^\pm w_{\mathrm{}}_1Ke^{\kappa t},K^+,t>0.$$

Since by (A2) we have $`1=M_0^+M_0^{}`$, this implies

$`w(t,)w_{\mathrm{}}_1`$ $`=w^+(t,)w^{}(t,)(M_0^+M_0^{})w_{\mathrm{}}_1`$
$`w^+(t,)M_0^+w_{\mathrm{}}_1+w^{}(t,)M_0^{}w_{\mathrm{}}_1`$
$`Ce^{\kappa t},C^+.`$
For entropy generators of the form (4.14), we get from (4.20)

$$w^\pm (t,)M_0^\pm w_{\mathrm{}}_{2,\delta }e^{\kappa t}w_0^\pm M_0^\pm w_{\mathrm{}}_{2,\delta }t>0,$$

where $`_{2,\delta }`$ is the $`L^2`$-norm with weight $`\delta =1/w_{\mathrm{}}`$, i.e.

$$f_{2,\delta }^2:=_{^d\times ^d}|f(x,\xi )|^2w_{\mathrm{}}^1(dx,d\xi ).$$

As above, the use of the triangle inequality allows us to conclude that exponential decay (with rate $`\kappa `$) of $`w`$ holds in this weighted $`L^2`$-norm. Since $`1/w_{\mathrm{}}1`$ pointwise, we have $`f_2f_{2,\delta }`$ (here we have used the Lindblad condition (A1) to show that the constant in front of $`1/w_{\mathrm{}}`$ is indeed greater than 1).

This implies exponential decay in the usual $`L^2`$-norm and thus, by interpolation, we obtain the exponential convergence of $`w`$ towards the steady state in all $`L^p`$-norms with $`1p2`$. ∎

It should be noted, that the main drawback of the above theorem is the fact that it only holds for initial data with positive and negative part bounded in relative entropy (A4). Furthermore one should have in mind that the above theorem fails, if $`\gamma =0`$ or $`\omega _0=0`$, since then $`A`$ is no longer uniformly convex. However we have already seen that in these cases no finite-mass steady state, different from zero, exists.

The precise value of the rate $`\kappa `$ can be obtained by the following transformation. Define a new potential $`\stackrel{~}{A}`$ by

$$\stackrel{~}{A}(\stackrel{~}{x},\stackrel{~}{\xi }):=A\left(\sqrt{D}(x,\xi )^{}\right),$$
(4.25)
where $`\sqrt{D}`$ is the square root of $`D`$ in the sense of positive definite matrices. Then, it holds that

$$\frac{^2\stackrel{~}{A}}{(\stackrel{~}{x},\stackrel{~}{\xi })^2}\kappa \mathrm{I},$$
(4.26)
which gives the optimal rate $`\kappa `$ as the smallest eigenvalue of the Hessian of $`\stackrel{~}{A}`$. For details, see \[AMTU\] or \[Ri\].

Remark. In contrast to our result, the convergence to the equilibrium state for classical kinetic FP equations (2.20) is, by similar entropy methods, obtained with a suboptimal rate of order $`O(t^{\mathrm{}})`$, see \[DeVi\]. This is due to the fact, that for the classical FP equation, in addition to the thermal equilibrium state (3.31), there exist other states with zero entropy dissipation. Indeed explicit calculations show that $`\kappa `$ depends on $`\overline{}h`$ and converges to $`0`$ as $`\overline{}h0`$.

We hope to extend the above results to more general, maybe nonlinear, potentials $`V`$ in forthcoming works. Although quite simple, the presented harmonic oscillator case is expected to be an essential prerequisite in this analysis.

Acknowledgment

The first author thanks A. Arnold and L. Neumann for helpful discussions.

This work has been supported by the FWF Wissenschaftskolleg ”*Differentialgleichungen*”, the *”Wittgenstein Award 2000”* of P.M. and the FWF-Project P$`14876`$, ”*Fokker-Planck und mittlere Feldgleichungen*”. Additional financial support has been obtained from the EU-TMR project, ”*Asymptotic Methods in Kinetic Theory*”, the Spanish DGES project PB98-1281 and the austrian-french program Amadeus $`02580`$ QF on ”*Partial differential equations for modelling semi-conductors.*