El método de mallado de Lagrange y aplicaciones en Física Cuántica

El método de mallado de Lagrange y aplicaciones en Física Cuántica Trabajo de Fin de Grado Doble Grado en Física y Matemáticas Ángel Lucena Pérez Tutorizado por D. José Miguel Arias Carrasco y D. José Antonio Lay Valera abstract 1 Abstract Solving the Schrödinger equation is essential for understanding quantum physical systems. With the exception of a few systems, most are described by Hamiltonians whose mathematical resolution does not provide analytical solutions. Therefore, the application of numerical methods is convenient for their study. This paper presents a method that allows us to solve the time-independent Schrödinger equation in matrix form. The formalism used is based on states located in configuration space (CLS). Its construction begins with a family of orthogonal polynomials. With these, we will calculate relevant physical quantities such as energy or wave functions. First, the results will be compared with potentials with analytical solution such as the Morse potential. Second, it will be applied to potentials of interest in Molecular Physics, such as the Kratzer-Fues, Deng-Fan, and Varshni potentials. Finally, the method will be extended to be applied to 3D systems and the ro-vibrational espectrum of 𝑂2 molecule will be studied. The following pages will describe the theoretical framework of the method and the obtained results with a code implemented in Matlab (which can be consulted in appendix A). Both the code used to implement the method and the visualization of the results were created entirely by the author. 2 Resumen La resolución de la ecuación de Schrödinger resulta fundamental para la comprensión de los sistemas físicos cuánticos. A excepción de algunos pocos sistemas, la mayoría se describen mediante hamiltonianos cuya resolución matemática no proporciona soluciones analíticas. Es por ello que resulta conveniente la aplicación de métodos numéricos para poder abordar el estudio de los mismos. En este trabajo se presenta un método que nos permite resolver de forma matricial la ecuación de Schrödinger independiente del tiempo. El formalismo usado se basa en estados localizados en el espacio de configuración, o CLS (siglas procedentes del inglés). Para su construcción, se partirá de una familia de polinomios ortogonales. Con ellos, calcularemos magnitudes relevantes físicas como la energía o las funciones de onda. En primer lugar, se compararán los resultados con potenciales de solución conocida como el de Morse. En segundo lugar, se aplicará a potenciales de interés en Física Molecular como el potencial de Kratzer-Fues, Deng-Fan y Varshni. Finalmente, el método se extenderá para su aplicación a sistemas en 3D y se estudiará el espectro vibro-rotacional de la molécula de 𝑂2. Las siguientes páginas describirán el marco teórico del método y los resultados obtenidos con un código implementado en Matlab (que puede consultarse en el apéndice A). Tanto el código utilizado para implementar el método como la visualización de los resultados fueron creados íntegramente por el autor. Índice general Abstract.............................................. 1 Resumen ............................................. 2 1. Nociones básicas de Física Cuántica. 5 1.1. Los postulados de la Física Cuántica. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2. La ecuación de Schrödinger independiente del tiempo. . . . . . . . . . . . . . . . . . 6 1.3. Propiedades de las funciones de onda. . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4. Introducción a la Física molecular. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4.1. Nociones básicas sobre la molécula. . . . . . . . . . . . . . . . . . . . . . . . 9 1.4.2. La aproximación de Born-Opphenheimer para moléculas diatómicas. . . . . 10 1.4.3. Elenlaceiónico................................... 12 1.4.4. Elenlacecovalente................................. 13 2. Herramientas matemáticas. 15 2.1. Familia de polinomios ortogonales. . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1.1. Relación de recurrencia. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.1.2. Fórmula de Christoffel-Darboux. . . . . . . . . . . . . . . . . . . . . . . . . 17 2.1.3. Relacióndiferencial................................ 18 2.1.4. Polinomios de Hermite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.1.5. Polinomios generalizados de Laguerre. . . . . . . . . . . . . . . . . . . . . . 20 2.2. Resolución algebraica de la ecuación de Schrödinger. . . . . . . . . . . . . . . . . . 21 3. El método del mallado de Lagrange. 23 3.1. CuadraturadeGauss.................................... 23 3.2. ElmétododeLagrange................................... 24 3.2.1. Mallado de Lagrange-Hermite. . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3. Configuration Localised States: CLS. . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3.1. Propiedades de los CLS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3.2. Aplicación. .................................... 27 3.4. CLS con el oscilador armónico unidimensional. . . . . . . . . . . . . . . . . . . . . 29 4 3.4.1. Validación del método: el potencial de Morse. . . . . . . . . . . . . . . . . . 30 3.5. Comparación método CLS y método de Lagrange-Hermite. . . . . . . . . . . . . . . 37 3.5.1. Conclusiones.................................... 39 3.6. Aplicación a diversos potenciales interatómicos. . . . . . . . . . . . . . . . . . . . . 39 3.6.1. El potencial de Kratzer-Fues. . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.6.2. El potencial de Deng-Fan. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.6.3. El potencial de Varshni. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.6.4. Conclusiones.................................... 46 4. CLS en 3D. 47 4.1. CLS con el oscilador armónico tridimensional. . . . . . . . . . . . . . . . . . . . . . 47 4.2. La molécula de 𝑂2...................................... 50 4.2.1. El potencial de 𝑂2. ................................ 50 4.3. Espectro vibro-rotacional. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.3.1. Estados vibracionales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.3.2. Estadosrotacionales................................ 59 5. Resumen y conclusiones. 67 A. Código en MATLAB. 71 A.1. CLSunidmensional..................................... 71 A.2. CLS3D............................................ 73 Capítulo 1 Nociones básicas de Física Cuántica. La Física Cuántica es la rama de la Física que se encarga del estudio de los sistemas y sus interacciones a escala microscópica. Debido a la naturaleza de los fenómenos que describe, es necesario un marco totalmente distinto al que sigue la Física Clásica. Es por ello que el objetivo de este primer capítulo es entender los fundamentos de dicha disciplina y la manera de estudiarla. 1.1. Los postulados de la Física Cuántica. Siguiendo la interpretación de Copenhague de la Mecánica Cuántica, las bases de dicha teoría se exponen a través de los siguientes postulados: Primer postulado: Todo sistema cuántico 𝑆tiene asociado un espacio de Hilbert , en el que cada vector |𝜓⟩∈representa un posible estado de 𝑆y tal que |||𝜓⟩||= 1. Segundo postulado: Cada observable del sistema tiene asociado un operador lineal y hermítico definido en ∗(espacio dual de ). Tercer postulado: La medida de un observable solo puede dar como resultado un autovalor del operador. Esto es :  𝐴𝜓 =𝑎𝜓 , siendo 𝑎un autovalor de  𝐴y𝜓un autovector asociado a dicho autovalor. Cuarto postulado: La probabilidad de obtener el autovalor 𝑎𝑛asociado a un operador 𝑎 viene dado por: 𝑃(𝑎𝑛) = |⟨𝜙𝑛|𝜓⟩|2, siendo 𝜙𝑛el autovector de  𝐴asociado al autovalor 𝑎𝑛. Quinto postulado: La medida de un sistema altera dicho sistema de forma que, si en la medida se obtiene el valor 𝑎𝑛, el sistema colpasa al autoestado asociado a dicho autovalor, 𝜙𝑛. Sexto postulado: La evolución temporal de la función de onda satisface la ecuación de Schrödinger: 𝑖ℏ𝜕𝜓(𝑟, 𝑡) 𝜕𝑡 = 𝐻(𝑟, 𝑡)𝜓(𝑟, 𝑡),(1.1) donde  𝐻representa el hamiltoniano del sistema y 𝜓la función de onda del sistema. 6 1.2. La ecuación de Schrödinger independiente del tiempo. La resolución de dicha ecuación es crucial para analizar cómo evoluciona el sistema. Desde el punto de vista matemático, es una ecuación compleja en derivadas parciales que no siempre presenta una solución analítica. Es por ello, que la finalidad de este trabajo es su resolución mediante métodos numéricos. El punto de partida será el hamiltoniano. El hamiltoniano del sistema representa la energía total del mismo. El operador hamiltoniano de forma general puede escribirse como una combinación lineal del operador energía cinética,  𝑇(𝑟, 𝑡), y del operador energía potencial,  𝑉(𝑟, 𝑡), es decir:  𝐻(𝑟, 𝑡) =  𝑇(𝑟, 𝑡) +  𝑉(𝑟, 𝑡).(1.2) La acción de dicho operador sobre la función de onda nos proporciona la ecuación de Schrödinger dependiente del tiempo: −ℏ2 2𝑚∇2𝜓(𝑟, 𝑡) + 𝑉(𝑟, 𝑡)𝜓(𝑟, 𝑡) = 𝑖ℏ𝜕𝜓(𝑟, 𝑡) 𝜕𝑡 .(1.3) Como se observa, el operador energía cinética no depende del tiempo de manera explícita. Además, existen múltiples casos de interés en los que el operador energía potencial tampoco depende del tiempo. Consideremos dicho caso: 𝑉(𝑟, 𝑡) = 𝑉(𝑟). De esta forma, el hamiltoniano pasa a ser separable y podemos aplicar separación de variables. 𝜓(𝑟, 𝑡) = 𝜙(𝑟)𝑇(𝑡).(1.4) Como 𝜓satisface la ecuación (1.1), se obtiene: 𝑖ℏ 1 𝑇(𝑡) 𝜕𝑇 (𝑡) 𝜕𝑡 =1 𝜙(𝑟)  𝐻(𝑟, 𝑡)𝜙(𝑟).(1.5) En la ecuación (1.5) suponemos que tanto 𝑇(𝑡)como 𝜙(𝑟)son no nulas. En caso contrario, obtendríamos la solución idénticamente nula, que no tiene ningún interés práctico. Como la parte derecha solo tiene dependencia en 𝑟y la parte izquierda sólo depende de 𝑡, la única solución posible es: 𝑖ℏ 1 𝑇(𝑡) 𝜕𝑇 (𝑡) 𝜕𝑡 =𝐸 , (1.6) 1 𝜙(𝑟)  𝐻(𝑟, 𝑡)𝜙(𝑟) = 𝐸 . (1.7) La constante 𝐸representa la energía total del sistema cuando se encuentra en el estado 𝜙(𝑟). Asimismo, la ecuación (1.6) tiene solución analítica y viene dada por: 𝑇(𝑡) = 𝑁𝑒−𝑖𝐸 ℏ𝑡,(1.8) donde 𝑁es una constante que usaremos para la correcta normalización de la función de onda. Por otro lado, de la ecuación (1.7), se obtiene:  𝐻𝜙(𝑟) = 𝐸𝜙(𝑟).(1.9) 1. nociones básicas de física cuántica. 7 La ecuación (1.9) depende del término de energía potencial y no es más que un problema de autovalores y autovectores asociado al operador hamiltoniano. Recuperando la ecuación (1.4), la solución completa está dada por: 𝜓(𝑟, 𝑡) = 𝜙(𝑟)𝑒−𝑖𝐸 ℏ𝑡,(1.10) teniendo en cuenta que la constante de normalización se ha incluido en 𝜙(𝑟). Lo que realmente tiene interés físico es la densidad de probabilidad. Para estados como (1.10) su expresión viene dada por: |𝜓(𝑟, 𝑡)|2=𝜙∗(𝑟)𝑒+𝑖𝐸 ℏ𝑡𝜙(𝑟)𝑒−𝑖𝐸 ℏ𝑡=𝜙(𝑟)∗𝜙(𝑟).(1.11) Los estados, cuya densidad de probabilidad no depende del tiempo (como es el caso anterior), se denominan estados estacionarios. El hamiltoniano tendrá autovalores 𝐸𝑛con autoestados 𝜙𝑛(𝑛∈ 1,2,3, ...). Dichos autoestados forman una base del espacio, por lo que cualquier solución a la ecuación de Schrödinger se podrá escribir como: 𝜓(𝑟, 𝑡) = ∑ 𝑛 𝑔(𝑛) ∑ 𝑗=1 𝑐𝑗 𝑛𝜙𝑗 𝑛(𝑟)𝑒−𝑖𝐸𝑛 ℏ𝑡,(1.12) con 𝑔(𝑛)la degeneración de dicho autovalor y 𝑐𝑛∈ℂ. 1.3. Propiedades de las funciones de onda. Debido a la interpretación probabilística del cuadrado de la función de onda, dichas funciones deberán cumplir una serie de propiedades. A consecuencia del primer postulado, las funciones de onda han de estar adecuadamente normalizadas, esto es: ∫𝑡𝑜𝑑𝑜 𝑒𝑙 𝑒𝑠𝑝𝑎𝑐𝑖𝑜 |𝜓(𝑟, 𝑡)|2𝑑3𝑟= 1 .(1.13) Para que la expresión (1.13) esté bien definida, la función de onda ha de ser monoevaluada, finita y de cuadrado integrable (para que exista la integral y tenga un valor finito único). Es más, el hecho de satisfacer la ecuación de Schrödinger implica que: La función de onda ha de ser continua, diferenciable y asintóticamente nula. Esta restricción se impone para que exista solución clásica de la ecuación de Schrödinger, o sea, para que la solución sea Ψ(𝑟, 𝑡) ∈ 2(ℝ3;ℝ). La suma de funciones de onda ha de ser una función de onda. Esto se basa en que el espacio de soluciones es un espacio vectorial. Como a lo largo de este trabajo vamos a trabajar al principio con potenciales en una dimensión, vamos a comentar una serie de propiedades exigibles a 𝜙(𝑥), aunque son propiedades extrapolables a mayores dimensiones. Para hablar de los estados del potencial será útil recordar el concepto de "estado ligado". Este hace referencia a un estado que cumple: 𝑉(−∞) < 𝐸 < 𝑉 (+∞). Los estados ligados en una dimensión, además de lo anterior, cumplen que: 14 necesario emplear potenciales como el de Morse o de Kratzer-Fues para reproducir tanto el mínimo como la región en la que los átomos están alejados entre sí. Capítulo 2 Herramientas matemáticas. Con vista a la resolución de la ecuación de Schrödinger para distintos tipos de potenciales, en este trabajo se va a desarrollar un método de resolución numérica fundamentalmente basado en estados localizados en el espacio de configuración. En este capítulo se va a desarrollar las herramientas matemáticas necesarias para, en el siguiente capítulo, poder abordar la construcción de dicho método. 2.1. Familia de polinomios ortogonales. Dado un espacio de Hilbert,, un producto escalar se define como una aplicación, < 𝑢|𝑣 >∶= ×→ℂ, que verifica: ⟨𝑢|𝑢⟩≥0, y < 𝑢|𝑢 >= 0 ⇐⇒ 𝑢= 0 ,∀𝑢∈. ⟨𝑢|𝑣⟩=⟨𝑣|𝑢⟩∗,∀𝑢, 𝑣 ∈. ⟨𝛼𝑢|𝑣⟩=𝛼⟨𝑢|𝑣⟩,∀𝑢, 𝑣 ∈. ⟨𝑢|𝑣+𝑤⟩=⟨𝑢|𝑣⟩+⟨𝑢|𝑤⟩,∀𝑢, 𝑣, 𝑤 ∈. Para lo que nos ocupa estaremos interesados en el siguiente producto escalar: sea 𝑤(𝑥) ∈ ℝtal que 𝑤(𝑥)>0en (𝑎, 𝑏). Entonces, dadas dos funciones 𝑓y𝑔, se puede definir el siguiente producto escalar: ⟨𝑓|𝑔⟩=∫𝑏 𝑎 𝑤(𝑥)𝑔∗(𝑥)𝑓(𝑥)𝑑𝑥 , sobre el espacio 𝐿2 𝑤([𝑎, 𝑏](−∞ ≤𝑎, 𝑏 ≤+∞), formado por las funciones que cumplen : ∫𝑏 𝑎 𝑤(𝑥)|𝑓(𝑥)|2𝑑𝑥 < +∞ . Diremos que dos funciones son ortogonales cuando el valor de su producto escalar es nulo. Se dirá que {𝑣}𝑣≥0forman un conjunto de polinomios ortogonales si se cumple: ∫𝑏 𝑎 𝑤(𝑥)𝑣𝑣′𝑑𝑥 =𝑁𝑣𝛿𝑣𝑣′. siendo 𝑁𝑣≥0. 16 En un conjunto de polinomios ortogonales {𝑣(𝑥); 𝑣≥0}, se cumple que el polinomio 𝑁(𝑥)de grado 𝑁tiene 𝑁raíces reales en el intervalo [𝑎, 𝑏]. Por tanto, si denotamos {𝑥𝑖}𝑁 𝑖=1 a dichas raíces y el polinomio es mónico, lo podemos escribir: 𝑁(𝑥) = 𝑁 ∏ 𝑖=1 (𝑥−𝑥𝑖).(2.1) 2.1.1. Relación de recurrencia. Un conjunto de polinomios ortogonales cumple una relación de recurrencia del tipo: 𝑁+1(𝑥)=(𝛽𝑁𝑥−𝛼𝑁)𝑁(𝑥) − 𝛾𝑁𝑁−1(𝑥).(2.2) Vamos a obtener esta expresión de forma explícita. Reordenando los términos se obtiene: 𝑁+1 − (𝛽𝑁𝑥−𝛼𝑁)𝑁(𝑥)=−𝛾𝑁𝑁−1 .(2.3) Puesto que 𝑁−1(𝑥)es un polinomio de grado 𝑁−1, los coeficientes de 𝑥𝑁+1 y𝑥𝑁en la expresión anterior deben ser nulos. Dichos coeficientes son: 𝑁+1(𝑥)−(𝛽𝑁𝑥−𝛼𝑁)𝑁(𝑥),(2.4) Si denotamos 𝑘𝑣al coeficiente de 𝑥𝑣en el polinomio 𝑣(𝑥), y centrándonos en el coeficiente 𝑥𝑁+1 tenemos: 𝑘𝑁+1 −𝛽𝑁𝑘𝑁= 0 ⇒𝛽𝑁=𝑘𝑁+1 𝑘𝑁 .(2.5) De esta forma, tomando dichos 𝛽𝑁, la expresión 𝑁+1(𝑥) − 𝛽𝑁𝑥𝑁(𝑥)es un polinomio de grado 𝑁que puede escribirse como una combinación lineal de la siguiente forma: 𝑁+1(𝑥) − 𝛽𝑁𝑥𝑁(𝑥) = 𝑁 ∑ 𝑣=0 𝑏𝑣𝑣(𝑥).(2.6) Como el coeficiente asociado a 𝑥𝑁debe ser nulo, si denotamos 𝑘′ 𝑣al coeficiente que acompaña a 𝑥𝑣−1, se cumple: 𝛼𝑁= −𝑏𝑁= − 1 𝑘𝑁 (𝑘′ 𝑁+1 −𝛽𝑁𝑘′ 𝑁) = 𝑘𝑁+1 𝑘2 𝑁 𝑘′ 𝑁−1 𝑘𝑁 𝑘′ 𝑁+1 .(2.7) Luego, el polinomio (2.4) es un polinomio de grado 𝑁− 1 que puede expresarse de la siguiente forma: 𝑁+1(𝑥) − 𝛽𝑁𝑥𝑁(𝑥) = 𝑁−1 ∑ 𝑣=0 𝑐𝑣𝑣(𝑥).(2.8) Usando las propiedades de ortogonalidad podemos determinar el resto de coeficientes. Multiplicamos la expresión por 𝑤(𝑥)𝑚(𝑚∈ [0, 𝑁 − 1]. 2. herramientas matemáticas. 17 El término de la derecha solo sobrevive el término 𝑚: ∫𝑏 𝑎 𝑑𝑥 𝜔(𝑥)𝑚(𝑥) 𝑁−1 ∑ 𝑣=0 𝑐𝑣𝑣(𝑥) = 𝑐𝑚∫𝑏 𝑎 𝑑𝑥 𝜔(𝑥)[𝑚(𝑥)]2=𝑐𝑚𝑚.(2.9) El término izquierdo puede descomponerse en tres integrales usando la linealidad de la integral: ∫𝑏 𝑎 𝑑𝑥 𝜔(𝑥)𝑚(𝑥)𝑁+1(𝑥) + 𝛼𝑁∫𝑏 𝑎 𝑑𝑥 𝜔(𝑥)𝑚(𝑥)𝑁(𝑥) − 𝛽𝑁∫𝑏 𝑎 𝑑𝑥 𝜔(𝑥)𝑚(𝑥)𝑥𝑁(𝑥).(2.10) Las dos primeras integrales se anulan para todo 𝑚≤𝑁− 1 debido a la ortogonalidad. Por otro lado, puesto que 𝑥𝑚(𝑥)es un polinomio de orden 𝑚+1, la tercera integral sólo será distinta de cero si 𝑚+1 = 𝑁. Por lo tanto, los coeficientes 𝑐𝑚deben ser nulos salvo para 𝑚=𝑁−1. De esta manera, obtenemos a partir de la ecuación (1.9) que: 𝛾𝑁= −𝑐𝑁−1 . Para obtener el valor concreto de 𝛾𝑁, estudiamos la última integral en (2.10) para 𝑚=𝑁− 1. Podemos desarrollar el polinomio 𝑥𝑁−1(𝑥)como: 𝑥𝑁−1(𝑥) = 𝑘𝑁−1𝑥𝑁+ 𝑁−2 ∑ 𝑛=0 𝑎𝑛𝑥𝑛+1 .(2.11) Dado que 𝑥𝑛+1 es ortogonal a 𝑁(𝑥)para 𝑛≤𝑁− 2, ninguno de estos términos del desarrollo contribuirá a la integral. El único término relevante será 𝑘𝑁−1𝑥𝑁: ∫𝑏 𝑎 𝑑𝑥 𝜔(𝑥)𝑥𝑁−1(𝑥)𝑁(𝑥) = 𝑘𝑁−1 ∫𝑏 𝑎 𝑑𝑥 𝜔(𝑥)𝑥𝑁𝑁(𝑥) = 𝑘𝑁−1 𝑘𝑁∫𝑏 𝑎 𝑑𝑥 𝜔(𝑥)[𝑁(𝑥)]2=𝑘𝑁−1 𝑘𝑁 𝑁. (2.12) Sustituyendo los valores de las integrales (2.9) y (2.12) en la igualdad, obtenemos el valor del coeficiente 𝑐𝑁−1 y por consiguiente: 𝛾𝑁= −𝑐𝑁−1 =𝑘𝑁+1𝑘𝑁−1 𝑘2 𝑁 𝑁 𝑁−1 .(2.13) Finalmente, la relación de recurrencia queda: 𝑁+1(𝑥) = [𝑘𝑁+1 𝑘𝑁 𝑥−(𝑘𝑁+1 𝑘2 𝑁 𝑘′ 𝑁−1 𝑘𝑁 𝑘′ 𝑁+1)]𝑁(𝑥) − 𝑘𝑁+1𝑘𝑁−1 𝑘2 𝑁 𝑁 𝑁−1 𝑁−1(𝑥).(2.14) 2.1.2. Fórmula de Christoffel-Darboux. A partir de la relación de recurrencia anterior puede obtenerse la fórmula de Christoffel-Darboux, expresión crucial en el desarrollo de los estados localizados en el espacio de la configuración (CLS). 18 Multiplicando (2.14) por 𝑁(𝑦): 𝑁+1(𝑥)𝑁(𝑦) = [𝑘𝑁+1 𝑘𝑁 𝑥−(𝑘𝑁+1 𝑘2 𝑁 𝑘′ 𝑁−1 𝑘𝑁 𝑘′ 𝑁+1)]𝑁(𝑥)𝑁(𝑦)−𝑘𝑁+1𝑘𝑁−1 𝑘2 𝑁 𝑁 𝑁−1 𝑁−1(𝑥)𝑁(𝑦). (2.15) Intercambiando 𝑥por 𝑦para restar ambas expresiones, se obtiene: 𝑁+1(𝑥)𝑁(𝑦) − 𝑁(𝑥)𝑁+1(𝑦) = 𝑘𝑁+1 𝑘𝑁 (𝑥−𝑦)𝑁(𝑥)𝑁(𝑦) + 𝑘𝑁+1𝑘𝑁−1 𝑘2 𝑁 𝑁 𝑁−1 [𝑁(𝑥)𝑁−1(𝑦) − 𝑁−1(𝑥)𝑁(𝑦)].(2.16) Podemos iterar esta expresión y desarrollar de la misma manera el último término: 𝑁(𝑥)𝑁−1(𝑦) − 𝑁−1(𝑥)𝑁(𝑦) = 𝑘𝑁 𝑘𝑁−1 (𝑥−𝑦)𝑁−1(𝑥)𝑁−1(𝑦) + 𝑘𝑁𝑘𝑁−2 𝑘2 𝑁−1 𝑁−1 𝑁−2 [𝑁−1(𝑥)𝑁−2(𝑦) − 𝑁−2(𝑥)𝑁−1(𝑦)].(2.17) Con sucesivas iteraciones se obtiene la fórmula de Christoffel-Darboux: 𝑘𝑁+1 𝑘𝑁 (𝑥−𝑦) 𝑁 ∑ 𝑣=0 𝑁 𝑣 𝑣(𝑥)𝑣(𝑦) = 𝑁+1(𝑥)𝑁(𝑦) − 𝑁(𝑥)𝑁+1(𝑦).(2.18) 2.1.3. Relación diferencial. Los polinomios ortogonales además cumplen la siguiente relación: 𝑔2(𝑥)𝑑 𝑑𝑥𝑁(𝑥) = 𝑔1(𝑥)𝑁(𝑥) + 𝑔0(𝑥)𝑁−1(𝑥).(2.19) Los valores de 𝑔0,𝑔1y𝑔2se pueden encontrar en la referencia [4]. 2.1.4. Polinomios de Hermite. Los polinomios de Hermite físicos 𝐻𝑛(𝑥)son una familia de polinomios ortogonales con gran utilidad en diversos problemas físicos. Dichos problemas suelen estar descritos por ecuaciones diferenciales en las que interviene el operador Laplaciano. Matemáticamente, son casos particulares del problema de Sturm-Liouville, en los que los polinomios de Hermite resultan útiles para construir soluciones de dichas ecuaciones. Son polinomios bien definidos en ℝy conforman una familia de polinomios ortogonales con la función peso 𝑤(𝑥) = 𝑒−𝑥2. Vienen dados de forma explícita por: 2. herramientas matemáticas. 19 𝐻2𝑚(𝑥) = (−1)𝑚 𝑚 ∑ 𝑝=0 (2𝑚)! (𝑚−𝑝)!(2𝑝)!(−4𝑥2)𝑝,(2.20) 𝐻2𝑚+1(𝑥) = (−1)𝑚2𝑥 𝑚 ∑ 𝑝=0 (2𝑚+ 1)! (𝑚−𝑝)!(2𝑝+ 1)!(−4𝑥2)𝑝.(2.21) El coeficiente líder es 𝑘𝑛= 2𝑛. De esta forma, la constante de normalización : 𝑛=∫∞ −∞ 𝑑𝑥 𝑒−𝑥2𝐻2 𝑛(𝑥)=2𝑛𝑛!√𝜋 . (2.22) Además, satisfacen la siguiente relación de recurrencia: 𝐻𝑁+1(𝑥)=2𝑥𝐻𝑁(𝑥)−2𝑁𝐻𝑁−1(𝑥),(2.23) Se puede comprobar fácilmente, usando las expresiones (2.21), que se cumple la siguiente relación diferencial: 𝑑 𝑑𝑥𝐻𝑁(𝑥)=2𝑁𝐻𝑁−1(𝑥).(2.24) Algunos de los polinomios de Hermite físicos se representan en la figura 2.1 y son: 𝐻0(𝑥) = 1 ; (2.25) 𝐻1(𝑥)=2𝑥;(2.26) 𝐻3(𝑥)=8𝑥3− 12𝑥;(2.27) 𝐻4(𝑥) = 16𝑥4− 48𝑥2+ 12 ; (2.28) 𝐻5(𝑥) = 32𝑥5− 160𝑥3+ 120𝑥 . (2.29) Figura 2.1: Representación de los polinomios de Hermite de grado 𝑣∈ {0, ..., 4} 20 2.1.5. Polinomios generalizados de Laguerre. Al igual que los polinomios de Hermite, los polinomios generalizados de Laguerre son una familia de polinomios ortogonales entre sí que resultan ser parte de soluciones de ecuaciones tan trascendentales en la Física como la ecuación de la parte de onda radial en el átomo de hidrógeno. Los polinomios 𝐿(𝛼) 𝑛(𝑥)son ortogonales con la función peso 𝑥𝛼𝑒−𝑥(𝛼 > −1) y dominio (0,∞). Estos se definen como siguen: 𝐿(𝛼) 𝑛(𝑥) = 𝑛 ∑ 𝑝=0 (−1)𝑝Γ(𝑛+𝛼+ 1) Γ(𝑝+𝛼+ 1)Γ(𝑛−𝑝+ 1)Γ(𝑝+ 1)𝑥𝑝.(2.30) De la expresión anterior se concluye que los coeficientes de los términos principales (𝑘𝑛,𝛼 coeficiente del término de grado n y 𝑘′ 𝑛,𝛼 coeficiente del término de grado n-1): 𝑘𝑛,𝛼 =(−1)𝑛 𝑛!, 𝑘′ 𝑛,𝛼 =(−1)𝑛−1(𝑛+𝛼) (𝑛− 1)! . Por otro lado, el valor de la integral es bien conocido: 𝑛,𝛼 =∫∞ 0 𝑑𝑥 𝑥𝛼𝑒−𝑥[𝐿(𝛼) 𝑛(𝑥)]2=Γ(𝑛+𝛼+ 1) Γ(𝑛+ 1) .(2.31) Finalmente, los polinomios de Laguerre generalizados satisfacen la siguiente ecuación diferencial: 𝐿(𝛼) 𝑁+1(𝑥)=− 1 𝑁+ 1 [𝑥− (2𝑁+𝛼+ 1)]𝐿(𝛼) 𝑁(𝑥) − 𝑁+𝛼 𝑁+ 1𝐿(𝛼) 𝑁−1(𝑥).(2.32) Relación diferencial: 𝑥𝑑 𝑑𝑥𝐿(𝛼) 𝑁(𝑥) = 𝑁𝐿(𝛼) 𝑁(𝑥)−(𝑁+𝛼)𝐿(𝛼) 𝑁−1(𝑥).(2.33) Algunos de los polinomios de Laguerre generalizados se representan en la figura 2.2 y son: 𝐿1 0(𝑥) = 1; (2.34) 𝐿1 1(𝑥) = 2 − 𝑥;(2.35) 𝐿1 2(𝑥) = 1 2𝑥2− 3𝑥+ 3; (2.36) 𝐿1 3(𝑥)=−1 6𝑥3+ 2𝑥2− 6𝑥+ 4; (2.37) 𝐿1 4(𝑥) = 1 24𝑥4−5 6𝑥3+ 5𝑥2− 10𝑥+ 5; (2.38) 2. herramientas matemáticas. 21 Figura 2.2: Representación de los polinomios de Laguerre, 𝐿(1) 𝑣,𝑣∈ {0, ..., 4} 2.2. Resolución algebraica de la ecuación de Schrödinger. En esta sección vamos a ilustrar que la resolución de la ecuación de Schrödinger independiente del tiempo es equivalente a diagonalizar el hamiltoniano en una base del espacio. Como hemos visto, la ecuación de Schrödinger independiente del tiempo (  𝐻𝜙𝑘=𝐸𝑘𝜙𝑘) es un problema de autovalores y autovectores asociado al operador hamiltoniano. Supongamos que conocemos una base del espacio de Hilbert: {𝜙𝑛}𝑛≥0. Al ser base, es un sistema generador del espacio, luego, cualquier solución del problema se puede expresar como combinación lineal de ellas. Esto es: 𝜑𝑘= ∞ ∑ 𝑛=0 𝑐(𝑘) 𝑛𝜙𝑛.(2.39) Introduciendo dicha expresión en la ecuación:  𝐻 ∞ ∑ 𝑛=0 𝑐(𝑘) 𝑛𝜙𝑛=𝐸𝑘 ∞ ∑ 𝑛=0 𝑐(𝑘) 𝑛𝜙𝑛(2.40) resulta evidente que conocer 𝑐(𝑘) 𝑛es análogo a conocer 𝜑𝑘.Para determinar dichos coeficientes vamos a hacer uso de la ortogonalidad (autoestados asociados a autovalores distintos son ortogonales entre sí): ∞ ∑ 𝑛=0 𝑐(𝑘) 𝑛⟨𝑚|𝐻|𝑛⟩=𝐸𝑘 ∞ ∑ 𝑛=0 𝑐(𝑘) 𝑛⟨𝑚|𝑛⟩=𝐸𝑘𝑐(𝑘) 𝑚,(2.41) 22 donde hemos usado la notación de Dirac para facilitar los cálculos. 𝜙𝑛→|𝑛⟩ ∫𝑑𝑉 𝜙∗ 𝑚𝜙∗ 𝑛→⟨𝑚|𝑛⟩ ∫𝑑𝑉 𝜙∗ 𝑚𝐻𝜙∗ 𝑛→⟨𝑚|𝐻|𝑛⟩ Usando el formalismo matricial, la ecuación (2.41) se puede escribir como: ⎛⎜⎜⎜⎝ 𝐻0,0𝐻0,1𝐻0,2⋯ 𝐻1,0𝐻1,1𝐻1,2⋯ 𝐻2,0𝐻2,1𝐻2,2⋯ ⋮ ⋮ ⋮ ⋱⎞⎟⎟⎟⎠⎛⎜⎜⎜⎜⎝ 𝑐(𝑘) 0 𝑐(𝑘) 1 𝑐(𝑘) 2 ⋮ ⎞⎟⎟⎟⎟⎠ =𝐸𝑘⎛⎜⎜⎜⎜⎝ 𝑐(𝑘) 0 𝑐(𝑘) 1 𝑐(𝑘) 2 ⋮ ⎞⎟⎟⎟⎟⎠ Por todo ello, para resolver la ecuación de Schrödinger independiente del tiempo, vamos a tratar de diagonalizar el hamiltoniano en una base del espacio de Hilbert solución. Esto es, obtener la matriz asociada al operador y resolver el problema de autovalores y autovectores. En la práctica, como las bases del espacio son de dimensión infinita, se usan bases truncadas. Claramente, los resultados que se obtendrán presentarán una desviación con respecto a los reales. Aún así, se pueden obtener resultados muy precisos y exactos. En el capítulo 3, profundizaremos en el método de mallado de Lagrange. Se abordará dicho método mediante las funciones de Lagrange y con las funciones de oscilador armónico a través del método de los CLS. Capítulo 3 El método del mallado de Lagrange. El método del mallado de Lagrange es un método variacional aproximado que resuelve las ecuaciones no en un continuo, sino en una malla usando la aproximación de cuadratura de Gauss. Se basa en polinomios ortogonales que presentan buenas propiedades tales como diferenciabilidad infinita. El método presenta una forma simple y, con el número de puntos adecuados, se puede lograr una gran precisión en los cálculos. Con dicho método es posible calcular las energías, las funciones de onda y distintos observables de estados ligados descritos por cualquier potencial. 3.1. Cuadratura de Gauss. La cuadratura de Gauss aproxima una integral mediante una suma ponderada de valores de la función en ciertos puntos específicos (raíces de un polinomio ortogonal). Esta técnica es exacta cuando la función integrada es un polinomio de grado hasta 2𝑁−1, donde 𝑁es el número de puntos de malla. Si la función no es lo suficientemente suave (por ejemplo, si tiene discontinuidades), el error puede ser significativo. Sea (𝑎, 𝑏)un intervalo donde está definida una familia de polinomios ortogonales {𝑘}con respecto a la función peso 𝑤(𝑥). Supongamos que los polinomios pueden escribirse de la siguiente forma: 𝑃𝑁(𝑥) = 𝑘𝑁𝑥𝑁+𝑘′ 𝑁𝑥𝑁−1 +... (3.1) Consideremos la siguiente integral: ∫𝑏 𝑎 𝑃𝑘(𝑥)𝑤(𝑥)𝑃𝑘′𝑑𝑥 =ℎ𝑘𝛿𝑘,𝑘′(3.2) con ℎ𝑘el módulo cuadrado de 𝑘. Una cuadratura estándar de Gauss consistiría en dividir el intervalo de integración en 𝑁puntos, siendo cada punto un cero del polinomio de grado 𝑁. Se tiene entonces un conjunto de puntos {𝑥𝑖}𝑁 𝑖=1 tales que 𝑁(𝑥𝑖) = 0 ∀𝑖= 1, ..., 𝑁. Sea 𝑤𝑖los pesos asociados a cada 𝑥𝑖,𝑖∈ {1, ..., 𝑁}. Por la aproximación de cuadratura: ∫𝑏 𝑎 𝑓(𝑥)𝑤(𝑥)𝑑𝑥 = 𝑁 ∑ 𝑘=1 𝑤𝑘𝑓(𝑥𝑘).(3.3) Para las aplicaciones en el método de mallado de Lagrange, se prefiere una forma modificada de la cuadratura que evita el uso explícito de la función peso 𝑤(𝑥). La integral se reescribe como una 30 donde los solapamientos son ⟨𝑗, 𝑣|𝐶𝐿𝑆;𝑗, 𝑠⟩=𝑤1∕2 𝑗𝑠 𝑁−1∕2 𝑗𝑣 𝐻𝑣(𝑥𝑠),(3.43) siendo 𝑁𝑗𝑣 la constante de normalización de las autofunciones de (3.41) y 𝑤−1 𝑗𝑠 =𝑗 𝑁𝑗𝑗−1 (𝐻𝑗−1(𝑥𝑠))2.(3.44) El hecho de que los solapamientos resulten ser funciones analíticas, nos permitirá evaluar los elementos de matriz de la ecuación (3.38) muy rápidamente desde el punto de vista computacional. Tomando una base formada por las 𝑗= 10, podemos definir los CLS del oscilador armónico unidimensional, representados en la figura 3.1. Figura 3.1: Representación de los estados CLS asociados al oscilador armónico unidimensional para una base de 𝑗= 10 estados. Se observa que para cada 𝑠∈ {1, ..,10}, se obtiene una función que se encuentra localizada en el espacio en torno al punto 𝑥=𝑥𝑠. Además, dicho polinomio se anula en todos los ceros de 𝐻10 excepto en 𝑥=𝑥𝑠. Asimismo, de la figura anterior se comprueba que si 𝑥𝑠es un cero del polinomio de Hermite de grado 10, entonces −𝑥𝑠también lo será. Eso provocará que al representar todos los CLS en una misma figura haya una simetría especular en torno al eje de ordenadas. Finalmente, comprobamos las buenas propiedades analíticas de los CLS. Se tratan de funciones suaves que además están normalizadas. 3.4.1. Validación del método: el potencial de Morse. Para evaluar la efectividad del método numérico que hemos construido basado en los CLS, vamos a obtener los autoestados y autofunciones del potencial de Morse y vamos a comparar la solución analítica con los resultados obtenidos. Como queremos comprobar la potencia del método, se ha 3. el método del mallado de lagrange. 31 optado por el uso del potencial de Morse adimensional (las longitudes están en unidades de longitud del oscilador y las energías en unidades de ℏ𝜔). De esta manera, el potencial de Morse adquiere la siguiente expresión: 𝑉𝑀(𝑥) ℏ𝜔 =𝐷 ℏ𝜔 [(1 − 𝑒−𝛽𝑥)2− 1].(3.45) Donde: −𝐷es el valor del potencial en el mínimo. La coordenada 𝑥=𝑟−𝑟𝑒 𝑎0 está en unidades de la longitud de oscilador 𝑎0; siendo 𝑟la coordenada física y 𝑟𝑒la posición de equilibrio. 𝛽es un parámetro adimensional relacionado con el alcance de la interacción.. Podemos caracterizar los valores de los parámetros del potencial a través de 𝑗𝑀. Este representa el número de estados ligados que soporta el potencial. Fijado dicho valor, las expresiones de 𝛽y𝐷 vienen dadas por: 𝛽=1 √𝑗𝑀+1 2 ,(3.46) 𝐷 ℏ𝜔 =𝛽2 2(𝑗𝑀+1 2)2 =1 2(𝑗𝑀+1 2).(3.47) Los estados del potencial se pueden denotar por dos números cuánticos que se denotan |𝑗𝑀𝑚⟩. La etiqueta 𝑚se relaciona con la energía del estado más ligado (correspondiente a 𝑚=𝑗). Sin embargo, en Física Molecular es más común usar un número alternativo: 𝑣=𝑗𝑀−𝑚. Dicho nuevo número no es más que el número de fonones anarmónicos y toma valores en {0, ..., 𝑗𝑀−1}. El estado fundamental se asocia a 𝑣= 0. Con esta notación, las energías (en unidades de ℏ𝜔) de los niveles del potencial de Morse son: 𝜖𝑗𝑀,𝜈 =𝐷 ℏ𝜔 −(𝜈+1 2)1 2𝑗𝑀+ 1 (𝜈+1 2)2 .(3.48) Los autoestados de este potencial son bien conocidos y tienen la siguiente expresión: 𝜓𝑗𝑀,𝑣(𝑥) = ⟨𝑥|𝑀;𝑗𝑀, 𝑣⟩=−1∕2 exp(−𝑦∕2) 𝑦𝑗𝑀−𝑣𝐿(2𝑗𝑀−2𝑣) 𝑣(𝑦),(3.49) donde se emplea la variable de Morse 𝑦= (2𝑗𝑀+1) exp(−𝛽𝑥), los 𝐿(𝑎) 𝑛son los polinomios de Laguerre generalizados, y la constante de normalización es: −1∕2 =√𝛽(2𝑗𝑀− 2𝑣)Γ(𝑣+ 1) Γ(2𝑗𝑀−𝑣+ 1) .(3.50) En la figura 3.2 representamos el potencial de Morse para 𝑗𝑀= 10. Como las energías y las autofunciones de dicho potencial son conocidas, hemos representado algunos de los autoestados sobre su energía correspondiente. Si la base utilizada en el método CLS fuera infinita, se hallarían las energías de forma exacta. Sin embargo, esto es inviable en la realidad. Cabe preguntarse si con una base finita se pueden obtener las energías de forma precisa. 32 Figura 3.2: Potencial de Morse con 𝑗𝑀= 10. Se representan algunas de las energías del potencial con su correspondiente función de onda asociada. Figura 3.3: Reprentación de la energía (en unidades de ℏ𝜔 ) frente al tamaño de la base usada de los CLS para distintos niveles del potencial de Morse adimensional con 𝑗𝑀= 10. Se repsenta en línea discontinua los valores de referencia. Cada color corresponde a un nivel distinto del potencial. Es por ello que en la figura 3.3 se han representado los valores de energía (en unidades de ℏ𝜔) obtenidos mediante el formalismo de los CLS para distintos tamaños de la base. Se aprecia que las energías convergen bastante rápido a las energías exactas. Asimismo, dicha convergencia es mayor cuanto más ligado está el estado. Esto es debido a que en dichos estados la base de oscilador armónico se asemeja mejor al potencial de Morse. 3. el método del mallado de lagrange. 33 Figura 3.4: Representación de la función de onda de onda del tercer estado excitado del potencial de Morse adimensional con 𝑗𝑀= 10 usando distintos tamaños de la base. En línea discontinua negra se representa la función de onda exacta Figura 3.5: Representación de la función de onda de onda del tercer estado excitado del potencial de Morse adimensional con 𝑗𝑀= 10 usando distintos tamaños de la base. En línea discontinua negra se representa la función de onda exacta 34 En cuanto a las autofunciones, ocurre algo similar que con la energía. A medida que se aumenta el número de elementos de la base, las funciones de onda se asemejan a las exactas. Para ver esto, se ha representado en la figura 3.4 y 3.5 el cuarto estado excitado para distintos tamaños de la base. Mientras que para una base de 𝑗= 50 elementos no se obtienen prácticamente desviaciones, para una base de 𝑗= 10 elementos las desviaciones son significativas. De hecho, la función obtenida para 𝑗= 10 presenta más ceros de los que debería tener el estado. Como se puede apreciar en los resultados, la energía necesita de un menor número de elementos de la base para que el valor numérico presente menores desviaciones del real. Si recordamos que: 𝐸𝑛=⟨𝐻⟩|𝜙𝑛⟩,(3.51) es decir, si la base utilizada abarca la mayor parte del dominio en el que está definida la función de onda, el error cometido es pequeño. Esto es porque la principal contribución a la energía procede de aquellas zonas donde la densidad de probabilidad es más elevada, o sea, la región en torno a la posición de equilibrio (es decir, donde el potencial alcanza su mínima energía). En el caso del potencial considerado, esta región está adecuadamente representada por los elementos de la base. Sin embargo, la función de onda experimental en cada punto del espacio depende de que la base sea lo suficientemente extensa para describir la zona, en otras palabras, que los ceros del polinomio de Hermite usado tomen valores en dicha región. De no ser así, la función obtenida puede anularse aunque no le corresponda o presentar oscilaciones que realmente no posee. Para ilustrar lo anterior, se va a elegir una base de 𝑗= 50 elementos, de forma que no se van a obtener todas las energías ni las funciones de onda de forma exacta, pero es lo suficientemente grande como para poder reproducir bien los primeros estados ligados. En la tabla 3.1 se recogen las energías exactas del potencial de Morse junto a las obtenidas mediante el formalismo de los CLS construido con una base de oscilador de 𝑗= 50. v𝐸𝑒𝑥𝑎𝑐𝑡𝑎 𝑣𝐸𝐶𝐿𝑆 𝑣 0 -4,7619 -4.7619 1 -3,8571 -3.8571 2 -3,0476 -3.0476 3 -2,333 -2.3333 4 -1.7143 -1.7143 5 -1.1905 -1.1902 6 -0.7619 -0.7518 7 -0.4286 -0.3381 8 -0.1905 0.1355 9 -0.0476 0.6880 Tabla 3.1: Energías exactas y energías obtenidas con una base de CLS con 𝑗= 50 estados para el potencial de Morse adimensional con 𝑗𝑀=10. De la observación directa de la tabla, se aprecia que el método es capaz de reproducir de forma fidedigna los estados energéticos de los primeros cuatro estados. A partir de ahí, se obtienen energías cercanas pero cometiendo errores relativos significativos (21 % para el estado 𝑣= 7). Finalmente, el formalismo de los CLS predice que los estados 𝑣= 8 y𝑣= 9 son estados no ligados cuando realmente sí lo son. Cuanto mayor número de elementos de la base de CLS, mayor será la precisión de las energías. Asimismo, como estamos usando una base de 𝑗= 50, los CLS que estamos tomando son funciones 3. el método del mallado de lagrange. 35 de onda localizadas en torno a los ceros de 𝐻50(𝑥). El último cero se alcanza en 𝑥≈ 9,18, pero como bien vemos en la figura 3.2 las funciones de onda asociadas a los últimos estados del potencial de Morse se extienden más allá de dicho punto. Por tanto, para estudiar esos estados se necesita una base mayor. En cuanto a las autofunciones, representadas en la figura 3.6, se comportan igual que lo ocurrido para las energías. Para los primeros estados el método es capaz de reproducirnos dichas curvas. A pesar de ello, a partir del estado 𝑣= 5 ya se aprecian desviaciones en la caída de la función. Finalmente, para 𝑣= 9 prácticamente son funciones distintas. De nuevo, la base usada se anula para 𝑥 > 9,18, por lo que es imposible que pueda estudiarse el comportamiento de la función en dicha región del espacio (𝑥 > 9,18). 36 Figura 3.6: Autonfunciones obtenidas por el método de los CLS (línea continua azul) para el potencial de Morse con 𝑗𝑀= 10 con una base 𝑗= 50. La línea discontinua representa las autofunciones exactas de dicho potencial. Diagrama de Birge-Sponer. Como ya se ha mencionado, el potencial de Morse surge como una forma más realista de explicar el enlace molecular. Principalmente, este refleja el efecto de la anarmonicidad. En una molécula real los niveles energéticos no se encuentran equidistantes sino que disminuyen a medida que aumentamos el estado vibracional ( incrementamos el valor de 𝑣). La energía del potencial de Morse viene dada por la expresión : 𝐸𝑣= −𝐷+(𝑣+1 2)−1 2𝑗𝑀+ 1(𝑣+1 2)2 , 𝑣 = 0,1, ..., 𝑗𝑀− 1 .(3.52) Por lo tanto, la diferencia de energía entre niveles contiguos puede calcularse como: Δ𝐸𝑣=𝐸𝑣+1 −𝐸𝑣= 1 − (𝑣+ 1) 2𝑗𝑀+ 1 .(3.53) El diagrama de Birge-Sponer no es más que la representación de Δ𝐸𝑣frente a 𝑣. De ella, podemos estimar el número máximo de estados ligados. Para ello, consideramos el máximo estado el punto de corte del diagrama con el eje de abscisas. No tendría sentido considerar Δ𝐸𝑣<0, pues implicaría que aumentar el estado vibracional conlleva disminuir la energía, lo cual es incoherente. 3. el método del mallado de lagrange. 37 En este caso, para el potencial de Morse adimensional, el máximo estado ligado es: 𝑣𝑚𝑎𝑥 ≤𝑗𝑀− 1 .(3.54) En la figura 3.7 se representa el diagrama de Birge-Sponer para el potencial de Morse con 𝑗𝑀= 10. El espaciado energético exacto sigue una clara dependencia lineal con 𝑣, pero los obtenidos por CLS con 𝑗= 50 sólo los primeros estados lo cumplen (recordemos que los CLS para los dos últimos estados ni siquiera nos confirman que son estados ligados, para esta base). Figura 3.7: Diagrama de Birge-Sponer para el potencial de Morse con 𝑗𝑀= 10. Se representa Δ𝐸𝑣= 𝐸𝑣+1 −𝐸𝑣(en unidades de ℏ𝑤) frente a 𝑣tanto para las energías obtenidas del método de los CLS para una base con 𝑗= 50 estados y para las exactas. El diagrama anterior nos pone de manifiesto que los CLS predicen mejor la energía para estados más ligados. Esto es porque en dichos estados la forma del potencial es más cercana a la del oscilador. Finalmente, hacer notar que el diagrama de Birge-Sponer para el Morse predice que el número máximo de estados ligados que predice el potencial es 𝑣𝑚𝑎𝑥 = 9, como corresponde a la elección 𝑗𝑀= 10. 3.5. Comparación método CLS y método de Lagrange-Hermite. Tanto el método de los CLS como el método de Lagrange-Hermite presentan la misma filosofía: resolución numérica de potenciales no analíticos. En este apartado, vamos a analizar en qué difieren y por qué se ha optado por el método de los CLS. En primer lugar, ambos métodos se basan en la discretización del espacio y el uso de cuadratura para la resolución algebraica de la ecuación de Schrödinger. La diferencia principal es la forma de abordar el término cinético en el hamiltoniano. 38 Para comparar ambos métodos, vamos a tomar como referencia las energías del potencial de Morse para 𝑗𝑀= 10. En la figura 3.8 se representa la energía obtenida por el método de los CLS y por el método de Lagrange-Hermite, tomando distintos tamaños de la base. Figura 3.8: Energía del estado fundamental y del primer estado excitado para el potencial de Morse con 𝑗𝑀= 10 en función del tamaño de la base de oscilador. Se incluyen las energías exactas para ver la convergencia del método. Con cls nos referimos al resultado obtenido por el formalismo de los CLS; y, LH, usando las funciones de Lagrange-Hermite. Ambos métodos convergen bastante rápido. Basta con una base con menos de 10 elementos para poder reproducir la energía del estado fundamental y del primer estado excitado. Bien es cierto que el método de los CLS presenta una mayor desviación para bases de muy pocos elementos. Sin embargo, dado que la convergencia es bastante rápida y que en realidad no resulta muy costoso incluso tomar bases mayores ( por ejemplo 𝑗= 200 ), esto no será un problema en la práctica. Realmente, las funciones de Lagrange-Hermite y los CLS son las mismas salvo signo. Veamos cómo aborda cada método el término de energía cinética. Por un lado, usando funciones de Lagrange se obtienen: 𝑇𝐺 𝑖≠𝑗= (−1)𝑖−𝑗2 (𝑥𝑖−𝑥𝑗)2,(3.55) 𝑇𝐺 𝑖𝑖 =1 3(2𝑁+1−𝑥2 𝑖).(3.56) Por otro lado, el método de los CLS introduce el término de potencial del oscilador, por lo que: 𝐻=𝑇+𝑉(𝑥) = 𝑇+𝑉𝑂𝐴(𝑥)+[𝑉(𝑥) − 𝑉𝑂𝐴(𝑥)] .(3.57) De esta forma, se evalúa englobando la parte de energía cinética a la energía total del oscilador armónico que está bien definida y es conocida, y en vez de evaluar el potencial que describe al sistema, evalúa otro potencial equivalente (la diferencia entre dichos potenciales), que desde el punto de vista computacional tiene el mismo costo. Sin embargo, cuando intentamos generalizar el método del mallado de Lagrange usando funciones de Lagrange a tres dimensiones, el término de energía cinética se vuelve más complejo: 𝑇𝐺 𝑖𝑗 =⟨𝑓𝑖|𝑝2 2𝑚|𝑓𝑗⟩,(3.58) 3. el método del mallado de lagrange. 39 y la evaluación de dichos elementos de matriz ya no es inmediato. De hecho, no hay un consenso de cuál es la mejor opción para abordar el método de mallado en este caso (puede verse esta cuestión más extendida en la referencia [1]). Aquí es cuando se evidencia la conveniencia del método de los CLS. Como la energía del oscilador armónico tridimensional está bien definida, se puede llevar a cabo la evaluación mostrada en la expresión (3.57) 3.5.1. Conclusiones. A lo largo del capítulo se ha desarrollado el formalismo de los CLS y se ha ilustrado cómo utilizarlo para obtener los autovalores y las autofunciones de un hamiltoniano. Usando el potencial de Morse, se ha podido profundizar en la calidad de los resultados. De todo ello podemos concluir: El método converge. Usando un tamaño de la base lo suficientemente grande se ha comprobado que las energías y autofunciones resultantes tienden, rápidamente, a las expresiones analíticas conocidas. Aumentar el tamaño de la base provoca que los puntos en torno a los cuales se localizan los CLS cubran mayor parte del espacio. Para obtener resultados precisos necesitamos tener cubierta toda la región en la que la función de onda sea no nula. El método que implica el uso de los CLS es más flexible para mayores dimensiones que el método de Lagrange. Por todo ello, en el resto del trabajo se ha optado por el método de los CLS. 3.6. Aplicación a diversos potenciales interatómicos. Como hemos visto en el capítulo 1, los potenciales interatómicos que describen los enlaces en moléculas diatómicas son curvas suaves que tienden a cero cuando la distancia entre los constituyentes que forman la molécula es elevada; y, divergen cuando las componentes se encuentran ínfimamente cerca, de forma que la curva entre dichos valores extremos es una curva regular que presenta un mínimo. Es por ello que el oscilador armónico no sería una buena elección de potencial para describir al sistema. De hecho, no podría explicar fenómenos tan comunes de las moléculas como la disociación, pues, ni siquiera existen estados no ligados. Por tanto, se recurre a potenciales como el de Morse, el de Kratzer-Fues, el de Deng-Fan o el de Varshni, que nos permiten reproducir de manera más fidedigna el comportamiento de la interacción molecular. 3.6.1. El potencial de Kratzer-Fues. El potencial de Kratzer-Fues es un potencial que pertenece a la familia de los potenciales que describen moléculas diatómicas. Su expresión es la siguiente: 𝑉(𝑞) = 𝐷(𝑞2 0 𝑥2−2𝑞0 𝑥).(3.59) El significado de las constantes es: 46 3.6.4. Conclusiones. El objetivo fundamental de esta sección ha sido no sólo ilustrar cómo el formalismo de los CLS puede ser usado para la resolución de la ecuación de Schrödinger independiente del tiempo para potenciales unidimensionales, sino también mostrar que dicho método es equivalente; incluso puede llegar a ser más eficiente que el método de Lagrange-Hermite. Se ha podido verificar su validez gracias a la ayuda de un potencial de solución analítica y se han explorado las limitaciones del método en función del número de elementos de la base elegida. Principalmente, el método tiene una excelente convergencia, siempre y cuando la base tenga un alcance suficiente como para cubrir los puntos del espacio en los que está definida la función de onda. En el siguiente capítulo, se va a mostrar que además es posible generalizar dicho formalismo para potenciales en 3D y se va a aplicar un sistema físico real. De esta forma, podremos llevar a cabo un estudio energético más rico: espectro molecular vibro-rotacional. Volveremos a comprobar la eficacia del método y abordaremos la dependencia de la elección de la longitud de oscilador para dotar a las soluciones de las unidades adecuadas. Capítulo 4 CLS en 3D. 4.1. CLS con el oscilador armónico tridimensional. Para el estudio de sistemas reales físicos, la variable 𝑟que denota la distancia internuclear, no toma valores en la recta real sino que 𝑟 > 0. Por ello, los CLS en una dimensión no resultan convenientes para su estudio al tomar valores en ℝ. Presentamos aquí el método de los CLS usando un oscilador armónico en 3D. Con él, solucionaremos el tema de la variable de longitud, y además, nos dará los grados de libertad suficientes para llevar a cabo un estudio más exhaustivo, en cuanto a las vibraciones y rotaciones de la molécula. El Hamiltoniano del oscilador armónico isótropo en 3D es:  𝐻𝑂𝐴 =1 2𝜇𝑝2+𝜇𝜔2 2𝑟2.(4.1) Siendo 𝜇la masa reducida del sistema y 𝜔la frecuencia. Las autofunciones de este potencial pueden separarse en dos componentes, una radial y una angular (armónicos esféricos): 𝜑𝑛𝓁𝑚(𝑟) = 𝑅𝑛𝓁(𝑟)𝑌𝓁𝑚(𝜃, 𝜙),(4.2) donde 𝑛= 1,2,…,𝓁= 0,…, 𝑛 − 1,𝑚= −𝓁,−𝓁+ 1,…,𝓁. La componente angular de las autofunciones, 𝑌𝓁𝑚(𝜃, 𝜙), son los armónicos esféricos; funciones conocidas y que están tabuladas. Además, esta parte angular es común a todos los potenciales centrales. Para la construcción de los CLS nos centraremos en las autofunciones radiales, en concreto en la función 𝑢𝑛𝓁(𝑟) = 𝑟𝑅𝑛𝓁(𝑟)cuya expresión explícita es: 𝑢𝑛𝓁(𝑟) = ⟨𝑟|𝑛, 𝓁⟩=−1∕2 𝑛𝓁𝑟𝓁+1 exp (−𝑟2 2)𝐿(𝓁+1∕2) 𝑛(𝑟2),(4.3) donde empleamos la coordenada 𝑟en unidades de la longitud del oscilador 𝑎0=√ℏ 𝜇𝜔 , los 𝐿𝛼 𝑛son los polinomios de Laguerre generalizados de grado 𝑛y parámetro 𝛼, y la constante de normalización es: 𝑛𝓁=√𝜋(2𝑛+ 2𝓁+ 1)!! 22𝑛+𝓁(2𝑛)!! .(4.4) 48 La constante de normalización se puede redefinir de la siguiente forma : 0𝓁=2𝓁+ 1 4Γ(𝓁+1 2),(4.5) 𝑛𝓁=(2𝑛+𝓁+ 1)(𝑛+𝓁)! 4𝑛𝛽(𝑛, 𝓁+3 2).(4.6) Haciendo el cambio de variable 𝑦=𝑟2, podemos reescribir estas autofunciones en la forma descrita en el formalismo : 𝑢𝑛𝓁(𝑟) = −1∕2 𝑛𝓁𝐹(𝑦)𝐿(𝓁+1∕2) 𝑛(𝑦),(4.7) tomando 𝐹(𝑦) = 𝑦𝑙+1 2𝑒−𝑦∕2 .(4.8) Particularizando el formalismo de la sección 3.3 para el caso del oscilador armónico cuántico unidimensional, las funciones CLS, en términos de los autoestados del oscilador se pueden escribir: |𝐶𝐿𝑆;𝑗, 𝑠, 𝑙⟩= 𝑗−1 ∑ 𝑣=0⟨𝑗, 𝑣, 𝑙|𝐶𝐿𝑆;𝑗, 𝑠, 𝑙⟩|𝑗, 𝑣, 𝑙⟩,(4.9) donde los solapamientos son ⟨𝑗, 𝑣, 𝑙|𝐶𝐿𝑆;𝑗, 𝑠, 𝑙⟩=𝑤1∕2 𝑗𝑠𝑙 𝑁−1∕2 𝑣𝑙 𝐿(𝑙+1∕2) 𝑣(𝑥𝑠),(4.10) siendo 𝑁𝑗𝑣 la constante de normalización de las autofunciones de (3.41) y 𝑤−1 𝑗𝑠𝑙 =−𝑗 𝑁𝑗𝑗−1 (𝑗+𝑙+1 2)(𝐿(𝑙+1∕2) 𝑗−1 (𝑥𝑠))2.(4.11) El hecho de que los solapamientos resulten ser funciones analíticas, nos permitirá evaluar los elementos de matriz de la ecuación (3.38) muy rápidamente desde el punto de vista computacional. Ahora 𝑙representa un grado de libertad que nos va a permitir obtener distintas familias de CLS. Fijado el valor de 𝑗, este parámetro a medida que aumenta provocará que los ceros se desplacen a la derecha. Como bien se observa en la figura 4.1, a medida que aumentamos el valor de 𝑙, las funciones de onda se van desplazando a la derecha. 4. cls en 3d. 49 Figura 4.1: CLS radiales en 3D con una base de oscilador de 𝑗= 10 y valores de 𝑙= 0,3,6 50 4.2. La molécula de 𝑂2. El oxígeno es un elemento químico no metálico caracterizado por su número atómico 𝑍= 8. Su configuración electrónica en el estado fundamental es (1𝑠)2(2𝑠)2(2𝑝)4. Debido a que su capa más externa no se encuentra completa, y siguiendo la regla del octeto, el oxígeno tenderá a formar enlaces con otros elementos para pasar a un estado más estable energéticamente. Entre una de las múltiples posibilidades, un oxígeno puede interaccionar con otro oxígeno y dar lugar a la molécula de 𝑂2. En dicha molécula, los dos oxígenos que la constituyen interaccionan mediante un doble enlace covalente. A temperatura ambiente, la molécula constituye un gas incoloro, inoloro e insípido que interviene en procesos vitales como la fotosíntesis o la respiración celular. Desde el punto de vista de la teoría de enlace molecular, los orbitales 2𝑝de ambos oxígenos solapan, de manera que se forman orbitales moleculares enlazantes y antienlazantes. Los electrones que se encuentran en los orbitales enlazantes están más bajos en energía y eso contribuye a la estabilidad de la molécula. Por el contrario, los electrones que se encuentran en los orbitales antienlazantes, más energéticos, tenderán a inestabilizar la molécula. No debemos olvidar que los sistemas físicos tienden a minimizar su energía. Como la contribución total neta a la energía resulta negativa, se forma el doble enlace. Aunque, a priori, pueda parecer un problema simple, este resulta ser bastante complejo, pues, un hamiltoniano que incluya tanto las interacciones de los electrones entre sí, el término de espín-órbita o las vibraciones y rotaciones moleculares, resulta en un problema de solución no analítica. Es por ello que se recurre a cálculos ab-initio para el estudio de la molécula de 𝑂2. Aunque en la actualidad, el desarrollo de las tecnologías computacionales resulte una gran herramienta para llevar a cabo el cálculo numérico de multitud de problemas, el elevado número de parámetros que suelen describir a dichos potenciales o la necesidad de tomar una base elevada para obtener una buena precisión conllevan un coste considerable en tiempo y almacenamiento. Por todo ello, en este capítulo vamos a partir de un potencial obtenido por cálculos ab-initio y le aplicaremos el formalismo de los CLS en 3D para poder estudiar el espectro vibro-rotacional. La gran ventaja de este método es que, al resolver de forma algebraica la ecuación de Schrödinger independiente del tiempo, no supone un gasto computacional demasiado elevado. 4.2.1. El potencial de 𝑂2. Los estados rotacionales y vibracionales de la molécula de dioxígeno han sido ampliamente estudiados y medidos experimentalmente. En la referencia [8] se obtiene la superficie de energía potencial con métodos ab-initio y el resultado se ajustó a una combinación de funciones gaussianas: 𝑉(𝑟) = 7 ∑ 𝑘=0 𝑎𝑘𝑒−𝛼𝛽𝑘𝑟2.(4.12) Los parámetros obtenidos en [8] están en la tabla 4.1 y el potencial se ha dibujado en la figura 4.2. 𝛼0.785 𝑎3154738.09175 𝛽1.307 𝑎4-215074.85646 𝑎0-2388.5641690 𝑎5214799.54567 𝑎118086.977116 𝑎6-148395.42850 𝑎2-71760.197585 𝑎773310.78145 Tabla 4.1: Parámetros extraídos de la referencia [8] para la obtención de energía potencial (PES). Las unidades de 𝑎𝑘son en miliHartree, 𝛽es adimensional y 𝛼en  𝐴 4. cls en 3d. 51 Figura 4.2: Superficie de energía potencial (PES) para 𝑂2extraído de la referencia [8]. La primera tarea que tenemos que abordar antes de ejecutar el formalismo desarrollado en el capítulo anterior es la elección adecuada de las unidades. Mientras que las aplicaciones mostradas en capítulos anteriores eran adimensionales, pues su finalidad era mostrar la potencia del método, ahora estamos tratando con un potencial que describe un sistema físico real y que, por tanto, debe ir acompañado de sus respectivas unidades. La expresión de la longitud de oscilador: 𝑎0=√ℏ 𝜇𝜔 =√(ℏ𝑐)2 𝜇𝑐2ℏ𝜔 ,(4.13) siendo en el caso de la molécula de 𝑂2,𝜇𝑂2= 15,9949146221∕2 𝑢.𝑚.𝑎. Figura 4.3: Representación del espaciado entre el nivel vibracional 𝑣= 1 y𝑣= 0, denotado por 𝐺1−𝐺0, frente a la longitud de oscilador empleado para la construcción de los CLS. En base a la representación anterior, se elige como longitud de oscilador aquella que proporcione el resultado más cercano al valor de la referencia. 52 En la figura 4.3 se ha representado el espacio entre los dos primeros niveles de energía vibracionales (usando el método de los CLS) frente a la longitud de oscilador usada. Como se puede apreciar, usando como 𝑎0= 0,1452  𝐴, se obtiene que el espaciado obtenido por el método es de 1564,06 𝑐𝑚−1, mientras que en la referencia marca que dicho valor debe ser de 1563,91 𝑐𝑚−1. Es decir, el error cometido es menor del 0,01 %. Consecuentemente, usaremos para el resto de cálculos ℏ𝜔 = 200 𝑐𝑚−1, que son las unidades de energía que nos proporcionan tal longitud del oscilador. 4.3. Espectro vibro-rotacional. En el capítulo 1, se explicó que el movimiento de las moléculas es más variado que el de los átomos al existir las rotaciones y las vibraciones. Asimismo, se estimó la energía asociada a cada contribución del movimiento. Es por ello que en el espectro vibro-rotacional de la molécula de 𝑂2 aparecerán sobre cada nivel vibracional distintos niveles rotacionales. En nuestro desarrollo, será el número cuántico 𝑣= 0,1,2,3, ... el que caracterice los niveles vibracionales, mientras que el número cuántico 𝐿= 0,1,2,3, ... representará las rotaciones. Las energías asociada a cada nivel vibro-rotacional de una molécula pueden ajustarse de forma general por: 𝐸𝑣,𝐿 =𝐸𝑣,0+𝐵𝑣𝐿(𝐿+ 1) − 𝐷𝑣𝐿2(𝐿+ 1)2,(4.14) donde 𝐸𝑣,0=ℏ𝜔𝑒(𝑣+1 2)−𝑥𝑒ℏ𝜔𝑒(𝑣+1 2)2 ,(4.15a) 𝐵𝑣=𝐵𝑒−𝛼𝑒(𝑣+1 2),(4.15b) 𝐷𝑣=𝐷𝑒+𝛽𝑒(𝑣+1 2).(4.15c) ℏ𝜔𝑒es la contribución armónica a las vibraciones y 𝑥𝑒ℏ𝜔𝑒da cuenta de la anarmonicidad. Por otra parte, se denomina constante rotacional a 𝐵𝑣y constante de deformación centrífuga a 𝐷𝑣. Los valores de 𝐷𝑒y𝐵𝑒son los valores de las constantes rotacionales en el estado vibracional fundamental (𝑣= 0 y despreciando el efecto de la rotación). Por último, 𝛼𝑒y𝛽𝑒son constantes positivas que representan los efectos de la anarmonicidad y la dispersión centrífuga en la molécula. 4.3.1. Estados vibracionales. En primer lugar, comenzaremos estudiando los estados vibracionales puros, esto es, los estados con 𝐿= 0. Para ello, tomaremos las autofunciones del oscilador armónico en 3D con 𝑙= 0 y usaremos, según hemos discutido anteriormente, una base de 𝑗= 200 elementos, una longitud de oscilador de 𝑎0= 0,1452  𝐴y, consecuentemente, ℏ𝜔 = 200 𝑐𝑚−1 . En la tabla 4.2, se exponen los espaciados energéticos entre estados consecutivos. Como se puede observar, los valores obtenidos por el método de los CLS coinciden prácticamente con los de la referencia. El error cometido es menor del 1 %, con lo cual, los datos obtenidos son bastante precisos y además refuerzan la robustez del método de los CLS. Examinando la tabla, los estados vibracionales que presentan un mayor error corresponden a 𝑣′𝑠más altas, lo cual es lógico, ya que, a medida que aumentamos el nivel vibracional, el potencial que describe a la molécula se aleja más de la forma de oscilador armónico usada, lo cual hace que la aproximación sea peor. 4. cls en 3d. 53 𝑣(𝐺𝑣,0−𝐺𝑣−1,0)𝑟𝑒𝑓 (𝑐𝑚−1) (𝐺𝑣,0−𝐺𝑣−1,0)𝐶𝐿𝑆 (𝑐𝑚−1)𝜖𝑟𝑒𝑙( %) 1 1.563,91 1.564,06 <0,01 2 1.538,50 1.538,65 <0,01 3 1.513,41 1.513,55 <0,01 4 1.488,65 1.488,79 <0,01 5 1.464,23 1.464,36 <0,01 6 1.440,14 1.440,27 <0,01 7 1.416,37 1.416,50 <0,01 8 1.392,91 1.393,03 <0,01 9 1.369,71 1.369,82 <0,01 10 1.346,74 1.346,85 <0,01 15 1.233,42 1.233,51 <0,01 20 1.116,17 1.116,23 <0,01 25 983,28 983,31 <0,01 30 815,54 815,50 <0,01 35 486,76 486,51 0,05 40 98,88 98,81 0,07 Tabla 4.2: Espaciado energético entre niveles consecutivos,𝑣y𝑣−1 (𝐺𝑣−𝐺𝑣−1) en unidades de 𝑐𝑚−1. Se comparan los resultados obtenidos con el formalismo de los CLS con los valores de la referencia [8] . Además, en la tabla 4.2 se observa el efecto de la anarmonicidad. Básicamente, la anarmonicidad es debida a que el enlace entre los dos núcleos no viene dado por un rotor rígido; más bien, viene descrito por una constante elástica 𝑘𝑒. La rotación en la molécula produce una fuerza centrífuga que genera un aumento en la distancia internuclear. Como consecuencia de ello, los niveles energéticos no se encuentran interespaciados, sino que, a medida que aumenta 𝑣, la separación energética entre los niveles disminuye; esto es porque la corrección anarmónica se vuelve más grande. Figura 4.4: Representación del potencial de 𝑂2junto con los espaciados energéticos obtenidos en la tabla 4.2. Manifiesta que la separación entre niveles no es constante. 54 En la figura 4.4 se han representado los datos obtenidos en la tabla 4.2. Como se puede apreciar, los primeros estados representados (𝑣= 0, ..., 10), prácticamente se encuentran equiespaciados. Esto es debido a que la corrección centrífuga para estados más cercanos al fondo del pozo es menor (el potencial molecular se asemeja al del oscilador armónico). Por el contrario, para los últimos estados estos se agrupan de tal manera que tienden a formar un cuasi-continuo. Asimismo, el nivel vibracional 𝑣= 40 se encuentra muy cercano a la energía de disociación de la molécula. Es por ello que existirá un estado vibracional límite, a partir del cual la molécula se rompe. Si representamos el diagrama de Birge-Sponer (figura 4.5), se observa que para los niveles más bajos de 𝑣, el espaciado energético sigue una tendencia lineal. Sin embargo, en torno al estado 𝑣= 25 se pierde dicha linealidad. Sumado a esto, si no existiera la corrección centrífuga y se mantuviera la tendencia lineal, se predeciría que existen 𝑣𝑚𝑎𝑥 = 67 estados ligados. Sin embargo, si representamos todos los datos de diferencias energéticas entre niveles consecutivos, a partir del estado 𝑣= 43 se pasa a aumentar el espaciado energético. Esto apunta a que, a partir de dicho estado, no estamos describiendo un estado ligado, por lo que podemos afirmar que la molécula de 𝑂2tiene, aproximadamente, 43 estados ligados vibracionales puros. La figura anterior manifiesta el hecho de que no se puede despreciar el efecto rotacional en las moléculas. Figura 4.5: Diagrama de Birge-Sponer para las energías vibracionales puras con una base CLS con 𝑗= 200 estados de la molécula de 𝑂2. Como para las vibraciones sabemos que Δ𝐸𝑣=𝐸𝑣+1 −𝐸𝑣=ℏ𝜔𝑒− 2𝑥𝑒ℏ𝜔𝑒(𝑣+ 1) ,(4.16) podemos llevar a cabo una regresión lineal de Δ𝐸𝑣frente a 𝑣+1 para así obtener ℏ𝜔𝑒y𝑥𝑒ℏ𝜔𝑒(véase figura 4.6). Los resultados obtenidos son: ℏ𝜔𝐶𝐿𝑆 𝑒= (1582,4±1,0) 𝑐𝑚−1 ,(4.17) 𝑥𝑒ℏ𝜔𝐶𝐿𝑆 𝑒= (11,66 ± 0,07) 𝑐𝑚−1 .(4.18) 4. cls en 3d. 55 Figura 4.6: Representación del espaciado energético, Δ𝐸en 𝑐𝑚−1 frente a 𝑣+ 1. Da cuenta de que el modelo vibracional describe bien los estados poco excitados vibracionalemente. Comparando con los datos de la referencia [10], ℏ𝜔𝑟𝑒𝑓 𝑒= 1580,4𝑐𝑚−1 y𝑥𝑒ℏ𝜔𝐶𝐿𝑆 𝑒= 12,07 𝑐𝑚−1. Los resultados tienen un error relativo menor al 3 % para el caso de 𝑥𝑒ℏ𝜔𝑒; y, menor del 0,1 % para ℏ𝜔𝑒. De lo anterior, podemos inferir que 𝑥𝑒≈ 7,39 ⋅10−3 (adimensional). Aunque la corrección centrífuga es pequeña, esta no es despreciable y como depende directamente de 𝑣, se hace mayor a medida que esta aumenta. Asimismo, podemos relacionar la frecuencia de vibración con la constante elástica: 𝑘𝑒=𝜇𝜔2 𝑒.(4.19) Los resultados obtenidos son: 𝑘𝐶𝐿𝑆 𝑒= 1182 𝑁∕𝑚;(4.20) 𝑘𝑟𝑒𝑓 𝑒= 1177 𝑁∕𝑚 , (4.21) que posee un error relativo menor al 0,5 %. Autofunciones. En la figura 4.7 se han representado los estados vibracionales puros usados en la tabla 4.2. Dichas funciones de onda cumplen las propiedades descritas en el capítulo 1. Por inspección sobre las curvas, se obtienen funciones continuas, sin puntos angulosos (diferenciables) y monovaluadas. En general, son curvas suaves. Se observa que a medida que aumentamos el nivel vibracional (incrementamos el estado excitado), la función de onda se anula en un punto más (aumenta en una unidad el número de ceros). Asimismo, la zona clásicamente permitida no contiene a ±∞, y las funciones son asintóticamente nulas (condición necesaria para su correcta normalización). La falta de simetría en las funciones de onda es claramente debida a la asimetría del potencial. 62 lejanos y como la subida para valores 𝑟 < 𝑟𝑣es mucho más abrupta, provoca que el valor esperado de la posición aumente. 𝑣 𝑟𝐶𝐿𝑆 𝑣( 𝐴)𝑟𝑟𝑒𝑓 𝑣( 𝐴) 0 1,211084 ±0,000004 1,210892 1 1,217855 ±0,000004 1,217580 2 1,224739 ±0,000004 1,224340 3 1,231732 ±0,000004 1,231188 4 1,238829 ±0,000005 1,237972 5 1,246027 ±0,000005 1,244824 6 1,253320 ±0,000005 1,252256 7 1,260715 ±0,000005 1,259647 8 1,268231 ±0,000005 1,267041 9 1,275843 ±0,000005 1,274512 10 1,283572 ±0,000005 1,282217 15 1,324414 ±0,000006 1,323152 20 1,371124 ±0,000006 1,369768 25 1,429563 ±0,000007 1,427141 30 1,514527 ±0,000008 1,513446 Tabla 4.5: Distancia de equilibrio de la molécula de 𝑂2, a partir de los valores de la constante rotacional obtenida mediante los valores de la referncia como el formalismo de los CLS. Reiters que el comportamiento de la molécula no debe basarse únicamente en el modelo de rotor-rígido. De hecho, como sabemos que a partir de 𝑣= 43 estamos en un estado no ligado, podemos estimar cuál sería la distancia máxima a la que los dos oxígenos podrían formar un estado ligado (antes de disociarse). Dicha distancia resulta ser 𝑟𝑚𝑎𝑥 ≈ 2,78  𝐴. En adición, de la pendiente de la regresión lineal (4.22) podemos obtener los valores de 𝐷𝑣, la constante centrífuga. Dicha constante surge como una corrección al modelo de rotor rígido, que tiene en cuenta que la molécula se deformará más cuanto más rápido gire. Los datos obtenidos se representan en la tabla 4.6. 𝑣 𝐷𝑟𝑒𝑓 𝑣(10−6𝑐𝑚−1)𝐷𝐶𝐿𝑆 𝑣(10−6𝑐𝑚−1)𝜖𝑟𝑒𝑙( %) 0 4,854 4,788 1 1 4,957 4,806 3 2 5,000 4,821 4 3 5,000 4,834 3 4 5,049 4,844 4 5 5,145 4,851 6 6 5,006 5,021 <1 7 4,953 5,153 4 8 4,990 4,867 2 9 5,012 4,873 3 10 4,990 4,879 2 15 5,134 4,971 3 20 5,557 5,283 5 25 6,648 6,093 8 30 8,400 8,246 2 Tabla 4.6: Constantes centrífugas en función del estado vibracional, tomando las energías procedentes de la referencia [8] y del método de los CLS. Dan cuenta de la importancia del efecto anarmónico en cada nivel vibracional. 4. cls en 3d. 63 Aunque las desviaciones con respecto a los datos experimentales de la referencia [8] no son muy altas (para la mayoría son menores a 5 %), son más altas que para el caso de las constantes rotacionales. De igual forma, como el error asociado al resultado obtenido es menor que el número de cifras significativas que la referencia tiene, se usan tantas como éstas tienen. Aunque la tendencia general de la constante centrífuga es aumentar a medida que aumenta el estado vibracional, esto no siempre se da (véase 𝑣= 4 y𝑣= 5). El aumento de dicha constante significa que la corrección centrífuga se hace más presente (coherente con el hecho de que nos encontramos con estados menos ligados que el fundamental) y que el enlace se vuelve un tanto más flexible. Sin embargo, la corrección centrífuga es mucho menor que la corrección de rotor no rígido y eso se refleja en los órdenes de magnitud de las constantes. Retomando la expresión (4.15c), podemos llevar a cabo un ajuste y obtener tanto el valor de 𝐷𝑒 como 𝛽𝑒(véase figura 4.9). Para la constante rotacional se obtiene: 𝐷𝐶𝐿𝑆 𝑒= (4,4±0,2) ⋅10−6𝑐𝑚−1 ,(4.28) 𝐷𝑟𝑒𝑓 𝑒= (4,5±0,2) ⋅10−6𝑐𝑚−1 .(4.29) Por otro lado, para 𝛽𝑒: 𝛽𝐶𝐿𝑆 𝑒= (0,082 ± 0,016) ⋅10−6𝑐𝑚−1 ,(4.30) 𝛽𝑟𝑒𝑓 𝑒= (0,089 ± 0,015) ⋅10−6𝑐𝑚−1 .(4.31) Figura 4.9: Valores de la constante de distorsión centrífuga 𝐷𝑣obtenida por el método de los CLS para una base de 𝑗= 200 estados y obtenidas de la referencia. Mediante un ajuste lineal podemos inferir de dichos datos tanto 𝐷𝑒como 𝛽𝑒. Se tiene en línea azul continua el ajuste lineal empleando los datos obtenidos por el formalismo de los CLS y en línea roja continua la regresión usando los datos de la referencia. Asimismo, se marca en en cruces azules los datos experimentales obtenidos por el formalismo de los CLS y con círculos rojos los datos de la referencia [8]. 64 Del mismo modo que ocurría con 𝛼𝑒y𝐵𝑒, el valor de 𝛽𝑒es bastante más pequeño que el valor de 𝐷𝑒. Igualmente, se tiene que 𝛽𝑒es una constante que tiene en cuenta que la distorsión centrífuga se hace más presente para estados vibracionales superiores. Como era de esperar, el valor de 𝐷𝑒es menor que el del resto, pues, 𝐷𝑒es la constante de deformación centrífuga que se obtendría en el ideal de molécula completamente rígida. Como dicha constante es inversamente proporcional a la constante elástica, se predice una mayor constante elástica en dicha situación ( el modelo no puede explicar que el enlace se vuelve más blando para estados vibracionales mayores). Los resultados anteriores deben entenderse más bien de forma cualitativa o como una primera referencia, puesto que el ajuste lineal realizado no presenta un buen coeficiente de determinación. Esto se ve claramente en que los puntos no siguen una tendencia lineal. Como podemos relacionar la frecuencia de oscilación con la constante centrífuga, podemos obtener una primera noción de la rigidez del enlace. Dicha expresión es: 𝐷𝑣=4𝐵3 𝑣 (ℏ𝜔𝑣)2.(4.32) En la tabla 4.7 encontramos dichas frecuencias de oscilación. Vemos que para estados de altos valores de 𝑣le corresponden frecuencias de oscilación muy pequeñas. Ahí, el enlace se ha vuelto más largo (los oxígenos se encuentran más separado) y por ende el enlace es más débil. LLegara un momento en el que el enlace sea tan débil que la molécula se disocie. Como primera aproximación podemos estimar dicha frecuencia calculando la que le correspondería al último estado ligado. Este valor es de ℏ𝜔𝑚𝑎𝑥 = 28,68 𝑐𝑚−1. vℏ𝜔𝐶𝐿𝑆 𝑣(𝑐𝑚−1)ℏ𝜔𝑟𝑒𝑓 𝑣(𝑐𝑚−1)𝜖( %) 0 1574.985536 1564.901655 1 1 1545.931771 1523.181323 1 2 1517.586926 1491.632917 2 3 1489.942771 1466.881987 2 4 1462.965323 1435.880185 2 5 1436.595498 1399.062902 3 6 1387.560368 1393.24942 <0.1 7 1345.828188 1376.171292 2 8 1360.240783 1347.197865 1 9 1335.296084 1320.736075 1 10 1310.459743 1299.925497 1 15 1181.818986 1166.26139 1 20 1033.158539 1010.398296 2 25 848.8190996 816.7847172 4 30 613.6386079 609.2742158 1 Tabla 4.7: Frecuencia de oscilación de las moléculas en distintos estados vibracionales. Los resultados mostrados son los obtenidos empleando el método de los CLS y los valores de la referencia [8]. Constituyen una noción básica del comportamiento vibracional de la molécula. 4. cls en 3d. 65 Autofunciones. De igual forma, podemos obtener las autofunciones asociadas a los estados rotacionales. Estas presentan la misma forma que la de los estados vibracionales puros, pero se encuentran desplazadas a la derecha (véase la figura 4.10). El origen de este desplazamiento reside en la rotación de la molécula. Para mayores valores de 𝐿, el efecto de la fuerza centrífuga es mayor, provocando que los oxígenos se separen. Figura 4.10: Autofunción del estado vibracional 𝑣= 41 con 𝐿= 0,𝐿= 5 y𝐿= 10. Se ha obtenido empleando el formalismo de los CLS tridimensional usando una base de 200 estados. Esclarece el efecto de la rotación en las moléculas. El hecho de que la función de onda tenga una apariencia prometedora refleja que con la base elegida hemos descrito todo el dominio de la función. Si recordamos, los puntos en torno al cual están localizados los estados CLS usados se construyen mediante los ceros del polinomio 𝐿(0) 200(𝑦), 𝐿(5) 200(𝑦)y𝐿(10)(𝑦) 200 . La coordenada 𝑦se relaciona con la coordenada radial del oscilador armónico 3D como 𝑦=𝑟2 𝑜𝑠𝑐, que además se encuentra en unidades de la longitud del oscilador. Por todo ello, los CLS están localizados en torno a los puntos 𝑟𝑠=𝑎0√𝑦𝑠. Lo que ocurre es que con la elección de la base empleada no estamos obteniendo el número de puntos adecuados para poder describir dicha región, y de ahí que adquiera esa forma la curva. Precisando el asunto anterior con cifras, el último cero de 𝐿(10) 200 (𝑦)toma el valor 𝑦≈ 788,38. Expresandolo en unidades de 𝑟, se obtiene que 𝑟= 4,08  𝐴, valor a partir del cual la función de onda se anula. Capítulo 5 Resumen y conclusiones. El propósito fundamental de este trabajo ha sido desarrollar un método numérico que nos permita resolver la ecuación de Schrödinger independiente del tiempo. Al tratarse de un sistema cuántico, se han repasado conceptos básicos de Física Cuántica tales como los Postulados de la Mecánica Cuántica o las propiedades que siguen las funciones de onda, centrándonos en aquellas que dependen únicamente de una variable, ya que los potenciales a tratar en principio son unidimensionales. Asimismo, hemos hecho énfasis en Física Molecular, estudiando de forma cualitativa la forma de los potenciales que describen a los enlaces. Una vez que adquirimos las nociones básicas del fundamento físico del problema que queremos tratar, pasamos a introducir los conceptos matemáticos para plantear el método. Entre ellos, cabe destacar el de familia de polinomios ortogonales, que será la base de la construcción de los CLS junto con la relación de Christoffel-Darboux. Asimismo, se muestra el caso de los polinomios de Hermite y Laguerre generalizado, dos familias de polinomios ortogonales relevantes en la resolución de sistemas cuánticos. Finalmente, se acaba el capítulo mostrando la resolución algebraica de la ecuación de Schrödinger, la esencia del método. En el capítulo 3, se aborda y se desarrolla el método que vamos a usar. Se definen nuevos conceptos como el de cuadratura de Gauss. Se presenta brevemente el método de mallado de Lagrange y se demuestra que los estados CLS son una alternativa ventajosa al método Lagrange-Hermite, prinicpalemnte por su clara generalización a dimensiones superiores. Por ello, se decide usar el método de mallado basado en los CLS. Se da un procedimiento explícito de cómo construir los CLS y usarlos en la resolución de la ecuación de Schrödinger independiente del tiempo. Se comienza con el potencial de Morse adimensional (cuya resolución analítica es conocida) y se comprueba la eficacia del método. Obtenemos así parámetros físicos de interés como la energía y las funciones de onda. En adición, en base a ellos, podemos explicar algunas de las limitaciones del modelo. Para finalizar el capítulo, presentamos resultados CLS para otros potenciales moleculares de interés sin solución analítica. Por último, en capítulo 4, primero, extendemos el método CLS a 3D y, a continuación, abordamos un problema realista: la molécula de 𝑂2. Empezamos describiendo el potencial y la elección de las unidades adecuadas para tratar el problema. De hecho, la aplicación del formalismo en 3D nos posibilitará estudiar el espectro vibro-rotacional de la molécula. Así, seremos capaces de determinar algunos parámetros cruciales en la rotación y vibración, como las constantes vibracionales ℏ𝜔𝑒y 𝑥𝑒ℏ𝜔𝑒y las rotacionales 𝐵𝑣y𝐷𝑣, entre otras. En adición, podremos ver el efecto del giro en las moléculas. Uno de los motivos por los que los resultados obtenidos son de tan alta calidad, además de la efectividad del método, si no que el potencial empleado reproduce de manera muy fiel el comportamiento de la molécula. Y esto es gracias a su derivación por complejos cálculos 𝑎𝑏 𝑖𝑛𝑖𝑡𝑖𝑜. De esta forma, una vez conocido el potencial, el método presentado es capaz de obtener con precisión los autoestados y las autofunciones de forma sencilla y rápida, empleando un tamaño razonable de 68 la base. Más aún, este trabajo se completa con un anexo en el que se incluye el código de MATLAB para su implementación. El desarrollo autónomo del código ha permitido una mayor comprensión del trabajo, así como una corroboración de las ventajas computacionales que presenta. Como conclusión, el formalismo de los CLS nos proporciona resultados muy precisos, siempre y cuando la elección del número de estados de la base sea la adecuada. De esta manera, al sustituir integrales por evaluaciones de funciones construidas en los puntos adecuados, nos ayuda a reducir el coste en cálculo computacional y tiempo. Bibliografía [1] Baye, D. (2015). The Lagrange-mesh method.Physics Reports, 565, 1. [2] Carvajal, M., Arias, J. M., & Gómez-Camacho, J. (1999). Configuration localized Morse wave functions: Application to vibrational transitions in anharmonic diatomic molecules. Physical Review A, 59, 1852. [3] Pérez-Bernal, F., Arias, J. M., Carvajal, M., & Gómez-Camacho, J. (2000). Configuration localized wave functions: General formalism and applications to vibrational spectroscopy of diatomic molecules. Physical Review A, 61, 042504. [4] Abramowitz, M., & Stegun, I. A. (Eds.). (1965). Handbook of mathematical functions with formulas, graphs, and mathematical tables (Vol. 55). Dover Publications. [5] Doman, B. G. S. (2015). The classical orthogonal polynomials. World Scientific. [6] Morse, P. M. (1929). Diatomic molecules according to the wave mechanics II. Vibrational levels. Physical Review, 34, 57. [7] Merrill, J. R. (1972). Introductory quantum mechanics with the computer. American Journal of Physics, 40, 138–143. [8] Bytautas, L., Matsunaga, N., & Ruedenberg, K. (2010). Accurate ab initio potential energy curve of O2. II. Core-valence correlations, relativistic contributions, and vibration-rotation spectrum. Journal of Chemical Physics, 132, 074307. [9] Herzberg, G. (1959). Spectra of Diatomic Molecules. Van Nostrand Reinhold. [10] Resnick, R., & Eisberg, R. (1978). Física cuántica: Átomos, moléculas, sólidos, núcleos y partículas (E. Flores, Trad.). Editorial Limusa. [11] Griffiths, D. J. (2005). Introduction to quantum mechanics (2nd ed.). Pearson Prentice Hall. [12] Bransden, B. H., & Joachain, C. J. (1995). Physics of atoms and molecules (2nd ed.). Longman. Apéndice A Código en MATLAB. En este anexo vamos a recoger los códigos principales que hemos usado en MATLAB para desarrollar este trabajo. Para la representación gráfica de los resultados se ha usado el comando plot. Asimismo, para el ajuste lineal en los diagramas de Birge-Sponer se ha usado la herramienta cftool. A.1. CLS unidmensional. Primero, obtenemos los ceros del polinomio de Hermite de grado n. Creamos otra función que representa a los pesos. Se ha nombrado como wjs.