Surface Diffusion by Means of Stochastic Wave Functions. The Ballistic Regime

Abstract

Stochastic wave function formalism is briefly introduced and applied to study the dynamics of open quantum systems; in particular, the diffusion of Xe atoms adsorbed on a Pt(111) surface. By starting from a Lindblad functional and within the microscopic Caldeira–Leggett model for linear dissipation, a stochastic differential equation (It$\widehat{o}$-type differential equation) is straightforwardly obtained. The so-called intermediate scattering function within the ballistic regime is obtained, which is observable in Helium spin echo experiments. An ideal two-dimensional gas has been observed in this regime, leading to this function behaving as a Gaussian function. The influence of surface–adsorbate interaction is also analyzed by using the potential of two interactions describing flat and corrugated surfaces. Very low surface coverages are considered and, therefore, the adsorbate–adsorbate interaction is safely neglected. Good agreement is observed when our numerical results are compared with the corresponding experimental results and previous standard Langevin simulations.

1. Introduction

Most physical and chemical systems in nature are open classical or quantum systems. The corresponding dynamics are usually well described by considering an environment with an infinite number of degrees of freedom [1,2]. The dynamics of open quantum systems [3] is a central issue in different areas of physics, such as materials science, atomic, molecular and statistical physics of complex systems. The interaction between a system and its environment often leads to dissipation, quantum fluctuations and the irreversible evolution of the system. These processes affect surface dynamics, surface tunneling, non-adiabatic effects, and so on. Many disciplines, from physics to biology, have developed increasingly powerful methods for modeling open quantum systems [4,5].

The adsorption of rare gases has been of interest from the pioneering work by Langmuir [6]. The diffusion of atoms/molecules on metal surfaces is usually considered in the quantum regime due to internal and/or external vibrations of adsorbates, as well as energy fluctuations at microscopic scale [7,8]. These studies provide valuable information on adsorbate–substrate and adsorbate–adsorbate interactions. Thus, it seems to be essential to study these phenomena using the theoretical formalism of open quantum systems in order to better understand the main microscopic and elementary processes occurring at surfaces. Moreover, these processes can be seen many times as a preliminary step in the study of more complicated phenomena as, for example, in heterogeneous catalysis [9], crystal growth, chemical vapor deposition [7], desorption [10], adhesion [11], photon–electron- induced surface reactions [12], etc.

Depending on the surface temperature, surface coverage and the adsorbate–substrate interaction, diffusion times can vary by orders of magnitude [13,14]. The experimental techniques used for the characterization of surface diffusion can be classified according to the diffusion time that they are capable of measuring, such as fast or slow diffusion. Scanning tunneling microscopy (STM) [15] and field ion microscopy (FIM) [16] are especially useful for slow diffusion, where the time between jumps is of the order of seconds. These measurements provide a direct observation of the diffusion process. Among the different experimental techniques used to study these processes and also extract information indirectly about interaction potentials, the so-called quasielastic He atom scattering (QHAS) [17,18,19,20] is considered as the surface science analogue of quasielastic neutron scattering, which is widely and successfully applied to analyze fast diffusion (time between jumps of the order of microseconds). In order to interpret QHAS results, a more elaborate theory is needed to have reliable adsorbate–adsorbate and adsorbate–substrate interactions. A more recent experimental and powerful technique is also available, the so-called He spin echo (HeSE), which is similar to the neutron spin echo with ultra-high energy resolution [21]. In the QHAS technique, the so-called dynamic structure factor or scattering law is directly measured, which provides us with line shapes for three elementary processes: the adsorbate diffusion process, low frequency adsorbate internal motions, and surface phonon excitations. The inverse Fourier transform in the frequency domain of the dynamic structure factor is the so-called intermediate scattering function which is an observable in the HeSE technique. This function provides us with the same information as the dynamic structure factor but in the time domain.

During the past few decades, many efforts have also been carried out to devise various phenomenological models [22] and more recently to derive via first-principles the dissipation functional from microscopic Hamiltonian [23]. Pertinent to surface dynamics, the effect of surface degrees of freedom into atom/molecule-surface dynamics within the density matrix approach is included [24,25,26]. In this formalism, microscopic events, such as quantum jumps, which exist in open individual systems, disappear by ensemble average. An alternative and less explored method consists of using a set of stochastic wave functions. The master equation for this method is a particular case of stochastic differential equations (SDE) known as the It$\widehat{o}$ equation [4]. This method reduces the dimensional problem from $N^2$ to N. In this sense, the SDE provides, aside from its statistical equivalent to the master equation, a more convenient description of individual systems. This is important for the study of quantum noise in mesoscopic [27] systems and individual atoms/molecules [28]. As far as we know, this is the first time that this theoretical analysis has been carried out in the surface diffusion context.

In particular, in this work, we have analyzed the QHAS experimental results by Ellis et al. [20] where, for the first time, experimental evidence for a two-dimensional (2D) ideal gas of Xe atoms on Pt(111) at low coverage ($\theta =0.017$), low incident helium atom energy ($E_i=10.15$ meV) and surface temperature $T=105$ K is observed. This observation implies that Xe atoms behave as independent particles on the Pt(111) surface during a certain time, the well-known ballistic regime in diffusion processes. This regime is established for times $t<<{\eta }^{-1}$, where $\eta$ is the friction coefficient. This study is carried out for a flat and corrugated surface and the corresponding results are compared with experimental and previous Langevin numerical results within the ballistic regime [14,20,29].

