Surface diffusion within the Caldeira–Leggett formalism

Abstract

Surface diffusion is described in terms of the intermediate scattering function in the time domain and reciprocal space. Two extreme time regimes are analyzed, ballistic (very short times) and Brownian or diffusive (very long times). This open dynamics is studied from the master equation for the reduced density matrix within the Caldeira–Leggett formalism. Several characteristic magnitudes in this decoherence process such as the coherence length, ensemble width and purity of the density matrix are analyzed. Furthermore, for flat surfaces, the surface diffusion is considered for the Schrödinger cat states and identical adsorbates or adparticles, bosons and fermions. The analytical results are compared with those issued from solving the Lindblad master equation through the stochastic wave function method. This numerical analysis is extended to be applied to corrugated surfaces.

---

1 Introduction

In 1954 van Hove¹ showed that the scattering cross section of probe particles (low energy neutrons) by a system of interacting particles can be expressed, within the Born approximation, in terms of the so-called generalized pair-distribution function or van Hove space-time correlation function G(r, t), depending on the space vector r and time t. Under this approximation, the scattering problem is reduced essentially to a problem of statistical mechanics where the nature of the scattered particles (neutron, light, atoms, etc.) and details of the interaction potential with the interacting system are irrelevant.² This G-function, which is a straightforward extension of the well known static pair or radial distribution function g(r) widely used in liquid theory,³ gives us the probability that given a particle at the origin and at time t = 0, any particle including the same one is to be found at the position r and at time t. A natural splitting of this G-function can be written as

G(r, t) = G_s(r, t) + G_d(r, t), (1)

where s refers to self and d to distinct, describing the correlations between positions of one and the same particle and pairs of different particles with time, respectively. The well-known properties of these distribution functions are^1,2,4where N is the number of particles in the interacting system and δ(r) is Dirac's function in space.

In this type of scattering, the sample response is linear and it is given by the function S(Δk, ω) which is proportional to the scattering cross section. This response function is often called the scattering law or dynamic structure factor (DSF) depending on the momentum transfer (or scattering vector) Δk = k_f − k_i and energy transfer ħω = E_f − E_i. The corresponding inverse Fourier transform in time, which is also a response function, is known as the intermediate scattering function (ISF), I(Δk, t). Finally, the double inverse Fourier transform (in momentum and energy) of S(Δk, ω) gives us the G(r, t)-function. Thus, one can write

where 〈·〉 stands for thermal (Boltzmann) averages or quantum expectation values depending on if we are considering classical or quantum scattering and r_l, r_j are the position operators of particle l and j, respectively, which do not commute at different times. The functions G_s(r, t) and G_d(r, t) give contributions of the same (l = j) and distinct (l ≠ j) particles, respectively. Due to the linearity of the Fourier transform, one can writeThe so-called coherent and incoherent scattering cross sections in neutron scattering are related to S(Δk, ω) and S_s(Δk, ω) and, therefore, from eqn (3), to G(r, t) and G_s(r, t), respectively. In classical mechanics, these G-functions are real and r_l(0) and r_j(t) represent the initial location of particle l and the classical trajectory of particle j, respectively.

Now, by considering that one particle is at the origin at time t = 0 and simultaneously a second particle is at r′, then the probability that, in an elapsed time t, the second particle goes from r′ to r with a net displacement given by r′ − r is denoted by P₀(r, r′, t). The so-called convolution approximation due to Vineyard⁴ is written as

where it is often assumed that P₀(r, r′, t) ≃ G_s(r − r′, t). This approximation is expected to be better at long times and distances where the correlation between two particles is faded out. Vineyard also derived eqn (5) from a quantum scenario. Within this approximation, one can obtain G_d asin order to have more or less good approximations to G(r, t), I(Δk, t) and S(Δk, ω). Thus, it is then primordial to know G_d to characterize the full scattering process.

In diffusion of adsorbates/adparticles on surfaces sampled by He atoms instead of neutrons, the scattering is fully coherent⁵, that is, only I(Δk, t) and S(Δk, ω) are observables; I_s(Δk, t) and S_s(Δk, ω) being relevant only when the surface coverage is really very small since correlations between any of two particles are almost negligible. In this context, the momentum transfer Δk as well as adsorbate trajectories R(t) are written in capital letters meaning that they are quantities parallel to the surface. Two well established surface experimental techniques are used, the quasi-elastic He atom scattering (QHAS)^6,7 and neutron scattering (QENS)⁸ which overlap in spatial and time resolution. The first one goes more or less from 10⁻¹³ to 10⁻⁸ seconds and from 10⁻¹¹ to 10⁻⁸ m, whereas the second one goes from 10⁻¹³ to 10⁻⁶ seconds and from 10⁻¹⁰ to 10⁻⁶ m.^9,10 QENS is, in general, more convenient for processes occurring in bulk and QHAS is essentially sensible for surfaces. More recently, these two techniques have been complemented by using spin-echo (SE), HeSE^9,10 and neutron spin-echo NSE.⁸ In any case, the experimental results issued from any scattering technique mentioned above require an adequate theoretical framework to properly process them and extract relevant information about the physical systems of interest. For neutrons, molecular dynamics calculations are generally used where a full description of the force fields (adsorbate–adsorbate and adsorbate–substrate interactions) involved is necessary whereas, for He atoms, the Langevin formalism or its generalization to include memory effects is widely applied starting from the well-known Caldeira–Leggett Hamiltonian.¹¹ In the second approach, the surface is usually well represented by a thermal bath consisting of an infinite number of harmonic oscillators. Friction and noise (white or color) appear after integrating over the degrees of freedom of the surface.

After eqn (3), the ISF can be calculated from the classical and quantum stochastic trajectories represented by R(t) for each adsorbate. Classical stochastic mechanics is usually sufficient except for very light adsorbates. One of the main goals is to provide simple analytical expressions to directly interpret the experimental results, shedding light on the interaction of adparticles with the surface and extracting friction and diffusion coefficients. As mentioned above, the corresponding observables or response functions are resolved in time (ISF) or energy (DSF). For such a goal, the dynamics of only one adsorbate is studied (assuming the surface coverage is small). Within the so-called Gaussian approximation, simple analytical expressions for I_s and S_s are easily obtained for flat and corrugated surfaces¹² or with memory friction.¹³ Numerical results from Langevin calculations are thus better interpreted. One way to also provide analytical expressions when the coverage is important and therefore the distinct contribution can not be neglected is through a model called interacting single adsorbate (ISA).^14,15 In the ISA model, the interaction among adsorbates is replaced by a shot noise by assuming that the adsorbates are also a thermal bath (the so-called two bath model). It is clear that when we are faced to strong correlated motions due to high surface coverage, the ISA model is no longer valid. For example, this is the case for sodium atoms diffusing on a Cu(111) surface.¹⁶ However, within the range of applicability of this approximation, one is able to provide analytical expressions for I and S, overcoming the explicit knowledge of I_d or S_d. Two extreme time regimes are well characterized in this open dynamics, the ballistic regime, at very short times, where the dynamics is friction and noise free and the Brownian or diffusive regime, at very long times, where the thermal equilibrium is reached with the surface. Analytical expressions for I_s and S_s are given by means of very simple functions. In the first regime both response functions are Gaussian functions and, in the second regime, I_s is an exponential function and S_s a Lorentzian function.

On the other hand, for very light adparticles, we have to resort to the quantum realm. The so-called quantum Brownian motion is the paradigm in this context. Here position operators do not commute at different times and some care needs to be taken into account. However, a quite straight forward extension of the methodology used for classical diffusion can be still used to provide reasonable analytical expressions.¹⁷ A second approach is to consider quantum stochastic trajectories within the Bohmian framework.¹⁸ An alternative and widely used approach to deal with open quantum systems is within the so-called system-plus-environment approach. A master equation is then derived for the evolution of this matrix which contains both frictional and thermal effects, the so-called Caldeira–Leggett (CL) master equation¹⁹ which is of Markovian type. The corresponding diagonal matrix elements give the quantum probabilities and the off-diagonal elements, the so-called coherences. Time evolution of coherences gives us an indication of how the decoherence process is gradually established leading to certain timescales of the system under study and exponential suppression of spatial interference terms.

Strictly speaking the splitting of the G-function into two parts (self and distinct) is only valid in classical statistical mechanics. Several quantum scenarios can be devised for this failure due to the presence of interference or cross terms which can drastically modify the probability density or, equivalently, the G-function: the so-called Schrödinger cat states for adsorbates which can be at both sides of two domains or shallow wells on a surface, and when several distinguishable or identical (interacting or noninteracting) adsorbates are considered in the diffusion dynamics. In this last case, the symmetry of the wave function (symmetric for bosons and anti-symmetric for fermions) adds a new ingredient to this open dynamics. In all of these example, one should only consider the G function, not anymore G_s and G_d separately. Interestingly enough, we have observed that in the ballistic regime the ISF and DSF no longer display a Gaussian function. However, in the Brownian or diffusion regime these two response functions still conserve their analytical shapes. In the interest of simplicity, this theoretical analysis is carried out for flat surfaces where an analytical expressions can be found by solving the CL master equation. On the other hand, a numerical technique has been developed to solve the stochastic Schrödinger differential equation for corrugated surfaces in order to analyze these two time regimes for low energy barriers at high temperatures. The analytical results are then compared with those issuing from solving the so-called Lindblad master equation^20,21 through the stochastic wave function method.^22,23 Numerical results will also be reported for model calculations on the Xe–Pt(111) system in order to check the application of this method in the surface diffusion context.

In brief, a new theoretical approach for the description of measurements of atomic diffusion by particle beams such as realized in helium atom scattering and neutron scattering is presented. This method should represent an important step on the development of quasi-elastic helium scattering and neutron spectroscopy. Due to the fact it can be applicable to a wide range of realistic systems, a high impact on a wide range of physicochemical research is expected.

This paper is organized as follows. In Section 2 the CL master equation in the coordinate representation is briefly introduced and applied in this context for the first time. Several characteristic magnitudes are then studied governing the coherence of the diffusion dynamics such as coherence length, ensemble width, purity and coherence time are presented and analyzed. In Sections 3 and 4, the Schrödinger cat state model and the diffusion of noninteracting bosons and fermions are discussed, respectively. In Sections 5 and 6, the Lindblad master equation through the stochastic wave function numerical analysis is carried out in order to compare with previous analytical results for flat and corrugated surfaces. Finally, in the last section, some conclusions and future perspective are discussed.

2 The Caldeira–Leggett formalism. The one particle problem

Many important problems in physics and chemistry are dealing with open quantum system.^24–26 A complete mathematical model for the system-plus-environment dynamics can be very complicated. The environment is usually well represented by a reservoir with an infinite number of degrees of freedom.^22,27 It turns out to be useful to formulate this dynamics using an appropriate equation of motion for the corresponding density matrix, known as quantum master equation, and employing some approximations.

Several techniques have been used to describe the dynamics of an exactly open quantum system in terms of the reduced density matrix. One of this techniques is the influence functional formalism used in the path integral representation where the environmental variables can be eliminated by using the Feynman–Vernon influence functional.²² This formalism is an efficient theoretical tool and has been applied for solving Langevin equation with damped harmonic oscillator,²⁸ non-Markovian systems,²⁹ and the spin-boson model^30,31 among others.

We will focus on discussing the simplest case, the motion of a quantum Brownian particle³² in the weak-coupling and high-temperature limits. Thus, for example, by assuming short environmental correlation times, the memory effects can be safely neglected in the reduced dynamics and a Markovian quantum master equation can be easily derived which is known as the CL master equation.³³ This formalism gives the evolution of the reduced density matrix in the coordinate representation by tracing out the environment degrees of freedom and contains both frictional and thermal effects due to the environment in analogy with the classical Fokker-Planck equation.²²

This pioneering CL work describes the dynamics of such a particle linearly coupled to an Ohmic environment³⁴ (linear dissipation). The CL master equation for the reduced density matrix in the position representation for one dimension ρ(x, x′, t), at the high temperature limit, takes the form³³

where V is the external interaction potential, m the particle mass, γ the friction coefficient and D the diffusion coefficient given by

D = 2mγk_BT, (8)

with k_B being Boltzmanns constant and T the surface temperature. This master eqn (7) for one particle coupled to an Ohmic environment (constant friction) with a generic potential does not have an analytic solution. However, for flat surfaces (V = 0), analytical solutions are obtained by considering, for example, a Gaussian ansatz³⁵ or using the Wigner representation.³⁶

Let us consider the initial state as a Gaussian wave packet representing a localized particle centered at the position x₀, with momentum p₀ and width σ₀,

The initial density matrix, given by

ρ(x, x′, 0) = ψ*(x′, 0)ψ(x, 0), (10)

can be extracted from the initial state (9). Replacing the initial wave packet (9) in (10), the density matrix can be expressed asUsing the Gaussian ansatzand inserting eqn (12) into the master eqn (7), a system of coupled differential equations for the time dependent coefficients A(t), B(t), C(t), D(t), E(t) and N(t) is obtained and shown in Appendix A. The diagonal elements give the probability to find the particle at position x and time t whereas the nondiagonal terms give the so-called coherences which provide the correlation of the particle at two different positions x and x′ at time t. With time, this correlation or coherence decays obviously to zero. This solution can also be obtained by using the Wigner representation of the reduced density matrix, see Appendix B.

The generalized pair-distribution function of van Hove G(x, t) can be associated with the reduced density matrix given by eqn (12) and by imposing the condition x′ = x leading to

where

x_t = x₀ + v₀f(t), (14)

Thus, G(x, t) is a Gaussian function with a width given by eqn (16) and a center moving along the classical trajectory (14). For one particle, the general properties of G or G_s given eqn (2) are fulfilled. As known,^37–39 the variance of a Gaussian wave packet in general is given byIn correspondence with eqn (17), the time dependent width given by eqn (16) has three contributions: (i) the initial width, (ii) the commutator through −[x (0), x(t)]²/4σ₀² and (iii) the mean square displacement (MSD). The expression for G(x, t) provides an analytical solution for the probability of finding the particle at position x at time t given that it was at x₀ at t = 0.

Once G(x, t) has been identified, the corresponding ISF, I(Δk, t), and the dynamic structure factor, S(Δk, ω), can be derived from eqn (3). Thus, one has that

which is Gaussian in the momentum transfer which is represented along this work as Δk. In diffusion processes, two well established time regimes are characterized, ballistic and Brownian or diffusive. In the ballistic regime, at short times, t ≪ γ⁻¹, and by considering the linear term of the exponential, eqn (14) and (16) reduce toand the corresponding ISF is given bythat is, a pure Gaussian function in time and momentum transfer is obtained for a zero initial position. The time exponential factor depends on the particle mass and the initial width. In the classical regime, this factor depends on the thermal velocity of the particle.¹² On the other hand, in the Brownian regime, we have just the opposite case, t ≫ γ⁻¹, andand the ISF isAs expected, in the Brownian regime, I(Δk, t) is an exponential function with time and is accompanied by the diffusion coefficient D which contains itself the surface temperature T and the friction coefficient γ according to eqn (8). Thus, the ISF gradually passes from a Gaussian to an exponential function when going from the short time to the long time regime.

In Fig. 1, the real components of the intermediate scattering function, eqn (18), are plotted for Xe (blue solid curve) and Na (red dashed curve) atoms. The parameters for the initial wave packet are: x₀ = 0, σ₀ = 0.003 a.u. for Xe, σ₀ = 0.02 a.u. for Na, p₀ = 0. The friction γ = 0.25 ps⁻¹, momentum transfer Δk = 0.12 Å⁻¹ and surface temperature T = 105 K have been extracted from previous works with this adsorbates.^40,41 For the Xe adsorbate, both the Gaussian function at short times and the exponential at long times give good fitting to the ISF. This result indicates that the adsorbate goes from a ballistic to diffusive regime. The corresponding DSF S(Δk, ω) is the time Fourier transform of eqn (18) and can not be obtained analytically due to the presence of the time dependent width. However, in both limiting regimes, analytical expressions are easily obtained. The dynamical structure factor in the ballistic regime is expressed as

and in the Brownian or diffusive regimeFrom eqn (23) and (24), it can be seen that the dynamical structure factor gradually passes from a Gaussian to a Lorentzian function with the frequency when going from the ballistic to Brownian regime. This gradual transition is the so-called motional narrowing effect.⁴²

