Transiciones de fase cuánticas en un modelo de dos niveles para la coexistencia átomo-diátomo

TRABAJO FIN DE GRADO Grado en Física Transiciones de fase cuánticas en un modelo de dos niveles para la coexistencia átomo-diátomo Alumno: Ignacio Baena Jiménez _________________________________________ Tutores: Arias Carrasco, José Miguel Rodríguez Gallardo, Manuela Índice 1. Introducción ...................................................................................................................................... 1 1.1 Fases y transiciones de fase ......................................................................................................... 1 1.2 Transiciones de fase clásicas ....................................................................................................... 1 1.3 Transiciones de fase cuánticas .................................................................................................... 4 1.3.1 Incompatibilidad con los conceptos clásicos ....................................................................... 4 1.3.2 Redefinición de fase y transición de fase ............................................................................. 5 2. Modelo de dos niveles ....................................................................................................................... 7 2.1 Introducción teórica a los operadores del hamiltoniano .............................................................. 7 2.2 Hamiltoniano del problema ......................................................................................................... 9 2.3 Solución exacta del modelo ...................................................................................................... 10 2.3.1 Energía ................................................................................................................................ 10 2.3.2 Número de partículas en cada nivel ................................................................................... 11 2.3.3 IPR ...................................................................................................................................... 12 2.3.4 Entropía de Rényi ............................................................................................................... 13 2.4 Teoría de Campo Medio ............................................................................................................ 15 3. Resultados ....................................................................................................................................... 20 3.1 Dependencia con 𝑁 ................................................................................................................... 20 3.2 Dependencia de 𝜆𝑐 .................................................................................................................... 23 3.3 Teoría de campo medio ............................................................................................................. 26 3.4 Entropía ..................................................................................................................................... 31 4. Conclusiones ................................................................................................................................... 34 5. Referencias ...................................................................................................................................... 36 6. Apéndices ........................................................................................................................................ 38 6.1 Apéndice A: resolución matricial del problema ........................................................................ 38 6.2 Apéndice B: autofunciones en el espacio de configuración ...................................................... 39 6.3 Apéndice C: Reescritura de la energía en TCM ........................................................................ 41 5.4 Apéndice D: cálculo de 𝜆𝑐 ........................................................................................................ 43 6.5 Apéndice E: Códigos de MATLAB .......................................................................................... 46 6.5.1 Programa principal ............................................................................................................. 46 6.5.2 Operadores de creación y destrucción .............................................................................. 53 6.5.3 Polinomios de Hermite ....................................................................................................... 54 6.5.4 Modelo de campo medio ................................................................................................... 54 1 1. Introducción El objetivo de este proyecto es realizar un estudio sobre las transiciones de fase cuánticas en un modelo de dos niveles. Se trata de describir por mecanismos cuánticos la coexistencia de un gas de átomo-diátomo a temperatura absoluta cero utilizando un modelo de dos niveles que se resolverá con el programa 𝑀𝐴𝑇𝐿𝐴𝐵 (apéndice E). 1.1 Fases y transiciones de fase En un sistema físico, podemos definir fase como aquella en la que el sistema posee unas características uniformes. Por tanto, se puede asociar “fase” con un cierto comportamiento homogéneo del sistema, esto es, no solo su composición química sino su estado físico. Sin embargo, si cambian las condiciones que rodean al sistema de forma adecuada, este comportamiento homogéneo puede cambiar drásticamente. Se dice entonces que tenemos una transición de fase. Nótese que esta definición es deliberadamente ambigua. Es sólo en casos concretos cuando se adquiere cierta rigurosidad. A continuación, particularizamos este concepto para diversos sistemas físicos. 1.2 Transiciones de fase clásicas En una descripción clásica; por ejemplo, en el caso de estados de agregación de la materia, es sencillo determinar qué es una fase. Las características uniformes vienen dadas por parámetros como la densidad o la capacidad calorífica y las condiciones que rodean al sistema vienen a ser parámetros como la presión y la temperatura. Se define como parámetro o variable de control de una transición de fase como aquella que “controla”, que determina, en qué fase se encuentra el sistema [1]. En este proceso, dada una presión y una temperatura, podemos saber en qué fase nos encontramos. Otros ejemplos de transiciones de fase comunes son los cambios de conductor a superconductor o de paramagnético a ferromagnético. La fase vendría dada por la conductividad y por la imanación espontánea respectivamente y el parámetro de control por la temperatura en ambos casos [2]. Ya que la idea física de una transición de fase es un cambio abrupto en sus propiedades, deberá haber alguna magnitud física, de la que estas dependan, que sufra también un cambio repentino. 2 Esta es la idea detrás de la teoría de Ehrenfest [3], que explica las transiciones de fase clásicas como reza a continuación. Las transiciones de fase clásicas son aquellas que ocurren entre dos estados termodinámicos y están gobernadas por una competición entre la energía libre del sistema y la entropía asociada a las fluctuaciones térmicas [4]. Podemos definir entropía desde la mecánica estadística mediante los conceptos de microestado y macroestado [5]. Un microestado es una configuración microscópica de un sistema termodinámico, es decir, un punto su espacio de fases. Por el contrario, un macroestado es una descripción macroscópica, en la que existe al menos una magnitud extensiva, compatible con una ecuación de estado. La entropía se define como una función de estado que caracteriza el número de microestados compatibles con el macroestado de equilibrio. En una descripción termodinámica en la línea de Ehrenfest, la definición de entropía [6] se realiza en base al calor intercambiado en un proceso reversible y en base a la temperatura a la que sucede este intercambio: 𝑑𝑆=(𝑑𝑄 𝑇)𝑝𝑟𝑜𝑐𝑒𝑠𝑜 𝑟𝑒𝑣𝑒𝑟𝑠𝑖𝑏𝑙𝑒. (1.1) En un sistema clásico, esta competición no existe en el cero absoluto, ya que no hay entropía definida en él. Paul Ehrenfest clasificó las transiciones de fase clásicas por primera vez introduciendo el concepto de orden de transición, que va ligado al concepto de energía libre de Gibbs. La energía libre de Gibbs, 𝐺, es un potencial termodinámico que se puede usar para calcular el trabajo máximo en un proceso reversible isotérmico e isobárico de un sistema. El orden de una transición para Ehrenfest es la derivada de menor orden de la energía libre de Gibbs respecto a una variable termodinámica que presenta discontinuidad. Así, una primera derivada discontinua es indicativo de una transición de primer orden y una derivada discontinua de segundo orden lo es de una transición de segundo orden [7]. En nuestros ejemplos, las transiciones de primer orden serían las dadas por cambios de estado de agregación de la materia, ya que hay cambios discontinuos en la densidad, que tiene que ver con el inverso de la derivada de la energía libre respecto a la presión. Las derivadas primeras relevantes son: 3 (𝜕𝐺 𝜕𝑇)𝑝=−𝑆, (𝜕𝐺 𝜕𝑃)𝑇=𝑉. (1.2) Para las transiciones de segundo orden estamos hablando de transiciones de fase ferromagnéticas. La primera derivada de la energía libre respecto al campo aplicado es la magnetización, que crece continuamente hasta alcanzar la temperatura de Curie, y decrece para temperaturas superiores. Sin embargo, su segunda derivada presenta una discontinuidad. Sea el calor específico a presión constante 𝐶𝑃, la compresibilidad isotérmica 𝛽𝑇 y, el coeficiente de compresibilidad isotermo 𝜅𝑇, podemos dar expresiones de las segundas derivadas que mencionábamos [6]: (𝜕2𝐺 𝜕𝑇2)𝑃=−𝐶𝑃 𝑇, (𝜕2𝐺 𝜕𝑇 𝜕𝑃)𝑃=𝑉 𝛽𝑇, (𝜕2𝐺 𝜕𝑃2)𝑇=−𝑉 𝜅𝑇. (1.3) Existe un parámetro introducido por Ehrenfest, llamado parámetro de orden, que se comporta como la energía libre de Gibbs. El parámetro de orden es una medida del grado de orden del sistema. En una transición de fase, normalmente varía entre 0 y un valor distinto de cero [8]. Aunque esta clasificación es muy útil, no sirve para clasificar por completo las transiciones de fase, puesto que no toma en consideración aquellas transiciones en las que la derivada de la energía libre diverge. Esto ocurre con nuestro tercer ejemplo, en las transiciones ferromagnéticas, cuya capacidad calorífica diverge. Este fenómeno también es relevante en transiciones de fase a superconductor, cuya importancia no ha hecho más que aumentar en las últimas décadas. [9] Una clasificación más moderna distingue entre dos tipos de transiciones, en honor a Ehrenfest, de primer orden o discontinuas y las transiciones continuas. Las transiciones de primer orden, o discontinuas, son las que involucran un calor latente [6]: 𝐿=𝑇 ∆𝑆. (1.4) En estas transiciones el sistema absorbe o emite una cantidad (normalmente grande) de calor por unidad de volumen. Durante este proceso, la temperatura del sistema se mantiene constante y el parámetro de orden cambia discontinuamente en el punto crítico. En cambio, en las continuas, aparece una discontinuidad en la entropía, por lo que no podemos definir un calor latente [10] y el parámetro de orden cambia continuamente. 4 1.3 Transiciones de fase cuánticas Las transiciones de fase cuánticas son aquellas que ocurren a una temperatura absoluta nula. Por lo tanto, las transiciones no ocurren por fluctuaciones térmicas, sino por fluctuaciones cuánticas. 1.3.1 Incompatibilidad con los conceptos clásicos Ya que la temperatura es nula, no se puede definir un calor latente que nos clasifique la transición. Es más, ni siquiera podemos definir una energía libre que analizar. Como consecuencia, ninguno de los dos criterios de clasificación es aplicable. A lo que sí tenemos acceso es al hamiltoniano cuantizado que describe el sistema. ¿Puede un hamiltoniano describir una transición de fase? En un problema clásico sencillo como puede ser un vaso de agua que se evapora, el hamiltoniano no puede describir la transición de fase de líquido a vapor. ¿Dónde está la temperatura en un hamiltoniano? Y es que un hamiltoniano no es una energía libre. En cierto modo, un hamiltoniano no es más que un modelo que cuenta interacciones. No obstante, a veces uno puede disfrazar una energía libre de hamiltoniano. Por ejemplo, para un sistema de partículas de posiciones relativas 𝑟𝑖𝑗, momentos lineales 𝑝𝑖 y cargas eléctricas 𝑞𝑖: 𝐻=∑𝑝𝑖2 2𝑚𝑖 𝑖+∑𝑞𝑖 𝑞𝑗 4𝜋𝜀 𝑟𝑖𝑗 𝑖≠𝑗 . (1.5) En el vacío 𝜀=𝜀𝑜, pero en la materia 𝜀=𝜀(𝑇). En cierto modo, resolver el problema para una permitividad eléctrica está incluyendo la información de la temperatura. Pero esto no siempre (casi nunca) es posible. Para terminar, estamos obviando el problema más grave: según el tercer principio de la termodinámica, el cero absoluto es inaccesible. 5 1.3.2 Redefinición de fase y transición de fase Todo esto quiere decir que las transiciones de fase cuánticas no serían un análogo de las clásicas. Esto no se trata de la cuantización de un problema clásico en el que la temperatura caiga a cero. Si 𝐾𝐵 la constante de Boltzmann y ℏ la constante reducida de Planck, lo que sucede es que las fluctuaciones térmicas, del orden de 𝐾𝐵𝑇, compiten con las fluctuaciones cuánticas, provenientes del principio de incertidumbre, del orden ℏ𝜔, donde 𝜔 es la frecuencia del oscilador cuántico equivalente [11]. Para temperaturas lo suficientemente bajas, los efectos cuánticos dominan a los termodinámicos y son los responsables del comportamiento del sistema. Se habla entonces de fases y transiciones de fase cuánticas. La teoría es la siguiente: sea un sistema descrito por la cuantización de un hamiltoniano, una transición de fase ocurre cuando hay una discontinuidad en la derivada del hamiltoniano respecto a parámetros del sistema, como pueden ser la presión, el potencial químico o el campo magnético. Existe un tercer dominio, y sería aquel denominado punto crítico cuántico, en el que las fluctuaciones cuánticas son del orden de las termodinámicas. En este caso la descripción cambia mucho en función de las condiciones del problema concreto y es en general el más complicado. Este tipo de transiciones de fase son útiles en astrofísica, por ejemplo, para la teoría de agujeros negros, o para ciertos resultados de la magnetohidrodinámica relativista. La idea física que se pretende corroborar en este texto es que, en ausencia de temperatura, un hamiltoniano cuántico describe por completo al sistema físico. Por esto, un cambio abrupto en las propiedades físicas de un sistema va acompañado por discontinuidades en derivadas del hamiltoniano, al igual que ocurría con la energía libre en el caso de transiciones de la física clásica. Esto tiene cierto sentido, ya que en física clásica son los conceptos termodinámicos los que definen las fluctuaciones termodinámicas y es el hamiltoniano cuántico el que nos puede dar información sobre las fluctuaciones cuánticas. Y son unas y otras fluctuaciones en cada caso las responsables del fenómeno de las transiciones de fase. Esto nos permite redefinir el concepto de fase como aquel delimitado por un subconjunto conexo en el espacio funcional de las magnitudes que dependen de la energía y sus derivadas: llamémoslo espacio energético. 6 Sea un sistema en un estado 𝐴 y sea el mismo sistema en un estado 𝐵, si existe al menos un proceso en el que puede evolucionar de 𝐴 a 𝐵 de forma que las funciones del espacio energético son analíticas, continuas y derivables, se dice que 𝐴 y 𝐵 están en la misma fase. Esta definición delimita un subconjunto conexo por caminos en el que no aparecen cambios de fase. Es decir, delimita una fase. Nótese que en ningún momento hemos pedido que el estado esté dado por variables macroscópicas, conceptos estadísticos o que sea un estado cuántico. Se entiende por estado de un sistema a lo que determina a dicho sistema, y ya esto será una cosa u otra según la rama de la física desde la que se realiza su estudio. 7 2. Modelo de dos niveles En física cuántica un sistema de dos niveles (o dos estados) es un sistema que es compatible con una función de onda que se puede escribir como cualquier superposición de dos estados cuánticos independientes y distinguibles físicamente [12]. El ejemplo más común es el de un sistema de partículas de espín ½, por ejemplo, electrones, que interaccionan con un campo magnético externo. El estado de este sistema es una combinación lineal de estados con proyección de espín +½ y proyección de espín −½. En esta sección vamos a introducir nuestro modelo de dos niveles desde un punto de vista teórico-práctico. Es decir, vamos a explicar qué operadores, magnitudes y método de resolución vamos a usar. 2.1 Introducción teórica a los operadores del hamiltoniano En primer lugar, vamos a introducir los operadores de creación y destrucción de bosones. Sean 𝑎† y 𝑎 dos operadores (uno el conjugado del otro) que cumplen la relación de conmutación [𝑎,𝑎†]=1. Definamos el problema de autovalores de su producto: 𝑎†𝑎 |𝑛𝑎⟩=𝑛𝑎 |𝑛𝑎⟩, (2.1) donde 𝑛𝑎 es el autovalor y |𝑛𝑎⟩ es el autovector. A partir de su conmutador: [𝑎,𝑎†]|𝑛𝑎⟩=𝑎𝑎†|𝑛𝑎⟩−𝑎†𝑎|𝑛𝑎⟩, 1 |𝑛𝑎⟩=𝑎𝑎†|𝑛𝑎⟩−𝑛𝑎 |𝑛𝑎⟩, 𝑎𝑎†|𝑛𝑎⟩=(𝑛𝑎+1)|𝑛𝑎⟩. (2.2) Como conclusión, 𝑎† y 𝑎 cambian el estado cuántico de forma contraria de un autoestado del operador 𝑎†𝑎 a otro del mismo operador. Esto es debido a que aplicar ambos operadores en cualquier orden dejan al estado (salvo constante) sin cambiar. Por tanto, hay constantes asociadas a cada operador. Si realizamos unas sencillas comprobaciones, nos damos cuenta de que sólo es consistente que: 𝑎†|𝑛𝑎⟩=√𝑛𝑎+1 |𝑛𝑎+1⟩, (2.3) 𝑎 |𝑛𝑎⟩=√𝑛𝑎 |𝑛𝑎−1⟩. (2.4) 14 𝐼(𝑝1𝑝2)=𝐼(𝑝1)𝐼(𝑝2). Esto es consistente con la entropía de Hartley, que consideró que 𝐼(𝑝)= −ln𝑝. Ahora bien, la información total será el valor esperado de esta variable. Asumamos que tenemos 𝑀 eventos con probabilidades 𝑃={𝑝1,𝑝2…𝑝𝑀}. Si 𝐼𝑘=𝐼(𝑝𝑘), 𝐼(𝑃)=∑𝑝𝑘 𝐼𝑘 𝑀 𝑘=1 . (2.23) Esto puede ser definido como la entropía de Shannon. En nuestro caso, estamos sumando sobre una variable continua, la posición, de modo que debemos extender el sumatorio como una integral. La probabilidad estaría dada por el módulo cuadrado de la función de onda. Por tanto, la entropía de Shannon, 𝑆 es: 𝑆=−∫𝜌(𝑞)ln𝜌𝑑𝑞, (2.24) donde 𝑞 representa a las coordenadas generalizadas del problema, en nuestro caso, las posiciones 𝑥 e 𝑦. La densidad de probabilidad está dada por 𝜌=𝜌𝑛=|𝜓𝑛(𝑥,𝑦)|2. Por tanto, 𝑆=∬|𝜓𝑛(𝑥,𝑦)|2ln|𝜓𝑛(𝑥,𝑦)|2 𝑑𝑥𝑑𝑦 . (2.25) Sin embargo, Rényi, razonó que había una suposición implícita en la ecuación de Shannon, y es que usa la media lineal, ∑𝑝𝑘 𝐼𝑘 𝑀 𝑘=1 , la cual no es la única posibilidad. En la teoría, la información debería cumplir una propiedad matemática que implica que existe una función 𝑔 invertible tal que: 𝐼(𝑃)=𝑔−1(∑𝑝𝑘 𝑔(𝐼𝑘) 𝑀 𝑘=1 ). (2.26) Lo cual tiene dos soluciones: 𝑔(𝑥)=𝑐 𝑥, (2.27) 𝑔(𝑥)=𝑐−2(1−𝛼)𝑥. (2.28) La primera, la ecuación (2.27), lleva a la forma lineal de Shannon mientras que la segunda, (2.28) nos conduce a: 𝐼𝛼(𝑃)=1 1−𝛼ln(∑𝑝𝑘𝛼 𝑀 𝑘=1 ), (2.29) 15 para 𝛼 positivo y distinto de la unidad. Haciendo el mismo paso a la integral, obtenemos la expresión para la entropía de Rényi [15]: 𝑅𝛼=1 1−𝛼ln∫𝜌𝛼(𝑞) 𝑑𝑞 ∀𝛼∈[0,1)∪(1,∞). (2.30) En nuestro problema, esto se traduce a: 𝑅𝑛𝛼=1 1−𝛼ln∬|𝜓𝑛(𝑥,𝑦)|2𝛼 𝑑𝑥𝑑𝑦 ∀𝛼∈[0,1)∪(1,∞). (2.31) Para resolver la integral, debemos calcular la representación en el espacio de configuración de los autovectores que vienen de resolver la ecuación de Schrödinger. Para pasar de un vector a una función dependiente de las coordenadas espaciales, vamos a escribir esta 𝜓 como combinación lineal de autofunciones del oscilador armónico. En el apéndice B detallamos cómo podemos hacer esto con los cálculos y explicaciones pertinentes. En resumidas cuentas, nuestra matriz estaba escrita en la base del oscilador armónico. Por tanto, la matriz de cambio de base que nos lleva de la base del oscilador armónico a la base propia, (la que diagonaliza al hamiltoniano), será la de autovectores puestos por columnas. Es la matriz que hemos calculado con MATLAB. Por tanto, las autofunciones del hamiltoniano son una combinación lineal de autofunciones del oscilador armónico (que sabemos cómo son en el espacio de configuración) y los coeficientes de esta combinación lineal son las componentes de los autovectores. Es decir, 𝜓0(𝑥,𝑦)=∑∑𝐶𝑛𝑎𝑛𝑏 𝒩𝑛𝑎 𝐻𝑛𝑎(𝑥)exp(−𝑥2 2)𝒩𝑛𝑏 𝐻𝑛𝑏(𝑦)exp(−𝑦2 2) 𝑛𝑏 𝑛𝑎. (2.32) Solo quedaría sustituir esta expresión en la entropía de Rényi para calcularla y resolver la integral aplicando los métodos correspondientes en MATLAB. 2.4 Teoría de Campo Medio En teoría de campo medio (TCM) se estudia el comportamiento de modelos estocásticos de un número elevado de grados de libertad. Para ello, resuelve un modelo más simple que el original en el que se reducen los grados de libertad. En la TCM, el efecto de cada uno de los elementos de un sistema es sustituido por la media de sus efectos. Pasaríamos de un problema de muchos cuerpos a un problema reducido, normalmente de 1 cuerpo o 2. 16 Modelar nuestro problema con la teoría de campo medio nos permite, reescribiendo el hamiltoniano, aproximar nuestro problema por dos osciladores armónicos. En el apéndice C está realizado el cálculo pertinente con más detalle. La idea fundamental es reescribir el hamiltoniano de forma que podamos tener términos del estilo del oscilador armónico con resultados aproximables para 𝑁 grande. En primer lugar, definimos los siguientes operadores: 𝐾+≡12 𝑎†𝑎†, (2.33) 𝐾−≡12 𝑎𝑎, (2.34) 𝐾0≡12 (𝑎†𝑎+12). (2.35) Introducimos también las relaciones de Holstein-Primakoff. Estas son: 𝐾+=𝐶†(12+𝐶†𝐶)1/2, (2.36) 𝐾−=(12+𝐶†𝐶)1/2𝐶, (2.37) 𝐾0=𝐶†𝐶+14. (2.38) Con estos cambios, y las reglas de conmutación demostradas en el apéndice C, obtenemos el hamiltoniano escrito de la siguiente forma: 𝐻=𝜔0 𝐶†𝐶+ 𝜔 𝑏†𝑏+√2 𝑁 𝜆 ((12+𝐶†𝐶)1/2𝑏† 𝐶+𝑏 𝐶†(12+𝐶†𝐶)1/2). (2.39) Como hacíamos en el caso del oscilador armónico, vamos a pasar al espacio de posiciones y momentos (adimensionales). Sean (𝑥,𝑝) posiciones y momentos de átomos libres y (𝑦,𝑞) los correspondientes a moléculas diatómicas: 𝐶=√𝑁 𝑥+𝑖𝑝 √2, 𝐶†=√𝑁 𝑥−𝑖𝑝 √2, 𝑏=√𝑁 𝑦+𝑖𝑞 √2, 𝑏†=√𝑁 𝑦−𝑖𝑞 √2. (2.40) Por tanto, el hamiltoniano quedaría como: 𝐻=𝐻0+𝐻𝑖𝑛𝑡 (1)+𝐻𝑖𝑛𝑡 (2). (2.41) 17 Donde estos operadores son: 𝐻0=𝜔0𝑁2 (𝑥2+𝑝2)+𝜔 𝑁2 (𝑦2+𝑞2), (2.42) 𝐻𝑖𝑛𝑡 (1)=√𝑁2 𝜆 √12+𝑁2 (𝑥2+𝑝2) (𝑥𝑦+𝑝𝑞+𝑖𝑝𝑦−𝑖𝑞𝑥), (2.43) 𝐻𝑖𝑛𝑡 (2)=√𝑁2 𝜆 (𝑥𝑦+𝑝𝑞−𝑖𝑝𝑦+𝑖𝑞𝑥)√12+𝑁2 (𝑥2+𝑝2). (2.44) Ahora bien. Nuestra idea es que la descripción cuántica tiende a la clásica para un número elevado de partículas. Es fácil ver en este supuesto que el término en (2.43) y (2.44) de √12+𝑁 2 (𝑥2+𝑝2)≈√𝑁 2 (𝑥2+𝑝2). Siendo esto igual para un número infinito de partículas. Con esta aproximación, podemos dividir por el número de partículas. Sean las cantidades por partículas las dadas por (2.45): ℎ=𝐻 𝑁, ℎ0=𝐻0 𝑁. (2.45) Podemos escribir la ecuación (2.41) dividida por el número de partículas como: ℎ=ℎ0+𝜆2(√𝑥2+𝑝2 (𝑥𝑦+𝑝𝑞+𝑖𝑝𝑦−𝑖𝑞𝑥)+ (𝑥𝑦+𝑝𝑞−𝑖𝑝𝑦+𝑖𝑞𝑥)√𝑥2+𝑝2). (2.46) Si bien cuánticamente no podemos conseguir energía igual a cero, podemos ver qué ocurre cuando no hay energía cinética y todo es energía potencial, es decir, sea 𝑉(𝑥,𝑦)= lim 𝑝,𝑞→0ℎ(𝑝,𝑞) , 𝑉(𝑥,𝑦)=𝜔0 2 𝑥2+ 𝜔2 𝑦2+𝜆 𝑥2𝑦. (2.47) Además, el operador que nos cuenta partículas, 𝑎†𝑎+2𝑏†𝑏 escrito con las mismas relaciones y dividido por el número de partículas totales 𝑁 da como resultado: 𝑥2+𝑝2+ 𝑦2+𝑞2=1. (2.48) Esta es una condición de conservación de la energía asociada a la conservación de las partículas. En nuestro caso de límite a energía potencial, 𝑥2+𝑦2=1. Es decir, podemos hacer un cambio de variables a polares. 18 𝑉(𝜃)=𝜔2+𝜔0−𝜔 2 cos2𝜃+𝜆 cos2𝜃sen𝜃. (2.49) Si derivamos respecto a 𝜃 y calculamos los mínimos de energía, obtenemos que según los valores de 𝜆, habrá un valor mínimo posible distinto. Como está demostrado en el apéndice D, el 𝜆 a partir del cual podemos encontrar un valor de mínimos de energía distinto será el valor crítico: 𝜆𝑐=𝜔0−𝜔 2. (2.50) Por tanto, la energía mínima en función de 𝜆 es: 𝑉𝑐(𝜆)=𝜔2+(𝜆𝑐+𝜆sen𝜃𝑐)cos2𝜃𝑐. (2.51) Donde, 𝜃𝑐(𝜆)= { (2𝑘+1)𝜋2 ∀𝑘∈ℤ 𝑝𝑎𝑟𝑎 𝜆≤𝜆𝑐 arcsen(−∆𝜔±√(∆𝜔)2+12𝜆2 6𝜆 ) 𝑝𝑎𝑟𝑎 𝜆>𝜆𝑐 (2.52) La fracción de partículas monoatómicas es 𝑥=cos𝜃. En el estado fundamental será 𝑥= cos𝜃𝑐. Es decir, 𝑥𝑐(𝜆)={0 𝑝𝑎𝑟𝑎 𝜆≤𝜆𝑐 cos(𝑎𝑟𝑐𝑠𝑒𝑛(−∆𝜔±√(∆𝜔)2+12𝜆2 6𝜆 )) 𝑝𝑎𝑟𝑎 𝜆>𝜆𝑐 (2.53) Podemos saber qué transición de fase tendríamos según la teoría de campo medio. Para ello, vamos a usar 𝜆𝑐=0.5. Para distinguir el orden de transición, lo que debemos hacer es calcular la primera y la segunda derivada de la energía y estudiar su continuidad. En la figura 2, representamos el valor de la primera derivada de la ecuación (2.51) frente a 𝜆. 19 Fig. 2. Primera derivada de la energía mínima frente a 𝜆 para 𝜆𝑐=0.5. No existe ninguna discontinuidad en la primera derivada. Por lo tanto, no debe ser una transición de primer orden. Pasamos a representar frente a 𝜆 la segunda derivada en la figura 3. Fig. 3. Segunda derivada de la energía mínima frente a 𝜆 𝑝𝑎𝑟𝑎 𝜆𝑐=0.5. Como aparece una discontinuidad en la segunda derivada, se trata de una transición de segundo orden. 20 3. Resultados Si este hamiltoniano puede modelar un sistema cuántico en el que hay una transición de fase, las magnitudes que hemos ido definiendo, deberían sufrir cambios drásticos a partir de ciertos valores del parámetro 𝜆. Sin embargo, 𝜆 no es lo único que afecta al valor de estas magnitudes. Tenemos otras variables como el número de partículas 𝑁, que era importante en la deducción de 𝜆𝑐, y el valor de 𝜔 y 𝜔0. Para analizar los resultados, vamos a fijar el valor de 𝜔 y 𝜔0, lo cual nos da un valor crítico del parámetro de control, 𝜆𝑐 si lo que acabamos de demostrar teóricamente es correcto. Estudiaremos el valor de las magnitudes definidas anteriormente en función de 𝜆 para distintos valores de 𝑁. Luego, para un valor “apropiado” de 𝑁 (veremos qué significa esto), probaremos a cambiar el valor crítico 𝜆𝑐 cambiando los valores de 𝜔 y 𝜔0. A partir de ahí concluiremos si los resultados son consistentes con la teoría. 3.1 Dependencia con 𝑁 En esta sección vamos a establecer 𝜔=1 y 𝜔0=2 y, por tanto, 𝜆𝑐=(𝜔0−𝜔)/2=0.5. En la figura 4, vamos a representar la energía del estado fundamental frente a 𝜆 para distintos 𝑁, desde un valor bajo a un valor alto. Fig. 4: Energía del estado fundamental 𝐸0 frente a 𝜆 para 𝜆𝑐=0.5. 21 Como cabía esperar, tenemos más energía cuantas más partículas hay. Además, efectivamente se observa que la energía sufre un cambio brusco a partir del valor que esperábamos de 𝜆𝑐. Mientras más partículas tenemos, más pronunciado es el cambio en la energía fundamental. En la figura 5, vamos a hacer lo mismo con el valor esperado del número de átomos libres. Es decir, vamos a representar 𝑁𝐴 en frente a 𝜆. Fig. 5: Valor medio de 𝑁𝐴 frente a 𝜆 para 𝜆𝑐=0.5 . Este valor es cero o prácticamente cero cuando el sistema se encuentra en una fase y su número se dispara hasta alcanzar el equilibrio con el número de moléculas diatómicas, a cuyo valor medio podemos llamarle (siendo coherentes con la notación) 𝑁𝐵 . Y es que el valor de 𝑁𝐴 tiende a ser dos tercios de 𝑁, mientras que 𝑁𝐵 tiende a ser un tercio del mismo número. Esto es visible en la figura 6 si representamos hasta 𝜆 mayores el valor de 𝑁𝐴 por partícula (fracción de átomos libres) para 𝑁 grande (𝑁=500). 22 Fig. 6: Valor medio de 𝑁𝐴/𝑁 frente a 𝜆, para 𝜆𝑐=0.5 y 𝑁=500. Como conclusión, habría el mismo número de núcleos sin enlazar, que enlazados. Una vez analizados los observables fundamentales, energía y valor medio de átomos libres 𝑁𝐴, vamos a analizar otra magnitud que esperamos que nos pueda marcar la presencia de una transición de fase. Vamos a continuar con el IPR, que definimos en la sección (2.3.3). Se recuerda al lector que se trata de un valor que pide la polarización de un estado, en este caso el fundamental, en la base. En la figura 7 se representa el IPR para el estado fundamental del sistema frente al parámetro de control, 𝜆. Fig. 7: IPR frente a 𝜆 para 𝜆𝑐=0.5. 23 Como preveíamos, el IPR comienza siendo cercano a la unidad, y crece drásticamente hasta saturar en un valor. Sin embargo, este no es el valor máximo que podría alcanzar. Esto es indicativo de que el sistema parte del equilibrio en el estado es el de más baja energía del oscilador armónico, y llega a otro estado de equilibrio que es una combinación de varios autoestados del hamiltoniano. A la vista de esto resultados, cuando 𝜆 es lo suficientemente pequeño, tenemos un condensado de bosones diatómicos. Por tanto, el número de átomos libres 𝑁𝐴 es prácticamente cero y el estado cuántico tiene pocas contribuciones de la base. Realmente sólo una significativa, la del estado |0 𝑁/2 ⟩ que es un estado compatible con 100% de diátomo. A partir del valor crítico de 𝜆, que en este caso es 0.5, la energía del cambio de átomo-diátomo se hace relevante. Entonces, el número de átomos libres aumenta (creándose una mezcla átomodiátomo), la energía disminuye y contribuyen los estados compatibles con la mezcla. Entonces, se tiende a un valor constante de átomos y del IPR. Estamos entonces en el caso de coexistencia átomo-diátomo en equilibrio. Este estado de equilibrio en coexistencia será combinación de un número de autoestados del hamiltoniano en el que podamos encontrar valores los números de partículas que corresponden. Es por esto que el IPR no llega a su valor máximo, pero sí satura. 3.2 Dependencia de 𝜆𝑐 Una vez hemos comprobado la dependencia con 𝑁, vamos a hacer uso de una magnitud más significativa del sistema, la energía del estado fundamental por partícula. Veamos si el valor de 𝜆𝑐 es consistente con lo predicho en la teoría. Vamos a darle varios valores a 𝜔0 y 𝜔. Siempre respetando que 𝜔0>𝜔, ya que es la forma en la que los dos niveles energéticos están definidos. Como el modelo que nos predice 𝜆𝑐 funciona mejor para un número grande de partículas, vamos a representar las magnitudes anteriores para distintos valores de 𝜔0− 𝜔, con 𝑁 constante y grande, 700 partículas, por ejemplo. Usamos unos valores recogidos en la Tabla 1. 30 Fig. 14: Desviación de 𝑁𝐴/𝑁 frente a 𝑁 discreto para 𝜆𝑐=0.5. Obtenemos en la figura 14 exactamente la misma dependencia con 𝑁 que en el caso de la energía. Como conclusión, hemos probado que el modelo del campo medio funciona razonablemente bien sólo cuando tenemos un número muy elevado de partículas. Es decir, si 𝑁 es lo suficientemente grande, con muy buena aproximación, podemos predecir el comportamiento del sistema sin resolver el problema exacto. Una vez hemos confirmado la validez de la teoría de campo medio para un número alto de partículas, podemos darles veracidad a las predicciones realizadas por la teoría de campo medio para este problema. Una de ellas era acerca del orden de la transición: no sólo tendríamos una transición de fase (cosa que de momento es indiscutible) sino que, además, sería de segundo orden. Sin embargo, aunque la TCM lo sugiere, no hemos demostrado en el problema exacto que la transición de fase sea compatible con una de segundo orden. Esto sería posible mediante el concepto de entropía desarrollado en la sección 2.3.4. Si la teoría de la información en física cuántica ha capturado bien el concepto de entropía, deberíamos ser capaces de determinar el orden de la transición estudiando su continuidad. 31 3.4 Entropía Para calcular la entropía sólo es necesario resolver las integrales (2.25) y (2.31). Esto, una vez más, lo debemos programar con MATLAB. La entropía ha sido calculada para 𝑁=100, puesto que una menor dimensión agiliza la integral a resolver. La entropía no es un concepto que hayamos predicho por la TCM, así que un número 𝑁 pequeño o grande es irrelevante para lo que tratamos de confirmar. Es más, lo que precisamente nos interesa es aquel comportamiento que difiere significativamente del modelo de campo medio. Lo que esperamos obtener si representamos una entropía frente a 𝜆, en un comportamiento que nos permita distinguir dos fases bien diferenciadas. Podemos comenzar con la primera entropía que hemos definido, la entropía de Shannon, representada en la figura 15 para 𝜆𝑐=0.5. Fig. 15: Entropía de Shannon frente a 𝜆 para 𝜆𝑐=0.5 y 𝑁=100. 32 En esta, vemos dos fases diferenciadas caracterizadas por dos valores de la entropía. Hay un máximo en la transición de fase debido al enorme cambio de probabilidades, al cambiar los autoestados relevantes del sistema. La entropía debería ser menor en la fase de condensado de moléculas diatómicas, por la localización del estado cuántico y mayor en la fase de coexistencia. Efectivamente, concuerda con el resultado obtenido. Sabemos que la entropía de Rényi, 𝑅𝛼 en (2.31) tiende a la de Shannon cuando 𝛼→1. Lo que tenemos que determinar es qué sucede para 𝛼>1 y para 𝛼<1. En primer lugar, podemos escoger un valor de 𝛼 entre 0 y 1, por ejemplo 𝛼=1/2. Si calculamos el valor de la entropía de Rényi para este valor, 𝑅1/2 en función de 𝜆, obtenemos la figura 16. Fig. 16 Entropía de Rényi para 𝛼=1/2, 𝑁=100 y 𝜆𝑐=0.5. Es una figura con cierto parecido a la de Shannon, con dos fases diferenciadas por diferentes valores de la entropía, menor y mayor para las fases mencionadas con anterioridad. Atendiendo a la ecuación (2.30) o a la (2.31), para 𝛼→0 la entropía pesa a todos los eventos de forma más igualitaria independientemente de su densidad de probabilidad, dada por |𝜓𝑛(𝑥,𝑦)|2. Por tanto, obtenemos un valor parecido al que obtenemos con la expresión de Shannon, pero la diferencia 33 entre las entropías a izquierda y derecha es menor. Por tanto, lo contrario debería ocurrir si aumentamos 𝛼. En la figura 17 representamos la entropía de Rényi para 𝛼=2. Fig. 17 Entropía de Rényi para 𝛼=2, 𝑁=100 y 𝜆𝑐=0.5. En este caso, podemos apreciar incluso más claramente un comportamiento de transición de fase. Además, con cambio bastante brusco, considerando el bajo valor de 𝑁. No sólo se cumple lo anterior, sino que la diferencia entre los valores de la entropía de cada fase es mucho mayor. Este comportamiento de la entropía es análogo al que se da en una transición de segundo orden al que nos apuntaba el modelo de campo medio. La clasificación de esta transición de fase está dada ahora tanto por la resolución exacta del problema como por su aproximación en el modelo de campo medio. 34 4. Conclusiones Hemos modelado un sistema de dos niveles mediante un hamiltoniano en segunda cuantización y hemos resuelto su problema de autovalores y autovectores. A partir de estos, hemos calculado magnitudes relevantes, definidas en la sección 2.3, que pudiesen poner de manifiesto que ocurre una transición de fase cuántica para el nivel fundamental. En las secciones 3.1, 3.2 y 3.4 hemos comprobado que para distintos valores de 𝑁 y de las energías de cada nivel, siempre existe para el estado fundamental un comportamiento diferenciado en dos fases: una fase con moléculas diatómicas y otra con coexistencia átomo-diátomo. En estas fases tenemos propiedades diferenciadas: la energía, el valor medio de las partículas en cada nivel energético, el IPR y la entropía. En la sección 3.3 hemos comprobado que este comportamiento tiende a lo predicho por el modelo de campo medio para un número lo suficientemente elevado de partículas. Por tanto, no sería necesario resolver el problema exacto bajo estas condiciones. Por tanto, hemos comprobado que efectivamente un sistema descrito por un hamiltoniano de estas características podría en su estado fundamental estar en dos fases cuánticas diferenciadas y que realizaría una transición de una a otra cambiando bruscamente sus propiedades. Es decir, que existiría una transición de fase cuántica tal y como hemos postulado en la introducción, para un sistema de dos niveles. También hemos comprobado por el modelo exacto y el modelo de campo medio que la transición de fase sería de segundo orden o continua. Cabe destacar que esto no prueba la existencia de transiciones de fase cuánticas. Como comentábamos en la sección 1.3.1, un hamiltoniano no es más que un modelo, como todo tratamiento en física. Es tarea de la disciplina experimental comprobar para qué sistema real, si alguno, se tiene un comportamiento que se corresponda con muy buena aproximación con un hamiltoniano de estas características. Es decir, lo que hemos demostrado no es que existan las transiciones de fase cuánticas, sino que para un sistema cuántico que se describe como en nuestro modelo, existiría una transición de fase cuántica. No obstante, dado que modelos tan idílicos y simples como el oscilador armónico, la partícula libre o el gas ideal funcionan razonablemente bien en determinadas circunstancias, sería lógico que se puedan dar nuestras hipótesis en un sistema real. Al fin y al cabo, sólo estamos considerando que el hamiltoniano de un sistema de dos niveles está dado por la energía debida a las partículas 35 que hay en cada nivel y por una interacción simple de paso de un nivel a otro. La única hipótesis adicional que puede chocar con un problema real es que nos hemos restringido a un sistema con un número par de partículas. Si bien es cierto que esto puede suponer diferencias en los resultados para un valor de 𝑁 muy pequeño, las propiedades de un sistema de muchas partículas no van a cambiar de forma significativa por añadir o eliminar una. Es decir, 𝑁±1≈𝑁. De hecho, el modelo de campo medio (que considera 𝑁→∞) no tiene este problema. Además, hemos demostrado que este modelo de TCM funciona de forma suficientemente parecida a resolver el problema exacto, que tiene su restricción de 𝑁 par. Por lo tanto, se confirma que 𝑁 par es sólo una condición para hacer más ventajoso el cálculo, sin que esto altere significativamente las propiedades del sistema [5]. Por todo esto, si hay un sistema real descrito por este hamiltoniano, se cumplirían las definiciones de fase y transición de fase en la sección 1.3.2. Es un hecho que en la naturaleza se observan cambios bruscos en las propiedades de los sistemas a partir de valores críticos de parámetros que gobiernan sus comportamientos. Y ya sea en física clásica o en física cuántica, estos cambios bruscos se traducen en comportamientos en ciertas magnitudes que se pueden modelar como discontinuos. En todos los fenómenos conocidos, dichas magnitudes estarían relacionadas con conceptos energéticos. Por lo tanto, damos por satisfactoria la redefinición de transición de fase y su descripción cuántica. 36 5. Referencias [1] Cracknell, A. P., Lorenc, J., & Przystawa, J. A. Landau’s theory of second-order phase transitions and its application to ferromagnetism. Journal of Physics C: Solid State Physics. 9(9), 1731–1758. https://doi.org/10.1088/0022-3719/9/9/015, (1976). [2] Yeomans, J. M., & Rudnick, J. Statistical Mechanics of Phase Transitions. Physics Today. 46(7) 80 https://doi.org/10.1063/1.2808979, (1993). [3] Jaeger, Gregg. "The Ehrenfest Classification of Phase Transitions: Introduction and Evolution". Archive for History of Exact Sciences. 53 (1): 51–81, (1998). [4] Fernández Tejero, C & Baus, M. Física Estadística del Equilibrio: fases de la materia. Aula documental de investigación, (2000). [5] Reif, F., & Scott, H. L, Fundamentals of Statistical and Thermal Physics. American Journal of Physics. https://doi.org/10.1119/1.19073, (1998). [6] Zamora, M Termo I: un estudio de los sistemas termodinámicos. Publicaciones de la Universidad de Sevilla, (1998) [7] Pokrovsky, V. L, Landau and modern physics. Physics-Uspekhi, 52(11), 1169–1176. https://doi.org/10.3367/ufne.0179.200911j.1237 (2009). [8] McNaught, A. D. & Wilkinson, A. Compendium of chemical terminologythe “Gold book.” In Blackwell Scientific. https://doi.org/10.1351/goldbook.I03352, (1997). [9] Sachdev, S. Quantum phase transitions. In Physics World (Vol. 12, Issue 4). https://doi.org/10.1088/2058-7058/12/4/23, (1999). [10] Blundell, S. J., & Blundell, K. M. Concepts in Thermal Physics. In Concepts in Thermal Physics. Oxford Scholarship Online. https://doi.org/10.1093/acprof:oso/9780199562091.001.0001, (2010). [11] Carr, L. Understanding quantum phase transitions. In Understanding Quantum Phase Transitions. CRC Press. https://doi.org/10.1201/b10273, (2010). [12] Dubbers, D., & Stöckmann, H.-J. Quantum Physics: The Bottom-Up Approach. In From the Simple Two-Level System to Irreducible Representations. Springer. (2013) 37 [13] Juan Gamito, Trabajo Fin de Máster "Análisis de un quantum quench en el modelo de Lipkin anarmónico", Universidad de Sevilla, 2019 y referencias en él [14] Sarah S Chehade and Anna Vershynina (2019), Scholarpedia, 14(2):53131. doi:10.4249/scholarpedia.53131 [15] Romera, E., del Real, R., Calixto, M., Nagy, S., & Nagy, Á. (2013). Rényi entropy of the U(3) vibron model. Journal of Mathematical Chemistry. https://doi.org/10.1007/s10910-012-0106-7 38 6. Apéndices A continuación, presentamos aquellos cálculos que se han omitido por no ser relevantes para el informe científico, pero que complementan la información que se da en el mismo. Así mismo, presentamos el código realizado en MATLAB con el que hemos conseguido toda la información relevante. 6.1 Apéndice A: resolución matricial del problema En esta sección vamos a ilustrar cómo se resolvería el problema exacto. Comenzamos representando las matrices: Al ser base propia de 𝐻𝑎 y 𝐻𝑏, estos son diagonales y en la diagonal aparece el número de bosones monoatómicos y diatómicos respectivamente que hay en cada estado. 𝐻𝑎𝑏 tiene ceros en la diagonal y términos distintos de cero con coordenadas 𝐻𝑖 𝑖±1 sin exceder los límites de la dimensión de la matriz. Para un sistema con 6 partículas la base viene dada por los siguientes vectores: {|0 3⟩,|2 2⟩,|4 1⟩,|6 0⟩}. 𝐻𝑎↔12𝜔0 ( ⟨0 3|𝑎†𝑎|0 3⟩ ⟨0 3|𝑎†𝑎|2 2⟩ ⟨2 2|𝑎†𝑎|0 3⟩ ⟨2 2|𝑎†𝑎|2 2⟩⟨0 3|𝑎†𝑎|4 1⟩ ⟨0 3|𝑎†𝑎|6 0⟩ ⟨2 2|𝑎†𝑎|4 1⟩ ⟨2 2|𝑎†𝑎|6 0⟩ ⟨4 1|𝑎†𝑎|0 3⟩ ⟨4 1|𝑎†𝑎|2 2⟩ ⟨6 0|𝑎†𝑎|0 3⟩ ⟨6 0|𝑎†𝑎|2 2⟩⟨4 1|𝑎†𝑎|4 1⟩ ⟨4 1|𝑎†𝑎|6 0⟩ ⟨6 0|𝑎†𝑎|4 1⟩ ⟨6 0|𝑎†𝑎|6 0⟩ ) . (6.𝐴.1) Como hemos introducido, el operador 𝑎†𝑎 cuenta el número de partículas 𝑛𝑎, y como los autovectores de un oscilador armónico son ortonormales: 𝐻𝑎↔12 𝜔0(0 0 0 2 0 0 0 0 0 0 0 0 4 0 0 6). (6.𝐴.2) De forma análoga: 𝐻𝑏↔𝜔 ( ⟨0 3|𝑏†𝑏|0 3⟩ ⟨0 3|𝑏†𝑏|2 2⟩ ⟨2 2|𝑏†𝑏|0 3⟩ ⟨2 2|𝑏†𝑏|2 2⟩ ⟨0 3|𝑏†𝑏|4 1⟩ ⟨0 3|𝑏†𝑏|6 0⟩ ⟨2 2|𝑏†𝑏|4 1⟩ ⟨2 2|𝑏†𝑏|6 0⟩ ⟨4 1|𝑏†𝑏|0 3⟩ ⟨4 1|𝑏†𝑏|2 2⟩ ⟨6 0|𝑏†𝑏|0 3⟩ ⟨6 0|𝑏†𝑏|2 2⟩ ⟨4 1|𝑏†𝑏|4 1⟩ ⟨4 1|𝑏†𝑏|6 0⟩ ⟨6 0|𝑏†𝑏|4 1⟩ ⟨6 0|𝑏†𝑏|6 0⟩ ) , 39 𝐻𝑏↔𝜔(3 0 0 2 0 0 0 0 0 0 0 0 1 0 0 0) (6.𝐴.3) Una vez más, aplicamos los operadores para 𝐻𝑎𝑏: 𝜆 𝐻𝑎𝑏↔𝜆 √2·6 ( 0√1 √2 √3 √1 √2 √3 0 0 0 √4 √2 √3 0 0 √4 √2 √3 0 0 0√1 √5 √6 √1 √5 √6 0 ) . 𝜆 𝐻𝑎𝑏↔𝜆 √12 ( 0√6 √6 0 0 0 2√6 0 0 2√6 0 0 0√30 √30 0 ) . (6.𝐴.4) Evidentemente, obtenemos una matriz 𝐻 hermítica. ________________________________________________ Usamos el símbolo ↔ para denotar que un operador “se corresponde” con una matriz. Lo cual es conceptualmente distinto a que el operador sea una matriz. Un operador existe como elemento abstracto en un cierto espacio vectorial y una matriz es una representación en una base concreta de dicho operador. 6.2 Apéndice B: autofunciones en el espacio de configuración En esta sección vamos a ver cómo se representan en el espacio de configuración los autovectores obtenidos de resolver la ecuación de Schrödinger (2.12). Recordamos que para construir el hamiltoniano en forma matricial habíamos escrito sus componentes en (2.14) 𝐻𝑖𝑗=⟨𝑖|𝐻|𝑗⟩. (2.14) Donde los vectores |𝑖⟩ y |𝑗⟩ son los autovectores de un problema del oscilador armónico. Al diagonalizar el hamiltoniano, lo que estamos haciendo es “girar” los vectores con una matriz 𝑉 definida como: 𝑉|𝑖⟩≡|𝑖′⟩. (6.𝐵.1) Es relevante comentar que operador se define por cómo actúa sobre una base completa. 46 𝜆′≥12 24=1 2 , (6.𝐷.10) es decir, el valor mínimo de 𝜆′ a partir del cual puede encontrarse un mínimo distinto en la energía del estado fundamental será 0.5. Este cambio en el comportamiento de la energía del estado fundamental nos marcaría un cambio de fase. Si deshacemos el cambio, 𝜆𝑐=∆𝜔 2=𝜔0−𝜔 2. (6.𝐷.11) 6.5 Apéndice E: Códigos de MATLAB En esta sección se escribirán los diversos códigos, complementarios entre sí, que hemos escrito para conseguir lo expuesto en la sección 3 (resultados). El código está comentado adecuadamente para ir siguiendo qué es lo que se hace en cada momento y por qué. 6.5.1 Programa principal %% Modelo Exacto % Tenemos el problema exacto, que vamos a resolver paso a paso % Existen varias formas de configurar el programa en función de lo que queramos concretamente % Esta es la versión más general y adaptable a cada necesidad. %% Parametros del problema % Si queremos calcular la entropía de Rényi Entropia_value=false; % Si queremos magnitudes por partícula Porparticula_value=true; % Podemos cambiar colores fácilmente con la siguiente variable stringer=['b' 'r' 'g' 'm' 'k']; % Debe tener al menos tanta longitud como N % Si queremos calcular las magnitudes´E0 y NA en el punto crítico critico_value=false; % Si queremos una forma aproximada de la ec de autovalores. No tiene máyor % interés ajustepolinomico_value=false; % En las primeras versiones del programa, se introducía "a mano" el valor % de N. La idea es no tener que ir a las primeras lineas continuamente para % ir probando distintos N. % Es posible sustituirlo simplemente por cualquier valor % Elegimos N. Puede ser un vector si queremos el mismo cálculo para varios % N. No es aconsejable hacerlo para la entropía de Rényi. 47 N=input('Elige número de partículas N= '); Nlength=length(N); for kn=1:Nlength % Admitimos que N debe ser real, par y positivo. % Introducimos condiciones para arreglar los posibles valores ilógicos que % el ejecutante del programa pueda introducir if imag(N(kn))~=0 % Correción para números con parte imaginaria warning('N debe ser real') N(kn)=real(N(kn)); fprintf('N(kn)= Parte real de N: %0.3f \n',N) end if round(N(kn))~=N(kn) % Correción para números con parte decimal warning('N debe ser entero') N(kn)=round(N(kn)); fprintf('N= parte entera de N: %0.0f \n',N(kn)) end if N(kn)<0 warning('N debe ser positivo') N(kn)=-N(kn); fprintf('N= valor absoluto de N: %0.0f \n',N(kn)) end if mod(N(kn),2)~=0 % Corrección para números impares warning('N=N+1 (debe ser par)') N(kn)=N(kn)+1; fprintf('Valor de N= %0.0f \n',N(kn)) end if N(kn)==0 N(kn)=2*round(200*rand); % Número aleatorio de partículas fprintf('Valor aleatorio de N= %0.0f \n',N(kn)) end end for kn=1:Nlength %% Contrucción de la base % La base está dada por v=[C Na Nb] que son: % C = constante que multiplica al ket. Si C=1, el estado está normalizado % Na = número de partículas atómicas % Nb = número de partículas diatómicas % Definimos los estados en los que puede estar el sistema, desde el nivel % fundamental [1 0 N], primero excitado: [1 1 N-2], etc. Nb=linspace(0,N(kn)/2,N(kn)/2+1); D=length(Nb); % Dimensión de la base Aux=zeros(1,D); 48 Na=N(kn)-Nb*2; % Vamos a almacenar cada vector de la base en una fila distinta % de una matriz. La dimensión de la matriz será D x 3. base=ones(D,3); for i=1:D base(D+1-i,2)=Na(i); base(D+1-i,3)=Nb(i); end % La base está ordenada. El autovalor con menor energía debería ser el % primero que aparezca en la matriz diagonal. Si no nos damos cuenta de % esto, el programa se encarga de encontrar este valor más adelante %% hamiltoniano adimensional actuando sobre vectores de la base % Sin considerar las constantes delante de cada operador. % Tenemos 4 términos. Son, por orden, los del anunciado. % Definimos H | base > = E Ea=zeros(D,3); Eb=zeros(D,3); Eab1=zeros(D,3); Eab2=zeros(D,3); % Las siguientes M-funciones están en otros ficheros. Son los operadores a, % a daga, b, b daga. for i=1:D Ea(i,:)=adaga(a(base(i,:))); Eb(i,:)=bdaga(b(base(i,:))); Eab2(i,:)=(bdaga(a(a(base(i,:))))); Eab1(i,:)=b(adaga(adaga(base(i,:)))); end % Aquí tenemos cómo actúa cada elemento adimensional del hamiltoniano sobre % cada elemento de la base, almacenado en matrices %% Valor medio del hamiltoniano (contribuciones adimensionales) % < base | hamiltoniano adimensional | base > % < base | Matrices Ea, Eb... > % < C Na Nb | C' Na' Nb'> (Donde C=1) % Esto será cero, {zeros(D,D)} a no ser que Na==Na' y Nb==Nb'. % Por tanto: < 1 Na Nb | C' Na' Nb'> = C' o 0. Ha=zeros(D,D); Hb=zeros(D,D); Hab1=zeros(D,D); Hab2=zeros(D,D); % Cada if le pregunta a un hamiltoniano adimensional sobre lo anterior for i=1:D 49 for j=1:D if base(i,2)==Ea(j,2) && base(i,3)==Ea(j,3) Ha(i,j)=Ea(j,1); end if base(i,2)==Eb(j,2) && base(i,3)==Eb(j,3) Hb(i,j)=Eb(j,1); end if base(i,2)==Eab1(j,2) && base(i,3)==Eab1(j,3) Hab1(i,j)=Eab1(j,1); end if base(i,2)==Eab2(j,2) && base(i,3)==Eab2(j,3) Hab2(i,j)=Eab2(j,1); end end end % Ahora solo tenemos que sumar las contribuciones con sus respectivas % constantes, omega, omega0 y lambda %% Cálculo de magnitudes omega=1; omega0=2; % Podríamos cambiar estos parámetros omega0red=omega0/2; critico=(omega0-omega)/2; % punto crítico. Lo calculamos para graficarlo también Cn=1/sqrt(2*N(kn)); % A continuación, calculamos la energía del estado fundamental para % distintos valores de lambda que contienen al punto de interés. if Entropia_value==true lambda=linspace(0,2*critico,80); else lambda=linspace(0,3*critico,1000); % lambda conteniendo el punto de interés end % el cálculo para rényi es mucho menor para poder aligerar el programa. L=length(lambda); E0=zeros(1,L); % Aquí almacenamos la energía del estado fundamental F_onda0=zeros(D,L); %Aquí almacenamos la F. de onda del estado fundamental ValorMedioHa=zeros(1,L); % Almacenamos aquí el valor medio IPR=zeros(1,L); % Almacenamos aquí los valores del IPR en función de lambda % en un único bucle calculamos IPR,E0, NA y funciones de onda 50 for i=1:L H=omega0red*Ha+omega*Hb+(Hab1+Hab2)*lambda(i)*Cn; % Matriz Hamiloniano [V,Diagonal]=eig(H); % Resolución problema autovalores y autovectores Autovalores=diag(Diagonal)'; E0(i)=Autovalores(1); % La base está en orden para ahorrarnos el comando min(Autovalores) F_onda0(:,i)=V(:,1); % Calculamos el valor medio ValorMedioHa(i)=F_onda0(:,i)'*Ha*F_onda0(:,i); % IPR IPR(i)=1/sum(V(:,1).^4); end if Porparticula_value==true E0=E0/N(kn); ValorMedioHa=ValorMedioHa/N(kn); end % V = matriz autovectores % Diagonal = hamiltoniano diagonalizado % Autovalores = vector fila formado por autovalores % E0=energía del estado fundamental para cada lambda % Representamos las energías del estado fundamental frente a lambda figure(1) plot(lambda,E0,stringer(kn)) % axis([0 3*critico min(E0) max(E0)+5]) axis([0 1.5 0 0.55]) title('Energía estado fundamental frente a {\lambda}','Fontsize',20) xlabel('{\lambda}','Fontsize',20) ylabel(['E_{' int2str(0) '}'], 'Fontsize',20); hold on % Para poner el * en el valor crítico. E0 para el lambda crítico: if critico_value==true H=omega0red*Ha+omega*Hb+(Hab1+Hab2)*critico*Cn; [V,Diagonal]=eig(H); Autovalores=diag(Diagonal)'; Ecritico=Autovalores(1); % Y lo representamos en rojo plot(critico,Ecritico,'*r') legend('Energía fundamental','Valor Crítico') end % Representamos ahora el IPR figure(2) plot(lambda,IPR,stringer(kn)) title('IPR frente a {\lambda} para N=500') xlabel('{\lambda}') ylabel('Coeficiente de participación inversa') %% NA 51 figure(3) % Lo representamos plot(lambda,ValorMedioHa,stringer(kn)) hold on subindex='A'; title('Valor medio de N_A frente a {\lambda}','Fontsize',20) xlabel('{\lambda}','Fontsize',20) ylabel('N_A','Fontsize',20) % Para poner el * en el valor crítico if critico_value==true F_onda0Critico=V(:,1); % En V está almacenado el valor crítico ValorMedioHaCritico=F_onda0Critico'*Ha*F_onda0Critico; plot(critico,ValorMedioHaCritico,'rx') legend('Valor medio de Ha', 'Critico') end %% Si quisiesemos la ecuación de autovalores dependiente de lambda if ajustepolinomico_value==true cE0=polyfit(lambda,E0,D); % Me da la E0 en función de lambda % % Comprobación gráfica E0pol=0; for i=1:length(cE0) E0pol=E0pol+cE0(i)*lambda.^(D-i+1); end figure(1) plot(lambda,E0pol,'r') end end %% Entropía de Renyi if Entropia_value==true alpha=1/2; C_Ralpha=1/(1-alpha); dosalpha=2*alpha; % los números cuánticos son: 52 na=base(:,2)'; nb=base(:,3)'; %% ENTROPIA DE RENYI % descomentar para Renyi: % limite=40; % pasos=0.1; % Shannon limite=13; % comentar para Shannon pasos=0.02; % igual x=-limite:pasos:limite; y=x; xlength=length(x); I=zeros(xlength,xlength); S=I; Ralpha=zeros(1,L); Shannon=Ralpha; n1=0:1:N(kn)/2; n2=N(kn)/2+2:2:N(kn); n=[n1 n2]; nlength=length(n); LordHermite=zeros(N(kn)+1,xlength); % Aquí calculamos los polinomios de hermite relevantes for i=1:nlength for j=1:xlength LordHermite(n(i)+1,j)=hermite2(n(i),x(j))*exp(- x(j)^2/2)/sqrt(2^(n(i)).*factorial(n(i))*sqrt(pi)); end end % Aquí calculamos las entropías de Shannon o Rényi para cada lambda for t=1:L % Aquí calculamos el integrando de las expresiones correspondientes for i=1:xlength for j=1:xlength IH=0; for k=1:D IH=IH+F_onda0(k,t)*LordHermite(na(k)+1,i)*LordHermite(nb(k)+1,j); end % I(i,j)=abs(IH)^dosalpha; % Rényi % Entropía de Shannon: 53 if IH~=0 IH=abs(IH); S(i,j)=IH*log(IH); end % Fin de entropía de Shannon end end % E integramos: % Ralpha(t)=C_Ralpha*log(trapz(y,trapz(x,I,2))); % Renyi Shannon(t)=-trapz(y,trapz(x,S,2)); % Shannon I=zeros(xlength,xlength); S=I; % Shannon end % Rényi % figure(4) % hold on % plot(lambda,Ralpha) % % xlabel('{\lambda}','Fontsize',20) % ylabel('R^{1/2}','Fontsize',18) % Shannon figure(5) plot(lambda,Shannon) xlabel('{\lambda}','Fontsize',20) ylabel('S','Fontsize',18) hold on end 6.5.2 Operadores de creación y destrucción En esta sección se exponen los 4 operadores relevantes. Estos son: Operador de creación de partícula libre function[B]=adaga(A) if A(1)==0 B=zeros(1,3); else B(1)=A(1)*sqrt(A(2)+1); B(2)=A(2)+1; B(3)=A(3); end end Operador de creación de partícula diatómica 54 function[B]=bdaga(A) if A(1)==0 B=zeros(1,3); else B(1)=A(1)*sqrt(A(3)+1); B(2)=A(2); B(3)=A(3)+1; end end Operador de destrucción de partícula libre: function[B]=a(A) if A(1)==0 B=zeros(1,3); else B(1)=A(1)*sqrt(A(2)); B(2)=A(2)-1; B(3)=A(3); end end Operador de destrucción de partícula diatómica: function[B]=b(A) if A(1)==0 B=zeros(1,3); else B(1)=A(1)*sqrt(A(3)); B(2)=A(2); B(3)=A(3)-1; end end 6.5.3 Polinomios de Hermite function [h] = hermite2 (n, x) % h = hermite(n, x) devuelve el polinomio de Hermite de grado n evaluado en x h = zeros(1,n+1); n_fact = factorial(n); m = 0:floor(n/2); h(2*m+1) = n_fact .* (-1).^m ./ (factorial(m) .* factorial(n-2.*m)) .* 2.^(n2.*m); if exist('x','var') h = polyval(h, x); end end 6.5.4 Modelo de campo medio % desviacion_value será true si queremos calcular la diferencia de lo % calculado anteriormente por el modelo de campo medio. desviacion_value=false; 55 %% Derivadas de la energía V=@(lam) 0.5+(0.5+lam*(-1-sqrt(1+12*lam^2))/(6*lam))*cos(asin((-1sqrt(1+12*lam^2))/(6*lam)))^2; Vtheta=@(theta,lam) omega0/2*cos(theta)^2+omega/2*sin(theta)^2+lam*cos(theta)^2*sin(theta); dV=@(lam) (2*lam*(((12*lam^2 + 1)^(1/2) + 1)^2/(36*lam^2) - 1))/(12*lam^2 + 1)^(1/2) - (((12*lam^2 + 1)^(1/2) + 1)^2/(18*lam^3) - (2*((12*lam^2 + 1)^(1/2) + 1))/(3*lam*(12*lam^2 + 1)^(1/2)))*((12*lam^2 + 1)^(1/2)/6 - 1/3); dV2=@(lam) (2*(((12*lam^2 + 1)^(1/2) + 1)^2/(36*lam^2) - 1))/(12*lam^2 + 1)^(1/2) - ((12*lam^2 + 1)^(1/2)/6 - 1/3)*((8*((12*lam^2 + 1)^(1/2) + 1))/(12*lam^2 + 1)^(3/2) - ((12*lam^2 + 1)^(1/2) + 1)^2/(6*lam^4) - 8/(12*lam^2 + 1) + (2*((12*lam^2 + 1)^(1/2) + 1))/(lam^2*(12*lam^2 + 1)^(1/2))) - (4*lam*(((12*lam^2 + 1)^(1/2) + 1)^2/(18*lam^3) - (2*((12*lam^2 + 1)^(1/2) + 1))/(3*lam*(12*lam^2 + 1)^(1/2))))/(12*lam^2 + 1)^(1/2) - (24*lam^2*(((12*lam^2 + 1)^(1/2) + 1)^2/(36*lam^2) - 1))/(12*lam^2 + 1)^(3/2); % dV=@(theta,lam) lam*cos(theta)^3 - cos(theta)*sin(theta) - 2*lam*cos(theta)*sin(theta)^2; % dV2=@(theta,lam) 2*lam*sin(theta).^3 - cos(theta).^2 + sin(theta).^2 - 7*lam*cos(theta).^2*sin(theta); thetacritico=[asin(1) asin(-1) asin(1/3)]; E0teo=zeros(1,L); dVc=E0teo; dV2c=E0teo; Nateo=E0teo; thetalambda=asin((-ones(1,L)-sqrt(ones(1,L)+12*lambda.^2))./(6*lambda)); thetalambdamal=asin((-ones(1,L)+sqrt(ones(1,L)+12*lambda.^2))./(6*lambda)); for i=1:L if lambda(i)<0.5 E0teo(i)=0.5; % Nateo(i)=0; else dVc(i)=dV(lambda(i)); dV2c(i)=dV2(lambda(i)); E0teo(i)=V(lambda(i)); Nateo(i)=cos(asin((-1-sqrt(1+12*lambda(i)^2))/(6*lambda(i))))^2; end end figure(1) hold on plot(lambda,E0teo,'--g') legend('N=50','N=700','Campo Medio') figure(3) hold on plot(lambda,Nateo,'--g')