The paper is organized as follows. In Section 2, we briefly present the Linblad formalism, as well as the equivalent stochastic differential equation for open quantum systems. In Section 3, numerical results are presented and discussed in the ballistic regime for flat and corrugated surfaces. Finally, Section 4 summarizes the conclusions that can be drawn from this work.

2. Theory

2.1. The Intermediate Scattering Function

A long time ago, Van Hove [30] described the differential cross section of the scattering of slow neutrons by a system of interacting particles in terms of the generalized pair distribution function, the so-called $G(\mathbf{r},t)$ function of van Hove (with $\mathbf{r}$ being a position vector and t a time interval). In the Born approximation, the scattering event is reduced mainly to a statistical mechanics problem [30,31,32], the nature of the scattered particles (neutrons, light, atoms, etc.) and details of the interaction potential with the system under study being irrelevant (the linear response theory). In this formalism, the linear response of the system implies that it is determined entirely by the properties exhibited by the system in the absence of probe particles. The differential cross section can be written in terms of the momentum transfer, $\hslash \Delta \mathbf{k}$, and energy transfer $\hslash \omega$ as [30]

$$\begin{array}{c}\hfill \frac{d^2\mathcal{R}(\Delta \mathbf{k},\omega )}{d\mathrm{\Omega }d\omega }\propto S(\Delta \mathbf{k},\omega )={\left(2\pi \right)}^{-1}N\int e^{i(\Delta \mathbf{k}.\mathbf{r}-\omega t)}G(\mathbf{r},t)d\mathbf{r}dt\end{array}$$

(1)

which gives the probability that the probe particles scattered from the diffusing system reach a certain solid angle $\mathrm{\Omega }$ in an interval of outgoing energy $\hslash \omega$. The response function $S(\Delta \mathbf{k},\omega )$ is also termed the scattering law or dynamic structure factor (DSF), N being the number of interacting particles in the system under study. As is well known, the spatial Fourier transform of the G-function

$$\begin{array}{c}\hfill I(\Delta \mathbf{k},t)=\int e^{i\Delta \mathbf{k}.\mathbf{r}}G(\mathbf{r},t)d\mathbf{r}\end{array}$$

(2)

is called the intermediate scattering function (ISF) and, therefore, S and I are related by the inverse Fourier transform in frequency.

For surface diffusion, Helium atoms are generally used for probing the dynamics of adsorbates or adparticles on surfaces. Due to the scattering, He atoms undergo an energy exchange $\hslash \omega =E_{final}-E_{initial}$ and a parallel (to the surface) wave vector transfer $\Delta \mathbf{K}={\mathbf{K}}_{final}-{\mathbf{K}}_{initial}$. Capital letters are used here for variables in the surface plane. The prominent peak displayed by the DSF around the zero energy transfer or quasi-elastic peak provides indirect information of the diffusion constant. Satellite weaker peaks at low energy transfers observed are attributed to low frequency motions of some adsorbates, as well as surface phonon excitations. Furthermore, long distance and time correlations in the interacting system are extracted from the scattering law when considering small values of $\Delta \mathbf{K}$ and $\hslash \omega$, respectively. The DSF dynamic is usually expressed as [30]

$$\begin{array}{c}\hfill S(\mathbf{K},\omega )={\left(2\pi \right)}^{-1}N\int e^{-i\omega t}I(\mathbf{K},t)dt,\end{array}$$

(3)

with

$$\begin{array}{c}\hfill I(\mathbf{K},t)\equiv \frac{1}{N}\langle \sum _{j,j^{\prime }}^N{e}^{-i\mathbf{K}·{\mathbf{R}}_j\left(0\right)}{e}^{i\mathbf{K}·{\mathbf{R}}_{j^{\prime }}\left(t\right)}\rangle \end{array}$$

(4)

where the brackets denote an ensemble average and ${\mathbf{R}}_j\left(t\right)$ the position vector of the j adparticle at time t on the surface. From a theoretical point of view, one of the main goals is to calculate the corresponding trajectories. For open quantum systems, there are several ways to tackle such a calculation. The standard and generalized Langevin equation formalism are widely used. However, as far as we know, the so-called stochastic wave function formalism has not been used in this context.

2.2. The Linblad Formalism

In the density matrix formalism, the time evolution of the reduced density matrix $\rho$ of an open quantum system is ruled by the equation of motion

$$\begin{array}{c}\hfill \frac{d\rho \left(t\right)}{dt}+\frac{i}{\hslash }[\widehat{H},\rho \left(t\right)]=\mathcal{L}\rho \left(t\right),\end{array}$$

(5)

where $\widehat{H}$ and $\mathcal{L}$ are the Hamiltonian and Liouville operators of the system. The derivation of the Liouville operator $\mathcal{L}$ from microscopic Hamiltonians usually leads to a quantum master equation. From the pioneering work of Caldeira and Leggett (CL) [23] for an oscillator linearly coupled to an Ohmic environment (linear dissipation), the master equation is written as

$$\begin{array}{c}\hfill \frac{d\rho }{dt}+\frac{i}{\hslash }[\widehat{H},\rho ]=-\frac{\eta }{\hslash }\left\{\frac{2m{k}_{B}T}{\hslash }\left[\widehat{x},[\widehat{x},\rho ]\right]+i\left[\widehat{x},{[\widehat{p},\rho ]}_{+}\right]\right\},\end{array}$$

(6)