Fig. 1 The real components of the intermediate scattering function given by eqn (18) are plotted for Xe (blue solid curve) and Na (red dashed curve) atoms. The parameters for the initial wave packet are: x₀ = 0, σ₀ = 0.003 a.u. for Xe, σ₀ = 0.02 a.u. for Na, and p₀ = 0. The fitting of a Gaussian function (black line) for short times and an exponential function (purple line) for long times are shown for the Xe adsorbate. It is assumed that γ = 0.25 ps⁻¹, Δk = 0.12 Å⁻¹ and T = 105 K.

2.1 Characteristic magnitudes in decoherence dynamics

The dynamics of adsorbates on a surface seen as a thermal bath leads gradually to decoherence or, in other words, to the quantum-classical transition. Several magnitudes are used to characterize this decoherence process in the position representation of the reduced density matrix, namely: the coherence length l(t), the ensemble width Δx(t) and the so-called purity δ(t).³⁵

One of the most important quantities in decoherence is the so-called coherence length l(t). This magnitude represents a measure of the distance over which quantum correlations are important (the width of ρ in the x–x′ direction, the nondiagonal matrix elements). In correspondence with,⁴³ the coherence length l(t) can be obtained from the Gaussian ansatz (12) and the density matrix (11),

It also measures the characteristic distance over which the system can exhibit spatial interference effects. The explicit time dependence is found in Appendix C.

In Fig. 2, the time dependence of the coherence length for three adsorbates with different masses is plotted: m_Xe = 131.293, m_Na = 22.989 and m_Li = 6.941. These calculations are carried out for an initial Gaussian wave packet (9) with the following initial conditions: x₀ = 0, σ₀ = 0.003(Xe), 0.02(Na) and 0.2(Li) a.u. and momentum p₀ = 0. The surface temperature is T = 105 K and the friction coefficient γ = 0.25 ps⁻¹. The color curves are under the presence of the environment (solid lines) with Xe (blue), Na (red) and Li (purple) whereas the dashed curves give the free evolution of Xe, Na and Li. At very short times (t ≪ 4 ps), l(t) displays a very pronounced linear increase identical to the free evolution but suddenly the presence of the environment starts reducing significantly the coherence length reaching an asymptotic value around 4 ps. The decoherence process is gradual established but it is going faster for heavier adsorbates according to eqn (120) from Appendix C.

Fig. 2 Time dependence of the coherence length l(t) for three adsorbates Xe, Li and Na. The initial parameters for the Gaussian wave packet (9) are: position x₀ = 0, width σ₀ = 0.003(Xe), 0.02(Na) and 0.2(Li) a.u., and momentum p₀ = 0, a temperature T = 105 K and a friction coefficient γ = 0.25 ps⁻¹. With environment (solid lines) and without (dashed line): Xe (blue), Na (red) and Li (purple).

The second key quantity is the so-called ensemble width, Δx(t), which represents the probability distribution size P(x, t) ≡ ρ(x, x, t) and can be found from the following expression

which coincides with the width of the probability distribution.

Fig. 3 shows the time dependence of the ensemble width, given by eqn (26), for the above three adsorbates Xe, Na and Li by assuming again the same initial conditions and environment parameters. The same color codes are used as in Fig. 2. At very short times (t ≪ 4 ps), there is no difference among the ensemble width with and without environment. This means that in the ballistic regime the adsorbates moves freely in agreement with the coherently spread out of a free Gaussian wave packet (126) from Appendix C. However, at long times (t ≫ 4 ps), the ensemble width of the three adsorbates increases more slowly and much more for heavier masses. This behavior indicates that the adsorbate can diffuse less freely as a result of the interaction with environment.

Fig. 3 Time dependence of the ensemble width Δx(t) for the three adsorbates Xe, Li and Na. The same initial conditions for the wave packet, surface temperature and friction are used as in Fig. 2 The color codes for the curves are also used in Fig. 2.

Finally, the third key quantity is a dimensionless measure of decoherence given by the ratio of the coherence length l(t) and the ensemble width Δx(t) called the purity of the reduced density matrix,

The purity of a normalized quantum state is also defined as δ(t) = Tr(ρ²). This magnitude gives information on how much a state is mixed. If the quantum state ρ(x, x′, t) represents a pure-state density matrix, ρ² = ρ and the purity is δ(t) = 1.

The explicit expression is again found in Appendix C.

In Fig. 4, the time dependence of the purity δ(t) is also analyzed for the same adsorbates and initial conditions. Without environment, the reduced density matrix purity remains constant equal to one since the free ensemble size eqn (126) and the free coherence length eqn (118) are equal. Under the presence of the environment, the purity decreases faster for Xe than for Na and Li. This decreasing occurs at very short times (t ≪ 4 ps). In other words, the gradual decoherence process is very fast for the initial conditions chosen in these examples.

Fig. 4 Time dependence of the purity δ(t) for the three adsorbates Xe, Li and Na. The same initial conditions for the wave packet, surface temperature and friction are used as in Fig. 2 The color codes for the curves are also used in Fig. 2 and without environment for the three adsorbates (black solid curve).

3 The Schrödinger cat states

Sometimes the adsorbate can be at two different positions/domains at the same time; for example, if there is a barrier and the tunneling process takes place, the adsorbate has a probability different from zero to be in both side at the same time.

Let us consider now, for a single adsorbate, an initial wave function expressed as the superposition of two equal weighted Gaussian wave packets (9) centered around x = ±x₀ at t = 0 and the same initial momentum p₀

where [scr N, script letter N] is the normalization constant. From eqn (28), the corresponding initial reduced density matrix is given byconsisting of four contributions: two terms associated with each wave packet separately and the interference terms. Due to the linearity of the master eqn (7), the time evolution of the reduced density matrix when V = 0 leads towhere the interference term is written asUsing again a Gaussian ansatzthe explicit expression for the coefficients can be found in Appendix D. The probability of finding the adsorbate at (x, t) can be extracted from the reduced density matrix eqn (30) by imposing the condition x = x′, that is,whereand analogously for P_−x0(x, t). Using the property of the hermiticity of the reduced density matrix, the interference term at x and at time t is given then bywhere is the modulus of ρ_interf(x, x, t) and Θ(x, t) its phase. Eqn (33) can be expressed as the typical interference pattern expression according to⁴⁴where Γ(t) is known as the decoherence function written asAfter some straightforward calculations, one obtains the explicit expressions for the normalization constantthe phasewithand the decoherence functionAs can be clearly seen from eqn (41), the decoherence function is negative since the width σ_t, given by eqn (16), is always greater than the initial width σ₀. Furthermore, when D = 0 or T = 0, the width σ_t isand the decoherence function vanishes. This is an expected behavior since the decoherence process⁴⁵ is represented by the last term of the CL master eqn (7).

In the zero dissipation limit, γ → 0, the decoherence function is expressed as

and when σ₀ ≪ 2x₀, the decoherence function can be written aswhere d = 2x₀ is the separation between initial wave packets and τ_D is known as the decoherence time. From eqn (44) one sees again that τ_D depends on the friction coefficient, temperature and the initial wave packets separation.⁴⁶

In Fig. 5, the decoherence function is plotted for Li adsorbates with different friction coefficients (left panel) and surface temperatures (right panel). The color codes are inside each plot. As expected, the decoherence function which is a negative function decreases with friction and temperature.

Fig. 5 Decoherence function for Li adsorbates with different friction coefficients (left panel) and surface temperatures (right panel) for an initial position x₀ = ±2 a.u., width σ₀ = 0.2 a.u. and momentum p₀ = 0. On the left panel, the friction coefficients are 0.02 ps⁻¹ (blue solid curve), 0.25 ps⁻¹ (red dashed curve) and 2.50 ps⁻¹ (purple dotted curve). On the right panel, the surface temperatures are 50 K (blue solid curve), 105 K (red dashed curve) and 200 K (purple dotted curve).

The ISF is then given by

I(Δk,t) = [scr N, script letter N]²(I(Δk,t) + I_inter(Δk,t)), (45)

whereIn the ballistic regime, t ≪ γ⁻¹, eqn (40) reduces toand the corresponding ISF isAs can be seen, this function is not longer a Gaussian function.

However, in the Brownian regime, t ≫ γ⁻¹, one has

and the ISF iswhere the same exponential function of time is still valid.

