Equation de Schrödinger et états comprimés

Ce TD a pour but de simuler numériquement le TD “Etats quantiques d’atomes de césium dans un piège harmonique” du cours de physique quantique.

On prendra \(m=1\), \(\hbar=1\). On considère de plus un piège harmonique de pulsation \(\omega = 2\pi\).

1. Ecrire (avec un papier et un crayon!) l’équation de Schrödinger.

2. Ecrire (avec un papier et un crayon) la fonction d’onde de l’état fondamental du piège harmonique.

On va échantillonner l’espace entre -5 et 5 avec un pas de 0.02

3. Ecrire une fonction Python qui calcule un état Gaussien d’écart type \(\sigma_x\) de vitesse \(v_0\) et centré en \(x_0\).

4. Ecrire une fonction Python qui calcule le Laplacien d’une fonction d’onde \(\psi\).

5. Ecrire la fonction Python psi_prime(t, psi, V) qui donne la dérivée de psi. Vérifier graphiquement que l’on a bien \(\frac{\partial \psi}{\partial t} = -i\frac{\omega}{2} \psi\) pour l’état fondamental du piège.

6. Résoudre numériquement sur un temps \(t=1.5\) l’équation de Schrödinger pour différents états initiaux : par exemple un paquet centré en \(-2\) initialement au repos et de largeur \(0.3\). On utilisera la méthode RK23.

7. Ecrire une fonction python moments_x(psi) qui retourne \(\mathbb{E}(X)\) (position moyenne) et \(\mathbb{E}((X - \mathbb{E}(X))^2)\) (carré de la largeur du paquet)

8. Tracer ces moments en fonction de t.

9. Question 7 et 8, mais avec l’impulsion. Dans le calcul de l’espérance, on suprimera les parties hautes impulsions qui sont dûes à du bruit numérique

10. Pendant 0.5 seconde, on laisse la paquet évoluer librement (V=0) puis pendant une durée \(t<0.5\), le paquet évolue avec le potentiel harmonique. Tracer la largeur en impulsion en fonction de cette durée.

---

Solution

Equation de Schrödinger :

(1)\[\begin{equation} i\hbar\frac{\partial\psi}{\partial t} = -\frac{\hbar^2}{2m}\frac{\partial^2\psi}{\partial x} + \frac{1}{2} m \omega^2 x^2 \end{equation}\]

Etat fondamental : on définit \(a_\mathrm{oh} = \sqrt{\frac{\hbar}{m\omega}}\). L’état fondamental est donné par :

(2)\[\begin{equation} \psi_0 = \frac{1}{\left( \pi a^2_\mathrm{oh}\right)^{\frac{1}{4}}} e^{-\frac{x^2}{2a_\mathrm{oh}^2}} \end{equation}\]

a_ho = np.sqrt(hbar/(m*omega)) psi_0 = gaussian_packet(a_ho, 0, 0)

import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML

import numpy as np
from scipy import integrate

from numba import jit, njit

#Question 3

L = 10.
dx = 0.02
N_x = int(round(L/dx))