where T and $\eta$ represent the bath temperature (here, the surface temperature) and friction coefficient between the oscillator and the Ohmic bath. Similar master equations have been obtained along this line for a particle in a more general environment [33] with linear and nonlinear couplings. The master Equation (6) is known to have a few drawbacks. First, it is only valid at high temperatures, or equivalently, in the classical limit. The quantum mechanical behavior is known to be important mainly in the low-temperature regime and, therefore, it is important to extend this model to such a regime. Additionally, second, it deals with the non-positivity of $\rho$ on a time scale proportional to $\tau \sim {\eta }^{-1}$, as shown in different works [33,34].

There have been several attempts to derive master equations, preserving the positive density evolution from microscopic Hamiltonians, and they can be found in the existing vast literature. Such derivations are mainly limited to the so-called weak-coupling [8] or high-temperature regimes [33]. However, the non-positivity problem is not found in the theory of quantum dynamical semigroups by Lindblad [35], Kossakowski [36] and Gorini and Sudarshan [2]. In particular, it was shown that the generator for a completely positive map should be governed by

$$\begin{array}{c}\hfill \mathcal{L}\rho =\sum _k\left\{[{\widehat{A}}_k,\rho {\widehat{A}}_k^{+}]+[{\widehat{A}}_k\rho ,{\widehat{A}}_k^{+}]\right\},\end{array}$$

(7)

where ${\widehat{A}}_k$ are known as the Lindblad dissipation operators. However, these operators are, in general, unknown, and the compact structure of the generator (7) does not generally assure system equilibrium with the bath.

2.3. Construction of the Lindblad Operators and Extension to Low Temperature

A single Lindblad operator $\widehat{A}$ is generally taken to be a linear combination of the operator position $\widehat{x}$ and the momentum $\widehat{p}$ as

$$\begin{array}{c}\hfill \widehat{A}=\mu \widehat{x}+i\nu \widehat{p},{\widehat{A}}^{+}=\mu \widehat{x}-i\nu \widehat{p}.\end{array}$$

(8)

where $\mu$ and $\nu$ are arbitrary complex numbers. This particular choice is motivated by considering (i) the form of the dissipation operator for a damped harmonic oscillator with mass m and frequency $\mathrm{\Omega }$ [36], that is, ${\widehat{a}}_d=\sqrt{m\mathrm{\Omega }/2\hslash }[\widehat{x}+(i/m\mathrm{\Omega })\widehat{p}]$, which is a linear combination of both operators, (ii) the dissipation terms of the master Equation (6), and (iii) Dekker’s constraints [37] for the diffusion coefficients of the master equation for a damped oscillator.

The Lindblad operators (8) in the generator (7) are replaced by the master Equation (5) according to

$$\begin{array}{cc}\hfill \frac{d\rho }{dt}+\frac{i}{\hslash }[{\widehat{H}}^{\prime },\rho ]& =-{\mu }^2\left[\widehat{x},[\widehat{x},\rho ]\right]-2i\mu \nu \left[\widehat{x},{[\widehat{p},\rho ]}_{+}\right]-{\nu }^2\left[\widehat{p},[\widehat{p},\rho ]\right],\hfill \\ & {\widehat{H}}^{\prime }=\widehat{H}-2\mu \nu \hslash \widehat{x}\widehat{p}.\hfill \end{array}$$

(9)

The first two terms on the right-hand side of (9) are essentially the same as those in the CL master Equation (6). This comparison suggests the choice of $\mu$ and $\nu$ as follows

$$\begin{array}{cc}{\mu }^2& =2kTm\eta /{\hslash }^2,\hfill \\ 2\mu \nu & =\eta /\hslash ,\hfill \\ {\nu }^2& =\eta /\left(8mkT\right)\underset{T\to \infty }{\overset{}{\to }}0.\hfill \end{array}$$

(10)

The coefficient ${\nu }^2$ is negligible at high-temperatures but not at low temperatures. The master Equation (9), together with the coefficients given by Equation (10) represent the simplest Lindblad master equation [33]. It provides a way to overcome the non-positivity problem of the CL model. Remember that this choice is only valid at high temperatures. In order to use the two coefficients given by Equation (10) at low temperatures, Equation (9) has to be closely analyzed. The two dissipation terms of (6) have a different physical origin. The second term on the right side corresponds to the dissipation induced by the bath and is temperature independent

$$\begin{array}{c}\hfill 2\mu \nu =\eta /\hslash .\end{array}$$

(11)

Unlike the second term, the first one depends on T and describes the environment-induced fluctuation. The origin of this term comes from the high temperature approximation.

Thus, in order to extend the master equation to the low-temperature regime, the CL master Equation (6) is derived within the path-integral approach [23] and theoretical physics methods [33]. The generator $\mathcal{L}$ is then written in terms of the so-called noise $D_1\left(\tau \right)$ and dissipation $D\left(\tau \right)$ kernels as

$$\begin{array}{c}\hfill \mathcal{L}\rho =-\frac{1}{{\hslash }^2}{\int }_0^{\infty }d\tau (iD\left(\tau \right)\left[\widehat{x},{[\widehat{x}(-\tau ),\rho ]}_{+}\right]-D_1(-\tau )\left[\widehat{x},[\widehat{x}(-\tau ),\rho ]\right]),\end{array}$$

(12)

with

$$\begin{array}{cc}\hfill D_1\left(\tau \right)& =\hslash {\int }_0^{{\mathrm{\Omega }}_c}J\left(\omega \right)cos\left(\omega \tau \right)coth(\hslash \omega /2k_{B}T)d\omega ,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill D\left(\tau \right)& =\hslash {\int }_0^{{\mathrm{\Omega }}_c}J\left(\omega \right)sin\left(\omega \tau \right)d\omega ,\hfill \end{array}$$