The Schrödinger cat state or the quantum superposition states become more likely for H particles than for Na or Xe. The ratio between the decoherence time τ_D and relaxation time τ_r for a superposition of two different positions at a distance Δx is τ_r/τ_D = (Δx/λ_dB)², where the thermal de Broglie wavelength is given by . With increasing mass, the decoherence time scale is much smaller than the relaxation time scale required to reach thermal equilibrium. It becomes more difficult to observe interference patterns for massive particles and therefore the creation of a spatial superposition.^47,48

In Fig. 6, the ISF given by eqn (45) is plotted (blue solid line) for an H adsorbate with width σ₀ = 0.2 a.u., positions x₀ = ±0.6 a.u. and momentum p̄₀ = −p₀ = 0.01 a.u., a friction coefficient γ = 0.25 ps⁻¹, a momentum transfer Δk = 0.12 Å⁻¹ and two values of temperature: T = 10 K (left panel) and T = 100 K (right panel). The insets show the interference term (47) together with the ISF (46) with a (dashed red line) and a (dotted purple line), respectively. The interference term breaks the ballistic behavior of the system. As expected, when the temperature increases, the ISF decreases faster leading to the typical exponential behavior with time.

Fig. 6 Time dependency of the real part of the intermediate scattering function for a H particle with the initial wave packet (28) with width σ₀ = 0.2 a.u., positions x₀ = ±0.6 a.u. and momentum p̄₀ = −p₀ = 0.01 a.u. for a friction coefficient γ = 0.25 ps⁻¹, a momentum transference Δk = 0.12 Å⁻¹ and two values of temperature: T = 10 K (left panel) and T = 100 K (right panel). The inside panel shows the interference term (47) and the self intermediate scattering function (46) with a (dashed red line) and a (dotted purple line) respectively and the (blue solid line) the intermediate scattering function (45).

Thus, the corresponding DSF S_s(Δk, ω), in the Brownian regime, is again a Lorentzian function with the frequency.

4 Identical adsorbates

For neutral atomic adsorbates, those with an even number of neutrons are bosons whereas those with an odd number of neutrons are fermions.⁴⁹ In this section, we are going to consider the diffusion of noninteracting identical adsorbates where the only interaction comes from the symmetry of the corresponding wave functions with respect to the interchange of particles, symmetric for bosons and anti-symmetric for fermions. In particular, we want to analyze how the well-known bunching and anti-bunching properties of bosons and fermions, respectively, are manifested in this open dynamics; in other words, in what extent the symmetry of the wave function is robust enough to keep it along time and how the gradual loss of being indistinguishable is established.

The initial wave function for the identical adsorbates (symmetric, +, and anti-symmetric, −) can be written as

Ψ_±(x₁,x₂,0) = [scr N, script letter N]_±{ψ(x₁,0)ϕ(x₂,0) ± ϕ(x₁,0)ψ(x₂,0)}, (52)

where ψ and ϕ are one-particle Gaussian wave packets (9) with parameters x₀, p₀, σ₀ and x̄₀, p̄₀, [small sigma, Greek, macron]₀, respectively. Under the Caldeira–Leggett eqn (7), the evolution of the two-particle system is given bywhereNote that ρ₁₁(x, x′, t) and ρ₂₂(x, x′, t) are the density matrix of one adsorbate. As seen before, the probabilities are given by the diagonal elements of the density matrix eqn (53),where P_ij(x, t) = ρ_ij(x, x, t). The last or cross term of eqn (55) is responsible for symmetry effects. For the single-particle density one has thatand the probability is written as

P_sp,±(x, t) = [scr N, script letter N]_±²{P₁₁(x, t) + P₂₂(x, t) ± 2Re[P₁₂(x, t)s(t)]}, (57)

where P₁₁(x,t) and P₂₂(x,t) are given by eqn (34), the overlapping integral s(t) beingand the normalization factor [scr N, script letter N]_±From eqn (58) it can be seen that s(t) is time independent and does not depend on environment parameters such as γ and T, only depends on the initial conditions. In analogy to the cat state problem, eqn (35), one can write

P_sp,±(x, t) = [scr N, script letter N]_±²{P₁₁(x, t) + P₂₂(x, t) ± 2|P₁₂(x, t)s(t)|cos Θ(x, t)}, (60)

whereand Γ₁₂(t) is the decoherence function for noninteracting identical adsorbate, given byThe probability P₁₂(x, t) is now written aswhere

andNote that for [small sigma, Greek, macron]₀ = σ₀, the width is b₂(t) = σ_t/2 and the explicit formula for the phase Θ(x, t) iswithFrom eqn (63), the interference term of the probability (60) can be written aswhereWhen x̄₀ = x₀, the probabilities P₁₁(x, t) and P₂₂(x, t) are given by (34) for the momentum p₀ and p̄₀, respectively. In this case, from eqn (62), one obtains the decoherence functionand the phasein correspondence with ref. 50.

Now, if we identify

G(x, t) ≡ P_sp,±(x, t) (71)

the ISF is found to be

I_sp,±(Δk, t) = [scr N, script letter N]_±²(I_x0(Δk, t) + I_x̄0(Δk, t) ± I_c(Δk, t)), (72)

whereand the cross term

In contrast to distinguishable adsorbates where the intermediate scattering function is given by I_x0(Δk, t) + I_x̄0(Δk, t), for noninteracting identical adsorbates the cross term (74) takes into account the symmetry of the wave function. From eqn (72), we have that for the ballistic regime I_sp,± is no longer a Gaussian function whereas for the Brownian regime the exponential behavior with time is still conserved.

In Fig. 7 the initial probability eqn (60) is plotted for distinguishable adsorbates (solid black curve), bosons (blue dashed curve) and fermions (red dotted curve) and for two different initial conditions of the wave packets: (i) different initial positions x₀ = 0 and x̄₀ = 14 a.u. with the same initial momenta p₀ = p̄₀ = 2 a.u. (left panel), and (ii) different initial momenta p₀ = 0.2 a.u. and p̄₀ = 0.3 a.u. with the same initial positions x₀ = x̄₀ = 0 a.u. (right panel). In both cases, the initial widths of the Gaussian functions are the same, σ₀ = [small sigma, Greek, macron]₀ = 2 a.u. The bunching and anti-bunching property of the bosons and fermions are already seen at t = 0, respectively. The distribution is a little bit narrower for the former and the distinguishable distribution is in-between. For fermions, the bimodal character of the initial distribution is kept for both sets of initial conditions. The same behavior is reported elsewhere.^51,52

Fig. 7 Initial spatial probability for a light adsorbate considered to be distinguishable (black solid curve), boson (blue dashed curve) and fermion (red dotted curve). Two different initial conditions of the wave packets are considered: (i) different initial positions x₀ = 0 and x̄₀ = 14 a.u. with the same initial momenta p₀ = p̄₀ = 2 a.u. (left panel), and (ii) different initial momenta p₀ = 0.2 a.u. and p̄₀ = 0.3 a.u. with the same initial positions x₀ = x̄₀ = 0 a.u. (right panel). In both cases, the initial widths of the Gaussian functions are the same, σ₀ = [small sigma, Greek, macron]₀ = 5.2 a.u.

In Fig. 8, the ISF for a light adsorbate considered to be distinguishable (black solid curves), boson (blue dashed curves) and fermion (red dotted curves) with a friction coefficient 0.25 ps⁻¹, momentum transfer Δk = 0.12 Å⁻¹ and surface temperature T = 105 K. In the four panels several initial conditions are considered: (a) x₀ = 0, x̄₀ = 14 a.u., p₀ = p̄₀ = 0.2 a.u., and σ₀ = [small sigma, Greek, macron]₀ = 5.2 a.u.; (b) same as (a) but σ₀ = [small sigma, Greek, macron]₀ = 4.8 a.u.; (c) x₀ = x̄₀ = 0, p₀ = 0.2 a.u., p̄₀ = 0.3 a.u. and σ₀ = [small sigma, Greek, macron]₀ = 5.2 a.u. and (d) same as (c) but σ₀ = [small sigma, Greek, macron]₀ = 4.8 a.u. Several interesting features are worth discussing. First, it is clear that unlike the black curve (distinguishable adsorbates), the ISF for fermions, for the ballistic regime, is no longer a Gaussian function. Second, in the Brownian or diffusion regime, the symmetry of the wave function is not robust enough to distinguish bosons and fermions since adsorbates tend to the same asymptotic exponential function. Third, the initial width makes the difference between the symmetry of the wave function more pronounced in the intermediate time region. And fourth, it should be possible to observe these differences in an experiment where the fermions decay slower than bosons in panels (c) and (d) (and bosons decay quite similar to distinguishable adsorbates); in other words, the decoherence process is much more gradual for fermions than bosons for certain initial widths.

