


@ 1,283 +1,190 @@





\documentclass{ian}





\documentclass[pra,twocolumn]{revtex42}










\usepackage{graphicx}





\usepackage{array}





\usepackage{color}





\usepackage{amssymb}





\usepackage{amsmath}










\begin{document}










\hbox{}





\hfil{\bf\LARGE





Exact solution of the Schr\"odinger equation\par





\hfil for photoemission from a metal\par





}





\vfill





\title{Exact solution of the 1D timedependent Schr\"odinger equation for the emission of quasifree electrons from a flat metal surface by a laser}










\hfil{\bf\large Ovidiu Costin}\par





\hfil{\it Department of Mathematics, The Ohio State University}\par





\hfil{\tt\color{blue}\href{mailto:costin.9@osu.edu}{costin.9@osu.edu}}\par





\vskip20pt





\author{Ovidiu Costin}





\author{Rodica Costin}





\affiliation{Ohio State University,Columbus, OH 43210, USA}










\hfil{\bf\large Rodica Costin}\par





\hfil{\it Department of Mathematics, The Ohio State University}\par





\hfil{\tt\color{blue}\href{mailto:costin.10@osu.edu}{costin.10@osu.edu}}\par





\vskip20pt





\author{Ian Jauslin}





\affiliation{Princeton University, Princeton, NJ 08544, USA}










\hfil{\bf\large Ian Jauslin}\par





\hfil{\it Department of Physics, Princeton University}\par





\hfil{\tt\color{blue}\href{mailto:ijauslin@princeton.edu}{ijauslin@princeton.edu}}\par





\vskip20pt





\author{Joel L. Lebowitz}





\affiliation{Rutgers University, Piscataway, NJ 08854, USA}










\hfil{\bf\large Joel L. Lebowitz}\par





\hfil{\it Department of Mathematics and Physics, Rutgers University}\par





\hfil{\tt\color{blue}\href{mailto:lebowitz@math.rutgers.edu}{lebowitz@math.rutgers.edu}}\par





\vskip20pt










\vfill















\hfil {\bf Abstract}\par





\medskip





{\small





We solve rigorously the time dependent Schr\"odinger equation describing electron emission from a metal surface by a laser field perpendicular to the surface.





We consider the system to be onedimensional, with the halfline $x<0$ corresponding to the bulk of the metal and $x>0$ to the vacuum.





The laser field is modeled as a classical electric field oscillating with frequency $\omega$, acting only at $x>0$.





We consider an initial condition which is a stationary state of the system without a field, and, at time $t=0$, the field is switched on.





We prove the existence of a solution $\psi(x,t)$ of the Schr\"odinger equation for $t>0$, and compute the surface current.





The current exhibits a complex oscillatory behavior, which is not captured by the ``simple'' three step scenario.





As $t\to\infty$, $\psi(x,t)$ converges with a rate $t^{\frac32}$ to a time periodic function with period $\frac{2\pi}{\omega}$ which coincides with that found by Faisal, Kami\'nski and Saczuk (Phys Rev A {\bf 72}, 023412, 2015).





However, for realistic values of the parameters, we have found that it can take quite a long time (over 50 laser periods) for the system to converge to its asymptote.





Of particular physical importance is the current averaged over a laser period $\frac{2\pi}\omega$, which exhibits a dramatic increase when $\hbar\omega$ becomes larger than the work function of the metal, which is consistent with the original photoelectric effect.





}










\vfill










\tableofcontents










\vfill





\eject










\setcounter{page}1





\pagestyle{plain}





\begin{abstract}





We solve exactly the onedimensional Schr\"odinger equation for $\psi(x,t)$ describing the emission of electrons from a flat metal surface, located at $x=0$, by a periodic electric field $E\cos(\omega t)$ at $x>0$, turned on at $t=0$.





We prove that for all physical initial conditions $\psi(x,0)$, the solution $\psi(x,t)$ exists, and converges for long times, at a rate $t^{\frac32}$, to a periodic solution considered by Faisal et al. (Phys. Rev. A {\bf 72}, 023412 (2005)).





Using the exact solution, we compute $\psi(x,t)$, for $t>0$, via an exponentially convergent algorithm, taking as an initial condition a generalized eigenfunction representing a stationary state for $E=0$.





We find, among other things, that: (i) the time it takes the current to reach its asymptotic state may be large compared to the period of the laser;





(ii) the current averaged over a period increases dramatically as $\hbar\omega$ becomes larger than the work function of the metal plus the ponderomotive energy in the field. For weak fields, the latter is negligible, and the transition is at the same frequency as in the Einstein photoelectric effect;





(iii) the current at the interface exhibits a complex oscillatory behavior, with the number of oscillations in one period increasing with the laser intensity and period.





These oscillations get damped strongly as $x$ increases.





\end{abstract}










\maketitle










\section{Introduction}





\indent In this note, we continue our investigation of the exact solution of the timedependent Schr\"odinger equation describing the emission of electrons from a metal surface by an electric field \cite{CCe18}.





We use the Sommerfeld model of quasifree electrons with a Fermi distribution of energies, confined by a step potential $U=\mathcal E_F+W$, where $\mathcal E_F$ is the Fermi energy and $W$ is the work function of the metal.





This setup is commonly used as a simple model for the process of emission, both for a constant field and an oscillating field \cite{FN28,YGR11,Fo16,Je17}.





In both cases one focuses attention on electrons, part of the Fermi sea, with energy $\frac12k^2$ (in atomic units), moving from the left towards the metal surface at $x=0$.





These are described by a wave function $e^{ikx}$, $k>0$, $x<0$.





When the electrons reach the metal surface they are either reflected or transmitted through the barrier.





\bigskip





\indent There have been many advances in recent years in the development and application of short intense laser pulses to produce femtosecond and even attosecond beams of electrons from metallic surfaces \cite{HKK06,SKH10,BGe10,KSH11,KSe12,THH12,HSe12,PPe12,HBe12,PSe14,HWR14,EHe15,BBe15,YSe16,FSe16,RLe16,FPe16,LJ16,HKe17,SSe17,PHe17,Je17,WKe17,KLe18,LCe18,SMe18}.





A full microscopic description of the shorttime behavior of the emission process is therefore highly desirable.










\indent The time evolution of the wave function of an electron in such a beam is described by the one dimensional Schr\"odinger equation,





\indent In this note, we present, for the first time, an exact solution for the timedependent Schr\"odinger equation describing the emission of electrons from a flat metal surface by an oscillating electric field.





We use the Sommerfeld model of quasifree electrons with a Fermi distribution of energies, confined by a step potential $U=\mathcal E_F+W$, where $\mathcal E_F$ is the Fermi energy and $W$ is the work function of the metal.





This setup was first used by Fowler and Nordheim \cite{FN28} in 1928 for a timeindependent field, and is commonly used as a model for the process of emission, both for a constant and an oscillating field \cite{FN28,FKS05,HKK06,KSH11,YGR11,KSe12b,PA12,YHe13,CPe14,ZL16,Fo16,Je17,KLe18}.





In both cases one imagines the metal occupies the half space $x<0$, and focuses attention on electrons, part of the Fermi sea, moving from the left towards the metal surface at $x=0$.





These are described by a wave function $e^{ikx}$, $k>0$, $x<0$ and have energy $\frac12k^2$ (in atomic units).





In the sequel, we shall generally consider values of $k$ such that $\frac12k^2=\mathcal E_F$.





The field is described classically.





\medskip










\indent The time evolution of the wave function of an electron in such a beam subjected to an oscillating field for $x\geqslant 0$, is described by the one dimensional Schr\"odinger equation: for $x\in\mathbb R$ and $t>0$,





\begin{equation}





i\partial_t\psi(x,t)=H(x,t)\psi(x,t)=\frac12\Delta\psi(x,t)+\Theta(x)(UEx\cos(\omega t))\psi(x,t)





i\partial_t\psi(x,t)=\frac12\Delta\psi(x,t)+\Theta(x)(UEx\cos(\omega t))\psi(x,t)





\label{schrodinger}





\end{equation}





where $\Theta(x)=0$ for $x<0$ and $\Theta(x)=1$ for $x>0$, $E$ is the electric field perpendicular to the surface, $\frac\omega{2\pi}$ is the frequency and we are using atomic units $\hbar=m=1$.





Here, we ignore the carrier wave and Shottky effect.





We also only focus on the direction that is orthogonal to the metal surface, ignoring the other two directions \cite{KSH11}.





We note that in experiments one usually applies the laser field to a sharp tip in order to enhance the strength of the field.





One also includes a carrier wave envelope.





Here we ignore these as well as the Shottky effect.





Including them would greatly complicate the problem.





We believe that the simpler model considered here already captures many of the relevant physical phenomena so we focus on its exact solution.





The values of the field we use in our computations are those generally used for the enhanced field at a sharp tip.





The short time behavior would be the same as if the field was cut off after some time $t_0$.










\indent In the absence of an external field, $E=0$, the Schr\"odinger equation (\ref{schrodinger}) has a ``stationary" solution $e^{i\frac12k^2t}\varphi_0(x)$ in which there is, for $k^2<2U$, a reflected beam of the same energy and intensity as the incoming beam $e^{ikx}$ and an evanescent, exponentially decaying tail on the right,





\indent In the absence of an external field, $E=0$, the Schr\"odinger equation (\ref{schrodinger}) has a ``stationary" solution $e^{i\frac12k^2t}\varphi_0(x)$ in which there is, for $k^2<2U$, a reflected beam of the same energy and intensity as the incoming beam $e^{ikx}$ and an evanescent, exponentially decaying tail on the right.





The requirement of continuity of $\psi$ and its spatial derivative at $x=0$ then gives \cite{Je17}





\begin{equation} \label{initial}





\varphi_0(x)=





\left\{\begin{array}{ll}





e^{ikx}+R_0e^{ikx}&\mathrm{for\ }x<0\\





T_0e^{ \alpha x}&\mathrm{for\ }x>0





\left\{\begin{array}{>\displaystyle ll}





e^{ikx}+\frac{ik+\sqrt{2Uk^2}}{ik\sqrt{2Uk^2}}e^{ikx}&\mathrm{for\ }x<0





\\[0.5cm]





\frac{2ik}{ik\sqrt{2Uk^2}}e^{\sqrt{2Uk^2} x}&\mathrm{for\ }x>0





\end{array}\right.





\end{equation}





Using the continuity of $\psi$ and its derivative at $x=0$ yields





\begin{equation}





\alpha:=\sqrt{2Uk^2}





,\quad





1+R_0=T_0=\frac{2ik}{ik\sqrt{2Uk^2}}





.





\end{equation}





The current,





\begin{equation}





j(x,t):=\mathcal Im(\psi^*(x,t)\partial_x\psi(x,t))





\label{current}





\end{equation}





is zero and no electrons leave the metal.





Note that $\varphi_0$ is not square integrable, which will complicate matters when it comes to solving the Schr\"odinger equation.





\medskip










\indent If we turn on a constant field, $E>0, \ \omega=0$, at $t=0$, the Schr\"odinger equation has a ``stationary'' solution ${\rm e}^{\frac 12 k^2 t}\varphi_{E}(x)$ with $\varphi_E$ satisfying the equation





\indent In \cite{CCe18} we solved the timedependent Schr\"odinger equation (\ref{schrodinger}) for any initial condition, including $\varphi_0(x)$, for the constant field.





This corresponds to setting $\omega=0$ in (\ref{schrodinger}).





We showed that $\psi(x,t)$ converges, as $t\to\infty$, to the well known FowlerNordheim (FN) solution \cite{FN28} for emission by a constant field.





FN assumed a solution of (\ref{schrodinger}), with $\omega=0$, of the form $e^{\frac12k^2t}\varphi_E(x)$, so that $\varphi_E(x)$ satisfies the equation





\begin{equation}





\frac12k^2\varphi_E=\frac12\frac{\partial^2}{\partial x^2}\varphi_E+\Theta(x)(UEx)\varphi_E





\frac12\Delta\varphi_E+\Theta(x)(UEx)\varphi_E=\frac12k^2\varphi_E





.





\end{equation}





with incident beam $\Theta(x)e^{ikx}$.





This equation was solved by Fowler and Nordheim \cite{FN28} in 1928 to describe the quantum tunneling of electrons through a triangular barrier.





Their solution, with some corrections, is commonly used to describe the stationary current produced by a field $E$ acting on the surface of a metal \cite{Je17}.





\bigskip





The solution $\varphi_E(x)$ has the form $\varphi_E(x)=e^{ikx}+R_Ee^{ikx}$ for $x<0$ and an Airy function expression for $x>0$.





The FN computation of the tunneling current from $\varphi_E(x)$ via (\ref{current}), is still the basic ingredient for the analysis of experiments at present \cite{Je17}.





The main modification is the use of the Schottky factor \cite{Je17,Fo16}, rounding off the barrier at $x=0$, which, as already noted by FN, is only important for $\frac12k^2\sim U$.










\indent In \cite{CCe18} we solved the timedependent Schr\"odinger equation (\ref{schrodinger}) with initial condition $\varphi_0(x)$ for the constant field.





We showed that $\psi(x,t)$ converges, as $t\to\infty$, to the FowlerNordheim (FN) solution ${\rm e}^{\frac 12 k^2 t}\varphi_{E}(x)$





The exact rate of convergence is $t^{\frac32}$.





The ``effective'' time of approach to the FN solution was found to be, for realistic values of the parameters, of the order of femtoseconds.





The rate of convergence of the solution of the timedependent Schr\"odinger equation (\ref{schrodinger}), with $\omega=0$ and an initial condition of the form (\ref{initial}), to the FN solution was shown in \cite{CCe18} to be like $t^{\frac32}$.





Surprisingly, the deviation of the current from the FN solution becomes quickly very small, so the ``effective'' time of approach to the FN solution was found to be, for realistic values of the parameters, of the order of femtoseconds.





It is therefore not significant for emission in constant fields acting over much longer times.










\indent Here, we investigate solutions of (\ref{schrodinger}) with $E>0$, $\omega>0$, $\psi(x,0)=\varphi_0(x)$ as given in (\ref{initial}).





Since we are interested in describing situations in which the laser is turned on for only a short time, with pulses as short as a few periods, understanding the role of the initial state is important.





This problem turns out to be considerably more difficult than the constant field case, $\omega=0$.





\bigskip





\indent Here, we investigate solutions of (\ref{schrodinger}) with $E>0$, $\omega>0$ and general $\psi(x,0)$.





This covers a wide range of physical situations, depending on $\omega$ and $E$, ranging from mechanically produced oscillating fields to those produced by lasers of high frequency.





As the Keldysh parameter $\gamma:=\frac{2\omega}E\sqrt W$ increases, the process goes from tunneling to multiphoton emission \cite{EHe16}.





In situations in which the laser is turned on for only a short time, with pulses as short as a few laser periods, knowing the early time behavior is important.





In fact, we shall see later that the asymptotic approach to the periodic state considered in Faisal et al. \cite{FKS05}, which we discuss in detail in Section \ref{longtime}, can be much longer than a laser period.





Solving (\ref{schrodinger}) for $\omega\neq0$ turns out to be much more difficult than the constant field case, since here, the solution cannot be written in terms of known functions.





The very existence of physical solutions, which are bounded at infinity, is not mathematically obvious.





\medskip










\indent There are several studies in the literature about the solution of the timedependent Schr\"odinger equation, see \cite{KLe18} and references therein.





The method most often used is the CrankNicolson \cite{CN47} numerical scheme.





In the case treated here, however, as we will discuss in more detail in the following, because of the discontinuity in the potential at $x=0$, the CrankNicolson algorithm is actually inaccurate for short times.





We will introduce an alternative scheme that does not suffer from the discontinuity of the potential.





\indent In addition to the work in \cite{FKS05} there have been very many studies of the emission process using various approximations \cite{Wo35,Ke45,CFe87,Tr93,HKK06,KSH11,YGR11,KSe12b,PA12,YHe13,CPe14,ZL16,KLe18}.





Of particular note is the work of Yalunin, Gulde and Ropers \cite{YGR11} who consider the same setup as we (except for a phase difference in the field).





They use various analytic approximations for obtaining solutions of the periodic type considered in \cite{FKS05}.





They also carry out numerical solutions using the CrankNicolson method.





This method is discussed in section \ref{Sec3} and compared to the exact result in Figure \ref{fig:cn}.





As shown there the CrankNicolson method is not correct for very short times.





In particular, the current does not tend to its initial value, zero, as time tends to 0, see \cite[Figure 5]{YGR11}.





\medskip










\indent Of particular note is the work by Yalunin, Gulden and Roper \cite{YGR11}, where the initial condition and potential are the same as in (\ref{schrodinger}) and (\ref{initial}).





Our result is quite different from theirs for short time, but the values we obtain at longer times seem to agree qualitatively with those obtained in \cite{YGR11} for short times.





Furthermore, in \cite{YGR11}, there does not seem to be a distinction between long and short times, whereas we find that for some important physical quantities, like the current averaged over a laser period, one has to wait many periods to be close to the asymptotic value.





\bigskip





\indent The outline of the rest of the paper is as follows.





In section \ref{sec:math}, we give a brief description of the method used to solve (\ref{schrodinger}).





In section \ref{Sec3}, we present the results for the initial state $\psi(x,0)=\varphi_0(x)$ in (\ref{initial}).





In section \ref{longtime}, we describe the asymptotic form of $\psi(x,t)$ as $t\to\infty$.





The appendix contains more information about the derivation of the solution of (\ref{schrodinger}).










\indent We look for an exact solution of (\ref{schrodinger}) for $t>0$ such that $\psi(x,t)$ and its derivative $\partial_x\psi(x,t)$ are continuous at $x=0$ and be bounded for $x\to\infty$.





To do this, we first obtain the solutions for $x<0$ and for $x>0$ with general initial conditions and unknown time dependent boundary conditions $\psi(0,t)=\psi_0(t)$ and $\partial_x\psi(x,t)_{x=0}=\psi_{x,0}(t)$.





A solution of (\ref{schrodinger}) with appropriate $x\to\pm\infty$ conditions will exist if and only if we can match these boundary conditions coming from the right and from the left.





This approach to solving such an equation as (\ref{schrodinger}) goes under the general name of initialboundaryvalue (IBV) problem, developed in particular by Fokas for these types of problems \cite{Fo97,Ry18}.










\indent We prove the existence of a solution of (\ref{schrodinger}) by showing that a necessary and sufficient condition for this is that $\psi_0(t)$ satisfies an integral equation.





The existence of the solution $\psi(x,t)$ for all $x$ and $t$ then follows from the proof that this integral equation has a unique solution, which we provide.





Expressions for $\psi$ and for the current $j(x,t)$ are given by explicit integrals involving $\psi_0(t)$ \cite{CCe}.





These can be evaluated numerically.





To do so accurately for both short and moderate times, we developed a numerical scheme which goes beyond the CrankNicolson algorithm \cite{CN47}, see section\,\ref{Sec3}.





\bigskip










\indent $\psi(x,t)$ has the long time form $e^{i\frac12k^2t}\mathcal U(x,t)$, where $\mathcal U$ is (up to a phase) the periodic in time solution of (\ref{schrodinger}) derived by Faisal et al \cite{FKS05,ZL16}.





We note that Faisal et al. did not prove the existence of a solution to their infinite set of equations.





Instead, they solved numerically a truncated set of these equations.





Our results prove the existence of such a solution.










\indent The approach to this asymptotic form $\mathcal U$ goes like $t^{\frac32}$.





$\mathcal{U}(x,t)$ consists of resonance terms which can be associated, as was done in \cite{FKS05}, with the absorption (emission) of different number of photons, see Section\,\ref{longtime}.





These resonances show up dramatically when considering the dependence of the average current on $\omega$ in the vicinity of $\omega_c=W$, as seen in early experiments on the photoelectric effect where light with frequency $\omega>\omega_c$ was necessary to get electron emission.















\section{Existence of solutions of the Schr\"odinger equation}





\indent Here we give an outline of the proof derived in \cite{CCe}. To simplify the notations we rescale: let





\begin{equation}\label{rescale}





\bar t=\omega t,\ \ \bar x=\sqrt{2\omega}x,\ \ \bar U=U\omega^{1},\ \ \bar E=E\omega^{3/2}2^{1/2}





\end{equation}





Equation (\ref{schrodinger}) becomes





\begin{equation}\label{eq12}





i\partial_{\bar t}\psi(\bar x,\bar t)=\partial^2_{\bar x}\psi(\bar x,\bar t)+\Theta(\bar x)\left( \bar U2\bar E \bar x\cos(\bar t)\right)\psi(\bar x,\bar t)





\end{equation}





As mentioned earlier, we will solve (\ref{eq12}) with the general initial condition





$\psi(\bar x,0):=f(\bar x)$ separately for $\bar x<0$, $\psi_(\bar x,\bar t)$ and for $\bar x>0$, $\psi_+(\bar x,\bar t)$, and then match the values of $\psi_\pm(0,\bar t)$ and $\partial_{\bar x}\psi_\pm(0,\bar t)$.





\bigskip










\indent To obtain these solutions we first go from $\bar x$ to Fourier space.





We write





\begin{equation}





\hat\psi(\xi,\bar t)=\frac1{\sqrt{2\pi}}\int_{\infty}^\infty d\bar x\ e^{i\bar x\xi}\psi(\bar x,\bar t)





=\hat\psi_(\xi,\bar t)+\hat\psi_+(\xi,\bar t)





\end{equation}





where $\hat\psi_\pm$ is the halfline Fourier transform along $\pm \bar x>0$.





\bigskip










\indent The equation for $\hat\psi_(\xi,\bar t)$ has the form





\begin{equation}\label{eq5}





i\frac{\partial \hat{\psi}_}{\partial\bar t}=\xi^2 \hat{\psi}_ \frac1{\sqrt{2\pi}} \psi_{\bar x,0}(\bar t)i\xi \frac1{\sqrt{2\pi}} \psi_0(\bar t)





\end{equation}





where





\begin{equation}





\psi_0(\bar t):=\psi(0,\bar t)





,\quad





\psi_{\bar x,0}(\bar t):=\partial_{\bar x}\psi(0,\bar t)





\end{equation}





are treated as unknown functions to be determined later by the matching conditions.





The solution of (\ref{eq5}) is given by





\begin{equation} \label{psim}





\hat{\psi}_(\xi,\bar t)= {\rm e}^{i\xi^2\bar t}\left\{ C_(\xi) + \frac1{\sqrt{2\pi}} \int_0^{\bar t} {\rm e}^{i\xi^2\bar s}\ \left[ i\psi_{\bar x,0}(\bar s)\xi \psi_0(\bar s) \right]\, d\bar s\right\}





\end{equation}





where





$C_(\xi)=\int_{\infty}^0{\rm e}^{i\bar x\xi}f(\bar x)\, d\bar x$.





Inverting the Fourier transform we obtain, for all $\bar x<0$,





\begin{equation}\label{psiminus}





{\psi}_(\bar x,\bar t)=h_(\bar x,\bar t)+ \frac i{{2\pi}} \int_0^{\bar t} d\bar s\, \psi_{\bar x,0}(\bar s) H_0(\bar t\bar s,\bar x)+ \frac i{{2\pi}} \int_0^{\bar t} d\bar s\, \psi_0(\bar s) \, \partial_{\bar x}H_0(\bar t\bar s,\bar x)





\end{equation}





where





$$





h_(\bar x,\bar t) =\frac 1{{2\pi}}\int_{\infty}^0 f(\bar y) H_0(\bar t,\bar x\bar y)\, d\bar y,





\quad





H_0(\bar t,\bar r)= \frac {\sqrt{\pi}}{\sqrt{i}{\sqrt{\bar t}}} {\rm e}^{\frac{i\bar r^2}{4\bar t}}





$$





\bigskip










\indent The equation for $\hat{\psi}_+(\xi,\bar t)$ is given by





$$i\frac{\partial \hat{\psi}_+}{\partial\bar t}=i 2\bar E \cos(\bar t)\, \frac{\partial \hat{\psi}_+}{\partial \xi} +\left( \xi^2+\bar U\right) \hat{\psi}_+(\xi,\bar t) + \frac1{\sqrt{2\pi}} \psi_{\bar x,0}(\bar t)+i\xi \frac1{\sqrt{2\pi}} \psi_0(\bar t) $$





Solving for $\hat\psi_+(\xi,\bar t)$ and inverting the Fourier transforms, we obtain





\begin{equation}\begin{array}{>\displaystyle c}\label{calcpsip}





\psi_+(\bar x,\bar t)=h_+(\bar x,\bar t)+\\





\frac1{{2\pi}} {\rm e}^{2i\bar E \bar x\sin(\bar t)}\int_{\infty}^{\infty} du\,{\rm e}^{i\bar xu\Phi(u,\bar t)} \int_0^{\bar t} \, {\rm e}^{\Phi(u,\bar s)}\left[ i \psi_{\bar x,0}(\bar s)+2\bar E\sin\bar s\, \psi_0(\bar s)\right]\, d\bar s\\





+ \frac1{{2\pi}} {\rm e}^{2i\bar E \bar x\sin(\bar t)}\int_{\infty}^{\infty} du\,{\rm e}^{i\bar xu\Phi(u,\bar t)} \int_0^{\bar t} \, {\rm e}^{\Phi(u,\bar s)}\,u\, \psi_0(\bar s)d\bar s





\end{array}\end{equation}





where





\begin{equation}\begin{array}{>\displaystyle c}





\Phi(u,\bar t)=i\left[ (u^2 +\bar U+2\bar E^2)\bar t  4\bar E u \cos(\bar t) \bar E^2 \sin(2\bar t) \right], \\





h_+(\bar x,\bar t)=\frac 1{2\sqrt{i \pi }}{\rm e }^{2i\bar x\bar E\sin(\bar t)i(\bar U+2\bar E^2)\bar t+i\bar E^2\sin2\bar t}\int_0^{\infty }d\bar y\, f(\bar y)\, \frac 1{\sqrt{\bar t}}\,{\rm e}^{\frac{i(\bar x\bar y4\bar E+4\bar E\cos(\bar t))^2}{4\bar t}}





\end{array}\end{equation}





\bigskip










\indent Taking now the limit $\bar x\to0$ from left and right in (\ref{psiminus}) and in (\ref{calcpsip}) and setting $\psi_(0,\bar t)=\psi_+(0,\bar t)$, and similarly for $\partial_{\bar x}\psi_\pm $ we obtain (in fact $\psi_(0,\bar t)=\psi_0(\bar t)=\psi_+(0,\bar t)$ implies that $\partial_{\bar x}\psi_(0,\bar t)=\partial_{\bar x}\psi_+(0,\bar t)$) one equation with the only unknown function $\psi_0(\bar t)$, which has to satisfy the integral equation





\begin{equation}\begin{array}{>\displaystyle c}\label{eqcontr}





\psi_0(\bar t)= h_+(0,\bar t)+ h_(0,\bar t) \frac{1}{\pi}\int_0^{\bar t}\left( \int_0^{\bar s} \frac{h_(0,\bar s')}{\sqrt{\bar s\bar s'}}\, d\bar s'\right) \, G_1(\bar s,\bar t)\,d\bar s\\





+\frac{1}{2\pi}\int_0^{\bar t} \left( \int_0^{\bar s} \frac{\psi_0(\bar s')}{\sqrt{\bar s\bar s'}}\, d\bar s'\right) \, G_1(\bar s,\bar t)\, d\bar s \\





+\frac{\bar E}{\sqrt{i\pi}}\int_0^{\bar t} d\bar s \, \,\psi_0(\bar s)\, \frac 1{\sqrt{\bar t\bar s}}\, \left( \sin\bar s\, + \frac{ \cos(\bar t) \cos(\bar s)}{\bar t\bar s}\right) \,{\rm e}^{iF_0(\bar s,\bar t)}





\end{array}\end{equation}





where





$$G_1(\bar s,\bar t)=\frac d{d\bar s}\frac{e^{iF_0(\bar s,\bar t)}1}{\sqrt{\bar t\bar s}}$$





with





\begin{equation}\label{Fzero}





F_0(s,\bar t)=\frac{4\bar E^2 (\cos(\bar t) \cos(\bar s))^2}{\bar t\bar s} (\bar U_2+2\bar E^2)(\bar t\bar s) +\bar E^2(\sin 2\bar t\sin 2\bar s)





\end{equation}





We then prove that the integral equation (\ref{eqcontr}) has a unique solution, which implies that (\ref{eq12}) does as well.





\bigskip










\indent Setting the initial condition $f(\bar x)=\varphi_1(\bar x)$ (obtained from (\ref{initial}) after the rescaling (\ref{rescale}))





we obtain explicit functions in equation (\ref{eqcontr}) for $\psi_0(\bar t)$,





\begin{equation}\label{hminspec}





h_(0,\bar t)= {{\rm e}^{i{\bar k}^{2}\bar t}}\,\frac12\, \left[ 1+{\rm erf} \left( i^{3/2}\bar k\sqrt {\bar t}\right) \right]





R_0\, {{\rm e}^{i{k}^{2}\bar t}}\,\frac12\, \left[ 1+{\rm erf} \left( i^{3/2} \bar k\sqrt {\bar t}\right) \right]





\section{Solution of the Schr\"odinger equation}\label{sec:math}





\indent We solve (\ref{schrodinger}) by using the one sided Fourier transforms $\hat{\psi}_(\xi,t)=\frac{1}{\sqrt{2\pi}}\int_{\infty}^0{\rm e}^{i\xi x}\psi(x,t)\, dx$ and $\hat{\psi}_+(\xi,t)=\frac{1}{\sqrt{2\pi}}\int_0^{\infty}{\rm e}^{i\xi x}\psi(x,t)\, dx$. These satisfy the equations





\begin{equation}\label{A}





i\frac{\partial \hat{\psi}_(\xi,t)}{\partial t}\frac{\xi^2}2 \hat{\psi}_





= \frac1{\sqrt{2\pi}} \frac{\partial\psi}{\partial y}(0,t)i\xi \frac1{\sqrt{2\pi}} \psi(0,t)





\end{equation}





and





\begin{equation}\begin{array}{>\displaystyle c}\label{hpluspec}





h_+(0,\bar t)= T_0\,\exp[{i(\bar U+2\bar E^2)\bar t+i\bar E^2\sin 2\bar t4\,\bar E \alpha\,\cos\bar t +i{\alpha}^{2}\bar t+4\,\bar E \alpha}] \times \\





\frac 1{2} \left[1 {\rm erf} \left(2\,i^{3/2}\bar E\frac{\cos\bar t1}{\sqrt {\bar t}}+ \alpha \sqrt{i\bar t} \right) \right]





\end{array}\end{equation}





where





$$\bar k=k/\sqrt{2\omega},\ \ \ \alpha=\sqrt{\bar U\bar k^2},\ \ \ 1+R_0=T_0=\frac{2i\bar k}{i\bar k\sqrt{\bar U\bar k^2}}.$$





\begin{multline}\label{B}





i\frac{\partial \hat{\psi}_+(\xi,t)}{\partial t}+ i \frac E2 \cos (\omega t)\, \frac{\partial \hat{\psi}_+}{\partial \xi}





\left(\frac{\xi^2}2+U\right) \hat{\psi}_+(\xi,t) \\=





\frac1{\sqrt{2\pi}}\frac{\partial\psi}{\partial x}(0,t) +i\xi \frac1{\sqrt{2\pi}} \psi(0,t)





\end{multline}





Both (\ref{A}) and (\ref{B}) admit explicit solutions for initial values $\hat\psi_\pm(\xi,0)$ and specified boundary values $\psi(0,t)$ and $\partial_x\psi(0,t)$, see the Appendix. The continuity conditions for $\psi$ and $\partial_x\psi$ at $x=0$ then lead to an integral equation for $\psi(0,t)$ of the form





\begin{equation}





\label{eq:inteq}





\psi(0,t)=h(t)+L\psi(0,t)





\end{equation}





where $L$ is some compact integral operator whose expression is rather involved, see the appendix, and $h(t)$ is a function of the initial condition $\psi(x,0)$. We prove the existence and uniqueness of a physical solution of (\ref{eq:inteq}) for all $t>0$, by showing that $L$ is a contraction.





Given that solution $\psi(0,t)$, we can obtain $\psi(x,t)$ for all $x$ by direct integration.





To evaluate the solution numerically for the initial condition (\ref{initial}), we expand $\psi(0,t)$ in a Chebyshev polynomial series and identify the coefficients of this expansion. The complexity of this integral operator $L$ results in complex behavior of its solutions as discussed in the sequel.





The full mathematical proof of the existence and uniqueness of the solution of (\ref{schrodinger}) will be presented separately \cite{CCe}.










\section{Short time behavior: approach to the stationary state}\label{Sec3}





\indent We carried out numerical solution of (\ref{eqcontr}) using (\ref{hminspec}) and (\ref{hpluspec}) for a variety of values of $E$.





To compare with the work in \cite{YGR11}, we take $\frac {\hbar^2k^2}{2m}=\mathcal E_F=4.5$\,eV, $W=5.5$\,eV.





%\section{Long time behavior: approach to the stationary state}\label{sec:longtime}





%$\psi(x,t)$ has the long time form $e^{i\frac12k^2t}\mathcal U(x,t)$, where $\mathcal U$ is (up to a phase) the periodic in time solution of (\ref{schrodinger}) derived by Faisal et al \cite{FKS05,ZL16}.





%Faisal et al. conjectured the existence of a solution to their infinite set of equations, and computed a solution for a truncated subset of these equations numerically.





%Our results prove the existence of such a solution.





%





%\indent The wavefunction approaches its asymptotic form $\mathcal U$ like $t^{\frac32}$.





%$\mathcal{U}(x,t)$ consists of resonance terms which can be associated, as was done in \cite{FKS05}, with the absorption (emission) of different number of photons.





%These resonances show up dramatically when considering the dependence of the average current on $\omega$ in the vicinity of $\omega_c=W$, as seen in Figure \ref{fig:omega} and in experiments on the photoelectric effect where light with frequency $\omega>\omega_c$ is necessary to get electron emission \cite{Mi16}.










\section{Short time behavior}\label{Sec3}





\indent We carried out the numerical solution of the integral equation for $\psi_0(t)$ (\ref{eq:inteq}) with initial condition $\varphi_0(x)$ given in (\ref{initial}), using an exponentially convergent algorithm.





The numerical computation is based on expressing the solution $\psi_0(t)$ of (\ref{eq:inteq}) in terms of Chebyshev polynomials.





Since this solution becomes periodic at long times, it is actually convenient to first split time into small intervals, and expand $\psi_0(t)$ into Chebyshev polynomials in each interval.





Then, to compute the right side of (\ref{eq:inteq}), we must compute integrals of $\psi_0(t)$, which we carry out using Gauss quadratures.





Once that is done, (\ref{eq:inteq}) is approximated (by truncating the Chebyshev expansion and the Gauss quadratures) by a finite linear system of equations, which can be solved easily.





One has to pay attention to make sure that this approximation is good.





In this work, we have striven to ensure that the approximation converges to the exact solution exponentially fast as the truncations of the Chebyshev expansion and Gaussquadratures are removed.





This is not entirely trivial, as both $\psi_0$ and $L$ have square root singularities, so the Chebyshev polynomial expansion and the Gauss quadratures have to be adjusted to take these into account.





\medskip










\indent We take $\frac {\hbar^2k^2}{2m}=\mathcal E_F=4.5$\,eV, $W=5.5$\,eV.





Unless otherwise specified, we also take $\omega=1.55\ \mathrm{eV}$.





The laser period $\tau=\frac{2\pi}\omega$ is then equal to $2.7\ \mathrm{fs}$.





\bigskip





\medskip










\indent The details of the numerical algorithms, which take into account the singular behavior of $\psi(x,t)$ at $x=0$ and $t=0$, will be given in a separate publication.





\bigskip










\indent Figure \ref{fig:wavefunction} shows the real and imaginary parts of $e^{i\frac{k^2}2t}\psi_0(t)$, as well as $\psi_0(t)^2$ and the normalized current $\frac1kj(0,t)$.





While $\psi_0^2$ seems more in phase with the real part of$e^{i\frac{k^2}2t}\psi_0$, the current looks more in phase with the imaginary part.





\indent Figure \ref{fig:wavefunction} shows the density at the interface $\psi_0(t)^2$.





The maxima and minima of the density are approximately in phase with the field.





Figure \ref{fig:current} shows the values of the current $j(0,t)$ passing through the origin as a function of time for different strengths of the field, all at $\omega=1.55$\, eV.





We see there a change of behavior as $E$ increases from $E=1$\,V/nm to $E=30$\,V/nm (the Keldysh parameter $\gamma=\frac{2\omega}E\sqrt W$ goes from 18.6 to 0.62).





For large values of $E$, fast oscillations appear, which become faster and larger as $E$ grows.





The oscillatory behavior within one period for strong fields is also observed in the approximate solution \cite{YGR11}, though the details vary.





We have no simple physical explanation for these fast oscillations.





They do not fit in with the ``simple three step'' scenario \cite{KSK92,Co93,KLe18}.





In Figure \ref{fig:current+}, the current is plotted for various values of $\omega$.





It is seen there that the fast oscillations appear only for small values of $\omega$.





It is also apparent that the frequency of these oscillations is not a function of just the Keldysh parameter.





In Figure \ref{fig:currentx}, we show a plot of the current for positive values of $x$, and see that the fast oscillatory behavior within one period is strongly damped as $x$ increases.





The fact that the electrons cross the surface at different phases of the field does not fit in with the ``simple man'' scenario \cite{KSK92,Co93,KLe18}, where electrons are ejected out of the metal only when the field is positive, or even only when the field is at its maximum.





These oscillations at $x=0$ are also observed in the approximate solution \cite{YGR11}, though the details vary.





We note that the height of the first maximum is linear in $E$, while its location is almost independent of $E$.





\bigskip





Note also that the rapid oscillations occur mostly when the field is increasing.





In the Conclusions we further discuss these oscillations and their possible link to the physics of this process.





\medskip










\begin{figure}





\includegraphics[width=8cm]{real.pdf}





\hfill





\includegraphics[width=8cm]{imag.pdf}





\par\penalty10000\bigskip\penalty10000





\includegraphics[width=8cm]{density.pdf}





\hfill





\includegraphics[width=8cm]{current.pdf}





\caption{





The real (a) and imaginary (b) parts of $e^{i\frac{k^2}2t}\psi_0(t)$ the density (c) and the current (d) as a function of $\frac{t\omega}{2\pi}$ for 3 periods, with $E=15\ \mathrm V\cdot\mathrm{nm}^{1}$, $\omega=1.55\ \mathrm{eV}$.





The dotted line represents the phase of the field: $\cos(\omega t)$.





The density $\psi_0^2$ (recall that $\psi_0(t)\equiv\psi(0,t)$ is the wavefunction at $x=0$) as a function of $\frac{t\omega}{2\pi}$ for 3 periods, with $E=15\ \mathrm V\cdot\mathrm{nm}^{1}$, $\omega=1.55\ \mathrm{eV}$.





The dotted line is the graph of $\cos(\omega t)$ (not to scale).





}





\label{fig:wavefunction}





\end{figure}




@ 285,116 +192,140 @@ We note that the height of the first maximum is linear in $E$, while its locatio





\begin{figure}





\includegraphics[width=8cm]{current01.pdf}





\hfill





\includegraphics[width=8cm]{current05.pdf}





\par\penalty10000\bigskip\penalty10000





\includegraphics[width=8cm]{current18.pdf}





\includegraphics[width=8cm]{current15.pdf}





\hfill





\includegraphics[width=8cm]{current30.pdf}





\caption{





The normalized current $\frac jk$ as a function of $\frac{t\omega}{2\pi}$ for $\omega=1.55\ \mathrm{eV}$ and various values of the electric field: $E=1,\ 5,\ 18.6,\ 30\ \mathrm{V}\cdot\mathrm{nm}^{1}$.





The Keldysh parameter $\gamma=\frac{2\omega}E\sqrt W$ for these fields is, respectively, $18.6$, $3.73$, $1.00$ and $0.621$.





The dotted line represents the phase of the electric field: $\cos(\omega t)$.





The normalized current $\frac jk$ at the interface (recall that we are using atomic units, so $\frac jk$ is dimensionless) as a function of $\frac{t\omega}{2\pi}$ for $\omega=1.55\ \mathrm{eV}$ and three values of the electric field: $E=1,\ 15,\ 30\ \mathrm{V}\cdot\mathrm{nm}^{1}$.





The Keldysh parameter $\gamma=\frac{2\omega}E\sqrt W$ for these fields is, respectively, $18.6$, $1.24$ and $0.621$.





The dotted line is the graph of $\cos(\omega t)$ (not to scale).





As the field increases, fast oscillations in the current appear.





These fast oscillations mostly occur while the field is increasing.





}





\label{fig:current}





\end{figure}










\begin{figure}





\includegraphics[width=8cm]{currentk621.pdf}





\hfill





\includegraphics[width=8cm]{currentk124.pdf}





\caption{





The normalized current $\frac jk$ at the interface as a function of $\frac{t\omega}{2\pi}$ for two values of the Keldysh parameter.





In (a), $\gamma=0.621$, the blue curve has $E=30\ \mathrm{V}\cdot\mathrm{nm}^{1}$ and $\omega=1.55\ \mathrm{eV}$, and the red curve has $E=15\ \mathrm{V}\cdot\mathrm{nm}^{1}$ and $\omega=0.755\ \mathrm{eV}$.





In (b), $\gamma=1.24$, the blue curve has $E=15\ \mathrm{V}\cdot\mathrm{nm}^{1}$ and $\omega=1.55\ \mathrm{eV}$, and the red curve has $E=7.5\ \mathrm{V}\cdot\mathrm{nm}^{1}$ and $\omega=0.755\ \mathrm{eV}$.





At fixed $E$, the frequency of the oscillations decreases with $\omega$ (compare the red curve in (a) with the blue curve in (b)).





However, this frequency does not only depend on the Keldysh parameter.





}





\label{fig:current+}





\end{figure}










\begin{figure}





\includegraphics[width=8cm]{currentx.pdf}





\caption{





The normalized current $\frac jk$ as a function of $\frac{t\omega}{2\pi}$ at positive $x$.





The parameters here are $E=30\ \mathrm{V}\cdot\mathrm{nm}^{1}$ and $\omega=1.55\ \mathrm{eV}$, and the values of $x$ are $0.12\ \mathrm{nm}$ (red), $0.24\ \mathrm{nm}$ (green) and $0.37\ \mathrm{nm}$ (purple).





The fast oscillations die down as $x$ gets larger.





}





\label{fig:currentx}





\end{figure}










\indent In Figure \ref{fig:cn}, we show a comparison of our solution with the results obtained from a direct solution of (\ref{schrodinger}) via the CrankNicolson algorithm.





The agreement, especially for the location of the maxima and minima after very early times is very good.





At short times, however, the agreement is not so good.





It is not so, however, for very short times.





This is expected since the CrankNicolson scheme is based on approximating derivatives by finite differences.





However, at short times, $\psi(0,t)\sim t^{\frac32}$, which has a singular second derivative at $0$, so $\partial_t\psi$ is poorly approximated at short times by finite differences.





However, at short times, $\psi(0,t)\sim t^{\frac32}$, which has a singular second derivative at $0$, so $\partial_t\psi$ is poorly approximated by finite differences.





The fact that the result of the CrankNicolson algorithm is different at short times affects the values at later times, as the initial error effectively changes the initial condition.





The agreement at the end of one period is still rather good, indicating that the behavior of the solution behaves weakly on the initial condition.










Note, also, that the CrankNicolson method requires a cutoff in $x$, restricting $\psi(x,t)$ to $x\in[a,a]$.





This causes distortions in $\psi$ due to reflections from these artificial boundaries, and so can only be used for short times (this can be avoided by using nonlocal boundary conditions, such as ``transparent boundary conditions'').





\bigskip





This causes distortions in $\psi$ due to reflections from these artificial boundaries, and so can only be used for short times (reflections can be avoided by using nonlocal boundary conditions, such as ``transparent boundary conditions'').





\medskip










\begin{figure}





\includegraphics[width=8cm]{cn.pdf}





\hfill





\includegraphics[width=8cm]{cn_short.pdf}





\caption{





The current computed with our method (blue, color online), compared with the CrankNicolson algorithm (red), for $\omega=1.55\ \mathrm{eV}$ and $E=15\ \mathrm{V}\cdot\mathrm{nm}^{1}$.





The maxima and minima seem to occur at the same time, and the agreement is somewhat good for $t>\frac\pi\omega$.





The graph on the right focuses on short times, for which the CrankNicolson algorithm produces a different, and quite surprising result: the current initially shoots down to negative values before rising back up, but the field is initially positive, which should drive the current up, not down.





The discrepancy is due to the fact that the wavefunction is singular at $t=0$ (it has a divergent second derivative), and that the CrankNicolson algorithm approximates derivatives using finite differences, which is good only for regular functions.





The maxima and minima seem to occur at the same time, and the agreement is pretty good for $t>\frac\pi\omega$.





The inset focuses on short times, for $\frac{t\omega}{2\pi}<0.0005$, for which the CrankNicolson algorithm produces a different, and unphysical result: the current initially shoots down to negative values before rising back up.





}





\label{fig:cn}





\end{figure}










\indent Also shown in the figures is the running average of the current










\indent In Figure \ref{fig:avgcurrent}, we plot the running average of the current





\begin{equation}





\left<j\right>_t:=\frac1\tau\int_{t\tau}^t j(0,t)





\left<j\right>_t:=\frac1\tau\int_{t\tau}^tds\ j(0,s)





,\quad





\tau:=\frac{2\pi}\omega





.





\end{equation}





It is seen there that while the current appears to approach its long time value rather rapidly it does not do so fully until much later.





As seen in Figure \ref{fig:avgcurrent}, where we plot $\left<j\right>$ for $\omega=6\ \mathrm{eV}$, $E=10\ \mathrm{V}\cdot\mathrm{nm}^{1}$, $\gamma=9.6$, $\left<j\right>_t$ is still not all that close to its asymptotic value even when $t\approx48\tau$.





\bigskip





We plot $\left<j\right>_t$ at $x=0$ for $\omega=6\ \mathrm{eV}$, $E=10\ \mathrm{V}\cdot\mathrm{nm}^{1}$, $\gamma=9.6$.





It is seen there that the relative deviation of $\left<j\right>_t$ from its constant asymptotic value, described in section \ref{longtime}, remains significant even when $t\approx48\tau$.





\medskip










\indent The rate of the decay of the average current to its asymptotic value is evaluated in Figure \ref{fig:decay}.





In order to compute this rate without having to guess the asymptotic value, we proceed as follows.





At the end of every laser period $t_n=\frac{2\pi}\omega n$, we compute the minimal value $\mu_n$ and maximal value $M_n$ of the average current in the period $(\frac{2\pi}\omega(n1),\frac{2\pi}\omega n]$.





The plot shoes $M_n\mu_n$ as a function of $t$, and shows that $\left<j\right>_t\approx\left<j\right>_\infty+g(t)t^{\frac32}$, where $g(t)$ is bounded and asymptotically constant.





\bigskip





The plot shows $M_n\mu_n$ as a function of $t$, and shows that $\left<j\right>_t\approx\left<j\right>_\infty+g(t)t^{\frac32}$, where $g(t)$ is bounded and asymptotically constant.





This is consistent with the exact result in section \ref{longtime}.





\medskip










\begin{figure}





\includegraphics[width=8cm]{avgcurrent.pdf}





\hfill





\includegraphics[width=8cm]{avgcurrent_zoom.pdf}





\hfil\includegraphics[width=8cm]{avgcurrent.pdf}





\caption{





The normalized average current $\frac1k\left<j\right>_t$ as a function if $\frac{t\omega}{2\pi}$ for $\omega=6\ \mathrm{eV}$ and $E=10\ \mathrm{V}\cdot\mathrm{nm}^{1}$.





The two plots show the same data, but the first one focuses on times between the 12th and 48th period.





The normalized average current $\frac1k\left<j\right>_t$ as a function of $\frac{t\omega}{2\pi}$ at $x=0$ (blue) and $x=0.37\ \mathrm{nm}$ (red) for $\omega=6\ \mathrm{eV}$ and $E=10\ \mathrm{V}\cdot\mathrm{nm}^{1}$.





The inset shows the same data, restricted to between the 12th and 48th period.





Here, $\omega$ is large enough that absorbing one photon suffices to overcome the work function.





Even after 48 periods, the average current has not converged to its asymptotic value.





Even after 48 periods the average current has not converged to its asymptotic value.





}





\label{fig:avgcurrent}





\end{figure}










\begin{figure}





\hfil\includegraphics[width=8cm]{decay.pdf}





\hfil\includegraphics[width=9cm]{decay.pdf}





\caption{





The convergence rate to the asymptote of the average current as a function of $\frac{t\omega}{2\pi}$ on a loglog plot for $\omega=6$ and $E=10\ \mathrm{V}\cdot\mathrm{nm}^{1}$.





The convergence rate to the asymptote of the average current as a function of $\frac{t\omega}{2\pi}$ on a loglog plot for $\omega=6\ \mathrm{eV}$ and $E=10\ \mathrm{V}\cdot\mathrm{nm}^{1}$.





The dots are computed at the end of each period $t_n=\frac{2\pi}\omega n$, and their value is the difference between the maximun $M_t$ and the minimum $\mu_t$ of the normalized average current $\frac1k\left<j\right>_t$ in the period immediately preceding $t_n$.





The red line is a plot of $0.0030\times(\frac{t\omega}{2\pi})^{\frac32}$, which fits the data rather well.





}





\label{fig:decay}





\end{figure}










\indent In Figure \ref{fig:omega}, we show an estimate of the asymptotic average current as a function of $\omega$ in the vicinity of the onephoton threshold $\omega_c=W$ for $E=10\ \mathrm{V}\cdot\mathrm{nm}^{1}$.





As we argued above, the average current converges slowly to its asymptotic value.





In order to estimate this asymptotic value without computing the current at very large times, we actually compute the average of the average current over a laser period:





\indent In Figure \ref{fig:omega}, we show an estimate of the asymptotic average current as a function of $\omega$ in the vicinity of the onephoton threshold $\omega_c=W+\frac{E^2}{4\omega^2}$ for $E=3,10,30\ \mathrm{V}\cdot\mathrm{nm}^{1}$.





In order to reduce the fluctuations and estimate the longtime average current, we took a second average, and computed the average over a laser period of the average current, defined as





\begin{equation}





\left<\left<j\right>\right>:=\frac1\tau\int_{T\tau}^T dt\ \left<j\right>_t





.





\end{equation}





By dividing the current by $\epsilon^2$, we see that the the average of the average of the current is proportional to $\epsilon^2$.





We see that there is a steep increase in $\left<\left<j\right>\right>$ as $\omega$ increases past $\omega_c$.





This is precisely what is observed in experiments on the photoelectric effect, where the emission of electrons from the metal surface has such a threshold \cite{Mi16}.





This became a key element in Einstein's ansatz of localized photons.





Here, in the semiclassical treatment of the laser field, this phenomenon appears as a resonance.





The ponderomotive term $\frac{E^2}{4\omega^2}$, which appears in the photon energy, see below, is negligible compared to the work function $W$.





\bigskip





Here, in the classical treatment of the laser field, this phenomenon appears as a consequence of the quantum treatment of the electrons.





It shows that, despite it simplicity, this model captures essential features of the physical phenomena.





\medskip










\begin{figure}





\hfil\includegraphics[width=8cm]{omega.pdf}





\caption{





The average of the average of the normalized current $\frac1k\left<\left<j\right>_t\right>$ after 12 periods as a function of $\omega$, for $E=10\ \mathrm{V}\cdot\mathrm{nm}^{1}$.





After 12 periods, the average current $\left<j\right>_t$ is still fluctuating quite a lot, so, in order to eliminate these fluctuations and get closer to the asymptotic value $\left<j\right>_\infty$, we average the average current over a laser period.





The average of the average of the current $\frac1{\epsilon^2}\left<\left<j\right>_t\right>$ after 12 periods as a function of $\omega\omega_c$, for various values of the field: $E=3\ \mathrm{V}\cdot\mathrm{nm}^{1}$ (blue), $E=10\ \mathrm{V}\cdot\mathrm{nm}^{1}$ (red), $E=30\ \mathrm{V}\cdot\mathrm{nm}^{1}$ (green).





We see a sharp transition as $\omega$ crosses $\omega_c$.





}





\label{fig:omega}





\end{figure}










\section{Long time behavior of $\psi(x,t)$}\label{longtime}





\indent The long time behavior of the system is given by the poles of the Laplace transform





$$\hat{\psi}(x,p)=\int_0^{\infty}{\rm e}^{pt} \psi(x,t)\, dt$$





\begin{equation}





\hat{\psi}(x,p)=\int_0^{\infty}dt\ e^{pt} \psi(x,t)





\end{equation}





on the imaginary axis, as can be seen by taking the inverse Laplace transform





\begin{equation}





\psi(x,t)=\int_{i\infty}^{i\infty}dp\ e^{pt}\hat\psi(x,p)





.





\end{equation}





Other poles, which necessarily have negative real parts, give rise to exponentially decaying terms while branch cuts generally contribute $t^{\frac32}$ terms to the approach to the asymptotic state.





As mentioned earlier, the time asymptotic of $\psi$ is of the form





where the path of integration can be deformed to get contributions only from poles and branch cuts in the negative real halfplane $\mathcal Re(p)\leqslant 0$.





The poles that have negative real parts would give rise to exponentially decaying terms while branch cuts generally contribute $t^{\frac32}$ terms to the approach to the asymptotic state.





The contribution from poles on the imaginary $p$axis then give the longtime asymptotics of $\psi$ is of the form





\begin{equation}\label{4p1}





\psi(x,t)\sim {\rm e}^{\frac 12 i k^2 t} \bar{\psi}(x,t)





\end{equation}




@ 403,20 +334,22 @@ This corresponds to the poles of $\hat{\psi}(p,t)$ on the imaginary axis being a





\begin{equation}\label{4p2}





p=\frac 12 i k^2+i\omega n,\ \ \ n\mathrm{\ \ integer}





\end{equation}





\bigskip





\medskip










\indent We further show that $\bar{\psi}(x,t)$ coincides, after the change of gauge, with the wave function $\mathcal{U}(x,t)$ computed by Faisal et al \cite{FKS05}.





In that work it was assumed (\ref{schrodinger}) (in the magnetic gauge) has a solution of the form (\ref{4p1}) with an incoming beam $e^{ikx}$ for $x<0$, without considering any initial conditions.





Computing the residues at the poles on the imaginary axis we show that





\begin{equation}\label{4p3}





\bar\psi(x,t)\sim \left\{\begin{array}{>\displaystyle ll}





e^{ikx}





+\sum_{m\in\mathbb Z}e^{im\omega t}e^{ix\sqrt{k^2+2m\omega}}\mathcal R_m





&\mathrm{for\ }x<0





\\[0.5cm]





e^{i\frac E\omega x\sin\omega t}\sum_{n,m\in\mathbb Z}e^{in\omega t}g_{nm}^{(\kappa_m)}e^{\kappa_mx}\mathcal T_m





&\mathrm{for\ }x>0





\end{array}\right.





\indent As shown in \cite{CCe}, $\bar{\psi}(x,t)$ coincides with the wave function $\mathcal{U}(x,t)$ computed by Faisal et al \cite{FKS05}.





In that work it was assumed that (\ref{schrodinger}) has a solution of the form (\ref{4p1}) with an incoming beam $e^{ikx}$ for $x<0$, without considering any initial conditions.





Computing the residues at the poles on the imaginary axis we show that, for $x<0$,





\begin{equation}





\bar\psi(x,t)\sim





e^{ikx}





+\sum_{m\in\mathbb Z}e^{im\omega t}e^{ix\sqrt{k^2+2m\omega}}\mathcal R_m





\label{asym1}





\end{equation}





and for $x>0$,





\begin{equation}





\bar\psi(x,t)\sim





e^{i\frac E\omega x\sin\omega t}\sum_{n,m\in\mathbb Z}e^{in\omega t}g_{nm}^{(\kappa_m)}e^{\kappa_mx}\mathcal T_m





\label{asym2}





\end{equation}





where





\begin{equation}




@ 428,41 +361,200 @@ and





g_{nm}^{(\kappa_m)}=\frac\omega{2\pi}\int_0^{\frac{2\pi}\omega}dt\ e^{i(nm)\omega t}e^{\frac{i\frac{E^2}{4\omega^2}}\omega\sin(2\omega t)+\kappa_m\frac{2E}{\omega^2}\cos(\omega t)}





.





\end{equation}





This is exactly of the form considered in \cite{FKS05}.





$\mathcal R_n$ and $\mathcal T_m$ are computed by matching boundary values of $\psi(0,t)$ and $\psi_{x,0}(t)$ at $x=0$.





\bigskip





This is exactly of the form obtained in \cite{FKS05}.





$\mathcal R_n$ and $\mathcal T_m$ are computed by matching boundary values of $\psi(x,t)$ and $\partial_x\psi(x,t)$ at $x=0$.





The phase $e^{i\frac E\omega x\sin\omega t}$ comes from a change of gauge with respect to \cite{FKS05} (we use the ``length'' gauge, instead of the ``magnetic gauge'' \cite{CFe87}).





\medskip










\indent A physical interpretation of (\ref{4p3}), see \cite{FKS05}, is that an electron in a beam coming from $\infty$ and moving in the positive $x$direction, ${\rm e}^{ikx},\ k>0$, absorbs or emits ``$m$ photons'' and is either reflected, transmitted or damped.





\indent A physical interpretation of (\ref{asym1})(\ref{asym2}), see \cite{FKS05}, is that an electron in a beam coming from $\infty$ and moving in the positive $x$direction, ${\rm e}^{ikx},\ k>0$, absorbs or emits ``$m$ photons'' and is either reflected, transmitted or damped.





Transmission occurs when $m\omega>U+\frac {E^2}{4\omega^2}\frac{k^2}2\equiv\omega_c$.





The $\frac{E^2}{4\omega^2}$ in (\ref{kappa}) corresponds to the ponderomotive energy of the electron in the laser field.





Damping occurs when the inequalities are in the opposite direction.





$\omega_c$ is the minimum frequency necessary to push the electron with kinetic energy $\frac12k^2$ (in the $x$direction) over the potential barrier of height $U\frac12k^2=W$.





For $x\gg0$, the current in the $m$photon channel will have kinetic energy $m\omega\omega_c$ and the current will be given by $\sqrt{m\omega\omega_c}$.





Since we are taking $\frac{k^2}2=\mathcal E_F$, we have that $U\frac{k^2}2=W$, the work function.





The $\frac{E^2}{4\omega^2}$ in (\ref{kappa}) term corresponds to the ponderomotive energy of the electron \cite{Wo35} in the laser field.





Damping occurs when $m\omega<\omega_c$.





$\omega_c$ is the minimum frequency necessary to let the electron with incoming kinetic energy $\frac12k^2$ (in the $x$direction) propagate to the right of the potential barrier.





For large $x>0$, the current in the $m$photon channel will have kinetic energy $m\omega\omega_c$ and the current will be given by $\sqrt{m\omega\omega_c}$.





This will also be equal to the average current at large $t$, which is independent of $x$.





This explains the picture in Figure 4 for $m=1$.





This explains the picture in Figure \ref{fig:omega} for $m=1$.





The larger $m$ values necessary for smaller $\omega$ are difficult to see.





\bigskip





\medskip










\indent The asymptotic form (\ref{4p3}) is true for all initial conditions of the form $\Theta(x)e^{ikx}+f_0(s)$ as long as $f_0(x)$ only contains terms which are square integrable.





\section{Concluding remarks}





\indent In this paper we presented, for the first time, the exact solution of the timedependent Schr\"odinger equation (\ref{schrodinger}).





This is the simplest physical model describing the emission of electrons from a flat metal surface by an oscillating field.





The model was first used by Fowler and Nordheim \cite{FN28} for emission by a constant electric field, $\omega=0$.





Their formula for the steady state current, obtained from the stationary solution of (\ref{schrodinger}), with $\omega=0$, still forms the basis of the interpretation of experiments at present time \cite{Je17}.





There are modifications due to the Shottky effect, but these are not expected to change the basic results.





The situation is different when the field acting on the metal surface is periodic in time.





Equation (\ref{schrodinger}) no longer has an explicit ``stationary'' (in this case, periodic) solution of the type considered by Faisal et al \cite{FKS05}.





In fact, even the existence of physical solutions of (\ref{schrodinger}), i.e. ones bounded for all $x$, is problematic from a mathematical point of view.





This is what we establish here by proving that the integral equation (\ref{eq:inteq}) indeed has solutions which give a physical $\psi(x,t)$ for all $t>0$.





We prove this for a very general class of initial conditions and carry out exact numerical solutions for the case of a particular physically motivated initial state.





The numerics are proven to give arbitrary accuracy for any fixed $t$ and specified bounded initial state.





\medskip










\indent Our results reveal, as shown in the figures, many new features of the exact solution of (\ref{schrodinger}), e.g. the slow convergence in time of the average of the current to its asymptotic value, and the rapid oscillations at the interface for strong fields and small $\omega$.















A more detailed examination of our solution shows that the rapid oscillations are confined to a very narrow region close to the metal surface. A time Fourier transform of the wave function  which corresponds to looking in energy space  indicates that these fast oscillations are due to energy absorptions, $E_n=n \hbar\omega$ for all $n$ such that $E_n$ exceeds the work function plus the ponderomotive energy. Farther away from the metal surface, due to the transition to a semiclassical behavior, energy absorption and hence the rapid oscillations, cease rapidly. The purely quantum processes occur in the tunneling region proximal to the surface.










The slow convergence of the average of the current indicates that different initial conditions may give different results in short pulse experiments.





On the other hand, we prove that there is indeed an asymptotic periodic state of the form assumed by Faisal et al. \cite{FKS05}.





The asymptotic form (\ref{asym1})(\ref{asym2}) is true for all initial conditions of the form $\Theta(x)e^{ikx}+f_0(x)$ as long as $f_0(x)$ only contains terms which are square integrable.





The additional terms in $\psi(x,t)$ which come from $f_0(x)$ go to zero as $t\to\infty$.





This follows from the fact, proven in \cite{CCe}, that the Floquet operator associated to (\ref{schrodinger}) has no point spectrum.





\bigskip





\medskip










\begin{acknowledgements}





We thank David Huse, Kevin Jensen and Donald Schiffler for very valuable discussions, as well as the anonymous referees who made invaluable comments and helped improve this paper. This work was supported by AFOSR Grant No. FA95001610037. O.C. was partially supported by the NSFDMS (Grant No. 1515755). I.J. was partially supported by the NSFDMS (Grants No. 31128155 and 1802170). J.L.L. thanks the Institute for Advanced Study for its hospitality.





\end{acknowledgements}










\bibliographystyle{apsrev42}





\bibliography{bibliography}










\vfill










\delimtitle{\bf Acknowledgements}





We thank David Huse, Kevin Jensen and Donald Schiffler for very valuable discussions. This work was supported by AFOSR Grant No. FA95001610037. O.C. was partially supported by the NSFDMS (Grant No. 1515755). I.J. was partially supported by the NSFDMS (Grants No. 31128155 and 1802170). J.L.L. thanks the Institute for Advanced Study for its hospitality.





\enddelim










\eject










\begin{thebibliography}{WWW99}





\small





\IfFileExists{bibliography/bibliography.tex}{\input bibliography/bibliography.tex}{}





\end{thebibliography}





\appendix





\section{The exact solution of \eqref{schrodinger}}










Let $\psi_0(t)=\psi(0,t)$ and $\partial_x\psi_{0}(t)=\partial_x\psi(x,t)_{x=0}$.





The operator $L$ in (\ref{eq:inteq}) is given by





\begin{equation}





\begin{array}{>\displaystyle l}





L\psi_0(t):=





\frac{E}{2\omega\sqrt{2i\pi}}\int_0^tds\ \psi_0(s)\frac{\alpha(s,t)}{\sqrt{ts}}e^{if(s,t)}





\\[0.5cm]\hfill





+\frac1{2\pi}\int_0^t du\ \psi_0(u)\int_u^t ds\ \frac{g(s,t)}{\sqrt{su}}





\end{array}





\end{equation}





where





\begin{equation}





\alpha(s,t):=\sin(\omega s)+\frac{\cos(\omega t)\cos(\omega s)}{\omega(ts)}





\label{alpha}





\end{equation}





\begin{equation}





\begin{array}{>\displaystyle l}





f(s,t):=





\frac{E^2(\cos(\omega t)\cos(\omega s))^2}{2\omega^4(ts)}





\\[0.5cm]\hfill





\left(V+\frac{E^2}{4\omega^2}\right)(ts)





+\frac{E^2}{8\omega^3}(\sin(2\omega t)\sin(2\omega s))





\end{array}





\end{equation}





and





\begin{equation}





g(s,t):=





\frac{e^{if(s,t)}1}{2(ts)^{\frac32}}+\frac{i\partial_sf(s,t)e^{if(s,t)}}{\sqrt{ts}}





.





\end{equation}










The function $h$ in (\ref{eq:inteq}) is given by





\begin{equation}





\begin{array}{r@{\ }>\displaystyle l}





h(t):=&h_+(0,t)+h_(0,t)





\\[0.5cm]&





\frac1{\pi}\int_0^tdu\ h_(u)\int_u^t ds\ \frac{g(s,t)}{\sqrt{su}}





\end{array}





\end{equation}





where





\begin{equation}





h_(x,t):=\frac{e^{\frac{i\pi}4}}{\sqrt{2\pi t}}\int_{\infty}^0 dy\ \varphi_0(y)e^{\frac i{2t}(xy)^2}





\end{equation}





and





\begin{equation}





\begin{array}{r@{\ }>\displaystyle l}





h_+(x,t):=&





\frac{e^{\frac{i\pi}4}}{\sqrt{2\pi t}}





e^{i\frac E\omega x\sin(\omega t)i(V+\frac{E^2}{4\omega^2})t+\frac{iE^2}{8\omega^3}\sin(2\omega t)}





\\[0.5cm]&\times





\int_0^\infty dy\ \varphi_0(y)





e^{\frac i{2t}(x+y+\frac{E}{\omega^2}(1\cos(\omega t)))^2}





\end{array}





\end{equation}





For our choice (\ref{initial}) of $\varphi_0$, these two functions are explicit:





\begin{equation}





\begin{array}{>\displaystyle l}





h_(x,t)=\frac{e^{\frac{ik^2}{2m}t}}2\Big(





e^{ikx}





\mathrm{erfc}({\textstyle e^{\frac{i\pi}4}(\sqrt{\frac t2}k+\frac1{\sqrt{2t}}x})





\\[0.5cm]\hfill





+\frac{ik+\sqrt{2Uk^2}}{ik\sqrt{2Uk^2}}





e^{ikx}





\mathrm{erfc}({\textstyle e^{\frac{i\pi}4}(\sqrt{\frac t2}k+\frac1{\sqrt{2t}}x)})





\Big)





.





\end{array}





\end{equation}





\medskip





and





\begin{equation}





\begin{array}{>\displaystyle l}





h_+(x,t)=





\frac{ik}{ik\sqrt{2Uk^2}}





e^{i\frac E\omega\sin(\omega t)x\sqrt{2Uk^2}x}





\\[0.5cm]\hfill\times





e^{\frac{E}{\omega^2}(1\cos(\omega t))\sqrt{2Uk^2}i(\frac{k^2}2+\frac{E^2}{4\omega^2})t+i\frac{E^2}{8\omega^3}\sin(2\omega t)}





\\[0.5cm]\hfill\times





\mathrm{erfc}(e^{\frac{i\pi}4}({\textstyle i\sqrt{\frac t2}\sqrt{2Uk^2}+\frac E{\omega^2}\frac{1\cos(\omega t)}{\sqrt{2t}}\frac1{\sqrt{2t}}x}))





.





\end{array}





\end{equation}





$\hat{\psi}_$ is obtained by explicitly solving (\ref{A}):





\begin{equation}\begin{array}{c}





\hat{\psi}_(\xi,t)= e^{i\frac{\xi^2}2t} \frac{1}{\sqrt{2\pi}}\int_{\infty}^0 e^{i\xi x}\varphi_0(x)\, dx





\\ \\





+\frac1{2\sqrt{2\pi}} \int_0^t e^{i\frac{\xi^2}2(ts)}\ \left[ i\partial_x\psi_0(s)\xi \psi_0(s) \right]\, ds





\end{array}\end{equation}





The PDE (\ref{B}) can be solved explicitly by characteristics to give $\hat{\psi}_+$:





\begin{equation}





\hat{\psi}_+(\xi,t) =G(\xi \frac E{\omega}\sin \omega t,t)





\end{equation}





where





\begin{multline}





G(u,t)=e^{i\Phi(u,t)} \frac{1}{\sqrt{2\pi}}\int_0^{\infty}e^{ikx}\varphi_0(x)\, dx +





\\





+\frac1{2\sqrt{2\pi}}\int_0^t ds\, e^{i(\Phi(u,t)\Phi(u,s))}\cdot\hfill\\





\cdot\left[ i \partial_x\psi_0(s)+(u+\frac E{\omega}\sin (\omega s)) \psi_0(s)\right]





\end{multline}





where





\begin{multline}





\Phi(u,t)=\\





=\left(\frac{u^2}2 +V+\frac{E^2}{4\omega^2}\right)t + u \frac{E}{\omega^2}(1\cos(\omega t))\frac{E^2}{8\omega^3} \sin(2\omega t)





.





\end{multline}





One can then check that





\begin{equation}





\label{eq:postconvo}





\partial_x\psi_0(t)=\frac{\sqrt2}{\sqrt{i \pi}} \frac d{dt}\left[ \psi_0(t)*t^{1/2} 2 h_(0,t)*t^{1/2}\right]





\end{equation}





where '*' denotes the Laplace convolution





\begin{equation}





[f*g](t)=\int_0^t ds\ f(s)g(ts)





\end{equation}





is continuous for $t>0$.





\bigskip










The solution $\psi(x,t)$ of the Schr\"odinger equation (\ref{schrodinger}) is, for $x<0$, the inverse Fourier transform of $\hat{\psi}_$, while for $x>0$ it equals the inverse Fourier transform of $\hat{\psi}_+$. Namely,





\begin{multline}





\label{eq:psimxt}





{\psi}_(x,t)=h_(x,t)+\\+ \frac{e^{\frac{i\pi}4}}{2\sqrt{2\pi}} \int_0^t ds\, \left(\partial_x\psi_0(s)+i\psi_0(s)\frac x{ts}\right)\frac{e^{\frac{ix^2}{2(ts)}}}{\sqrt{ts}}\\





\end{multline}





and





\begin{multline}





\psi_+(x,t)=h_+(x,t)\\





e^{i\frac E\omega\sin(\omega t)x}\frac{e^{\frac{3i\pi}4}}{2\sqrt{2\pi}}\int_0^t ds\frac{\Gamma_+(s,x,t)}{\sqrt{ts}}e^{iF(x,s,t)}





\end{multline}





where





\begin{multline}





\label{eq:h0}





\Gamma_+(s,x,t):=





i\partial_x\psi_0(s)+\\





+\left(\frac E\omega\sin(\omega s)+\frac E{\omega^2}\frac{\cos(\omega t)\cos(\omega s)}{ts}+\frac x{ts}\right)\psi_0(s)





\end{multline}





and





\begin{equation}





\label{eq:defx}





F(x,s,t)=f(s,t)+x\frac{E}{\omega^2}\frac{\cos(\omega t)\cos(\omega s)}{ts}





+\frac{ix^2}{2(ts)}





.





\end{equation}





\end{document}














