Estudio de la reacción de captura radiativa 14C(n,g)15C

Estudio de la reacción de captura radiativa 14 C (n, γ)15 C Autor: Domingo Parrales de la Cruz Tutores: Jesús Casal Berbel Antonio Matias Moro Muñoz Julio 2023 Índice 1. Introducción 4 1.1. Núcleoshalo................................... 5 1.2. Interés astrofísico del 15 C............................ 9 1.3. El núcleo de 15 C................................. 10 1.4. Objetivos .................................... 10 1.5. Metodología................................... 11 1.6. Organización de la memoria . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2. Fundamento teórico 12 2.1. Solución de la ecuación de Schrödinger . . . . . . . . . . . . . . . . . . . . 12 2.2. Probabilidad de transición eléctrica . . . . . . . . . . . . . . . . . . . . . . 13 2.2.1. Momento multipolar eléctrico . . . . . . . . . . . . . . . . . . . . . 14 2.3. TransicionesE1................................. 16 2.3.1. Solución radial. Normalización y aproximación asintótica. . . . . . . 16 2.3.2. Cálculo analítico de la probabilidad de excitación reducida . . . . . 19 2.4. Sección ecaz de captura y tasa de reacción . . . . . . . . . . . . . . . . . 19 2.5. Modelodeestructura.............................. 20 3. Resultados 22 3.1. Constante de normalización asintótica . . . . . . . . . . . . . . . . . . . . . 22 3.2. Comparación con el cálculo numérico . . . . . . . . . . . . . . . . . . . . . 22 3.3. Comprobación del modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3.1. Comprobación de las aproximaciones . . . . . . . . . . . . . . . . . 27 3.4. Sección ecaz: ruptura y captura radiativa . . . . . . . . . . . . . . . . . . 28 3.4.1. Sección ecaz de captura: primer estado excitado . . . . . . . . . . 29 3.5. Tasadereacción ................................ 33 4. Resumen y conclusiones 34 A. Funciones de onda: integral analítica 36 B. THOx 37 B.1.Input....................................... 38 1 B.2.Output...................................... 41 C. Matlab 41 C.1.Lecturadedatos ................................ 42 C.2.Coecientesdeajuste.............................. 43 C.3.Cálculos..................................... 43 Referencias 45 2 Resumen En este trabajo estudiamos la reacción de captura 14 C (n, γ)15 C, que es relevante en varios sistemas astrofísicos. Para ello, nos basamos en las propiedades de núcleo halo que presenta el 15 C, que nos permite simplicar el problema a un modelo de dos cuerpos ( n+14 C). Después, resolvemos la ecuación de Schrödinger numéricamente, con lo que se consideran las funciones de onda exactas para los estados inicial y nal, y mediante un desarrollo analítico, para el que consideramos la aproximación asintótica de los estados que sean ligados y la de ondas planas para los del continuo. A partir de estos resultados, se calculan la sección ecaz de captura radiativa y su tasa de reacción correspondiente Abstract In this work we focus on the 14 C (n, γ)15 C capture reaction, with astrophysical relevance. To do that, we take into account the halo properties of 15 C nuclei, which allow as to simplify the problem to a two-body system ( n+14 C). Then we solve Schrödinger's equation numerically, with wich we obtain the exact wave functions, and with an analitic model, considering only the asymptotic behaviour of the bound states and plain wave solutions for the continuum states. With those, we calculate the radiative cross section and its reaction rate. 3 1. Introducción Con el desarrollo de los aceleradores de haces radiactivos se ha vuelto posible el estudio de núcleos inestables, aquellos con semivida corta, que, además de presentar algunas propiedades diferentes a las de los núcleos estables, resultan ser de interés en áreas como la astrofísica, debido a su presencia en los procesos de formación de núcleos pesados y en sistemas, como estrellas o incluso supernovas. En la Figura 1, se muestra un diagrama en el que cada punto corresponde a un núcleo y se representa el número de protones que tiene frente a su número de neutrones. Los núcleos inestables serán aquellos que se alejen de la línea central de núcleos estables, representada con cuadrados negros. Vemos que podemos alejarnos de esta aumentando el número de protones, Z , o el número de neutrones, N , pero conforme lo hagamos estaremos reduciendo la energía mínima necesaria para separar un nucleón N del núcleo, SN , hasta que, llegado cierto punto, se tenga que SN= 0 , provocando que no podamos añadir ningún otro nucleón del tipo N al núcleo. Es por esto que en la Figura 1 vemos que solo podemos medir núcleos en el entorno de la línea de núcleos estables hasta cier- to límite, siendo estos bordes las líneas de evaporación. Por tanto, podemos distinguir entre la línea de evaporación de neutrones (límite inferior) y la línea de evaporación de protones (línea superior), y en general será cerca de estas donde se tengan sistemas nucleares que presenten propiedades inusuales, a los que nos referimos como núcleos exóticos. Una consecuencia de la disminución de SN será un aumento del radio nuclear, R , que para núcleos estables obedece aproximadamente la regla R=rA1/3 (1) donde r= 1,2 fm se ha ajustado empíricamente y A=Z+N es la suma del número de protones y neutrones respectivamente. Sin embargo, existen evidencias de núcleos ligeros cerca de la línea de evaporación, típicamente de neutrones, y con una energía de separación del último neutrón baja, Sn , en los que se mide un radio superior al que se espera con la ecuación (1), lo que provoca un incremento notable en la sección ecaz de ruptura, que se ha interpretado como la formación de una nube o halo en torno al resto de nucleones, a los que nos referiremos como core [11]. 4 Figura 1: Carta de Segrè de NuDat [12]. Respecto de los núcleos estables, puntos negros, distinguimos principalmente dos zonas de núcleos inestables: los que están por encima tienen un exceso de protones y los que están por debajo un exceso de neutrones. Los límites de estas zonas son las líneas de evaporación. 1.1. Núcleos halo La categoría de núcleo exótico engloba núcleos con diversas propiedades, como puede verse en la Figura 2. Por ejemplo, tenemos los núcleos halo, representados con fondo azul cuando este halo los forman protones y en rojo cuando lo forman neutrones. Este halo puede estar compuesto por más de una partícula, y principalmente suelen ser neutrones, debido al connamiento que genera la barrera coulombiana por la interacción con otros protones. No obstante, otro efecto que limita la formación del halo, incluso para el de neutrones, se debe al potencial centrífugo en la supercie nuclear, que viene dado por L(L+ 1)/2A2/3 MeV [16], siendo L el momento angular orbital del nivel en el que se encuentran los nucleones del halo. Por tanto, se verá favorecida la formación de halos de neutrones en niveles con l= 0 (ondas s). Los núcleos halo más estudiados son 6 He, 11 Li y 11 Be, pero también se tienen otros como el 14 Be, 14 B, 15 C y 19 C. En este trabajo nos centramos en el 15 C, que tiene un halo de un neutrón, al que solemos referirnos como neutrón de valencia, y, por tanto, un core 5 Figura 2: Representación de uno de los sectores de la Figura 1 en la que tenemos un número de nucleones pequeño [15]. Cada núcleo se representa por un cuadrado con un número en su interior, que es el número de nucleones, A . La línea de núcleos estables se representa en negro y los núcleos halo se representan con fondo azul (halo de protones) y en rojo (halo de neutrones), con un círculo si se tiene un sistema nucleón+ core y con dos círculos si se tiene un sistema nucleón+nucleón+ core , que solemos llamar núcleos borromeo . de 14 C que podemos suponer inerte, porque la energía mínima necesaria para que el 14 C se excite es de 6093 keV [26], mientras que para el 15 C se tiene que Sn= 1218 keV. Por tanto, podemos considerar que el 15 C es un sistema de dos cuerpos (neutrón+ core ) en el que no hay interacción coulombiana (porque el neutrón tiene carga eléctrica nula) y que el 14 C es inerte. Podemos suponer entonces que el neutrón de valencia se encuentra sometido a un potencial producido por el core . Por ejemplo, podemos suponer un potencial de Woods- Saxon, como el de la Figura 3, en donde tendremos un número nito de soluciones para ciertas energías que no superen el valor del potencial en las paredes (estados ligados), e innitas soluciones para cada valor de la energía por encima de la de las paredes del pozo (estados del continuo). Además para los estados ligados podremos distinguir: Zona clásicamente permitida (solución oscilante) : cuando la energía correspondiente a la función de onda es superior al potencial. 6 Figura 3: Función de onda (línea azul) en un potencial de Woods-Saxon, en el que se distingue la zona clásicamente permitida (solución oscilante) y la clásicamente prohibida (decaimiento exponencial). Zona clásicamente prohibida (decaimiento exponencial) : cuando la energía de la función de onda es inferior al potencial. En esta zona la función de onda presenta un decaimiento exponencial que será más rápido cuanto mayor sea la diferencia de energía respecto del límite superior de la barrera, en este caso, cero. Entonces, en el caso de un núcleo halo, como la energía de separación es menor que la de los núcleos ordinarios, se espera que haya un lento decaimiento en la zona clásicamente prohibida, lo que implica que es más probable que el nucleón menos ligado penetre en la región clásicamente prohibida [3]. Esto provocará que el radio que se calcule para un núcleo halo sea mayor que el esperado, aunque no podemos medirlo directamente. No obstante, es posible inferir ese radio a partir de otros observables, por ejemplo: Sección ecaz de interacción con otros núcleos: hacemos incidir un haz de los núcleos que queremos estudiar sobre un blanco, que será un núcleo conocido. Este actuará como una barrera y mediremos la intensidad, I , del haz que la atraviesa en relación a la intensidad incidente, I0 , según I=I0exp(−σIt) (2) siendo t el espesor de la barrera y σI la sección ecaz de interacción, que guarda relación con los radios de los núcleos usados como proyectil y como blanco. Por tanto, 7 11 Li 6 He 11 Be 15 C 8,75 ms 806,7 ms 13,76 s 2,449 s Cuadro 1: Comparación de la semivida de núcleos halo bien estudiados, como son el 11 Li, el 6 He y el 11 Be, con la del 15 C. podemos estimar el del proyectil ( 15 C) si conocemos el del blanco, que podemos tomarlo por ejemplo de 12 C. Sección ecaz de la dispersión Rutherford: podemos ver que se comporta como un sistema de dos cuerpos midiendo su ruptura ya que, debido a su baja Sn respecto de otros isótopos de C, debería romperse el sistema de dos cuerpos y dejar de presentar dispersión elástica antes que los otros isótopos [13]. Dada su estructura, en el caso de blancos pesados ( Z grande) deberíamos ver favorecida una ruptura análoga a la del dipolo eléctrico, es decir, se favorecen las transiciones del neutrón de valencia de un estado ligado a uno del continuo (ruptura) del tipo E1 (dipolo eléctrico). Momento del core : mediante reacciones knockout 1 hacemos perder un neutrón al 15 C y medimos el momento del core de 14 C, ~p , resultante, que será cercano al momento antes de la interacción. Para núcleos halo esperamos que la dispersión en el valor de ~p sea estrecha en comparación con la de otros isótopos de C, por lo que según el principio de indeterminación, deberá presentar una dispersión en la coordenada radial r notablemente superior. Otro factor a tener en cuenta es que el 15 C tiene una semivida grande, T1/2= 2,449 s, comparable a la de los núcleos más estudiados (Cuadro 1), con lo que vemos que se encuentra entre el 6 He y el 11 Be, y será posible usar el 15 C en reacciones nucleares para estudiar sus propiedades. 1 Las reacciones knockout consisten en hacer incidir sobre un blanco un haz del núcleo que queramos estudiar a alta velocidad, de forma que se pierda uno o dos nucleones del proyectil sin que haya tiempo para otro tipo de reacciones. 8 pero teniendo en cuenta las condiciones anteriores, podemos simplicar notablemente esta expresión: h(lf,1/2)kJf||Yλµ(ˆr)||(li,1/2)Jii= (−1)λ+lf+Ji+1/2(−1)Ji+Jf+λ ×hλ0, li0|lf0iˆ Jiˆ liˆ Jiˆ Jfr(2λ+ 1) 4π   li1/2Ji Jfλ lf     Ji0Ji Jfλ Jf   . La estructura del último coeciente 6- j es análoga a la de [4]    abc 0c b    =(−1)a+b+c ˆ bˆc, (18) debido a las relaciones de permutación de los coecientes 6- j    Ji0Ji Jfλ Jf   =   0JiJi λ JfJf   =   λ JfJi 0JiJf   =(−1)λ+Jf+Ji ˆ Jiˆ Jf . (19) Entonces h(lfs)kJf||Yλµ(ˆr)||(lis)Jii= (−1)λ+lf+Ji+1/2ˆ Jiˆ lf   li1/2Ji Jfλ lf  hlf||Yλ||lii, (20) siendo hlf||Yλ||lii=hλ0, li0|lf0iˆ li ˆ lfr(2λ+ 1) 4π. (21) Finalmente, podemos dar la distribución de probabilidad combinando las ecuaciones (16), (20) y (11) dB d (Eλ, Ji→kJf) = 2Jf+ 1 2Ji+ 1 µk (2π)3~2× X lf  eZ(λ) eff 4π k(−i)lfI(lf,1/2)Jf (li,1/2)Ji(λ, r)(−1)λ+lf+Ji+sˆ Jiˆ lf   li1/2Ji Jfλ lf  hlf||Yλ||lii 2 , y, usando la ecuación (21), dB d (Eλ, Ji→kJf) = 2Jf+ 1 2Ji+ 1 µk (2π)3~2eZ(λ) eff 2× X lf  4π kr(2λ+ 1)(2Ji+ 1)(2li+ 1) 4π   li1/2Ji Jfλ lf  hλ0, li0|lf0i 2 I(lf,1/2)Jf (li,1/2)Ji(λ, r) 2. (22) 15 2.3. Transiciones E1 Ahora vamos a centrarnos en el caso de excitaciones dipolares ( λ= 1 ) para sistemas tipo halo con un neutrón de valencia ( Zv= 0 ). Por tanto, según la ecuación (14), tenemos que Z(λ=1) eff (Zv= 0) ≡Zeff =Zc−mv mc+mv, (23) luego dB d (E1, Ji→kJf) = 2Jf+ 1 2Ji+ 1 µk (2π)3~2(eZeff )2 4π kr3(2Ji+ 1)(2li+ 1) 4π!2 × X lf    JfJi1 lilf1/2  h10, li0|lf0iI(lf,1/2)Jf (li,1/2)Ji(1, r) 2 . (24) Por tanto, la mayor parte del cálculo consistirá en resolver la integral I(lf,1/2)Jf (li,1/2)Ji(1, r) , para lo que necesitamos la parte radial, u(r) , de las soluciones de la ecuación de Schrödinger. Esta dependerá de si la función de onda representa estados del continuo o estados ligados. Si tomamos como cero de energía el valor del potencial en r→ ∞ tendremos que los estados con energía negativa son ligados y los de energía positiva son del continuo. 2.3.1. Solución radial. Normalización y aproximación asintótica. Ahora vamos a plantear simplicaciones sobre las soluciones radiales que nos permitan resolver la integral I(lf,1/2)Jf (li,1/2)Ji(1, r) de forma sencilla. Para ello nos centraremos en la región asintótica ( r >> ), ya que, como vemos en la ecuación (17), su valor depende fundamentalmente de la región de separaciones grandes. Estados ligados: en este caso tenemos Er<0 , por lo que las funciones radiales están normalizadas Z∞ 0 dr |f(lisi)Ji(r)|2= 1, (25) lo que signica que son funciones de cuadrado integrable. El comportamiento asintótico de estas funciones viene dado por una función de Whittaker (ecuación (19) de [2]) f(lisi)Ji(r)→ClsjW−η,l+1/2(2kBr), (26) donde tenemos el coeciente de normalización Clsj , el parámetro de Sommerfeld η y el número de onda del estado ligado k2 B= 2µSn/~2 , siendo Sn=|Er| . Sin embargo, en 16 nuestro caso estamos tratando un sistema netrón + core , por lo que el parámetro de Sommerfeld (que depende del producto de las cargas) es nulo. Esto nos permite reducir la complejidad del problema, ya que la función de Whittaker se reduce a una función de Hankel de argumento complejo f(li,1/2)Ji(r)→Dr h(1) li(ikBr), (27) siendo h(1) li(r) la función de Hankel esférica de orden 1, li el momento angular del estado y D una constante de normalización. Estados del continuo: en este caso Er>0 , por lo que no son normalizables, es decir la integral diverge Z∞ 0 dr |g(lfsf)Jf(k, r)|2=∞. (28) Sin embargo, una normalización es necesaria para poder usarla, y es por esto que la expresamos mediante deltas de Dirac [10] Z∞ 0 dr g∗ (lfsf)Jf(k, r)g(lfsf)Jf(k0, r) = π 2δ(k−k0), (29) donde el factor π/2 se añade debido a la denición de g(lfsf)Jf(k, r) , ya que cuando estudiamos su comportamiento asintótico se obtienen las funciones de Coulomb [14] g(lfsf)Jf(k, r)→exp i(σlf+δlf)×(cos(δlf)F(kr) + sin(δlf)G(kr)), (30) donde la fase no juega ningún papel relevante en estos cálculos, ya que estamos trabajando con el módulo al cuadrado de la integral. Por otro lado, k2= 2µ/~2 , F y G son las funciones regular e irregular de Coulomb respectivamente, y δlf es el corrimiento de fase, que almacena la información de la interacción en el estado nal. Por tanto, si despreciamos el efecto del potencial en los estados nales, tendremos que δlf= 0 . Así, podremos reescribir g(lfsf)Jf(k, r)→F(kr) , que podemos desarrollar en términos de las funciones esféricas de Bessel [14] Flf(η, ρ) = (2lf+ 1)!! ρ Clf(η) ∞ X k=lf bkjk(ρ), (31) siendo jk la función de Bessel esférica, y los coecientes bk y Clf dependen de η según blf= 1, bk>lf∝η= 0, Clf(η) = 2lf e −πη/2|Γ(lf+1+iη)| Γ(2lf+ 2) =2lfΓ(lf+ 1) Γ(2lf+ 2) =1 (2lf+ 1)!!, 17 con lo que nos acaban quedado ondas planas g(lfsf)Jf(k, r)→Flf(η= 0, kr) = krjlf(kr) = rπ kr 2Jlf+1/2(kr), (32) donde Jl+1/2 es la función de Bessel, y el factor pπ/2 que lleva multiplicando es el responsable de que aparezca el factor π/2 al hacer el producto g∗ (lfsf)Jf(k, r)g(lfsf)Jf(k0, r) . Por tanto, si el neutrón se encuentra en la región en la que la energía de ligadura es cercana a cero (o equivalentemente, la energía de separación Sn es cercana a cero) podremos tomar el límite asintótico de la función de Hankel [6] y reescribir el estado ligado f(li,1/2)Ji(r)→Dr lim r→∞ h1 li(ikBr) = Dr (−i)li+1 e−kBr ikBr. (33) Para el estado del continuo vamos a jar lf para facilitar los cálculos. Esto lo haremos teniendo en cuenta que las transiciones más probables serán aquellas que partan del estado fundamental 1/2+ , que es un estado con li= 0 (ondas s), y, en consecuencia, al pasar al continuo mediante transiciones E1 debemos acabar en estados con lf= 1 (ondas p), que serán en nuestro caso el nivel 1/2− y el 3/2− . Por tanto, nos queda g(lfsf)Jf(k, r)→g(r)≡krjlf=1(kr) = sin(kr)−kr cos(kr) kr . (34) Por otro lado, la función del estado ligado quedará f(li,1/2)Ji(r)→f(r)≡Dr e−kBr kBr, (35) donde hemos incluido el signo negativo en la constante D . Con esto ya podemos hacer el cálculo analítico de la parte radial (desarrollo en el Anexo A) I ≡ I(1,1/2)Jf (0,1/2)1/2(1, r) = DZ∞ 0 drsin(kr)−kr cos(kr) kr re−kBr=2Dk2 (k2+k2 B), (36) válido tanto para Jf= 1/2 como para Jf= 3/2 . Y, si además tomamos las deniciones de kB y k de las ecuaciones (27) y (30), podemos dar la parte radial en función de las energías I= 2D~2 2µ  (2+S2 n)2. (37) 18 2.3.2. Cálculo analítico de la probabilidad de excitación reducida Sustituyendo en la ecuación (24) los momentos angulares elegidos para el cálculo anterior dB d (E1, Ji= 1/2→kJf) = (2Jf+ 1) 12πµc2 (2π)3(~c)2 1 k(eZeff )2   Jf1/2 1 0 1 1/2   2 |I|2, donde se ha tomado h10,00|10i= 1 , y Jf depende del estado del continuo al que se excite. Como por transiciones E1 desde 1/2+ tenemos dos opciones ( 1/2− y 3/2− ), siguiendo la ecuación (12), sumamos la aportación de cada uno para calcular la probabilidad total dB d (E1) = dB d (E1, Ji= 1/2→kJf= 1/2) + dB d (E1, Ji= 1/2→kJf= 3/2). (38) donde lf= 1 en ambos casos, debido a que tenemos transiciones E1 y partimos de li= 0 . Sin embargo, la forma del coeciente 6-j es idéntico al de la ecuación (18) y, al tomar su módulo, este término se vuelve independiente de Jf , de forma que la probabilidad de excitación reducida solo depende de Jf a través del factor (2Jf+ 1) . Por tanto: dB d (E1) = 21 2+ 1+23 2+ 1 12πµc2 (2π)3(~c)2 1 k(eZeff )2 1 √6 2 |I|2, (39) y sustituyendo la ecuación (37) llegamos a que dB d (E1) = 3 √2C(eZeff )2 2π2(~c)2 µc23/23/2 (+Sn)4, (40) donde se ha denido C≡ |D|2 . Como podemos ver, esta expresión coincide con la ecuación (5) del trabajo de Nagarajan [5] salvo por la constante de normalización, ya que, en su caso, lo ha calculado con la ecuación (25), obteniendo C= 2kB . En nuestro caso, vamos a deducirla a partir de una ajuste de f(r) en la región asintótica, como se hace en el artículo de Typel y Baur [2]. Por otro lado, se observa que si Sn disminuye, el máximo de la distribución se dará a una energía más baja, y si nos jamos en los valores límite de Sn (cero e innito) también se aprecia que el valor de B(E1) debe crecer conforme Sn disminuya. 2.4. Sección ecaz de captura y tasa de reacción El cálculo de σrc lo hacemos partiendo de la reacción inversa, la reacción de fotodisociación 15 C (γ, n)14 C. La expresión de σγ la extraemos de [8] en la que se sigue un 19 procedimiento análogo para hallar σrc de 7 B (p, γ)8 B , σ(λ) γ(γ) = (2π)3(λ+ 1) λ[(2λ+ 1)!!]2γ ~c2λ−1dB d (πλ), (41) donde γ es la energía del fotón γ=+Sn necesario para provocar la transición desde el estado fundamental al estado del continuo con energía  . Por el principio de balance detallado, podemos dar la de la reacción inversa, en este caso la de captura a+b→c+γ σ(λ) rc () = 2(2Ic+ 1) (2Ia+ 1)(2Ib+ 1) kγ k2 σ(λ) γ(γ), (42) siendo Ia , Ib e Ic los espines de las partículas involucradas y kγ=γ/(~c) el número de onda del fotón. En nuestro caso solo consideramos que λ= 1 , ya que, si nos jamos en la ecuación (41) vemos que cuanto mayor sea λ menor serán los términos que acompañan a dB/d(E1) . También, a partir de este resultado podemos calcular la tasa de reacción, R(T) , que representa el ritmo con el que ocurre la reacción en un entorno astrofísico. Lo calculamos promediando σrc con una distribución de Maxwell-Boltzmann, y suele escribirse como [1] R(T) = NA8 µπ(kBT)31/2Zσrc()exp − kBTd, (43) siendo NA el número de Avogadro, T la temperatura y kB la constante de Boltzmann. 2.5. Modelo de estructura Para comprobar hasta qué punto la aproximación asintótica y la fórmula analítica son válidas, vamos a calcular también dB/d(E1) resolviendo la ecuación de Schrödinger directamente con el programa THOx [27]. Para ello, vamos a considerar que tenemos un sistema neutrón de valencia + core de 14 C. Este último producirá un potencial efectivo que, siguiendo la elección habitual, tendrá un termino central, con la forma de Woods-Saxon, y un término no central de espín-órbita, como el empleado para describir los estados de 15 C en el artículo de Capel et al. [9]: Vcf (~r) = V0(r)+4~ L·~ IVso(r), (44) donde Vcf representa las interacciones internas entre el core y el neutrón de valencia, ~ L es el momento angular orbital relativo e I es el espín del neutrón. El primer sumando representa el término central V0(r) = −Vlf(r, R0, a), (45) 20 con el factor de forma de Woods-Saxon f(r, R0, a) = 1 1 + exp r−R0 a, (46) y además, el término del acoplamiento espín-órbita vendrá dado por Vso(r) = VLS 1 r d drf(r, R0, a), (47) donde R0 es el radio reducido y viene dado por la ecuación (1); y a , que lo hemos tomado como 0,6 fm, es un parámetro relacionado con la difusividad del campo. Otros parámetros que tenemos que introducir en el programa son Vl y VLS . Estos se eligen para que el sistema reproduzca los estados ligados y una de las resonancias del 15 C [9]: Vl=0 = 52,814 MeV y Vl>0= 51,3 MeV. También tenemos que especicar características como el mallado en r y en energías para las funciones del continuo, las características del neutrón y del core , y los estados iniciales y nales en el programa, pero estos comentan en detalle en el Anexo B. De esta forma es posible calcular el elemento de matriz I numéricamente, y, en consecuencia, dB/d (E1), σγ , σrc y R(T) , a partir de soluciones más generales para la parte radial de las funciones de onda. 21 3. Resultados 3.1. Constante de normalización asintótica Como se comentó al obtener la expresión analítica de dB/d(E1) a partir de la aproximación asintótica, ecuación (40), vamos a deducir la constante de normalización D ajustando la función que marca nuestro modelo, es decir, una de la forma f(r;D) = D e −kBr , a la función de onda del estado ligado generada numéricamente por los tutores con el programa FRESCO [28] 2 . Si nos jamos en la Figura 5 podemos ver la función de onda numérica, en línea continua, junto con varios ajustes de la función analítica para distintos valores de D , que se obtienen utilizando la función t de Matlab (Anexo C). Cada ajuste comienza en un valor de r diferente, por ejemplo, el primero se hace desde rmin = 0 , es decir, se toman todos los puntos, por lo que este resultados será equivalente a la normalización usual, descrita en la ecuación (25). Es evidente que este ajuste no es correcto para la región asintótica, pero conforme se toma un rmin mayor el ajuste mejora, hasta que al pasar del mínimo se tiene un comportamiento aproximadamente exponencial. De esta forma, evaluando el coeciente tomando varios valores de rmin &5 fm se obtiene un valor para el coeciente D=−1,37 fm 1/2 , que podemos vericar que se comporta correctamente representando |f(r;D=−1,37 fm 1/2)| en escala logarítmica junto al valor absoluto de las funciones ligadas generadas numéricamente, como puede verse a la derecha de la Figura 5. 3.2. Comparación con el cálculo numérico Al realizar el cálculo numérico de la probabilidad de excitación lo que en primer lugar se espera observar es que la transición al estado 3/2− sea mucho más probable que la del estado 1/2− por la degeneración de los estados, g(1/2−) = 2 y g(3/2−) = 4 . Esto se debe a que, como se vio en la ecuación (39), la diferencia entre la dB/d(E1) de los dos estados viene dada por un factor (2Jf+ 1) ≡ deg (Jf) . En la gura 6 se observa este comportamiento y además se ve que efectivamente el área bajo la curva correspondiente 2 El motivo de usar FRESCO en lugar de THOx es que este último genera las funciones de onda radiales mediante la normalización del hamiltoniano en una base de estados normalizables. Esto, en general, describirá bien la función de onda en el interior, pero no asintóticamente. Sin embargo, con FRESCO se obtienen por integración directa de la ecuación de Schrödinger. 22 Figura 5: En ambas grácas tenemos que la línea continua es la función de onda del estado ligado generada numéricamente con FRESCO y el resto son los ajustes f(r;D) = Dexp(−kBr) , con D el parámetro de ajuste. En la de la izquierda tomamos el valor mínimo de r para el ajuste ( rmin ) cada vez más lejos de r= 0 , es decir, cada vez se toma rmin más grande. En la gráca derecha tenemos en escala logarítmica el ajuste exponencial para el valor de D calculado. Figura 6: Cálculo numérico de la probabilidad de excitación para transiciones E1 para la energía de separación Sn(15C) = 1,22 MeV. Las líneas discontinuas son las transiciones desde el estado fundamental 1/2+ a los del continuo 1/2− y 3/2− . La línea continua es la probabilidad total de transición, que se calcula sumando las dos contribuciones anteriores. 23 Figura 7: La línea continua representa dB/d(E1) calculada numéricamente y el resto de líneas se reeren al cálculo analítico de dB/d(E1) a partir de la ecuación (40). La curva negra en línea de puntos y rayas se ha hecho para Sn(15C) = 1,22 MeV, y las otras se han tomado para algunos valores en torno a él. al 3/2− es aproximadamente el doble. Sumando la aportación de cada transición podemos dar la dB/d(E1) total, que es con la que compararemos el resultado analítico. En la Figura 7 se representa en línea discontinua los valores de dB/d(E1) obtenidos con la expresión analítica (40) para distintos valores de Sn . También hemos representado en línea continua los resultados del cálculo numérico de dB/d(E1) , y, del mismo color, la curva analítica para la misma Sn= 1,22 MeV. Al compararlas vemos que hay un acuerdo razonable en la forma, aunque esta última presenta una magnitud un 29 % mayor. 3.3. Comprobación del modelo Hemos visto que, para Sn≈1,22 MeV, hay una diferencia apreciable entre el cálculo analítico y el numérico. No obstante, si reducimos el valor de Sn (o, equivalentemente, si aumentamos la distancia media entre el neutrón y el core ), esperamos que el acuerdo de la B(E1) calculada analíticamente y a partir de las funciones de onda generadas por métodos numéricos mejore, debido a que los cálculo numéricos se harán bajo unas condiciones 24 Figura 11: A la izquierda se representa en línea continua la sección ecaz de captura total, que es la suma de la del estado ligado (línea discontinua) y la del primer estado excitado (línea de puntos y rayas). A la derecha tenemos una zona ampliada de la de la izquierda, centrándonos en energías bajas, en la que se representan datos experimentales de Nakamura [24] (ruptura coulombiana) y de Reifarth [25] (medida directa). A la derecha de la Figura 10 tenemos la sección ecaz de captura, en la que vemos que la principal aportación viene de la transición del 7/2− . Esto es esperable, ya que la B(E1) de este estado es notablemente superior a la de los otros dos, pero además se tienen otros factores que lo hacen incluso más predominante: Por un lado, en la ecuación (42) hay un factor (2Jf+ 1) , por lo que el estado con mayor degeneración, el 7/2− , se ve más favorecido que los otros dos. Por otro lado, en la ecuación (41) hay un factor (γ/~c)2 , que será mayor conforme aumente  . Como la dB/d(E1) de la transición al 7/2− abarca una zona con energías más altas, este término también favorecerá más esta transición que la del estado 5/2− . Por último, podemos ver en la Figura 11 (izquierda) una comparación entre σrc obtenida para el estado fundamental, que es la más relevante a energías bajas, y para el estado excitado, que empieza a ser más importante que la del estado fundamental a partir de los 31 Figura 12: Tasa de reacción como función de la temperatura representada en escala logarítmica en ambos ejes. La línea discontinua es la tasa de reacción para transiciones al estado fundamental, la línea de puntos y rayas es de transiciones al primer estado excitado y en línea continua se tiene la suma de ambas. La zona sombreada marca el rango de temperatura de los sistemas astrofísicos relevantes. 4,5 MeV aproximadamente, y alcanza unos valores más elevados que la del estado fundamental. Además, en la gráca derecha de la Figura 11, se ha representado el resultado obtenido junto con datos experimentales de σrc [23] del trabajo de Nakamura [24] y el de Reifarth [25]. Los datos de Nakamura se obtienen mediante la reacción de disociación coulombiana de 15 C en un blanco de Pb, a partir del que se extrae información de la contribución de las transiciones E1 desde el estado fundamental a la reacción de captura 14 C (n, γ)15 C. Por otro lado, los datos de Reifarth se han tomado por medida directa [25], por lo que medirá fotones de la reacción 14 C (n, γ)15 C tanto del estado ligado como del primer excitado. Lo que se observa de ambos conjuntos de datos es que, efectivamente, sus máximos se encuentran en torno a 1 MeV y siguen una tendencia similar a la que marcan las curvas del cálculo numérico. Además, para energías bajas los datos de Reifarth coinciden con los nuestros, mientras los de Nakamura en general se encuentran bajo la curva de σrc del estado fundamental. 32 3.5. Tasa de reacción Por último, podemos usar la sección ecaz de captura para el cálculo de la tasa de reacción, ecuación (43), que representamos en la Figura 12. Lo más importante de este resultado es el hecho de que la captura al estado fundamental es la reacción principal hasta, aproximadamente, los 10 GK. Para resaltar este hecho, se ha sombreado la zona de la gráca correspondiente al rango de temperaturas relevantes para los sistemas astrofísicos expuestos. En concreto, la temperatura en el núcleo de estrellas AGB está entorno a los T∼3·108 K para las estrellas AGB poco masivas y T∼3,5·108 K para las más masivas [29], y también, en el modelo BBN el rango de temperatura de interés es T∼0,2−1,2 GK [1]. Equivalentemente, el rango de energías de interés estará por debajo de max = 1 GK ×kB≈0,09 MeV, siendo kB la constante de Boltzmann. Otro factor llamativo es que aun siendo σrc del estado excitado mayor que la del estado fundamental, la tasa de reacción favorece a este último. Esto podemos explicarlo basándonos en la ecuación (43), ya que vemos que σrc está ponderada por una distribución de Maxwell-Boltzmann. Concretamente, la exponencial variará en función de la energía característica kBT , de forma que si la temperatura es baja , la exponencial tenderá a cero rápidamente, y solo afectará a la integral los valores de σrc para energías pequeñas, que es justamente donde predomina la del estado fundamental. Sin embargo, cuando T es alta la exponencial tiene argumento próximo a cero para un amplio rango de energías, por lo que el valor de σrc para energías más altas también afectará al cálculo, y la captura al estado excitado comenzará a ser relevante 3 . 3 Al hablar de temperatura baja o alta lo hacemos respecto de la escala de la Figura 12. Por tanto, una temperatura pequeña será inferior a 1 GK, mientras que una temperatura alta será superior a los 10 GK. 33 4. Resumen y conclusiones En este trabajo, hemos estudiado el proceso de captura radiativa 14 C (n, γ)15 C. Este proceso toma importancia en sistemas con una alta densidad de neutrones, como las estrellas AGB o el sistema propuesto en la teoría BBN. En este contexto, la reacción de captura 14 C (n, γ)15 C es relevante para el estudio de la producción de energía a través de ciclos CNO y para explicar el origen de núcleos con A > 20 [1]. Otro punto de interés es que el 15 C es un núcleo con halo de un neutrón, de forma que podemos describir dicho sistema de forma aproximada mediante un modelo de dos cuerpos (neutrón + core ) en ausencia de interacción coulombiana. Después, para calcular las soluciones analíticas, hemos considerado la aproximación asintótica para el estado fundamental y de ondas planas para los estados del continuo, lo que nos ha permitido calcular analíticamente la probabilidad de excitación reducida dB/d( E1 ) partiendo del estado fundamental ( 1/2+ ) a estados de energía positiva p3/2 y p1/2 . Para comprobar la validez de la aproximación asintótica, hemos comparado los resultados con el cálculo numérico de estas transiciones obtenido usando el programa THOx . Para ello se ha tomado un potencial de Woods-Saxon de un estudio previo [9], con el que se han generado las funciones de onda ligadas y del continuo. Además, se ha estudiado el efecto de la aproximación asintótica en los estados ligados y los del continuo por separado, concluyendo que la mayor diferencia proviene de la descripción del estado ligado y que la aproximación asintótica funciona mejor cuanto menor sea la energía de separación Sn . Con el cálculo numérico de dB/d( E1 ) hemos calculado la sección ecaz de fotodisociación ( σγ ), y, a partir de ella, la de captura radiativa ( σrc ). Finalmente, hemos calculado la tasa de reacción como función de la temperatura ( R(T) ) para el estado fundamental y el estado excitado ( 5/2+ ). Este último presenta transiciones E1 de un estado ligado d ( l= 2 ) a estados del continuo p3/2 , f5/2 y f7/2 . Podemos observar al compararla con la del estado fundamental que solo es relevante para energías un orden de magnitud más altas. En consecuencia, en el rango de temperatura de los sistemas astrofísicos relevantes, como las estrellas AGB, será despreciable frente a la contribución del estado fundamental. Finalmente podemos comparar las resultados de σrc a partir de datos experimentales del artículo de Nakamura [24] (ruptura coulombiana) y el de Reifarth [25] (medida direc- 34 ta), encontrando un acuerdo razonable. Cabe destacar que este procedimiento seguido con el cálculo numérico es general para núcleos halo de un neutrón, como son el 19 C o el 11 Be, mientras que el cálculo analítico depende de la estructura de capas de cada caso y de las transiciones posibles. Por ejemplo, para el primer estado excitado del 15 C el cálculo requiere hacer dos integrales análogas a la del Anexo A, una para la transición d→p y otra para las transiciones d→f , lo que nos deja con funciones de Hankel esféricas y de Bessel esféricas de orden superior a la que se obtienen para el estado fundamental. 35 Anexos A. Funciones de onda: integral analítica La parte radial de la probabilidad de excitación reducida toma la forma hf|r|ii=Z∞ 0 dr 1 kr( sin(kr)−kr ·cos(kr) )×r×D e −kBr =D kZ∞ 0 dr [sin(kr)−kr ·cos(kr))] e −kBr =D kZ∞ 0 dr sin(kr) e −kBr−DZ∞ 0 dr cos(kr)r e −kBr, que podemos reescribir como hf|r|ii=D1 kIA−IB, (52) donde hemos denido IA=Z∞ 0 dr sin(kr) e −kBr, IB=Z∞ 0 dr cos(kr)r e −kBr. Integral A: tiene estructura de integral cíclica IA=Z∞ 0 dr sin(kr) e −kBr=Z∞ 0 dr sin(kr)d dr −1 kB e −kBr =−1 kBsin(kr) e −kBr +∞ 0−kZ∞ 0 dr cos(kr) e −kBr =k kBZ∞ 0 dr cos(kr)d dr − e −kBr kB =−k k2 Bcos(kr) e −kBr +∞ 0+Z∞ 0 dr (ksin(kr)) e −kBr =k k2 B{1−kIA}. Por tanto IA=k k2 B+k2. (53) 36 Integral B: podemos relacionarla con IA . IB=Z∞ 0 dr cos(kr)r e −kBr=Z∞ 0 dr d dr sin(kr) kr e −kBr =1 ksin(kr)r e −kBr +∞ 0−Z∞ 0 dr sin(kr) e kBr−kBr e −kBr =−1 kIA+kB kZ∞ 0 dr sin(kr)r e −kBr =−1 kIA+kB kZ∞ 0 dr d dr −1 kcos(kr)r e −kBr =−1 kIA+cos(kr)r e −kBr +∞ 0−Z∞ 0 dr cos(kr) e −kBr−kBr e −kBr =−1 kIA+kB k2Z∞ 0 dr cos(kr) e −kBr−kBIB. También podemos relacionar la integral que nos queda con IA Z∞ 0 dr cos(kr) e −kBr=Z∞ 0 dr d dr 1 ksin(kr) e −kBr =1 ksin(kr) e −kBr ∞ 0−Z∞ 0 dr sin(kr)−kB e −kBr =kB kZ∞ 0 dr sin(kr) e −kBr=kB kIA, y sustituyendo IB=−1 kIA+k2 B k3IA−k2 B k2IB⇔IB=k2 B−k2 k2 B+k2 1 kIA. (54) Finalmente, a partir de las ecuaciones (52), (53) y (54) obtenemos la parte radial hf|r|ii=2Dk2 (k2+k2 B)2. (55) B. THOx El programa THOx nos permite calcular numéricamente las soluciones del sistema, por lo que podemos comparar estas soluciones más generales con la solución analítica (con aproximación asintótica). Para trabajar con él usamos el terminal de Linux y un archivo de entrada, o input , de la forma archivo.in , en el que asignamos las características del sistema, el mallado en r de las funciones de onda generadas, el tipo de potencial, etc. 37 Figura 13: Características del sistema. En el apartado & SYSTEM se especican las propiedades básicas del sistema neutrón+ core y en & corestates el estado del core ( 14 C). Para trabajar de forma ordenada, hemos optado por hacer un input con toda la información que vamos a usar, de forma que para hacer una transición concreta comentamos (símbolo !) o descomentamos la información relevante. Además, hay que tener en cuenta que las unidades en las que trabaja el programa son: L→ fm , E → MeV , m → uma . B.1. Input Podemos dividir las partes que modicamos del input en cuatro: Características del sistema. Características del cálculo numérico. Estados del sistema. Características de las funciones de onda del continuo. Características del sistema: en la Figura 13 vemos que en esta parte especicamos Z y A tanto del core como del nucleón de valencia (en nuestro caso será 14 C y un neutrón). Además, también jamos los espines de cada uno, y en el caso del core también su paridad. Para el neutrón tenemos: Z= 0 , A= 1 y su espín sn = 1/2 . Como el core es de 14 C, Z= 6 y A= 14 , y como tenemos todos los niveles llenos, su espín será 0+ . 38 Figura 14: Características del cálculo numérico. En & GRID se especica el mallado en r de las funciones de onda radiales generadas, y en & POTENTIAL el tipo y las propiedades del potencial. Características del cálculo numérico: en la Figura 14 vemos que en esta parte especicamos por un lado el mallado, que será ( rmin , rmax ) con pasos dr , y las propiedades del potencial. En nuestro caso hemos tomado un potencial de Woods-Saxon ( ptype=1 ), del que tenemos que especicar la difusividad a0 , el radio reducido r0 (ecuación (1)) y el valor de vl0 que, como se comenta en el trabajo, extraemos de la bibliografía y en general varía con l . Estados del sistema: en la Figura 15 vemos que en esta parte se especican todos los estados relevantes para el estudio de las transiciones E1 del estado fundamental ( 1/2+ ) y el primer excitado ( 5/2+ ). Para ello tenemos que especicar fundamentalmente: Propiedades: espín del estado JTOT , su paridad PARITY y el momento angular orbital lmax . Rango de búsqueda: dadas las propiedades de los estados se buscan en un cierto rango de energía ( exmin , exmax ). Es por esto que cuando buscamos estados del continuo tomamos exmin =0, ya que no existen estados del continuo por debajo, y seguimos la misma lógica con los estados ligados, en los que tomamos emax =0. Características de las funciones de los estados del continuo: en la Figura 16 vemos que en esta parte se especican las energías de los estados del continuo que se toman para el cálculo. Esto es necesario debido a que hay innitos con las mismas características, por lo que sería inviable realizar el cálculo si no acotamos la elección. Para ello 39 Figura 15: Estados del sistema. Se especifícan todos los estados relevantes para el estudio, primero los ligados y más abajo los del continuo. La propiedades fundamentales que se especican para cada uno son el espín ( JTOT ), la paridad ( PARITY ) y el momento angular orbital ( lmax ). Figura 16: Características de las funciones de los estados del continuo. En & SCATWF se especican las energías del continuo que se consideran para imprimir las funciones de onda, y con & belambda las energías para imprimir otros resultados, como B(E1) , σγ y σrc . 40 [27] "THOx : A few-body coupled-channels code with core and target excitations", A. M. Moro, J. A. Lay, M. Gómez-Ramos y R. de Diego, Universidad de Sevilla. Sin publicar. [28] Computer code fresco , I. J. Thompson,Comp. Phys. Rep. 7, 167(1988) [29] Busso, M., Gallino, R., & Wasserburg, G. J. (1999). Nucleosynthesis in asymptotic giant branch stars: Relevance for Galactic Enrichment and Solar System Formation. Annual Review of Astronomy and Astrophysics, 37, 239309. https://doi.org/10.1146/annurev.astro.37.1.239 47