Numerical simulation of Alfven waves in non uniform plasmas

University of Seville Degree in Physics Numerical simulation of Alfv´ en waves in non uniform plasmas Daniil Kabirov Leontieva Supervisors: Eleonora Viezzer and Haotian Chen Departamento de F´ ısica At´ omica, Molecular y Nuclear Facultad de F´ ısica Plasma Science and Fusion Technology Abstract Nuclear fusion can provide a sustainable and nearly unlimited source of energy, but for it to become a viable alternative, the fuel must reach very high temperatures, which inevitably leads to the creation of a plasma medium inside a fusion reactor. Thus, a thorough understanding of this state of matter is essential for this goal. One of the best known models of the plasma is the magneto-hydrodinamic model, or MHD, which describes it as particular fluid, subjected to equations analogous to Naver-Stokes’ equations, where the Maxwell equations also play an important role in defining the behaviour of the plasma. Using this model, it is possible to predict certain instabilities that hinder the containment of the plasma in a nuclear fusion reactor. One of the most relevant instabilities of this kind is the Shear Alfv´ en Wave, or SAW, the object of study of this thesis. For this project a simple plasma system is assumed, from which a compact equation for the Shear Alfv´ en Wave is deduced. However, the equation features a singular point, which is treated performing a Fourier Transform and solving it through numerical methods, returning a profile for the SAW in real space. Contents 1 Introduction 2 2 Theoretical background 5 2.1 Introduction to plasma physics and MHD model . . . . . . . . . . . . . . . . . 5 2.2 SAW in homogeneous plasma . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Linear equations for a non-homogeneous plasma . . . . . . . . . . . . . . . . 15 2.4 Alfv´ en-Interchange-Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.5 Simplifications of the deduced equation . . . . . . . . . . . . . . . . . . . . . 24 3 Numerical implementation 31 3.1 Obtaining a low order approximation . . . . . . . . . . . . . . . . . . . . . . . 31 3.2 Boundaryconditions ............................... 33 3.3 Solvercode .................................... 36 3.3.1 Summary of the problem . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.3.2 Differential equation solver . . . . . . . . . . . . . . . . . . . . . . . 38 3.3.3 Solving the eigenvalue problem . . . . . . . . . . . . . . . . . . . . . 39 4 Results 42 5 Summary and discussion 46 6 Appendix 48 6.1 Vector algebra and calculus identities . . . . . . . . . . . . . . . . . . . . . . 48 Bibliography 49 1 Chapter 1 Introduction Our modern society, despite all its advancements in our general quality of life, is faced with a grave challenge: a constantly increasing energy demand for an ever growing human population, which is currently supplied, in large part, by fossil fuels. This energy source is not only threatening the long term survival of Earth’s ecosphere, but can also lead to undesirable geopolitical tensions due to its in inhomogeneous distribution around the globe. Thus, humanity is in dire need of alternative energy sources, one of them being the main motivator behind this thesis: nuclear fusion. The nuclei that make up all the elements around us, as we know, are composed of neutrons and protons. However, it is well known that the mass of a certain nucleus isn’t exactly the sum of the masses its components: there is a certain mass difference associated in the whole nucleus which is due to the strong interactions between the nucleons. Using Einstein’s mass energy-equivalence E=mc2, this mass excess can be associated to a binding energy of a system, which is the sum of the rest energy all its components minus its total rest energy BE=Pmic2−Mc2. Thus, by managing to transmute a system of protons and neutrons bound in a certain configuration, into another configuration of the same protons and neutrons which is more bound (less massive), then there is an energy excess which translates into kinetic energy of the final system that can be eventually repurposed for energy generation. 2 3 As it is shown in the following graph (Figure 1.1), one of those energetically viable processes is to combine two light nuclei into a heavier nucleus which is more bound than its original components. Figure 1.1: Chart of binding energies per nucleon for stable nuclei. From: https://www.researchgate.net/figure/Binding-energy-per-nucleon-B-Z-N-A-as-a-function-of-the-mass-number fig1 234065870 However, to make that possible the two nuclei must overcome the Coulomb barrier due to the electrostatic interaction of the protons, which requires a certain amount of kinetic energy. In the Sun’s core, where the pressure is in the order of hundreds of billions of atmospheres, this process starts at proton level, eventually producing 4 2He. Meanwhile, on Earth it is impossible to reach such pressures, so a different reaction is used to achieve fusion: 2 1D + 3 1T−−→ 4 2He + 1 0n + 17.6 MeV Note that other reactions are also possible in a nuclear fusion reactor, but are considerably less frequent. In this reaction, 2 1Dand 3 1Tdenote deuterium and tritium, respectively, which can be easily and cheaply obtained from Earth’s oceans: directly, in the case of deuterium, or indirectly, by obtaining lithium and transmuting it into tritium in the same nuclear fusion reactor. However, in order to achieve this nuclear reaction, a temperature in the order of 100 millions of degrees Celsius is needed, which will lead to the complete ionization of the atoms in the reactor: 4CHAPTER 1. INTRODUCTION thus a plasma is created. In addition, due to the extremely high temperatures, the particles have very large velocities (of the order of 107m/s), easily escaping containment. One of the proposed solutions to this problem are Tokamak fusion reactors, consisting of a toroidal field coil that generates a a magnetic field perpendicular to the solenoid surfaces. Due to the fact that the plasma is composed of charged particles, they experience a Lorentz force which makes them follow the magnetic field lines in a helicoidal orbit. However, due to the weakening of the magnetic field the further we go from the rotational axis of symmetry, there is a drift of the particle towards the outer walls of the solenoid. Thus, a secondary magnetic field in the poloidal direction is created, twisting the total magnetic field in a way that the particles are stuck in closed orbits due to the conservation of energy and magnetic momentum. A schematic representation of a Tokamak is shown in the next figure: Figure 1.2: A Tokamak with its magnetic fields. From: https://www.energy.gov/science/doe-explainsTokamaks There is an additional challenge for containment: the plasma is subject to instabilities in a macroscopic scale, which are hard to control phenomena that eventually lead to the loss of containment of the plasma. These disruptions can be characterized using the MHD (Magnetohydrodinamic) model. As it name implies, it assumes that the plasma can be modeled as a fluid governed Maxwell equations. As it will be seen soon, one of the most well studied instabilities is the Shear Alfv´ en Wave, which will be the main topic of this thesis. Chapter 2 Theoretical background 2.1 Introduction to plasma physics and MHD model A plasma is defined formally as an ionized gas which is is quasi-neutral and displays collective effects. In more detail: • Quasi-neutrality: despite the fact that the plasma is composed of positive (ions) and negative particles (electrons), for scales much larger than to a coefficient called Debye Length, λD, the positive charge and negative charge densities can be considered the same in modulus, thus cancelling out. This coefficient has the following expression: λD=sϵ0kbTe e2n0 (2.1) where kbis the Boltzmann constant, n0the mean density of electrons and Tethe temperature associated to electrons, which is generally different from the temperature associated to the ions. Note that these temperatures aren’t defined as in Thermodynamics, as that would require a system in equilibrium. 5 12 CHAPTER 2. THEORETICAL BACKGROUND For this system to be solvable, its determinant must be zero. Thus, 3 different roots for ω2are obtained: ω2=k2 ∥v2 A(2.22a) ω2=1 2k2(v2 A+v2 S)(1 ±√1−α2); α= 4k2 ∥ k2 v2 Av2 S (v2 A+v2 S)2(2.22b) The first root (2.22a) corresponds to the Shear Alfv´ en Wave (SAW). Substituting it into the system of equations (2.21), v1xbecomes undetermined, while the remaining two equations are trivial, giving us v1y=v1x= 0. Therefore, v1xis non-zero. It can be easily seen that the Shear Alfv´ en Wave is incompressible, as ∇ · v1= 0. It can be proven that when an arbitrary plasma displacement compresses the fluid, it must perform work, decreasing its kinetic energy and increasing the potential energy, stabilizing the perturbation. As this effect is not present in SAW waves, they tend to be more unstable, making their study to be of special interest. Recalling the definition of the group velocity of a 3D wave: vg=∇kω, the group velocity of the SAW is the following: vG=vA B0 B0 (2.23) Thus, the collective motion of the perturbation is parallel to the magnetic field. However, from Ohm’s Law (2.29d) and Faraday’s law (2.29e), it is deduced that B1is parallel to the perturbed velocity V1and perpendicular to the equilibrium magnetic field B0: v1=±vA B1 B0 (2.24) Recalling that the perturbed quantities have an oscillatory element to them f1(r, t) = ˆ fe−i(ωt−k·r), the modulus of B1oscillates along z. That means that there is a field line bending, that leads to a perturbed current J1due to Ampere’s Law (2.29f), which in turn leads to a time variation in v1 due to the momentum equation (2.29b). Thus, the SAW is transverse wave with an oscillation mechanism analogous to a string vibrating under tension. A figure representing the bending can be seen in the following page. 2.2. SAW IN HOMOGENEOUS PLASMA 13 Figure 2.1: Diagram of the Shear Alfv´ en Wave in a homogeneus plasma. From [7]. The second mode, which is the positive root in (2.22b), is called the fast magneto-acoustic wave, as it can easily be seen that its characteristic frequency, given by the equation in (2.22b), is greater than the Alfv´ en frequency. There is an interesting limit which is relevant to a plasma in a Tokamak. It is given by the βparameter, defined as the ration between the plasma scalar pressure Pwe have defined in the previous section, and the magnetic pressure Pmag =B2 0 2µ0 associated to the magnetic field, which behaves similarly. β=P Pmag =2µP B2 0 (2.25) Usually Tokamaks function in the low βlimit, that is, the pressure associated to the magnetic field is much bigger than the fluid pressure, which means that that v2 s/v2 a≪1, Then, the fast magneto-acoustic wave, which is then called fast compressional wave, has the following dispersion relation: ω2=k2v2 A(2.26) The group velocity would be: vG=vA k k(2.27) 14 CHAPTER 2. THEORETICAL BACKGROUND The perturbed fluid velocity v1components are now v1x= 0 and v1y, v1z= 0, thus ∇·v1= 0: in this mode of oscillation the plasma is compressible, and has a longitudinal and transversal component. This leads to this mode being described by a magnetic tension force, as in the last case, plus a magnetic compression force, as in a sound wave. For the negative root in (2.22b), the slow magneto-sonic wave appears, whose frequency is smaller that the SAW frequency. In the low beta limit of interest: ω2=k2 ∥v2 S(2.28) In this mode the results for the perturbed velocity are the same, but in the low beta limit case: v1y/v1z≫1, which means that the main force at play here is the plasma compression. As a final note, let’s make another analogy with the elastic string under tension for SAW and fast compressional waves: recalling the relationship between the elastic constant of a spring and its frequency, it turns out that the terms k∥vAfor SAW and k vAfor the fast compressional wave act as the analogue for said elastic constants in the plasma. Usually, the bigger the elastic constant, the more kinetic energy is needed to produce a perturbation of the same amplitude. This is another argument by which compressional waves are more stable than SAW. Energetically speaking, a k∥vAelastic constant implies kinetic energy absorption through the bending of the magnetic field lines, while a k⊥vAelastic constant corresponds to compressional energy absorption. For SAW waves, the modes of interest are the more dangerous and unstable ones that can get the plasma out of confinement, so when studying SAW, |k∥|is made as small as possible, of the order of the inverse of the typical length of the Tokamak. Moreover, there is another energy absorbing mechanism we want to minimize, which happens for large |k∥|: damping through wave-particle resonance, where the wave energy is transferred to the resonant particles, also called Landau Damping. Essentially, what happens is that particles moving with a similar speed an direction as an MHD wave, get accelerated if they are moving slightly slower than the wave and decelerated in the opposite case. Thus, if the number of particles that go slower is greater than the ones that go faster, the wave cedes energy to the slow particles, losing energy. 2.3. LINEAR EQUATIONS FOR A NON-HOMOGENEOUS PLASMA 15 This mechanism can’t be explained by the MHD fluid model, these phenomena won’t play a role in this case. Still, it is another justification for the approximation |k∥| ≪ |k⊥|that can be often found in the bibliography, as in references [3] and [2] 2.3 Linear equations for a non-homogeneous plasma The main topic of this thesis is the characterization of the SAW in a non-homogeneous plasma slab, where the equilibrium quantities may depend on a single coordinate. Before performing the linearization, some simplifications can be made: due to the quasi-neutraility condition, E0= 0and it will be assumed that ∇·v1= 0, which is equivalent to saying that the perturbed mode will be in-compressible. There are two reasons behind this simplification: • There are multiple types of perturbations, or plasma modes, that can be obtained in these equations. While the Shear Alfv´ en Mode is the one which is in-compressible, it isn’t the case for other modes. Thus, it is a way to isolate the study of SAW. Still, fast compressional modes will still play a part in explaining the the behaviour of the wave later. • The fact that a perturbation is compressible has been proven to make the system more stable. The phenomena of interest are the ones that could destabilize the system, as they are of relevance for plasma containment in Tokamaks. Thus, the linearized MHD equations are the following: ∂ρ1 ∂t +∇ρo·v1= 0 (2.29a) ρo ∂v1 ∂t =−∇P1+1 cJ1×B0+1 cJ0×B1(2.29b) ∂P1 ∂t +v1·∇P0= 0 (2.29c) E1+1 cv1×B0= 0 (2.29d) ∇×E1=−1 c ∂B1 ∂t (2.29e) 16 CHAPTER 2. THEORETICAL BACKGROUND ∇×B1=4π cJ1(2.29f) ∇·B1= 0 (2.30) 2.4 Alfv´ en-Interchange-Equation Having obtained the linearized MHD equations, the next objective is to obtaining a single equation that can describe the SAW. Thus, all the perturbed quantities will be expressed in terms of the perturbed potential ϕ1, which makes it easier to add the aformentioned wave-particle interaction phenomenon for a future project. The first step is multiplying vectorially the equation (2.29b) B0from the right side: B0×ρo ∂v1 ∂t =−B0×∇P1+1 cB0×(J1×B0) + 1 cB0×(J0×B1)(2.31) However, B0doesn’t depend explicitly on t, so it can be inserted into the time derivative. Invoking (2.29d) along with the vector property A1, from Appendix A, leads to: ρoc∂E1 ∂t =−B0×∇P1+1 cB0×(J1×B0) + 1 cB0×(J0×B1)(2.32) The second right-hand side term (RHT) can be transformed by applying (2.29f) and vector property A2: B0×(J1×B0) = cB2 0 4π[(∇×B1)−(∇×B1)·bb]≡cB2 0 4π(∇×B1)⊥(2.33) where b≡B0 |B0|is the unitary vector parallel to the magnetic field line, bb is the tensor diad that verifies v·bb ≡(v·b)b,vbeing a vector in 3D space. With this notation, v−v·bb ≡v⊥ is defined as the projection of the vector vperpendicular to the magnetic field line. Substituting the new expression (2.33) in (2.32): ρoc∂E1 ∂t =−B0×∇P1+B2 0 4π(∇×B1)⊥+1 cB0×(J0×B1)(2.34) 2.4. ALFV ´ EN-INTERCHANGE-EQUATION 17 Differentiating the equation by ∂ ∂t and using the relations given by (2.29c) and (2.29f) leads to the following expression: ρoc∂2E1 ∂t2=B0×∇(∇P0·v1)−cB2 0 4π(∇×∇×E1)⊥+B0×[J0×(∇×E1)] (2.35) Assuming that v1⊥B0, which is true for SAW waves in an homogeneous plasma, Ohm’s Law (2.29d) can be inverted: v1=−c B2 0 E1×B0(2.36) This result can be used in the first RHS term in (2.35): ρo ∂2E1 ∂t2=B0×∇"(b×∇P0)·E1 B0#−B2 0 4π(∇×∇×E1)⊥+1 cB0×[J0×(∇×E1)] (2.37) where the triple product has been reorganized and the cconstant has been cleared. Now it is necessary to review the electromagnetic potential formulation, by which the perturbed electromagnetic fields can be expressed using the scalar and vector perturbed potentials: E1=−∇ϕ1−1 c ∂A1 ∂t (2.38a) B1=∇×A1(2.38b) Coulomb’s Gauge has also been used: ∇·A1= 0. Taking into account that the equation is not only linear for any of the three Cartesian coordinates, but also that the equilibrium properties of the plasma will vary along one of those coordinates, which will be x. Therefore, as in the last section, every perturbed quantity can be expressed as f1=˜ f1(x)e−i(kyy+kzz). Moreover, the del operator into two parts: one parallel and one perpendicular to the magnetic field lines: ∇=∇∥+∇⊥. Using Coulomb’s Gauge along with said decomposition: ∇·A1=∇∥·A1+∇⊥·A1=k∥v1∥+k⊥v1⊥ ! = 0 (2.39) 18 CHAPTER 2. THEORETICAL BACKGROUND As it has been assumed that |k∥|≪|k⊥|, the following relation holds true: |A1∥|≫|A1⊥|=⇒A1≈A1∥b(2.40) However, going back to Ohm’s law (2.29d), the perturbed electric field is perpendicular to the magnetic equilibrium field due to the properties of the vector product, so E1∥= 0. Using the electrostatic potential equation (2.38a) for said parallel component: ∇∥ϕ1=−1 c ∂A1∥ ∂t (2.41) By inserting that condition into the E1(2.38a), taking the curl of the equation with the approximation (2.40): ∇×E1=((((( ( −∇×∇ϕ1−1 c∇× ∂A1∥ ∂t b!=∇×(∇∥ϕ1b)(2.42) However, the approximation (2.40) also leads to: E1=E1⊥=∇⊥ϕ1(2.43) These results can be inserted into the previous equation (2.37), leading to: ρ∂2E1 ∂2t=B0×∇"(b×∇P0)·∇⊥ϕ1 B0# | {z } I +B2 0 4π[∇×∇×(∇∥ϕ1b)]⊥ | {z } II + +1 cB0×J0×[∇×(∇∥ϕ1b)] | {z } III (2.44) The equation is now going to acquire a lot of terms, so each of the terms has been labeled with a roman numeral as they will be treated separately. 2.4. ALFV ´ EN-INTERCHANGE-EQUATION 19 Term I: B0×∇"(b×∇P0)·∇⊥ϕ1 B0#(2.45) In order to continue, the following expression must be proven: ∇× b B0 =4π cB2 0 (J∥−J⊥) + 2 Bb×κ(2.46) where κ≡(b·∇)bis the curvature vector of the magnetic field. Proof (2.46): First, the LHT (left hand term) is expanded by virtue of the successive application of vector calculus property A4: ∇× b B0 =1 B0∇×b−1 B2 0∇B0×b=2 B0∇×b−1 B2 0∇×B0(2.47) For the first RHT in the original equation (2.46), (2.29f) and the definition of the the vectors parallel v∥and perpendicular v⊥to the equilibrium magnetic field are used, and the obtained expression is simplified by using A4 where B0=B0b: 4π cB2 0 (J∥−J⊥) = 2 B2 0 (∇×b)bb −1 B2 0∇×B0(2.48) For the second RHT in (2.46), the alternative definition of κwhich is κ=−b×(∇ × b)is applied: 2 B0 b×κ=−2 B0 b×b×(∇×b) = −2 B0 (∇×b)bb +2 B0 (∇×b)(2.49) By adding (2.48) to (2.49) we get (2.47), thus proving (2.46).End of proof. 20 CHAPTER 2. THEORETICAL BACKGROUND Using curl chain rule A4 in the term I (2.45) multiple times, while writing B0as B0=B2 0 b B0 B0×∇(b×∇P0)·∇⊥ϕ1 B0= =B2 0−∇×(b×∇P0)·∇⊥ϕ1 B2 0 b+∇×b B0(b×∇P0)·∇⊥ϕ1 B0 (2.50) This way, applying the proof at (2.46), term 1 (2.45) is decomposed into the following expansion: B0×∇(b×∇P0)·∇⊥ϕ1 B0=4π cJ∥(b×∇P0)·∇⊥ϕ1 B0 | {z } I.A −B2 0∇×(b×∇P0)·∇⊥ϕ1 B2 0 b | {z } I.B −4π c[J⊥−B0×κ](b×∇P0)·∇⊥ϕ1 B0 | {z } I.C (2.51) Term II: Ignoring the prefactor, the term II would be: [∇×∇×(∇∥ϕ1b)]⊥≡[∇×∇×(∇∥ϕ1b)] −[∇×∇×(∇∥ϕ1b)]bb (2.52) In order to perform a similar expansion of the terms, the following expression must be proven: [∇×∇×(∇∥ϕ1b)]bb ={(∇∥ϕ1)(∇·κ)−∇·[∇⊥(∇∥ϕ1)] + (∇∥ϕ1)|∇×b|2}b(2.53) Proof (2.53):For (∇∥ϕ1), which is a scalar, the notation (∇∥ϕ1)≡Fwill be used. First, use κ=−b×(∇×b)on the first RHT: F∇·κ=−F∇·[b×(∇×b)] = −F|∇×b|2+F(∇×∇×b)·b(2.54) Then, applying A4 to the second RHT of the immediately previous equation (2.54) two times: F(∇×∇×b)·b={∇×[∇×(Fb)] −∇×(∇F×b)−∇F×(∇×b)}·b(2.55) 2.4. ALFV ´ EN-INTERCHANGE-EQUATION 21 Now, by using A5 on the second RHT of the original equation (2.53): ∇·[∇⊥F] = −∇·{b×(b×∇F)}= = [∇×(b×∇F)] ·b−(b×∇F)·(∇×b) (2.56) Now (2.54), (2.55), (2.56) have been substituted in (2.53). Some terms will go away, resulting in: F∇·κ−∇·(∇⊥F) + F|∇×b|2= = (b×∇F)·(∇×b) + {∇×[∇×(Fb)] −∇F×(∇×b)}·b (2.57) By using the cyclic property of the triple product in the third RHT it goes away with the first term: F∇·κ−∇·(∇⊥F) + F|∇×b|2={∇×[∇×(Fb)]}·b(2.58) The equation we get is precisely (2.53), but in scalar form. Multiplying by bon both sides gives the original equation.End of proof. Now, the term II (2.52) can be expressed in the following way: B2 0 4π[∇×∇×(∇∥ϕ1b)]⊥≡B2 0 4π[∇×∇×(∇∥ϕ1b)] | {z } II.A −B2 0 4π[∇×∇×(∇∥ϕ1b)]bb; [∇×∇×(∇∥ϕ1b)]bb ={(∇∥ϕ1)(∇·κ) | {z } II.B −∇·[∇⊥(∇∥ϕ1)] | {z } II.C + (∇∥ϕ1)|∇×b|2 | {z } II.D }b (2.59) Note that the terms II.B, II.C, II.D also include the band the B2 0 4πprefactors. Term III: The term in question is: B0×J0×[∇×(∇∥ϕ1b)] (2.60) 28 CHAPTER 2. THEORETICAL BACKGROUND In the kink term, we have the perpendicular component a vector, which points perpendicularly to the magnetic field, thus the divergence acts as ∇⊥: ∇·"4πJ0∥ cB0∇×(b∇∥ϕ1)#⊥∼k⊥k2 ∥ ϵ a∼ϵ3 a4(2.89) The final term has a mix of both parallel and perpendicular components of J0, which verify J0∥J0⊥. Therefore: ∇·"16π2 c2B2 0 (J∥0−J0⊥)(J∥0−J0⊥)·∇ϕ1#∼ϵ4 a4(2.90) After evaluating the order of each term, leaving only the lowest order terms leads to: ∇·"1 v2 a∇⊥ ∂2 ∂t2ϕ1#−∇·hb∇·∇⊥(∇∥ϕ1)i= 0 (2.91) It is possible to perform a discrete Fourier expansion for ϕ1in the yand zcoordinates, for the only dependence of the involved variables is on the x axis, as well as a forward Laplace time transform: ϕ1= ∞ X ky,kz=∞ ˆ ϕ(x, ky, kz, p)ei(kyy+kzz);ky, kz∈R(2.92) ˆ ϕ≡ L(ϕ1) = Z∞ 0dt ϕ1(x, t)e−pt;p∈C(2.93) Due to the fact that there is no dependence in ynor z, there is a separate equation for each ˆ ϕ term, in which the complex exponential is also omitted. Also, the following properties of the forward Laplace transform will be used: L ∂ϕ ∂t !=pˆ ϕ−ϕ0;ϕ0=ϕ1(r, t = 0) (2.94a) L ∂2ϕ ∂t2!=p2ˆ ϕ−pϕ0−˙ ϕ0;˙ ϕ0≡∂ϕ1 ∂t r,t=0 (2.94b) 2.5. SIMPLIFICATIONS OF THE DEDUCED EQUATION 29 Substituting these identities into (2.91): ∇·"1 v2 a∇⊥(p2ˆ ϕ−pϕ0−˙ ϕ0)#−∇·hb∇·∇⊥(∇∥ˆ ϕ)i= 0 (2.95) Now define an equilibrium magnetic field which verifies all the previous assumptions: B0=B0 ez+α(x)ey q1 + α2(x)(2.96) in which α(x)≪1. Thus: ∇∥f≡ikz+kyα √1 + α2 ez+αey √1 + α2f≡ik∥fb(2.97) ∇⊥f≡ ∇f−∇∥f=∂f ∂xex−if (ky−αkz)(ey−αez) 1 + α2(2.98) Applying these operators to the first term of (2.95): ∇·"1 v2 a∇⊥(p2ˆ ϕ−pϕ0−˙ ϕ0)#=1 v2 a ∂2 ∂x2−ˆ k!(p2ˆ ϕ−pϕ0−˙ ϕ0)(2.99) ˆ k≡(ky−αkz)2 1 + α2 In order to simplify the equation, the initial condition will be such that it verifies: ∂2 ∂x2−ˆ k!(pϕ0+˙ ϕ0) = 0 (2.100) Now, the second term of (2.95) would be: ∇·hb∇·∇⊥(∇∥ˆ ϕ)i=−k∥(∂2 ∂x2−ˆ k)(k∥ˆ ϕ) = ˆ kk2 ∥ˆ ϕ−k∥ ∂2k∥ ∂x2ˆ ϕ+ 2∂k∥ ∂x ∂ˆ ϕ ∂x +k∥ ∂2ˆ ϕ ∂x2  (2.101) 30 CHAPTER 2. THEORETICAL BACKGROUND Putting (2.99) and (2.101) together and grouping the derivatives of the same order: − p2 v2 A +k2 ∥!ˆ kˆ ϕ+k∥ ∂2k∥ ∂x2ˆ ϕ+ p2 v2 A +k2 ∥!∂2ˆ ϕ ∂x2+ 2k∥ ∂k∥ ∂x ∂ˆ ϕ ∂x = 0 (2.102) Assuming that α(x)is small as well as all its derivatives (very smooth change in direction), so that the second order derivative of k∥goes away. In addition, the last two terms can be joined together: ∂ ∂x   p2 v2 A +k2 ∥!∂ˆ ϕ ∂x − p2 v2 A +k2 ∥!ˆ kˆ ϕ= 0 (2.103) The derivative of vA(x)2can also be considered negligible compared to the derivative of ˆ ϕThus, multiplying by v2 Aand simplifiying some terms leads to a very simple equation: ∂xP ∂xˆ ϕ−Qˆ ϕ= 0 (2.104) P≡p2+k2 ∥v2 A;Q≡ˆ k P This equation has a singular point at p2=−k2 ∥v2 Awhere P= 0. In the next section, a Fourier transform for the remaining spatial variable xwill be performed to get rid of it. Chapter 3 Numerical implementation 3.1 Obtaining a low order approximation While the previous equation has a simple expression, if all the functions are written explicitly in terms of x, the whole equation can get quite complex. Thus, more simplifications must be made: •k0is a function of x due to the α(x)term. We assume for now that α(x)=0so k0(x) = k2 y, simplifying the calculations. • The function is even around x= 0. This consideration doesn’t impact the physics of the problem and it creates additional restriction which allows obtaining useful eigenvalues out of our eigenvalue equation. • Then the Alfv´ en frequency is defined as k2 ∥v2 A≡ω2 A(x), which is assumed analytical and can be expanded in a Taylor series up to second order around x= 0. Performing the aforementioned expansion: ω2 A(x) = ω2 A(0) + ω2 A ′(0)x+1 2ω2 A ′′(0)x2(3.1) 31 32 CHAPTER 3. NUMERICAL IMPLEMENTATION In order to simplify the notation, define λ=ω2 A(0) and λnas it’s successive derivatives: ω2 A(x) = ˜ λ+λ1x+1 2λ2x2;˜ λ=p2+λ(3.2) Substituting in the previous equation: (λ1+λ2x)∂xˆ ϕ+ (˜ λ+λ1x+1 2λ2x2)∂2 xˆ ϕ−k2 y(˜ λ+λ1x+1 2λ2x2)ˆ ϕ= 0 (3.3) Now, assuming that ˆ ϕis well behaved, a Fourier transform can be performed: Φ(k)≡Z∞ −∞ dx ˆ ϕ(x)e−ikx;k∈R(3.4) kmust be real by definition of the Fourier transform, later on, it will have to be to extended analytically to the complex plane due to convergence issues. Let f(x)be a well behaved function, fn(x)its n-th derivative, and F[ ] the Fourier transform operator. Some useful properties of the Fourier transform are: F[fn(x)] = (ik)nF[f(x)] (3.5a) F[xnf(x)] = (i∂k)nF[f(x)] (3.5b) F[xnfm(x)] = (i∂k)m{(ik)nF[f(x)]}(3.5c) These properties can be deduced by integration by parts assuming f(x)vanishes at infinity. Now these properties can be substituted in the equation, in which Φ≡F(ˆ ϕ)and Φn= (∂k)nF(ˆ ϕ): ikλ1Φ−λ2(kΦ′+ Φ) −k2˜ λΦ−iλ1(k2Φ′+ 2kΦ) + 1 2λ2(k2Φ′′ + 4kΦ′+ 2Φ)−(3.6) −k2 y(−1 2λ2Φ′′ +iλ1Φ′+˜ λΦ) = 0 3.2. BOUNDARY CONDITIONS 33 Regrouping all the derivatives of Φntogether: −[(k2+k2 y)˜ λ+ikλ1]Φ + [kλ2−i(k2+k2 y)λ1]Φ′+1 2(k2+k2 y)λ2Φ′′ = 0 (3.7) Due to the assumption ω2 Ais even around 0, so is λ1= 0: −(k2+k2 y)˜ λΦ + kλ2Φ′+1 2(k2+k2 y)λ2Φ′′ = 0 (3.8) Simplifying this equation even more by taking common factor k2 yand 1 2λ2leads to: n∂θh(1 + θ2)∂θi+ Ω2(1 + θ2)oΦ = 0 (3.9) where θ≡k kyand Ω2≡ −2˜ λ λ2. Thus, a second order equation in k-space is obtained, with a arbitrary complex parameter Ω2. Now appropriate boundary conditions have to be imposed in order to solve it. 3.2 Boundary conditions By writing the equation explicitly: (1 + θ2)∂2 θΦ+2θ ∂θΦ+Ω2(1 + θ2) Φ = 0 (3.10) Taking the limit θ≫1leads to the following simplified form: θ2∂2 θΦ+2θ ∂θΦ+Ω2θ2Φ=0 (3.11) Now a new variable Ψ = θΦis defined, and after some algebra: ∂2 θΨ = −Ω2Ψ(3.12) 34 CHAPTER 3. NUMERICAL IMPLEMENTATION The two possible solutions for this trivial equations are Φ±=C e±iΩθ, with Cbeing a complex scalar. Undoing the variable change, the asymptotic boundary condition appears: Φ±=Ce±iΩθ θ(3.13) The question now is which one of the solutions to select. For that goal the causality condition, which will now be explained, is used. In order to undo the Fourier transform in xand the Laplace transform in t, the following integral is done: ˆ ϕ(x, t) = Klim T→∞ Zη+iT η−iT dΩZ∞ −∞ dθ Φ(θ, Ω)ei(θ x−Ωt)(3.14) Where Kis a constant complex constant due to the variable change from kto θand from pto Ω, while ηis a real number which must be chosen in a way that the function converges through all the complex integration path. This way the instability can be regarded as a group of plane waves propagating in a certain direction. Assuming it is possible to perform the following Taylor expansion up to first order: Ω(θ) = Ω0+dΩ dθ θ=0 (θ−θ0)(3.15) What is left is the the classical expression for a wave group modulated by a perfect monochromatic wave: ˆ ϕ(x, t) = Kei(θ0x−Ω0t)ZΩdΩZθdθ Φ(θ, Ω)ei(θ−θ0)(x−dΩ dθ t)(3.16) The exponent in the integral will oscillate periodically in the domain of θ, as the (x−dΩ dθ t) term is constant when integrating for all θ. Thus, the final result will be really small unless the exponent is approximately zero. Putting the time origin at t= 0, and due to symmetry, the perturbation can be considered to propagate along the right half of the plane. Thus: x t≈dΩ dθ θ=0 >0(3.17) 3.2. BOUNDARY CONDITIONS 35 However, instead of Ω, it is possible to expand θ: θ(Ω) = θ0+dθ dΩθ=0 (Ω −Ω0)(3.18) This way, the integral becomes: ˆ ϕ(x, t) = Kei(θ0x−Ω0t)ZΩdΩZθdθ Φ(θ, p)ei(Ω−Ω0)( dθ dΩx−t)(3.19) However, Ωis complex, so the exponential will not simply oscillate, but it can increase or decrease in modulus depending on the sign of the real exponential. Thus, the integral will have a divergent or convergent factor e−Im(Ω) dθ dΩx. Assuming Ω(θ)and θ(Ω) are continuous and differentiable functions, the Inverse Function Theorem along with (3.17), can be used: dθ dΩ= dΩ dθ !−1 >0,for x > 0(3.20) There are certain phenomena in plasma which are unstable: they start in a certain point (let’s assume x= 0 in this case), expand from that point and then dissipate when they go further in time and space from the starting point. Therefore, the plane waves that compose it must be spatially decreasing in t= 0 and Im(Ω) >0, otherwise the real exponential factor leads to a phenomenon that propagates at infinite speed through all space when it appears at a certain point: a violation of causality. However, the integral is performed for all values of θpositive and negative, so unless some condition on the plane wave amplitude Φ(θ, Φ) is put, part of the integral would diverge for large values of θ. Thus, it is imposed that: Φ(θ, Φ) ≈Cesign(θ)iΩθ θ(3.21) Nonetheless, the function must be even in Fourier space, as the original function is even and real in x space as well. Therefore, another boundary condition must be ∂θΦ(θ= 0) = 0. 36 CHAPTER 3. NUMERICAL IMPLEMENTATION 3.3 Solver code 3.3.1 Summary of the problem The problem that needs to be solved is the following: the equation (3.9) is a second order differential ordinary equation for Φ(θ)in the complex plane where the independent variable θ has been assumed to be real, spanning the whole set of R. However, due to the parameter Ω, the equation becomes part of an eigenvalue problem, where the solution verifies the following boundary conditions: 1. The solution must be even, both in the Fourier Space and Real Space. This can be accomplished by imposing: ∂θΦ(θ= 0) = 0 2. For large θ∈R, the equation must tend asymptotically to: Φas(θ)≡Cesign(θ)iΩθ θ Numerically, there are no strict infinities or strict zeros, so the boundary conditions must be modified in the following way: 1. The even parity condition can be approximated as Φ′(θ= 0) = ϵ≪1. However, a problem arises when the value isn’t exactly zero: the value Φ(θ= 0) tends to matter. For example, the numerical derivative of f(x) = Acos(x)at x≈0is f′(x≈0) = Aϵ: the amplitude of the function is affecting the value of ϵ. Those, a normalized derivative is defined: F≡∂kΦ(θ= 0) |Φ(θ= 0)|=ϵ≪1(3.22) 2. A large value of |θ| ≫ 1is used instead of infinity, in this case θ= 100. Thus, the boundary condition is transformed into a initial condition for the value of the solution 3.3. SOLVER CODE 37 and its derivative at θ= 100 : Φ(θ= 100) = Φas(θ= 100) = eiΩθ θθ= 100 (3.23a) ∂θΦ(θ= 100) = ∂θΦas(θ= 100) = (iΩθ+ 1) e−iΩθ θ2θ= 100 (3.23b) There are three three different boundary conditions for a second order ODE. That is because Ω is undetermined, so while the ODE is solved with the asymptotic initial conditions (3.23a) and (3.23b), we are going to determine Ωby using (3.22). Thus, the proposed problem becomes eigenvalue problem. Another restriction to the solution is that stable plasma modes are considered: Im[Ω] <0. This is due to the fact that in the current model, without any kinetic effects, there are no dissipation mechanisms for the Shear Alfv´ en Wave, thus it won’t dissipate for large values of |x|. This means that the imposed boundary conditions in θspace are divergent. This is problematic, as computing a solution for divergent function gives rise to numerical complications. However, by extending θto the complex plane and solving the problem along a straight line in the complex plane, passing through the origin, which due to the Stokes phenomenon, will give us a convergent solution. The Stokes phenomenon is the mathematical property of complex functions by which their asymptotic behaviour changes for different regions of the complex space. For example, f(z) = ez, diverges in modulus for Re[z]→+∞, but tends to zero for Re[z]→ −∞. Thus, the complex plane is separated by the imaginary axis, where the function is purely oscillatory, into two halves, in which the asymptotic behaviour of this function is completely different. The line which separates two different regions of different convergence is called the Anti-Stokes line. In the relevant case, the function tends asymptotically to eiΩθin the real axis. However by assuming θis a complex number, for Re[θ]≫1,Φ(θ)→eiΩθ. In this case, the Anti-Stokes line, assuming that we are far away from the origin, would be given by a complex line that is parallel to Re[iΩθ]=0line. Thus, if we assume θgoes along a straight line in complex space, it could eventually cross the Anti-Stokes line and boundary conditions would become convergent. 44 CHAPTER 4. RESULTS Indeed, it is easily seen that the eigenfunction is even in Fourier space. However, recalling that the Fourier Transform and Inverse Fourier Transform are only defined for real values of θ, the eigenfunction with complex θcannot be used directly for obtaining the solution in real space. Thus, the equation will be solved again for θ∈ R while using the complex eigenvalue Ωcomputed before. The modulus of the resulting solution is shown in the following figure: Figure 4.3: Representation of modulus of the eigenfunction in real space when assuming θ∈ R. In the first place, the function in real space is even, which seems to indicate that the obtained eigenvalue is also valid for the real space solution. However, solving the ODE for other arbitrary Ωalso yields even functions in real space, and while not present in the figure, the solution is actually a complex function: ˆ ϕ∈ C. These results are due to the fact that while the ODE has some intrinsic parity in real space, the parity condition ∂θΦ(θ= 0) = 0 in θ∈ R Fourier space for our solution is violated. For a function in Fourier space to be even and real, Φ(θ) = Φ(−θ) = Φ(θ), the function in real space must also be even and real ˆ ϕ(x) = ˆ ϕ(−x)∈ R, which isn’t the case for this eigenfunction, where ∂θΦ(θ=0) Φ(θ=0) ≈ −3. 45 The proof of the the even-parity property is the following: ˆ ϕ(x) = ˆ ϕ(−x)∈R=⇒Φ(x) = Z∞ −∞ dθ Re[ϕ(x)eiθx]! =Z∞ 0dθ ϕ(θ) cos(θx)∈ R (4.4) Still, by taking the modulus of the function, the solution correctly describes the behaviour of the Shear Alfv´ en Wave in the magnetohydrodinamic model: two peaks for the amplitude of the perturbed potential ϕ1can be seen at x≈ ±0.68, which correspond to the the singular points where P= 0 in (2.104), as there is a considerable energy accumulation compared to the rest of the plasma. This can be explained by the fact that the SAW has a local complex resonance frequency given by p2=−k2 ∥(x)v2 A(x)at a certain point in the plasma. In equation (3.3), assuming λ1= 0 as the function is even, there is a singularity in the real part of the ODE if x2=−2˜ λ λ2≡Ω2=⇒x=±Re[Ω], which corresponds with the maxima of the numerical results, but instead of getting a singular point, due to a the discrete nature of the numerical method, a sharp maximum is obtained. The harmonic oscillator analogue serves again as a good analogy: if we have a undampened harmonic oscillator of a certain frequency, and excite it with a force at the same frequency and phase, it will start absorbing energy indefinitely. This is a simple explanation of what is happening at the singularity point. In this case, the driving force analogue responsible for the energy accumulation are the fast compressional waves. Even-though the incompressible plasma assumption has been made, there is still some coupling between the SAW waves and compressional waves, which have velocities perpendicular to one another if |k∥|≪|k⊥|. Thus, in the xaxis, in which the SAW propagates, there is a surface which can accept an unlimited amount of energy, which is given by a perturbation coming from the yand zaxes. If Figure 4.2 can be generalized as the y= 0 plane in a two-dimensional plot, the compressional waves would propagate along the direction perpendicular to the paper. Chapter 5 Summary and discussion Starting from the MHD equations, the Alfv´ en-Interchange equation was successfully obtained, which allows the study the different perturbations in a low-pressure plasma system. After performing some geometrical and physical simplifications in order to isolate the physical phenomenon of interest, it was shown that behaviour is explained by the magnetohydrodinamic model: the Shear Alfv´ en Wave, due to its spectrum, which depends on the x-spatial coordinate, has a resonance point where we can observe an energy build up. The same numerical simulation was performed in reference [1], where the equation was also resolved for θ∈ C on the same straight line, and while the perturbation in the the real plane isn’t shown, the eigenvalue obtained for the same boundary conditions was, up to three significant digits, Ωbib = 0.681 −0.616i. Compared to the obtained value Ω = 0.6814 −0.6178, there is a very small deviation in its imaginary value of less than 1%. This deviation is expected, due to the fact that our boundary conditions at the origin and infinity are approximated, and that the numerical algorithm has some error due to its discrete nature as well. Thus, there is a certain arbitrariety in the choice of the approximated boundary conditions and algorithm which can explain this discrepancy. 46 47 Still, the obtained eigenvalue has been found to give rise to a overall complex solution when assuming θ∈R, but it doesn’t pose a problem for the physical interpretation of the phenomenon. However, from references [3] and [2] it can be seen that the empirical observations don’t correspond to the numerical predictions presented here. This is due to the fact that we have assume a plasma with no dissipation processes, which allows an infinite energy build up, as in a nondampened harmonic oscillator which is driven by an external force, which with time oscillates more and more, reaching infinite amplitude. In reality, as it can be read in [3], there are dampening mechanisms, which can be explained by the interaction of the collective MHD perturbation with the plasma particles, such as the aforementioned Landau Damping, by which the Shear Alfv´ en Wave to accelerates particles that move slightly slower than the wave. This phenomena creates an energy transfer from the wave to the particles, which carry away the energy built up at the singularity, making it finite. The incorporation of the terms responsible for these dissipation mechanisms is beyond the scope of this thesis, as it requires combining this fluid model with a discrete particle model. Nonetheless, the formulation of the MHD equations in terms of the perturbed potential leads nicely to these corrections, which could give place to an expansion of this thesis later on. Chapter 6 Appendix 6.1 Vector algebra and calculus identities Let A,Band Cbe differentiable vector fields and fa differentiable scalar field. A×B=−B×A(A1) A×B×C= (A·C)B−(A·B)C(A2) (A×B)·(C×D) = (A·C)(B·D)−(A·D)(B·C)(A3) ∇×(fv) = ∇f×v+f∇×v(A4) ∇·(v×u) = u·(∇×v)−v·(∇×u)(A5) ∇⊥f=−b×b×∇f(A6) 48 Bibliography [1] C.Z Cheng, Liu Chen, and M.S. Chance. “High-n Ideal and Resistive Shear Alfv´ en Waves in Tokamaks”. In: Annals of Physics (1985). [2] L. Chen and F. Zonca. “Theory of Shear Alfven Waves in Toroidal Plasmas”. In: Physica Scripta T60 (1995). [3] G. Vlad, L. Chen, and F. Zonca. “Dynamics of Alfv´ en Waves in Tokkamaks”. In: IL NUOVO CIMENTO (2003). [4] J.P. Goedbloed and S. Poedts. Principles of Magnetohydrodinamics: With Applications to Laboratory and Astrophysical Plasmas. Cambridge University Press, 2004. Chap. 10. [5] P.M. Bellan. Fundamentals of Plamsa Physics. Cambridge University Press, 2006. [6] R.B. White. The Theory of Toroidally Confined Plasmas. 2nd ed. Imperial College Press, 2006. Chap. 9. [7] J.P. Freidberg. Ideal MHD. Cambridge University Press, 2014. [8] L. Chen and F. Zonca. “Physics of Alfven waves and energetic particles in burning plasmas”. In: Reviews of Modern Physics 88 (2016). [9] L. Chen and H. Chen. “Gyrokinetic theory of low-frequency electromagnetic waves in finite-b anisotropic plasmas”. In: Phys. Plasmas 28 (2021). [10] Haotian Chen. “Notes on SAW in ideal MHD”. In: (2022). [11] P. Ricci. Lecture Notes in Plasma Physics and Application to Fusion Energy, Astrophysics and Industry. 49