Fig. 8 Intermediate scattering function for a light adsorbate considered to be distinguishable (black solid curves), boson (blue dashed curves) and fermion (red dotted curves) with a friction coefficient 0.25 ps⁻¹, momentum transfer Δk = 0.12 Å⁻¹ and surface temperature T = 105 K. Panels: (a) x₀ = 0, x̄₀ = 14 a.u., p₀ = p̄₀ = 0.2 a.u., and σ₀ = [small sigma, Greek, macron]₀ = 5.2 a.u.; (b) same as (a) but σ₀ = [small sigma, Greek, macron]₀ = 4.8 a.u.; (c) x₀ = x̄₀ = 0, p₀ = 0.2 a.u., p̄₀ = 0.3 a.u. and σ₀ = [small sigma, Greek, macron]₀ = 5.2 a.u. and (d) same as (c) but σ₀ = [small sigma, Greek, macron]₀ = 4.8 a.u.

5 The Îto stochastic differential equation

The dynamics of an open quantum systems so far discussed has been based on the CL or reduced density matrix formalism which is valid for weak coupling and at high temperatures. In the density matrix formalism, the Liouville functional [script L] governing the dissipative evolution of the reduced density matrix ρwhere Ĥ is the Hamiltonian of the system.

As demonstrated in a number of works,^53,54 the Caldeira–Leggett master eqn (7) does not preserve the positivity of ρ(t). Derivation of master equations from microscopic Hamiltonians preserving the positivity of the density matrix is a fundamental problem in previous works.^20,55

The Lindblad theory of quantum dynamical semigroups²⁰ together with Kossakowski et al.²¹ showed that the generator for a completely positive map should be

where Â_k are known as Lindblad dissipation operators. The structure of these operators Â_k are unknown and does not generally assure equilibrium between the system and the bath.

One might ask the question of whether the CL eqn (7) can be brought into Lindblad form. The answer is positive since it can be written in Lindblad form by just adding an additional term which is small in the high temperature limit.

Rewriting CL master eqn (7) in the operator form,

employing the minimal invasive modification, an additional term of the form γ[p̂, [p̂, ρ]]/8mk_BT can be added to the master eqn (77)

One can associate the CL equation with the Lindblad form if one Lindblad operator  given by a linear combination of the position x̂ and momentum p̂ operators is taken

 = μx̂ + iνp̂, Â⁺ = μx̂ − iνp̂, (79)

where the coefficients μ and ν are parameters to be determined below. Replacing the generator (76) and the Lindblad operator (79) in the Liouville eqn (75), a Markovian master equation is obtainedIn the coordinate space, the master eqn (80) takes the following formwhereWhile the master equation can be directly solved in a double-space (two-dimensional) representation (81), it can also be solved by using a set of stochastic wave functions {|ψ〉}, where each function obeys the following differential equationwhere x̄ = 〈ψ|x̂|ψ〉, p̄ = 〈ψ|p̂|ψ〉 and dξ are complex stochastic variables fulfilling the properties of a Wiener process: , and M[dξ_kdξ_n] = 0⁵⁶ (see Appendix E).

This formulation has been used in previous works^23,55 for describing the dynamics of adsorbates on metal surfaces with some constraints.

6 Numerical results

The stochastic wave functions ψ_i(x, t) are obtained by solving numerically the Îto stochastic differential eqn (83). By using the variant described in ref. 57, the split operator method has been used to solve (83). The terms with only x̂ or p̂ are propagated in the coordinate or momentum space, respectively. In correspondence with previous works,²² the temperature has been incorporated in the system dynamics through the Lindblad operator coefficients and .

The initial wave function is a Gaussian wave packet is chosen to be

where the initial position is chosen to be x₀ = 0, initial width σ₀ = 0.1 a.u., momentum p₀ is distributed according to the Maxwell–Boltzmann distribution for a given temperature T and an adsorbate with mass m.

The numerical stochastic wave function ψ (x, t) is evaluated in equally spaced points belonging to the interval [x₁, x_Ns], with N_s being a power of 2 in order to use the fast Fourier transform and for a number of realizations N. At each time t, the discrete wave function is represented by

ψ(x, t) ≡ {ψ(x₁, t);ψ(x₂, t);ψ(x₃, t);…;ψ(x_Ns, t)}. (85)

The discrete wave function normalization is calculated aswhere Δx is the length between two consecutive points of the grid. The normalized stochastic wave function is thenand the mean value of any operator Ĉ can be expressed as

The numerical code was parallelized to solve N times the stochastic differential eqn (83). The number of realization is chosen in order to reach the numerical stability for the mean value of each observable. The unidimensional space x is chosen to be in the interval x ∈ [x_i, x_f], where x_i = −120 Å and x_f = 120 Å and the number of points N_s = 4096. The increment Δx can be calculated from Δx = (x_f − x_i)/N_s and by using these parameters, Δx = 0.05 Å.

The parameters that describe the simulation of the system were selected by taking into account the results presented elsewhere.⁵⁸ For this goal, the temporal integration step is Δt = 66.1462 a.u.= 1.6 fs, γ = γ₀/2m with γ₀ = 5 × 10⁻² ps⁻¹ and the temperature T = 121 K.

6.1 Flat surface

It is important to analyze how the width of the wave packet σ_t changes with time. By fitting every 10 ps the stochastic wave function ψ(x, t) to a Gaussian wave packet, exp{−(x − x_t)²/(2σ_t²)}, the width values σ_t in this time interval can be extracted as seen in Table 1. In the first column, the analytical values are reported whereas in the second and third columns the width values are listed for N = 3200 and N = 10 000. This comparison shows the importance of the statistics in this type of simulations. The temperature dependence of the coefficients μ and ν fulfill the condition reported in ref. 55. In Fig. 9–11 the time dependence of the wave packet widths up to 50 ps are displayed for different surface temperatures, frictions and adsorbate masses, respectively. The numerical (points) and analytical (solid curves) results are shown in each plot. In Fig. 9, the friction is given by γ = γ₀/2m where γ₀ = 5 × 10⁻² ps⁻¹ with different values of T and initial widths σ₀ for Xe adsorbates: T = 170 K and σ = 0.08 a.u. (black line), T = 121 K and σ = 0.1 a.u. (red line) and T = 70 K and σ = 0.14 a.u. (blue line).

Table 1 Comparison between the numerical and analytic width for the dynamics of a Xe adsorbate after a time of 50 ps, a Gaussian initial wave packet with width σ₀ = 0.1 a.u., friction γ = γ₀/2m, where γ₀ = 5 × 10⁻² ps⁻¹ and surface temperature T = 121 K

σ 0 = 0.049 Å	σ 0 = 0.049 Å	σ 0 = 0.049 Å
Analytic^50	N = 3200	N = 10 000
0.858	0.858	0.858
0.062	0.061	0.061
0.046	0.045	0.046
0,039	0.039	0.039
0.035	0.034	0.035
0.033	0.031	0.032

---

Fig. 9 Time dependence of the wave packet width for different temperatures and widths σ for Xe adsorbates up to 50 ps with friction γ = γ₀/2m, where γ₀ = 5 × 10⁻² ps⁻¹. The numerical and analytic results are shown by points and solid curves, respectively: T = 170 K and σ = 0.08 a.u. (black line), T = 121 K and σ = 0.1 a.u. (red line) and T = 70 K and σ = 0.14 a.u. (blue line).

Fig. 10 Time dependence of the wave packet width for Xe adsorbates and different values of γ and σ up to 50 ps. The surface temperature is T = 121 K. The numerical and analytic results are shown by points and solid curves: γ = 1 ps⁻¹ and σ = 0.7 a.u. (green line), γ = 5 ps⁻¹ and σ = 0.1 a.u. (blue line) and γ = 50 ps⁻¹ and σ = 0.25 a.u. (red line).

Fig. 11 Time dependence of the wave packet width for three different adsorbates Xe, Li and Na and different values of σ with the same friction coefficient γ = γ₀/2m, where γ₀ = 5 × 10⁻² ps⁻¹ and temperature T = 121 K after a time of 50 ps. The numerical and analytic results are shown by points and solid curves: γ = 50 ps⁻¹ and σ = 0.07 a.u. (red line), γ = 1 ps⁻¹ and σ = 0.08 a.u. (green line) and γ = 5 ps⁻¹ and σ = 0.1 a.u. (blue line).

