diff options
Diffstat (limited to 'Costin_Costin_Jauslin_Lebowitz_2019.tex')
| -rw-r--r-- | Costin_Costin_Jauslin_Lebowitz_2019.tex | 468 |
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} + |