TREBALL DE FI DE MÀSTER
DETERMINACIÓ DE LES CARACTERÍSTIQUES DE MODES NORMALS A PARTIR DE SIMULACIONS NUMÈRIQUES MITJANÇANT L’ANÀLISI CEOF: ABSORCIÓ RESSONANT
Biel Enric Castell Cladera
Màster Universitari Física Avançada i Matemàtica Aplicada (Especialitat/Itinerari d’Astrofísica i Relativitat)
Centre d’Estudis de Postgrau
Any Acadèmic 2019-20
DETERMINACIÓ DE LES CARACTERÍSTIQUES DE MODES NORMALS A PARTIR DE
SIMULACIONS NUMÈRIQUES MITJANÇANT L’ANÀLISI CEOF: ABSORCIÓ RESSONANT Biel Enric Castell Cladera
Treball de Fi de Màster
Centre d’Estudis de Postgrau
Universitat de les Illes Balears
Any Acadèmic 2019-20
Paraules clau del treball:
Física solar, bucle coronal, oscil·lacions transversals, modes normals, anàlisi CEOF.
Nom Tutor/Tutora del Treball Ramón Julio Oliver Herrero
Índex
1 Introducció 3
1.1 El sol i l’atmosfera solar . . . 3
1.2 Bucles coronals i les seves oscil·lacions transversals . . . 5
1.3 Equacions de la magnetohidrodinàmica. . . 7
1.3.1 Deducció de les equacions de la MHD . . . 7
1.3.2 Equacions de la MHD ideal . . . 10
1.3.3 Aproximacióβ = 0 . . . 11
1.3.4 Ones de la MHD ideal . . . 11
1.4 Estructura de la tesis. . . 12
2 Model d’equilibri, equacions linealitzades i modes normals 13 2.1 Equilibri . . . 13
2.2 Equacions linealitzades de la MHD ideal . . . 14
2.3 Modes normals . . . 15
2.3.1 Perfil de densitat discontinu . . . 16
2.3.2 Perfil de densitat continu . . . 18
3 Simulacions numèriques d’evolució temporal 20 3.1 Perfil de densitat discontinu . . . 20
3.2 Perfil de densitat continu . . . 22
4 Anàlisi "Complex Empirical Orthogonal Functions" 24 4.1 Introducció a l’anàlisi CEOF . . . 24
4.2 Exemple (Terradas et al. (2004)) [9] . . . 25
4.3 Aplicació del CEOF a les simulacions numèriques d’evolució temporal . . . 27
4.3.1 Perfil de densitat discontinu . . . 27
4.3.2 Perfil de densitat continu . . . 30
5 Mètode iteratiu 33 5.1 Perfil de densitat discontinu . . . 34
5.2 Perfil de densitat continu . . . 36
5.3 Temps d’esmorteïment . . . 40
6 Errors i ritme de convergència 41 6.1 Errors de les funcions pròpies . . . 41
6.2 Ritme de convergència del mètode iteratiu . . . 42
7 Conclusions i treball futur 44
Resum
Aquest TFM és una continuació de Rial et al. (2019, ApJ, 876, 86) [1]. En aquest article s’han extret les propietats d’alguns modes normals d’un sistema a partir de simulacions numèriques dependents del temps. El sistema consisteix en una llosa uniforme de plasma, és a dir, una làmina fina uniforme de plasma il·limitada en dues direccions ortogonals i incrustada en un entorn uniforme. Les simulacions numèriques resolen les equacions de la magnetohidrodinàmica i descriuen l’evolució temporal del sistema després de provocar una pertorbació inicial. Aquesta evolució temporal és igual a la suma de les evolucions temporals dels modes normals excitats per la pertorbació inicial. El mètode utilitzat per determinar les funcions dels modes normals és l’anàlisi de les funcions ortogonals empíriques complexes (CEOF) que, una vegada aplicat a les simulacions, proporciona una aproximació de les funcions pròpies i de la freqüència pròpia dels modes normals excitats. Per millorar aquesta aproximació, es selecciona un mode normal i s’utilitza com a condició inicial per una nova simulació numèrica, a la qual es pot tornar a aplicar l’anàlisi CEOF.
El fet de repetir aquest procés condueix a una millora de les característiques del mode normal obtingut.
En aquest treball substituïm el perfil de densitat discontinu a les fronteres de la làmina per una transició suau des de la llosa cap a l’entorn. Aquesta modificació fa aparèixer l’absorció ressonant, a partir de la qual l’energia dels modes normals amb velocitat transversal és transferida als modes normals amb velocitat paral·lela a les fronteres de la llosa. L’objectiu d’aquest TFM és provar el mètode de Rial et al.
(2019) [1] en aquest problema amb un espectre de modes normals més ric. A més a més, ara els modes normals no posseeixen una amplitud constant com a Rial et al. (2019) [1]. Hem pogut comprovar que el mètode funciona per al perfil de densitat continu i ens serveix per trobar bones aproximacions a les funcions pròpies i valors propis, la qual cosa inclou tant la freqüència com el ritme d’esmorteïment. El ritme de convergència dels errors és inferior al de Rial et al. (2019) [1]. Per aquest motiu, per obtenir una precisió determinada de les funcions pròpies són necessàries més iteracions.
1 | Introducció
1.1 El sol i l’atmosfera solar
El Sol és un objecte de gran bellesa i fascinació que ha estat estudiat amb interès durant milers d’anys.
Va néixer a partir d’un núvol interestel·lar en contracció fa uns 4600 milions d’anys. Aquesta protoestrella va començar en un estat on la gravetat era lleugerament superior al gradient de pressió, fet que provocava una contracció lenta i contínua que pot a poc encalentia el plasma. A un cert punt, la temperatura del nucli va ser suficientment alta per fusionar hidrogen a heli i aturar la contracció. Actualment el Sol es troba en aquest estat d’equilibri, aproximadament a la meitat del camí de la seva vida. D’aquí uns 5000 milions d’anys, quan l’hidrogen del nucli s’acabi, continuarà la fusió d’hidrogen al voltant del nucli d’heli, la qual cosa provocarà una expansió de les capes més externes del Sol i aquest es convertirà en una gegant vermella. Finalment, el nucli nu d’aquesta gegant vermella quedarà com una nana blanca, contenent la major part de la massa original en un volum similar al de la Terra.
El Sol és l’estrella més propera a la Terra i el cos principal del Sistema Solar, contenent un 99.8%
de la seva massa. Té una temperatura superficial aproximadament de 5500 K i és l’encarregat de proveir energia a la Terra i per tant permetre la vida tal i com la coneixem. És una esfera quasi perfecta, lleuger- ament el·lipsoidal degut a la seva rotació, que actualment conté un 74% d’hidrogen, un 25 % d’heli i petites quantitats d’altres elements més pesats. A la Taula 1 es poden veure alguns paràmetres típics del Sol.
Paràmetres Valors
Radi (R) 6.96×108 m
Massa (M) 1.99×1030 kg
Densitat mitjana 1.408 kg m−3 Gravetat superficial 273.95 m s−2 Lluminositat (L) 3.827×1026 W Període de rotació equatorial 25.38 dies Distància a la Terra (1 AU) 1.496×1011 m
Taula 1.1: Paràmetres físics del Sol
L’interior del Sol està estructurat en capes concèntriques. La frontera física i les característiques químiques entre les diferents capes són difícils d’establir degut a que no podem observar directament l’interior, però els models actuals aconsegueixen explicar satisfactòriament la major part dels processos observats. Les regions que podem diferenciar a l’interior del Sol són:
• Nucli: S’estén des del centre fins aproximadament 0.25 R. La densitat del nucli és de 1.5× 105kg m−3 i la temperatura és de 15×106K. És la zona on es genera l’energia a partir de la fusió d’àtoms d’hidrogen per formar àtoms d’heli. Aquesta energia posteriorment es transfereix a les capes superiors. Els fotons generats al nucli presenten un recorregut lliure mitjà molt més curt que el radi del nucli solar. Aquests fotons són absorbits i reemesos moltes vegades abans de sortir a través de la superfície solar. Per aquest motiu no ens arriba informació de com va ser creat aquest fotó ni de quines característiques tenia el medi.
• Zona Radiativa: S’estén des de r= 0.25 R fins aproximadament r = 0.7 R. En aquesta zona la densitat i la temperatura decreixen cap a l’exterior. L’energia és transportada per radiació. Els fotons són absorbits i reemesos contínuament per protons o nuclis d’Heli. Durant el viatge dels
fotons cap a l’exterior, de aproximadament106 anys, aquests canvien la seva longitud d’ona des de la zona de rajos x cap a zones de longitud d’ona més altes.
• Zona Convectiva: S’estén des der= 0.7R fins a la superfície solar. El gas de la base de la zona convectiva s’encalanteix, passa a ser menys dens i puja. Es troba amb gas més fred, intercanvia energia, es refreda, i torna a baixar. L’energia tèrmica és transportada cap a dalt a través de moviments de massa (convecció). Aquests moviments de material s’organitzen en forma de cel·les convectives.
Figura 1.1: Fotosfera solar. Esquerra: Imatge del disc solar complet. Font: https://sohowww.
nascom.nasa.gov/gallery/images/mdi20031028.html. Dreta: Imatge d’un grup de taques solars, on la granulació també es visible, de la regió activa 10030, feta pel 1m-Swedish Solar Telescope. Font:
https://www.nasa.gov/mission_pages/stereo/multimedia/TRACE_images.html.
Figura 1.2: Corona solar. Esquerra: Imatge del Sol durant un eclipsi. Es pot veure de color blanc la corona. Font: https://skyandtelescope.org/online-gallery/solar-corona-2/. Dreta: Imatge presa des del satèl·lit espacial SOHO amb un filtre centrat a 171 Å on es poden distingir zones brillants anomenades regions actives, que són regions que es troben damunt grups de taques solars i a les quals tenen lloc fenòmens com explosions violentes i ejeccions de plasma coronal. Font: https://sci.esa.
int/web/soho/-/47707-active-regions-and-magnetic-loops-in-the-sun.
Per sobre de la zona convectiva comença l’atmosfera solar. Degut a que el material es troba en estat de plasma, és difícil diferenciar exactament on comença. Podem diferenciar-hi quatre zones:
• Fotosfera: És la zona visible del Sol i la primera capa de l’atmosfera solar (veure figura1.1, esquerra).
Es situa per sobre de la zona convectiva i té una amplada d’uns 550 km i una temperatura de 5500 K. En el rang visible es poden diferenciar petits grànuls en continu moviment, els quals són la part superior de les cel·les convectives (figura 1.1, dreta). Per altra part també podem trobar taques
solars amb una temperatura d’uns 4000 K (més freda que la de la fotosfera). Aquestes taques són fenòmens magnètics causats per l’emergència de flux magnètic de l’interior del Sol. Al centre d’aquestes taques solars el camp magnètic pot prendre valors de fins a 3000 G. A la part superior de la fotosfera la temperatura baixa fins als 4000 K.
• Cromosfera: Damunt la fotosfera trobam una regió menys densa i més transparent anomenada Cromosfera, amb un grossor d’uns 2000 km. En aquesta part de l’atmosfera la temperatura creix amb l’altura, arribant a un valor d’uns 10000 K.
• Zona de transició: Per sobre de la cromosfera existeix una capa molt fina on la temperatura passa sobtadament de 10000 K a106K. La causa d’aquest fenomen és desconeguda. A la zona de transició la temperatura és prou elevada com per que el plasma estigui completament ionitzat, mentre que a la fotosfera i cromosfera està parcialment ionitzat.
• Corona: És la capa més externa de l’atmosfera solar. La llum provinent de la fotosfera fa que sigui difícil veure la corona. Només durant un eclipsi (figura1.2, esquerra) o mitjançant tècniques coronogràfiques és possible veure-la a ull nu. La densitat mitjana d’electrons és de aproximadament 1014m−3però pot ser de fins a 20 vegades superior en algunes zones. Aquesta densitat disminueix ràpidament a la superfície, arribant a un valor de 1010 m−3 a 10R. El plasma, completament ionitzat, pot arribar a temperatures de fins a 107 K. Degut a aquesta alta temperatura la corona emet sobretot radiació ultravioleta (figura1.2, dreta).
La figura1.3mostra la variació de la temperatura i la densitat en un model simplificat d’atmosfera solar.
Es pot apreciar com la densitat decreix amb l’altura i que la temperatura, de manera contraintuïtiva, comença a créixer a partir d’uns 500 km d’altura.
Figura 1.3: Temperatura i densitat a l’atmosfera solar. La línia contínua representa la temperatura mentre que la discontínua representa la densitat. Podem veure com la densitat disminueix amb l’altura, mentre que la temperatura presenta un augment sobtat de dos ordres de magnitud a la regió de transició, entre la cromosfera i la corona. Font: Priest (2013) [3].
1.2 Bucles coronals i les seves oscil·lacions transversals
Els bucles coronals són estructures brillants i corbades que apareixen com a arcs a la corona solar (figura 1.4). Els bucles coronals es troben principalment a regions actives, però també poden connectar una regió activa amb una altra. Als bucles coronals el plasma flueix seguint les fortes línies de camp magnètic,
Figura 1.4: Bucles coronals des del satèl·lit TRACE/NASA. Font: https://aasnova.org/2019/05/21/
its-getting-hot-hot-hot-with-coronal-loops/.
donant a lloc les seves formes característiques.
Normalment, però no sempre, els bucles coronals estan associats a les taques solars. Aquestes taques solars són en certa manera manifestacions visibles de zones on camps magnètics molt forts “atravessen” la fotosfera (superfície) estenent-se cap a l’atmosfera del Sol. Generalment trobam parells de taques solars amb pols magnètics oposats a les quals estan arrelats aquests arcs, que es poden arribar a estendre milers de kilòmetres per sobre la fotosfera cap a la corona i arribar a unes temperatures de més d’un milió de graus.
Els bucles coronals són més comuns al màxim del cicle solar, quan el camp magnètic del Sol està molt alterat i les taques solars són numeroses. Són més fàcils de veure en les regions de rajos X i ultravioleta de l’espectre electromagnètic.
Un dels principals problemes que els astrofísics intenten resoldre actualment és per què l’atmosfera del Sol és molt més calenta que la superfície del Sol (Figura1.3). Una de les possibles explicacions és degut a que d’alguna manera l’energia atrapada en els camps magnètics causa aquest escalfament. Els bucles coronals són probablement una de les claus en aquesta investigació.
L’any 1998 el satèl·lit Transition Region and Coronal Explorer (TRACE) va registrar per primera veg- ada oscil·lacions transversals de bucles coronals en una zona activa (Figura 1.5). L’amplitud més gran
Figura 1.5: Oscil·lacions transversals de bucles coronals. Esquerra: representació esquemàtica de l’amplitud d’oscilació transversal d’un bucle coronal. La component dominant del camp magnètic és en la direcció del bucle mentre que les oscil·lacions són perpendiculars a la direcció del camp magnètic.
Dreta: evolució temporal del desplaçament de l’àpex del bucle començant el dia 14 de Juliol de 1998 a les 13:13:51 UT. S’observa un esmorteïment de les oscil·lacions després d’uns pocs períodes d’oscil·lació.
Font: Nakariakov et al. (1999) [4].
d’oscil·lació es va veure a l’àpex (vèrtex de l’arc), mentre que els peus del bucle es mantengueren fixes.
Aquests tipus d’oscil·lacions són pràcticament harmòniques i fortament esmorteïdes, amb uns períodes típics d’oscil·lació de l’ordre dels minuts i un esmorteïment que dura normalment menys de 10 períodes.
Les oscil·lacions transversals de bucles coronals es poden interpretar com ones magnetohidrodinàmiques atrapades dins d’una estructura de plasma. Un dels principals models per a representar els bucles coro- nals ha estat fins el moment un tub cilíndric amb una densitat interiorρi (plasma que forma el bucle) i una densitat exteriorρe (medi que enrevolta el bucle). Adoptar aquest model implica ignorar l’efecte de la curvatura de les línies de camp magnètic, que ha resultat ser una aproximació bastant bona. Així i tot, també s’ha emprat sovint el model de llosa, en geometria cartesiana, que consisteix en una làmina infinita (el bucle coronal) submergida en un medi infinit (la corona al voltant del bucle). Aquest segon model és el que emprarem en aquest treball. Tot i ser més simple que el tub cilíndric ens permetrà explicar el mecanisme físic que en aquests moments es creu que causa l’esmorteïment: l’absorció ressonant. Al capítol 2 explicarem el model i en què consisteix l’absorció ressonant. Abans, però, convé introduir les equacions de les ones magnetohidrodinàmiques.
1.3 Equacions de la magnetohidrodinàmica
La magnetohidrodinàmica (MHD) és la teoria de la física que descriu el comportament macroscòpic d’un plasma. Les equacions de la MHD es poden entendre com una combinació de les equacions de la dinàmica de fluids juntament amb les equacions de Maxwell. Degut a que el Sol es troba en un estat de plasma, tenim un fluid (el gas) que està ionitzat i per tant l’electromagnetisme juga un paper important.
1.3.1 Deducció de les equacions de la MHD
Per fer la derivació de les equacions de la MHD en aquesta secció seguim https://www.docsity.com/
es/deduccion-de-ecuaciones-mhd/4169393/. Consideram una cel·la en una dimensió d’amplada∆x.
La massa a l’interior de la cel·la pot canviar degut a la diferència de la massa que entra i surt de la cel·la, que anomenam F(x). La massa a l’interior de la cel·la és
m=
Z x+∆x x
ρdx, (1.1)
onρ(x, t)és la densitat de la cel·la. Si volem saber com varia la massa en el temps
∂
∂t
Z x+∆x x
ρdx
!
=F(x)−F(x+ ∆x), (1.2)
∂ρ
∂t = lim
∆x→0
F(x)−F(x+ ∆x)
∆x , (1.3)
∂ρ
∂t +∂F(x)
∂x = 0. (1.4)
El flux de massa per segon a les fronteres de la cel·la és
F(x, t) =ρ(x, t)vx(x, t). (1.5)
En 3 dimensions podem generalitzar l’equació (1.4) a
∂ρ
∂t +∇ ·(ρ~v) = 0. (1.6)
Aquesta equació es coneix com a l’equació de continuïtat. Es pot aplicar un procediment anàleg per a qualsevol magnitud que es conservi, com per exemple, a la densitat de moment o la densitat d’energia.
Si R
~
u dxes conserva, llavors
∂~u
∂t +∇ ·F~ = 0. (1.7)
A la dinàmica de fluids, la relació entre la derivada total i la derivada parcial és d
dt = ∂
∂t+~v· ∇, (1.8)
on d/dt ens indica com varia una quantitat en un punt que es mou amb el fluid i ∂/∂t ens indica com varia una quantitat en un punt fixe de l’espai.
Si no tenim intercanvis de calor per conducció o transport, els canvis en la pressió i el volum d’un element de fluid (movent-se amb el fluid) són adiabàtics. En aquest cas es compleix que
pVγ =ct, (1.9)
onγ és el coeficient de dilatació adiabàtica. D’aquesta forma d
dt(pVγ) = 0. (1.10)
Si ens movem sobre un element de fluid la massa es conserva, i per tantV ∝ρ−1 d
dt p
ργ
= 0. (1.11)
Aquesta equació es coneix com a l’equació d’energia adiabàtica.
Procedint de forma anàloga a com hem fet per l’equació de continuïtat, el moment dins de la cel·la unidimensional és ρvx∆x. El moment total de la cel·la canvia degut al gradient de pressió. D’aquesta forma, l’equació que ens queda és
∂
∂t(ρvx∆x) =F(x)−F(x+ ∆x) +p(x)−p(x+ ∆x), (1.12) on hem suposat que altres forces, com per exemple la gravitatòria, són menyspreables. Ara F(x) és el flux de moment per segonF =ρvxvx. Fent el límit∆x→0 arribam a
∂
∂t(ρvx) +∂F
∂x =−∇p. (1.13)
Amb l’ajuda de l’equació de conservació de massa arribam a ρ∂vx
∂t +ρvx∂vx
∂x =−∇p. (1.14)
Utilitzant la regla de la cadena
dvx(x, t) dt = ∂vx
∂t +∂x
∂t
∂vx
∂x, (1.15)
ρdvx
dt =−∇p. (1.16)
Ampliant a tres dimensions ens queda
ρd~v
dt =−∇p. (1.17)
A un plasma, però, tenim partícules carregades. La força de Lorentz ens diu que la força sobre una partícula carregada en un camp electromagnètic és:
F~ =q(E~ +~v×B).~ (1.18)
La força electromagnètica total per unitat de volum sobre els electrons del plasma és−nee(E~ +~ve×B)~ i sobre els ions del plasma nie(E~ +~vi×B), on~ ne i ni són les densitats del nombre d’electrons i ions respectivament. Per tant, la força electromagnètica total per unitat de volum és
F~ =e(ni−ne)E~ + (eni~vi−ene~ve)×B.~ (1.19) Si el plasma és quasi-neutral (ni'ne) podem aproximar l’expressió anterior com
F~ =en(~vi−~ve)×B~ =~j×B,~ (1.20)
on~jés la densitat de corrent. D’aquesta forma obtenim l’anomenada equació de moment, afegint la força de Lorentz al mebre de la dreta de l’equació (1.17)
ρd~v
dt =−∇p+~j×B.~ (1.21)
Cal notar que de moment l’únic canvi en les equacions de la magnetohidrodinàmica respecte a les de fluids és el terme~j×B~ a l’equació de moment. Ara necessitam una equació per a la densitat de corrent i el camp magnètic per tancar el sistema. Per trobar-les partim de les equacions de Maxwell
∇ ·E~ = ρ 0
, (1.22)
∇ ·B~ = 0, (1.23)
∇ ×E~ =−∂ ~B
∂t , (1.24)
∇ ×B~ =µ0~j+0µ0
∂ ~E
∂t . (1.25)
L’equació (1.25) és l’anomenada llei d’Ampère. Començam fent un anàlisi dimensional per dos dels tres termes d’aquesta equació: el rotacional del camp magnètic i el terme de desplaçament de corrent.
∇ ×B~
∼ B0
L , (1.26)
0µ0
∂ ~E
∂t
∼0µ0
E0
T ∼ 1 c2
E0
T . (1.27)
on L, T, B0 i E0 corresponen a les magnituds característiques del plasma de longitud, temps, camp magnètic i camp elèctric, i c2 és la velocitat de la llum al quadrat. Més endavant mostrarem que a la magnetohidrodinàmicaE0/B0∼V, onV correspon a la magnitud de velocitat. Comparant aquests dos termes
0µ0∂ ~E/∂t
∇ ×B~ ∼ L T
V c2 ∼ V2
c2. (1.28)
Aquest anàlisi ens condueix a que per velocitats baixes o feqüències baixes podem ignorar el terme de desplaçament de corrent. Així, l’equació (1.25) queda
∇ ×B~ =µ0~j. (1.29)
A partir de l’equació de continuïtat, podem escriure les següents expressions per als ions i els electrons
∂ni
∂t +∇ ·(ni~vi) = 0, (1.30)
∂ne
∂t +∇ ·(ne~ve) = 0. (1.31)
Multiplicant cada una d’aquestes equacions per la càrrega corresponent i sumant obtenim
∂σ
∂t +∇ ·~j= 0, (1.32)
onσés la densitat de càrrega i~jés la densitat de corrent. Si consideram freqüències baixes, multiplicant cada costat de l’equació (1.29) per l’operador de divergència
∇ ·(∇ ×B) =~ µ0∇ ·~j, (1.33)
∇ ·~j= 0→ ∂σ
∂t = 0. (1.34)
Per tant, per processos de baixa freqüènciani'ne. Això es coneix com quasi-neutralitat.
Si consideram un fluid conductor, com el plasma, que es mou a una velocitat ~v en un camp magnètic
B, el moviment relatiu del fluid indueix un camp elèctric~ E~ que exerceix una força elèctrica sobre les partícules carregades, donant a lloc un corrent elèctric~j. L’equació que descriu aquesta situació és la coneguda llei d’Ohm generalitzada
mene
d~ve
dt =−nee ~E+nemeν(~vi−~ve)−ene~ve×B,~ (1.35) ne,mei~vesón la càrrega, la massa i la velocitat dels electrons respectivament, iνés la freqüència de xoc entre els electrons i els ions. Degut a que els electrons tenen una massa molt més petita que la dels ions podem menysprear el terme de l’esquerra de l’equació. Per altra part, com que els electrons es mouen molt més ràpid que els ions (~ve~vi) l’equació anterior ens queda
0 =−nee ~E+nemeν~ve−ene~ve×B,~ (1.36) 0 =−nee(E~ +~ve×B) +~ meν
e ~j, (1.37)
E~ +~v×B~ = 1
σ~j, (1.38)
on σ = nee2/νme és la conductivitat elèctrica. Escrivint l’equació anterior en funció de la resistivitat elèctricaη=σ−1obtenim
E~ +~v×B~ =η~j. (1.39)
Amb tot això, el sistema d’equacions de la magnetohidrodinàmica resultant és el següent:
∂ρ
∂t +∇ ·(ρ~v) = 0, (1.40)
ρd~v
dt =−∇p+~j×B,~ (1.41)
d dt
p ργ
= 0, (1.42)
∂ ~B
∂t =−∇ ×E,~ (1.43)
~j= 1
µ0∇ ×B,~ (1.44)
E~ +~v×B~ =η~j. (1.45)
Cal recordar que degut a que hem assumit quasi-neutralitat hem de tenir freqüències baixes i velocitats molt menors a la velocitat de la llum. Per altra banda, comentar que no tenim en compte viscositat, transferència de calor, conducció, radiació, gravetat, rotació, ni ionització.
1.3.2 Equacions de la MHD ideal
Les equacions generals de la magnetohidrodinàmica (1.40) – (1.45) que hem vist a la secció anterior es simplifiquen a les equacions ideals de la magnetohidrodinàmica si consideram un sistema sense resistivitat (η = 0). En aquest cas podem substituir el camp elèctric de l’equació (1.45) dins l’equació (1.43), i el corrent elèctric (equació (1.44)) dins l’equació (1.41). Si a més a més reescrivim les derivades parcials en forma de derivades totals obtenim
dρ
dt =−ρ∇ ·~v, (1.46)
ρd~v
dt =−∇p+ 1
µ(∇ ×B)~ ×B,~ (1.47)
dp dt =γp
ρ dρ
dt, (1.48)
∂ ~B
∂t =∇ ×(~v×B),~ (1.49)
p=ρRT
˜
µ. (1.50)
Aquestes equacions juntament amb la condició solenoidal ∇ ·B~ = 0 formen les equacions ideals de la MHD. Les variablesρ,T,p,~v,B,~ R,µ,µ˜ són respectivament la densitat, la temperatura, la pressió del gas, la velocitat, el camp magnètic, la constant de gas ideal, la permitivitat magnètica i el pes atòmic mitjà.
1.3.3 Aproximació β = 0
Una gran varietat de les estructures solars que observam pareixen romandre sense moviment durant llargs períodes de temps. Per això, el seu estudi es pot realitzar mitjançant solucions d’equilibri magneto- hidrostàtic a les equacions de la MHD. Exemples d’aquests equilibris magnetohidrostàtics són les arcades coronals, les taques solars, les protuberàncies solars, o com en el cas d’aquest treball, els bucles coronals.
Partim de l’equació (1.47) i vem que el terme de l’esquerra pot ser menyspreat sempre que la velocitat del flux sigui suficientment baixa. Fent un anàlisi adimensional a aquesta equació tenim
ρV T = p
L+B2 µL → V
T = p Lρ+ B2
µρL →V2= p ρ+B2
µρ. (1.51)
Definim la velocitat del so del plasma com vs= (γp/ρ)1/2 i la velocitat d’Alfvén comvA= (B2/µρ)1/2. D’aquesta manera
v2= vs2
γ +v2A. (1.52)
En els bucles coronalsvsivAprenen valors típics de l’ordre de 200 km s−1i 1000 km s−1, respectivament, mentre que les velocitats del plasma són generalment molt menors. Sempre que la velocitat del flux sigui molt més petita que la velocitat del so i que la velocitat d’Alfvén podem menysprear el terme de l’esquerra.
En aquest punt podem pensar en comparar els termes del gradient de pressió i de la força de Lorentz.
Definim β (beta del plasma) com el quocient entre la pressió tèrmica i la pressió magnètica.
β' |∇p|
~j×B~
= |∇p|
jB =
p L B
µLB = p
B2/µ. (1.53)
Si β =µp/B21qualsevol gradient de pressió és dominat per la força de Lorentz i podem considerar que aquesta és l’única força que queda a l’equació de continuïtat del moment.
~j×B~ = 0. (1.54)
A partir de les definicions de les velocitats del so i d’Alfvén veim que la beta del plasma se pot escriure β =v2s/γvA2.Pels valors típics d’aquestes velocitats mencionats abans,β ∼0.1i per tant en aquest treball utilitzarem aquesta aproximació (β = 0). Aquesta aproximació és bona per a aplicacions a la corona solar, on la pressió magnètica domina sobre la pressió del gas. La seva implicació és que podem triar lliurement els perfils de densitat del nostre model i per tant de velocitat d’Alfvén si el camp magnètic és conegut.
També cal dir que aquesta aproximació ens simplificarà de forma significativa l’anàlisi matemàtic de les oscil·lacions transversals del nostre model.
Aplicant l’aproximació β = 0, és a dir, menyspreant la pressió del gas, el sistema anterior es redueix
a ∂ρ
∂t =−∇ ·(ρ~v), (1.55)
ρ∂~v
∂t =−ρ(~v· ∇)~v+1
µ(∇ ×B)~ ×B,~ (1.56)
∂ ~B
∂t =∇ ×(~v×B).~ (1.57)
1.3.4 Ones de la MHD ideal
Si pensam en un medi homogeni que és pertorbat a una regió i a un instant donat, ens apareixeran diferents tipus d’ones que es propagaran en totes direccions. En absència de camp magnètic, tenim les ones acústiques clàssiques que es propaguen isotròpicament en el plasma a la velocitat del so del medi vs com a variacions de pressió i altres quantitats termodinàmiques. Però degut a que tenim un camp magnètic, les ones també es veuen sotmeses a forces magnètiques a part de a la pressió del plasma.
Aquestes tipus d’ones són les ones magnetohidrodinàmiques i en podem tenir de dos tipus: ones d’Alfvén i ones magnetoacústiques. Les ones d’Alfvén són incompressibles i es propaguen únicament degut a la tensió magnètica. Són ones no dispersives i viatgen en la direcció de les línies de camp magnètic a una velocitat constant vA anomenada velocitat d’Alfvén. Per altra banda les ones magnetoacústiques són ones compressibles que sorgeixen degut a la combinació de l’efecte del gradient de pressió i les forces
magnètiques. Podem tenir dos tipus d’ones magnetoacústiques: les ràpides i les lentes.
L’interior del Sol està en continu moviment degut als moviments convectius. A la superfície solar, on la β 1, aquests moviments del plasma arrosseguen els camps magnètics, canviant la seva complexi- tat topològica per damunt la superfície. Com a conseqüència, el camp magnètic de l’atmosfera del Sol emmagatzema energia contínuament fins arribar a un punt al qual l’excés d’energia és alliberat en esde- veniments catastròfics, com les flamarades solars, fortes explosions normalment acompanyades d’ejeccions de massa coronal. Aquests esdeveniments ocorren sobretot a zones properes a regions actives i per tant sovint aprop de bucles coronals. Tant el medi com les estructures que puguin enrevoltar les flamarades solars són pertorbats i es generen ones, com per exemple les oscil·lacions transversals de la secció1.2. A la següent secció descriurem el model que emprarem per al bucle coronal i linealitzarem les equacions per tal de poder introduir després una pertorbació inicial i veure com reacciona.
1.4 Estructura de la tesis
Al capítol2descriurem el model de bucle coronal que emprarem i linealitzarem les equacions de la MHD per tal d’estudiar els modes normals del bucle i per poder introduir les solucions analítiques del prob- lema. També veurem que segons com triem el perfil de densitat a l’equilibri podem tenir una resposta oscil·latòria periòdica o una resposta oscil·latoria esmorteïda degut al fenomen d’absorció ressonant. Al capítol3 resoldrem les equacions linealitzades a partir de simulacions numèriques emprant una condició inicial gaussiana, que farà la funció de depositar una energia inicial a la llosa, la qual cosa provocarà una excitació dels diferents modes normals d’oscil·lació. Podrem analitzar quins modes s’han excitat realitzant un espectre de potències per les freqüències. Al capítol 4 introduirem l’anàlisi CEOF, una tècnica multivariant que aplicada a les simulacions numèriques ens permetrà extreure una aproximació de la distribució espacial de diversos modes normals del sistema (en concret estudiarem el mode fona- mental). Al capítol 5explicarem en què consisteix el mètode iteratiu que empram en aquesta tesis. Les aproximacions a les funcions pròpies del mode fonamental extretes de l’anàlisi CEOF ens serviran com a condició inicial per a una nova simulació numèrica, provocant que el mode fonamental s’exciti amb més energia que a la primera simulació numèrica i per tant la resta de modes excitats ho facin amb menys energia. Podem pensar en repetir aquest procés tantes vegades com vulguem per tal d’obtenir una millor aproximació al mode fonamental després de cada iteració. Al capítol 6 calcularem els errors en funció de cada iteració entre les aproximacions extretes del CEOF per al mode fonamental d’oscil·lació i les solucions analítiques trobades al capítol 2(error global), i entre dues aproximacions successives extretes del CEOF per al mode fonamental (diferència global), per tal de poder calcular el ritme de convergència del procés iteratiu i comparar-lo amb els resultats obtinguts a Rial et al. (2019) [1].
2 | Model d’equilibri, equacions linea- litzades i modes normals
En aquest capítol presentam els dos models d’equilibri que emprarem per representar un bucle coronal, les equacions linalitzades de la MHD ideal en el límit β = 0per la nostra configuració i discutirem els seus modes normals. En el desenvolupament seguirem Arregui et al. (2007) [2] i Rial et al. (2019) [1].
2.1 Equilibri
Per estudiar les oscil·lacions dels bucles coronals es sol menysprear la seva curvatura, el que condueix a models basats en estructures cil·líndriques. Aquestes configuracions se poden simplificar encara més anant a estructures pla-paral·leles invariants en una direcció. Així, el model del bucle coronal que farem servir consisteix en una llosa amb un gruix de2aen la direccióx, infinit en la direccióyi amb una altura L en direcció z (veure figura2.1). El significat dels dos plansz =±L/2 s’explicarà a la secció 2.2. El sistema és invariant en direcció ya l’equilibri.
Figura 2.1: Esquema del model d’equilibri, on l’àrea ombrejada representa el bucle coronal limitat pels plansx=±aiz=±L/2. La configuració és infinita en direcció y. Tenim un camp magnètic uniforme a tot l’espai dirigit en la direcció z i una densitatρi a l’interior del bucle i ρe al medi que l’envolta. La figura de l’esquerra està agafada de Rial et al. (2019).
Tot i que aquest model és una simplificació excessiva i està lluny de ser una representació precisa de bucles coronals reals, la seva senzillesa ens permetrà un estudi detallat del nostre problema. Consider- arem dues situacions d’equilibri: la primera correspondrà a un perfil de densitat constant on la densitat no dependrà explícitament de x, sinó que tendrem una densitat interior ρi corresponent al plasma en repòs del bucle coronal, i una densitat exterior ρe que representa el medi que envolta al bucle (ambdues constants), de tal forma que:
ρ0(x) =
ρi , |x| ≤a, ρe , |x|> a.
(2.1)
Figura 2.2: Perfil de densitat continu amb dues capes de transició d’amplada l que connecten la densitat interior, ρi, amb la densitat exterior,ρe.
La segona situació consisteix en afegir dues capes de transició d’amplada l en la direcció x, que ens connectaran la densitat a l’interior de la llosa, ρi, amb la densitat a l’exterior, ρe, de forma contínua i suau (veure figura2.2). En aquesta situació tendrem el següent perfil de densitat:
ρ0(x) =
ρi , |x| ≤a−l/2, f(x) , a−l/2≤ |x| ≤a+l/2, ρe , |x| ≥a+l/2,
(2.2)
on
f(x) =ρi 2
1 +ρe
ρi
−
1−ρe ρi
sin
π|x−a|
l
. (2.3)
La idea és mostrar com amb un perfil de densitat continu es produeix l’efecte d’absorció ressonant (que explicarem més endavant) i obtenim oscil·lacions transversals esmorteïdes del bucle coronal, mentre que amb un perfil de densitat discontiu obtenim oscil·lacions transversals periòdiques sense esmorteïment.
Per altra banda tenim un camp magnètic uniforme a tot l’espai B~0 = B0ˆez en la direcció longitudi- nal del bucle, és a dir, en direcció z. Definim les velocitats d’Alfvén a l’interior i a l’exterior de la llosa com
vA(x) =
vAi≡ √Bµρ0
i , |x| ≤a, vAe≡ √Bµρ0
e , |x| ≥a.
(2.4)
2.2 Equacions linealitzades de la MHD ideal
Una vegada introduït l’estat d’equilibri (un camp magnètic uniforme a tot l’espai i els dos perfils de densitat que es tendran en compte, discontinu i continu en la direcció xi ambdós independents de z) passam a la linealització de les equacions de la MHD. Considerarem l’aproximació β = 0, donada per les equacions (1.55) – (1.57). En el nostre estudi ens volem concentrar en les ones de petita amplitud i per això escriurem la densitat, les components de la velocitat i les components del camp magnètic com la suma del valor que prenen a l’estat d’equilibri més una pertorbació. Degut a que a l’equilibri tenim el plasma quiet, les components de la velocitat són zero. D’aquesta forma podem escriure les variables com
ρ(~r) =ρ0(x) + ˆρ(~r), (2.5)
~v(~r) =~0 + ˆ~v(~r), (2.6)
B(~~ r) =B0ˆk+~b(~ˆ r). (2.7)
Substituint les expressions (2.5) – (2.7) dins les equacions (1.55) – (1.57) i eliminant els termes de segon ordre o superior de les variables pertorbades obtenim el següent sistema d’equacions:
∂ρˆ
∂t =−ρ0
∂vˆx
∂x −ikyvˆy
, (2.8)
∂ˆvx
∂t = B0
µρ0
∂ˆbx
∂z −∂ˆbz
∂x
!
, (2.9)
∂vˆy
∂t = B0
µρ0
∂ˆby
∂z +ikyˆbz
!
, (2.10)
∂vˆz
∂t = 0, (2.11)
∂ˆbx
∂t =B0∂ˆvx
∂z , (2.12)
∂ˆby
∂t =B0
∂ˆvy
∂z , (2.13)
∂ˆbz
∂t =−B0
∂ˆvx
∂t −ikyˆvy
. (2.14)
Cal comentar que s’ha assumit que les pertorbacions es propaguen en direccióy amb nombre d’onaky i per tant la dependència de les variables en aquesta direcció és del tipus exp(−ikyy). En aquest sistema d’equacions podem veure com la densitat pertorbada, ρ, només apareix a l’equació (2.8) i que pot serˆ calculada una vegada es coneguinˆvxivˆy. També veim que l’absència de forces en la direccióz condueix a que la velocitat en aquesta direcció sigui zero. En general, en el cas β 6= 0 aquesta component de la velocitat no s’anul·la . Per això, les equacions que ens interessa resoldre són les (2.9), (2.10), (2.12), (2.13) i (2.14). Recordam també queρ0 =ρ0(x) (equacions (2.1) i (2.2)) i que les variables són funcions que ara depenen de x,z it(ja que hem imposat com és la dependència amb y).
Les equacions diferencials que acabam de presentar s’han de complementar amb condicions de contorn.
Imposarem que la velocitat als plans z = ±L/2, que són anàlegs als peus del bucle coronal, és zero.
Aquesta condició ve imposada per les observacions, que mostren que els peus dels bucles estan fixos du- rant les oscil·laccions transversals. D’altra banda, el nostre sistema és infinit en la direccióxi per aquest motiu imposarem que les pertorbacions s’anul·len quanx→ ±∞.
2.3 Modes normals
Considerant que les funcions pròpiesvˆx,vˆy,ˆbx,ˆby iˆbz són funcions separables de la formaf(x)g(z)h(t), les equacions (2.9), (2.10), (2.12), (2.13) i (2.14) ens revelen que la dependència en z devˆx, vˆy iˆbz és cos(kzz)mentre que perˆbxiˆby éssin(kzz). Per satisfer les condicions de contorn als extrems de la llosa kz = (n+ 1)π/L, amb n= 0per al mode fonamental,n= 1per al primer mode excitat, etc. Per a un estudi dels modes normals imposam una dependència entde la formaexp(iwt). Per tant, les components pertorbades per a la velocitat i el camp magnètic tenen la següent forma
ˆ
vx(x, z, t) = ˆvx(x) cos(kzz)eiwt, (2.15) ˆ
vy(x, z, t) = ˆvy(x) cos(kzz)eiwt, (2.16) ˆbx(x, z, t) = ˆbx(x) sin(kzz)eiwt, (2.17) ˆby(x, z, t) = ˆby(x) sin(kzz)eiwt, (2.18) ˆbz(x, z, t) = ˆbz(x) cos(kzz)eiwt. (2.19) Subsituïnt les expressions (2.15) – (2.19) dins les equacions (2.9), (2.10), (2.12), (2.13) i (2.14) obtenim el següent sistema d’equacions:
wˆvx=−B0 µρ0
kz(iˆbx)− d dx(iˆbz)
, (2.20)
w(iˆvy) = B0
µρ0
h
kzˆby+ky(iˆbz)i
, (2.21)
w(iˆbx) =−B0kzvˆx, (2.22)
wˆby=B0kz(iˆvy), (2.23)
w(iˆbz) =−B0 dˆvx
dx −ky(iˆvy)
. (2.24)
Això no és més que un problema de valors propis i funcions pròpies on els valors propis són les freqüències wi les funcions pròpies sónvˆx(x),iˆvy(x),iˆbx(x),ˆby(x)iiˆbz(x). El fet de que algunes d’aquestes funcions pròpies vagin acompanyades del nombre imaginarii és únicament perquè aquestes siguin reals en el cas d’una densitat d’equilibri discontínua a les fronteres de la llosa (equació (2.1)). Cal destacar que les funcions pròpiesiˆbxiˆby són proporcionals avˆxiiˆvy respectivament, per tant no serà necessari calcular- les ja que es poden obtenir directament multiplicant pel factor corresponent. Recombinant les equacions anteriors es pot arribar a una equació diferencial de segon ordre on només apareix la funció pròpiavˆx:
vA2(x)(w2−k2zvA2(x)) w2−(k2y+k2z)vA2(x)
d2vˆx(x)
dx2 − ky2vA2(x)w2dρdx0(x) [w2−(k2y+k2z)vA2(x)]2ρ0(x)
dˆvx(x) dx +(w2−k2zvA2(x))ˆvx(x) = 0.
(2.25)
Com hem explicat a la secció 2.1, en aquest treball considerarem dos perfils diferents de densitat. A continuació analitzarem l’equació (2.25) per a les dues situacions que podem tenir a l’equilibri.
2.3.1 Perfil de densitat discontinu
La primera situació correspon a ρ0(x) = ct. a cada una de les tres regions del sistema (equació (2.1)).
En aquest cas la densitat presenta una discontinuïtat a x = ±a, concretament passam de tenir una densitatρia l’interior de la llosa a una densitatρea l’exterior. Així i tot, aquestes dues densitats interior i exterior són constants, de tal forma que dρ0(x)/dx = 0. Tenint això en compte, l’equació (2.25) es simplifica considerablement:
d2vˆ1x
dx2 =m2vˆ1x, (2.26)
on m2 = k2y+kz2−w2/vA2. La quantitat m2 pren els valors m2i i m2e a l’interior i exterior de la llosa, respectivament, ja que les velocitats d’Alfvén depenen de la densitat i tenen valors diferents a l’interior i a l’exterior de la llosa: vAi ivAe. Una vegada resolta l’equació (2.26) podem obteniriˆvy iiˆbz a partir de
iˆvy = ky
m2 dˆvx
dx, (2.27)
iˆbz B0
=−1 w
κ2 m2
dˆvx
dx, (2.28)
onκ2=k2z−w2/vA2 =m2−k2y. Altra vegada tenim dos possibles valors perκ,κi iκedepenent de si ens trobem a l’interior o a l’exterior de la llosa.
Per resoldre l’equació (2.26) hem d’imposar en la direccióxcondicions de contorn als extrems del sistema, concretament ˆvx → 0 quan x→ ±∞. A més, imposarem condicions de continuïtat al salt de densitat que tenim a x=±a, que segons Goedbloed et al. (2004) [5] són la continuïtat de la componentvˆx i de la pressió total, que alhora implica la continuïtat deiˆbz. Obtenim la següent solució:
ˆ vx(x) =
Cexp(mex) , x <−a, Acosh(mix) , −a≤x≤a, Cexp(−mex) , x > a.
(2.29)
Imposant la condició de continuïtat als límits de la llosa,vˆx(±a−) = ˆvx(±a+), obtenim
C=Aexp(mea) cosh(mia). (2.30)
Ara que ja coneixem l’expressió per avˆx, a partir de l’equació (2.28) trobam:
iˆbz(x) B0
=−1 w
Cexp(mex)κ2e/me , x <−a, Asinh(mix)κ2i/mi , −a≤x≤a,
−Cexp(−mex)κ2e/me , x > a.
(2.31)
Imposant la condició de continuïtat als límits de la llosaˆbz(±a−) = ˆbz(±a+), obtenim:
C=−Aexp(mea) sinh(mia)κ2i κ2e
me
mi. (2.32)
La constant A es pot triar arbitràriament perquè ens trobam en el límit ideal, així que triam A = 1.
D’aquesta manera queden obtingudes analíticament les funcions pròpies ˆvx, iˆvy iiˆbz. Recordam queiˆbx
iˆbysón proporcionals a les funcions pròpies anteriors (equacions (2.22) i (2.23)). Per aquest motiu no es mostraran a les gràfiques aquestes dues funcions pròpies durant tot el treball. A la figura 2.3 es poden veure les solucions de les funcions pròpies per al mode fonamental (n= 0). Els valors dels paràmetres que s’han emprat sónρi/ρe= 10,L= 50a, ikya= 0.5. Els valors adimensionals s’obtenen amb l’ajuda de la longitud a, la velocitatvAi i el tempsτAi=a/vAi. Aquests valors dels paràmetres es mantindran durant tot el treball.
Figura 2.3: Cas amb perfil de densitat discontinu. Representació de les funcions pròpiesvˆx,iˆvy iiˆbz del mode fonamental lateralment evanescent. Les línies discontínues representen els límits de la llosa. Font:
Rial et al. (2019)[1].
En quant als valors propis, a partir de les equacions (2.30) i (2.32) trobam la relació de dispersió tanh(mia) =−κ2e
κ2i mi
me
, (2.33)
la qual ens permet trobar els valors propis del problema, és a dir, les freqüències pròpies. La freqüència per al mode fonamental és w/τAi = 0.1011, per al primer excitat tenim w/τAi = 0.2989 i per al segon excitat w/τAi= 0.4852.
Amb el perfil de densitat discontinu a l’equilibri hem pogut trobar solucions analítiques al problema.
Això serà de gran utilitat més endavant ja que podrem comparar els resultats obtinguts a les simulacions numèriques tant amb les funcions pròpies com amb els valors propis. Un detall important és que a l’hora d’obtenir els resultats numèricament emprarem un domini finit, i fins ara hem imposat la condició de que les funcions pròpies siguin evanescents, és a dir que se’n vagin a zero a x = ±∞. Més endavant explicarem com serà el tractament numèric, però de moment avançam que el domini de computació en direccióxserà−20a≤x≤20a . El fet de tenir aquests límits finits ens farà aparèixer uns nous modes d’oscil·lació no evanescents, que anomenarem confinats. Aquests modes, tot i no ser una solució física si no una solució deguda al domini computacional, també s’excitaran amb la pertorbació inicial i conviuran durant un temps amb els modes evanescents descrits abans. Repetint el procés analític anterior però amb la condició de que les funcions pròpies s’anul·len a x= ±20a enlloc de a x = ±∞, obtenim uns nous valors propis i funcions pròpies per al problema. Ens interessa sobretot el valor propi del mode fonamental lateralment confinat: w/τAi= 1.676. Per a més informació, veure Rial et al. (2019) [1].
2.3.2 Perfil de densitat continu
De moment hem analitzat el cas amb perfil de densitat discontinu a l’equilibri (equació (2.1)) a partir del qual l’equació (2.25) s’ha simplificat considerablement i hem pogut trobar solucions analítiques al problema de valors propis i funcions pròpies. De manera similar, en el cas del perfil de densitat continu (equació (2.2)), tot i que no podem trobar solucions analítiques per a les funcions pròpies i els valors propis, sí que podrem resoldre el sistema numèricament i trobar les funcions pròpies i valors propis cor- responents (veure Arregui et al (2007) [2]).
En aquest segon cas és interessant introduir el concepte d’absorció ressonant. Si ens fixem en l’equació (2.25), els denominadors que acompanyen els termes amb derivades es fan zero quan
w2(x) = (ky2+k2z)v2A(x). (2.34) Definim la freqüència d’Alfvén com wA(x) = (ky2+k2z)1/2vA(x). Degut a que la densitat ρ0(x) depèn de la coordenada x, la velocitat d’Alfvén vA(x) (equació (2.4)) i la frequència d’Alfvén wA(x) també varien ambx. Això presenta una diferència important amb el cas de densitat discontínua on la freqüència d’Alfvén era una funció a trossos constant que botava abruptament del valor intern al valor extern. Ara la freqüència d’Alfvén és una funció contínua que connecta les freqüències d’Alfvén interiors i exteriors (veure figura2.4).
Una solució interessant de la relació de dispersió (equació (2.33)) en el cas de densitat discontinua es pot obtenir en el límit de propagació quasi-perpendicular,ky kz, que implica k2y(k2z−w2/v2A). En aquesta situació, la part de l’esquerra de la relació de dispersió es pot aproximar com
tanh(mia)'1. (2.35)
Definim wk com la freqüència en aquest límit (anomenada també com a freqüència del mode kink):
wk =kz
s
ρiv2Ai+ρevAe2 ρi+ρe
. (2.36)
El mode kink és l’únic que permet desplaçar la llosa lateralment i per tant és l’única aproximació consis- tent amb les observacions de oscil·lacions transversals de bucles coronals. Tornam al cas d’un perfil de densitat continu i anem a suposar de moment que la llosa oscil·la transversalment amb una freqüència aproximada a la freqüència kink, wk. Aquest valor sempre ens caurà dins el continu de la freqüència d’Alfvén, tal i com es pot veure a la figura 2.4. D’aquesta forma, tenim dos punts dins les capes de transició onwk =wA. Al voltant d’aquests punts (ja quewk és una aproximació) és on l’equació (2.25) presenta una singularitat i es produeix la ressonància d’Alfvén. La presència d’aquesta ressonància causa l’esmorteïment de l’ona MHD kink, ja que l’energia d’aquest mode és transferida als modes d’Alfvén al voltant del punt de ressonància.
Figura 2.4: Perfil de freqüències on podem observar com tenim dos punts on la freqüència d’Alfvén i la freqüència del mode kink s’igualen.
Tot i que no podem trobar solucions analítiques per al cas del perfil de densitat continu, podem resoldre el sistema d’equacions de forma numèrica i trobar les funcions pròpies per aquest cas. Cal comentar que tant les funcions pròpies com el valor propi (freqüència) són complexes. A més a més, la part imaginària de la freqüència és positiva, fet que ens indica que tendrem esmorteïment. La freqüència corresponents al mode fonamental ésRe(w) = 0.1027(prou similar a la del perfil de densitat continu) iIm(w) = 4.285×10−3. A la figura2.5hem representat les parts real i imaginària d’aquestes funcions pròpies. Com podem veure les parts reals són bastant similar al a les de la figura 2.3, llevat que a la funció pròpia iˆvy tenim una contribució prou important prop de les fronteres de la llosa.
Figura 2.5: Cas amb perfil de densitat continu. Representació de la part real (gràfiques superiors) i imaginària (gràfiques inferiors) de les funcions pròpies vˆx, iˆvy i iˆbz del mode fonamental lateralment evanescent. Les línies discontínues representen els límits de la llosa, la regió a la qual la densitat passa del valor de la llosa al valor de l’entorn coronal. Les funcions pròpies d’aquesta figura han estat calculades amb el programa emprat per Arregui et al. (2007) [2].
Ara que ja sabem quina forma tenen les funcions pròpies per als dos casos de densitat i hem intro- duït els conceptes bàsics, la idea és comprovar mitjançant simulacions numèriques que amb el perfil de densitat discontinu tenim oscil·lacions periòdiques, mentre que el fet d’introduir el perfil de densitat continu fa aparèixer el fenomen d’absorció ressonant i les oscil·lacions s’esmorteeixen.
3 | Simulacions numèriques d’evolució temporal
En aquest capítol presentam els resultats de dues simulacions numèriques de l’evolució temporal de la llosa quan es pertorbada amb un desplaçament inicial transversal. La secció 3.1 3.1correspon al cas de densitat no contínua a les fronteres x=±ai conté els resultats de Rial et al. (2019) [1]. La secció 3.2 correspon al cas de densitat que varia suaument entre la corona i la llosa i conté resultats nous.
A continuació resolem numèricament les equacions (2.8) – (2.14) a la regió−Lx≤x≤Lx,−Lz≤z≤Lz. La pertorbació inicial és el producte de dues gaussianes en les direccionsxizque fa la funció de depositar una energia inicial al centre de la llosa pertorbant la component transversal de la velocitat:
vx(x, z, t= 0) =v0exp
−x2 a2
exp
−z2 a2
. (3.1)
La resta de variables són inicialment zero.
Els valors que s’han utilitzat per al domini espacial són Lx = 20a i Lz = 50a i una xarxa uniforme de4001×51punts en les direccionsxiz. El fet de discretitzar la direccióxamb 4001 punts és degut a que necessitam suficient precisió per captar bé les petites escales espacials que, com mostra la figura2.3, ens poden aparèixer aprop dels límits de la llosa (x=±a). La simulació s’atura després det= 280τAi, que equival a ∼4.5 i∼75 períodes del mode fonamental lateralment evanescent i del mode fonamental lateralment confinat, respectivament. El pas temporal és∆t= 0.704τAi.
El mètode emprat per resoldre les equacions lineals és el mètode de línies, que permet tractar inde- pendentment la part temporal de la part espacial. Per la part temporal s’empra un Runge-Kutta de quart ordre, mentre que per la part espacial s’empra un mètode de diferències finites de quart ordre.
S’inclou dissipació artificial per tal d’evitar les oscil·lacions d’escala de la xarxa. El codi, anomenat MoLMHD, es descriu a Bona et al. (2009) [10].
3.1 Perfil de densitat discontinu
A les animacions de la simulació numèrica es pot veure com la condició inicial gaussiana excita un gran nombre d’harmònics tant evanescents com confinats en la direcció x, juntament amb ones que viatgen cap a l’exterior de la llosa. A partir de t≈50τAi aquestes ones que viatgen cap a l’exterior abandonen el sistema a les fronteres x=±Lxi queden únicament oscil·lant els modes evanescents i confinats. A la figura 3.1es mostren les variablesˆvx,iˆvy iˆbz per a quatre instants de temps representatius. At= 0τAi tenim la gaussiana per a la variableˆvxmentre que les altres dues variables són zero. Al cap de poc temps, quan t= 4.9τAi, podem veure com aquesta gaussiana fa ha excitat tant els modes normals lateralment evanescents com els modes normals lateralment confinats, cada un d’ells amb la seva freqüència pròpia, donant a lloc una forma on s’aprecia una superposició de tots ells. A partir det= 31.0τAija es comença a apreciar una forma de les variables semblant a la de les funcions pròpies, mostrada a la figura2.3. En aquest instant l’amplitud d’oscil·lació per a la funció pròpia vˆx és màxima. Ens trobem aproximada- ment a la meitat del període d’oscil·lació per al mode fonamental lateralment evanescent. Finalment, per at= 64.8τAi tenim completat aproximadament el primer període d’oscil·lació corresponent al mode fonamental lateralment evanescent. A partir d’aquest instant les variables continuen oscil·lant com una superposició dels modes lateralment evanescents i modes lateralment confinats.
Figura 3.1: Cas amb perfil de densitat discontinu. Variables vx, ivy ibz per a quatre instants de temps corresponents a la simulació numèrica. Una animació d’aquestes imatges es pot descarregar a https:
//drive.google.com/file/d/1YEZi54vwTecPb7o5cf0ViVyitO4C-BSU/view?usp=sharing.
L’animació dels resultats de la simulació numèrica mostra que per t .50τAi el comportament del sis- tema està dominat per un transitori, com hem mencionat més amunt. A partir d’aquest instant es veu clarament la combinació del mode fonamental lateralment evanescent, amb un període llarg, i del mode fonamental lateralment confinat, amb un període molt més curt.
Per tal d’analitzar els modes normals excitats en aquesta simulació és convenient construir un espec- tre de freqüències. A la figura 3.2 podem veure els periodogrames de Lomb-Scargle fets a partir de la simulació numèrica per a les tres variables. Per a vˆx s’ha triat el punt x = 0, z = 0, mentre que per iˆvy i ˆbz s’ha triat el punt x = a, z = 0 degut a que són els punts de major amplitud. Les línies discontínues vermelles corresponen a les freqüències dels modes lateralment evanescents trobades analíticament a la secció 2.3 mentre que la línia verda correspon a la freqüència del mode fonamental lateralment confinat trobat també a la secció2.3. Com es pot veure tenim una sèrie de pics que s’ajusten als resultats analítics. Aquests pics corresponen a les freqüències ν = 0.01601/τAi (w = 0.1006/τAi), ν = 0.04721/τAi (w = 0.2966/τAi), i ν = 0.07681/τAi (w = 0.4826/τAi) per als modes lateralment evanescents; i ν = 0.2544/τAi (w = 1.598/τAi) per al mode fonamental lateralment confinat. Aquests valors són molt semblants als resultats analítics.
Figura 3.2: Cas amb perfil de densitat discontinu. Periodogrames de Lomb-Scargle per a vˆx, iˆvy iˆbz a partir de la simulació numèrica. Les dades per un tempst <50τAi no s’han tengut en compte per fer els periodogrames. Les línies discontínues vermelles i verda corresponen a les freqüències dels primers modes normals evanescents i el mode fonamental confinat, respectivament.
3.2 Perfil de densitat continu
Repetim el procés per al perfil de densitat continu (equació (2.2)). Com podem veure a la figura3.3, els resultats obtinguts de moment són bastant semblants als del perfil de densitat discontinu (figura 3.1).
Cal destacar, però, principalment, dues diferències. En primer lloc, si comparam les amplituds de la component transversal de la velocitat en la direcció x,vˆx, per al tempst= 64.8τAi veim com per al cas del perfil de densitat continu el valor màxim (corresponent als punts x= ±a) és lleugerament inferior a 0.4, mentre que per al perfil de densitat discontinu el valor és lleugerament superior a 0.5. Veim, per tant, que l’amplitud comença a esmorteir-se amb el temps. Això és conseqüència de la segona diferència que es pot veure a simple vista: a les gràfiques per a la component en direcció y de la velocitat per als temps t = 31.0τAi it = 64.8τAi de la figura 3.3ens apareixen uns pics aproximadament a x=±a, mentre que a la figura 3.1 no ens apareixen aquests pics. La presència d’aquests pics és una evidència clara de la transferència d’energia de la component transversal de la velocitat a la componenty. Aquesta és, com vàrem dir al capítol 2, una característica típica de l’absorció ressonant a una llosa. L’animació dels resultats de la figura3.3 reforça aquestes dues conclusions, és a dir, l’esmorteïment temporal de les oscil·lacions i l’efecte de l’absorció ressonant.
A la figura3.4veim el periodograma de Lomb-Scargle per aquest cas. Veim com també s’exciten tant els modes normals lateralment evanescents com els modes normals lateralment confinats, però en aquest segon cas els pics de potència són més dispersos al voltant dels valors propis. Per altra part, el periodograma corresponent a la funció pròpiaˆbz presenta una menor excitació dels modes lateralment evanescents que en el cas de densitat discontinua. Les freqüències presents a aquesta segona iteració sónν = 0.01601/τAi (w= 0.1006/τAi), ν= 0.04747/τAi (w= 0.2982/τAi) iν = 0.07761/τAi (w= 0.4876/τAi) per als modes lateralment evanescents; iν= 0.2533/τAi(w= 1.5918/τAi) per al mode fonamental lateralment confinat.
Aquests valors són pràcticament idèntics als de la figura3.2, el que es pot interpretar com una variació de la freqüència menyspreable quan la densitat passa de tenir un bot a les fronteres de la llosa a ser contínua. Aquesta és precisament la suposició que vàrem fer a la discussió de l’absorció ressonant a la secció 2.3.2.
De moment tenim les dades de la simulació numèrica temporal en un domini espacial x−z per a les funcions pròpies vˆx, iˆvy iˆbz. Aquestes dades contenen la informació de l’evolució temporal de tots els modes normals tant lateralment evanescents com confinats que s’han excitat a partir de la pertorbació inicial gaussiana. El que interessa a partir d’ara és obtenir la distribució espacial en direcció x de les variables per a cada un dels diferents modes per separat. El problema és que a les figures 3.1 i 3.3 les funcions pròpies que veim són una superposició dels diferents modes que s’han excitat amb la pertorbació inicial. A la següent secció explicam en què consisteix un anàlisi CEOF. Analitzant les dades de la simu- lació numèrica aconseguirem extreure tres funcions pròpies (pervˆx,iˆvy iˆbz) que s’aproximaran a les del mode que nosaltres seleccionem; en el nostre cas estudiarem el mode fonamental lateralment evanescent.
Ja que coneixem la distribució espacial de les funcions pròpies (figures 2.3 i 2.5) podrem comparar els resultats extrets del CEOF amb les funcions pròpies.
Figura 3.3: Igual que la figura 3.1 però per al cas del perfil de densitat continu. Una animació d’aquestes imatges es pot descarregar a https://drive.google.com/file/d/
1xuaSdv8dVuFqW5r38X7Rid1pIhmSdhFg/view?usp=sharing.
Figura 3.4: Igual que la figura 3.2però per al cas del perfil de densitat continu.
