Initial commitv0.0

1 files changed, 468 insertions, 0 deletions

diff --git a/Costin_Costin_Jauslin_Lebowitz_2019.tex b/Costin_Costin_Jauslin_Lebowitz_2019.tex
new file mode 100644
index 0000000..5f4c367
--- /dev/null
+++ b/Costin_Costin_Jauslin_Lebowitz_2019.tex
@@ -0,0 +1,468 @@
+\documentclass{ian}
+
+\usepackage{graphicx}
+\usepackage{array}
+
+\begin{document}
+
+\hbox{}
+\hfil{\bf\LARGE
+Exact solution of the Schr\"odinger equation\par
+\hfil for photoemission from a metal\par
+}
+\vfill
+
+\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
+
+\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
+
+\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
+
+\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 one-dimensional, with the half-line $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}
+
+
+\section{Introduction}
+\indent In this note, we continue our investigation of the exact solution of the time-dependent Schr\"odinger equation describing the emission of electrons from a metal surface by an electric field \cite{CCe18}.
+We use the Sommerfeld model of quasi-free 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 The time evolution of the wave function of an electron in such a beam is described by the one dimensional Schr\"odinger equation,
+\begin{equation}
+  i\partial_t\psi(x,t)=H(x,t)\psi(x,t)=-\frac12\Delta\psi(x,t)+\Theta(x)(U-Ex\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}.
+
+\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,
+\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
+  \end{array}\right.
+\end{equation}
+Using the continuity of $\psi$ and its derivative at $x=0$ yields
+\begin{equation}
+  \alpha:=\sqrt{2U-k^2}
+  ,\quad
+  1+R_0=T_0=\frac{2ik}{ik-\sqrt{2U-k^2}}
+  .
+\end{equation}
+The current,
+\begin{equation}
+  j(x,t):=\mathcal Im(\psi^*(x,t)\partial_x\psi(x,t))
+\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.
+
+\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
+\begin{equation}
+  \frac12k^2\varphi_E=-\frac12\frac{\partial^2}{\partial x^2}\varphi_E+\Theta(x)(U-Ex)\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
+
+\indent In \cite{CCe18} we solved the time-dependent 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 Fowler-Nordheim (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.
+
+\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 There are several studies in the literature about the solution of the time-dependent Schr\"odinger equation, see \cite{KLe18} and references therein.
+The method most often used is the Crank-Nicolson \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 Crank-Nicolson algorithm is actually inaccurate for short times.
+We will introduce an alternative scheme that does not suffer from the discontinuity of the potential.
+
+\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 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 initial-boundary-value (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 Crank-Nicolson 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 U-2\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 half-line 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 y-4\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]
+\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 t-4\,\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 t-1}{\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}}.$$
+
+\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.
+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
+
+\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.
+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}.
+We note that the height of the first maximum is linear in $E$, while its location is almost independent of $E$.
+\bigskip
+
+\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)$.
+  }
+  \label{fig:wavefunction}
+\end{figure}
+
+\begin{figure}
+  \includegraphics[width=8cm]{current01.pdf}
+  \hfill
+  \includegraphics[width=8cm]{current05.pdf}
+  \par\penalty10000\bigskip\penalty10000
+  \includegraphics[width=8cm]{current18.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)$.
+    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}
+
+\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 Crank-Nicolson 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.
+This is expected since the Crank-Nicolson 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.
+Note, also, that the Crank-Nicolson method requires a cut-off 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 non-local boundary conditions, such as ``transparent boundary conditions'').
+\bigskip
+
+\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 Crank-Nicolson 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 Crank-Nicolson 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 Crank-Nicolson algorithm approximates derivatives using finite differences, which is good only for regular functions.
+  }
+  \label{fig:cn}
+\end{figure}
+
+\indent Also shown in the figures is the running average of the current
+\begin{equation}
+  \left<j\right>_t:=\frac1\tau\int_{t-\tau}^t j(0,t)
+  ,\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
+
+\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(n-1),\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
+
+\begin{figure}
+  \includegraphics[width=8cm]{avgcurrent.pdf}
+  \hfill
+  \includegraphics[width=8cm]{avgcurrent_zoom.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.
+    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.
+  }
+  \label{fig:avgcurrent}
+\end{figure}
+
+\begin{figure}
+  \hfil\includegraphics[width=8cm]{decay.pdf}
+  \caption{
+    The convergence rate to the asymptote of the average current as a function of $\frac{t\omega}{2\pi}$ on a log-log plot for $\omega=6$ 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 one-photon 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:
+\begin{equation}
+  \left<\left<j\right>\right>:=\frac1\tau\int_{T-\tau}^T dt\ \left<j\right>_t
+  .
+\end{equation}
+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 semi-classical 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
+
+\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.
+  }
+  \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$$
+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 
+\begin{equation}\label{4p1}
+\psi(x,t)\sim {\rm e}^{\frac 12 i k^2 t} \bar{\psi}(x,t)
+\end{equation}
+where $\bar{\psi}(x,t)$ is periodic in $t$ with period $2\pi/\omega$ and $\frac 12  k^2$ is the energy of the electrons in the incoming beam.
+This corresponds to the poles of $\hat{\psi}(p,t)$ on the imaginary axis being at 
+\begin{equation}\label{4p2}
+p=-\frac 12  i k^2+i\omega n,\ \ \ n\mathrm{\ \ integer}
+\end{equation}
+\bigskip
+
+\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_{n-m}^{(\kappa_m)}e^{-\kappa_mx}\mathcal T_m
+    &\mathrm{for\ }x>0
+  \end{array}\right.
+\end{equation}
+where
+\begin{equation}
+  \kappa_m=\sqrt{2U-k^2+\frac{E^2}{2\omega^2}-2m\omega}
+  \label{kappa}
+\end{equation}
+and
+\begin{equation}
+  g_{n-m}^{(\kappa_m)}=\frac\omega{2\pi}\int_0^{\frac{2\pi}\omega}dt\ e^{-i(n-m)\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
+
+\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.
+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}$.
+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$.
+The larger $m$ values necessary for smaller $\omega$ are difficult to see.
+\bigskip
+
+\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.
+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
+
+
+
+\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. FA9500-16-1-0037. O.C. was partially supported by the NSF-DMS (Grant No. 1515755). I.J. was partially supported by the NSF-DMS (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}
+
+
+\end{document}
+