Plasmon effects in photoemission

We develop the concept of scattering matrix and we use it to perform stable numerical calculations of photo-emission from nano-tips. Electrons move in an external space and time dependent nonperturbative electric field. We apply our algorithm for different strengths and spatial configurations of the field.


I. INTRODUCTION
The aim of this paper is to investigate some particular quantum processes taking place in an arbitrary space-dependent scalar potential and a time-and spacedependent vector potential. Vector potential is periodic in time and describes a laser field. Its space-dependence results from the interaction of the laser field with electrons in solids. Such conditions are met for example in semiconductor nanostructures [1][2][3] (like quantum wires or wells), photoemission from a metal tip [4,5], carbon nanotubes or graphene [6][7][8] or in surface physics [9][10][11][12][13][14][15]. To make our presentation as clear as possible we shall restrict ourselves to the one-space-dimensional case, although extension of the presented method to systems of higher dimensionality is possible (see, e.g. [16]). We shall apply our method to investigation of the photo-emission process.
This paper is organized as follows. In Sec. II the most general solution of the Schrödinger equation is introduced. The transfer-matrix method and matching conditions are analyzed in Sec. III, whereas reflection and transition probabilities are introduced in Sec. IV. These probabilities must sum up to 1, which puts a very strong check for the accuracy of numerical calculations. The most important part of this paper, i.e. the concept of the scattering-matrix method, is discussed in the next section, where it is shown why the scattering-matrix algorithm has to be introduced, instead of a much simpler transfer-matrix algorithm. Numerical illustrations of the applicability of this algorithm are presented in Sec. VI, and are followed by short conclusions.
In our numerical illustrations we use atomic units, if otherwise stated.

II. SOLUTION OF THE SCHRÖDINGER EQUATION
Let us start with one-dimensional Schrödinger equation of the form [17], i∂ t ψ(x, t) = 1 2 +V (x) ψ(x, t). (1) Space-dependent mass m(x), scalar potential V (x) and vector potential A(x, t) are spatially constant in finite intervals. Their values in any interval (x i−1 , x i ) will be denoted as m i , V i and A i (t). We require also that the function A(x, t) is periodic in time, that is where T = 2π/ω and ω is the frequency of the oscillating in time electric field. Defining in a standard way the probability density ρ(x, t), and the probability current j(x, t), we show using Eq. (1) that the conservation of probability condition is satisfied. Indeed, assuming the above definitions, we get the continuity equation, Space dependence of mass in Eq. (1) forces one to impose non-standard continuity conditions on any solution of this equation. It is now the wavefunction ψ(x, t) and the quantity that have to be continuous at points of discontinuity of mass m(x) and both potentials V (x) and A(x, t) [17][18][19][20]. Before passing to a general solution ψ(x, t) of Eq. (1) in any given interval (x i−1 , x i ), which we shall denote as ψ i (x, t), let us note that due to time periodicity of the Hamiltonian, ψ i (x, t) can be chosen such that the Floquet condition, is satisfied, where E is the so-called quasienergy. A general solution ψ i (x, t) of Eq. (1) in any interval (x i−1 , x i ) takes then the following form [22,23], where C σ iN are arbitrary complex numbers to be determined and with U i = e 2 A 2 i (t) /2m i being the ponderomotive energy, where A 2 i (t) means the time-average of A 2 i (t) over the laser-field oscillation. Components for which p iN are purely imaginary are called closed channels. These channels are not observed for a particle in initial or final states, but they have to be taken into account in order to satisfy the unitary condition of the time evolution. In a general case, the B M−N (σp iN ) functions are components of the following Fourier expansion, provided that the vector potential A(x, t) is periodic in time. Functions Φ σ iN (t) are defined as follows: It is easily seen from the above equation that the B M−N (σp iN ) functions depend on the form of the vector potential A(x, t), that is on the laser field applied.

