<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Computational Quantum Dynamics
(Lorenzo Biasi)
Project
I upload some library and initialize settings
Step1:  1. Superexchange in a three-level system.
(a)
For calculating the occupation probability we can diagonalize the hamiltonian and compute $\Psi(t) as$
Step2:  (b)
One can plot the given functions over the computation for $|\Psi_{1, 2}(t)|^2$ for values of $V = 1$ and $E_2  = 20$. For being able to distinguish the approximation from $|\Psi_{1, 2}(t)|^2$ we plot and enlargment for the the first coordinate.
Step3:  2. The one-dimensional soft-core potential.
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.
Step4:  3. Ionization from a one-dimensional soft-core potential.
(a)
Step5:  (b)
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)$. 
I will first  do (d) and then (c) and (e) together.
(d)
We can approximate our system to an electron in a sinusoidal electric field in the classical regime. Newton's equation then is 
Step6:  (f)
I run the previous program for different values of $E_0$.
Step7:  We see as expected that the lower the amplitude of the electric field the lower the final ionization probability will be.
<ASSISTANT_TASK:>
Python Code:
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)
Explanation: Computational Quantum Dynamics
(Lorenzo Biasi)
Project
I upload some library and initialize settings
End of explanation
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()
Explanation: 1. Superexchange in a three-level system.
(a)
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.
End of explanation
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]);
Explanation: (b)
One can plot the given functions over the computation for $|\Psi_{1, 2}(t)|^2$ for values of $V = 1$ and $E_2  = 20$. For being able to distinguish the approximation from $|\Psi_{1, 2}(t)|^2$ we plot and enlargment for the the first coordinate.
End of explanation
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')
Explanation: 2. The one-dimensional soft-core potential.
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.
End of explanation
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)))
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')
Explanation: 3. Ionization from a one-dimensional soft-core potential.
(a)
End of explanation
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)
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))
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
Explanation: (b)
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)$. 
I will first  do (d) and then (c) and (e) together.
(d)
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}
(c) and  (e)
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.
End of explanation
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)
Explanation: (f)
I run the previous program for different values of $E_0$.
End of explanation
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)
Explanation: We see as expected that the lower the amplitude of the electric field the lower the final ionization probability will be.
End of explanation
<END_TASK>