\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}