(14)

$J\left(\omega \right)={\sum }_{n=1}^{N_1}(k_n^2/2m_n{\omega }_n)\delta (\omega -{\omega }_n)$ being the spectral function of the phonon bath and ${\mathrm{\Omega }}_c$ the bandwidth cutoff frequency. For an Ohmic bath, $J\left(\omega \right)$ is usually approximated by $J\left(\omega \right)=\eta \omega /\pi$. A different approximation to the noise kernel can be used at low temperatures by observing that $J\left(\omega \right)coth(\hslash \omega /2k_{B}T)$ is a smooth function of $\omega$, while $cos\left(\omega \tau \right)$ is fast oscillating, then

$$\begin{array}{cc}\hfill D_1\left(\tau \right)& \approx \frac{4\hslash m\eta }{\pi }{\omega }_{c}coth\left(\frac{\hslash {\omega }_c}{2k_{B}T}\right){\int }_0^{{\mathrm{\Omega }}_c}cos\left(\omega \tau \right)d\omega \hfill \\ & =4m\eta \hslash {\omega }_{c}coth\left(\frac{\hslash {\omega }_c}{2k_{B}T}\right)\tilde{\delta }\left(\tau \right),\hfill \end{array}$$

(15)

where $\tilde{\delta }\left(t\right)=1/\pi {\int }_0^{{\mathrm{\Omega }}_c}cos\left(\omega t\right)d\omega$ with ${\omega }_c$ being the center of the $J\left(\omega \right)$ band. The Markovian limit is easily recovered when ${\mathrm{\Omega }}_c\to \infty$ and then $\tilde{\delta }\left(t\right)\to \delta \left(t\right)$. Thus, $2k_{B}T\to \hslash {\omega }_{c}coth\left(\hslash {\omega }_c/2k_{B}T\right)$ in the first diffusion term of (6) and ${\mu }^2$ is given by

$$\begin{array}{c}\hfill {\mu }^2\left(T\right)=\frac{2\eta m{\omega }_c}{\hslash }coth\left(\frac{\hslash {\omega }_c}{2k_{B}T}\right).\end{array}$$

(16)

The parameter ${\omega }_c$ can be uniquely determined by a harmonic oscillator approximation at $T=0$. The Lindblad operator $\widehat{A}$ should then be degenerated with the annihilation operator of the harmonic oscillator, ${\widehat{a}}_d=\sqrt{m\mathrm{\Omega }/2\hslash }[\widehat{x}+(i/m\mathrm{\Omega })\widehat{p}]$, implying $\mu /\nu =m\mathrm{\Omega }$. Furthermore, the relation $2\mu \nu =\eta /\hslash$, valid for all temperatures, and (16) lead to $\mu /\nu =2m{\omega }_c$ at $T=0$. The harmonic approximation thus gives ${\omega }_c=\mathrm{\Omega }/2$. Then, the temperature dependence of the two coefficients reads as

$$\begin{array}{c}\hfill {\mu }^2\left(T\right)=\frac{\eta m\mathrm{\Omega }}{\hslash }coth\left(\frac{\hslash \mathrm{\Omega }}{4k_{B}T}\right),\end{array}$$

(17)

$$\begin{array}{c}\hfill {\nu }^2\left(T\right)=\frac{\eta }{2\hslash m\mathrm{\Omega }}tanh\left(\frac{\hslash \mathrm{\Omega }}{4k_{B}T}\right),\end{array}$$

(18)

with $2\mu /\nu =\gamma /\hslash$. Both expressions reduce to (10) in the high-temperature regime and are based on the fluctuation-dissipation theorem. Therefore, they bring the equilibrium behavior into the Lindblad formalism through their temperature dependence. Equations (17) and (18) fulfill Dekker’s constraints at any temperature.

In coordinate space, the master Equation (9) has the following form -4.6cm0cm

$$\begin{array}{cc}\hfill \frac{\partial \rho (x,x^{{}^{\prime }},t)}{\partial t}+\frac{i}{\hslash }\left[\tilde{H}\left(x\right)-{\tilde{H}}^{*}\left(x^{\prime }\right)\right]\rho (x,x^{\prime },t)=& -\{{\mu }^2\left(T\right){(x-x^{\prime })}^2+\eta (x-x^{\prime })\left(\frac{\partial }{\partial x}-\frac{\partial }{\partial x^{\prime }}\right)\hfill \\ & -{\nu }^2\left(T\right){\hslash }^2{\left(\frac{\partial }{\partial x}+\frac{\partial }{\partial x^{\prime }}\right)}^2\}\rho (x,x^{\prime },t),\end{array}$$

(19)

with

$$\begin{array}{c}\hfill \tilde{H}\left(x\right)=H\left(x\right)+i\hslash \eta x\frac{\partial }{\partial x}+i\frac{\hslash \eta }{2}.\end{array}$$

(20)

This master equation presents two important aspects: (i) a renormalized Hamiltonian with an extra imaginary term involving a frictional force; (ii) an extra diffusion term (the ${\nu }^2$ term) with temperature-dependent diffusion coefficient. In this way, the corresponding equation also guarantees the positivity of the density matrix at short times.

The master Equation (19) can be directly solved in a double-space (two-dimensional) representation [38] but it can also alternatively be solved using a set of stochastic wave functions. This approach can be a possible advantage of the Lindblad scheme, when the dimensionality is a problem, to look for a direct solution.

2.4. Solution of the Master Equation with Stochastic Wave Functions

