I upload some library and initialize settings
In [1]:
from pylab import *
from copy import deepcopy
from matplotlib import animation, rc
from IPython.display import HTML
%matplotlib inline
rc('text', usetex=True)
font = {'family' : 'normal',
'weight' : 'bold',
'size' : 15}
matplotlib.rc('font', **font)
For calculating the occupation probability we can diagonalize the hamiltonian and compute $\Psi(t) as$:
\begin{equation*} \Psi(t) = D^\dagger e^{-i\Lambda t /\hbar} D \Psi(0) \end{equation*}Where $e^{-i\Lambda t /\hbar}$ is a diagonal matrix with on the diagonal $e^{-i\lambda_i t /\hbar}$ for each eigenvalue.
In [2]:
E1, E2, E3 = 0., 20., 0.
V12, V23 = 1., 1.
psi0 = array([1, 0, 0], dtype='complex')
Nt = int(1e4)
psi = zeros((Nt, 3), dtype='complex')
psi[0, :] = psi0
for E2, tf in zip(arange(4) * 20, [20, 200, 200, 200]):
times = linspace(0, tf, Nt)
H = array([[E1, V12, 0],
[V12, E2, V23],
[0, V23, E3]])
lambd, Q = eigh(H)
Q_inv = Q.T.conj()
for it in range(1, Nt):
psi[it, :] = Q_inv @ psi0
psi[it, :] = diag(np.exp(-1j * lambd * times[it])) @ psi[it, :]
psi[it, :] = Q @ psi[it, :]
plot(times, abs(psi) ** 2)
ylabel(r'$\|\Psi(t)\|^2$')
xlabel(r'$t$')
legend(['$\|\Psi(t)_1\|^2$', '$\|\Psi(t)_2\|^2$', '$\|\Psi(t)_3\|^2$'], loc=1)
figure()
In [3]:
y = cos(V12 ** 2 / E2 * times) ** 2
plot(times, y)
y = sin(V12 ** 2 / E2 * times) ** 2
plot(times, y)
plot(times, abs(psi[:, 0]) ** 2, label='$\|\Psi(t)_1\|^2$')
plot(times, abs(psi[:, 2]) ** 2, label='$\|\Psi(t)_3\|^2$')
ylabel(r'$\|\Psi(t)\|^2$')
xlabel(r'$t$')
legend(loc=1)
figure()
plot(times, abs(psi[:, 0]) ** 2, label='$\|\Psi(t)_1\|^2$')
y = cos(V12 ** 2 / E2 * times) ** 2
plot(times, y, label=r'$\cos(V_{12}^2 / (E_2 t)^2)$')
ylabel(r'$\|\Psi(t)\|^2$')
xlabel(r'$t$')
legend(loc=1)
ylim([.99, 1.01])
xlim([-.3, 3]);
We can discretize the hamiltonian with the finite difference scheme and compute the eigenvalues and eigenvector. I plotted the eigenfunctions for the 3 lowest energy eigenstates and in the second plot you can see the different bound energies compered to the potential.
In [4]:
def V(x, Z=1):
return -Z / sqrt(2 / Z ** 2 + x ** 2)
N = 2 ** 10
x0, x1 = -25, 25
x = linspace(x0, x1, N)
dx = (x1 - x0) / (N - 1)
H = diag(ones(N - 1), -1) - 2 * diag(ones(N)) + diag(ones(N - 1), 1)
H *= -1 / (2 * dx**2)
H += diag(V(x))
E, Psi_tot = eigh(H)
E_bound=E[E<0]
for k, E_ in enumerate(sorted(E_bound)[:3]):
print('E_{' + str(k) + '} = ' + "{:1.4f}".format(E_))
plot(x, Psi_tot[:, 0] / sqrt(dx), label=r'$\Psi_0(x)$')
plot(x, Psi_tot[:, 1] / sqrt(dx), label=r'$\Psi_1(x)$')
plot(x, Psi_tot[:, 2] / sqrt(dx), label=r'$\Psi_2(x)$')
legend(loc=1)
xlabel('x')
ylabel('$\Psi(t)$')
figure()
plot(x, V(x))
plot(x, E_bound * ones_like(x)[:, newaxis])
legend([r'$V(x)$', r'$E_0$', r'$E_1$', r'$E_2$'])
xlabel('x')
ylabel('Energy')
Out[4]:
In [5]:
def E(t, E0, omega, n):
t_ = maximum(omega * t, 0)
t_ = minimum(t_, 2 * np.pi * n)
return E0 * sin(t_) * sin(t_ / (2 * n))
def A(t, E0, omega, n):
pref = -E0 / omega
t_ = maximum(omega * t, 0.)
t_ = minimum(t_, 2 * np.pi * n)
return pref * (cos(t_) * (n * n * cos(t_ / n) - n * n + 1) + n * sin(t_) *
sin(t_ / n) - 1) / (2 * (n * n - 1))
def vanish(V0, x, x0, x1):
V0 *= 2
xs, xe = x[0], x[-1]
potential = np.maximum(0, (V0 * (x - x0) / (xs - x0)))
return np.maximum(potential, (V0 * (x - x1) / (xe - x1)))
In [6]:
omega = .02
n = 3 / 2
E0 = .05
Z = 1
x0, x1 = -15, 15
dx = .1
x = arange(x0, x1, dx)
N = len(x)
p = fftfreq(N, d=dx / (2 * pi))
dt = 0.5
ts = np.arange(- pi / omega, 2 * np.pi * (n + .5) / omega, dt)
plot(ts, E(ts, E0, omega, n))
title('Electric field for n = 3/2')
xlabel('t')
ylabel('E(t)')
figure()
t_star = np.arange(-pi / omega, 2 * np.pi * (5 + .5) / omega, 0.01)
plot(t_star, E(t_star, E0, omega, 5))
xlabel('t')
ylabel('E(t)')
title('Electric field for n = 5')
figure()
plot(ts, A(ts, E0, omega, n))
xlabel('t')
ylabel('A(t)')
title('Magnetic potential field for n = 3/2')
figure()
plot(t_star, A(t_star, E0, omega, 5))
xlabel('t')
ylabel('A(t)')
title('Magnetic potential field for n = 5')
Out[6]:
We know that in our simulation we work with a system with low ionation probability, this suggests that the wave function does not move much from the initial position. All of this would indicate that the kinetic momentum is low, so the optimal gauge is the length gauge, as the velocity gauge would shift the momentum by $A(t)$.
We can approximate our system to an electron in a sinusoidal electric field in the classical regime. Newton's equation then is : \begin{equation*} m\ddot{x} = E_0 \sin(\omega t) sin^2(\frac{\omega t}{2 n}) \end{equation*} by integrating this equation we obtain: \begin{equation*} m\dot{x} = \frac{E_0}{\omega} \frac{\cos(\omega t) (n^2 \cos(\frac{\omega t}{n}) - n^2 + 1) + n \sin(\omega t) \sin(\frac{\omega t}{n}) -1}{2 (n^2 - 1)} \end{equation*} looking at the possible maxima of this function one can put an upper bound on the momentum in the case of $n \leq 1$: \begin{equation*} \tilde{p} \leq (n^2 E_0) / ((n^2 -1) \omega) \end{equation*} if a general bound is needed one can use $((n^2 + 1) E_0)$ in place of the numerator. Obviously this is just an upper bound on the momentum and not a maximum, one can use a better value depending on the parameters, by calculating the maximum numerically.
Having an approximation for the maximum momentum we estimate the spacing of the spacial grid and temporal step with: \begin{equation*} \Delta x \leq \frac{\hbar\pi}{\tilde{p}} \quad \text{and} \quad \Delta t \leq \frac{\hbar\pi}{\tilde{E}} = \frac{2\pi \hbar}{\tilde{p}^2} \end{equation*} For the simulation we need also a maximum and minimum for $x$. Since we do not expect our wave function to move dramatically we can look and the ground state of the soft-core potential and look when the function is sufficiently close to 0. In the simulation we opted for -15 a.u., 15 a.u. as boundery.
Lastly we will choose the mean of the vanishing potential $\tilde{V}$ based around the size of the damping region $d$ and the average momentum $p$, by using: \begin{equation*} \frac{\hbar p}{4m_0d} \ll \tilde{V} \ll \frac{d p^3}{2 m_0 \hbar} \end{equation*}
In this first part we are going to set up the various parameters needed for the simulation. For $\tilde p$ we could also decide to use the maximum momentum obtained numerically, but we decided to be consistent with what was said previously.
In [7]:
omega = .02
n = 3 / 2
E0 = .05
Z = 1
x0, x1 = -15, 15
xl, xr = -10, 10
d = x1 - xr
t_temp = np.linspace(0, 2 * np.pi * (n + .5) / omega, 1000)
A_max = max(A(t_temp, E0, omega, n)) # the maximum momentum is equal to the
# maximum value of the magnetic potential
p_tilde = n**2 * E0 /(n**2 - 1) / omega
print('dx using the approximation ',\
"{:1.4f}".format(pi / p_tilde), 'a.u.')
print('dx using the maximum of the momentum calculated numerically',\
"{:1.4f}".format(pi / A_max), 'a.u.')
print('dt using the approximation ',\
"{:1.4f}".format(2 * pi / p_tilde ** 2), 'a.u.')
print('dt using the maximum of the momentum calculated numerically',\
"{:1.4f}".format(2 * pi / A_max ** 2), 'a.u.')
print("{:1.4f}".format(p_tilde / (8 * d)), 'a.u. < tilde_V <',\
"{:1.4f}".format(p_tilde ** 3 / 2 ** 4), 'a.u.')
V_tilde = 5.
dx = pi / p_tilde
x = arange(x0, x1, dx)
N = len(x)
p = fftfreq(N, d=dx / (2 * pi))
dt = 2 * pi / p_tilde ** 2
ts = np.arange(0, 2 * np.pi * (n + .5) / omega, dt)
In [8]:
H = diag(ones(N - 1), -1) - 2 * diag(ones(N)) + diag(ones(N - 1), 1)
H *= -1 / (2 * dx ** 2)
H += diag(V(x, Z))
U_2 = exp(-1j * 0.5 * p ** 2 * dt)
_, Psi_tot = eigh(H)
Psi = Psi_tot[:, 0].astype('complex')
Psi /= np.sqrt(sum(abs(Psi) ** 2 * dx))
psi0 = deepcopy(Psi)
norm = np.zeros(len(ts))
overlap = np.zeros(len(ts))
In [9]:
for k, t in enumerate(ts):
U_1 = exp(-0.5 * 1j * (V(x, 1) - 1j *
vanish(V_tilde, x, xl, xr) - x * E(t, E0, omega, n)) * dt)
Psi *= U_1
Psi = fft(Psi)
Psi *= U_2
Psi = ifft(Psi) # go to real space
Psi *= U_1
norm[k] = sum(abs(Psi) ** 2 * dx)
overlap[k] = abs(vdot(Psi, psi0)) * dx
In [11]:
N_e = 20
ionizs = np.zeros(N_e)
norms = np.zeros((N_e, len(ts)))
for j, E0 in enumerate(np.linspace(0, .05, N_e)):
Psi = deepcopy(psi0)
for k, t in enumerate(ts):
U_1 = exp(-0.5 * 1j * (V(x, 1) - 1j * vanish(V_tilde, x, -10, 10) - x * E(t, E0, omega, n)) * dt)
Psi *= U_1
Psi = fft(Psi)
Psi *= U_2
Psi = ifft(Psi) # go to real space
Psi *= U_1
norms[j, k] = sum(abs(Psi) ** 2 * dx)
ionizs[j] = 1 - sum(abs(Psi) ** 2 * dx)
We see as expected that the lower the amplitude of the electric field the lower the final ionization probability will be.
In [18]:
title('Ionization probabilies in time')
plot(ts, 1 - norms.T[:, ::-1])
legend([r'$E_0 = 0.05$ a.u.'])
xlabel(r't')
ylabel(r'$1 - |<\Psi(t)|\Psi(t)>|^2$')
figure()
title(r'ionization probabilies at $t_{end}$')
ylabel(r'$1 - |<\Psi(t_{end})|\Psi(t_{end})>|^2$')
xlabel(r'$E_0$')
plot(np.linspace(0, .05, N_e), ionizs)
Out[18]:
In [ ]: