Movimiento de una partícula cargada en el campo de un dipolo magnético

Computa numéricamente las órbitas de un electrón en un campo magnético externo de forma dipolar.

El Lagrangiano de la partícula libre en presencia de un campo con simetría asimutal es:

\begin{equation} \mathcal{L} = -mc^2\sqrt{1-v^2/c^2} + \frac{e}{c}A(r,\theta)r\sin{\theta} v^{\phi}, \end{equation}

donde $A$ es la componente del co-vector campo de gauge. Recordemos que en estas coordenadas tenemos $v^2 = \dot{r}^2 + r^2 \dot{\theta}^2 + r^2\sin^2(\theta) \dot{\phi}^2$. Es conveniente redefinir $A \to A/c$ $m=1$, $e \to e/m$, ya que las cantidades solo dependen de estos parámetros.

La independencia del lagrangiano con respecto al tiempo implica que la energía es conservada, en este caso (sin potencial electrostático), viene dada por:

\begin{equation} H = \frac{c^2}{ \sqrt{1-v^2/c^2}} = \gamma c^2 , \end{equation} lo que implica que también el módulo de la velocidad es conservado. La invariancia asimutal nos dice que el momento angular asociado es también conservado,

\begin{equation} P := \frac{\partial \mathcal{L}}{\partial \dot{\phi}} = \frac{\dot{\phi} r^2 \sin^2{\theta} }{\sqrt{1-v^2/c^2} }+ e A(r,\theta) r \sin{\theta} = \gamma \dot{\phi} r^2 \sin^2{\theta} + e A(r,\theta) r \sin{\theta} \end{equation}

En coordenadas esféricas el sistema completo es:

\begin{eqnarray*} \frac{\partial \phi}{dt} &=& \frac{h}{\gamma r\sin{\theta}} \;\;\;\;\;\;\; h = [\frac{P}{r\sin{\theta}} - e A] \\ \frac{\partial \dot{r}}{dt} &=& e \frac{1}{ \gamma}\frac{\partial(rA)}{\partial r} \sin{\theta} \dot{\phi} + \frac{v^2 - \dot{r}^2}{r} \\ \frac{\partial \dot{\theta}}{dt} &=& [\frac{e}{\gamma} \frac{1}{r}\frac{\partial(\sin{\theta}A)}{\partial\theta} + \sin{\theta}\cos{\theta} \dot{\phi}]\dot{\phi} - 2\frac{\dot{r}}{r}\dot{\theta} \end{eqnarray*}

La conservación de $H$ también implica la conservación de $v^2$, usando la expresión para el momento angular, obtenemos,

\begin{equation} v^2 = \dot{r}^2 + r^2 \dot{\theta}^2 + r^2 \sin^2{\theta} \dot{\phi}^2 \;\;\;\;\; \dot{r}^2 + r^2 \dot{\theta}^2 = v^2 - \frac{h^2}{\gamma^2} . \end{equation}

Por lo tanto el movimiento es solo posible cuando $h^2 = [\frac{P}{r\sin{\theta}} - e A] < v^2\gamma^2$

Para un multipolo genérico con simetría axial tenemos,

\begin{equation} A(r,\theta) = \frac{a_n}{r^{n+1}}P^1_n(\cos{\theta}) \end{equation}

Caso dipolar

En particular tenemos, $P^1_1(\cos{\theta}) = -\sin(\theta)$.

Definiendo $e a_1 = a_0$,

\begin{eqnarray*} \frac{\partial \phi}{dt} &=& \frac{h}{\gamma r\sin{\theta}} \;\;\;\;\;\;\; h = [\frac{P}{r\sin{\theta}} + \frac{a_0}{r^2}\sin{\theta} ] \\ \frac{\partial \dot{r}}{dt} &=& \frac{1}{\gamma}\frac{a_0}{r^2} \sin^2{\theta} \dot{\phi} + \frac{v^2 - \dot{r}^2}{r} \\ \frac{\partial \dot{\theta}}{dt} &=& [- \frac{2a_0}{\gamma r^3} + \dot{\phi}]\dot{\phi}\sin{\theta}\cos{\theta} - 2\frac{\dot{r}}{r}\dot{\theta} \end{eqnarray*}

Donde \begin{equation} v^2 = \frac{\dot{r}^2 + r^2 \dot{\theta}^2 + h^2}{1+h^2}, \end{equation} es constante a lo largo del movimiento.

In [1]:
import numpy as np
from scipy.integrate import ode
#%matplotlib
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
#import matplotlib.animation as animation
from matplotlib import animation, rc
from IPython.display import HTML
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D
Vendor:  Continuum Analytics, Inc.
Package: mkl
Message: trial mode expires in 30 days
In [2]:
p=-0.1 # momento angular (reescaleado con m)
a = 0.1 # amplitud del campo (reescaleado con e, m y c)
In [3]:
N=5 # número de variables
In [4]:
# Primero las condiciones iniciales U=(rho,rho_t,theta, theta_t,phi)
x0 = 0. # tiempo inicial
rho0 = 1. 
#rho_t0 = -0.02
rho_t0 = -0.0
#pert = 0.000000000000000
pert = 0.01
theta0 = np.pi/2.0
theta_t0 = pert
phi0 = 0