After Gisin and Percival [39], a master equation of Lindblad class can be equivalently solved by a set of wave functions $\left\{\right|\psi \rangle \}$ by using the representation of the density operator as

$$\begin{array}{c}\hfill \rho =\sum _j|{\psi }_j\rangle \langle {\psi }_j|=M|\psi \rangle \langle \psi |.\end{array}$$

(21)

An N-dimensional state space can be expressed in many ways as a mean M over the distribution of normalized pure-state projection operators. Within the stochastic equation framework, the variation $|d\psi \rangle$ at time $dt$ is given by the Itô form

$$\begin{array}{c}\hfill |d\psi \rangle =|v\rangle dt+\sum _{j=1}^{\tilde{N}}|u_j\rangle d{\xi }_j,\end{array}$$

(22)

where $|v\rangle dt$ is the drift term and the differential stochastic fluctuations are represented by a sum over independent Wiener process $d{\xi }_j$, which are complex with equal and independent fluctuations in their real and imaginary parts. This complex form has simple invariance properties under unitary transformations. As is known, the mean and variance of fluctuations are zero and $d{t}^{1/2}$, respectively, with

$$\begin{array}{c}\hfill M\left[d{\xi }_{k}d{\xi }_n^{*}\right]=2dt{\delta }_{kn}M\left[d{\xi }_k^{*}d{\xi }_n^{*}\right]=0M\left[d{\xi }_{k}d{\xi }_n\right]=0,\end{array}$$

(23)

where M is also used to represent a mean over fluctuations due to the stochastic process.

In order to preserve the normalization of the state vector $|\psi \rangle$, the fluctuations $|u_j\rangle$ in the state must be orthogonal to that state

$$\begin{array}{c}\hfill \langle \psi |u_j\rangle =0.\end{array}$$

(24)

By carrying out means over $|d\psi \rangle$ (22) and $|d\psi \rangle \langle d\psi |$, we have that

$$\begin{array}{c}\hfill M|d\psi \rangle =|v\rangle dt,M|d\psi \rangle \langle d\psi |=2\sum _{j=1}^{\tilde{N}}|u_j\rangle \langle u_j|dt.\end{array}$$

(25)

Now from (25), the change in $\rho$ is given by

$$\begin{array}{c}\hfill d\rho =M\left(|\psi \rangle \langle d\psi |+|d\psi \rangle \langle \psi |+|d\psi \rangle \langle d\psi |\right),\end{array}$$

(26)

and then the equation of motion can be written as