In Fig. 10, the analysis is carried out for Xe adsorbates and T = 121 K with different values of γ and initial widths. And in Fig. 11 different adsorbates such Xe, Na and Li and different values of σ₀ with the same friction coefficient and surface temperature T = 121 K are plotted. As expected, with temperature, larger value of the wave function widths are obtained. With friction and mass, the corresponding widths reach an asymptotic value at short times. In all of cases studied, the agreement between numerical and analytical results is fairly good indicating that the number of realizations used is good enough to obtain reliable results (in fact, with 2000–4000 realizations the numerical stability is reached).

6.2 Corrugated surface

In order to know how works the stochastic wave function method (83) for more realistic systems, several numerical calculations are compared with previous Langevin simulations and experimental results from Ellis et al.⁴¹

The potential energy surface is obtained from the pairwise additive potential approximation between the Xe and a square network of Pt atoms with a distance equal to the unit cell length a = 3.93 Å.²³ The metal surface is considered to be large enough so that the boundary effects do not influence the interaction potential. For this goal, the number of unit cells N_c of the surface in the XY plane is increased until the variation in the potential is as small as a given threshold, say 10⁻³ meV. A cosine corrugation function is assumed to represent the potential energy surface of the Xe–Pt(111) system given by , where V₁ = −13.64 meV and V₂ = −12.32 meV. These values are obtained by fitting the cosine function to the numerical result of the potential energy surface. Both interaction potentials are very similar because its amplitude is around −28 meV and the periodicity match with the unit cell length of the Pt in both cases. The unidimensional corrugated potential used in the simulations is . The amplitude V₀ is taken so that the energy barrier for this potential (E_b = V₀) coincides with the Langevin molecular simulations.⁴¹ The diffusion of the Xe adsorbate on the metal surface of Pt is studied along the (100) direction.

First, in the ballistic regime, the broadening (Γ) is analyzed in terms of the surface parallel momentum transfer Δk. In this regime, the ISF is a Gaussian function and has been fitted to exp(−t²/2σ_t). As we have seen previously,

For an ideal gas, the width of the Gaussian function is . By using eqn (89), the quasielastic broadening isThe unidimensional space x is chosen to be in the interval x ∈ [−800, 800] Å, the number of points N_s = 8192, the position step Δx = 0.2 Å, the time step Δt = 66.1462 a.u. = 1.6 fs and the surface temperature T = 105 K. The number of realizations is N = 2000. This comparison is shown in Fig. 12 for Xe adsorbates on a Pt(111) surface at two different friction coefficients, γ = 0.25, 2.00 ps⁻¹. In panel (a) the following results are plotted: numerical simulations (NS) for a flat surface with 0.25 ps⁻¹ (blue dots) and 2.00 ps⁻¹ (red dots), Langevin dynamics simulations for 0.25 ps⁻¹ (orange star), 2.00 ps⁻¹ (purple star), theoretical results for an ideal gas (90) (black line) and QHAS experimental results at two incident He beam energies E_i = 10.15 meV (circles) and E_i = 26.85 meV (squares). In panel (b), numerical results are plotted for two different potential energy barriers E_b = 9.6 meV (blue dot) and E_b = 23.6 meV (red dot) with a friction coefficient of 0.05 ps⁻¹. The Langevin simulations are given for orange stars and purple, respectively. Experimental values are the same as before. In general, the overall agreement between our method, Langevin simulations and experimental results is fairly good.

Fig. 12 Γ versus Δk for Xe adsorbates on a Pt(111) surface in the ballistic regime. QHAS experimental results at two incident He beam energies E_i = 10.15 meV (circles) and E_i = 26.85 meV (squares) are shown and the ideal gas (black line) given by (90). In panel (a), the following results at two different friction coefficients are plotted: numerical simulations (NS) with 0.25 ps⁻¹ (blue dots), 2.00 ps⁻¹ (red dots); Langevin simulations (LS) for 0.25 ps⁻¹ (orange star) and 2.00 ps⁻¹ (purple star). In panel (b), results for different potential energy barriers are shown: NS with 9.6 meV (blue dot), 23.6 meV (red dot); LS with 9.6 meV (orange star) and 23.6 meV (purple star).

Second, the next step for checking our numerical method is within the Brownian regime. For this goal, the ISF I(Δk, t) for the same system with an energy barrier E_b = 23.6 meV, surface temperature T = 105 K and initial velocity v₀ = 115 m s⁻¹ is calculated. As has previously seen, in this regime, this function is an exponential function and has been fitted to exp(−αt). In Fig. 13, the decay rate α is plotted versus Δk. The analysis has to be carried out in terms of the so-called Chudley–Elliott model⁵⁹ which describes this jump diffusion by means of a rate equation according to

where τ_j is the average time between successive jumps over the one-dimensional vector j and the summation runs over all lattice vectors. In this model, it is assumed that the time for a simple jump is very short compared with the time τ between successive jumps. Thus, the total jump rate is , with τ_j = τ_−j. Now, due to the linearity property of the Fourier transformthe solution of this first order differential equation beingthat is, the effect to include a simple Bravais lattice is to add an exponential factor ruling the jumps in the diffusion process. Thus,Numerical results are the red dots and eqn (94) is plotted for j = ±1 (black curve) and j = ±6 (blue curve). As expected, at small values of Δk, one observes a quadratic behavior. This is well reproduced by considering j = ±1. At higher values of Δk more and more values of j are needed. This is explained because a very small friction coefficient has been considered.

Fig. 13 The decay rate α versus Δk for the Xe–Pt(111) system is displayed for a corrugated surface with an energy barrier E_b = 23.6 meV, surface temperature T = 105 K and initial adsorbate velocity v₀ = 115 m s⁻¹. The corresponding numerical results are analyzed in terms of the so-called Chudley–Elliott model given by eqn (94). Numerical results (red dots) and Chudley–Elliott curves for j = ±1 (black curve) and six j = ±6 (blue curve) are plotted.

7 Conclusions

In this work, we have used the CL formalism as well as the stochastic wave function method for the first time in surface diffusion. As has been discussed along this work, the splitting of the G-function in G_s and G_d is no longer valid when considering coherent quantum surface corrugation. In particular, for such a goal, we have analyzed the cat states and identical adsorbates, bosons and fermions. The ISF function is not Gaussian in the ballistic regime but it keeps the exponential decay in the Brownian regime. This study has been carried out first for flat surfaces where analytical results are obtained. The stochastic wave function numerical method has also been used to compare with analytical results. In order to show that this numerical method can be an alternative to Langevin simulations for corrugated surfaces, model calculations have been carried out which have been analyzed in terms of the so-called Chudley–Elliott model. The agreement in all of cases studied here is fairly good. This should serve us to extend this method to different realistic systems such as for example Na–Cu(111) and Li–Cu(111) where available experimental results also exist.^42,58,60 Furthermore, it should be very interesting to extend and analyze systems to light adsorbates where surface tunneling is present. The stochastic wave function method should be then modified to consider low surface temperatures. Work in this this direction is now in progress.

Author contributions

All authors have contributed equally to this work.

Conflicts of interest

There are no conflicts to declare.

Appendix A: Gaussian ansatz coefficients for one particle

From the initial wave packet (9), the initial density matrix can be extractedIf one replacing the Gaussian ansatz (12) into the master eqn (7), a system of differential equations for the coefficients A(t), B(t), C(t), D(t), E(t) and N(t) we have thatFrom the initial density matrix (95), the initial conditions for the coefficients can be extractedWith the initial conditions (101), the solutions of the differential equations for the coefficients (96) are

Appendix B: Wigner representation

This representation is often used as an alternative to the density matrix for systems described by a continuous degree of freedom. This alternative way of the reduced density matrix used by Unruh and Zurek⁶¹ allows us to carry out this analysis in the (k, Δ) space. In this representation, the reduced density matrix is the characteristic function associated with the Wigner function^62,63 and is defined by

ρ(k, Δ) = Tr{ρ exp i(kx + Δp/ħ)}, (103)

where the Wigner function is given byand both representations are connected through a Fourier transform according towhere Δ = x − x′ measures the distance from the diagonal in positions and is the wave number in the direction parallel to the diagonal⁶⁴. The CL master eqn (7) is then written aswhose solution is