III. MATCHING CONDITIONS AND TRANSFER MATRIX
Continuity conditions discussed above and applied to a general solution (8) of the Schrödinger equation (1) lead to an infinite chain of equations connecting constants C σ iN in the neighboring domains. These matching conditions can be written in the matrix form, where C ± iN = [C ± i ] N are the components of the columns C ± i . The matrices B(i, x) and C i are defined as follows, (13) The elements of B(i, x) can be computed in the following way.
For an arbitrary function A(x, t), periodic in time with the period T we have where ω = 2π/T . In the interval (x i−1 , x i ) coefficients b n (x) assume constant values, which we shall denote as b i,n . Using the condition of the continuity of the wavefunction ψ i (x, t) at the point x i−1 , we compute the elements of the matrices B + and B − , On the other hand elements of the B ′ matrix can be evaluated by substituting a general solution (8) to the expression (6) and applying the continuity condition to it at x i−1 . After some algebraic manipulations we obtain finally the expression for the B ′ -matrices, and a set of equations for vectors C i , where These relations allow to connect a solution in a given domain where T ji is the so-called transfer matrix [19,21,22].

IV. REFLECTION AND TRANSITION PROBABILITIES
It is clear now that on the basis of Eq.(19) we can connect solutions in the boundary domains (−∞, x 0 ) and (x L−1 , ∞). Values of mass m(x), scalar potential V (x) and vector potential A(x, t) in these domains will be denoted as m 0 , V 0 , A 0 (t) and m L , V L , A L (t), respectively. We can then write down solutions of (1) for each of these domains. These solutions represent incident (ψ inc ), reflected (ψ ref ) and transmitted (ψ tr ) waves, and take the following form, where Constants C − 0,N and C + L,N will be denoted from now on as R N and T N , respectively. Using continuity conditions for functions defined above, we get the probability conservation equation for reflection and transition amplitudes, R N and T N , where summations are over such N for which p N and q N are real, i.e., over the open channels. This equation permits us to interpret and as reflection and transition probabilities for a tunneling process in which absorption (N > 0) or emission (N < 0) of energy N ω by electrons occurred [20,22]. In case of a monochromatic laser field this process can be interpreted as absorption or emission of N photons from the laser field. The unitary condition (24) can be also interpreted as the conservation of electric charge. To this end, let us define the quantities proportional to the density of electric currents, Then Eq. (24) adopts the form of the first Kirchhoff low, Using (19) we can calculate constants C − 0,N = R N and C + L,N = T N appearing in equations (20) - (22). Indeed, since where transfer matrix T = T L0 , and because T , C 0 and C L adopt the following block forms, we arrive at where R and T denote the columns of R N i T N , and [C + 0 ] N = δ 0,N . Thus, after some algebraic manipulations, we have, which allows us to determine the quantities R N and T N for a given transfer matrix T . For open channels, these quantities are the amplitudes of reflection (R N ) and transition (T N ) probabilities, from which one can compute reflection and transition probabilities using equations (25) and (26).

V. THE SCATTERING MATRIX
We note from equations (15) and (16) that each of the B i matrices that constitute the transfer matrix T ji contain elements exp(±ip i,N x i ) that depend on the x i coordinates at which the discontinuities appear. For closed channels, that is when the p i,N momenta are purely imaginary, these numbers are real and may assume arbitrary values, very large or very small, depending again on the x i coordinates. Number of the B i matrices is equal to the number of discontinuity points, that is it depends on how we divide the space into short intervals in order to make our potential tractable by our algorithm. It may therefore turn out that in order to compute the transfer matrix T ji , we have to multiply a large number of the B i matrices, each containing both very small and very large numbers. It is clear that such a procedure is numerically unstable. We have to find a way to modify our method of calculations in order to compute the elements of each B i matrix at the same point x = 0 independently of where the 'real' x i is. This would eliminate "dangerous" exp(±ip i,N x i ) elements (turning them to 1), however at the cost of appearing somewhere else. We shall see later that these 'left-overs' of the shift into x = 0 appear only as differences x i+1 − x i and therefore do not cause any harmful side-effects. We shall see now that such a modification is possible and the price we pay for it is worth the effort.
It follows from Eq. (19) that in the neighboring domains, (x i−2 , x i−1 ) and (x i−1 , x i ), we have,  (iσp iN δ). These constants can be included in coefficients C σ iN . In this way we get a new set of constants which we shall denote asC σ iN , We shall interpret these constants as coefficients in solution (8), given by the continuity conditions at point Eq. (36) written in the matrix form becomes, where and In the equation above P σ i (δ) is a diagonal matrix, whereas iN . It follows from the form of the matrix P i (δ) that the following relations are satisfied: Let us notice also that translation of the system defined above modifies the transfer matrix T i,i−1 . We have and we can write it down as Matrix elements denoted with the tilde symbol refer to the translated system. Using the method defined above and the relation (19), we can connect now the solution in the domain (−∞, x 0 ) with the solution in any other domain (x i−1 , x i ). In this way the elements of the transfer matrix, which have been computed until now at the points of discontinuity x 0 . . . x i−1 , are computed now each time at the same point x = 0. Let us illustrate this method for a special case of i = 3 Equation (47) connects constants C 0 and C 3 using the matrices T 0 j,j−1 all computed at x = 0 independently of j, and diagonal matrices P j (δ j ), given by the relations (38) and (40), where δ j = (x j − x j−1 ). Edge matrices P 0 (x 0 ) and P −1 3 (x 2 ) in the equation (47) can be omitted while computing the transmission and reflection probability amplitudes since their only role is to multiply the amplitudes by phase quotients which disappear while computing the probabilities. Although these matrices lead to significant modifications of the closed channels in the domains of x < x 0 and x > x 3 in this particular case, these channels do not influence the reflection and transition amplitudes. Transmission and reflection probabilities can thus be computed using a modified transfer matrix, The matrices T 0 i,i−1 are equal to the matrices B i in Eq. (18) calculated however for x i−1 = 0. This fact speeds up numerical calculations since now matrix B(i, x = 0) in Eq. (18) have to be inverted only once. Further on we shall omit the superscript 0 in T and the tilde over C in order to simplify notation.
The method presented above is still numerically unstable. The reason for this instability lies in the existence of large numerical values of elements of the P − i (δ) matrix for imaginary momenta p iN . In other words, for the source of numerical instabilities are matrix elements T −− i,i−1 that contain large numbers. There is however a chance for improving the stability, if only its reverse will be used, T −− i,i−1 −1 . It appears that it is possible provided that in our numerical algorithm only the so-called scattering matrix will be applied. For this reason we will show below how to compute the scattering matrix, S j,i , using only elements of the transfer matrix, T j,i . For the transfer matrix T j,i we have, Thus, On the basis of (51) we now want to compute the elements of the S j,i matrix. This matrix is supposed to connect the coefficients C ± i and C ± j in the following way, (52) Using the set of linear equations (51), we easily compute the coefficients C − i and C + j on the left-hand side of equation (52), as functions of the coefficients C − j and C + i . We get then the following relations, Finally we compute the elements of the matrix S j,i , As expected, the matrix S j,i contains only numerically stable elements (T −− j,i ) −1 . It follows from Eq. (19) that the transfer matrix T j,i can be written as the product of two transfer matrices, T j,k and T k,i (i < k < j), where matrices T j,k and T k,i are defined as follows, Applying the method presented above, for each of the transfer matrices T j,k and T k,i we can construct a scattering matrix, S j,k and S k,i respectively. Elements of the scattering matrix S j,i can be computed using only elements of the matrices S j,k and S k,i . Using the notation above, we obtain the following expressions for the elements of the S j,i matrix, It is clear from the above that the S j,i matrix is not merely a product of two matrices S j,k and S k,i , but rather a complicated nonlinear composition of them. It is important however to note that despite its evident complexity, such a construction of the scattering matrix is numerically stable, as opposed to the transfer matrix method which fails if a system with a large number of discontinuity points x i is considered. Stability of such an algorithm has been proven in our numerical investigations by checking that the condition (24) is satisfied with an error smaller than 10 −14 . Such an accuracy can never be achieved for systems with a large number of discontinuity points if the transfer matrix is applied.