$$\begin{array}{c}\hfill \dot{\rho }=(|\psi \rangle \langle v|+|v\rangle \langle \psi |+2\sum _{j=1}^{\tilde{N}}|u_j\rangle \langle u_j|(.\end{array}$$

(27)

The stochastic terms $|u_j\rangle$ are determined by the component of $\dot{\rho }$ in the space orthogonal to $|\psi \rangle$

$$\begin{array}{c}\hfill 2\sum _{j=1}^{\tilde{N}}|u_j\rangle \langle u_j|=({\widehat{I}}_S-|\psi \rangle \langle \psi \left|\right)\dot{\rho }({\widehat{I}}_S-|\psi \rangle \langle \psi \left|\right),\end{array}$$

(28)

where ${\widehat{I}}_S$ is the identity operator. Note that, although the $|u_j\rangle$ are not uniquely determined by (28), the operator given by the sum over their projections does. This is enough to determine the diffusion process uniquely, since the first and second moments of $|d\psi \rangle$ given by (22) are the same for any $|u_j\rangle$ satisfying (28). The remaining moments can be neglected for an It$\widehat{o}$ process. By multiplying (27) for the state vector $|\psi \rangle$

$$\begin{array}{c}\hfill \dot{\rho }|\psi \rangle =|\psi \rangle \langle v|\psi \rangle +|v\rangle ,\end{array}$$

(29)

and applying the scalar product properties

$$\begin{array}{c}\hfill \langle \psi \left|\dot{\rho }\right|\psi \rangle =2Re\langle \psi |v\rangle ,\end{array}$$

(30)

the drift term is given by

$$\begin{array}{c}\hfill |v\rangle ={\dot{\rho }}_S|\psi \rangle -\left(\frac{1}{2}\langle \psi |{\dot{\rho }}_S|\psi \rangle +ic\right)|\psi \rangle ,\end{array}$$

(31)

where $ic$ is a non-physical imaginary phase change constant, which is determined (by convention) to agree with the usual Schrödinger equation in the absence of interaction with the environment.

The above theory and Equations (28) and (31) apply to any linear differential equation for the density operator which is first order in time. Thus [38]

$$\begin{array}{c}\hfill \dot{\rho }=-i[\widehat{H},\rho ]+\sum _{k=1}^{N^2-1}{\tilde{\eta }}_k\left(2{\widehat{A}}_k\rho {\widehat{A}}_k^{+}-{\widehat{A}}_k^{+}{\widehat{A}}_k\rho -\rho {\widehat{A}}_k^{+}{\widehat{A}}_k\right),\end{array}$$

(32)

where the coefficients ${\tilde{\eta }}_k$ are given by ${\tilde{\eta }}_k={\eta }_k/2$. From (32) and (31), we have that

$$\begin{array}{c}\hfill |v\rangle =-\frac{i}{\hslash }\widehat{H}|\psi \rangle +\sum _{k=1}^{N^2-1}{\tilde{\eta }}_k\left(2\langle {\widehat{A}}_k^{+}\rangle {\widehat{A}}_k-{\widehat{A}}_k^{+}{\widehat{A}}_k-\langle {\widehat{A}}_k^{+}\rangle \langle {\widehat{A}}_k\rangle \right)|\psi \rangle .\end{array}$$

(33)

where $\langle {\widehat{A}}_k\rangle =\langle \psi |{\widehat{A}}_k|\psi \rangle$ represent the mean value of the Lindblad operator ${\widehat{A}}_k$. By replacing (32) in (28), the stochastic terms are given by

$$\begin{array}{c}\hfill |u_k\rangle =\left({\widehat{A}}_k-\langle {\widehat{A}}_k\rangle \right)|\psi \rangle .\end{array}$$

(34)

Now, from (33) and (34), the resulting stochastic differential equation for the state vector (22) is -4.6cm0cm

$$\begin{array}{c}\hfill d|\psi \rangle =-\frac{i}{\hslash }\tilde{H}|\psi \rangle +\sum _{k=1}^{N^2-1}{\tilde{\eta }}_k\left(2\langle {\widehat{A}}_k^{+}\rangle {\widehat{A}}_k-{\widehat{A}}_k^{+}{\widehat{A}}_k-\langle {\widehat{A}}_k^{+}\rangle \langle {\widehat{A}}_k\rangle \right)|\psi \rangle dt+\left({\widehat{A}}_k-\langle {\widehat{A}}_k\rangle \right)|\psi \rangle d\xi .\end{array}$$

(35)

3. Results and Discussion

As discussed above, the ISF is a key observable in the surface diffusion context. If the surface coverage $\theta$ is low, the adsorbates can be considered as independent particles, that is, the adsorbate–adsorbate interaction can be safely neglected. Thus, the ISF now reads as

$$\begin{array}{c}\hfill I(\Delta \mathbf{K},t)=\langle e^{-i\Delta \mathbf{K}·\widehat{\mathbf{R}}\left(0\right)}{e}^{i\Delta \mathbf{K}·\widehat{\mathbf{R}}\left(t\right)}\rangle .\end{array}$$

(36)

For our problem, $\Delta \mathbf{K}=\Delta k_x\overrightarrow{i}$ is the wave vector transfer parallel to the surface, $\widehat{\mathbf{R}}\left(t\right)=\widehat{x}\left(t\right)\overrightarrow{i}$ and $\widehat{\mathbf{R}}\left(0\right)=\widehat{x}\left(0\right)\overrightarrow{i}$ are the position operators of the adatom/adsorbate in the unidimensional space and the initial position operator, respectively. The x-dimension can also be considered as one of the symmetry directions of the surface, where the projection on the wave vector transfer is carried out.

The stochastic wave functions ${\psi }_i\left(t\right)$ are obtained by numerically solving (35) and using the splitting operator method together with the Fast Fourier transform (FFT). The terms with only $\widehat{x}$ or $\widehat{p}$ operators are propagated in the coordinate or momentum space, respectively. All the initial conditions for the position are chosen to be the same, $x_0=0$, and the propagation of the wave function moves towards negative values of x since the initial velocity of the adsorbate is assumed negative, $v_0<0$. Based on previous works [8,18,29], the initial momentum $p_0$ satisfies the Boltzmann relation: $p_0=\sqrt{2m{k}_{B}T}$, where T is the surface temperature and m is the adsorbate mass (Xe). The initial wave function $\psi (x,0)$ is a Gaussian function given by $\psi (x,0)=\sqrt[4]{\frac{m\omega }{\hslash \pi }}exp\left[-\frac{m\omega }{2\hslash }{(x-x_0)}^2+\frac{i}{\hslash }{p}_0(x-x_0)\right]$ where $\omega$ is the corresponding width.

The one dimensional discrete stochastic wave function is evaluated at N equally spaced points belonging to the interval $[x_1,x_N]$, with N being a power of 2 in order to use the FFT method. The discrete wave function is given at each time t by

$$\begin{array}{c}\hfill \psi (x,t)\equiv \{\psi (x_1,t);\psi (x_2,t);\psi (x_3,t);...;\psi (x_N,t)\}.\end{array}$$

(37)

The normalization of the numerical wave function is calculated as

$$\begin{array}{c}\hfill F_{norm}=\sqrt{{\sum }_{i=1}^N{\psi }^{*}(x_i,t)\psi (x_i,t)\Delta x}\end{array}$$

(38)

where $\Delta x$ is the length between two consecutive points of the grid. The normalized stochastic wave function is thus calculated as

$$\begin{array}{c}\hfill {\psi }_n(x_i,t)=\frac{\psi (x,t)}{F_{norm}}.\end{array}$$

(39)

The mean value of any operator $\widehat{B}$ can be expressed as

$$\begin{array}{c}\hfill <\widehat{B}>=\sum _{i=1}^N{\psi }_n^{*}\left(x_i\right)\widehat{B}{\psi }_n\left(x_i\right)\Delta x.\end{array}$$

(40)

From the definition of scalar product, the ISF (36) is now written as

$$\begin{array}{c}\hfill I(\Delta k_x,t)=\frac{1}{N_s}\sum _{j=1}^{N_s}\langle {\psi }_{n_j}(x,t)|e^{-i\Delta k_x[\widehat{x}\left(t\right)-\widehat{x}\left(0\right)]}|{\psi }_{n_j}(x,t)\rangle .\end{array}$$

(41)

where the mean value of the operator $e^{-i\Delta k_x[\widehat{x}\left(t\right)-\widehat{x}\left(0\right)]}$ and $\widehat{x}$ can be obtained from (40). Thus, the ISF $I(\Delta k_x,t)$ is then calculated according to

$$\begin{array}{cc}\hfill I(\Delta k_x,t)& =\frac{1}{N_s}\sum _{j=1}^{N_s}\left(\sum _{i=1}^N\left({\psi }_{n_j}^{*}(x_i,t)e^{-i\Delta k_x[x_i\left(t\right)-x\left(0\right)]}{\psi }_{n_j}(x_i,t)\Delta x\right)\right)\hfill \\ & =\frac{1}{N_s}\sum _{j=1}^{N_s}{I}_j(\Delta k_x,t),\hfill \end{array}$$

(42)

where $I_j(\Delta k_x,t)$ is the ISF for each stochastic realization. These values are calculated independently and summed up to have the average value $I(\Delta k_x,t)$. The stochastic differential Equation (35) is solved $N_s$ times. For a number of realizations of $N_s=1000$, the numerical stability for the mean value $I(\Delta k_x,t)$ is achieved. In particular, the unidimensional space x is chosen to be in the interval $x\in [x_i,x_f]$, where $x_i=-150\AA$ and $x_f=30\AA$. The number of points chosen is $N=8192$. The increment $\Delta x$ can be calculated from $\Delta x=(x_f-x_i)/N$. For these parameters, the value obtained for this increment is $\Delta x=0.02\AA$. For simplicity, from now on, $\Delta k_x$ is replaced by $\Delta k$

3.1. Flat Surface

It is important to first analyze the diffusion on a flat surface, $V=0$. This example is representative of low corrugated surfaces where the role of the activation barrier is negligible since the adsorbate thermal energy is much greater than the barrier. The adsorbate-surface interaction is through the friction coefficient $\eta$. In the limiting case of small friction and flat potential, the results of the simulations have been compared with the case of a two-dimensional (2D) ideal gas. The ISF (36) for the 2D ideal gas is a Gaussian function

$$\begin{array}{c}\hfill {\displaystyle I(\Delta k,t)=e^{-\frac{\Delta k^2{v}_0^2{t}^2}{4}},}\end{array}$$

(43)

where $v_0=\sqrt{2k_{B}T/m}$ is the initial thermal adsorbate velocity.

In Figure 1, and in order to first check our propagation method, $I(\Delta k_x,t)$ is plotted versus the propagation time for the surface parallel wave vector transfer $\Delta k=0.40{\AA }^{-1}$, friction coefficient $\eta =0$, flat potential $V=0$, surface temperature $T=105$ K and initial adsorbate velocity $v_0=115$ m/s with negative direction. The propagation time t of the stochastic wave function ${\psi }_i\left(t\right)$ and the time step $\Delta t$ are chosen to be equal to $t=14$ ps and $\Delta t=0.01$ ps, respectively. These values correspond to a number of time steps of $N_t=1400$. In Figure 1, our numerical simulation (SN) in red is also plotted showing a very good agreement with the analytic function (43), valid for a 2D ideal gas. This result indicates that the Xe-Pt(111) system under these conditions behave as a true 2D ideal gas where the diffusion of the adsorbate Xe is completely free on the Pt(111) surface. This behavior observed in Figure 1 is also reproduced for different $\Delta k$ values in the interval $[0.05,2.5]{\AA }^{-1}$ analyzed. This free particle motion is kept at very long times.

The next step is now to study the behavior of this system in the ballistic regime at different friction coefficients for a flat surface and compare with experimental results and previous standard Langevin simulations [20]. As mentioned above, this regime is only valid at times $t<<{\eta }^{-1}$. For such a goal, our purpose is to see how good the Gaussian function, which describes the ballistic regime, is. The Fourier transform in time of the Gaussian ISF is again a Gaussian function which corresponds to the DSF. The full width at half maximum (FWHM) of the SDF, $\Gamma$, is extracted and plotted at different frictions for different momentum transfer values $\Delta k\in [0.05,2.5]$${\AA }^{-1}$. The friction coefficients $\eta$ used in these simulations are equal to the values used in a previous work [20]. By fitting the ISF to a Gaussian function with width ${\sigma }_t$, $exp(-t^2/2{\sigma }_t)$, the values of $\Gamma$ can be obtained from ${\sigma }_t$ through

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\frac{2\hslash \sqrt{2ln2}}{{\sigma }_t}.\end{array}$$

(44)

The width in time space ${\sigma }_t$ of the ISF for a 2D ideal gas (43) is ${\sigma }_t=\sqrt{2}/v_0\Delta k$. Using (44), the dependence of $\mathrm{\Gamma }$ with the parallel wave vector transfer $\Delta k$ is given by

$$\begin{array}{c}\hfill \mathrm{\Gamma }=2\sqrt{ln2}\hslash v_0\Delta k.\end{array}$$

(45)

Figure 2 shows a comparison between the numerical results of $\mathrm{\Gamma }$ as a function of the parallel wave vector transfer $\Delta k$ at $T=105$ K for two different values of the friction coefficient, $\eta =0.25,2.0$ ps${}^{-1}$, and by assuming the corrugation of the Pt(111) surface to be negligible. Black lines correspond to Equation (45). Left top panel (a) displays our numerical simulations. Langevin molecular dynamics simulations are displayed on the right top panel (b) [20] for the same conditions. In panel (c), the fitting to a Gaussian function of the ISF at very short times and small friction is also plotted. The numerical results are in good agreement with those obtained in Ref. [20]. The slope of $\mathrm{\Gamma }$ decreases as the friction coefficient increases. For a fixed value of $\Delta k$, the value of $\mathrm{\Gamma }$ is higher for $\eta =0.25$ ps${}^{-1}$, so the width in the time domain is smaller and the function $I(\Delta k,t)$ decays faster. Thus, at higher values of $\Delta k$, the fitting has to be carried out at smaller times. In panel (c), the ballistic regime ($t<<4$ ps) seems to be well described around 0.8 ps. For $\eta =2.00$ ps${}^{-1}$, the friction is high enough to drastically reduce the ballistic regime at very, very short times. In this case, the interaction with the substrate is very important and the diffusive regime is reached at times $t>>{\eta }^{-1}$. This regime is not studied here.

3.2. Corrugated Potential

The potential energy surface (PES) is simulated using 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\AA$. To ensure that the metal surface is large enough so that the boundary effects do not influence the interaction potential, the number of 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}^{-3}$ meV.

A cosine corrugation function is proposed to represent the potential energy surface of the Xe-Pt(111) system given by $V_a(x,y)=V_{1}cos\left(\frac{2\pi }{a}x\right)cos\left(\frac{2\pi }{a}y\right)+V_2$, where $V_1=-13.64$ meV and $V_2=-12.32$ meV. These values are obtained by fitting the cosine function to the PES numerical result. Both numerical and analytic PSE are very similar because the amplitude of the interaction potential takes an approximated value equal to $-28$ meV and the periodicity match with the unit cell length of the Pt for both cases. The unidimensional corrugated potential is then used in the simulations with the form $V\left(x\right)=V_a(x,0)=V_{0}cos\left(\frac{2\pi }{a}x\right)$. The amplitude $V_0$ is taken so that the energy barrier for this potential ($E_b=V_0$) coincides with the Langevin molecular simulations [20]. In Figure 3, the directions (010) and (100) of the metallic surface of Pt are plotted for the corrugated potential $V_a(0,y)$ and $V_a(x,0)$, respectively, in the interval $x\in [90,110]\AA$. The shape of the potential cutoff along the (010) direction is analogous to the (100). Thus, the diffusion of the Xe adsorbate on the metal surface of Pt is studied using the (100) direction for the interaction potential.

In Figure 4, the influence of surface corrugation on the Xe-Pt(111) system for a simulation with small friction coefficient $\eta =0.03$ ps${}^{-1}$ and surface temperature $T=105$ K is analyzed. Figure 4 shows the quasielastic broadening $\mathrm{\Gamma }$ as a function of $\Delta k$. The slope of $\mathrm{\Gamma }$ decreases as $E_b$ increases. The adsorbate thermal energy can be obtained from $E_T=k_{B}T$. For a temperature of $T=105$ K, the thermal energy is $E_T=9.04$ meV. If the energy barrier is less than $9.04$ meV, the adsorbate can diffuse more or less freely over the metal surface of Pt. For energy barriers belonging to this range, the system is in a ballistic regime at very short times. If $E_b=9.6$ meV, the difference between both energies is small and the system is at the limit where some Xe gas particles can be trapped, but others can diffuse over the surface. However, if $E_b>9.6$ meV the thermal energy is not enough to overcome the barrier and almost all the adsorbate particles are trapped in the potential valley and there is no diffusion. In panel (c), the propagation time is much smaller than 1 ps in order to be in the ballisitc regime, $t<<{\eta }^{-1}$. In this case, the chosen energy barrier is $E_b=9.6$ meV and the fitting to a Gaussian function is worse. The system under these conditions behaves as a free 2D gas within a unit cell for very short times, particles are trapped inside the potential valley. If a lower barrier energy is chosen, the corresponding fitting is better and the propagation time in this regime is longer. Our numerical results are in good agreement with those obtained by Langevin molecular simulations [20] and experimental data marked by triangles and squares at low $E_b$ values.

4. Conclusions

The stochastic wave function formalism has been briefly reviewed and applied to the diffusion of a Xe gas on a Pt(111) surface. The Lindblad master equation, when it is solved by means of stochastic wave functions, provides a good physical description of the ballistic regime in this diffusion process. It provides an efficient numerical scheme to study the quantum dynamics of open systems. As far as we know, this is the first time this theoretical formalism has been applied to the surface diffusion context, since most of works use the Langevin formalism (in its generalized or standard version). This is particularly encouraging since recent experiments with better spatial and time resolutions (HeSE technique) tend to reveal the underlying properties of individual systems and their jump events directly in the time domain. New results for low and moderate values of the surface coverage for the Xe/Pt(111) system have been recently obtained. Thus, thanks to this type of study with quite satisfactory results, the next step will be to extend this theoretical formalism to deal with the diffusion regime, $t>>{\eta }^{-1}$. This would allow us to consider more physical systems where recent experimental results are waiting for a better interpretation. Adsorbates, such as $Na$, $Li$, $H_2$, etc., on metallic surfaces have been studied by using the HeSE technique. For light adsorbates, quantum methods are more convenient than classical ones. In the surface diffusion context, this extension, which would be carried out for the first time, is again very promising.

Perspectives

The dependence with the surface temperature and coverage can be studied for a better understanding of microscopic processes. It is suggested to also extend and improve this numerical method for times greater than 20 ps with periodic two-dimensional interaction potential in order to reach the Brownian regime. This method turns out to be a valuable alternative for dealing with more diffusive systems studied via QHAS or HeSE such as, for example, Na-Cu(001), Ar-MgO(100), etc., which have been previously considered within the Langevin formalism [14,18,40].