ρ(k, Δ, t) = exp{−c₁(t)k² − c₂(t)kΔ − c₃(t)Δ² − ic₄(t)k − ic₅(t)Δ − c₆(t)}. (107)

A system of linear differential equations for the coefficients c_i(t) is obtained by substituting eqn (107) into eqn (106)

The initial conditions for this coefficients can be extracted from the Fourier transform of the density matrix for the initial wave packet (9) written in the center of mass and relative coordinatesBy carrying out the space Fourier transform of both argument r and R to the variables Δ and kand comparing (110) with the density matrix (107) the following initial conditions for the coefficients c_i(t) areFrom eqn (111), the solution of the differential eqn (108) areThe coefficients c₁, c₂, c₃, c₄, c₅ and c₆ are related to the coefficients of eqn (12) through

Appendix C: limit cases for the characteristic magnitudes in decoherence dynamics

C.1 Coherence length

In the absence of environmental interactions (D = 0 and γ = 0), the coherence length becomesThis expression describes the free coherent spread of the wave packet for a free particle. However, the time evolution of the coherence length (25) for times t ≪ γ⁻¹ isDue to the influence of environmental scattering, at short times, the coherence length has a constant value equal to the initial width σ₀ and later on increases linearly with time. At long times t ≫ γ⁻¹, the coherence length has an asymptomatic behaviorindependently on the initial width σ₀. This value is inversely proportional to the diffusion coefficient D, which includes the particle mass and temperature (8), and it is directly proportional to the friction coefficient γ.

C.2 Ensemble width

The probability distribution is given bythe first two moments beingandhence the variance in position is entirely given by the function C asthenwhich quantifies the total size of the position-space ensemble. In the absence of any environmental interactions (D = 0 and γ = 0), the free ensemble width isFor short times (t ≪ γ⁻¹), the time evolution of the ensemble width (26) isThere is no difference between the ensemble width with or without environment. However, at long times t ≫ γ⁻¹, the ensemble width has an asymptomatic behavior

Appendix D: Gaussian ansatz coefficients for Schrödinger cat states

For the initial wave packet (28), the density matrix is given by (29) and can be expressed as a sum of four contributions, since the Markovian master eqn (7) is linear,where each term is the solution of (7) for the corresponding initial density matrix.

Using a Gaussian ansatz (32) for each term and inserting (129) into the equation of motion eqn (7), the coupled differential equations for the coefficients A_i(t), B_i(t), C_i(t), D_i(t), E_i(t) and N_i(t) with the initial conditions (29) can be solved. The coefficients of each contribution are

where i = 1, 2, 3, 4 are the corresponding subscripts 1p_x0, 1p_−x0, 2p₁₂, 2p₂₁. The coefficients A_i(t), B_i(t) and C_i(t) are the same in the four terms. However, the coefficients D_i(t) are differentThe same happens for the E_i(t) coefficientsand N_i(t):

Appendix E: Itô stochastic differential equation

Given a specific base of wave functions, the density matrix can be expressed as a linear combination of normalized wave function |ψ〉 from the basewhere M represent the average over the number of wave functions in the Hilbert space [script letter H]_S.

For an open system, the environmental interaction gives a non-deterministic evolution to the system density matrix. The probabilistic nature of the interaction turns the master equation for the density matrix into a stochastic differential equation. After a time dt, the variation of the wave function |dψ〉 is given by

where |v〉dt is a diffusive term and |u_j〉dξ_j a stochastic term given by the Wiener independent processes Ñ.⁵⁶ This expression is known as the Itô differential equation. In order to preserve the normalization of the state vector |ψ〉, it should be orthogonal to the fluctuations |u_j〉

〈ψ|u_j〉 = 0. (145)

From eqn (144) and the Wiener propertiesthe diffusive and stochastic terms can be found from (143) to be

dρ_S = M(|ψ〉〈dψ| + |dψ〉〈ψ| + |dψ〉〈dψ|), (147)

and by using eqn (146), the time evolution of the density matrix is given byThe vector |u_j〉 is the component of [small rho, Greek, dot above]_S in the orthogonal space to |ψ〉where Î_S is the identity operator. Now, from

[small rho, Greek, dot above]_S|ψ〉 = |ψ〉〈v|ψ〉 + |v〉, (150)

and applying the scalar product properties

〈ψ|[small rho, Greek, dot above]_S|ψ〉 = 2 Re〈ψ|v〉, (151)

the diffusive term can be written aswhere ic is an imaginary constant representing a phase change and his value is determined in such a way that, in the absence of interaction with the environment, the Schrödinger equation is recovered.

If the Lindblad master equation²²

is replaced in eqn (152), the diffusive term can be obtainedwhere 〈Â〉 = 〈ψ|Â|ψ〉 represent the operator average value  and [small eta, Greek, tilde]_k are non-negative quantities. Once the diffusive term has been found, the stochastic term |u_k〉 of the differential wave function eqn (144) can be found by replacing (153) in (149) and using the scalar product properties, leading to the following expression

|u_k〉 = (Â_k − 〈Â_k〉)|ψ〉. (155)

Thus, one has thatand replacing the Lindblad operator (79), the time evolution of the state vector |ψ〉 can be found by solving the stochastic differential eqn (83), given by a deterministic part and a stochastic term.

Acknowledgements

E. E. T.-M. and S. M. A. would like to thank Fundación Humanismo y Ciencia for financial support and Spanish Project PID2021-125735NB-I00. G. R.-L. and J. R.-S. would also thank the Advanced Computational Team at Higher Institute for Technologies and Applied Sciences (Havana University) for the support provided during the realization of this work and to Project NA223LH-INSTEC-003. Finally, we would like to thank the two anonymous reviewers for very useful comments and suggestions to improve this work.

Notes and references

