Valencia_Gomez_Jorge

Vista previa del material en texto

TRABAJO DE FIN DE GRADO
PARTÍCULAS Y CAMPOS EN UN
CONDENSADO DE BOSE-EINSTEIN
Jorge Valencia Gómez
Grado de Física
Facultad de Ciencias
Año académico 2021-22
PARTÍCULAS Y CAMPOS EN UN
CONDENSADO DE BOSE-EINSTEIN
Jorge Valencia Gómez
Trabajo de Fin de Grado
Facultad de Ciencias
Universidad de las Illes Balears
Año académico 2021-22
Palabras clave del trabajo:
condensado de Bose-Einstein, Ecuación Gross-Pitaevskii, impureza, hidrodinámica, vórtice
Nombre Tutor/Tutora del Trabajo: Cristóbal López Sánchez
Se autoriza la Universidad a incluir este trabajo en el Repositorio
Institucional para su consulta en acceso abierto y difusión en línea,
con fines exclusivamente académicos y de investigación
Autor/a Tutor/a
Sí No Sí No
� � � �
Resumen
En este trabajo, en primer lugar, estudiamos el movimiento de una impureza o partícula de
prueba inmersa en un condensado de Bose-Einstein (BEC) a temperatura cero. Se utiliza un
potencial repulsivo Gaussiano para modelar la interacción entre la impureza y las partículas
que forman el BEC, las cuales están descritas por una única función de onda macroscópica que
satisface la ecuación de Gross-Pitaevskii (GPE). La fuerza que ejerce el condensado sobre la
impureza se calcula mediante el teorema de Ehrenfest. Además, la llamada transformación de
Madelung permite escribir GPE en su versión hidrodinámica y estudiar el BEC como un fluido
cuántico. De esta manera, se comparan la fuerza ejercida, por parte del BEC, sobre la impureza
con la fuerza que actuaría sobre una partícula en el caso de la dinámica de fluidos clásica. Para
esta parte del trabajo se analizan dos casos simples: uno despreciando la presión cuántica y otro
considerándola.
La segunda parte del trabajo consiste en estudiar el flujo que se crea en el condensado
mediante un haz láser a temperatura diferente de cero. Para ello se resuelve numéricamente la
ecuación de Gross-Pitaevskii en dos casos: con potencial externo nulo y con potencial armónico.
Para el primero de los casos, se adopta una visión lagrangiana del fluido y se integra la posición
de una muestra de partículas del condensado obteniendo así sus trayectorias. Finalmente, a
partir de estas trayectorias se discute la existencia de un régimen caótico mediante el cálculo del
exponente de Lyapunov.
Resum
En aquest treball, primer estudiem el moviment d’una impuresa o partícula de prova immersa
en un condensat de Bose-Einstein (BEC) a temperatura zero. S’utilitza un potencial repulsiu
Gaussià per modelar la interacció entre la impuresa i les partícules que formen el BEC, les quals
estan descrites per una única funció d’ona macroscòpica que satisfà l’equació de Gross-Pitaevskii
(GPE). La força que exerceix el condensat sobre la impuresa es calcula mitjançant el teorema
d’Ehrenfest. A més, l’anomenada transformació de Madelung permet escriure GPE en la seva
versió hidrodinàmica i estudiar el BEC com a un fluid quàntic. D’aquesta manera, es comparen
la força exercida, per part del BEC, sobre la impuresa amb la força que hi actuaria sobre una
partícula en el cas de la dinàmica de fluids clàssica. Per a aquesta part del treball s’analitzen dos
casos simples: un menyspreant la pressió quàntica i un altre considerant-la.
La segona part del treball consisteix a estudiar el flux que es crea al condensat mitjançant
un feix làser a temperatura diferent de zero. Per això es resol numèricament l’equació de Gross-
Pitaevskii en dos casos: amb potencial extern nul i amb potencial harmònic. Per al primer dels
casos, s’adopta una visió lagrangiana del fluid i s’hi integra la posició d’una mostra de partícules
del condensat obtenint així les seves trajectòries. Finalment, a partir d’aquestes trajectòries es
discuteix l’existència d’un règim caòtic mitjançant el càlcul de l’exponent de Lyapunov.
Abstract
In this work, firstly, we study the motion of an impurity or test particle immersed in a Bose-
Einstein condensate (BEC) at zero temperature. A Gaussian repulsive potential is used to model
the interaction between the impurity and the particles that form the BEC, which are described
by a single macroscopic wave function that satisfies the Gross-Pitaevskii equation (GPE). The
force exerted by the condensate on the impurity is calculated using Ehrenfest’s theorem. In
addition, the so-called Madelung transformation allows writing GPE in its hydrodynamic version
and studying the BEC as a quantum fluid. In this way, the force exerted by the BEC on the
impurity is compared with the force that would act on a particle in the case of classical fluid
3
dynamics. For this part of the work, two simple cases are analyzed: one neglecting quantum
pressure and the other considering it.
The second part of the work consists of studying the flow that is created in the condensate
by means of a laser beam at a non-zero temperature. To do this, the Gross-Pitaevskii equation is
solved numerically in two cases: without external potential and with an harmonic-profile external
potential. For the first case, a Lagrangian view of the fluid is adopted and the position of a sample
of condensate particles is integrated, thus their trajectories are obtained. Finally, based on these
trajectories, the existence of a chaotic regime is discussed by calculating the Lyapunov exponent.
4
Índice general
1 Introducción 6
1.1 Descripción lagrangiana y euleriana de un fluido . . . . . . . . . . . . . . . . . . 7
1.2 Esquema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2 Fundamento teórico: descripción hidrodinámica 9
2.1 Ecuación de Gross-Pitaevskii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Hidrodinámica de un BEC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3 Fuerza sobre una impureza 12
3.1 Teorema de Ehrenfest . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.2 Potencial gravitatorio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.2.1 Despreciando la presión cuántica . . . . . . . . . . . . . . . . . . . . . . . 13
3.2.2 Considerando la presión cuántica . . . . . . . . . . . . . . . . . . . . . . . 14
4 Flujos complejos 17
4.1 Ecuación disipativa de Gross-Pitaevskii . . . . . . . . . . . . . . . . . . . . . . . 17
4.2 Vórtices cuánticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.2.1 Densidad de partículas en un vórtice . . . . . . . . . . . . . . . . . . . . . 19
4.2.2 Fase de la función de onda en un vórtice . . . . . . . . . . . . . . . . . . . 20
4.3 Potencial armónico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.3.1 Relajación de la función de onda . . . . . . . . . . . . . . . . . . . . . . . 21
4.3.2 Aproximación de Thomas-Fermi . . . . . . . . . . . . . . . . . . . . . . . 22
4.4 Sin potencial externo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.4.1 Relajación de la función de onda . . . . . . . . . . . . . . . . . . . . . . . 23
4.4.2 Regímenes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.4.3 Exponente de Lyapunov . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
5 Conclusiones 30
A Algoritmo numérico 32
A.1 Simulación dGPE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
A.2 Obtención de trayectorias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
B Deducción de la relación entre fase y velocidad 34
C Adimensionalización dGPE 35
Bibliografía 37
5
Capı́tulo 1
Introducción
Cuando los átomos de un gas se enfrían, su energía disminuye y también su velocidad, la cual
queda determinada por un rango de valores cada vez más acotado a medida que se rebaja la
temperatura. De esta forma, se conoce con menor precisión la posición de las partículas, tal y
como establece el principio de incertidumbre de Heisenberg. Si la temperatura sigue descendiendo,
Figura 1.1: Primera detección del con-
densado de Bose-Einstein. Se muestra
la distribución de momentos de los áto-mos de 87Rb. A medida que disminuye
la temperatura (de izquierda a dere-
cha) las partículas tienden a reducir
su energía y velocidad (tonos azules),
colapsando al mismo estado cuántico.
Imagen extraída de [1].
también lo hace el módulo del momento lineal, p, de
los átomos por lo que, según la relación de DeBroglie
λdB = h/p, donde h es la constante de Planck, el au-
mento en la longitud de onda es tal que las partículas
comienzan a convertirse en no localizadas. Es decir, la
naturaleza cuántica de los átomos predomina sobre la
naturaleza corpuscular. Por debajo de una temperatura
crítica, del orden de 500 - 800 nK, los átomos del gas
colapsan en el mismo estado cuántico (generalmente el
estado fundamental), donde todos ellos tienen la misma
energía, y forman un nuevo estado de la materia llamado
condensado de Bose-Einstein (BEC, por sus siglas en
inglés)[1]. Este nuevo estado solo puede estar formado
por bosones (partículas con espín entero en unidades de
~) debido al principio de exclusión de Pauli, que evita
que dos o más fermiones (cuyo espín es semientero en
unidades de ~) permanezcan en el mismo estado cuántico
simultáneamente. La primera predicción teórica de este
fenómeno fue debida a Einstein en 1924, continuando
con los estudios de Bose sobre la estadística de fotones
[2], y la primera observación experimental se produjo en
1995 mediante átomos de 87Rb [3], ver Fig. 1.1.
En las últimas décadas se han conseguido producir condensados en sistemas muy variados
y con innumerables aplicaciones en campos como superconductividad, superfluidez, sistemas
no lineales y sensores, entre otros [1]. Además, en la actualidad están presentes en una de las
principales y más novedosas ramas de investigación, la computación cuántica. Trabajos que datan
del 2004 [4] estudian esquemas de computación cuántica en los cuales se pueden implementar todo
un conjunto de puertas cuánticas (el análogo a las puertas lógicas en el mundo cuántico) mediante
dos condensados de Bose-Einstein formados por átomos distintos. Gracias a la superposición,
los estados fundamentales de dichos condensados pueden construir un estado entrelazado que
se comporta de manera muy similar a los qubits, lo que permite establecer una relación directa
entre los qubits y los BEC-qubits [5]. Asimismo, los condensados resultan ser de gran utilidad en
la creación de estados entrelazados porque estos minimizan la decoherencia.
6
1.1. Descripción lagrangiana y euleriana de un fluido
1.1 Descripción lagrangiana y euleriana de un fluido
A una temperatura igual a cero, y en un condensado en equilibrio, todas las partículas se
encuentran en el nivel de menor energía, con momento lineal nulo, formando el estado fundamental
del BEC. A medida que aumenta la temperatura o se abandona el equilibrio mediante fuerzas
externas, los bosones comienzan a poblar los niveles energéticamente superiores, dando lugar
a estados excitados, los cuales también están ocupados por un gran número de partículas [6].
Un análisis exhaustivo de estos estados (ya sea fundamental o excitado) pasaría por conocer la
función de onda formada por todos los constituyentes del sistema. Sin embargo, considerando
un gas muy diluido y debido a que, en este nuevo estado de la materia, todas las partículas se
encuentran en el mismo estado cuántico a causa de la deslocalización infinita, resulta razonable
aproximar el sistema por una gran onda de materia y describirlo con una única función de
onda efectiva [1]. Las partículas forman un estado colectivo cuya función de onda macroscópica
está modelada por la llamada ecuación de Gross-Pitaevskii (GPE, por sus siglas en inglés), la
cual puede escribirse en su versión hidrodinámica en términos de la densidad de partículas y la
velocidad del condensado. De esta forma, se puede tratar el BEC como un fluido cuántico.
Existen dos formas de caracterizar el movimiento de un fluido. Una de ellas es la descripción
Euleriana, que está ligada con la evolución temporal de las propiedades del fluido en un punto
espacial o en una región fija del dominio por donde circula [7]. Esta descripción resulta de gran
utilidad para analizar la distribución de densidad y velocidad en función de la posición y el
tiempo. El hecho de estudiar una región espacial y no las partículas del fluido conlleva una cierta
complicación matemática a la hora de calcular derivadas, ya que en regiones fijas el número
de partículas no es constante. Consecuentemente, se debe definir la derivada temporal de una
propiedad F (r, t) como la suma de un término no estacionario ∂F/∂t y un término advectivo
v ·∇F que tiene en cuenta el flujo de partículas que atraviesan la frontera del volumen de interés.
La otra forma de caracterizar el fluido es mediante la descripción lagrangiana, que está
relacionada con el estudio de las partículas individuales y es equivalente a la cinemática de una
partícula [7]. Estas partículas se etiquetan por sus posiciones iniciales ~r0 en un tiempo inicial
t0 y se siguen mientras se van moviendo por el fluido. Bajo este formalismo, la posición de la
partícula “i”, etiquetada por r0i y t0i , se escribe como ri(t; r0i , t0i) y cualquier propiedad, F , de
esta partícula será una función de la trayectoria seguida, es decir Fi = Fi[ri(t; r0i , t0i), t]. Para
los condensados que se estudian en este proyecto, la velocidad de las partículas del fluido es la
magnitud de mayor interés. La descripción lagrangiana, para este trabajo, será describir con
trayectorias las partículas de fluido del BEC, donde la velocidad de la partícula “i” viene dada
por vi = dri(t; r0i , t0i)/dt.
Tal y como las imágenes de Schrödinger y Heisenberg en la mecánica cuántica describen el
mismo sistema físico, las descripciones lagrangiana y euleriana, dada una magnitud del fluido F ,
deben representar la misma evolución temporal de esta.
1.2 Esquema
El objetivo del trabajo es, en primer lugar, estudiar las fuerzas hidrodinámicas que actúan
sobre una impureza sumergida en el condensado de Bose-Einstein y compararlas con las que
dicta la dinámica de fluidos clásica. En segundo lugar, es analizar la trayectoria seguida por las
partículas del fluido en flujos turbulentos y determinar si existe un régimen caótico. La memoria
se organiza de la siguiente manera:
Capítulo 2: se presenta la ecuación de Gross-Pitaevskii para un condensado de Bose-Einstein
ideal y su versión hidrodinámica.
Capítulo 3: A partir de las ecuaciones hidrodinámicas del condensado se estudian analíticamente
las fuerzas que actúan sobre la impureza considerando flujos sencillos, sin turbulencia. Los
7
1.2. Esquema
resultados se comparan con los de dinámica de fluidos clásica, discutidos por Maxey &
Riley en [8].
Capítulo 4: se introduce la ecuación amortiguada de Gross-Pitaevskii para condensados de
Bose-Einstein no ideales. Por otro lado, mediante métodos lagrangianos, se analiza el
movimiento de una partícula de fluido (no una impureza, sino una propia partícula de
fluido) en un flujo turbulento, el cual es creado por la rotación de un haz láser. A diferencia
del Capítulo 3, la complejidad de las ecuaciones conduce a la necesidad de un tratamiento
numérico. Además, se discute si la dinámica de las partículas es o no caótica mediante el
llamado exponente de Lyapunov.
Capítulo 5: se resumen y discuten los resultados obtenidos, así como futuros trabajos sobre
este tema.
Apéndice A: se describe el algoritmo numérico utilizado para la resolución de la ecuación de
Gross-Pitaevskii. También, se explica brevemente la modificación que se ha realizado al
algoritmo de Runge-Kutta para integrar ecuaciones diferenciales con matrices de datos en
lugar de funciones analíticas.
Apéndice B: se deriva la relación entre la fase de la función de onda y la velocidad del
condensado.
Apéndice C: se adimensionaliza la ecuación de Gross-Pitaevskii incluyendo el término repulsivo
propio de la interacción BEC-impureza.
8
Capı́tulo 2
Fundamento teórico: descripción
hidrodinámica
La fracción del número de partículas en el estado fundamentalde un condensado de Bose-
Einstein en función de la temperatura viene dado por [9]
N0
N
= 1−
(
T
Tc
)α
, Tc =
1
k
(
N
CαΓ(α)ζ(α)
)1/α
, (2.1)
con N0 el número de partículas en el estado fundamental, N el número de partículas total, T la
temperatura del sistema, Tc la temperatura crítica a partir de la cual las partículas comienzan
a poblar macroscópicamente el estado fundamental y α un parámetro real que depende de la
geometría del volumen confinante. En la segunda igualdad, k es la constante de Boltzmann, Cα
es una constante real que depende del volumen y, por su parte, Γ y ζ son las funciones Gamma
y Zeta de Riemann, respectivamente. De la Ec. (2.1) se desprende que para temperatura cero
todas las partículas del sistema se encontrarán en el estado fundamental formando el condensado.
A este tipo de condensados de Bose-Einstein se les llama ideales y la función de onda que los
describe satisface la ecuación de Gross-Pitaevskii. En este Capítulo introduciremos esta ecuación
y sus versiones hidrodinámicas, útiles para estudiar el condensado desde el punto de vista de la
dinámica de fluidos.
2.1 Ecuación de Gross-Pitaevskii
A temperatura igual a cero, un condensado de Bose-Einstein está bien descrito por una
función de onda ψ(r, t) tal que verifica la ecuación de Gross-Pitaevskii [10]:
i~
∂ψ
∂t
=
(
− ~
2
2m∇
2 + gB|ψ|2 − µ+ V
)
ψ. (2.2)
Esta tiene la forma de una ecuación de Schrödinger con un potencial efectivo que contiene un
término no lineal, gB|ψ|2, que modela el campo medio producido por los bosones del condensado
[9]. Por su parte, ψ está normalizada al número de partículas,
∫
|ψ|2d3r = N . En la Ec. (2.2), µ
es el potencial químico; m, la masa de las partículas del condensado; V (r, t), el potencial externo
y el parámetro gB está relacionado con la interacción entre los átomos que forman el BEC y tiene
unidades de energía multiplicada por volumen. Para un gas diluido las interacciones entre más
de dos cuerpos son despreciables frente a la interacción a pares, la cual viene dada por fuerzas
de Van der Waals [1]. La interacción descrita por gB puede ser repulsiva o atractiva según si su
valor es positivo o negativo, respectivamente. Experimentalmente se puede controlar tanto su
9
2.2. Hidrodinámica de un BEC
signo como su valor mediante campos magnéticos [1]. Para este trabajo se estudia el caso gB > 0
donde la interacción a pares es repulsiva. En este régimen gB viene dado por
gB =
4π~2as
m
, (2.3)
siendo as la longitud de dispersión en onda s 1. Para un condensado con interacciones débiles,
as � λdB, siendo λdB la longitud de onda de DeBroglie.
El sistema descrito por la ecuación de Gross-Pitaevskii puede modificarse introduciendo
una impureza que se mueve por el BEC. Para modelar la interacción entre la impureza y
las partículas del condensado se utiliza un potencial Gaussiano repulsivo en 3 dimensiones
Up(r, rp, t) = µ/(
√
2πσ)3e−(r−rp)2/(2σ2), donde σ es el tamaño efectivo de la impureza y rp, la
posición de su centro de masas. Así, la Ec. (2.2) toma la forma
i~
∂ψ
∂t
=
(
− ~
2
2m∇
2 + gB|ψ|2 − µ+ V + gpUp
)
ψ. (2.4)
donde gp es una constante (con unidades de volumen) que, para este trabajo, es positiva, ya que
caracteriza la repulsión entre la impureza y las partículas del BEC.
Para simplificar los futuros cálculos y siguiendo con el procedimiento de [10], se puede
adimensionalizar GPE introduciendo las variables ψ̃ = ψ/ψ0, t̃ = t/t0 y r̃ = r/ξ. Donde
ψ0 =
√
µ/gB es la solución homogénea de GPE y t0 = ~/µ, la unidad característica de tiempo.
Por su parte, ξ = ~/√mµ es la longitud de coherencia (healing length, en inglés), que es la mínima
distancia sobre la cual la función de onda puede variar. La deducción de las expresiones para ψ0,
t0 y ξ se encuentra en el Apéndice C. Definiendo
Ṽ ≡ V (r̃, t̃)/µ, g̃pŨp ≡ gpUp(r̃, r̃p, t̃)/µ
con g̃p = gp/ξ3, Up(r̃, r̃p, t̃) = µ/(
√
2πσ)3e−(r̃−r̃p)2/(2a2) y siendo a = σ/ξ, la Ec. (2.4) se puede
escribir como:
∂ψ̃
∂t̃
= i
(1
2∇̃
2 − |ψ̃|2 + 1− Ṽ − g̃pŨp
)
ψ̃. (2.5)
Para lo que resta de trabajo no escribiremos “ ∼ ” sobre las variables por tal de simplificar la
notación. Se sobreentiende que son variables adimensionales, a no ser que se indique lo contrario.
2.2 Hidrodinámica de un BEC
Con la finalidad de entender el condensado desde el punto de vista hidrodinámico resulta útil
derivar la ecuación de continuidad. Para ello, multiplicando por ψ∗(r, t) la Ec. (2.5) y restándole
el complejo conjugado se llega a
∂|ψ|2
∂t
+ ∇ ·
[ 1
2i (ψ
∗∇ψ − ψ∇ψ∗)
]
= 0, (2.6)
donde se ha considerado que los potenciales V y gpUp de GPE son reales. Definiendo |ψ|2 como
la densidad de partículas del condensado, n, y comparando con la ecuación de continuidad propia
de la mecánica de fluidos clásica
∂ρ
∂t
+ ∇ · (ρv) = 0, (2.7)
siendo ρ(r, t) = mn(r, t) la densidad de masa, se observa que la Ec. (2.6) es una ecuación de
continuidad para n. Esta comparación conduce a la definición de la velocidad del condensado
como:
v = 12ni (ψ
∗∇ψ − ψ∇ψ∗) = 1
n
=(ψ∗∇ψ), (2.8)
1El valor característico de as para átomos de 87Rb es del orden de 5.8 nm [1]
10
2.2. Hidrodinámica de un BEC
donde =(ψ∗∇ψ) es la parte imaginaria de ψ∗∇ψ. Por tanto, la ecuación de continuidad para un
condensado de Bose-Einstein se escribe
∂n
∂t
+ ∇ · (nv) = 0, (2.9)
que refleja la conservación del número de partículas (consistente con la normalización de ψ de la
Sección 2.1).
La función de onda que describe el condensado es compleja y como tal puede expresarse en
términos de su módulo y fase:
ψ(r, t) =
√
n(r, t)eiS(r,t). (2.10)
De tal forma que se verifica
v(r, t) = ∇S(r, t), (2.11)
como se demuestra en el Apéndice B. De esta manera, la función de onda queda expresada en
términos hidrodinámicos, lo que se conoce como transformación de Madelung [1]. De la Ec. (2.11)
se concluye que cuanto mayor sea la variación espacial de la fase, mayor será la velocidad de las
partículas [11]. Además, el flujo de partículas del condensado es irrotacional, ya que la velocidad
es el gradiente de una función escalar. Es decir, la vorticidad, ω = ∇× v = 0 siempre y cuando
S no sea singular2 [9]. Por otro lado, la divergencia del campo de velocidades generalmente no es
nula, ∇ · v = ∇2S 6= 0. Por lo que el condensado, además de irrotacional, representa también un
fluido compresible.
La segunda ecuación hidrodinámica del condensado se obtiene al combinar GPE con la
transformación de Madelung, llegando a
∂v
∂t
= −∇
(
1
2v
2 + n− 12
∇2
√
n√
n
+ V + gpUp
)
. (2.12)
En la Ec. (2.12), el término −∇(v2/2) tiene su origen en la energía cinética de las partículas del
condensado, mientras que −∇n es la fuerza debida a la existencia de un gradiente de densidad en
el fluido (esta contribución se anula en el caso de un fluido homogéneo). La fracción ∇2
√
n/
√
n
se conoce como el término de presión cuántica y −∇(V + gpUp) son las fuerzas producidas sobre
el fluido originadas por el potencial externo y la interacción BEC-impureza, respectivamente.
Tras alguna manipulación algebraica la Ec. (2.12) puede reescribirse como [1]
n
(
∂v
∂t
+ (v ·∇)v
)
= −∇(P + PQ)− n∇(V + gpUp), (2.13)
que tiene la forma de la ecuación de Euler para fluidos no viscosos (consultar Capítulo 4 de [7])
con una presión efectiva formada por dos términos. El primero de ellos es la presión del fluido y
viene dada por P = n2/2. Puesto que la presión es únicamente función de la densidad, se dice
que el fluido es barotrópico. Por otro lado, el término PQ = n∇2(lnn)/4 representa la presión
cuántica3. Comparando los órdenes de magnitud entre ambas contribuciones de la presión se
concluye que PQ � P cuando n varía en una escala similar a la longitud de coherencia, es decir,
cerca de las fronteras del sistema; mientras que PQ � P en el resto de regiones espaciales [1, 9].
En resumen, un condensado de Bose-Einstein suficientemente diluido puede tratarse de forma
general como un fluido inhomogéneo, compresible, barotrópico, con un flujo irrotacional y que
satisface unas ecuaciones que tienen la misma formaque las de la dinámica de fluidos clásica
pero con modificaciones puramente cuánticas.
2En el Capítulo 4 veremos que esta última condición no se cumple en el centro de un vórtice cuántico, ya que
en tales puntos S no está definida y, por tanto, tampoco lo estará la vorticidad.
3Se puede demostrar matemáticamente que ∇2
√
n/(2
√
n) = n∇2(lnn)/4, por lo que ambos términos representan
la presión cuántica.
11
Capı́tulo 3
Fuerza sobre una impureza
En este Capítulo, consideraremos un BEC ideal confinado en un cilindro de altura L y cuyo
radio haremos tender a infinito en aras de simplificación. Posteriormente, introduciremos una
impureza en el condensado y estudiaremos la fuerza a la que se ve sometida desde un punto de
vista hidrodinámico.
3.1 Teorema de Ehrenfest
El teorema de Ehrenfest establece el análogo de la segunda ley de Newton en el ámbito
cuántico y viene dado por [11]
m
d2
dt2 〈x〉 =
d
dt 〈p〉 = −〈∇V (r)〉 , (3.1)
siendo 〈x〉, 〈p〉 y 〈∇V (r)〉 los valores esperados de la posición, momento lineal y gradiente del
potencial, respectivamente, en el estado del condensado. Esta fórmula resulta útil para calcular
la fuerza que sufre una impureza sumergida en un condensado de Bose-Einstein. Sustituyendo
V (r) por el potencial de interacción gpUp en la Ec. (3.1) se obtiene FB, la fuerza que ejerce
dicha impureza sobre el condensado. Puesto que nuestro interés reside en la fuerza que sufre la
partícula de prueba, Fp, podemos hacer uso de la tercera ley de Newton para escribir
Fp(t) = −FB(t) = gp
∫
d3r ∇Up(r, rp, t) n(r, t)
i.p.p= −gp
∫
d3rUp(r, rp, t)∇n(r, t), (3.2)
en unidades adimensionales [10]. Para la última igualdad se ha integrado por partes y se ha
tenido en cuenta que la densidad de partículas se anula en las fronteras del volumen confinante.
Pese a que la integral de la Ec. (3.2) se extiende para todo el volumen del condensado, cuando
|r| � |rp| la contribución tiende a cero debido al corto rango de interacción de Up. También se
nota que si la densidad de partículas es homogénea (∇n = 0) dentro de este rango de interacción
de Up, entonces la fuerza sobre la impureza es nula.
3.2 Potencial gravitatorio
En esta sección estudiamos las fuerzas que actúan sobre una partícula de prueba inmersa en
un potencial gravitatorio y localizada en un condensado de Bose-Einstein estacionario (|v| = 0) el
cual está confinado en un cilindro orientado en la dirección del eje z. Las dimensiones del sistema
vienen caracterizadas por la altura del cilindro, que supondremos L = 10ξ (L = 10, en unidades
de la longitud de correlación) y su radio, el cual hacemos tender a infinito R→∞. Por su parte,
el potencial V de la ecuación de Gross-Pitaevskii toma la forma de potencial gravitatorio mgz,
12
3.2. Potencial gravitatorio
siendo m la masa de las partículas del condensado y g la aceleración gravitatoria. Para analizar
este sistema primero se ignora el término de la presión cuántica en la Ec. (2.12) obteniendo así
un perfil de densidad aproximado y, posteriormente, se mejoran los resultados considerando la
totalidad de la ecuación de Euler.
3.2.1 Despreciando la presión cuántica
Antes de introducir una impureza en el condensado, estudiamos cómo se comportan las
partículas de éste bajo un potencial gravitatorio. Eliminando los términos ∇2
√
n/
√
n, gpUp y
teniendo en cuenta un BEC estacionario se resuelve la Ec. (2.12) obteniendo un perfil de densidad1
n(z) = C − Kz, (3.3)
que representa la distribución de partículas que forman el condensado en función de la altura.
En la Ec. (3.3) aparecen dos constantes adimensionales: K = mgξ/µ que proviene de la adimen-
sionalización del potencial externo y C, que surge de la integración de la ecuación de Euler y se
determina con la condición de normalización
∫
n(z)d3r = N . Para este Capítulo se escoge un N
tal que la densidad de partículas en el centro del cilindro es 1. Aplicando esta última condición
se obtiene
n(z) = 1 + 12KL−Kz, (3.4)
de tal forma que en ausencia de potencial gravitatorio,K = 0, se recupera la densidad de partículas
para el caso homogéneo n(z) = 1, que en variables dimensionales sería n(z) = n0 = |ψ0|2 = µ/gB.
Para describir un sistema físico n(z) debería anularse en la base del cilindro, z = 0, y no
diverger para grandes alturas. Sin embargo, el perfil de densidad dado en la Ec. (3.4) describe un
condensado ficticio, pues no cumple estas condiciones de contorno. Pese a ello, esta expresión
puede usarse como una primera aproximación al comportamiento real de n(z) en zonas alejadas
de las paredes y la base del recipiente, donde la presión cuántica es despreciable (tal y como se
comentó en la Sección 2.2). En esa región, la densidad de partículas decae linealmente con la
altura.
Una vez calculado n(z) se puede introducir una impureza en el sistema y calcular la fuerza
que las partículas del condensado ejercen sobre ella mediante el teorema de Ehrenfest. Para ello,
se calcula la integral de la Ec. (3.2) únicamente en las regiones donde n(z) es válida, es decir,
suficientemente lejos de las fronteras del volumen confinante. De esta forma, para un entorno de
la impureza podemos suponer que que la base y la tapa del cilindro se encuentran infinitamente
lejos, despreciando así los efectos de borde. Por tanto la Ec. (3.2) da lugar a
Fp(t) = gpK ẑ, (3.5)
en unidades adimensionales. Para lograr comparar con los resultados de la teoría clásica de fluidos
debemos restaurar las dimensiones. Teniendo en cuenta que gp = ξ3g̃p y que Fp(t) = (µ/ξ)F̃p(t):
Fp(t) = gpρg ẑ, (3.6)
donde ρ ≡ m/ξ3 es una constante con unidades de densidad que difiere de la densidad másica
del BEC ya que esta última depende de la posición, ρBEC = mn(z).
La expresión de la fuerza de la Ec. (3.6) puede compararse con la fuerza de flotabilidad de
dinámica clásica de fluidos [8], por la cual una partícula inmersa en un fluido experimenta una
fuerza en sentido opuesto a su peso y proporcional al volumen de fluido desalojado. Para el caso
de una impureza esférica de radio σ y totalmente sumergida, la fuerza de flotabilidad es
Fp(t) =
4
3πσ
3ρfg ẑ, (3.7)
1Para este cálculo primero se ha escrito el potencial gravitatorio V (r, t) = mgz en unidades adimensionales.
Para ello: Ṽ ≡ V (r̃, t̃)/µ = (mgξ/µ)z̃.
13
3.2. Potencial gravitatorio
siendo 4πσ3/3 el volumen de la impureza, ρf la densidad constante del fluido y g la gravedad.
Este resultado es equivalente al de la Ec. (3.6) donde gp hace el papel del volumen de la partícula
de prueba y ρ, pese a no representar la densidad del fluido cuántico, tiene unidades de densidad al
igual que ρf. Esta analogía es válida solo si las inhomogeneidades introducidas por la presencia de
la impureza en el fluido son despreciables y los cambios en la velocidad de este varían débilmente
[10].
Las pequeñas modificaciones que causa la impureza se pueden tratar como perturbaciones del
estado homogéneo del fluido [10]. Para ello, podemos realizar un análisis perturbativo a primer
orden sustituyendo
n(r, t) = nh + gpδn(r, t)
v(r, t) = gpδv(r, t)
(3.8)
en las Ecs. (2.9) y (2.12). En la Ec. (3.8), nh describe el perfil de densidad sin impureza, dado por
la Ec. (3.4). Por su parte, cuando no hay ninguna impureza en el condensado, este se encuentra
en reposo, |vh| = 0, mientras que al introducirla se genera un pequeña modificación δv. Se llega
así a una ecuación de onda para las inhomogeneidades en la densidad del fluido:
∂2
∂t2
δn = ∇ · [nh∇(δn+ Up)], (3.9)
para el caso con V = 0. La solución son ondas planas emitidas en la cercanía de la impureza
cuando esta se introduce en el sistema o se mueve por él. Dichas ondas cumplen con la relación de
dispersión de Bogoliubov ω(k) = k
√
1 + k2/4 (con k el módulo del vector de onda), a partir de la
cual se puede calcular la velocidad de fase, vφ = w/k, y la velocidad de grupo, vg = ∂ω/∂k [10].
En el límite k → 0 o longitud de onda larga, propio de las ondas de sonido, ambas velocidades
son idénticas y valen c = 1 (en unidades adimensionales2). Además, pese a estar fuera del alcance
del trabajo,se podría calcular la fuerza sobre la impureza con la Ec. (3.2), lo que daría lugar a
otras correcciones añadidas a la fuerza de flotabilidad, que podrían comparase con su análogo
clásico de la Ref. [8].
3.2.2 Considerando la presión cuántica
Cuando se tiene en cuenta el término de presión cuántica y considerando un condensado con
las mismas características que en la sección anterior, la Ec. (2.12) adopta la forma
∂2ζ
∂z2
= 2
[
(Kz − η)ζ + ζ3
]
. (3.10)
En la ecuación anterior se ha definido ζ2 ≡ n y η como una constante adimensional que asegure la
conservación del número de partículas en el sistema. Por tanto, η se calcula mediante la condición
de normalización, que puede escribirse como
∫
d3rζ2(z) = N .
Las soluciones para la Ec. (3.10) en un volumen cilíndrico de altura L = 10, para diferentes
valores de K y η, se muestran en el panel (A) de la Fig. 3.1. En el método numérico utilizado
para resolver la ecuación diferencial, se he impuesto que la función de onda debe anularse en
las paredes impenetrables del sistema. Consecuentemente, n es cero tanto en la base como en la
parte superior del cilindro, en contraposición con el comportamiento de los perfiles de densidad
dados por la Ec. (3.4). Por su parte, la simetría de este cilindro de radio infinito nos permite
simplificar el sistema a una sola dimensión (eje z). Por este motivo, el área bajo cada una de las
curvas en el panel (A) de la Fig. 3.1 representa el número de partículas por unidad de área, el
cual es constante. Se observa como a medida que aumenta la intensidad del potencial gravitatorio,
es decir, K incrementa, la densidad de partículas, n, crece cerca de la base del cilindro. Este
comportamiento es el esperado, pues un valor mayor de K = mgξ/µ conlleva un incremento en
2La velocidad se adimensionaliza como c = c0c̃ siendo c0 = ξ/t0. Usando las relaciones de la Ec. (C.4) se llega a
c0 =
√
µ/m, que es la velocidad de propagación para las ondas de sonido emitidas por la impureza.
14
3.2. Potencial gravitatorio
0
1n
A
Densidad de partículas
0.5
0.0
F p
1e 2
B
Fuerza sobre la impureza
Parámetros adimensionales
= 0, 1.00
= 0.01, 1.05
= 0.05, 1.25
= 0.1, 1.49
= 0.15, 1.72
= 0.2, 1.94
= 0.25, 2.15
= 0.3, 2.39
0 1 2 3 4 5 6 7 8 9 10
z
0.0
0.5
E
C
Energía
Figura 3.1: Representación, mediante variables adimensionales, de varias magnitudes relevantes
en función de la distancia de las partículas a la base del cilindro, z (en unidades de ξ). El
panel (A) muestra el resultado de la integración numérica de la Ec. (3.10), donde la densidad
está normalizada por n0 ∼ 1020 m−3. En el panel (B) se dibuja la fuerza experimentada por
la impureza sumergida en el condensado, la cual ha sido calculada mediante la Ec. (3.2) con
gp = 0.01 y a = 1. Las líneas discontinuas reflejan la fuerza para el caso en el cual se desprecia
la presión cuántica. Finalmente, el panel (C) representa los perfiles de energía, calculados con
F = −∇E. Para cada panel se utilizan un amplio rango de valores para los parámetros K y η. La
fuerza y la energía se normalizan por µ/ξ ∼ 7.35×10−25 N y µ ∼ 7.35×10−31 J, respectivamente,
donde ξ ∼ 1 µm.
la gravedad g de nuestro sistema. Asimismo, otra consecuencia del aumento de K es la reducción
de la longitud de coherencia, puesto que la distancia entre el origen, donde ψ = 0, y el máximo
de la función de onda va acortándose.
La presión cuántica es únicamente relevante para escalas espaciales del orden o inferiores a ξ,
por lo que lejos de la base y la parte superior del cilindro puede despreciarse. De esta forma, en
tales regiones el perfil de la densidad de partículas tiene un comportamiento lineal con z, como
se esperaba de la Ec. (3.4), ver el panel (A) de Fig. 3.1. Los valores del gradiente (la derivada
con respecto a z) de los perfiles de n en estos puntos son cercanos a los de la pendiente mostrada
en la Ec. (3.4). En la Fig. 3.2 se muestran los perfiles lineales de densidad cuando no se considera
la presión cuántica, que coinciden con los perfiles no lineales de densidad en zonas ubicadas lejos
de la base y la parte superior del cilindro para valores de K hasta 0.1 (ver panel (A) de la Fig.
3.2). Cuando K, es decir, la gravedad, aumenta: 0.15, 0.2, 0.25 y 0.3 (panel (B) de la Fig. 3.2), la
curvatura del perfil de densidad es cada vez mayor, alejándose del comportamiento lineal de la
aproximación. Entonces, el término de la presión cuántica (que depende de la segunda derivada
de la densidad) deja de ser despreciable en zonas cada vez más alejadas de las paredes.
Si una impureza se sumerge en el condensado a una cierta altura, la fuerza que sufrirá, que
podemos llamar fuerza de flotabilidad, se representa en el panel (B) de la Fig. 3.1. A una altura
determinada, diferente para cada K, se produce un cambio de signo de la fuerza que sitúa un
15
3.2. Potencial gravitatorio
0 1 2 3 4 5 6 7 8 9 10
z
0
1
2
n
A
0 1 2 3 4 5 6 7 8 9 10
z
B
Densidad de partículas
Parámetros adimensionales
= 0, 1.00
= 0.01, 1.05
= 0.05, 1.25
= 0.1, 1.49
= 0.15, 1.72
= 0.2, 1.94
= 0.25, 2.15
= 0.3, 2.39
Figura 3.2: El eje y se comparte entre los dos gráficos, los cuales representan la densidad de
partículas, n, en función de z (ambas adimensionales). A diferencia de la Fig. (3.1), se muestran
los perfiles de densidad para el caso sin presión cuántica (líneas discontinuas) para diferentes
valores de K y η.
punto crítico por debajo del cual la partícula se hundirá hasta llegar a la base del cilindro. Justo
en el punto crítico, la fuerza de flotación ejercida por el BEC sobre la partícula es cero. Por
encima del punto crítico, la partícula flotará siempre que la fuerza de flotación sea mayor que su
peso. Ese comportamiento se puede observar fácilmente en el panel (C) de la Fig. 3.1 donde se
muestra el perfil de energía, E, calculado a partir de F = −∇E. El punto crítico se encuentra en
el máximo de cada curva. La zona de hundimiento está a la izquierda de ese máximo y la zona
de flotabilidad (si el peso de la partícula lo permite) está a la derecha. Para pequeños K como 0
o 0.01, pequeñas perturbaciones alrededor de z ∼ 5 desencadenarán el hundimiento o la flotación
de la partícula.
Además, la fuerza ejercida sobre una partícula de prueba lejos de las paredes del cilindro
debe ser aproximadamente constante, como se muestra en el panel (B) de la Fig. 3.1 entre
z ∼ 3 y z ∼ 7. Numéricamente, esta fuerza coincide con la Ec. (3.5) (líneas discontinuas) para
k = 0, 0.01, 0.05, 0.1 y 0.15 y comienza a discrepar para K superiores, a partir de donde los efectos
de la presión cuántica empiezan a ser notables.
16
Capı́tulo 4
Flujos complejos
Cuando la temperatura del sistema es lo suficientemente grande, las partículas pueden adquirir
la energía necesaria como para abandonar el estado de condensado y ocupar otros estados. Estos
forman un baño térmico [6, 10] que se eleva sobre el gran número de partículas que continúan
en el estado de condensado y, por tanto, se describen mediante GPE. Por tanto, a diferencia
de los capítulos anteriores, ahora el sistema físico no solo está formado por el condensado
sino que también se debe tener en cuenta el baño térmico a su alrededor. Este hecho conlleva
nuevas interacciones condensado-baño que modifican ligeramente las ecuaciones del Capítulo 2
introduciendo términos disipativos. En este Capítulo trataremos condensados bidimensionales
cuyos flujos, causados por el movimiento de un haz láser, son no estacionarios, inhomogéneos
y localmente irrotacionales. En primer lugar, presentaremos las modificaciones introducidas al
considerar un BEC con baño térmico y se justificará matemáticamente la aparición de vórtices
en el flujo. Posteriormente, obtendremos los campos de velocidades dependientes del tiempo, los
cuales se integrarán para, mediante la descripción lagrangiana del fluido, calcular las trayectorias
de algunas partículas de este. Finalmente, se discutirá sobre la existencia de caos en ausencia de
potencialexterno.
4.1 Ecuación disipativa de Gross-Pitaevskii
Considerar un sistema a una temperatura apreciable da lugar a la aparición de una fuerza
disipativa debido a la interacción de las partículas que forman el condensado con las del baño
térmico. Esta fuerza es comparable con la fuerza de arrastre de la dinámica clásica de fluidos
[8, 10] y debe, de alguna manera, tenerse en cuenta a la hora de formular la ecuación que gobierna
el condensado. Bajo estas circunstancias y teniendo en cuenta una impureza sumergida en el
BEC, la función de onda efectiva satisface [13]:
i~
∂ψ
∂t
= (1− iγ)
(
− ~
2
2m∇
2 + gB|ψ|2 − µ+ V + gpUp
)
ψ, (4.1)
que es la llamada ecuación disipativa de Gross-Pitaevskii adimensional. Donde γ es el coeficiente
disipativo o de arrastre térmico, un parámetro definido positivo que cuantifica el intercambio de
partículas entre el condensado y el baño mediante procesos de colisión [10].
Por tal de producir flujos turbulentos en el condensado se utiliza un haz láser que se desplaza
siguiendo un movimiento circular uniforme. El haz se puede interpretar como una impureza que
va “removiendo” el condensado, por lo que la interacción entre ambos está descrita por gpUp.
Así, se puede escribir gpUp(x, y, t) = V0 exp
[
− (x−xp(t))
2+(y−yp(t))2
2σ2
]
[12], siendo V0 ≡ gpµ/(2πσ2)
donde, para el caso bidimensional, gp tiene unidades de área. La posición del haz del láser viene
especificada por xp(t) = s cos(v`t/s) e yp(t) = s sin(v`t/s) con v` el módulo de la velocidad lineal
17
4.2. Vórtices cuánticos
de este y s el radio de la órbita circular que describe. La posición inicial del haz es (s, 0), a partir
de donde comienza un movimiento circular de periodo P = 2πs/v`. Por otro lado, consideramos
un potencial externo armónico V (x, y, t) = mω2ρ(x2 + y2)/2 donde ω2ρ = ω2x + ω2y es la frecuencia
de confinamiento en la dirección eρ (en coordenadas cilíndricas), siendo ωx = ωy debido a la
isotropía del sistema bajo estudio.
Sustituyendo gpUp y V de los párrafos anteriores en la Ec. (4.1) y escribiendo esta última con
variables adimensionales se llega a,
∂ψ̃
∂t̃
= (i+ γ)
{
1
2∇̃
2 − |ψ̃|2 + 1− 12Ω
2
ρ(x̃2 + ỹ2)− V0 exp
[
−(x̃− x̃p(t̃))
2 + (ỹ − ỹp(t̃))2
2a2
]}
ψ̃.
(4.2)
cuya adimensionalización se encuentra en el Apéndice C. En la Ec. (4.2), Ωρ ≡ ωρt0 y V0 ≡ V0/µ
son parámetros adimensionales que caracterizan la fuerza de la interacción del potencial confinante
y del potencial de interacción entre el haz láser y las partículas del condensado, respectivamente.
Como es costumbre, no escribiremos “ ∼ ” sobre las variables adimensionales por tal de simplificar
la notación.
Análogamente al sistema con temperatura cero, para este caso con T 6= 0, la transformación
de Madelung da paso a [13]:
∂n
∂t
+ ∇ · (nv) = 2nγ(1− Ueff) (4.3)
y
∂v
∂t
= −∇
(
Ueff −
γ
2n∇ · (nv)
)
, (4.4)
donde se ha definido el potencial efectivo Ueff ≡ 12v
2 + n− 12
∇2
√
n√
n
+ V + gpUp. Las Ecs. (4.3) y
(4.4) son una clara generalización de la ecuación de continuidad y de la ecuación de Euler del
Capítulo 2, respectivamente.
Si se considera un condensado homogéneo y estático (|v| = 0) y asumiendo Ueff = n, la Ec.
(4.3) indica que n = C, ∀C ∈ R, es solución para sistemas sin baño térmico (γ = 0), mientras
que en el caso contrario (γ 6= 0) solo n = 1 es solución. Es decir, la disipación provocada por el
baño térmico implica la existencia de un valor estacionario de la densidad en n = 1. Sin embargo,
esta densidad de equilibrio puede tomar cualquier valor en condensados ideales.
4.2 Vórtices cuánticos
Dos conceptos clave para comprender los vórtices son la vorticidad y la circulación, identificados
con ω y Γ, respectivamente. El primero de ellos cuantifica la capacidad de los elementos de
volumen de un fluido para rotar alrededor de sí mismos y se define como el rotacional del campo
de velocidad: ω = ∇× v. Por su parte, la circulación se calcula como la integral de línea de v a
lo largo de un contorno cerrado C, es decir: Γ =
∮
C v · d`. Ambos conceptos guardan una estrecha
relación pues se puede demostrar utilizando el teorema de Stokes que la vorticidad no es más
que la circulación por unidad de área. A continuación, se harán uso de estas magnitudes para
describir el flujo de un condensado de Bose-Einstein.
Tal y como se desprende de la Ec. (2.11), v es un campo conservativo ya que se calcula como
el gradiente de un campo escalar, siendo este último la fase de la función de onda. Entonces, como
cualquier otro campo conservativo, su integral de línea será idénticamente cero, por lo que Γ = 0,
siempre y cuando S no sea singular [9]. Como la circulación y la vorticidad son proporcionales,
un valor nulo para la circulación se traduce en ω = 0 en cualquier punto dentro del contorno
considerado, es decir, el flujo es irrotacional en esa región espacial. Sin embargo, si en el área que
encierra el camino C existe algún punto donde S sea singular, la función de onda debe seguir
18
4.2. Vórtices cuánticos
siendo univaluada [9], por lo que:
Γ =
∮
C
v · d` =
∮
C
∇S · d` = ∆S, (4.5)
donde ∆S ≡ S(2) − S(1) = 2πl, con l ∈ Z. Es decir, la circulación es diferente de cero si el
contorno considerado encierra una singularidad de la fase [9], dando lugar a la existencia de
puntos dentro de C donde la vorticidad no es cero. Además, a la vista del resultado de la Ec.
(4.5) y a diferencia de los vórtices clásicos donde la circulación puede tomar cualquier valor, los
vórtices cuánticos tienen una circulación cuantizada en intervalos discretos de 2π.
4.2.1 Densidad de partículas en un vórtice
Para entender los vórtices en el régimen cuántico consideraremos el caso más sencillo: un
condensado estacionario (∂tψ = 0), sin potencial externo (Ωρ = 0) y sin haz láser (V0 = 0) en el
que existe un vórtice en la posición (x, y) = (0, 0). Bajo estas suposiciones, la Ec. (4.2) queda:
0 =
(
1 + 12∇
2 − |ψ|2
)
ψ, (4.6)
donde se ha tenido en cuenta que (1 +γ) 6= 0 ya que γ > 0. Sabiendo que la función de onda tiene
una parte real y otra imaginaria, podemos utilizar el Ansatz: ψ(ρ, ϕ) = f(ρ)eig(ϕ) en coordenadas
polares, siendo f y g funciones reales y f(ρ) ≥ 0 ∀ρ. Con este Ansatz la densidad de partículas
del condensado viene dada por f2(ρ) y la velocidad por ∇g(ϕ). Se debe tener en cuenta que en
coordenadas polares:
∇2 ≡ 1
ρ
∂ρ (ρ∂ρ)+
1
ρ2
∂2ϕ, ∇2ψ =
{
∂2ρf(ρ) +
1
ρ
∂ρf(ρ)−
1
ρ2
[
(∂ϕg(ϕ))2 − i∂2ϕg(ϕ)
]
f(ρ)
}
eig(ϕ).
(4.7)
Sustituyendo el resultado anterior en la Ec. (4.6) y separando la parte real e imaginaria se llega
a dos ecuaciones:
1
2∂
2
ρf(ρ) +
1
2ρ∂ρf(ρ) +
[
1− (∂ϕg(ϕ))
2
2ρ2
]
f(ρ)− f3(ρ) = 0, (4.8)
1
2ρ2 f(ρ)∂
2
ϕg(ϕ) = 0. (4.9)
De la última ecuación vemos como g(ϕ) = a1ϕ+a2 siendo a1 y a2 constantes reales (despreciamos
la solución f(ρ) = 0 ∀ρ, ya que anularía la función de onda en cualquier punto del sistema). La
cuantización de la circulación impone a1 = l, mientras que no hay ninguna restricción para a2,
por simplicidad tomaremos a2 = 0. Así:
g(ϕ) = lϕ, l ∈ Z. (4.10)
Con esta solución podemos reescribir la Ec. (4.8) como:
∂2ρf(ρ) +
1
ρ
∂ρf(ρ) + 2
[
1− l
2
2ρ2
]
f(ρ)− 2f3(ρ) = 0. (4.11)
A continuación, se realiza un análisis cualitativo de la Ec. (4.11) con el fin de conocer el
comportamiento de la densidad de partículas en la presencia de vórtices. Para puntos cercanos al
núcleo del vórtice, ρ� 1, se pueden despreciar los términos 2f(ρ) y 2f3(ρ) frente a l2
ρ2 f(ρ) [9].
De esta manera, se llega a una ecuación diferencial de Cauchy-Euler para f(ρ):
ρ2∂2ρf(ρ) + ρ∂ρf(ρ)− l2f(ρ) = 0, (4.12)
19
4.2. Vórtices cuánticos
la cual se resuelve con el Ansatz f(ρ) = ρα y se obtiene que α = ±l con l un número natural. Así:
f(ρ) = c1ρl + c2ρ−l. (4.13)
Por tal de obtener una solución finita en el origen, f(ρ = 0) 9∞, debemos imponer c2 = 0, por
lo que:
f(ρ� 1) ∝ ρl , l ∈ N (4.14)
que establece que la densidad de partículas en el núcleo del vórtice se anula.
Por otro lado, para regiones lo suficiente lejanas del vórtice, ρ� 1, se espera que la función
de onda tienda a la soluciónhomogénea ψ = 1 (o ψ = ψ0 con variables dimensionales); lo que
implicaría f(ρ) = 1 (o f(ρ) = f0(ρ) =
√
µ/gB utilizando variables dimensionales). Al contrario
que en el caso anterior, para distancias tales que ρ� 1, los términos 2f(ρ) y 2f3(ρ) dominan
frente a l2
ρ2 f(ρ) en la Ec. (4.11), la cual se puede reescribir como[1
2∇
2 + 1− |f(ρ)|2
]
f(ρ) = 0, (4.15)
que no es más que la propia ecuación de Gross-Pitaevskii para la amplitud de la función de onda.
Tiene como solución:
f(ρ� 1) = 1, (4.16)
es decir, la densidad de partículas se anula en el centro del vórtice pero al alejarse lo suficiente,
la función de onda, y por tanto la propia densidad, recupera su valor homogéneo. Este aumento
progresivo de la densidad desde cero hasta su valor homogéneo se produce en distancias del orden
de la longitud de coherencia, ξ. Por lo tanto se deduce que el radio de los vórtices debe ser de
orden ξ.
Por lo que respecta a la velocidad del condensado, esta se calcula como:
v(ρ, ϕ) = ∇g(ϕ) = eρ ∂ρg(ϕ)︸ ︷︷ ︸
vρ=0
+eϕ
1
ρ
∂ϕg(ϕ)︸ ︷︷ ︸
vϕ
= eϕ
1
ρ
∂ϕ(lϕ) =
l
ρ
eϕ. (4.17)
Es decir, en la vecindad de un vórtice la única componente no nula del campo de velocidad
es la tangencial/azimutal. Esta velocidad es inversamente proporcional a la distancia al núcleo
del vórtice, por lo que justo en ρ = 0 la velocidad sería infinita y, consecuentemente, la energía
cinética también. Este comportamiento obliga a la densidad de partículas a anularse en el centro
de un vórtice por tal que no haya ninguna partícula capaz de adquirir esa energía infinita, es
decir: f2(ρ→ 0)→ 0. Condición que se cumple automáticamente dada la Ec. (4.14).
Una vez calculada la velocidad del BEC, se puede calcular su rotacional para obtener la
vorticidad:
ω = ∇× v = 1
ρ
∂ρ(ρvϕ)eϕ =
1
ρ
∂ρ
(
ρ
l
ρ
)
eϕ, (4.18)
que es igual a cero ∀ρ 6= 0. Es decir, el flujo del condensado es irrotacional excepto en el núcleo
de un vórtice, donde la vorticidad no está definida (tal y como se adelantó en la Sección 2.2). A
este tipo de vórtices se les llama irrotacionales.
4.2.2 Fase de la función de onda en un vórtice
De la Sección anterior se puede escribir ψ(ρ, ϕ) = f(ρ)eilϕ siendo lϕ ≡ S la fase de la función
de onda del condensado:
S = arctan
[=(ψ)
<(ψ)
]
, (4.19)
20
4.3. Potencial armónico
representando <(ψ) y =(ψ) la parte real e imaginaria, respectivamente, de la función de onda.
En el núcleo de los vórtices la densidad de partículas se anula, es decir f(ρ = 0) = 0. Entonces
<[ψ(ρ = 0, ϕ)] = f(ρ = 0) cos(S) = 0
=[ψ(ρ = 0, ϕ)] = f(ρ = 0) sin(S) = 0
y, por tanto, según la Ec. (4.19) la fase de la función de onda no está definida. Este hecho es
equivalente a afirmar que la vorticidad en el núcleo de un vórtice no está definida puesto que
ω = ∇× v = ∇×∇S donde S tampoco lo está. Este último argumento reafirma el obtenido a
partir de la Ec. (4.18).
4.3 Potencial armónico
Consideraremos un condensado de Bose-Einstein en el seno de un potencial armónico con-
finante, V (x, y) = 12Ω
2
ρ(x2 + y2) (en unidades adimensionales). En esta Sección se estudiará el
estado transitorio por el que pasa la función de onda hasta alcanzar el régimen estacionario, en
el cual adopta la forma parabólica del potencial.
4.3.1 Relajación de la función de onda
Para un condensado bidimensional sometido a un potencial armónico, sin haz láser (V0 = 0) e
inhomogéneo (∇2ψ 6= 0), se resuelve numéricamente la Ec. (4.2). Tal y como muestra la Fig. 4.1a,
el potencial armónico causa oscilaciones de los perfiles de la densidad y la fase durante un régimen
transitorio. Tras t ≈ 400 (en unidades de ξ/c) cesan las oscilaciones y el régimen pasa a ser
estacionario, donde ambas magnitudes son constantes e iguales a sus valores de equilibrio: n = 1
y S = constante arbitraria. El perfil de densidad adopta una forma de paraboloide invertido, al
contrario que el potencial armónico ya que las partículas tienden a desplazarse hacia zonas de
equilibrio, es decir, de menor potencial (ver Fig. 4.1b).
0.8
1.0
1.2
1.4
n
(0
,0
,t
)
0 100 200 300 400 500 600
t
−π
−π
2
0
π
2
π
S
(0
,0
,t
)
Evolución de la densidad y la fase
(a)
x
60 30 0 30 60
y
60
30
0
30
60
0
0.5
1
Densidad en régimen estacionario
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
n
(r
,t
)
(b)
Figura 4.1: (a) Relajación de la densidad de partículas del BEC y la fase de la función de
onda macroscópica calculados en el punto (0, 0). (b) Representación de la densidad en el estado
estacionario. También se grafican las proyecciones 2D y 1D. Para la simulación se ha utilizado
Ωρ = 0.025 y n(r, t = 0) = 0.85.
21
4.4. Sin potencial externo
4.3.2 Aproximación de Thomas-Fermi
Para el sistema de la Sección anterior se puede obtener un resultado analítico aproximado de
la densidad de partículas en el estado estacionario.
Con esta finalidad, asumiremos que el flujo es aproximadamente homogéneo (∇2ψ ≈ 0), es
decir, el término de presión cuántica es despreciable respecto a los otros términos. El objetivo es
conocer la forma que presentará n(x, y, t) en el estado estacionario, por lo que ∂tψ = 0. Para este
sistema, la Ec. (4.2) se reduce a
0 = (i+ γ)
[
−|ψ|2 + 1− 12Ω
2
ρ(x2 + y2)
]
ψ. (4.20)
Con la relación |ψ(x, y, t)|2 = n(x, y, t) y evitando la solución trivial ψ = 0, se obtiene
n(x, y, t) = 1− 12Ω
2
ρ(x2 + y2), (4.21)
60 40 20 0 20 40 60
x
0.0
0.2
0.4
0.6
0.8
1.0
V
(x
),
n
(x
),
n
T
F
(x
)
Aprox. de Thomas-Fermi
n(x)
nTF(x)
V(x)
Figura 4.2: Potencial armónico (ver-
de), densidad calculada numéricamen-
te (azul) y densidad aproximada (rojo)
en función de la coordenada x. Se ha
utilizado y = 0 y Ωρ = 0.025.
que tiene forma de paraboloide invertido. Los puntos en
los que la densidad de partículas se vuelve nula crean
una circunferencia cuyo radio es conocido como el radio
de Thomas-Fermi (RTF). Imponiendo n(x, y, t) = 0 en
la Ec. (4.21) se llega a
RTF2 ≡ x2 + y2 = 2Ω2ρ
⇒ RTF =
√
2
Ωρ
, (4.22)
que muestra la relación de proporcionalidad inversa entre
RTF y la frecuencia confinante.
Cuando 12Ω
2
ρ(x2 + y2) > 1 la Ec. (4.21) da lugar
a valores negativos para la densidad de partículas, los
cuales carecen de sentido físico. Por tal de corregir estos
valores se toma:
nTF (x, y, t) =
{
1− 12Ω
2
ρ(x2 + y2) si x2 + y2 < RTF,
0 si x2 + y2 ≥ RTF,
(4.23)
con RTF dado por la Ec. (4.22) y nTF indica la densidad
de partículas con la aproximación de Thomas-Fermi. Es-
te perfil de densidad indica que en el estado estacionario,
todas las partículas que forman el condensado se encuen-
tran confinadas en un círculo de radio aproximadamente
igual al radio de Thomas-Fermi.
En aras de simplificación, la Fig. 4.2 muestra la de-
pendencia de V , n y nTF con x para y = 0 (corresponde
con una de las proyecciones unidimensionales de la Fig. 4.1b). La densidad de partículas calculada
numéricamente y la obtenida a través de la Ec. (4.23) son casi idénticas en todo el dominio de x
excepto cerca de x = ±56.57 que equivale al valor de RTF para Ωρ = 0.025, donde n presenta un
perfil suave que tiende a cero.
4.4 Sin potencial externo
En esta Sección estudiaremos el caso de un condensado de Bose-Einstein sin potencial externo.
Para ello impondremos que la frecuencia de confinamiento adimensional Ωρ es igual a cero
en la Ec. (4.2). En primer lugar, se analizará la relajación de la función de onda inicial a su
22
4.4. Sin potencial externo
valor homogéneo en ausencia de haz láser. Seguidamente, tras introducir una impureza en el
condensado, se resolverá numéricamente la ecuación de Gross-Pitaevskii disipativa dando lugar
a perfiles de densidad, fase y campos de velocidad en cuatro regímenes diferentes. Finalmente,
para uno de esos regímenes, mediante la descripción lagrangiana, se calculará la trayectoria de
un conjunto de pares de partículas del fluido y se estudiará si la dinámica es caótica.
4.4.1 Relajación de la función de onda
Se considera un condensado totalmente homogéneo (∇2ψ = 0), por lo tanto sin ningún haz
láser que produzca vórtices (V0 = 0) y sin potencial externo (Ωρ = 0). Con estas condicionesla
Ec. (4.2) queda
∂ψ
∂t
= (i+ γ)
(
1− |ψ|2
)
ψ. (4.24)
Además, se asume que inicialmente la función de onda es real y constante:
ψ(r, t = 0) =
√
α, α > 0 & α 6= 1, con r = (x, y), (4.25)
de tal forma que la densidad inicial de partículas es n(r, t = 0) = α. Puesto que el objetivo de
esta sección es analizar la relajación de la función de onda a sus valores homogéneos, se excluye
α = 1, donde no habría relajación.
Haciendo uso de la transformación de Madelung, Ec. (2.10), en la Ec. (4.24) y separando en
parte real e imaginaria se llega a
1
2n
∂n
∂t
= γ(1− n), (4.26)
∂S
∂t
= 1− n, (4.27)
que son las ecuaciones de evolución temporal para la densidad y la fase del condensado, respecti-
vamente. A medida que n tiende al valor homogéneo, n→ 1, la evolución de S tiende a cero. Es
decir, pasado un cierto tiempo, tanto la densidad como la fase del condensado tienden a valores
constantes. Estos valores constantes se obtienen solucionando las dos ecuaciones anteriores y
tomando el límite para tiempos tendiendo a infinito. En primer lugar, las soluciones de las Ecs.
(4.26) y (4.27) son
n = χe
2γt
1 + χe2γt , (4.28)
S = t− 12γ ln
∣∣∣∣∣1 + χe2γt1 + χ
∣∣∣∣∣ , (4.29)
donde se ha definido χ ≡ α/(1− α). Tomando el límite t→∞:
ĺım
t→∞
n = 1, (4.30)
ĺım
t→∞
S = − 12γ ln
∣∣∣∣ χ1 + χ
∣∣∣∣ , (4.31)
que determinan los valores asintóticos a los que tienden la densidad de partículas y la fase de
la función de onda del condensado. En el caso en que 0 < α < 1 la dinámica debería hacer que
n creciese hasta 1. Sin embargo, si α > 1, entonces n debería disminuir hasta alcanzar el valor
estacionario.
El comportamiento descrito por las Ecs. (4.28)-(4.31) se ilustra en la Fig. 4.3. Para los
gráficos se han escogido dos condiciones iniciales para n: una por encima del valor de equilibrio,
n(r, t = 0) = 1.15, y otra por debajo, n(r, t = 0) = 0.85 (en unidades de n0 = |ψ0|2 = µ/gB). La
evolución de ambos perfiles de densidad sigue la Ec. (4.26), de donde se infiere que si inicialmente
n < 1 entonces ∂tn > 0 y la densidad crecerá hasta alcanzar n = 1, valor para el cual ∂tn = 0.
23
4.4. Sin potencial externo
0.85
0.95
1.05
1.15
n
(r
,t
)
A
Evolución de la densidad y la fase
0 10 20 30 40 50 60 70 80
t
−π
−π
2
0
π
2
π
S
(r
,t
)
B
n(r, t= 0) = 1.15
n(r, t= 0) = 0.85
Figura 4.3: Representación del perfil de densidad de partículas n(r, t) (Panel A) y de la fase S(r, t)
(Panel B) en función del tiempo para dos condiciones iniciales distintas. Las líneas continuas
representan los perfiles analíticos dados por las Ecs. (4.28) y (4.29), mientras que las estrellas
y triángulos se han obtenido a partir de la resolución numérica de la Ec. (4.24). También se
ilustran con líneas discontinuas los valores asintóticos a los que tienden n y S, dados por las Ecs.
(4.30) y (4.31). Se ha utilizado γ = 0.03.
Y lo contrario sucede si inicialmente la densidad excede el valor n = 1. Por su parte, la fase
comienza de un valor inicial nulo y evoluciona siguiendo la Ec. (4.27). En t = 0 la fase del
condensado es cero ya que se ha elegido una condición inicial para la función de onda tal que
ψ(r, t = 0) ∈ R, por lo que =[ψ(r, t = 0)] = 0 y, consecuentemente, la arcotangente del cociente
de la parte imaginaria y real de ψ es idénticamente cero. Pasado t ≈ 60 ambos perfiles de fase
tienden a sus valores estacionarios dados por la Ec. (4.31).
En la siguiente sección, se estudiarán flujos más complejos, con la aparición de vórtices que
afectarán notablemente a los perfiles de densidad y fase, los cuales ya no serán homogéneos por
todo el sistema.
4.4.2 Regímenes
En esta Sección se estudia un condensado con un haz láser de radio σ siguiendo un movimiento
circular de radio s y con una velocidad lineal de módulo v`. De esta forma se crean inhomegenidades
(∇2ψ 6= 0) y vórtices en el fluido. Igual que en la Sección anterior, asumimos un potencial externo
nulo (Ωρ = 0).
La velocidad del láser resulta, intiuitivamente, un parámetro directamente relacionado con
la creación de inhomogeneidades en el flujo. Es decir, a mayor velocidad, se espera un mayor
número de vórtices. Por este motivo, se ha simulado la ecuación disipativa de Gross-Pitaevskii
variando v` en un amplio rango de velocidades. Al igual que en la Ref. [12], tras analizar los
resultados de dichas simulaciones se distinguen cuatro regímenes de turbulencia en función del
número de vórtices emitidos por el haz, mostrados en las Figs. 4.4-4.7.
El régimen de emisión cero se alcanza para v` . 0.5 (en unidades de c, la velocidad del
sonido en el medio), donde el haz no tiene suficiente velocidad como para inducir vórtices en
el condensado, ver Figura 4.4. Este régimen de flujo casi laminar no solo se obtiene con bajas
velocidades del haz, sino que también con valores de V0 no lo suficientemente grandes como para
anular la densidad de partículas en el núcleo de un hipotético vórtice. Es decir, la intensidad de
24
4.4. Sin potencial externo
Figura 4.4: Densidad de partículas (1a fila) y fase (2a fila) en régimen de emisión cero para
distintos tiempos. Para la simulación se ha utilizado: s = 10, v` = 0.3, a = 1, V0 = 10 y Ωρ = 0.
Se ha tomado una ventana de área (20ξ)2 pese a que la simulación se hace en una de (128ξ)2.
La dirección del campo de velocidades se indica con flechas normalizadas, cuyo módulo se lee
de la barra de colores situada en la parte inferior de cada gráfico. El tiempo y la velocidad en
variables dimensionales se pueden obtener multiplicando t y |v| por t0 = 100 µs y c = 2.26 mm/s,
respectivamente.
Figura 4.5: Densidad de partículas (1a fila) y fase (2a fila) en régimen dipolar para distintos
tiempos. En la simulación, se han utilizado los mismos parámetros que para el régimen dipolar
excepto por v` = 0.5. La dirección del campo de velocidades se indica con flechas normalizadas,
cuyo módulo se lee de la barra de colores situada en la parte inferior de cada gráfico. El tiempo
y la velocidad en variables dimensionales se pueden obtener multiplicando t y |v| por t0 = 100 µs
y c = 2.26 mm/s, respectivamente.
25
4.4. Sin potencial externo
Figura 4.6: Densidad de partículas (1a fila) y fase (2a fila) en régimen de clústers para distintos
tiempos. Se ha utilizado v` = 0.7. La dirección del campo de velocidades se indica con flechas
normalizadas, cuyo módulo se lee de la barra de colores situada en la parte inferior de cada
gráfico. El tiempo y la velocidad en variables dimensionales se pueden obtener multiplicando t y
|v| por t0 = 100 µs y c = 2.26 mm/s, respectivamente.
Figura 4.7: Densidad de partículas (1a fila) y fase (2a fila) en régimen de solitones oblicuos
para distintos tiempos. En la simulación, se han utilizado los mismos parámetros que para los
regímenes anteriores excepto por v` = 1.5. La dirección del campo de velocidades se indica con
flechas normalizadas, cuyo módulo se lee de la barra de colores situada en la parte inferior de
cada gráfico. El tiempo y la velocidad en variables dimensionales se pueden obtener multiplicando
t y |v| por t0 = 100 µs y c = 2.26 mm/s, respectivamente.
26
4.4. Sin potencial externo
la interacción entre el haz y las partículas que forman el BEC es débil y el potencial no es lo
suficientemente repulsivo como para crear regiones donde n = 0.
Alrededor de v` w 0.5 y con V0 ≈ 10, Fig. 4.5, comienzan a aparecer pares de vórtices de
manera periódica. Este es el régimen dipolar, donde se emiten parejas de vórtices aisladas, las
cuales prácticamente no interactúan unas con otras. Por otro lado, los constituyentes de un
mismo par son creados con circulaciones opuestas, lo que conlleva a una atracción mutua hasta
su posterior aniquilación.
Si se sigue incrementando la velocidad del haz láser, el condensado entra en el régimen de
clústers, Fig. 4.6. Para velocidades en el intervalo 0.5 . v` . 1 la frecuencia de emisión de los
pares de vórtices aumenta y, consecuentemente, la interacción entre ellos. Cuando dos parejas se
aproximan pueden intercambiar uno de sus constituyentes, modificar sus trayectoriaso colisionar
frontalmente.
Finalmente, para velocidades superiores a la velocidad del sonido en el medio, es decir v` & 1,
se observa una cadena de vórtices tras el haz y un gran frente de onda en la parte delantera,
ver Fig. 4.7. Este es el llamado régimen de solitones oblicuo. Su nombre se debe a la creación
de solitones oscuros oblicuos que, en un condensado bidimensional, se vuelven dinámicamente
inestables conduciendo a cadenas de pares de vórtices [12, 14]. En este régimen la vida de los
pares de vórtices es muy corta ya que la aniquilación se produce pocos instantes después de la
creación de estos.
En las Figs. 4.4-4.7 también se muestra el campo de velocidad, calculado con la Ec. (2.8), del
fluido para distintos tiempos (normalizados a t0 ∼ 100 µs). Las flechas indican la dirección del
campo, mientras que el módulo se indica mediante el color de estas, mostrado en la parte inferior
de cada gráfico. Como se espera de la Ec. (4.17), la velocidad en el entorno de los vórtices es
tangencial a ellos y cuyo módulo aumenta a medida que se reduce la distancia al núcleo. En el
resto del condensado, es decir lejos de los vórtices, el fluido es aproximadamente homogéneo,
por lo que la densidad de partículas tiende a su valor estacionario n = 1, como se derivó en la
Sección 4.4.1. Velocidades supersónicas pueden originarse en los núcleos de los vórtices, donde no
hay partículas; sin embargo, en las posiciones restantes la velocidad no excede el límite c = 1.
Por otro lado, los gráficos de la fase del condensado verifican que en el centro de un vórtice S
no está definida, observándose un continuo de colores (en concordancia con lo predicho en la
Sección 4.2.2).
4.4.3 Exponente de Lyapunov
El caos mide lo sensible que es un sistema a cambios en las condiciones iniciales [15]. Muchos
procesos disipativos de la física estadística y no lineal están descritos por una dinámica caótica
[16]. Una de las maneras de cuantificar el caos es mediante el exponente de Lyapunov, que
mide, en promedio, la divergencia o convergencia [16] de las órbitas de dos partículas separadas
inicialmente por una distancia idealmente infinitesimal. Un exponente positivo es un claro indicio
de un sistema caótico.
Consideremos dos partículas de un BEC bidimensional separadas inicialmente por una
distancia ‖δ(t = 0)‖. Es decir:
‖δ(t = 0)‖ =
√
[x02 − x01 ]2 + [y02 − y01 ]2 (4.32)
donde (x0i , y0i), con i = 1, 2, son las coordenadas iniciales de la partícula “i” y ‖ · ‖ es la norma
Euclideana. Para calcular el exponente de Lyapunov se supone que la separación entre ambas
partículas evoluciona como
‖δ(t)‖ ≈ eλt‖δ(0)‖, (4.33)
siendo λ el exponente de Lyapunov y donde la validez de la expresión requiere que ‖δ(0)‖ sea
suficientemente pequeño. Si λ < 0 ambas trayectorias convergen a una misma pero si λ > 0 la
separación entre ellas se acentúa exponencialmente, dando paso a un régimen caótico donde el
27
4.4. Sin potencial externo
comportamiento del sistema es extremadamente alterado por un sutil cambio en las condiciones
iniciales.
En esta Sección trataremos de calcular el exponente de Lyapunov para un condensado de
Bose-Einstein sin potencial externo y en el régimen de clústers. Para ello se ha adoptado la
descripción lagrangiana del fluido cuántico para calcular las trayectorias de 50 pares de partículas
a modo de “tracking". La separación inicial de las partículas de cada par satisfacen la Ec. (4.32) y
se han colocado de manera aleatoria por el sistema. Por tal de describir una situación estacionaria,
se comienza a recoger datos de las posiciones de las partículas cuando el número de vórtices es
aproximadamente constante, es decir cuando el número de vórtices creados es similar al número de
vórtices aniquilados. A partir de sus respectivas posiciones iniciales, la posición de cada partícula
se calcula mediante
d
dtri(t; r0i , t0i) = v[ri(t; r0i , t0i), t], (4.34)
siguiendo con la notación introducida en la Sección 1.1. En el régimen de clústers, para tiempos
mayores a 3P , siendo P el periodo de la órbita del haz láser, el número de vórtices es aparentemente
constante por lo que se establece como tiempo de referencia t0i = 3P ∀i.
La Ec. (4.34) se soluciona con una modificación al algoritmo de Runge-Kutta para funciones
discretas (ver Sección A.2 del Apéndice A). Para cada paso de integración se calcula la distancia
entre los constituyentes del par y, posteriormente, se realiza el promedio de los resultados de
todas las parejas ‖δ(t)‖. Para obtener λ, se lleva a cabo un ajuste lineal de los datos
log
(
‖δ(t)‖
‖δ(0)‖
)
= at+ b, (4.35)
donde la pendiente de la recta, a, se identifica con el exponente de Lyapunov y b es la ordenada
en el origen.
El perfil de la distancia entre las partículas de una pareja, en promedio, está caracterizado por
tres zonas: transitoria, crecimiento exponencial y saturación [17]. En la primera de ellas, ‖δ(t)‖
permanece aproximadamente constante. Para nuestro caso, esta zona corresponde al tiempo que
tarda el haz láser o un vórtice en pasar cerca del par de partículas y modificar la distancia entre
ellas. En el periodo de crecimiento exponencial el par se separa cada vez más y es en esta zona
donde el ajuste lineal de la Ec. (4.35) es válido. Es decir, la separación determina el valor de la
pendiente a de la regresión y, por tanto, determina λ. Finalmente, la distancia tiende a un valor
finito. Este efecto se produce cuando las partículas se encuentran en posiciones más allá del radio
de la órbita del haz láser, donde la concentración de vórtices se reduce notablemente debido a
procesos de aniquilación y la separación entre pares permanece aproximadamente constante.
Los resultados numéricos se encuentran en la Fig. 4.8. El ajuste lineal se realiza entre t ≈ 305
y t ≈ 365, que corresponde a la zona de crecimiento exponencial para el régimen de clústers. El
valor del exponente de Lyapunov obtenido es λ = 0.019, por lo que se puede concluir que, en
promedio, la dinámica lagrangiana del sistema es caótica.
28
4.4. Sin potencial externo
260 280 300 320 340 360 380 400 420
t
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
lo
g[
δ(
t)
/δ
(0
)]
R 2 = 0.99334
log[δ(t)/δ(0)] = − 5.609 + 0.019t
Exponente de Lyapunov promedio (N=100)
Ajuste lineal
Numérico
Figura 4.8: Promedio de la distancia entre los constituyentes de 50 pares de partículas moviéndose
en el régimen de clústers de un BEC,N = 100 partículas del fluido. Se ha utilizado una separación
inicial δ(0) = 1 = 4dx (en unidades de ξ) entre las parejas y una distribución inicial aleatoria en
un cuadrado de dimensiones (10ξ)2.
29
Capı́tulo 5
Conclusiones
En este trabajo se han analizado las fuerzas que actúan sobre una impureza, así como el
campo de velocidad y flujo turbulento producido por el movimiento circular de un haz láser en
un condensado de Bose-Einstein. Para ello, en el Capítulo 2 se ha descrito el BEC como un fluido
cuántico mediante la transformación de Madelung y, a partir de la ecuación de Gross-Pitaevskii
se han derivado los análogos de las ecuaciones de continuidad y de Euler de la dinámica de fluidos
clásica en el marco cuántico.
Para un condensado ideal, cilíndricamente confinado y sometido a un potencial gravitatorio
externo, existe un gradiente de densidad que provoca una fuerza de flotación sobre cualquier
impureza que se introduzca en él. Esta fuerza de flotación cuántica, Ec. (3.6), se compara con
la fuerza de Arquímedes, Ec. (3.7), y se deduce que gp hace el mismo papel que el volumen
de la partícula clásica sumergido en el fluido. Para este sistema, el perfil de densidad, n, y de
fuerza, Fp, presentan no linealidades cerca de la base y de la parte superior del cilindro. En
cambio, para el caso simple donde se desprecia la presión cuántica, n es lineal y da lugar a una
fuerza de flotabilidad constante con z. Pese a no ser este un resultado físicamente posible, ya
que la densidad no se anula en las fronteras del volumen, resulta una muy buena aproximación
en posiciones lo suficientemente alejadas de los bordes,tal y como se ilustra en las Figs. 3.1 y
3.2. También se concluye que a medida que aumenta el valor de la gravedad (proporcional al
parámetro K) aumenta la curvatura de n, ya que el término de la presión cuántica en la Ec. (2.12)
comienza a ser apreciable, por lo que la aproximación deja de ser válida a partir de K ≈ 0.15.
En el Capítulo 4 se demuestra que los vórtices originados por el haz láser son puntos del
condensado donde la densidad de partículas es nula y la fase no está definida, como sucede con
los vórtices irrotacionales de la dinámica de fluidos clásica. Consecuentemente, son las únicas
zonas donde el fluido cuántico no puede considerarse irrotacional. Posteriormente, se estudia
la relajación del perfil de densidad de un condensado no ideal bajo la acción de un potencial
armónico y sin potencial externo. En ambas situaciones se concluye que, gracias a la interacción
con un baño térmico, existe un valor estacionario para la densidad dado por n = 1.
Para el caso sin potencial externo se resuelve numéricamente la Ec. (4.2) con Ωρ = 0.025
y V0 = 10 dando lugar a cuatro regímenes en función de v` y del número de vórtices creados:
emisión cero (v` . 0.5), dipolar (v` ' 0.5), de clústers (0.5 . v` . 1) y de solitones oblicuos
(v` & 1). Para el mismo potencial externo y en el régimen de clústers, también se calcula el
exponente de Lyapunov concluyendo con la existencia de una dinámica caótica en promedio para
el movimiento de las partículas del fluido.
Finalmente y como parte de un futuro trabajo, se podría integrar numéricamente la Ec. (3.9)
por tal de obtener correcciones a primer orden de la fuerza que sufre una impureza inmersa en
un BEC. Además, la misma clasificación en cuatro regímenes que se ha llevado a cabo para
el caso sin potencial externo podría elaborarse con el potencial armónico. De esta manera se
podría concluir si la forma del perfil de densidad (plano o como un paraboloide invertido para
30
la primera y segunda situación, respectivamente) influye en el número de vórtices creado. Por
último, comparar los distintos exponentes de Lyapunov para los tres regímenes y analizar cuáles
son los parámetros cruciales en su determinación, sería otro de los posibles trabajos de futuro.
Agradecimientos
Querría dar las gracias tanto a Cristóbal como a Emilio por sus innumerables explicaciones,
su paciencia y su tiempo.
31
Apéndice A
Algoritmo numérico
En este Apéndice se explica brevemente el método numérico utilizado para la simulación de
la ecuación disipativa de Gross-Pitaevskii. También, se comenta la modificación del algoritmo
de Runge-Kutta para poder aplicarlo a un conjunto discreto de datos en lugar de a funciones
analíticas.
A.1 Simulación dGPE
Las soluciones numéricas de la Ec. (4.2) se simulan en un sistema de tamaño (128ξ)2 con un
paso de integración espacial dx = dy = 0.25ξ y temporal de dt = 0.01t0. Se usa s = 10 como
radio de la trayectoria circular del haz láser, por lo que la dinámica del condensado se concentra
en una ventana de aproximadamente (20ξ)2 y se pueden despreciar los efectos de borde. El resto
de parámetros característicos de la Ec. (4.2) utilizados son: γ = 0.03, V0 = 10 y a = 1, mientras
que Ωρ y v` varían dependiendo del sistema que se quiera analizar.
La ecuación que se quiere resolver es
∂ψ(x, y, t)
∂t
= −(i+ γ)Hψ(x, y, t), (A.1)
dondeH ≡ T+V con T = 12∇
2 y V = −|ψ|2+1−Ω2ρ(x2+y2)/2−V0e
(x−xp)2+(y−yp)2
2a2 . Aproximando
el Hamiltoniano como independiente del tiempo1 e integrando entre t y t+ ∆t [1]:
ψ(x, y, t+ ∆t) ≈ e−(i+γ)∆tHψ(x, y, t). (A.2)
Puesto que T y V no conmutan e−(i+γ)∆tH 6= e−(i+γ)∆tT e−(i+γ)∆tV . Consecuentemente, se utiliza
la aproximación
e−(i+γ)∆tHψ(x, y, t) ≈ e−
1
2 (i+γ)∆tV e−(i+γ)∆tT e−
1
2 (i+γ)∆tV ψ(x, y, t), (A.3)
cuyo error es O(∆t3). La acción de las exponenciales de los operadores de energía potencial,
V , y de energía cinérica, T , sobre la función de onda puede ser computacionalmente costoso.
Uno de los algoritmos numéricos que simplifica notablemente la simulación de la Ec. (4.2) es
el llamado Beam Propagation Method [18], el cual consiste en hacer actuar V sobre la función
de onda en el espacio de posiciones y T en el espacio de momentos, donde son diagonales
respectivamente. De esta manera, la acción de V simplemente consiste en multiplicar ψ por la
fase que aparece en la Ec. (A.3) evaluando V (x, y) en cada punto de la malla bidimensional, es
decir e−
1
2 (i+γ)∆tV (x,y)ψ(x, y, t). Por su parte, para T se debe primero pasar a la representación
1en realidad no es cierto porque xp e yp dependen del tiempo. Sin embargo, con ∆t = 0.01t0 se cumple que
∆t� P , siendo P el periodo del láser, y así el movimiento del haz es despreciable en ese intervalo de tiempo.
32
A.2. Obtención de trayectorias
de momentos por lo que ∇2 → (ik)2 y ψ(x, y, t)→ F [ψ(x, y, t)], donde k es el vector de onda en
dos dimensiones y la letra F indica la transformada de Fourier en dos dimensiones. Al igual que
para V , en el espacio de momentos la acción de T sobre ψ se reduce a e−
1
2 (i+γ)∆t(ik)
2F [ψ(x, y, t)].
Así, el algoritmo se expresa como
ψ(x, y, t+ ∆t) ≈ e−
1
2 (i+γ)∆tV (x,y)F−1
[
e−
1
2 (i+γ)∆t(ik)
2F [e−
1
2 (i+γ)∆tV (x,y)ψ(x, y, t)]
]
, (A.4)
que permite solucionar la Ec. (4.2) dada una condición inicial para la función de onda. Siempre y
cuando se cumpla que (T + V )dt � 2π el algoritmo es convergente y, además, se conserva la
norma de ψ ya que V es real [18]. Para calcular F y F−1 en dos dimensiones se han utilizado las
funciones de Python numpy.fft.fft2 y numpy.fft.ifft2, respectivamente.
A.2 Obtención de trayectorias
Para integrar los campos de velocidad en la Sección 4.4.3 y obtener las trayectorias de algunas
partículas del fluido, se ha solucionado la Ec. (4.34) mediante el algoritmo de Runge-Kutta en
2D combinado con una interpolación en 3D (2D espacial + 1D temporal).
Figura A.1: Representación esque-
mática de la matriz de datos tridi-
mensional para v.
La Ec. (4.34) se podría resolver sin complicaciones con
el conocido algoritmo de Runge-Kutta en dos dimensiones
(RK2D) si la velocidad v fuese una función continua y ana-
lítica. Sin embargo, v se calcula numéricamente a partir de
dGPE, por lo que toma valores discretos sobre cada uno
de los puntos de una malla computacional. En concreto,
la malla utilizada es tridimensional y tiene un tamaño de
1024×1024×len(t), siendo len(t) la longitud del vector tempo-
ral, que da paso a la tercera dimensión, tal y como se observa
en el esquema de la Fig. A.1. Consecuentemente es necesaria
una modificación de RK2D si lo que se quiere es solucionar
la Ec. (4.34).
El objetivo es tratar de aproximar una matriz de datos
tridimensional por una función, f , de tres variables x, y, t de
tal forma que para cualquier punto (x0, y0, t0), pese a no coincidir con un punto de la malla
computacional, obtenga una interpolación de los datos almacenados en la matriz de velocidades.
Para ello se ha usado scipy.interpolate.RegularGridInterpolator una función de la librería scientific
python, que es capaz de interpolar linealmente sobre una malla computacional en dimensión
arbitraria [19]. Esta función juntamente con el algoritmo de Runge-Kutta de cuarto orden, con
un paso de integración h = dt = 0.01 (en unidades adimensionales), permiten solucionar la Ec.
(4.34).
33
Apéndice B
Deducción de la relación entre fase y
velocidad
A continuación se deriva la Ec. (2.11) partiendo de la transformación de Madelung, Ec. (2.10).
Despejando la fase, S, de la Ec. (2.10) se obtiene:
S = 1
i
ln
(
ψ√
n
)
. (B.1)
Utilizando la definición de la densidad de partículas se tiene que
√
n = |ψ| =
√
ψψ∗, por tanto la
Ec. (B.1) se reescribe de la siguiente manera
S = 12i (lnψ − lnψ
∗) . (B.2)
El siguiente paso consiste en aplicar el operador ∇ a la ecuación anterior y operar algebraicamente
como sigue:
∇S = 12i (∇ lnψ −∇ lnψ
∗)
= 12i
(
∂ lnψ
∂ψ
∇ψ − ∂ lnψ
∗
∂ψ∗
∇ψ∗
)
= 12i
(∇ψ
ψ
− ∇ψ
∗
ψ∗
)
= 12i
ψ∗∇ψ − ψ∇ψ∗
ψψ∗
= 1
n
=(ψ∗∇ψ) .
(B.3)
Resultado que es igual a la velocidad