Relatividad general con Matlab

Universidad de Sevilla Facultad de Física Relatividad general con Matlab Trabajo Fin de Grado Autor: Mario Misas Arcos Tutores: Alberto Tomás Pérez Izquierdo y Carlos Soria del Hoyo Junio 2021 Resumen El objetivo de este proyecto es la resolución numérica de las ecuaciones de las geodésicas características de distintos espacio-tiempos en el marco de la Teoría General de la Relatividad. Para ello se hará uso del software MATLAB, en primer lugar para la construcción de dichas ecuaciones utilizando herramientas de cálculo simbólico y en segundo lugar para su resolución mediante métodos numéricos incluidos en el software. Posteriormente se realizará una comparación de los resultados obtenidos con diversas predicciones teóricas y algunos de los tests clásicos más relevantes de la relatividad general. Agradecimientos Debo en primer lugar expresar mi gratitud hacia aquellos que han hecho posible el proyecto que se expone en esta memoria. A mis tutores, Carlos Soria del Hoyo y Alberto Tomás Pérez Izquierdo, por su inestimable dirección durante el desarrollo del proyecto y todo el material que me han facilitado para su desarrollo. En especial a este último, quien durante el tercer curso del grado, mientras ejercía como alumno interno bajo su tutela, con su dedicación y entrega contribuyó a que experimentara mi formación con una motivación que me acompañará durante el resto de mis estudios y carrera profesional. A mis amigos y amigas, en especial a Beatriz, Laura, Alejandro, Miguel, Andrés, Rafael, Javier y María, quienes han sido un imprescindible apoyo sin el cual no podría asegurar haber llegado a donde me encuentro hoy. Y por último a mi familia, por nunca dejar de creer en mí y animarme a seguir recorriendo el camino con cada vez más fuerza. Sin todos vosotros, la realización de este proyecto habría resultado imposible. Gracias. Índice 1. Introducción 1 2. Fundamento teórico 3 2.1. Variedad, métrica y tensor métrico . . . . . . . . . . . . . . . . . . . . . . . 3 2.2. Laecuacióngeodésica............................... 5 2.3. El principio de covariancia y las ecuaciones de campo de Einstein . . . . . . 7 2.4. La métrica de Schwarzschild . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.5. LamétricadeGödel ............................... 12 3. Metodología 16 3.1. Construcción de la ecuación geodésica . . . . . . . . . . . . . . . . . . . . . . 16 3.2. Condicionesiniciales ............................... 21 3.3. Obtención de resultados: tests de relatividad general . . . . . . . . . . . . . . 24 3.4. Obtención de resultados: universo de Gödel . . . . . . . . . . . . . . . . . . . 30 4. Resultados 37 4.1. Tests de relatividad general . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2. UniversodeGödel ................................ 42 5. Discusión de resultados 45 5.1. Tests de relatividad general . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.2. UniversodeGödel ................................ 46 6. Conclusiones 50 1 INTRODUCCIÓN 1. Introducción El camino hacia la relatividad general y motivación para el uso de MATLAB La interacción gravitatoria es la suma gobernante de los fenómenos a gran escala en nuestro Universo. Como tal, su presencia en la vida cotidiana se hace palpable en casi cualquier escenario imaginable, siendo a la vez responsable de sucesos tan a priori desligados como lo pudieran parecer el ciclo día-noche o la caída de los frutos de un árbol. No fue sin embargo hasta 1687 cuando aparecería, en una publicación de título Philosophiæ naturalis principia mathematica , la formulación de cierta ley empírica que unicaría la dinámica de cuerpos celestes y terrestres sometidos a esta interacción. Esta no es otra que la célebre Ley de Gravitación Universal de Isaac Newton, cuya archiconocida expresión proporciona la fuerza que se ejercen dos cuerpos de masas m1 y m2 separados una distancia r , debido a interacción gravitatoria: F=Gm1m2 r2. (1.1) Aunque de incalculable valor para el desarrollo de la ciencia, tecnología, y comprensión de la naturaleza, proporcionando predicciones sobre el movimiento de los cuerpos del Sistema Solar con gran exactitud, esta ley no está exenta de limitaciones. En el ámbito empírico, se encontraron ciertas discrepancias con las observaciones, siendo la más conocida de ellas el error en la predicción de la precesión del perihelio de Mercurio, para la cual se pueden observar 43 segundos de arco por siglo no predichos por la gravitación newtoniana. El defecto más relevante de la Ley de Gravitación Universal es sin embargo de carácter matemático, y es su incompatibilidad con la relatividad especial, de lo cual Einstein se percató en 1907. Como se puede observar, la ecuación 1.1 no presenta ningún tipo de dependencia con el tiempo. Esto implica que la interacción gravitatoria se transmite de manera instantánea, teniendo que si por ejemplo, pudiésemos variar el valor de una de las masas m1 o m2 , la otra notaría el cambio en la fuerza atractiva que sufre al instante. No sólo entra esto en conicto con la imposibilidad de la transmisión de información a velocidades mayores que la de la luz, si no que además, el cambio en la masa que hemos imaginado y la variación en la fuerza que siente el otro cuerpo deberían ser simultáneos, y sabemos que en relatividad especial la simultaneidad de dos sucesos depende del observador. El mismo problema sufría la ley de Coulomb para la electrostática, que fue resuelto al utilizar las ecuaciones de Maxwell para describir el electromagnetismo, quedando la ley de Coulomb como un caso límite de una de estas ecuaciones. Bajo este marco, Einstein se propone en 1907 desarrollar una teoría dinámica de la interacción 1 1 INTRODUCCIÓN gravitatoria consistente con la relatividad especial tal y como las ecuaciones de Maxwell describían la interacción electromagnética guardando consistencia con esta teoría. En ese mismo año daría con el primer paso hacia la reconciliación de la gravedad y la relatividad especial: el principio de equivalencia; y terminaría culminando en 1915 con el nacimiento de la relatividad general. El principio de equivalencia fue, en en un primer momento, formulado como la idea de que un observador en caída libre (sobre el cual sólo actúa la fuerza gravitatoria) es completamente indistinguible de un observador inercial. Es decir: cualquier experimento realizado por dicho observador debía reproducir los resultados que se esperarían en el marco de la relatividad especial. Sin embargó, Einstein se percató de que bajo un campo gravitatorio no uniforme, se podrían llevar a cabo ciertos experimentos para detectar la presencia de este. Por ejemplo, un observador cayendo radialmente hacia la Tierra podría colocar dos cuerpos separados una cierta distancia tal que, a medida que se adentran en el campo gravitatorio, irían tendiendo a acercarse debido a que siguen trayectorias hacia el centro del planeta (no paralelas). Por lo tanto, el principio de equivalencia fue reformulado para pasar a enunciar que a escala local (es decir, a una escala sucientemente pequeña como para que las inhomogeneidades de los campos gravitatorios presentes se hagan despreciables) un observador en caída libre es completamente indistinguible de uno inercial, y por lo tanto los resultados de los experimentos realizados en su marco de referencia local deben ser acordes con la relatividad especial. Esto sugiere que podemos describir el espacio-tiempo como una unión de regiones sucientemente pequeñas de espacio-tiempo plano (es decir, de Minkowski) pero que da como resultado una unión que globalmente no es plana, sino que presenta curvatura. Einstein interpretó así la interacción gravitatoria como una manifestación de esta curvatura: a grandes rasgos, lo que Einstein expresó en las llamadas ecuaciones de campo de la relatividad general, es cómo la presencia de una determinada distribución de materia-energía induce una curvatura acorde en la geometría del espacio-tiempo, que se propagaría a la velocidad de la luz. Esta curvatura afectaría entonces a las trayectorias que seguirían la luz y las partículas presentes en dicho espacio-tiempo, dando lugar al efecto que clásicamente se ha interpretado como fuerza gravitatoria, o, en palabras de John Archibald Wheeler: El espacio-tiempo le dice a la materia cómo moverse; la materia le dice al espacio-tiempo cómo curvarse. Así, cuando una partícula o fotón no se viese afectado por ninguna otra interacción salvo la gravedad, su trayectoria a lo largo del espacio-tiempo que habita sería el análogo más cercano a una línea recta en el espacio de Minkowski. Estas trayectorias reciben el nombre de geodésicas , concepto que detallaremos en la siguiente sección. Por ahora, basta con hacer énfasis en la importancia de este tipo de trayectorias desde el punto de vista físico: al estar dominada la dinámica de los cuerpos celestes por la interacción gravitatoria, estudiar las geodésicas características de un cierto espacio-tiempo nos permitirá conocer propiedades físicamente relevantes del mismo y producir predicciones teóricas de altísima precisión en caso de que trabajemos con un espacio-tiempo que modele alguna región del que habitamos. La dicultad está en que las geodésicas vienen determinadas por un sistema de cuatro ecuacio2 2 FUNDAMENTO TEÓRICO nes diferenciales de segundo orden (distinto para cada espacio-tiempo) que ademas presentan una gran complejidad formal, lo que diculta encontrar una solución general en la mayoría de los casos. De hecho, suele ser necesario recurrir a diversas herramientas como el cálculo variacional sumado a argumentos de simetría para obtener leyes de conservación que simpliquen el problema [1]. Aquí es donde entra en juego MATLAB. Gracias al software Symbolic Math Toolbox podemos construir con sencillez las ecuaciones diferenciales que determinan las geodésicas de un espacio-tiempo a partir de su tensor métrico, que contiene información sobre la geometría de dicho espacio-tiempo. Una vez planteado el sistema, podemos utilizar una de las herramientas que incluye MATLAB para resolver ecuaciones diferenciales mediante métodos numéricos, obteniendo así una solución numérica a las ecuaciones que describirá trayectorias de, según el caso, luz o partículas en el espacio-tiempo en el que estemos trabajando. Así, a lo largo de este proyecto estudiaremos la ecacia de resolver la ecuación geodésica para distintos espacio-tiempos mediante el método descrito utilizando MATLAB. En primer lugar trataremos de reproducir los tests clásicos más conocidos de la relatividad general, como la predicción de la precesión del perihelio de Mercurio previamente mencionada, y comprobaremos si los resultados obtenidos se ajustan a las predicciones teóricas y observaciones. Posteriormente, estudiaremos las geodésicas características del universo de Gödel, cuyas propiedades peculiares que detallamos en la siguiente sección trataremos de reproducir y visualizar mediante MATLAB. Cabe destacar que durante todo el proyecto se trabajará en unidades naturales, es decir, expresando las velocidades en unidades de la velocidad de la luz en el vacío c , tal que c= 1. (1.2) En caso de que momentáneamente se cambie de sistema de unidades, se expresará explícitamente de antemano. Además, para la obtención de resultados se ha utilizado la versión R2021a de MATLAB. 2. Fundamento teórico En esta sección vamos a realizar una breve introducción a los conceptos teóricos en los que se fundamenta el proyecto. 2.1. Variedad, métrica y tensor métrico En la sección anterior hemos discutido como el principio de equivalencia conduce hacia la descripción del espacio-tiempo como la unión de regiones pequeñas de espacio-tiempo plano que resulta en uno con curvatura. Desde el punto de vista matemático, a un espacio global3 2 FUNDAMENTO TEÓRICO mente curvado que localmente tiende a uno plano se le conoce como variedad ([2], capítulo 5.4), y el espacio plano al que se aproxima en el entorno local de un punto P cualquiera de dicha variedad se le llama espacio tangente en el punto P ([3], capítulo 1.4). Ahora, necesitamos una manera de trabajar con escalares, vectores y tensores en esta variedad. Nuestra teoría física necesitará que sepamos cómo calcular cantidades como distancias entre puntos, módulos de vectores o ángulos formados por dos curvas en dicha variedad. En el espacio euclídeo En esto no entraña dicultad. Por ejemplo, con n= 2 , la distancia ds que separa un punto P1 de coordenadas cartesianas (x1, y1) y otro innitesimalmente cerca P2 de coordenadas (x1+dx, y1+dy) se calcula simplemente como: ds =pdx2+dy2. (2.1) Lo verdaderamente relevante de esta expresión es que la distancia innitesimal entre los puntos P1 y P2 no depende del valor de las coordenadas (x1, y1) en el espacio euclídeo. Esto no se cumple en general para una variedad n-dimensional, y en concreto no lo hace para la variedad tetradimensional con la que identicamos el espacio-tiempo. Por ello, se introduce una estructura llamada tensor métrico gµν , un tensor de orden 2 que en general será función de las coordenadas que utilicemos para describir el espacio-tiempo. Con él, se construye la ecuación métrica , que describe la distancia innitesimal entre puntos, llamada elemento de línea, en función de las coordenadas de la siguiente manera: ds2=gµνdxνdxµ. (2.2) Como claramente dxνdxµ=dxµdxν para cualquier conjunto de coordenadas, gµν debe ser un tensor simétrico. Cabe destacar además que el cuadrado del elemento de línea no depende del conjunto de coordenadas escogido, o lo que es lo mismo: es invariante ante transformaciones de coordenadas, lo que físicamente se traducirá a que ds2 será invariante ante transformaciones entre sistemas de referencia inerciales o no. Esto es un paso conceptualmente gigante con respecto a la relatividad especial, donde el invariante ds2 sólo lo era ante transformaciones entre sistemas de referencia inerciales. Esto es fácilmente demostrable, teniendo que un tensor contravariante (izquierda) y uno covariante (derecha) de orden n se transforman como: T0µ1...µn= i=n Y i=1 ∂x0µi ∂xνi!Tν1...νnT0 µ1...µn= i=n Y i=1 ∂xνi ∂x0µi!Tν1...νn. (2.3) Así, tendremos: ds02=g0 µνdx0µdx0ν=∂xα ∂x0µ ∂xβ ∂x0νgαβ ∂x0µ ∂xαdxα∂x0ν ∂xβdxβ (2.4) 4 2 FUNDAMENTO TEÓRICO y reordenando: ds02=∂xα ∂x0µ ∂x0µ ∂xα ∂xβ ∂x0ν ∂x0ν ∂xβgαβdxµdxν=gαβdxαdxβ=ds2 (2.5) Como ejemplo de algunos tensores métricos tenemos el del espacio euclídeo bidimensional: aµν =1 0 0 1  (2.6) o el del espacio de Minkowski, el cual se puede expresar de dos maneras diferentes según el convenio: ηµν =    −1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1     ηµν =    1 0 0 0 0−1 0 0 0 0 −1 0 0 0 0 −1     . (2.7) En nuestro caso, trabajaremos tomando la primera de las formas expresadas. Por último, una propiedad muy importante del tensor métrico es su uso para generar los llamados tensores asociados ([4], capítulo 5.3) mediante contracción, teniéndose que: Aν1...νn=gν1µ1gν2µ2. . . gνnµnAµ1...µn (2.8) Y también resulta relevante el tensor asociado al propio tensor métrico, cuya derivación podemos encontrar en el capítulo 5.3 de [4] y cumple que: gαµgµν =δν α (2.9) lo cual se traduce, cuando trabajemos con matrices, en que la forma matricial del tensor métrico será la inversa del asociado. 2.2. La ecuación geodésica Como hemos mencionado, la ecuación geodésica es fundamental para el estudio de las trayectorias de los cuerpos en relatividad general. Para obtenerla, debemos profundizar en lo mencionado anteriormente: que las trayectorias que siguen los cuerpos en caída libre son, geométricamente, el camino más recto que puede seguir dicho cuerpo en el espacio-tiempo en el que se encuentra. Intuitivamente, podríamos pensar que lo más parecido a una trayectoria recta en una variedad curvad se corresponde con el camino más corto entre dos puntos, tal y como en el espacio euclídeo el camino más corto entre dos puntos es una línea recta. Sin embargo, en el espacio de Minkowski, por ejemplo, una trayectoria recta entre dos puntos (eventos en el espacio-tiempo) no se corresponde con el camino más corto posible, si no con el más largo ([4], capítulo 6.6). 5 2 FUNDAMENTO TEÓRICO Así, una denición rigurosa de geodésica que cubre ambos casos sería aquella trayectoria que hace la distancia entre dos puntos estacionaria , tal que pequeñas desviaciones de la trayectoria producen un cambio nulo en la longitud S de esta, lo que se expresa como: δS = 0 (2.10) con S=Zds =Zpgµνdxµdxν (2.11) por la ecuación 2.2. Si consideramos ahora que las trayectorias en el espacio-tiempo constituyen curvas parametrizadas por un parámetro afín arbitrario λ tal que xµ=xµ(λ) , utilizando la regla de la cadena obtendremos: S=Zrgµν dxµ dλ dxν dλdλ (2.12) Las ecuaciones 2.10 y 2.12 constituyen un problema de cálculo variacional que, deniendo la función lagrangiana del sistema como L=rgµν dxµ dλ dxν dλ (2.13) y aplicando las ecuaciones de Euler-Lagrange (omitimos el desarrollo, que se encuentra en el capítulo 1.2.5 de [3]) conduce a la ecuación ∂2xν ∂λ2+ Γν αβ ∂xα ∂λ ∂xβ ∂λ = 0 (2.14) donde Γν αβ es la llamada conexión afín, también conocida como símbolos de Christoel de segunda especie, y su expresión es la siguiente: Γν βα =1 2gνρ ∂gαρ ∂xβ+∂gβρ ∂xα−∂gβα ∂xρ=1 2gνρNραβ (2.15) La ecuación 2.14 es la llamada ecuación geodésica, que como se puede observar, es en realidad un sistema de cuatro ecuaciones diferenciales entrelazadas, para ν= 0,1,2,3 . Su solución da la trayectoria xµ(λ) en función del parámetro afín escogido. Para cualquier métrica que describa un espacio-tiempo, existen tres tipos de geodésicas, dependiendo del signo del cuadrado de la distancia innitesimal (invariante ds2 ) que une dos puntos innitesimalmente cercanos de la trayectoria . Así, tendremos: ds2=     −1 para geodésicas temporales. 0 para geodésicas nulas. 1 para geodésicas espaciales. (2.16) 6 2 FUNDAMENTO TEÓRICO El modelo de universo descrito por Gödel es uno con cuatro características clave ([8]): Anisotropía: la materia estaría distribuida en una serie de capas; planos con curvatura negativa perpendiculares a una cierta dirección. Homogeneidad: en cada uno de los planos, la distribución de materia es totalmente homogénea, de manera que el universo estaría compuesto por planos de materia pulverulenta de densidad ρ . Carácter estático: el universo no presenta expansión, por lo que Λ6= 0 . Rotación: el universo en su totalidad presentaría una rotación en la dirección perpendicular a las capas de materia, con una velocidad angular ΩG= 2√πGρ . Figura 3: Representación de dos observadores A y B en el universo de Gödel. Interpretar la rotación en torno a cada uno de ellos con la imagen intuitiva lleva a conclusiones aparentemente paradójicas ([9], capítulo 4.2.5). El hecho de que exista una rotación global del universo parece problemático, pues indicaría que existe un punto privilegiado por el cual pasaría el eje de rotación. Puede sin embargo demostrarse mediante un análisis del modelo Newtoniano de este universo (el cual se encuentra en [8]) que todo observador puede considerarse a sí mismo en el eje de rotación: desde su sistema de referencia, todo el universo giraría en torno a él. La homogeneidad de la distribución de la materia sumado a este hecho lleva a la conclusión de que en el universo de Gödel todos los puntos son equivalentes, lo cual trataremos de comprobar mediante MATLAB. Gödel realizó por primera vez la descripción de este universo en [10], donde propuso la forma de la métrica en coordenadas cartesianas: ds2=a2dx2 0−dx2 1+e2x1/2dx2 2−dx2 3+ 2ex1dx0dx2 (2.43) donde a es una constante con dimensiones de longitud cuyo valor se deduce de imponer que la métrica cumpla las ecuaciones de Einstein ([10]), tal que se llega a: a=1 √8πGρ =1 √2ΩG . (2.44) Sin embargo, para nuestra simulación en MATLAB resulta más práctico trabajar con coordenadas cilíndricas, efectuando primero la transformación a las coordenadas cilíndricas adimensionales que aparece en [10] y posteriormente la que aparece en [11] (sección 3.1, ecuación 16), tal que la forma de la métrica queda (consistente con [6]): ds2=−dt2+dr2 1+[r/(2a)]2+r21−r 2a2dϕ2+dz2−2r2 √2adtdϕ (2.45) 13 2 FUNDAMENTO TEÓRICO 2.5.1. El horizonte óptico del universo de Gödel Realizando un análisis de las geodésicas nulas del universo de Gödel, es posible demostrar que, dada una fuente luminosa situada en el origen de coordenadas de un sistema de referencia, los rayos de luz que abandonen la fuente recorrerán una trayectoria cerrada, llegando hasta un determinado límite y regresando a la fuente. El desarrollo se encuentra en la sección 3.5 de [11]. Así se puede encontrar que para una geodésica nula que parta del origen, con velocidad radial distinta de cero, y contenida en un plano z=cte. , la coordenada radial toma la expresión (en función del parámetro afín de la trayectoria λ ): r(λ)=2a sin 1 2ηλ donde η≡√2u0(0)ΩG, (2.46) siendo u0(0) la componente temporal del cuadrivector velocidad inicial. Así, la máxima distancia que puede alejarse un rayo de luz de su fuente emisora es 2a . Un observador situado a una distancia r > 2a de la fuente, por lo tanto, sería incapaz de verla. Este límite es el llamado radio de Gödel rG= 2a . 2.5.2. Geodésicas circulares De entre las posibles soluciones de la ecuación geodésica para la métrica de Gödel, existe una clase de especial interés: las geodésicas circulares. La solución general a la ecuación geodésica y su particularización para geodésicas circulares se discute en [12], mientras que en [13] se puede encontrar un análisis de los tipos de trayectoria según el momento angular de la misma ( γ ) y las condiciones para obtener geodésicas circulares. Lo que hace llamativas a este tipo de trayectorias es su relación con la existencia de curvas temporales cerradas, uno de los fenómenos más destacados de entre los que se dan en el universo de Gödel. La discusión acerca de estas se puede encontrar en [8], y su análisis queda fuera de los objetivos del proyecto. Lo que sí resulta interesante es que estas geodésicas nos facilitarán el estudio de la homogeneidad del universo de Gödel, como detallaremos en la sección 3.4.4. En concreto, el análisis seguido en [12] proporciona que, con a= 1 , o lo que es lo mismo, trabajando con las distancias en unidades de a una geodésica circular tiene un radio de ([12], ecuación 56): r≤2 sinh 1 2arccosh √2. (2.47) Donde se ha aplicado la transformación de coordenadas que aparece en la ecuación 16 de [11] para pasar de las coordenadas cilíndricas adimensionales de [10] (en las que se trabaja en [12]) a las que se usan para expresar la ecuación métrica 2.45. Además, continuando con el análisis que se realiza en [12] sobre este tipo de trayectorias, se demuestra que al imponer las condiciones estudiadas para que la geodésica sea circular y 14 2 FUNDAMENTO TEÓRICO evaluar ds2 sustituyendo el valor de r por rm´ax = 2 sinh 1 2arccosh √2, (2.48) se obtiene ds2= 0 . Es decir: la geodésica circular de radio rm´ax es una geodésica nula, lo que trataremos de comprobar mediante MATLAB. Por último, también en el capítulo 4.2.6 de [9] y en [8] podemos encontrar una descripción de cómo la inclinación y la forma de los conos de luz de las partículas que viajan en trayectorias circulares va incrementando, lo cual trataremos también de visualizar grácamente en MATLAB. 2.5.3. Homogeneidad del universo de Gödel Figura 4: Inclinación y forma de los conos de luz para observadores en trayectorias de distinto radio. El cilindro sombreado representa aquellas trayectorias que son geodésicas. La última característica que vamos a tratar de reproducir mediante MATLAB es la homogeneidad del universo de Gödel. Esta se puede manifestar de dos maneras. En primer lugar, a través del tensor energía-momento, que describe la distribución de materia/energía planteada anteriormente, que, con la denición del elemento µ, ν del tensor energía-momento (adjunto) que damos en la sección 2.3.1, debe quedar ([9], capítulo 4.2.5) Tµν =    ρ0 0 0 0 0 0 0 0 0 0 0 0 0 0 0     . (2.49) Todo ujo de energía y momento es nulo. La única componente no nula, como vemos, es la densidad de energía T00 que (en unidades naturales) se corresponde con la densidad de la distribución de materia del universo de Gödel. Utilizaremos las herramientas de cálculo simbólico de MATLAB para tratar de reproducr este resultado. En segundo lugar, la estructura de las geodésicas debería reproducir la homogeneidad: aunque jemos un observador en r= 0 , tal que percibe la rotación en torno a él, es de esperar que todo punto r6= 0 pueda actuar como el origen de su propio set de geodésicas, cuya estructura debería ser topológicamente idéntica a aquellas calculadas tomando r= 0 como el origen ([13]). Sin embargo, debido a la dependencia de las componentes del tensor métrico con la coordenada radial r , como se puede ver en la ecuación 2.45, un observador con origen en un punto 15 3 METODOLOGÍA O observará una deformación en las geodésicas con origen en un punto O0 situado a una distancia r , y viceversa ([13]). Por ejemplo, una geodésica circular alrededor de O será vista desde el observador en O0 como una elipse, y una geodésica circular centrada en este último punto será vista desde O también como una elipse, siendo la deformación proporcional a la distancia r que separa los observadores. El análisis que realiza [13] a cerca del tipo de geodésica en función del parámetro de momento angular γ resulta muy útil para entender esto último, ya que es sencillo darse cuenta de que en un sistema de referencia con origen en O0 este parámetro se transformará a γ0 y por lo tanto la forma de la geodésica observada por este, que depende del parámetro de momento angular, debe variar. En MATLAB se llevará a cabo la visualización de este fenómeno construyendo geodésicas circulares y tratando de reproducirlas para orígenes distintos de r= 0 . 3. Metodología A lo largo de esta sección realizará una descripción de los scripts MATLAB construidos y se detallará cómo estos nos conducen a los resultados deseados. Los desarrollos explicativos serán acompañados de las líneas de código más relevantes de los scripts utilizados con el objetivo de ilustrar el procedimiento. 3.1. Construcción de la ecuación geodésica En primer lugar vamos a construir de forma simbólica la ecuación geodésica 2.14 y posteriormente transformar las variables simbólicas en numéricas para poder utilizar ode45 . Para ello necesitaremos primero obtener las componentes de la conexión afín Γν αβ . 3.1.1. Planteamiento y obtención de la conexión afín Comenzamos deniendo las variables simbólicas que vamos a utilizar: las coordenadas y alguna constante que necesitemos visualizar, como el semirradio de Gödel a o la constante de gravitación universal G según el caso. Para ello utilizamos el comando syms seguido de las variables que queremos denir. La dependencia de las coordenadas con el parámetro afín debe expresarse explícitamente. Como el tipo de geodésica a obtener (nula o temporal) se determinará con las condiciones iniciales para la resolución de la ecuación, llamaremos al parámetro s y le daremos la interpretación que corresponda según el caso. Por ejemplo, para trabajar con la métrica de Schwarzschild, el primer paso sería: 1 syms t(s) r(s) theta(s) phi(s) 16 3 METODOLOGÍA A continuación denimos el cuadrivector espacio-temporal. Sin embargo, existe un detalle a tener en cuenta: al depender las coordenadas de un parámetro, estas y las estructuras construidas a partir de ellas van a ser una variable del tipo symfun , y no sym o complex sym . Esto lleva a que, por ejemplo, dado un vector A(s) de tipo symfun , MATLAB interpretará que A(i) es el vector A(s) evaluado en s=i , y no el elemento i de A(s) . Para solventar esto introducimos un elemento intermedio donde especicamos la dependencia con el parámetro explícitamente. Continuando con el ejemplo anterior: 1 Spacetime=[t,r,theta,phi]; 2 x=Spacetime(s); El vector x ahora sí será de tipo sym y nos permitirá trabajar con cualquier elemento i del vector como x(i) . Ahora debemos construir el tensor métrico gµν . Lo más sencillo es usar el comando diag pues construye una matriz diagonal con los elementos especicados. De nuevo debemos utilizar una variable intermedia para que el resultado sea de tipo sym y no symfun . Para el caso de Schwarzschild, recordando que trabajaremos con rs= 1 : 1 metric = diag([-(1-1/r)*1^2,(1-1/r)^(-1),r^(2),r^(2)*sin(theta)^(2)]); 2 gdmatrix=metric(s); Mientras que para la métrica de Gödel, introducimos los elementos no diagonales manualmente (habiendo denido previamente las correspondientes coordenadas) y damos un valor numérico a ag , que se corresponde con el semirradio de Gödel a . Según el caso, también podemos trabajar con esta variable como una simbólica. Sin embargo, salvo para la construcción del tensor energía-momento, escogeremos a= 1 para poder comparar con la bibliografía que trabaja con coordenadas en unidades de este parámetro, como [8]; y así se mantendrá salvo que se indique lo contrario. 1 % Escogemos ag = 1: 2 ag = 1; 3 % Componentes del tensor metrico: 4 gtt = -1; 5 grr = (1+(r/(2*ag))^2)^(-1); 6 gpp = r^2*(1-(r/(2*ag))^2); 7 gzz = 1; 8 gtp = -r^2/(ag*sqrt(2)); 9 metric = diag([gtt,grr,gpp,gzz]); 10 gdmatrix=metric(s); 11 gdmatrix(1,3) = -r(s)^2/(ag*sqrt(2)); 12 gdmatrix(3,1) = gdmatrix(1,3); 17 3 METODOLOGÍA Para construir ahora la conexión afín vamos a comenzar por Nραβ , cuya expresión tal y como aparece en la ecuación 2.15 es: Nραβ =∂gαρ ∂xβ+∂gβρ ∂xα−∂gβα ∂xρ (3.1) Como tenemos tres índices distintos, podemos utilizar tres bucles for que vayan evaluando la expresión componente a componente. Para las derivadas utilizamos el comando diff(f,x, n) , donde f es la expresión a derivar (en nuestro caso el tensor métrico), x es la variable con respecto a la cual derivamos (en nuestro caso la componente correspondiente del cuadrivector xµ ) y n el orden de dicha derivada (que para este caso será uno o directamente puede omitirse pues el orden por defecto es uno). Así, independientemente de la métrica, tendríamos: 1 for p = 1:4 2 for b = 1:4 3 for a = 1:4 4 N(p,a,b) = diff(gdmatrix(a,p),x(b)) + diff(gdmatrix(b,p),x(a)) ... 5 - diff(gdmatrix(b,a),x(p)); 6 end 7 end 8 end Por último, queda efectuar la contracción 1 2gνρNραβ . Para esto resultan de suma utilidad las funciones basicTensor y tensorContraction , desarrolladas en marzo de 2021 por el Profesor Titular de la Universidad de Sevilla Carlos Soria del Hoyo. La primera nos permite construir tensores especicando la matriz o estructura que recoge las componentes del tensor ( compMatrix ) y el carácter (covariante o contravariante) de cada índice. Para ello, se introducen como argumentos de basicTensor la estructura mencionada y un vector la que contiene un +1 en las posiciones de los índices covariantes y un -1 en las de los contravariantes ( idxMask ). Así, para denir 1 2Nραβ y el tensor métrico inverso gµν como variables del tipo basicTensor , utilizando inv para calcular la inversa de la matriz que contiene las componentes de gµν : 1 dw = -1; 2 up = 1; 3 gu = basicTensor(inv(gdmatrix), [up up]); 4 Nt = basicTensor((1/2)*N,[dw dw dw]); La segunda efectúa la contracción de un par de tensores (variables del tipo basicTensor ) indicando con ayuda de una matriz de tres las qué indices se contraen y cuáles permanecen 18 3 METODOLOGÍA sin contraer; identicando cada índice con un valor numérico e introduciendo en la primera de las las los índices que tendrá el tensor nal contraído y en las siguientes los índices de los tensores a contraer. Si uno de los tensores es de menor orden que el resto, se debe introducir un 0 en las columnas correspondientes a los índices que este no tiene. Así, efectuamos la contracción: 1 alfa = 1; 2 beta = 2; 3 nu = 3; 4 rho = 4; 5 NU = 0; 6 F = tensorContraction(Nt,gu,[nu alfa beta;rho alfa beta;nu rho NU]); 7 Gamma = basicTensor(simplify(F.compMatrix),F.idxMask); Tanto F como Gamma son variables del tipo basicTensor que se corresponden con la conexión afín Γν αβ (aunque esta no sea realmente un tensor), pero hemos utilizado el comando simplify para visualizar el resultado de manera más sencilla. Así, ya podemos pasar a la construcción de la ecuación geodésica. 3.1.2. Ecuación geodésica Ahora vamos a implementar otro bucle for que evalue las derivadas cruzadas que aparecen en 2.14 para efectuar la contracción: Γν αβ dxα ds dxβ ds (3.2) por lo que usaremos de nuevo basicTensor para dar a la variable estructura de tensor: 1 for i = 1:4 2 for j = 1:4 3 dadb(i,j) = diff(x(i),s)*diff(x(j),s); 4 end 5 end 6 DADB = basicTensor(dadb,[up up]); Continuaremos efectuando la contracción y creando una nueva variable que contenga la compMatrix (matriz de componentes) del tensor de orden uno resultante. Obtendremos así un vector de cuatro componentes, siendo cada una un miembro de cada una de las componentes de la ecuación geodésica. Cabe destacar que el resultado sigue siendo una variable del tipo sym , pero ahora necesitamos que pase a ser del clase symfun pues posteriormente utilizaremos la función odeToVectorField que requiere una entrada de dicha clase. Para efectuar esta conversión resulta útil la función symfun , que cambia la clase de la variable de 19 3 METODOLOGÍA entrada a la deseada, especicando el parámetro del que dependen las variables simbólicas. Así: 1 GDS = tensorContraction(Gamma,DADB,[nu NU NU;nu alfa beta ;alfa beta NU]); 2 gds = GDS.compMatrix; 3 gds = symfun(gds,s); El otro miembro lo constituyen simplemente las derivadas de segundo orden de las coordenadas: d2xν ds2, (3.3) por lo que aplicando diff con n=2 al cuadrivector espacio-temporal y acoplando con la variable vectorial recién creada, tendremos un vector de salida cuyas componentes son cada una de las cuatro de la ecuación geodésica: 1 eqn = diff(Spacetime,s,s) == -gds; 3.1.3. Implementación en ode45 Una vez construida la ecuación geodésica de manera simbólica, debemos conseguir que ode45 , nuestra herramienta para resolver ecuaciones diferenciales numéricamente, sea capaz de interpretarla. Para ello son esenciales dos funciones. La primera es odeToVectorField , que es la que requería entrada de clase symfun . La utilidad de esta función reside en que transforma sistemas de ecuaciones diferenciales de orden superior a uno en sistemas de primer orden, creando nuevas variables simbólicas para las derivadas de las variables y aumentando el número de ecuaciones consecuentemente. En nuestro caso, teniendo un sistema de cuatro ecuaciones diferenciales de segundo orden, odeToVectorField producirá un sistema de ocho EDOs de primer orden creando cuatro variables nuevas de la forma Dr , Dt , etc., correspondientes con las derivadas primeras con respecto a s de las variables r , t , etc. Así, las cuatro nuevas ecuaciones simplemente indicarán la relación entre las nuevas variables (derivadas) y las originales (coordenadas). Como salida produce, además del nuevo sistema de ecuaciones, un vector al que llamaremos AUX que almacena la totalidad de variables simbólicas (coordenadas más derivadas) en el orden que serán tratadas a partir de ahora, al que debemos prestar atención a la hora de especicar las condiciones iniciales en ode45 y al trabajar con los resultados. La segunda es matlabFunction . Su función es transformar funciones simbólicas (variables clase symfun ) en funciones anónimas de MATLAB, es decir, variables de clase function_handle . Esto es imprescindible, pues ode45 trabaja con entradas de esta última clase. Además de nuestro sistema de ecuaciones simbólico, como entrada especicaremos las variables de la función anónima: un vector que contendrá las ocho que aparecen en las 20 3 METODOLOGÍA ecuaciones (coordenadas y derivadas) al que podemos llamar, por ejemplo, Y y el parámetro del que dependen (para el que hasta ahora hemos utilizado s ). Todo este proceso quedaría como: 1 [V,AUX] = odeToVectorField(eqn); 2 GEO = matlabFunction(V,'vars',{'s','Y'}); Esta última variable, en el caso ejemplo GEO , es la función anónima que será la primera entrada de ode45 . 3.2. Condiciones iniciales En este apartado detallaremos cómo escoger las condiciones iniciales sobre las coordenadas y sus derivadas que requiere ode45 como entrada para así obtener resultados físicamente relevantes. 3.2.1. Tetradas estándar Como ya comentamos en las secciones 1 y 2.1, en el entorno local de un punto de las variedades a través de las cuales describiremos cada espacio-tiempo, el espacio-tiempo toma forma plana o de Minkowski. La base ortonormal que vamos a utilizar del espacio tangente que se corresponde con esta región local de espacio-tiempo plano se denomina tetrada estándar e(i) . Esta se puede construir como una combinación lineal de los cuadrivectores unitarios en dirección de las coordenadas espacio-temporales ∂µ 1 : e(i)=eµ i∂µ (3.4) Las componentes eµ i se pueden obtener a partir de la condición de ortonormalidad [14]: gµνeµ (i)ev (j)=ηij (3.5) con ηij el primero de los expresado en la ecuación 2.7 (siguiendo nuestro convenio). La utilidad de las tetradas estándar es que imponiendo que cumplan la condición de ortonormalidad 3.5, al expresar la cuadrivelocidad en función de sus componentes se cumplen automáticamente la normalización de la cuadrivelocidad , que se expresa: [14]: gµv dxµ dλ dxν dλ =κ (3.6) con κ= 0 para geodésicas nulas y κ=−1 para geodésicas temporales. 1 ∂µ se corresponde con el cuadrivector unitario que apunta en la dirección de la componente µ de xµ . 21 3 METODOLOGÍA Así, sea vµ=dxµ dτ para geodésicas temporales y vµ=dxµ dλ para geodésicas nulas, la cuadrivelocidad inicial se puede expresar como ([14]): v=vµ∂µ=v(i)e(i)=v(0)e(0) +n (3.7) donde se tiene que: n=ψsin χcos ξe(1) + sin χsin ξe(2) + cos χe(3) (3.8) siendo ξ y χ los ángulos que, en esféricas, determinan la dirección inicial del fotón o la partícula, tal y como se puede ver en la gura ?? . El parámetro ψ es el módulo de la parte espacial del cuadrivector velocidad de la partícula o fotón con respecto al origen de coordenadas. Por lo tanto, para geodésicas nulas tendremos ψ= 1 , mientras que para geodésicas temporales ψ=γβ , con: γ=1 √1−v2;β=v∀v < 1 (3.9) donde v es el módulo de la velocidad (ordinaria) de la partícula. La componente v(0) se corresponde, en cada caso, con la componente 0 del cuadrivector velocidad de la partícula /fotón, quedando así que v0=±1 para geodésicas nulas y v0=±γ para temporales. El signo de esta componente determinará si la dirección temporal apunta hacia el futuro ( + ) o hacia el pasado ( − ). En este proyecto tomaremos por defecto el signo + . Figura 5: Representación de la dirección local en función de los ángulos ξ y χ ([14]). Por último, quedaría obtener las componentes de cada tetrada estándar e(i) . Estas se pueden deducir a partir de la ecuación 3.5 aunque omitimos el desarrollo explícito para evitar sobrecargar la memoria del proyecto. Las componentes de las tetradas estándar para las métricas de Schwarzschild y Gödel, que se pueden encontrar en [6], son: Métrica de Schwarzschild: e(t)=1 p1−1/r∂t,e(r)=r1−1 r∂r,e(θ)=1 r∂θ,e(ϕ)=1 rsin θ∂ϕ. (3.10) Métrica de Gödel: e(t)= Γ (∂t+ζ∂ϕ),e(r)=p1+[r/(2a)]2∂r, e(ϕ)= ∆Γ (A∂t+B∂ϕ),e(z)=∂z (3.11) donde: A= + (gtϕ +ζgϕϕ), B =−(gtt +ζgtϕ), Γ = 1 p−(gtt + 2ζgtϕ +ζ2gϕϕ),∆ = 1 qg2 tϕ −gttgϕϕ (3.12) Y la expresión de ζ se detalla en la sección 3.2.3. 22 3 METODOLOGÍA de radio el escogido: vc=rGM r0 . (3.24) Variaremos ligeramente el valor de la velocidad para obtener una órbita elíptica y no esférica. Mientras mayor sea la variación, más excentricidad presentará la órbita, pero debemos tener cuidado de no aumentar demasiado la velocidad u obtendremos órbitas parabólicas o hiperbólicas. Por ello, se ha escogido aumentar un 15% la velocidad vc . El período lo estimaremos aproximando nuestra órbita por una circular: aunque parezca una aproximación muy tosca, no resulta demasiado relevante pues a esta distancia tampoco será una buena aproximación tomar el tiempo propio τ como el tiempo ordinario, y por lo tanto sólo usamos el valor obtenido del período como algo orientativo. Como valor de referencia, podemos tomar 10 períodos. T=2πr0 vc (3.25) Resulta interesante realizar una comparación del valor calculado de la precesión del perihelio frente al que da la expresión 2.35, para comprobar como la validez de esta expresión disminuye a medida que reducimos la distancia al Sol. 3.3.3. Deexión de la trayectoria de la luz Ahora vamos a tratar de reproducir la deexión que provoca el Sol sobre la trayectoria de la luz. Las condiciones iniciales sobre la cuadrivelocidad que debemos escoger ahora vienen determinadas por 3.13. Para calcular la deexión que sufre la trayectoria, lo más sencillo es trazarla desde el punto de máximo acercamiento, es decir, integramos sobre media trayectoria (comparando con la gura 2 que muestra la trayectoria completa). Para ello debemos imponer que la velocidad radial inicial sea nula, por lo que debemos seleccionar cos ξ= 0 . También simplica el problema trabajar en el plano z= 0 , por lo que escogemos tanto θ(λ= 0) = π/2 como cos χ= 0 . El valor inicial de r debe corresponderse con la distancia de máximo acercamiento, b . En nuestro caso, como queremos calcular la deexión que sufre la trayectoria de un fotón al pasar tangente al Sol, debemos introducir r(λ= 0) = b=R que podemos consultar en el apéndice A. Ahora no tenemos ninguna referencia para escoger los límites de integración, pues el parámetro λ con respecto al cual integramos no tiene ningún signicado físico. Sin embargo, podemos tomar como referencia la longitud de la trayectoria que queremos simular para así ir variando el límite de integración en función de la precisión del resultado, ya que es importante que los últimos puntos de la trayectoria simulada se encuentren lo sucientemente lejos del origen como para que no sufran deexión. Por ejemplo, elegimos una integración sobre cien veces la distancia inicial. Además, ahora no requerimos de ningún evento, luego podemos integrar directamente tras especicar una tolerancia adecuada. 29 3 METODOLOGÍA La representación gráca se realiza de manera exactamente igual al caso de la órbita de Mercurio. Para obtener el ángulo de deexión de la trayectoria, debemos calcular el ángulo que forma la dirección inicial de la trayectoria del fotón con la dirección nal, y multiplicando por dos ya que integramos sólo sobre media trayectoria, obtendremos el ángulo deseado. La dirección inicial se corresponde con la de la velocidad inicial, que la calcularemos en cartesianas mediante la relación entre vectores unitarios que encontramos en el apéndice de [5] a partir de los valores de las componentes r , θ y ϕ de la cuadrivelocidad inicial. Para la dirección nal, suponiendo que hayamos integrado sobre un intervalo sucientemente grande como para que la trayectoria nal se aproxime lo máximo posible a una línea recta, puede obtenerse calculando el vector que une los puntos (en cartesianas) último y penúltimo de la trayectoria integrada. 3.3.4. Precisión del cálculo de la deexión en función de la distancia de máximo acercamiento A continuación, vamos a estudiar la validez de la expresión aproximada 2.41 comparando el ángulo de deexión obtenido al utilizar dicha expresión y el calculado mediante MATLAB para un rango de valores de la distancia b de máximo acercamiento. En concreto, se ha estudiado la relación entre ángulo teórico-aproximado y calculado para valores de b entre 20 rs y 2000 rs . Para ello, se ha implementado un bucle for dentro del cual se introduce la integración mediante ode45 y el cálculo del ángulo de deexión, para que se realice sucesivamente para cada valor de b . Almacenando en una variable vectorial los valores del ángulo de deexión calculados y en otra los dados por 2.41, podremos posteriormente compararlos grácamente. También, podemos realizar una representación gráca del ángulo calculado frente al teóricoaproximado, la cual sólo realizaremos hasta b= 100rs , pues la diferencia se hace muy poco apreciable para distancias mayores; y una representación de la razón entre ambos para comprobar que tiende asintóticamente a uno a medida que aumenta la distancia b , como esperaríamos dada la condición 2.34. Como complemento de este ejercicio, es interesante representar grácamente algunas trayectorias de la luz para distancias de máximo acercamiento b de órdenes cercanos al radio de Schwarzschild ([16] sección 4.1.1). En lugar del Sol, es más representativo representar un agujero negro de Schwarzschild con la masa del Sol como el responsable de la interacción gravitatoria. 3.4. Obtención de resultados: universo de Gödel 30 3 METODOLOGÍA 3.4.1. Horizonte óptico del universo de Gödel Para comprobar el límite de visión de un observador en un punto cualquiera del universo de Gödel, debemos construir geodésicas nulas que partan del origen en direcciones arbitrarias. Para ello establecemos unas condiciones iniciales adecuadas siguiendo la expresión dada en 3.16. Para simplicar, escogeremos trabajar en el LNRF , tal que: ζ=ω=−gtϕ gϕϕ (3.26) Para facilitar la tarea de establecer las condiciones iniciales, vamos a calcular por separado el valor inicial de A , B , ∆ y Γ . Para ello, damos un valor inicial a la coordenada r , de la que dependen las componentes de gµν y estas variables, escribimos su expresión en función de las componentes del tensor métrico siguiendo 3.12 y utilizamos subs para sustituir el valor inicial escogido y en las expresiones simbólicas que dependan de r . Además, utilizamos la función double para transformar estas variables de simbólicas a numéricas una vez realizada la sustitución. 1 w = -gtp/gpp; 2 F = 1/sqrt(-(gtt+2*w*gtp + w^2*gpp)); 3 A = gtp + w*gpp; 4 B = gtt + w*gtp; 5 Dl = 1/sqrt(gtp^2-gtt*gpp); 6 % Efectuamos la sustitucion, por ejemplo: 7 w = double(subs(w,r,rini)); 8 % Y asi sucesivamente con A, B y Dl Sin embargo, debemos tener en cuenta que las coordenadas cilíndricas no están denidas en el origen, luego no podemos establecer r = 0 como punto inicial. Sin embargo, resulta equivalente escoger un punto muy cercano al origen, tal que r(λ= 0) << 1 y así evitamos el problema. Por ejemplo, podríamos escoger r(λ= 0) = 10−6 . Ahora continuamos deniendo el resto de condiciones iniciales: es muy importante que la trayectoria esté contenida en un plano z= cte . ya que el universo de Gödel no es isótropo, por lo que debemos imponer cos χ= 0 . En cuanto a ξ , podemos escoger, por ejemplo, ξ= 0 para simular las trayectorias de un fotón que parte en dirección radial desde el origen (aunque la única imposición es que la velocidad radial sea no nula). Así ya podemos denir las componentes de la cuadrivelocidad como lo indica 3.16. Con el objetivo de visualizar la simetría rotacional del universo de Gödel alrededor del eje que constituye el observador que se encuentra en el origen de coordenadas, puede resultar interesante representar varias trayectorias simultáneamente para distintos valores de ϕ . Para ello utilizaremos un bucle for deniendo en un vector un conjunto de valores iniciales distintos para la ϕ . Además de la trayectoria podemos también representar grácamente el radio de Gödel, como una circunferencia de radio rG= 2a 31 3 METODOLOGÍA Por último quedaría comprobar si la trayectoria de la luz se ajusta a la ecuación 2.46. Para ello, debemos identicar el valor del parámetro afín s para el cual la coordenada radial se hace cero, que debe quedar almacenado en el vector sv , y comprobar si efectivamente cumple las ecuaciones mencionadas. El valor de u0(0) se corresponderá con el valor inicial que tome la componente temporal de la cuadrivelocidad. 3.4.2. Geodésicas circulares en el universo de Gödel Para simular geodésicas circulares en la métrica de Gödel, tal y como están descritas en la sección 2.5.2, podemos hacer uso de la herramienta dsolve de MATLAB que es capaz de resolver simbólicamente ecuaciones diferenciales. Esto es debido a que dada la libertad de elección en las condiciones iniciales, debemos encontrar aquellas que han de cumplirse tal que la geodésica resultante sea una circunferencia en torno al origen. Para ello, ya que necesitamos que el valor de r sea constante, debemos realizar el procedimiento especicado hasta la construcción de la ecuación geodésica, aunque con dos diferencias: en primer lugar, eliminamos la dependencia de r y z con el parámetro escogido s , ya que buscamos geodésicas circulares y además contenidas en un plano z= cte. Si comprobamos la forma explícita de la ecuación geodésica, se puede ver que ha pasado a ser el siguiente sistema de ecuaciones diferenciales: d2t ds2= 0 ,d2z ds2= 0 ,d2ϕ ds2= 0 1 8rdϕ ds2 r4+ 2r2−8+√2rr2 4+ 1dϕ ds dt ds= 0. (3.27) Este es un sistema con más ecuaciones que incógnitas. Despreciando la segunda de ellas pues el valor de z resulta irrelevante, quedan tres ecuaciones de las cuales debemos escoger dos, y el resultado no debe depender de dicha selección. Claramente, la ecuación en la segunda línea será siempre físicamente relevante pues es la única que contiene la coordenada r , por lo que debemos trabajar con ella. Ahora, tanto si escogemos la primera de las ecuaciones como si elegimos la tercera se llega a la misma condición tal y como era esperado. Esto se demuestra mediante el script MATLAB adjunto en el apéndice B. Por ahora, nos restringimos a resolver el sistema formado por la segunda y tercera ecuación de 3.27. Almacenando ambas ecuaciones en una variable vectorial, siendo cada ecuación un elemento del vector, y utilizándola como argumento de entrada de dsolve , se obtiene la solución general, que es: ϕ=C2+C1·s , t =C3+√2s·(2C1−C1r2) 4 (3.28) donde C1 , C2 y C3 son constantes a determinar a partir de las condiciones iniciales. En primer lugar, con t(0) = 0 y escogiendo ϕ(0) = 0 , se tiene que: C3=C2= 0 (3.29) 32 3 METODOLOGÍA Ahora, suponiendo que las geodésicas circulares que vamos a obtener son geodésicas temporales, podemos escoger β= 0 , lo cual simplica mucho la imposición de condiciones iniciales sobre las derivadas primeras de las variables. Así, el parámetro s pasará a identicarse con τ . De la ecuación 3.17 con β= 0 , se obtiene que: dϕ dτ(τ= 0) = Γζ;dt dτ(τ= 0) = Γ (3.30) dado que con β= 0 , γ= 1 . Ahora, imponiendo estas condiciones en nuestro sistema, se llega a: Γζ=C1,Γ = √2·(2C1−C1r2) 4 (3.31) donde podemos introducir el valor de C1 dado por la primera ecuación en la segunda, para despejar ζ y obtener: ζ=2√2 2−r2. (3.32) Es decir: toda geodésica temporal que cumpla la condición anterior será una geodésica circular con radio r , siempre que el valor del radio esté limitado tal y como marca la ecuación 2.47. Así, donde antes habíamos escogido ζ=ω para trabajar en el LNRF , ahora la deniremos explícitamente de acuerdo con la expresión 3.32. En cuanto a las condiciones iniciales para la cuadrivelocidad, debemos guardar coherencia con lo impuesto anteriormente: β= 0 . Los valores de χ y ξ son irrelevantes pues al imponer β= 0 todos los términos que dependen de los ángulos se hacen nulos. Para la geodésica nula circular de radio r=rm´ax , podemos basarnos en el desarrollo que se realiza en la sección 2.2 de [11], en la que, adaptado a nuestra notación, se deduce que debemos escoger para construir la tetrada estándar correspondiente: ζ=ω±rω2−gtt gϕϕ (3.33) con ω la denida en la ecuación 3.19. Es decir, debemos escoger, o bien ζm´ax o bien ζm´ın dependiendo del sentido en que busquemos que se recorra la circunferencia. Sin embargo, existe un inconveniente: se puede comprobar que el parámetro Γ no está denido para los valores extremos de ζ . Si evaluemos el denominador al sustituir ζm´ax o ζm´ın en su expresión (3.12): gtt + 2 ω±rω2−gtt gϕϕ gtϕ +ω±rω2−gtt gϕϕ 2 gϕϕ = =gtt −2g2 tϕ gϕϕ ±2gtϕsg2 tϕ g2 ϕϕ −gtt gϕϕ +2 g2 tϕ gϕϕ −gtt ∓2gtϕsg2 tϕ g2 ϕϕ −gtt gϕϕ = 0 (3.34) 33 3 METODOLOGÍA todos los términos se anulan y por lo tanto no podremos expresar las condiciones iniciales en función de la tetrada estándar con ζ=ζm´ax o ζm´ın . Sin embargo, podemos aprovechar que estamos resolviendo la ecuación geodésica de manera numérica por lo que, para una visualización gráca, será equivalente construir una geodésica que se aproxime a la teórica. Para ello, podemos, por ejemplo, escoger: ζ=ω±0.999 rω2−gtt gϕϕ (3.35) y así Γ continuará estando denido permitiéndonos establecer condiciones iniciales sobre la cuadrivelocidad utilizando la tetrada estándar. Así, podemos escoger ξ=χ=π/2 para anular las componentes vertical y radial de la velocidad, y como radio inicial el valor de rm´ax dado en 2.48. Podemos utilizar hold on para visualizar al mismo tiempo geodésicas circulares nulas y temporales si utilizamos distintos scripts. Por último, compararemos la orientación de los conos de luz para las partículas que recorren las geodésicas circulares obtenidas. Para ello, representaremos, en lugar de las coordenadas x , y y z , las coordenadas x , y y el tiempo t en el eje vertical. Para representar los conos una opción es utilizar mesh ; centrándolos en el punto inicial de cada trayectoria. 3.4.3. Tensor energía-momento del universo de Gödel Para este apartado trabajaremos exclusivamente con herramientas de cálculo simbólico. Hasta la obtención de la conexión afín, el procedimiento es exactamente el mismo que en los casos anteriores con la métrica de Gödel salvo por que esta vez denimos el semirradio de Gödel a como una variable simbólica, llamándola, por ejemplo, de nuevo ag . También es conveniente denir como variable simbólica la constante de gravitación universal G , que hemos llamado Gn . A partir de la conexión afín, debemos construir el tensor de Riemann. Para ello, creamos una nueva variable que almacene las componentes de la conexión afín en forma matricial, e implementamos un bucle for que evalúe las derivadas espaciales de la ecuación 2.24. Para las contracciones, podemos construir otro bucle for expresando de forma explícita la suma sobre el índice repetido, ya que resulta complicado utilizar tensorContraction cuando la salida es un tensor de orden 4. 1 % Componentes de la conexion afin: 2 CR = Gamma.compMatrix; 3 % Derivadas: 4 for a2 = 1:4 5 for b2 = 1:4 6 for v2 = 1:4 7 for p2 = 1:4 8 R12(a2,b2,v2,p2) = diff(CR(a2,b2,p2),x(v2))- diff(CR(a2,b2,v2),x(p2)); 34 3 METODOLOGÍA 9 end 10 end 11 end 12 end 13 % Contracciones: 14 for a2 = 1:4 15 for b2 = 1:4 16 for v2 = 1:4 17 for p2 = 1:4 18 R3(a2,b2,v2,p2) = CR(a2,1,v2)*CR(1,b2,p2) + CR(a2,2,v2)*CR(2,b2,p2) + ... 19 CR(a2,3,v2)*CR(3,b2,p2) + CR(a2,4,v2)*CR(4,b2,p2); 20 R4(a2,b2,v2,p2) = CR(a2,1,p2)*CR(1,b2,v2) + CR(a2,2,p2)*CR(2,b2,v2) + ... 21 CR(a2,3,p2)*CR(3,b2,v2) + CR(a2,4,p2)*CR(4,b2,v2); 22 end 23 end 24 end 25 end 26 % Creamos el tensor: 27 Rtens = basicTensor(simplify(R12 + R3 - R4),[up dw dw dw]) A continuación, construimos el tensor asociado covariante de Riemann contrayendo con gµν mediante tensorContraction . Contrayendo ahora con gµν ( gu ) llegamos al tensor de Ricci (2.26) ( Ricci ), y realizando una última contracción con gµν obtenemos el escalar de Ricci (2.28) ( Resc ). Ya tenemos casi todo lo necesario para construir el tensor energía-momento. Sólo falta la constante cosmológica, que podemos encontrar en [6] (introducimos un − manualmente pues debemos tener en cuenta que en [6] se utiliza un convenio distinto al utilizado en este proyecto tal que Λ aparece en la ecuación 2.29 con el signo cambiado): Λ = −R 2. (3.36) El siguiente paso sería, teniendo todos los términos calculados de la ecuación 2.27, despejar Tµν de la ecuación 2.29. 1 % Constante cosmologica: 2 Cosmo = -Resc/2; 3 % Componentes del tensor T_\mu \nu: 4 StEn = (1/(8*sym(pi)*Gn))*simplify(Ricci - 0.5*gdmatrix*Resc - Cosmo* gdmatrix); 5 % Damos estructura de tensor 6 Ttens = basicTensor(StEn,[dw dw]); 35 3 METODOLOGÍA Por último, quedaría realizar la contracción correspondiente para pasar de Tµν a Tµν : Tµν =gµαgνβTαβ (3.37) para lo cual utilizaremos de nuevo tensorContraction dos veces con gu , nuestro tensor métrico asociado y ya podemos recuperamos la matriz de componentes del tensor resultante para visualizar el resultado. 3.4.4. Homogeneidad del universo de Gödel a través de las geodésicas Como estudio nal de las propiedades del universo de Gödel, nos interesa tratar de obtener geodésicas en torno a un origen distinto a r= 0 , tal y como hemos descrito en la sección 2.5.3. Para simplicar la tarea, vamos a tratar de reobtener geodésicas circulares temporales en torno a otro origen, ya que en primer lugar son más útiles para visualizar la localización del nuevo origen y en segundo lugar resulta más sencillo inferir las condiciones que debemos imponer. Recordemos que, por lo que discutimos en la sección 2.5.3, al trabajar con el origen de coordenadas jado en r= 0 , en realidad estamos trabajando desde el sistema de referencia de un observador en r= 0 (punto O ), por lo que no deberíamos observar geodésicas circulares si no elípticas con excentricidad mayor a medida que se alejen del observador en el origen. Supongamos que buscamos obtener geodésicas temporales en torno a un origen situado sobre el eje x a una distancia r0 de O . El nuevo origen O0 tendrá como coordenadas r=r0 y ϕ= 0 al estar sobre el eje x . Sea ahora % el radio de la circunferencia centrada en O0 , realicemos una analogía con las geodésicas circulares obtenidas en torno a O . Escogiendo β= 0 , las ecuaciones 3.28 y 3.31, junto a la condición inicial ϕ(0) = 0 , llevan a ϕ= Γζ·s. (3.38) Eso quiere decir que una geodésica circular en torno a O0 (recordemos, circular para dicho observador ), podría tener una coordenada ϕ0 dada por: ϕ0= Γ0ζ0·s (3.39) donde ζ0 debería tomar forma análoga a 3.32 pero para el radio % seleccionado: ζ0=2√2 2−%2 (3.40) y la denición de Γ0 también guardaría analogía con la de Γ . Sustituyendo en 3.12 las componentes de gµν , para O0 se tendría: Γ0=1 q1 + ζ0%2√2/a −ζ02%2(1 −[%/(2a)]2) . (3.41) 36 4 RESULTADOS Queda ahora preguntarnos cómo relacionar nuestras coordenadas del sistema O con las de O0 . El desarrollo se encuentra adjunto en la gura del apéndice apéndice C. Se puede comprobar de manera sencilla que: ϕ= atan %sin ϕ0 r0+%cos ϕ0, r =r0+%cos ϕ0 cos ϕ. (3.42) Para hallar las condiciones iniciales a imponer, sólo queda derivar y evaluar en s= 0 , donde ϕ0(0) = 0 dada la ecuación 3.40. Esto se ha realizado mediante cálculo simbólico utilizando un script MATLAB que se adjunta en el apéndice B. dϕ ds(s= 0) = %ζ0Γ0 %+r0 ,dr ds(s= 0) = 0. (3.43) Por último, debemos establecer las condiciones iniciales adecuadas para la componente temporal de la cuadrivelocidad tal que se cumpla 3.6 con κ=−1 . Ya tenemos las condiciones a imponer sobre el resto de componentes de la cuadrivelocidad, luego planteando la ecuación 3.6 y evaluando las componentes del tensor métrico para el valor inicial del radio, podemos despejar la condición inicial sobre d t/ d s . Como ϕ0(0) = 0 , sustituyendo en 3.42 llegamos a: ϕ(s= 0) = 0 , r (s= 0) = r0+%. (3.44) Para despejar el valor inicial de d t/ d s se la ha denido como variable simbólica en MATLAB se ha utilizado la función solve para despejarla de 3.6. Adjuntamos en el apéndice B el script mediante el cual se realiza esta operación. Obtenemos dos valores posibles, tal y como al utilizar la tetrada estándar teníamos libertad de elección en el ± que acompañaba a a la componente temporal (3.17), aunque la elección sin embargo no inuirá en la visualización de la geodésica que queremos obtener. Denotando dϕ ds(s= 0) := ˙ϕ0 dt ds(s= 0) = ±q4 ˙ϕ2 0(r0+%)2+ ˙ϕ2 0(r0+%)4+ 4 2−√2 ˙ϕ0(r0+%)2 2 (3.45) Así, escogemos uno de ellos arbitrariamente y ya podemos integrar con ode45 y representar grácamente. 4. Resultados 4.1. Tests de relatividad general 4.1.1. Precesión del perihelio de Mercurio En este apartado presentaremos los resultados obtenidos para la simulación de la órbita de Mercurio y el cálculo de la precesión del perihelio. En la tabla 1 se presentan los distintos 37 4 RESULTADOS resultados para el cálculo de la precesión del perihelio frente al valor escogido en cada caso de 'RelTol' , siendo 'RelTol' = 2.22045 ·10−14 el mínimo que soporta MATLAB para esta integración. El error relativo para cada caso se ha calculado tomando como referencia el valor obtenido mediante la expresión teórico-aproximada que calculamos en 2.36. Posteriormente, en las guras presentamos las guras 6a y 6b obtenidas al representar grá- camente la trayectoria. Tolerancia relativa 'RelTol' Precesión del perihelio ( 00 /siglo) Error relativo 1·10−7 401.445992 8.34 ·102 % 1·10−8 77.617655 80.58% 1·10−9 48.182753 12.10% 1·10−10 43.235970 0.59% 1·10−11 42.998856 3,64 ·10−2 % 1·10−12 42.984375 2.67 ·10−3 % 1·10−13 42.982699 1.23 ·10−3 % 2.22045 ·10−14 42.982579 1.51 ·10−3 % Cuadro 1: Resultados de la precesión del perihelio de Mercurio para distintos valores de tolerancia relativa. (a) Vista cenital de la trayectoria. (b) Vista tridimensional de la trayectoria. Figura 6: Trayectoria de Mercurio a lo largo de cien órbitas. 4.1.2. Precesión del perihelio en una órbita cercana al Sol Se han simulado tres órbitas para tres valores distintos del radio inicial r0 : 100rs , 60rs y 20rs . A continuación presentamos la representación gráca de la órbita para cada caso y una tabla 38 5 DISCUSIÓN DE RESULTADOS 5. Discusión de resultados 5.1. Tests de relatividad general 5.1.1. Precesión del perihelio de Mercurio Los valores presentados en la tabla 1 dejan claro que existe un umbral de tolerancia relativa por encima del cual hemos de trabajar para obtener resultados precisos en el estudio de la precesión del perihelio de Mercurio. En concreto, los errores relativos sugieren que debemos trabajar con valores de 'RelTol' inferiores o iguales a 1·10−11 . La pérdida de precisión al aumentar la tolerancia no es lineal, pues podemos observar cómo se produce un salto brusco entre los valores de la precesión para los órdenes de magnitud más grandes, con un error relativo del 834% para 'RelTol' = 1·10−7 . Es cierto que el error relativo obtenido con, por ejemplo, 'RelTol' = 1·10−10 podría parecer lo sucientemente bajo, pero debemos tener en cuenta que en este proyecto estamos realizando simulaciones numéricas. Un error relativo del 0.56% resulta aceptable para una medida experimental, pero para un resultado producto de una simulación numérica debemos exigir más precisión. Así, se comprueba que a partir de 'RelTol' = 1·10−11 el valor de la precesión obtenido mediante MATLAB se ajusta con gran precisión al predicho a través de la expresión 2.35: el método utilizado resulta satisfactorio. Además, también encontramos que la precisión del cálculo crece en gran medida para los valores más exigentes de 'RelTol' . Sin embargo para el valor mínimo obtenemos un resultado de la precesión ligeramente menos preciso que para 'RelTol' = 1·10−13 : esto podría indicar que existe alguna otra fuente de error numérico tal que no podríamos llegar a una precisión innita con una tolerancia arbitrariamente baja. Sin embargo, resulta irrelevante pues la precisión obtenida es suciente para satisfacer los objetivos del proyecto. 5.1.2. Precesión del perihelio en una órbita cercana al Sol Tal y como esperábamos, las representaciones grácas mostradas en las guras 7 y 8 y los resultados recogidos en la tabla 2 revelan que la precesión se hace más acusada cuanto más cercana es la órbita al Sol. Esto es debido a la dependencia con el inverso del semieje mayor que encontramos en 2.35. Otra observación relevante es la pérdida de validez de la expresión 2.35. Dada la condición 2.34, era esperado que para órbitas cuyo perihelio se sitúe a distancias de órdenes de magnitud cercanos a rs 2.35 no devuelva una predicción precisa de la precesión. Aun así, para el rango de distancias estudiado, sigue siendo válida para obtener una estimación vaga de la precesión, 45 5 DISCUSIÓN DE RESULTADOS si por ejemplo sólo nos interesase conocer el orden de magnitud de este dato. 5.1.3. Deexión de la trayectoria de la luz La gura 2 apenas muestra desviación en la trayectoria de la luz, como era de esperar. Teniendo un ángulo de deexión de un orden de magnitud tan bajo (2.42), para que fuese apreciable en una representación gráca tendríamos que tomar una escala inmensamente grande. En la tabla 3 se puede apreciar como en esta ocasión el rango de tolerancias con el que podemos trabajar no es tan exigente como en el caso de la precesión del perihelio de Mercurio. Basta con que 'RelTol' ≤1·10−8 para obtener resultados con una precisión suciente como para que considerarlos satisfactorios. A partir de una tolerancia de 1·10−10 , los resultados se ajustan con una exactitud importante al valor obtenido en 2.42, con errores relativos del orden de 10−4% . 5.1.4. Precisión del cálculo de la deexión en función de la distancia de máximo acercamiento En primer lugar, en la gura 10a se puede observar como los valores del ángulo de deexión calculados mediante MATLAB y utilizando 2.41 divergen a medida que la distancia de máximo acercamiento disminuye. Sin embargo, la gura 10b resulta más práctica para visualizar la discrepancia a medida que disminuye el radio: se observa como la razón entre ángulos de deexión se aleja exponencialmente de la unidad al trabajar con distancias de máximo acercamiento b menores a 200rs aproximadamente. De nuevo, este resultado es el esperado al ser válida la expresión 2.41 en el rango de valores que cumplen 2.34. Por último, de las trayectorias representadas en la gura 11 la más interesante resulta la correspondiente a b= (3/2)rs . Esta trayectoria se corresponde con la llamada esfera de fotones , que como se demuestra en el ejemplo 3.1 de [3], es la órbita circular que sigue la luz cuando parte en dirección perpendicular a la radial de una distancia b=3 2rs del cuerpo masivo que orbita, en el caso de la gura un agujero negro de Schwarzschild. 5.2. Universo de Gödel 5.2.1. Horizonte óptico En primera instancia, en la representación gráca de la gura 12a se observa el fenómeno esperado y discutido en la sección 2.5.1: las trayectorias que sigue un rayo de luz emitido desde el origen de coordenadas son cerradas y están delimitadas por el radio de Gödel de 46 5 DISCUSIÓN DE RESULTADOS valor rg= 2a . Además, el análisis de la periodicidad de la trayectoria de la luz también lleva al resultado esperado: se cumple que el valor del parámetro λ , en nuestro código llamado s , es tal que la fase del seno de la ecuación 2.46 es 3.1416, una muy buena aproximación a π que sería el valor esperado tal que r= 0 según 2.46. Por último, la gura 12b reproduce con alta delidad una representación gráca de la ecuación 2.46. Utilizando plot para representar la función f(λ) = 2  sin λ 2 (5.1) donde se ha tenido en cuenta que con a= 1 , η= 1 al tener también u0 0= 1 como se comprobó en la sección 4.2.1, se obtiene lo representado en la gura 15, que guarda una forma prácticamente idéntica con la representación de la gura 12b, como cabía esperar. Figura 15: Representación gráca de la función expresada en la ecuación 5.1 La implicación física de este límite en las trayectorias de la luz es que un observador no podrá comunicarse con nada que se encuentre a una distancia mayor a rG= 2a , dado que las señales electromagnéticas que intentase mandarle regresarían a él sin avanzar más allá de esta distancia. El observador entonces sólo podría guardar relación causal con cualquier evento que sucediera dentro de dicho radio. 5.2.2. Geodésicas circulares La existencia de geodésicas circulares resulta más útil para el análisis de la sección 5.2.4, pero resulta interesante comprobar como a pesar de tener que realizar la aproximación descrita en la ecuación 3.35 la geodésica nula obtenida, representada en la gura 13a es, al menos a la escala representada, indistinguible de la que se describe teóricamente. Por otra parte, la gura 13b revela la discutida por [8] y [9] inclinación de los conos de luz. Cabe mencionar que además de variar en orientación, los conos de luz de una partícula siguiendo una geodésica circular también presentan un ángulo de apertura variable con el radio de la geodésica como se discute en [8]. Sin embargo, como el análisis en profundidad de la forma de los conos de luz en el universo de Gödel y sus implicaciones con la causalidad y curvas temporales cerradas escapan a los objetivos del proyecto, se ha optado por representar sólo la inclinación de estos a modo de visualización cualitativa del efecto. Por último, es interesante observar como la línea de universo de una partícula u observador situado en el origen es una línea recta vertical en el eje temporal: se corresponde con la de un observador estático (en el sistema de referencia escogido) en relatividad especial. 47 5 DISCUSIÓN DE RESULTADOS 5.2.3. Cálculo del Tensor energía-momento El tensor energía-momento coincide con el dado en la ecuación 2.49, como cabía esperar y su relación con la homogeneidad del universo de Gödel ya fue discutida en la sección 2.5.3. El aspecto más destacable de este resultado es que prueba la robustez de las herramientas de cálculo simbólico de MATLAB. La construcción del tensor de Riemann no es un procedimiento sencillo, como se detalla en la sección 3.4.3, ya que físicamente se corresponde con un tensor de orden cuatro. Asimismo, las herramientas de manejo y contracción de tensores, basicTensor y tensorContraction desarrolladas por el profesor Carlos Soria del Hoyo cumplen a la perfección, dada la cantidad de contracciones que se realizan para la construcción del tensor energía-momento. Para proyectos futuros que involucren álgebra tensorial probablemente resulten de gran ayuda. 5.2.4. Homogeneidad del universo de Gödel a través de las geodésicas En la gura 13 se aprecia el efecto discutido en la sección 2.5.3: a medida que los observadores se alejan de O , las formas de las geodésicas que percibe este último desde su sistema de referencia se alejan más de las circunferencias que perciben los observadores A , B , C y D en su respectivo sistema de referencia. Así, mientras que la geodésica centrada en el observador A se percibe casi como una circunferenica, la centrada en D adopta una forma de elipse con alta excentricidad vista desde O . Es interesante notar cómo para O los observadores se sitúan en uno de los focos de las elipses percibidas por este. De entre los cuatro observadores alejados de O , el más interesante resulta C . Esto es debido a que parte de la geodésica circular centrada en C cruza el radio de Gödel de O . Imaginemos así que el observador C decide poner en la órbita representada una antena receptora. Si O decidiera comunicarse con él mediante señales electromagnéticas, debería calcular cuándo la antena va a cruzar su radio de Gödel. De no encontrarse dentro de este, la señal jamás podría llegarle por lo discutido en la sección 5.2.1. 48 5 DISCUSIÓN DE RESULTADOS I1 I2 Figura 16: Hipotética comunicación entre O y C . Para el caso ejemplicado en la gura 16, si el observador O emite una señal luminosa que siga una de las geodésicas discutidas en 2.5.1, la señal debe ser emitida tal que se cruce con la antena en las intersecciones I1 o I2 para que esta sea recibida por la antena del observador C ([9]). Podríamos ahora preguntarnos qué ocurre desde el punto de vista de C . Claramente, la geodésica en la que orbita la antena para él es una circunferencia. ¾Cómo es posible entonces que cruce el radio de Gödel de O ? La respuesta a esta pregunta yace en el hecho de que claramente, al ser todos los puntos equivalentes en el universo de Gödel, cualquier observador puede delinear a su alrededor un radio de Gödel de rG= 2a . Sin embargo, al igual que ocurre con las geodésicas, la forma del radio de Gödel de cada observador no se mantendrá vista desde el sistema de referencia de otro. La discusión acerca de este hecho se encuentra en el capítulo 5.3.2 de [9]. Figura 17: Radios de Gödel al rededor de tres observadores O , A y B vistos desde el sistema de referencia de O ([9]). Así, para C , será el radio de Gödel de O el que con su deformación se interne en la geodésica circular que para él sigue la antena, tal y como ocurre en la gura 17, donde el observador 49 6 CONCLUSIONES allí llamado O percibe un radio de Gödel deformado para los observadores A y B . 6. Conclusiones La relatividad general presenta una diferencia clave con otras disciplinas de la física como pueden ser la mecánica cuántica o la electrodinámica: el número de experimentos que podemos reproducir en un laboratorio es muy limitado. Exceptuando experimentos como los relacionados con el redshift gravitacional ([4], capítulo 2), la mayoría de estudios empíricos relacionados con la relatividad general deben realizarse a través de observaciones (ya sea mediante instrumentos ópticos o de otra clase como los detectores de ondas gravitacionales). Así, en un panorama cientíco en el que la simulación cobra cada vez más importancia de acuerdo al desarrollo de las tecnologías adecuadas, la posibilidad de realizar estudios sobre relatividad general a través de simulaciones supone una oportunidad única para explorar esta teoría de maneras que antes del desarrollo de la computación habrían resultado imposibles. En nuestro caso, en vista de los resultados obtenidos tanto para los tests de relatividad general como a la hora de reproducir las características del universo de Gödel, podemos conrmar que MATLAB es una herramienta muy a tener en cuenta a la hora de realizar estudios relacionados con la relatividad general. Como se comentaba en la sección 1, la ecuación geodésica presentan una complejidad formal tan grande que, salvo mediante la utilización de técnicas de simplicación, encontrar una solución general puede suponer un esfuerzo inmenso. Hemos comprobado que ode45 , pese a no dar una solución analítica a dichas ecuaciones si no numérica, cumple con creces las expectativas a la hora de resolver estas ecuaciones. No sólo resulta clave el hecho de que mediante integración en ode45 podamos reproducir predicciones teóricas con elevadísima precisión, como el caso de los tests de relatividad general; también es fundamental el potencial de Symbolic Math Toolbox . Sin esta herramienta, habríamos tenido primero que calcular de forma manual las componentes de la ecuación geodésica para cada métrica y después implementarlas como function_handle en ode45 . Evitar esto ahorra, además de tiempo y esfuerzo, posibilad de cometer errores durante el procedimiento, y esta no es su única utilidad. Las herramientas solve y dsolve nos han permitido, por ejemplo, estudiar las condiciones que deben darse para que una geodésica en el universo de Gödel sea circular, o determinar qué condiciones iniciales imponer sobre la componente temporal de la cuadrivelocidad cuando manualmente escogíamos el resto (sección 3.4.4). En denitiva, MATLAB es una herramienta ecaz y versátil para llevar a cabo estudios sobre las geodésicas de un espacio-tiempo que además resulta intuitiva para usuarios no familiarizados con el cálculo simbólico. Lo explorado a lo largo de este proyecto no constituye si no una fracción de su verdadero potencial para esta materia y para aquellos con interés en el campo de la relatividad general puede convertirse en un valiosísimo aliado a la hora de llevar a cabo sus estudios. 50 Apéndice A: Valores de las constantes y magnitudes utilizadas Magnitud Símbolo Valor Unidades Referencia Velocidad de la luz en el vacío c 299792458 m/s [17] Constante de gravitación universal G6.67430 ·10−11 m 3 /(kg · s 2 ) [18] Semieje mayor de Mercurio a57.909 ·109 m [7] Excentricidad de la órbita de Mercurio e 0.2056 - [7] Masa de Mercurio MMERC 0.33011 ·1024 kg [7] Afelio de Mercurio ra69.817 ·109 m [7] Velocidad de Mercurio en el afelio va38.86 ·103 m/s [7] Período orbital de Mercurio en días TMERC 87.969 d [7] Período de la Tierra en días - 365.256 d [19] Masa del Sol M1988500 ·1024 kg [20] Radio del Sol R6.957 ·108 m [20] Apéndice B: Scripts MATLAB auxiliares para la sección 3.4 Geodésicas circulares en el universo de Gödel 1 % ------------------------------------------------------------------------ 2 % Resolucion simbolica de las ecuaciones geodesicas para la metrica de Godel 3 % imponiendo que la geodesica resultante sea circular. 4 % Script realizado por Mario Misas Arcos para el TFG: 5 % 'Relatividad General con MATLAB' tutorizado por Alberto 6 % Tomas Perez Izquierdo y Carlos Soria del hoyo, Junio de 2021. 7 % Dpto. de Electronica y Electromagnetismo, Facultad de Fisica, 8 % Universidad de Sevilla. 9 %------------------------------------------------------------------------ 10 11 %-------------------------------------------------------------------------- 12 % Inicio 13 %-------------------------------------------------------------------------- 14 clear all % Debemos desactivar clear all si vamos a utilizar este script despues de otro 15 close 16 clc 17 % Definimos las coordenadas a utilizar: Esta vez, anulamos la dependencia 18 % de r y z 19 syms t(s) r phi(s) z 20 % Utilizamos la funcion integrada en basicTensor que asigna +1 y -1 a up y 21 % dw automaticamente 22 [up, dw] = basicTensor.defineUpDw; 23 % Cuadrivector espacio-tiempo 24 Spacetime=[t,r,phi,z]; 25 % Seleccion del semirradio de Godel, esocogido uno para trabajar en 26 % distancias en unidades de este. 27 ag = 1; 28 %-------------------------------------------------------------------------- 29 % Construccion del tensor metrico 30 %-------------------------------------------------------------------------- 31 % Pasamos de symfun a sym el 4vector espacio-tiempo 32 x=Spacetime(s); 33 % Componentes del tensor metrico 34 gtt = -1; 35 grr = (1+(r/(2*ag))^2)^(-1); 36 gpp = r^2*(1-(r/(2*ag))^2); 37 gzz = 1; 38 gtp = -r^2/(ag*sqrt(2)); 39 % Construccion de la matriz. Como el tensor metrico solo depende de r que 40 % ya no depende de s, ya es sym y no symfun. 41 metric = diag([gtt,grr,gpp,gzz]); 42 gdmatrix=metric; 43 gdmatrix(1,3) = gtp; 44 gdmatrix(3,1) = gdmatrix(1,3); 45 % Estructura de tensor para las contracciones: 46 gd = basicTensor(gdmatrix,[dw dw]); 47 % Tensor adjunto: 48 gu = basicTensor(inv(gd.compMatrix), [up up]); 49 %-------------------------------------------------------------------------- 50 % Calculo de los simbolos de Christoffel 51 %-------------------------------------------------------------------------- 52 % Calculo de los simbolos de Christoffel: Creamos el bucle for que calcula 53 % el tensor de orden 3 a contraer con el adjunto del metrico. 54 for p = 1:4 55 for b = 1:4 56 for a = 1:4 57 GM(p,a,b) = diff(gd.compMatrix(a,p),x(b)) + diff(gd.compMatrix(b,p),x(a)) ... 58 - diff(gd.compMatrix(b,a),x(p)); 59 end 60 end 61 end 62 % Asignamos un valor a los caracteres griegos para trabajar de forma mas 63 % visual con tensorContraction 64 alfa = 1; 65 beta = 2; 66 nu = 3; 67 rho = 4; 68 NU = 0; 69 % Contraemos y simplificamos: 70 G=basicTensor((1/2)*GM,[dw dw dw]); 71 F = tensorContraction(G,gu,[nu alfa beta;rho alfa beta;nu rho NU]); 72 Gamma = basicTensor(simplify(F.compMatrix),F.idxMask); 73 %-------------------------------------------------------------------------- 74 % Construccion de las ecuaciones geodesicas 75 %-------------------------------------------------------------------------- 76 % Bucle que construye el termino de las derivadas de las coordenadas: 77 for i = 1:4 78 for j = 1:4 79 dadb(i,j) = diff(x(i),s)*diff(x(j),s); 80 end 81 end 82 % Le damos estructura de tensor para contraerlo 83 DADB = basicTensor(dadb,[up up]); 84 % Contraccion de las derivadas primeras con los S.Christoffel: 85 GDS = tensorContraction(Gamma,DADB,[nu NU NU;nu alfa beta ;alfa beta NU]); 86 gds = GDS.compMatrix; 87 % Conversion a symfun para ser utlizado en odeToVectorField y construccion 88 % de las ecuaciones; 89 gds=symfun(gds,s); 90 eqn = diff(Spacetime,s,s) == -gds; 91 %-------------------------------------------------------------------------- 92 % Resolucion de las ecuaciones resultantes 93 %-------------------------------------------------------------------------- 94 % En la memoria se ilustra la resolucion para la ecuacion en phi. Si 95 % escogemos la ecuacion en t, debemos llegar a lo mismo: para ello 96 % definimos como variables simbolicas las magnitudes usadas en la tetrada 97 % estandar. 98 syms zeta F 99 % Como vamos a escoger beta = 0, sabemos que: 100 % -. dt/dtau(0) = Gamma (Llamada F) 101 % -. dphi/dtau(0) = zeta Gamma 102 % Por lo que rescatamos las dos primeras ecuaciones que da eqn, que son 103 % las que nos interesan. Para ello pasamos de symfun a sym: 104 EQ = eqn(s); 105 AQ = [EQ(1) EQ(2)]; 106 % Ahora vamos a aplicar dsolve: 107 SOL = dsolve(AQ); 108 % Que nos da dos soluciones para phi, una constante (que claramente no es 109 % la solucion correcta pues buscamos una geodesica circular) y una para t. 110 % Asi seleccionamos la correcta en phi y cualquiera en t (dsolve da dos 111 % soluciones que son la misma): 112 phis = SOL.phi(2); 113 ts = SOL.t(1); 114 % Y como hemos impuesto que dt/dtau(0) = Gamma (Llamada F), tendremos que 115 % C1 = Gamma, y escogiendo t(0) = 0; C2 = 0. Para phi sabemos que C4 sera 116 % cero esocogiendo phi(0) = 0. Para encontrar zeta debemos declarar C1 117 % como variable simbolica: 118 syms C1 119 % Para poder sustituirla por Gamma y despues imponer que la derivada de phi 120 % en cero sea igual a Gamma z: 121 phis = subs(phis,C1,F); 122 % Por ultimo, usamos solve para despejar zeta: 123 % Derivada de phi en s = 0: 124 Dp0 = subs(diff(phis,s),s,0); 125 % Y llegamos a que zeta es: 126 zeta = solve(Dp0 == F*zeta) 127 % El resultado obtenido deberia ser el mismo que el visto en la memoria Homogeneidad del universo de Gödel a través de las geodésicas 1 % ------------------------------------------------------------------------ 2 % Script auxiliar para hallar las condiciones inciales sobre las derivadas 3 % de r y phi a la hora de construir geodesicas circulares para otros 4 % observadores. 5 % Script realizado por Mario Misas Arcos para el TFG: 6 % 'Relatividad General con MATLAB' tutorizado por Alberto 7 % Tomas Perez Izquierdo y Carlos Soria del hoyo, Junio de 2021. 8 % Dpto. de Electronica y Electromagnetismo, Facultad de Fisica, 9 % Universidad de Sevilla. 10 %------------------------------------------------------------------------ 11 12 %------------------------------------------------------------------------ 13 % Planteamiento 14 %------------------------------------------------------------------------ 15 % Variables simbolicas: parametro s, Gamma prima, zeta prima y los 16 % parametros espaciales (posicion, radio) de las geodesicas