VI. PHOTO-EMISSION
In our model investigations, we concentrate on some essential features of the solid-vacuum interface, as exemplified by the Sommerfeld model, in which the band structure is neglected. This simplification allows us to consider a quite general form of the laser field. To be more specific, the solid surface is described by a continuous step potential with The parameter w 0 determines the skin depth of a surface. For w 0 = 0, the surface potential represents the step function, commonly used in the Sommerfeld model. In our illustrations, we put w 0 = 5. We apply our theory to the gold surface and assume that the electron effective mass is close to the free electron mass. The work function and the Fermi energy for the gold metal are equal to 5.1 and 5.53 eV, respectively. This means that the constant V 0 above (as the sum of the work function and the Fermi energy) equals 10.63 eV. The surface potential described above can be generalized further to meet conditions suitable for other solids. In particular, one can take into account the spacedependent effective mass of electrons in semiconductor heterostructures or metals with effective masses different from the free electron mass.
On the other hand, the form of the laser field is assumed to depend on both space and time coordinates. Since, for laser pulses of duration ∼ 30 fs and the 800 nm wavelength, the monochromatic approximation works well, we therefore adopt the following form for the laser electric field: and similarly The parameter ǫ defines the plasmon-enhanced part of the laser field. For the incident laser beam, we choose the Ti:sapphire laser beam of frequency ω = 1.5498 eV (λ = 800 nm). This means that inside the solid the laser field intensity averaged over the time period decays exponentially, On the other hand, in vacuum, it stays constant close to the surface, and then again decays exponentially. In this way, we can mimic a real physical situation in which the radiation-filled space is finite. In our illustrations, we take ζ L = 40, which means that the penetration depth of the laser field intensity equals ζ L /2 = 20. The parameter a L ζ L describes the distance in a solid at which the intensity is not reduced substantially. On the other hand, b L µ L corresponds to the laser focus diameter in vacuum, whereas µ L alone determines the intensity reduction rate outside the focus. Similar parameters with the subscript P refer to the plasmon-enhanced part of the laser field. The remaining parameters have been chosen as follows: a L = 3, b L = 20, µ L = 100, a P = 1, ζ P = 8, b P = 4, µ P = 20, and ǫ = 0, . . . , 5. All dimensional parameters are in atomic units. In our discussions presented below, the laser field intensity is characterized by the dimensionless parameter ξ = U p /ω, where U p = E 2 0 /(4ω 2 ) is the ponderomotive energy of electrons in the monochromatic electromagnetic plane wave of frequency ω; hence E 0 = 2ω √ ωξ. In Fig. 1 we draw the space-dependence of the continuous step potential V (x) and the electric field amplitude E 0 (x) for ǫ = 5 and ξ = 0.1.
The total photo-emission probability is equal to We plot it in Fig. 2 as a function of the electron kinetic energy for ξ = 0.1 and for six values of ǫ. We clearly see the multi-photon structure in this distribution, i.e., the total probability jumps sometimes by a few orders of magnitude if the smaller number of laser photons is sufficient for photo-emission. As expected, the plasmon effect usually increases the photo-emission probability. Moreover, the energy of the multi-photon channel opening increases with increasing ǫ, which is due to the increase of the space-dependent ponderomotive energy of the laser field. The significance of this effect for the tunneling phenomena is going to be discussed in due course.

VII. CONCLUSIONS
As mentioned above, our algorithm is convergent provided that a sufficient number of discretization points is introduced. For systems considered here, this number should not be smaller than 100. If the laser field is very weak, this does not create significant numerical problems, except that calculations become longer. However, when the laser field is sufficiently intense, the algorithm based on the transfer matrix is unstable. This instability is due to the existence of closed channels, which introduce into numerical calculations very small and very large numbers at the same time. Augmenting precisions significantly slows down the calculation and does not diminish the problem. We have found that it is possible to make this algorithm numerically stable by just applying nonlinear matrix transformations, without introducing higher precisions.
Illustrations presented in this paper show that photoemission of electrons can be changed significantly by applied nonperturbative oscillating in time and spacedependent electric field. The efficiency of the algorithm presented in this contribution opens up the possibility of investigating surface phenomena in the presence of more realistic laser pulses that gradually decrease within solids and extend on a mesoscopic scale in vacuum.