x = (np.arange(N_x)-N_x//2)*dx


m = 1
hbar = 1

def norm(psi):
    return np.sqrt(np.sum(dx*np.abs(psi)**2))


def gaussian_packet(sigma, x_0=0, v_0=0):
    k_x = m*v_0/hbar
    psi = np.exp(-(x-x_0)**2 / (2.0 * sigma**2)) * np.exp(1j * k_x * x)
    psi = psi / norm(psi)
    return psi

T = 1
omega = 2*np.pi / T

V = 0.5 * m * omega**2 * x**2

# Question 4

@jit
def Laplace(psi, dx=dx):
    output = np.zeros_like(psi)
    for i in range(1, N_x-1):
        output[i] = (psi[i-1] + psi[i+1]) - 2*psi[i]
    output[0] = psi[1] - 2*psi[0]
    output[-1] = psi[-2] - 2*psi[-1]
    return output/(dx**2)

psi = gaussian_packet(1.5, 10)
_ = Laplace(psi)

%timeit Laplace(psi)

1.78 µs ± 31.1 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

# Question 5
def psi_prime(t, psi, V):
    return -1j * (- 0.5 * hbar / m * Laplace(psi) + V / hbar * psi)
psi_prime(0.0, psi, V);

#%timeit psi_prime(0.0, psi)

# On vérifie que l'on a bien un état propre

a_ho = np.sqrt(hbar/(m*omega))
psi_0 = gaussian_packet(a_ho, 0, 0)

plt.plot(x, np.abs(-1J*omega/2*psi_0)**2)
plt.plot(x, np.abs(psi_prime(0, psi_0, V))**2)
plt.xlabel('x')
plt.ylabel(r'$|\psi|^2$')

Text(0, 0.5, '$|\\psi|^2$')

# Question 6

t0 = 0.0
tf = 1.5
t_eval = np.linspace(t0, tf, 101)

psi_0 = gaussian_packet(.3, -2, 0)

sol = integrate.solve_ivp(psi_prime, t_span = [t0, tf], y0 = psi_0, t_eval = t_eval, method="RK23", args=(V,))

plt.plot(x, np.abs(psi_0)**2)
plt.plot(x, np.abs(sol.y[:, ::5])**2);

# Question 7
def moments_x(psi):
    proba = np.abs(psi)**2
    norm = np.sum(proba)
    center = np.sum(proba*x)/norm
    var = np.sum(proba*(x-center)**2)/norm
    return center, var

centers, width = np.array([moments_x(psi) for psi in sol.y.T]).T

# Question 8
fig = plt.figure()

ax1, ax2 = fig.subplots(2, 1, sharex=True)

ax1.plot(t_eval, centers)
ax2.plot(t_eval, width)

ax2.set_ylabel('Width')
ax1.set_ylabel('Position')
ax2.set_xlabel('t')

Text(0.5, 0, 't')

# Question 9

p = np.fft.fftshift(hbar*2*np.pi*np.fft.fftfreq(len(x), d=dx))

def psi_tilde(psi):
    out = np.fft.fftshift(np.fft.fft(psi))
    return out


def moments_p(psi):
    mask = np.abs(p)<40
    proba = np.abs(psi_tilde(psi))**2
    norm = np.sum(proba[mask])
    center = np.sum(proba[mask]*p[mask])/norm
    var = np.sum(proba[mask]*(p[mask]-center)**2)/norm
    return center, var

momentum, temperature = np.array([moments_p(psi) for psi in sol.y.T]).T

fig = plt.figure()

ax1, ax2 = fig.subplots(2, 1, sharex=True)

ax1.plot(t_eval, momentum)
ax2.plot(t_eval, temperature)

ax1.set_ylabel('momentum')
ax2.set_ylabel('temperature')
ax2.set_xlabel('t')

Text(0.5, 0, 't')

# Question 10

t0 = 0.0
tf = .5
t_eval = np.linspace(t0, tf, 101)

a_ho = np.sqrt(hbar/(m*omega))
psi_0 = gaussian_packet(a_ho, 0, 0)

sol = integrate.solve_ivp(psi_prime, t_span = [t0, tf], y0 = psi_0, 
                          t_eval = t_eval, method="RK23", args=(0*V,))

psi_1 = sol.y[:,-1]

sol2 = integrate.solve_ivp(psi_prime, t_span = [t0, tf], y0 = psi_1, 
                          t_eval = t_eval, method="RK23", args=(V,))

momentum, temperature = np.array([moments_p(psi) for psi in sol2.y.T]).T

fig = plt.figure()

ax = fig.subplots(1, 1, sharex=True)

ax.plot(t_eval, np.sqrt(temperature))

ax.set_ylabel(r'$\Delta p$')
ax.set_xlabel('t')
ax.set_ylim(0, None)
ax.set_xlim(0, None)
ax.grid()