Ian Jauslin
summaryrefslogtreecommitdiff
blob: e0ccaabdcadb4d4bde3b1b468bf13c5f92997a79 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
\documentclass{ian}

\usepackage{graphicx}
\usepackage{array}

\begin{document}

\hbox{}
\hfil{\bf\LARGE
Solution of the time dependent Schr\"odinger equation\par
\vskip10pt
\hfil leading to Fowler-Nordheim field emission\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{\it School of Mathematics, Institute for Advanced Study}\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
 We solve the time-dependent Schr\"odinger equation describing the emission of electrons from a metal surface by an external electric field $E$, turned on at $t=0$. Starting with a wave function $\psi(x,0)$, representing a generalized eigenfunction when $E=0$, we find $\psi(x,t)$ and show that it approaches, as $t\to\infty$, the Fowler-Nordheim tunneling wavefunction $\psi_E$. The deviation of $\psi$ from $\psi_E$ decays asymptotically as a power law $t^{-\frac32}$. The time scales involved for typical metals and fields of several V/nm are of the order of femtoseconds. We plot the short-time evolution of the current and density.

\vfill

\tableofcontents

\vfill
\eject

\setcounter{page}1
\pagestyle{plain}


\section{Introduction}

\indent The emission of electrons from a cold metal surface subjected to a constant (or oscillating) electric field is a subject of great practical and theoretical interest\-~\cite{Je03,Fo16,Je17}. The microscopic theory of such emissions by a constant field was developed by Fowler and Nordheim (FN) in the early days of quantum mechanics\-~\cite{FN28} (referred to then as the ``new mechanics''). They  considered an idealized situation in which the electrons  in the conduction band are treated, a la Sommerfeld, as free independent particles. Their energies are described by a Fermi distribution with maximum energy $E_F=\hbar^2k_{\mathrm F}^2/2m$; the deviation from this zero-temperature distribution is negligible at room temperatures. In the absence of an external field the electrons are confined by an external potential (caused by the positive ions) of magnitude $U=E_F+W$, where $W$ is the work function, i.e. the energy necessary to extract an electron from the metal.
\bigskip

\indent Considering emissions perpendicular to a flat surface at $x=0$, obtained when applying an external field $E$ for $x\geqslant 0$, assuming that the metal occupies all space $x<0$, leads to a one-dimensional tunneling problem in a triangular potential, see Fig.\-~\ref{fig:potential}. The one-dimensional Schr\"odinger equation describing an electron moving in this potential is then given by
\begin{equation}
  \label{schrodinger}
  i\partial_t\psi(x,t)=(-{\textstyle\frac12}\partial_x^2+V(x))\psi(x,t)
\end{equation}
(we write $\partial_x\equiv\frac\partial{\partial x}$) where
\begin{equation}
  V(x)=\left\{ \begin{array}{ll} 0, & x<0\\
  U-E x, & x>0 \end{array}\right.  
  \label{V}
\end{equation}
in atomic units $(\hbar=m=|e|=1)$.
\bigskip

\begin{figure}
  \hfil\includegraphics[width=8cm]{potential.pdf}\hfil\hbox{}\par
  \caption{The shape of the potential $V(x)$.}
  \label{fig:potential}
\end{figure}

\indent When $E=0$, the potential is, simply, a step function. The Schr\"odinger equation\-~(\ref{schrodinger}) with $E=0$ has stationary solutions with energies $k^2/2<U$, $\psi(x,t)=e^{-i\frac12k^2t}\psi_0(x)$, with $k>0$ and
\begin{equation}
  \psi_0(x)=\left\{ \begin{array}{ll} e^{i k x} +R_0 e^{-i k x} & x<0\\
  T_0e^{-\sqrt{2U-k^2}x} & x>0 \end{array}\right.  
  \label{psi0}
\end{equation} 
in which $R_0$ and $T_0$ are the {\it reflection} and {\it transmission} coefficients (we use a normalization in which the amplitude of the incoming wave with $k>0$ is 1):
\begin{equation}
  R_0=\frac{ik+\sqrt{2U-k^2}}{ik-\sqrt{2U-k^2}}
  ,\quad
  T_0=\frac{2ik}{ik-\sqrt{2U-k^2}}
  .
  \label{R0T0}
\end{equation}
These constants ensure that $\psi_0(x)$ and $\partial_x\psi_0(x)$ are continuous at $x=0$. Note that, in this state, the current vanishes:
\begin{equation}
  j_0(x)=i(\psi_0\partial_x\psi_0^*-\psi_0^*\partial_x\psi_0)=0
  .
\end{equation}
\bigskip

\indent When $E>0$, there is the possibility for an electron moving in the $+x$ direction, with kinetic energy $k^2/2<U$, to tunnel through the potential barrier and be emitted. This will then produce an electron current in the $+x$-direction. To obtain the probability of tunneling, FN computed the stationary solutions $\psi(x,t)=e^{-i\frac12k^2t}\psi_E(x)$ by solving
\begin{equation}
  (-{\textstyle\frac12}\partial_x^2+\Theta(x)(U-Ex)-{\textstyle\frac12}k^2)\psi_E(x)=0
  \label{eqpsiE}
\end{equation}
($\Theta(x)$ is the Heaviside function, which is equal to 1 if $x\geqslant 0$ and $0$ otherwise) whose solution is
\begin{equation}
  \psi_E(x)=
  \left\{ \begin{array}{ll}
    e^{i k x} +R_E e^{-i k x} & x<0\\
    T_E\Phi(x) & x>0.
  \end{array}\right.  
  ,\quad
  k>0
  \label{psiE}
\end{equation} 
in which $\Phi(x)$ is proportional to the Airy  function Ai$(x)$ (or the equivalent expression in terms of Hankel or Bessel functions), which decays when $x\to\infty$, and yet has a constant positive current for all $x$. This solution, see also\-~\cite{Ro11,Je03}, yielded the tunneling probability $D(k)=1-|R_E|^2$ of the electron as a function of $k,\,U$ and $E$. Integrating $kD(k)$ over the ``supply function'' corresponding to the density of electrons in the Fermi sea moving in the $+x$ direction with energy $k^2/2$, leads to an expression for the total steady state current $j_E$ in a static field $E$. An approximate expression for $j_E$ is~\-\cite{Fo08b,Je17}
\begin{equation}
  j_E\approx c_1E^2e^{-\frac{c_2}E}
  .
\end{equation}
\bigskip

\indent The FN formula for $j_E$, with various corrections for the idealizations made, e.g. flat surface, independent electrons, neglecting the Schottky effect, {\it etc.}, serves as the backbone of cold electron emission theory and experiment.  There is a vast literature on the subject (the original FN paper\-~\cite{FN28} has more than 6000 citations). We cite here only a few\-~\cite{Ro11,Fo08} and refer the reader for more information to the recent book by Jensen\-~\cite{Je17} and references therein.
\bigskip

\indent In this note we shall be concerned with a different problem, which, as far as we know, has not been investigated fully before. As an initial condition, we take a stationary solution of the Schr\"odinger equation at $E=0$, $\psi_0(x)$ in\-~(\ref{psi0}), and, at $t=0$, we turn the field on, and study the time evolution. In particular, we will investigate how long it will take, if ever, for the initial state $\psi(x,0)$ to approach the stationary state  $\psi_E(x)$ in\-~(\ref{psiE}). Of course, turning on $E$ instantaneously is an idealization, which we shall accept here. (In\-~\cite{YGR11}, this initial condition is considered, but the analysis then focuses mostly on the stationary solution.)
\bigskip

\indent In what follows, we shall prove that, for $\psi(x,0)=\psi_0(x)$, $\psi (x,t)$ approaches, for long times, the $\psi_E(x)$ of\-~(\ref{psiE}), i.e.,
\begin{equation}
  \label{psi_lim}
  \psi (x,t)\sim e^{-i\frac12k^2t}\psi_E(x)
\end{equation}
In fact, this holds for a wider class of initial conditions, in which the initial incident wave is $e^{ikx}$ and the initial reflected and transmitted waves are arbitrary. The deviation $\psi(x,t)-\psi_E(x)$ decays asymptotically as $t^{-\frac32}$. The actual time dependence, of course, depends on the exact form of $\psi(x,0)$. We shall calculate this for the $\psi(x,0)=\psi_0(x)$ given in\-~(\ref{psi0}) for different values of the parameters.
\bigskip

\indent Roughly speaking we find that for $U\approx9\ \mathrm{eV}$, $\hbar^2k_{\mathrm F}^2/2m=E_{\mathrm F}\approx 4.5\ \mathrm{eV}$ and $E\approx 4$-$8\ \mathrm{V}\cdot\mathrm{nm}^{-1}$, the time for the density $|\psi|^2$ and the current $j(t)$ to approach its final FN value is of the order of femtoseconds. The exact value depends on the position $x$ where we measure the current: for larger $x$, the time it takes for the current to stabilize is larger, see Fig.\-~\ref{current_density}. Such time scales are of practical relevance for short pulses of the order of femtoseconds or less. These are now common for oscillating laser fields for which the initial value problem will be considered in a later paper. (The ``steady state'' solution for laser fields was investigated in detail by Faisal et al\-~\cite{FKS05}; see also \cite{ZL16}.)

\section{Solution of the initial value problem}

\indent In order to emphasize the role of each term in the initial condition, we will split $\psi(x,0)$ into three terms: an incoming, a reflected, and a transmitted wave.
\begin{equation}
  \psi(x,0)=\psi^{(\mathrm I)}(x,0)+\psi^{(\mathrm R)}(x,0)+\psi^{(\mathrm T)}(x,0)
  \label{psi_sum}
\end{equation}
with
\begin{equation}
  \psi^{(\mathrm I)}(x,0)=\Theta(-x)e^{ikx}
  ,\quad
  \psi^{(\mathrm R)}(x,0)=R_0\Theta(-x)e^{-ikx}
  ,\quad
  \psi^{(\mathrm T)}(x,0)=T_0\Theta(x)e^{-\sqrt{2U-k^2}x}
   ,\quad
   k>0
\end{equation}
(recall that $\Theta(x)$ is the Heaviside function, which is equal to 1 if $x\geqslant 0$ and $0$ otherwise). Since the Schr\"odinger equation is linear, its solution will be the sum of the solutions for each term in $\psi(x,0)$.
\bigskip

\indent To obtain $\psi(x,t)$ we solve for $\hat\psi_p(x)$, the Laplace transform of $\psi(x,t)$,
\begin{equation}
  \hat\psi_p(x):=\int_0^\infty dt\ e^{-pt}\psi(x,t)
  \label{psip}
\end{equation}
which we obtain in closed form. We then compute, by inverting the Laplace transform, the long time asymptotics analytically, and the short time behavior numerically. This method provides an integral representation of the solution which can be evaluated numerically. It is thus better for our purposes than direct computations of the solution of\-~(\ref{schrodinger}). The latter requires cutoffs for the non-square integrable functions we are dealing with and cannot be used for long times. The Laplace transform of $\psi$ satisfies the equation
\begin{equation}
  (-{\textstyle\frac12}\partial_x^2+\Theta(x)(U-Ex)-ip)\hat\psi_p(x)=-i\psi(x,0)
  .
  \label{schrodinger_laplace}
\end{equation}
The physical solution to this equation is
\begin{equation} \label{sol11}
  \hat\psi_p(x)=
  \left\{\begin{array}{>\displaystyle ll}
    C_1(p)e^{\sqrt{-2ip}x}+F^{(\mathrm I)}_p(x)+R_0F^{(\mathrm R)}_p(x)
    &\mathrm{if\ }x<0\\[0.5cm]
    C_2(p)\varphi_p(x)+T_0F^{(\mathrm T)}_p(x)
    &\mathrm{if\ }x> 0
  \end{array}\right.
\end{equation}
where $R_0$ and $T_0$ are given in\-~(\ref{R0T0}),
\begin{equation}\label{2p6}
  F^{(\mathrm I)}_p(x):=-\frac{2ie^{ikx}}{-2ip+k^2}
  ,\quad
  F^{(\mathrm R)}_p(x):=-\frac{2ie^{-ikx}}{-2ip+k^2}
\end{equation}
\begin{equation}
  F^{(\mathrm T)}_p(x):=
  \frac{4\pi}{(2E)^{\frac13}}\left(\varphi_p(x)\int_0^x dy\ \eta_p(y)e^{-\sqrt{2U-k^2}y}
  +\eta_p(x)\int_x^\infty dy\ \varphi_p(y)e^{-\sqrt{2U-k^2}y}\right)
\end{equation}
and
\begin{equation}
  \varphi_p(x)=\mathrm{Ai}\left(2^{\frac13}e^{-\frac{i\pi}3}\left(E^{\frac13}x-E^{-\frac23}(U-ip)\right)\right)
\end{equation}
\begin{equation}
  \eta_p(x)=e^{-\frac{i\pi}3}\mathrm{Ai}\left(-2^{\frac13}\left(E^{\frac13}x-E^{-\frac23}(U-ip)\right)\right)
  \label{varphi}
\end{equation}
are two independent solutions of $(-\frac12\partial_x^2+U-Ex-ip)f=0$. The phases $e^{-\frac{i\pi}3}$ and $-1$ are cube roots of $-1$. The constants $C_1(p)$ and $C_2(p)$ are set so that $\hat\psi_p$ and $\partial_x\hat\psi_p$ are continuous at $x=0$:
\begin{equation}
  C_1(p)=-\frac{2iT_0}{\sqrt{-2ip}\varphi_p(0)-\partial\varphi_p(0)}\left(
    \frac{\sqrt{2U-k^2}\varphi_p(0)+\partial\varphi_p(0)}{-2ip+k^2}
    +\int_0^\infty dy\ \varphi_p(y)e^{-\sqrt{2U-k^2}y}
  \right)
\end{equation}
and
\begin{equation}\label{formC2}
  \begin{array}{>\displaystyle r@{\ }>\displaystyle l}
    C_2(p)=-\frac{2iT_0}{\sqrt{-2ip}\varphi_p(0)-\partial\varphi_p(0)}&\Bigg(
      \frac{\sqrt{2U-k^2}+\sqrt{-2ip}}{-2ip+k^2}
      \\\hfill&\hskip10pt
      -\frac{2i\pi}{(2E)^{\frac13}}(\sqrt{-2ip}\eta_p(0)-\partial\eta_p(0))\int_0^\infty dy\ \varphi_p(y)e^{-\sqrt{2U-k^2}y}
    \Bigg)
  \end{array}
\end{equation}
where $\partial\varphi_p(0)\equiv\frac{\partial\varphi_p(x)}{\partial x}\big|_{x=0}$ and similarly for $\partial\eta_p(0)$. The square root is defined with a branch cut along the positive imaginary axis, in such a way that $\sqrt{-2ip}$ has a branch cut along the real negative axis.
\bigskip

\indent A simple calculation shows that, as expected,
\begin{equation}
  \lim_{\displaystyle\mathop{\scriptstyle|p|\to\infty}_{\mathcal Re(p)>0}}p\hat\psi_p(x)=\psi(x,0)
  \label{laplace_init}
\end{equation}
which confirms that $\hat\psi_p(x)$ is, indeed, the Laplace transform of a function whose initial condition is $\psi(x,0)$.
\bigskip

\indent We then invert the Laplace transform:
\begin{equation}
  \psi(x,t)=\frac1{2i\pi}\int_{\gamma-i\infty}^{\gamma+i\infty}dp\ e^{pt}\hat\psi_p(x)
  \label{inv_laplace}
\end{equation}
in which $\gamma>0$ is an arbitrary small parameter taken close to $0$. 
\bigskip

\indent As is well known the integral on the right hand side of\-~(\ref{inv_laplace}) can be computed deforming the integration contour as in Fig.\-~\ref{fig:contour}, and studying the singularities, poles and branch points of $\hat\psi_p(x)$, lying in the half plane $\mathcal Re(p)\leq 0$. In particular, the only terms which do not decay as $t\to\infty$ come from poles on the imaginary $p$-axis. Analyzing\-~(\ref{sol11})-(\ref{formC2}) we find that the singularities of $\hat\psi_p(x)$ are, for $k>0$,
\begin{itemize}
  \item a pole on the imaginary axis, located at $-ik^2/2$, coming from  (\ref{2p6}), (\ref{formC2}) and (\ref{laplace_init})
  \item  poles with strictly negative real parts corresponding to the roots of $\sqrt{-2ip}\varphi_p(0)-\partial\varphi_p(0)$ appearing in the denominators of $C_1$ and $C_2$,
  \item a branch cut along the negative real axis coming from $\sqrt{-2ip}$.
\end{itemize}
\bigskip

\point{\bf Long time behavior.} The residue at $-ik^2/2$ yields the only term which does not decay in time: by an explicit computation, we find that the residue is equal to
\begin{equation}
  e^{-i\frac12k^2t}\psi_E(x)
\end{equation}
where $\psi_E$ is the FN solution\-~(\ref{psiE}).
\bigskip

\indent The residues of the poles with a negative real part decay exponentially in time (because of the factor $e^{pt}$ in~\-(\ref{inv_laplace})).
\bigskip

\indent The integral along the branch cut decays algebraically, as $t^{-\frac32}$: we define, for $p-i\epsilon\in\mathbb R_-$,
\begin{equation}
  \alpha:=e^{\frac{i\pi}4}\sqrt{-ip}
  ,\quad
  f(\alpha):=\hat\psi_{\alpha^2}(x)
\end{equation}
(recall the definition of $\hat\psi_p$ in~\-(\ref{psip})) and write the integral along the branch cut as
\begin{equation}
  \psi^{(\mathrm{BC})}(x,t)
  :=
  \int_{-\infty-i\epsilon}^{-i\epsilon} dp\ e^{pt}\hat\psi_p(x)
  +
  \int_{i\epsilon}^{-\infty+i\epsilon} dp\ e^{pt}\hat\psi_p(x)
  =
  2\int_0^\infty d\alpha\ e^{-\alpha^2t}\alpha(f(\alpha)-f(-\alpha))
  .
\end{equation}
By Taylor expansion, (in this context, this technique is usually called Watson's lemma)
\begin{equation}
  \psi^{(\mathrm{BC})}(x,t)
  =
  4\int_0^\infty d\alpha\ e^{-\alpha^2t}\alpha^2\partial f(0)
  +O(t^{-\frac52})
  =\left(\frac{t}{\tau_E(x)}\right)^{-\frac32}+O(t^{-\frac52})
  .
\end{equation}
with
\begin{equation}
  \tau_E(x)=\left\{\begin{array}{>\displaystyle ll}
    \left(c_E(\varphi_0(0)+x\partial\varphi_0(0))\right)^{\frac23}
    &\mathrm{if\ }x<0
    \\[0.5cm]
    \left(c_E\varphi_0(x)\right)^{\frac23}
    &\mathrm{if\ }x>0
  \end{array}\right.
\end{equation}
and
\begin{equation}
  c_E=
  -\frac{\sqrt2T_0e^{\frac{i\pi}4}}{\sqrt\pi(\partial\varphi_0(0))^2}
  \left(
    \frac{\sqrt{2U-k^2}\varphi_0(0)+\partial\varphi_0(0)}{k^2}
    +\int_0^\infty dy\ \varphi_0(y)e^{-\sqrt{2U-k^2}y}
  \right)
  .
\end{equation}

\indent All in all, we find that
\begin{equation}
  \psi(x,t)
  =e^{-i\frac12k^2t}\psi_E(x)+\left(\frac{t}{\tau_E(x)}\right)^{-\frac32}+O(t^{-\frac52})
  .
\end{equation}
Therefore, the wave function tends to the Fowler-Nordheim solution, with a rate $t^{-\frac32}$.
\bigskip

\begin{figure}
  \hfil\includegraphics[width=8cm]{contour.pdf}\hfil\hbox{}\par
  \caption{The deformed integration contour goes around the poles (one of which is on the imaginary axis, at $-ik^2/2$, while the others are in the negative real half-plane) and goes along the branch cut on the real negative axis.}
  \label{fig:contour}
\end{figure}

\point{\bf Short time behavior.} The behavior of $\psi(x,t)$ for small $t$ is more difficult to study analytically, but the inverse Laplace transform\-~(\ref{inv_laplace}) yields an integral formula that can be efficiently approximated numerically using fast Fourier transforms.
\bigskip

In Fig.\-~\ref{current_density} we have plotted the density $|\psi(x,t)|^2$, current
\begin{equation}
  j_k(x,t):=i(\psi\partial_x\psi^*-\psi^*\partial_x\psi)
\end{equation}
and integrated current (the current integrated over the supply function at 0 temperature)
\begin{equation}
  J_{k_{\mathrm F}}(x,t):=\int_0^{k_{\mathrm F}}dk\ j_k(x,t)
\end{equation}
as a function of time at two different values of $x$: $x_0:=\frac{2U-k_{\mathrm F}^2}{2E}\approx 11\ \mathrm{nm}$ and $10x_0$ ($x_0$ is the point at which $V(x_0)=\frac{k_{\mathrm F}^2}2$), and at two different values of $E$: $4$ and $8\ \mathrm{V}\cdot\mathrm{nm}^{-1}$. We have normalized the current $j$ by $2k$, which is the current of the incoming wave $e^{ikx}$, and the integrated current $J$ by $k_{\mathrm F}^2$, which is the current of the incoming wave integrated over the supply function. We find that there is a transient regime that lasts a few femtoseconds before the system stabilizes to the FN value. Note that the approach to the FN regime has some ripples, which come from the imaginary parts of the poles in the $p$-plane (see Fig.\-~\ref{fig:contour}). There is a delay before the signal reaches $x_0$, and between $x_0$ and $10x_0$. As expected, the asymptotic value of the current is independent of $x$. Note that the current and density depend strongly on the field $E$.
\bigskip

\begin{figure}
  \begin{tabular}{c|c}
    \includegraphics[width=7cm]{density-4.pdf}           {\scriptsize({\bf a})}&
    \includegraphics[width=7cm]{density-8.pdf}           {\scriptsize({\bf b})}\\
    \hline
    \includegraphics[width=7cm]{currents-4.pdf}          {\scriptsize({\bf c})}&
    \includegraphics[width=7cm]{currents-8.pdf}          {\scriptsize({\bf d})}\\
    \hline
    \includegraphics[width=7cm]{integrated_current-4.pdf}{\scriptsize({\bf e})}&
    \includegraphics[width=7cm]{integrated_current-8.pdf}{\scriptsize({\bf f})}
  \end{tabular}
  \caption{The density ({\bf a}),({\bf b}), current ({\bf c}),({\bf d}) and integrated current ({\bf e}),({\bf f}) as a function of time at $x=x_0\equiv\frac{2U-k_{\mathrm F}^2}{2E}$ and $x=10x_0$. We have taken $U=9\ \mathrm{eV}$ and $k^2/2=k_{\mathrm F}^2/2\equiv E_{\mathrm F}=4.5\ \mathrm{eV}$. In ({\bf a}),({\bf c}),({\bf e}), the field is $E=4\ \mathrm{V}\cdot \mathrm{nm}^{-1}$ and $x_0\approx11\ \mathrm{nm}$. In ({\bf b}),({\bf d}),({\bf f}), $E=8\ \mathrm{V}\cdot\mathrm{nm}^{-1}$ and $x_0\approx22\ \mathrm{nm}$. In ({\bf a}),({\bf c}),({\bf e}), the plots seem to indicate that the curves converge to 0, but they actually tend to small finite values.}
  \label{current_density}
\end{figure}

\noindent{\bf Remark}: While the time scale of the approach to the FN solution is clearly of order of femtoseconds we have not attempted to compute a ``tunneling time''. This is, as is well known, a tricky business, with many possible definitions, see, e.g. \cite{LM94,LK15}. Defining such a time in terms of the approach of the initial state to some steady state was investigated by McDonald et al.\-~\cite{MOe13}. Pfeifer and Fr\"ohlich\-~\cite{PF95} have computed rigorous bounds on the lifetimes of spatially confined states.

\section{Possible generalizations}  

\point{\bf Initial conditions.} As is shown by the computation described above, the long time asymptotic behavior of the wave function is independent of the initial reflected and transmitted waves. That is, the initial condition $\psi(x,0)=\psi^{(\mathrm I)}(x,0)=\Theta(-x)e^{ikx}$ leads to the same asymptotic formula:
\begin{equation}
  \psi(x,t)\sim e^{-i\frac12 k^2t}\psi_E(x).
\end{equation}
Indeed, the reflected and transmitted initial conditions do not actually give rise to any poles on the imaginary axis and therefore their contributions decay in time.
\bigskip

This leaves open the possibility to consider much more general initial conditions than\-~(\ref{psi_sum}): one can change the coefficients of the reflected and transmitted waves, add such waves with different wave vectors, or remove them altogether, without changing the asymptotic formula. Only the incoming wave $e^{ikx}$ affects it. In addition, one can add any square-integrable function to the initial condition without changing the long-time behavior. This is a consequence of the RAGE theorem \cite{Ru69,AG73,En78}, which states that whenever the Hamiltonian has absolutely continuous spectrum (as is the case here), the solution of the Schr\"odinger equation with a square-integrable initial condition vanishes point-wise as $t\to\infty$.
\bigskip

With this fact in mind, one can make an easy argument why the asymptotic behavior of the wave function must coincide with the stationary solution. Indeed, once we drop the initial reflected and transmitted waves,
 the Laplace transform of the wave function is of the form
\begin{equation}
  \hat\psi_p(x)=\frac{-2i}{-2ip+k^2}\left(\left(e^{ikx}+Re^{-\sqrt{-2ip}x}\right)\Theta(-x)+T\varphi_p(x)\Theta(x)\right)
\end{equation}
in which $\varphi_p$ is the solution of
\begin{equation}
  \left(-\frac12\partial_x^2+V(x)-ip\right)\varphi_p(x)=0
  ,\quad
  x>0
  .
  \label{eqvarphi}
\end{equation}
When $p=-ik^2/2$, (\ref{eqvarphi}) coincides with the equation\-~(\ref{eqpsiE}) for $\psi_E$. Assuming that $R$ and $T$ do not introduce any new poles on the imaginary axis (as we showed is the case), this implies that $\psi(x,t)$ converges to $e^{-i\frac12k^2t}\psi_E(x)$ as $t\to\infty$.
\bigskip

\point{\bf Potentials.} The exact form of the potential $V(x)=U-Ex$ was not really used in much of the computation above, so it can be carried out in very much the same way for many other $V(x)$. For instance, one could round off the triangular barrier as occurs in the Schottky effect\-~\cite{Fo08}. We could also consider a square barrier. The only real constraint on the potential is that it does not introduce bound states. To make this into a precise statement, one would also have to put constraints on the regularity and asymptotic properties of $V$, which we will not do here.
\bigskip

This leaves open the possibility of studying trains of pulses, in which the field is turned on and off repeatedly. The regime in which the field is off corresponds to a potential $V(x)=U\Theta(x)$, which can be studied using the method described above. Provided the time between the pulses is long enough, the system would stabilize to the stationary state in the time between each field switching.
\bigskip

\point{\bf Time-dependent fields.} It would be very interesting to consider a similar question in the case of an oscillating laser field $E=e_0\cos(\omega t)$. The stationary state of this problem was studied by Faisal et al.\-~\cite{FKS05}, and we are currently working on showing that the solutions of the initial value problem converge to this solution, and studying the short-time behavior.
\bigskip


\vfill

\delimtitle{\bf Acknowledgements}
\noindent This material is based upon work supported by the AFOSR under the award number FA9500-16-1-0037. OC was partially supported by the NSF-DMS grant 1515755. IJ was partially supported by the NSF-DMS grant 1128155. JLL thanks Kevin Jensen and Don Shiffler for useful discussions and the IAS for hospitality during part of this work.  
\enddelim

\vfill
\eject

\begin{thebibliography}{WWW99}
\small
\IfFileExists{bibliography/bibliography.tex}{\input bibliography/bibliography.tex}{}
\end{thebibliography}


\end{document}