1. L. Van Hove, Correlations in space and time and Born approximation scattering in systems of interacting particles, Phys. Rev., 1954, 95(1), 249.
2. S. W. Lovesey, Theory of neutron scattering from condensed matter, 1984, vol. 2.
3. D. A. McQuarrie, Statistical Mechanics, 1976.
4. G. H. Vineyard, Scattering of slow neutrons by a liquid, Phys. Rev., 1958, 110(5), 999.
5. D. J. Ward, A. Raghavan, A. Tamtögl, A. P. Jardine, E. Bahn, J. Ellis, S. Miret-Artés and W. Allison, Inter-adsorbate forces and coherent scattering in helium spin-echo experiments, Phys. Chem. Chem. Phys., 2021, 23(13), 7799–7805.
6. F. Hofmann and J. P. Peter Toennies, High-resolution helium atom time-of-flight spectroscopy of low-frequency vibrations of adsorbates, Chem. Rev., 1996, 96(4), 1307–1326 , PMID: 11848791.
7. A. P. Graham, The low energy dynamics of adsorbates on metal surfaces investigated with helium atom scattering, Surf. Sci. Rep., 2003, 49(4), 115–168.
8. I. Calvo-Almazán and P. Fouquet, The application of quasi-elastic neutron scattering techniques (QENS) in surface diffusion studies, Eur. Phys. J. Spec. Top., 2012, 213(4), 149–163.
9. A. P. Jardine, G. Alexandrowicz, H. Hedgeland, W. Allison and J. Ellis, Studying the microscopic nature of diffusion with helium-3 spin-echo, Phys. Chem. Chem. Phys., 2009, 11, 3355–3374.
10. A. P. Jardine, H. Hedgeland, G. Alexandrowicz, W. Allison and J. Ellis, Helium-3 spin-echo: Principles and application to dynamics at surfaces, Prog. Surf. Sci., 2009, 84(11), 323–379.
11. S. Miret-Artés and E. Pollak, The dynamics of activated surface diffusion, J. Phys.: Condens. Matter, 2005, 17(49), S4133–S4150.
12. J. L. Vega, R. Guantes and S. Miret-Artés, Quasielastic and low vibrational lineshapes in atom–surface diffusion, J. Phys.: Condens. Matter, 2004, 16(29), S2879–S2894.
13. P. S. M. Townsend and D. J. Ward, The intermediate scattering function for quasi-elastic scattering in the presence of memory friction, J. Phys. Commun., 2018, 2(7), 075011.
14. R. Martínez-Casado, J. L. Vega, A. S. Sanz and S. Miret-Artés, Line shape broadening in surface diffusion of interacting adsorbates with quasielastic He atom scattering, Phys. Rev. Lett., 2007, 98, 216102.
15. R. Martínez-Casado, J. L. Vega, A. S. Sanz and S. Miret-Artés, Surface diffusion and low vibrational motion with interacting adsorbates: A shot noise description, Phys. Rev. E, 2007, 75, 051128.
16. D. J. Ward, A. Raghavan, A. Tamtögl, A. P. Jardine, E. Bahn, J. Ellis, S. Miret-Artés and W. Allison, Inter-adsorbate forces and coherent scattering in helium spin-echo experiments, Phys. Chem. Chem. Phys., 2021, 23, 7799–7805.
17. R. Martínez-Casado, A. S. Sanz, J. L. Vega, G. Rojas-Lorenzo and S. Miret-Artés, Linear response theory of activated surface diffusion with interacting adsorbates, Chem. Phys., 2010, 370(1), 180–193.
18. S. Miret-Artés, Quantum surface diffusion in Bohmian mechanics, J. Phys. Commun., 2018, 2(9), 095020.
19. A. O. Caldeira and A. J. Leggett, Path integral approach to quantum Brownian motion, Phys. A, 1983, 121, 374.
20. G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys., 1976, 48(2), 119–130.
21. A. Kossakowski, On quantum statistical mechanics of non-Hamiltonian systems, Rep. Math. Phys., 1972, 3, 247.
22. H. P. Breuer and F. Petruccione, The theory of open quantum systems, 2002.
23. E. E. Torres-Miyares, G. Rojas-Lorenzo, J. Rubayo-Soneira and S. Miret-Artés, Surface diffusion by means of stochastic wave functions. The ballistic regime, Mathematics, 2021, 9(4), 362.
24. C. W. Gardiner and P. Zoeller, Quantum noise, Berlin, Heidelberg, New York, 1991.
25. A. O. Caldeira and A. J. Leggett, Influence of dissipation on quantum tunneling in macroscopic systems, Phys. Rev. Lett., 1981, 46(4), 211.
26. C. Chan, Theory and application of open quantum systems, 2012.
27. S. F. Huelga, et al., Open quantum systems: An introduction, 2012.
28. H. Grabert, P. Schramm and G. L. Ingold, Quantum Brownian motion: The functional integral approach, Phys. Rep., 1988, 168(3), 115.
29. D. Suess, A. Eisfeld and W. T. Strunz, Hierarchy of stochastic pure states for open quantum system dynamics, Phys. Rev. Lett., 2014, 113(15), 150403.
30. Y. Yan, F. Yang, Y. Liu and J. Shao, Hierarchical approach based on stochastic decoupling to dissipative systems, Chem. Phys. Lett., 2004, 395(4-6), 216–221.
31. M. A. Lane, D. Matos, I. J. Ford and L. Kantorovich, Exactly thermalized quantum dynamics of the spin-boson model coupled to a dissipative environment, Phys. Rev. B, 2020, 101(22), 224306.
32. R. Feynman, R. B. Leighton and M. Sands, The Brownian movement, Feynman Lect. Phys., 1964, 1, 41.
33. A. O. Caldeira and A. J. Leggett, Quantum tunneling in a dissipative system, Ann. Phys., 1983, 149(2), 374.
34. U. Weiss, Quantum Dissipative Systems, World Scientific Publishing, Singapore, 2nd edn, 1999.
35. M. A. Schlosshauer, Decoherence and the quantum-to-classical transition, 2007.
36. V. I. Tatarskiĭ, The Wigner representation of quantum mechanics, Soviet Physics Uspekhi, 1983, 26(4), 311.
37. G. W. Ford, J. T. Lewis and R. F. O'Connell, Quantum measurement and decoherence, Phys. Rev. A, 2001, 64, 032101.
38. G. W. Ford and R. F. O'Connell, Decoherence at zero temperature, J. Opt. B: Quantum Semiclassical Opt., 2003, 5(6), S609–S612.
39. H. C. Peñate-Rodríguez, R. Martínez-Casado, G. Rojas-Lorenzo, A. S. Sanz and S. Miret-Artés, Quantum zeno and anti-zeno effects in surface diffusion of interacting adsorbates, J. Phys.: Condens. Matter, 2012, 24(10), 104013.
40. D. J. Ward, A study of spin-echo lineshapes in helium atom scattering from adsorbates, 2013.
41. J. Ellis, A. P. Graham and J. P. Toennies, Quasielastic helium atom scattering from a two-dimensional gas of Xe atoms on Pt(111), Phys. Rev. Lett., 1999, 82(25), 5072.
42. R. Martínez-Casado, J. L. Vega, A. S. Sanz and S. Miret-Artés, Quasielastic He atom scattering from surfaces: A stochastic description of the dynamics of interacting adsorbates, J. Phys.: Condens. Matter, 2007, 19(30), 305002.
43. E. Joos and H. D. Zeh, The emergence of classical properties through interaction with the environment, Z. Phys. B: Condens. Matter, 1985, 59(2), 223–243.
44. S. V. Mousavi and S. Miret-Artés, On some unexplored decoherence aspects in the Caldeira-Leggett formalism: arrival time distributions, identical particles and diffraction in time, Eur. Phys. J. Plus, 2022, 137, 1–14.
45. W. H. Zurek, S. Habib and J. P. Paz, Coherent states via decoherence, Phys. Rev. Lett., 1993, 70(9), 1187.
46. T. Qureshi and A. Venugopalan, Decoherence and matter wave interferometry, Int. J. Mod. Phys. B, 2008, 22(08), 981–990.
47. W. H. Zurek, Pointer basis of quantum apparatus: Into what mixture does the wave packet collapse?, Phys. Rev. D, 1981, 24(6), 1516.
48. W. H. Zurek, Environment-induced superselection rules, Phys. Rev. D, 1982, 26(8), 1862.
49. S. Treiman, The odd quantum, in The Odd Quantum, Princeton University Press, 2002.
50. S. V. Mousavi and S. Miret-Artés, Dissipative quantum backflow, Eur. Phys. J. Plus, 2020, 135(3), 1–18.
51. S. V. Mousavi and S. Miret-Artés, Quantum-classical comparison: arrival times and statistics, Phys. Scr., 2015, 90(2), 025001.
52. S. Miret-Artés, R. S. Dumont, T. Rivlin and E. Pollak, The influence of the symmetry of identical particles on flight times, Entropy, 2021, 23(12), 1675.
53. V. Ambegaokar, Quantum Brownian motion and its classical limit, Ber. Bunseng. Phys. Chem., 1991, 95(3), 400–404.
54. A. Tameshtit and J. E. Sipe, Positive quantum Brownian evolution, Phys. Rev. Lett., 1996, 77(13), 2600.
55. S. Gao, Lindblad approach to quantum dynamics of open systems, Phys. Rev. B, 1998, 57(8), 4509.
56. R. Durrett, Probability: theory and examples, 2019, vol. 49.
57. J. Halliwell and A. Zoupas, Quantum state diffusion, density matrix diagonalization, and decoherent histories: A model, Phys. Rev. D, 1995, 52(12), 7294.
58. D. J. Ward, A study of spin-echo lineshapes in helium atom scattering from adsorbates, 2013.
59. R. J. Elliot and C. T. Chudley, Neutron scattering from a liquid on a jump diffusion model, Proc. Phys. Soc., 1961, 77(2), 353.
60. J. G. Mantecón, S. Miret-Artés and G. Rojas-Lorenzo, Diffusion of Ar atoms on MgO(100) surfaces, Revista Cubana de Física, 2012, 29(1), 8–13.
61. W. G. Unruh and W. H. Zurek, Reduction of a wave packet in quantum Brownian motion, Phys. Rev. D, 1989, 40(4), 1071.
62. M. Hillery, R. F. O'Connell, M. O. Scully and E. P. Wigner, Distribution functions in physics: Fundamentals, Phys. Rep., 1984, 106(3), 121–167.
63. C. M. Savage and D. F. Walls, Damping of quantum coherence: the master-equation approach, Phys. Rev. A, 1985, 32(4), 2316.
64. E. Joos, H. D. Zeh, C. Kiefer, D. J. W. Giulini, J. Kupsch and I. O. Stamatescu, Decoherence and the appearance of a classical world in quantum theory, 2013.

---