h0 = (p/rho0/np.sin(theta0) + a*np.sin(theta0)/rho0/rho0)
h20 = h0*h0
v2 = (rho_t0*rho_t0 + rho0*rho0*(theta_t0*theta_t0) + h20) / (1.+h20)

if (v2 > 0.99):
                                print('velocidad demasiado alta! v2=%f', v2)
                                exit()
gamma = np.sqrt(1.-v2)                        

u0 = [rho0, rho_t0, theta0, theta_t0, phi0]

Definimos la funcion del integrador, o sea el sistema de ecuaciones:

In [5]:
def f(x, u):
    sint = np.sin(u[2])
    cost = np.cos(u[2])
    phi_t = (p/u[0]/sint + a/u[0]/u[0]*sint)/u[0]/sint/gamma
    
    return [u[1], a/u[0]/u[0]*sint*sint*phi_t/gamma + (v2 - u[1]*u[1])/u[0],u[3],
            (-2.*a/u[0]/u[0]/u[0]/gamma + phi_t)*phi_t*sint*cost - 2.*u[1]*u[3]/u[0],phi_t]

Definimos cual integrador usar. Damos distintos tipos, para este problema que en algunos casos es muy sensible, es recomendable usar la dop853, pero pruebe también la dopri5

In [6]:
#r = ode(f).set_integrator('zvode', method='bdf')
#r = ode(f).set_integrator('dopri5')  # Runge Kutta de orden 4/5
r = ode(f).set_integrator('dop853')   # Runge Kutta de orden 8
r.rtol = 1.e-6

Damos los valores iniciales

In [7]:
 r.set_initial_value(u0, x0)
Out[7]:
<scipy.integrate._ode.ode at 0x10b5e8390>

El valor final de integracion y dividimos el intervalo en trozos iguales a los fines de graficar puntos intermedios.

In [8]:
#x1 = 10000 # 
x1 = 5000.
# Esto es solo a los fines de graficar, si da un error: out of bounds poner K mas grande
dx = 0.1
K=50001
K_anim = 5000 # puede ser igual a K o KK, pero en el notebook tarda mucho en graficar.
#K=100000
#dx=1.
#K=10000
KK=0 # para tener solo los puntos justos en las  gráficas.
U = np.zeros((K,N))
X = np.zeros(K)
In [9]:
# Integramos K veces separadas (para graficar)
i = 0
while r.successful() and r.t < x1 and r.y[0] < rho0*12.:
    r.integrate(r.t+dx)
    U[i,:] = r.y
    X[i] = r.t
    i=i+1
    KK=KK+1
In [10]:
# Finalmente graficamos la solucion    
plt.ioff()
ax = plt.subplot(111, projection='3d')
ax.grid(True)
ax.set_title("general orbits", va='bottom')

# cambiar para ajustar mejor el gráfico

ax.set_xlim((-1.5, 1.5))
ax.set_ylim((-1.5, 1.5))
ax.set_zlim((-0.5, 0.5))

ax.plot(U[0:KK,0]*np.sin(U[0:KK,4])*np.sin(U[0:KK,2]),
        U[0:KK,0]*np.cos(U[0:KK,4])*np.sin(U[0:KK,2]),
        U[0:KK,0]*np.cos(U[0:KK,2]), color='r', linewidth=1)





plt.show()

Ahora hacemos una animación del gráfico anterior. Suele demorar, paciencia!

In [11]:
fig = plt.figure()
#ax1 = fig.add_subplot(111, projection='3d')
ax1 = fig.gca(projection='3d')
ax1.grid(True)
line, = ax1.plot([], [], [], color='r' ,lw=1)
time_template = 'time = %.1fs'
time_text1 = ax1.text(0.01, 0.01, 0.95, '', transform=ax1.transAxes)
ax1.set_xlim((-1.5, 1.5))
ax1.set_ylim((-1.5, 1.5))
ax1.set_zlim((-1.5, 1.5))

def init():
    line.set_data([], [])
    line.set_3d_properties([])
    time_text1.set_text('')
    return line, time_text1

def animate(i):
    line.set_data(U[0:i,0]*np.sin(U[0:i,4])*np.sin(U[0:i,2]),
        U[0:i,0]*np.cos(U[0:i,4])*np.sin(U[0:i,2]))
    line.set_3d_properties(U[0:i,0]*np.cos(U[0:i,2]))
    time_text1.set_text(time_template % (i*dx))
    return line, time_text1

ani = animation.FuncAnimation(fig, animate, K_anim,
                              interval=15, blit=True, init_func=init, repeat=False)

#rc('animation', html='html5')
HTML(ani.to_html5_video())
#ani.save('3d-dipole.mp4', fps=15)
Out[11]: