Searching for Solitons in theMagnetosphere of a MagnetarbyMelody Candace WongA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFBACHELOR OF SCIENCEinThe Faculty of Undergraduate Studies(Physics and Astronomy)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)April 2013c© Melody Candace Wong 2013AbstractWe describe a non-pertrubative method of searching for solitary waves travelling through themagnetosphere of a magnetar. The wave will be supported by the combined effects of nonlinearityand dispersion caused by the presence of a quantum electrodynamic (QED) vacuum in a stronglymagnetized field and by a strongly magnetized plasma of the magnetosphere. Using this method,we have found a soliton in the form of an infinite current sheet. The method and the results ofthis paper could be conducive to research studying the emission of strongly magnetized stars.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Pulsar Emission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Method for Finding Nonlinear Wave Solutions in the Magnetosphere . . . . . . . . 41.3 Project Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1 Soliton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1.1 Nonperturbative Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Magnetosphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2.1 Phase and Group Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2.2 Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2.3 Nonlinearity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Mathematical Detail . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83.1 Usefulness of S = x− vt parameter . . . . . . . . . . . . . . . . . . . . . . . . . . 83.2 Storing Variables in Terms of Magnetic Vector Potential and Electric Potential . . 93.3 Setting Up Equations for Computation . . . . . . . . . . . . . . . . . . . . . . . . 13iiiTable of Contents3.3.1 Writing Equations Explicitly in Ay, Az, and φ . . . . . . . . . . . . . . . . 133.3.2 Turning on Ay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163.3.3 Turning on φ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193.3.4 Turning on Ay and Az . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214.1 Turning on Ay with Phase Velocity Greater Than C . . . . . . . . . . . . . . . . . 214.2 Turning on Ay with Phase Velocity Less Than C . . . . . . . . . . . . . . . . . . . 234.2.1 Soliton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 254.3 Turning on Ay and Az Values, Circumpolar . . . . . . . . . . . . . . . . . . . . . . 265 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29AppendicesA Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30ivList of Figures1.1 The intensity emissions of 200 pulses. Each horizontal line in the square corre-sponds to a pulse. The pulses are lined up by their central beam intensity. Noticethat the intensity substructures on either side of the central intensity appears insome pulses, while disappearing in others.[4] . . . . . . . . . . . . . . . . . . . . . . 21.2 Intensity emission from Pulsar B1857-26 presented in carousel configuration. It isthe reconstructed intensity emission of a slice of the pulsar beam, cut perpendicularto the beam axis. Black line indicates the position of one pulse. A series of pulsesare needed to form this carousel configuration. [4] . . . . . . . . . . . . . . . . . . 34.1 Phase diagram of Ay with a vp > c and a small input. Notice how circular it iswhen the input is small. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214.2 Phase diagram of Ay with a vp > c and a large input. Notice that when input islarge, it is an ellipse and it is more squared out. . . . . . . . . . . . . . . . . . . . . 224.3 Phase diagram of Ay with a vp < c and a small input. . . . . . . . . . . . . . . . . 234.4 Phase diagram of Ay with a vp < c and a large input. . . . . . . . . . . . . . . . . 234.5 Soliton is outlined in red. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 244.6 The discontinuity of A′y. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 254.7 Discontinuity of the magnetic field. . . . . . . . . . . . . . . . . . . . . . . . . . . . 264.8 Circularly polarized wave, where vp > c. Small initial input. . . . . . . . . . . . . . 264.9 Circularly polarized wave, where vp > c. Large initial input . . . . . . . . . . . . . 27vAcknowledgementsI am grateful to my supervisor, Dr. Jeremy Heyl, for the opportunity to have worked on thisresearch project with him. His expectional patience and guidance have made this paper possible.viChapter 1Introduction1.1 Pulsar EmissionMitra and Rankin observed the emission from Pulsar B1857-26 [4]. Each time the pulsar dida sweep along the observer’s line of sight, a pulse emission was recorded. In a data set of 200pulses, where the pulses were aligned by the emission’s central beam intensity, they noticed thatsubstructures to the left and the right of the main intensity emission would “disappear” in onepulse and then “reappear” in another (see Figure 1.1). This was called the nulling effect. Mitraand Rankin then presented these pulses in a carousel configuration (see Figure 1.2). By assumingthat that the pulsar beam slowly rotated around its beam axis, each beam sweep across theEarth’s line of sight, provided data along a slice of the pulsar beam cone. By artificially placingthe pulse intensities to where they may have originiated across the pulsar beam, Mitra and Rankinconstructed a picture of a slice of the pulsar beam, cut perpendicular to its beam axis - carouselconfiguration. The high intensity, central structure of the emission corresponded to the area nearthe main axis of the pulsar beam, while the intensity substructures corresponded to a ring of 20smaller beamlets surrounding the central beam.It is unclear why these beamlets are so regularly spaced apart from one another or why theyexist. We suggest that perhaps these beamlets are in actuality all one and the same. That is,a single, solitary electromagnetic wave is rotating around the pulsar beam axis. As the wavetravels, its presence can alter the local electric and magnetic field, and apply a force on the localions. By affecting the locations of the ions, the wave alters the locations of where radiation isemitted. If there is a soliton rotating around the pulsar beam axis, as it travels, it will go in andout of our line of sight. Hence if an intensity substructure is a soliton, the nulling effect wouldbe the soliton leaving our line of sight.11.1. Pulsar EmissionFigure 1.1: The intensity emissions of 200 pulses. Each horizontal line in the square correspondsto a pulse. The pulses are lined up by their central beam intensity. Notice that the intensitysubstructures on either side of the central intensity appears in some pulses, while disappearing inothers.[4]21.1. Pulsar EmissionFigure 1.2: Intensity emission from Pulsar B1857-26 presented in carousel configuration. It isthe reconstructed intensity emission of a slice of the pulsar beam, cut perpendicular to the beamaxis. Black line indicates the position of one pulse. A series of pulses are needed to form thiscarousel configuration. [4]31.2. Method for Finding Nonlinear Wave Solutions in the Magnetosphere1.2 Method for Finding Nonlinear Wave Solutions in theMagnetosphereIn a study presented by Mazur and Heyl, they were interested in stabilizing electromagneticwaves in the magnetosphere. They were able to find stabilized pulse trains in the magnetosphereby integrating over 15 coupled nonlinear ordinary differential equations (ODE) [3]. We are alsointerested in stabilizing waves in the magnetosphere, in particular we are interested in solitons.We will present a computationally simpler method of finding stabilized waves by integrating overthree nonlinear ODEs.1.3 Project OverviewIn this study, we are searching for mathematical solutions that could indicate the possible exis-tence of solitons in the magnetosphere. Our strategy will be to imagine an already formed wavetravelling through this medium and then try to figure out how to best maintain the shape of thistravelling wave. We decided that we could stabilize this wave using nonlinearity and dispersioneffects. That is, the effect of the strongly magnetized plasma and the effect of the QED vacuumin a strongly magnetized field, respectively. As we will further discuss in Sections 2.2.2 and 2.2.3,these effects when taken individually, would cause a wave to change shape. However, when theseeffects are put together, there may be cases where the effects can counter one another and endup supporting a wave - this is what we are looking for.Our study will make an ansantz that the travelling wave is a plane wave with the parameterS = x− vt. The variable x refers to position, as v and t refer to velocity and time, respectively.We will assume that our solitary wave is propagating in a homogenous medium, where we willneglect the effects of gravity and the magnetic field forces on the plasma medium. In addition, wewill describe our travelling electromagnetic wave in terms of the magnetic vector potential andthe electric potential in order to simplify computations. We are doing this in an effort to explorethe possibility of a single, solitary electromagnetic wave surviving the magnetosphere.4Chapter 2Theory2.1 SolitonIt is a single, solitary wave that has a finite extent in space and time. The wave has a start andan end, with an amplitude that approaches zero far from the peak of the wave. As it travelsthrough a medium, the soliton will maintain its shape.2.1.1 Nonperturbative MethodWe are making an ansatz that our soliton is a plane wave with a parameter S = v − xt. Thisparameter choice will help to simplify computations (see Section 3.1) and make the solutions thatwe derive, nonperturbative.Since the electromagnetic wave is travelling through a nonlinear medium (see 2.2.3), it is keythat we use a nonperturbative method of representing our solution. That is, we will be solving oursolutions directly, rather than using perturbative methods of linearly approximating the solution.Perturbative methods, such as using Fourier Theory to approximate a wave, will not provide ageneralized solution to nonlinear problems. They will only work for limited cases as we will needto solve these nonlinear equations computationally and can only integrate over a finite numberof Fourier modes.2.2 MagnetosphereA magnetar is a pulsar with an exceptionally strong magnetic field on the scale of 1011 Teslas [5].Due to the quickly changing magnetic fields near the spinning magnetar, a strong electric field isinduced on to the surface of the star. This electric field works against gravity, and pulls electrons52.2. Magnetosphereup from along the poles, forming the magnetosphere.2.2.1 Phase and Group VelocityIf we imagine a set of harmonic waves, each travelling at different velocities, these velocities willbe referred to as phase velocities, vp. If we add up the amplitudes of these waves while they aretravelling, they form a wave packet that will travel at its own charactersitic velocity. This is thegroup velocity, vg. This is not surprising as the position at which the wave crests coincide is itselfa location that alters with time [6].vp =ωk(2.1)vg =∂ω∂k(2.2)Note that group velocity cannot exceed the speed of light, while phase velocity can. This isbecause phase velocity cannot carry information, while group velocity can.2.2.2 DispersionAs will later be seen in the Section 3.2, the equations will use the plasma frequency, ωp. If someof the plasma gets displaced, the Coulomb force will cause the particles in the plasma to oscillateat this frequency ωp. Since this value is set based on the medium, then from the definition ofthe phase velocity (see Equation 2.1), waves of different spatial frequencies will travel at differentspeeds. If we imagine that our solitary wave is made up of a combination of different sinusoidalwaves, it is easy to imagine that the wave will disperse as different spatial frequency componentspropagate at different speeds.2.2.3 NonlinearityNonlinearity occurs when the dipoles of the medium no longer respond the alternating electricfield of our wave in a linear fashion. This nonlinearity is generated by the radiative correctionsof QED. In this case, the polarization has a nonlinear relationship with field.62.2. MagnetosphereThe significance is more obvious by seeing the relationship,D = E + P(E,B) (2.3)H = B−M(E,B) (2.4)In the magnetosphere, the medium that causes this nonlinearity is the QED vacuum in thepresence of a strong magnetic field. The magnetic field causes the QED vacuum to becomemagnetized and polarized[1]. This effect can become so powerful that it can alter the shape ofa wave as it is travelling through this magnetized vacuum. More specifically, this medium addsnonlinearity to the wave such that different amplitudes of the wave will travel at different speeds,much like a surfer’s wave where the top of the wave overtakes the base the wave. This nonlinearitywill cause the wave to steepen and shock.In our equations, we will account for this nonlinearity in our vacuum dielectric and inversemagnetic permeability tensors, ij and µ−1ij respectively. These tensors were derived from thefollowing electromagnetic Lagrangian relationships [1].D = E + P =∂L∂E(2.5)H = B−M = −∂L∂B(2.6)Note that rather than using the full, generalized equation, we are using the weak field ap-proximation, where the field is less than Bk = 4.413 × 1013 G. The reason for this is that thegeneralized equation will lead to complicated tensors that will bring difficulties in the computa-tion. So in this study we will test our computational method with the weak field approximation.Futher studies can then apply this computational method to the generalized tensors.7Chapter 3Mathematical DetailHere we will be rearranging variables in Maxwell’s equations into terms of magnetic vector po-tential and electric potential. The purpose of this is to minimize the number of variables that weneed to keep track of in order to simplify our compuations later on. After redefining the variablesin Maxwell’s equations into values of potentials, we can manipulate the equations and solve forthe values of magnetic vector potential and electric potential. From these values, we can thenfigure out the other variables in Maxwell’s equations.We will begin by presenting the mathematical benefit of using S = x − vt as its propertieswill be applied to our derivation.3.1 Usefulness of S = x− vt parameterOur choice in using the parameter S = x − vt allows us to solve ordinary differential equations(ODE) instead of partial differential equations (PDE). To demonstrate,∂f(S)∂x=dfdSdSdx=dfdx(3.1)∂f(S)∂t=dfdSdSdt= −vdfdx(3.2)Here we have transformed the PDEs that had partial derivatives with respect to x and to t intoequations of ODE. Note that the derivatives with respect to y and to z will go to zero. Wewill frequently use these parameter S properties as we manipulate Maxwell’s equations in thefollowing section 3.2.83.2. Storing Variables in Terms of Magnetic Vector Potential and Electric Potential3.2 Storing Variables in Terms of Magnetic Vector Potentialand Electric PotentialSince we are interested in an electromagnetic wave, we begin our derivations from Maxwell’sequations in matter, in Gaussian units. We are using Gaussian units so that the units for electricand magnetic fields will be the same.Recall that the variables are in terms of the parameter S = x − vt, and as such, the partialderivatives with respect to x and t are outlined in equations 3.1 and 3.2. Also note that in ournotation, primes indicate derivatives with respect to S (or as we now have shown in equation 3.1,with respect to x, as well).∇ ·B = 0 (3.3)∇×E = −1c∂B∂t(3.4)∇×H =4picJf +1c∂D∂t(3.5)∇ ·D = 4piρf (3.6)We will now proceed to describe the Maxwell variables in terms of potentials. Taking the firstMaxwell equation, 3.3, we make use of the magnetic field divergence being zero.B = ∇×A + BoB = −A′zy +A′yz + Bo (3.7)Taking the second Maxwell’s Equation 3.4 and applying 3.7,∇×E = −1c∂∂t(∇×A)∇×E = ∇× (−1c∂∂tA)0 = ∇× (E +1c∂∂tA) (3.8)93.2. Storing Variables in Terms of Magnetic Vector Potential and Electric PotentialSince the curl of E + 1cδδtA is zero, then we can set the term to equal the negative gradient of ascalar field φ. We will choose a scalar field such that Ax = 0 by Coloumb Gauge theory.−∇φ = E +1c∂A∂S∂S∂tE =vc∂A∂x−∇φE =vcA′ − φ′x (3.9)For the third Maxwell Equation 3.5,0 = ∇×H−1c∂D∂S∂S∂t−4picJ0 = −H ′zy +H′yz +vcD′ −4picJ4picJ = −H ′zy +H′yz +vcD′ (3.10)Focusing on the x, y and z components of equation 3.10 individually, we find that4picJx =vcD′x (3.11)4picJy = −H′z +vcD′y (3.12)4picJz = H′y +vcD′z (3.13)As will be shown in equation 3.16, there is a relationship between the current density J and theelectric field. Utilizing equation 3.9, we can then define the components of the current density interms of magnetic potential..In the final Maxwell Equation 3.6, we apply equation 3.11 and findρ =D′x4piρ =14pi4piJxvJx = vρ (3.14)103.2. Storing Variables in Terms of Magnetic Vector Potential and Electric PotentialIn order to close the system, we will also apply Lorentz Force Equation. First, we start withthe change in current density with respect to time and use the relationship J = ρv. The variablesJ, ρ, and v represent volume current density, charge density, and velocity respectively.∂J∂t=∂J∂S∂S∂t∂(ρv)∂t= −v∂J∂Sρ∂v∂t= −vJ′ (3.15)Now we will use the Lorentz Force Equation, F = q(E− v ×B). By assuming that the velocityof the charges in the plasma are close to zero, v ×B → 0, we find that ma = qE. Returning toequation 3.15,−vJ′ = ρ(qEm)= nq(qEm)=nmq2E=ω2p4piE (3.16)Above, we used definition of the charge density, ρ = nq, where q and n are the charge and theparticle density, respectively. We then applied the relationship nq2m =ω2p4pi , where ωp is the plasmafrequency.Here we take the x-component of the equation 3.16 and plug in the electric field from 3.9. We113.2. Storing Variables in Terms of Magnetic Vector Potential and Electric Potentialmake use of the fact that we have set Ax = 0, and then we apply equation 3.11.−vJ ′x =ω2p4pi(vcA′x − φ′)−vJ ′x = −ω2p4piφ′vJx =ω2p4piφv24piD′x =ω2p4piφv2D′x = ω2pφ (3.17)By applying the electric field relationship in equation 3.9 onto equation 3.16, we get −vJ′ =ω2p4pi (vcA′−φ′x). From this, we can find the x, y, and z components of J. In addition, if we includethe relationships 3.12, and 3.13, we find−H ′z +vcD′y =4picJy = −ω2pc2Ay (3.18)H ′y +vcD′z =4picJz = −ω2pc2Az (3.19)Finally, we will use the definition of the fields.Di = ijEj (3.20)Hi = µ−1ij Bj (3.21)For this study, we will apply the weak field limit of ij and µ−1ij (recovered from Klein andNigam’s paper [2]). A more generalized field can be found in a paper by Heyl and Hernquist [1].ij = δij +145piαB2k[2(E2 −B2]δij + 7BiBj ] (3.22)µ−1ij = δij +145piαB2k[2(E2 −B2)δij − 7EiEj ] (3.23)123.3. Setting Up Equations for ComputationWe have now redefined all of the necessary variables into terms of the magnetic vector andelectric potentials. More specifically, into values of magnetic potential components Ay and Az,and the electric potential φ.3.3 Setting Up Equations for ComputationIn this section we will be expanding our previous equations (from 3.2) and will solve for Ay, Az,and φ numerically with the Runge Kutta 4 method. A copy of this numerical solver is attachedA.As we are dealing with many different terms, Maple will be used to do the algebraic manip-ulations. Also note that we will be making the equations dimensionless by setting identities forthe length and the field strength.3.3.1 Writing Equations Explicitly in Ay, Az, and φHere we will rewrite the components of H’ and D’ into terms of Ay, Az, and φ. This is necessaryso that on both sides of the equation, for the equations 3.18, 3.19, and 3.17, will be written asvalues of the potentials. Later, these rewritten versions of 3.18, 3.19, and 3.17 will be used tonumerically solve for Ay, Az, and φ.From the field definitions in 3.20 and 3.21, we see that the components of H’ and D’ aredefined by the fields B and E, and the matrices ij and µ−1ij . Recall from 3.7 and 3.8 thatB = 〈B0,x, B0,y − A′z, B0,z + A′y〉 and E = 〈−φ′, vcA′y,vcA′z〉. The extended form of the matricesij and µ−1ij are shown below (taken from [2]).ij =1 + α45piB2k(2(E2 −B2)+ 7B2x)α45piB2k(7BxBy) α45piB2k(7BxBz)α45piB2k(7BxBy) 1 + α45piB2k(2(E2 −B2)+ 7B2y)α45piB2k(7ByBz)α45piB2k(7BxBz) α45piB2k(7ByBz) 1 + α45piB2k(2(E2 −B2)+ 7B2z)(3.24)133.3. Setting Up Equations for Computationµ−1ij =1 + α45piB2k(2(E2 −B2)− 7E2x)− α45piB2k(7ExEy) − α45piB2k(7ExEz)− α45piB2k(7ExEy) 1 + α45piB2k(2(E2 −B2)− 7E2y)− α45piB2k(7EyEz)− α45piB2k(7ExEz) − α45piB2k(7EyEz) 1 + α45piB2k(2(E2 −B2)− 7E2z)(3.25)The extended form of the matrices are still not explicitly written in terms of the potentials.Here we will resolve that by taking the products of E and B relevant to the matrices and findingtheir values in terms of potentials.From dot product of 3.8E2 = φ′2 +(vc)2A′2y +(vc)2A′2z (3.26)From dot product of 3.7B2 = B20,x +(B0,y −A′z)2+(B0,z +A′y)2(3.27)Subtracting Equation 3.27 from 3.26E2−B2 = φ′2+[(vc)2− 1]A′2y +[(vc)2− 1]A′2z −B20,x−B20,y−B20,z+2B0,yA′z+2B0,zA′y (3.28)Taking the derivative of 3.28 with respect to x[E2 −B2]′= 2φ′φ′′ + 2[(vc)2− 1]A′yA′′y + 2[(vc)2− 1]A′zA′′z + 2B0,yA′′z + 2B0,zA′′y (3.29)Now that B,E, ij , and µ−1ij are explicitly defined by Ay, Az, and φ, we can find the componentsof D and H using relationships 3.20 and 3.21.Note that the equations below are split into three rows in order to increase readability. Eachrow in the Di equation represents a cell in the matrix product ijEj . Likewise, each Hi equation143.3. Setting Up Equations for Computationis split into three rows so that each row corresponds to a cell µ−1ij Bj .D′x =− φ′′ −α45piB2k(2[E2 −B2]′φ′ + 2[E2 −B2]φ′′ + 7B20,xφ′′)+7α45piB2kvcB0,x[−A′′zA′y + (B0,y −A′z)A′′y]+7α45piB2kvcB0,x[−A′′yA′z + (B0,z −A′y)A′′z](3.30)D′y =7α45piB2k[B0,xA′′zφ′ +(B0,xA′z −B0,xB0,y)φ′′]+vcA′′y +α45piB2kvc[2[E2 −B2]′A′y + 2[E2 −B2]A′′y − 14(B0,y −A′z)A′′zA′y + 7(B0,y −A′z)2A′′y]+7α45piB2kvc[(B0,yA′′z − 2A′zA′′z) (B0,z +A′y)+(B0,yA′z −A′2z)A′′y](3.31)D′z =−7α45piB2kB0,x[φ′′(B0,z +A′y) + φ′A′′y]+7α45piB2kvc[−A′′z(B0,zA′y +A′2y ) + (B0,y −A′z)(B0,zA′′y + 2A′yA′′y)]+vcA′′z +α45piB2kvc[2[E2 −B2]A′′z + 2[E2 −B2]′A′z + 14(B0,z +A′y)A′′yA′z + 7(B0,z +A′y)2A′′z](3.32)H ′x =α45piB2k(2[E2 −B2]′ − 14(vc)2φ′φ′′)B0,x+7α45piB2kvc[φ′′(B0,yA′y −A′yA′z) + φ′(B0,yA′′y −A′′yA′z −A′yA′′z)]+7α45piB2kvc[φ′′(B0,zA′z −A′yA′z) + φ′(B0,zA′′z −A′′yA′z −A′yA′′z)](3.33)H ′y =7α45piB2kvc[φ′′A′y + φ′A′′y]B0,x−A′′z +α45piB2k[2[E2 −B2]′(B0,y −A′z)− 2[E2 −B2]A′′z − 14(vc)2A′yA′′y(B0,y −A′z) + 7(vc)2A′2y A′′z]−7α45piB2k(vc)2 [A′′z(B0,zA′y +A′2y ) +A′z(B0,zA′′y + 2A′yA′′y)](3.34)153.3. Setting Up Equations for ComputationH ′z =7α45piB2kvc[φ′′A′z + φ′A′′z]B0,x−7α45piB2k(vc)2 [A′′yA′zB0,y +A′yA′′zB0,y −A′′yA′2z − 2A′yA′zA′′z]+A′′y +α45piB2k[(2[E2 −B2]− 7(vc)2A′2z)A′′y +(2[E2 −B2]′− 14(vc)2A′zA′′z)(B0,z +A′y)](3.35)Now that H and D are explicitly defined by the magnetic and electric potentials, we canwrite the equations of 3.18, 3.19, and 3.17 entirely in terms of potentials. In this way we canmanipulate the equations so that we can solve for Ay, Az, and φ, using nonlinear ODEs.As our equations are nonlinear and have many different terms, finding a solution can bedifficult. So instead, we will be turning on our variables one at a time. That way, we can slowlybut steadily build our way up to the full solution.3.3.2 Turning on AyIn this section, we are only concerned with Ay, its derivatives and the constants in the equation.Our goal here is to find a general solution to our nonlinear Ay equation. With this generalsolution, we will try to find a case where a soliton might exist. Plugging in 3.35 and 3.31 intoequation 3.18, we find that−(ωpc)2Ay =−A′′y −α45piB2k[2(E2 −B2)A′′y + 2(E2 −B2)′A′y]+(vc)2A′′y +α45piB2k(vc)2[2(E2 −B2)′A′y + 2(E2 −B2)A′′y] (3.36)If we only focus on the Ay variables in 3.28 and 3.29, we find E2 − B2 =[(vc)2− 1]A′2y and2(E2 −B2)′= 2[(vc)2− 1]A′yA′′y. Therefore, equation 3.36 becomes−(ωpc)2Ay =[(vc)2− 1]A′′y +6α45piB2k[(vc)2− 1]2(A′y)2A′′y−(ωpc)2∣∣∣∣(vc)2− 1∣∣∣∣−1Ay =A′′y +6α45piB2k∣∣∣∣(vc)2− 1∣∣∣∣ (A′y)2A′′y (3.37)We can simplify the equation by making it dimensionless. We will set the term 2α45piB2k3∣∣∣(vc)2− 1∣∣∣,163.3. Setting Up Equations for Computationwhich is a measurement of field strength, to 1. Likewise, we will set(ωpc)2∣∣∣(vc)2− 1∣∣∣−1, which isin units of length−2, also to 1.Note that we are currently taking the field strength and length−2 terms as absolute values.In order to account for all cases of the equation 3.37, we need to consider when these terms arepostive and when they are negative. The sign of these terms will depend on whether the phasevelocity v is greater or less than the speed of light in a vacuum, c. The case when v = c will makethe equation a trivial problem as the field strength and length−2 terms will be zero.For now, we will just focus on the case when the phase velocity is greater than c. In particular,that the terms are set to 1.−Ay = A′′y +A′2y A′′y−Ay = A′′y(1 +A′2y )A′′y = −Ay1 +A′2y(3.38)In this case, with just one variable Ay and its derivatives, we can solve for Ay analytically.We shall do so right now.Letx = Ay,dxdt= A′y = yy = A′y,dydt= A′′y = y′ (3.39)173.3. Setting Up Equations for ComputationThen we can rewrite equation 3.38 as the expressiondydx=dydtdtdxdydx= −x1 + y21ydydx= −xy + y3(3.40)∫(y + y3)dy = −∫xdxy22+y44= −x22(3.41)Removing the substitions from 3.39, we find(A′y)2 +(A′y)42= −A2y (3.42)From the equation 3.40, we can use the information to make a phase diagram with A′y andAy for the case when the phase velocity is greater than c.dA′ydAy= −AyA′y + (A′y)3(3.43)In the case of when the phase velocity is less than c, the field strength and the length−2 termsbecome negative. As such, we can test this case by setting these terms as -1. We find thatdA′ydAy=AyA′y − (A′y)3(3.44)(A′y)2 −(A′y)42= A2y±[(A′y)2 −(A′y)42] 12= Ay (3.45)Due to the symmetry between Ay and Az, this method is also applicable to Az.183.3. Setting Up Equations for Computation3.3.3 Turning on φSimilar to section 3.3.2, we will only look at φ and its derivatives; the other variables will be setto zero. We will solve for φ using equation 3.17 and 3.30. Putting these equations together, wefind−ω2pv2φ = φ′′ +6α45piB2kφ′2φ′′ (3.46)If we set the constants −ω2pv2 and6α45piB2kto one, we see that this equation will be identical to theexpression in 3.38. Hence the solution found in section 3.3.2 will also be applicable to φ.3.3.4 Turning on Ay and AzHere we are allowing mixing between the variables Ay and Az, and its derivatives. Using equations3.18 and 3.19, and applying 3.34, 3.35, 3.31, and 3.32, we find−(wpc)2Ay =[(vc)2− 1]A′′y +M1A′′yA′2z + 2M1A′yA′zA′′z + 3M1A′2y A′′y (3.47)−(wpc)2Az =((vc)2− 1)A′′z +M1A′2y A′′z + 2M1A′yA′′yA′z + 3M1A′2z A′′z (3.48)where M1 isM1 = −4αv2c2 − 2αc4 − 2v4α45c4piB2k(3.49)We see that if we switch the y and z subscripts of one of the equations, then the equationswill become identical. As such, we see that the equations are perfectly symmetic.Here we take the equation 3.47 and divide both sides by ‖(vc)2− 1‖, being careful with theabsolute values. Recall in Section 3.3.2 that we set the term(wpc)2‖(vc)2− 1‖−1 = 1, here wewill do it again. Below, we will pull MA′′yA′2z + 2MA′yA′zA′′z together as [MA′yA′2z ]′ and factor thenumerator of M1.−Ay =v2 − c2‖v2 − c2‖A′′y + 3M2A′2y A′′y + [M2A′yA′2z ]′ (3.50)M2 =2α45c2piB2k(v2 − c2)2‖v2 − c2‖(3.51)193.3. Setting Up Equations for ComputationSince Ay and Az are symmetric in this section, the equation for −Az is just 3.50 with the yand z subscripts switched.20Chapter 4ResultsUsing the equations from Chapter 3 and plugging them into our nonlinear ODE solver (see A),we were able to produce phase diagrams of our equations. From these diagrams we can infer thecharacteristics of a wave and determine whether a soliton can be found with particular solution.Note that based on making the equations dimensionless (see Section 3.3.2), we realize thatwe need to consider both the case when the phase velocity, v, is less than c and the case when itis greater than c. As well, recall in Chapter 3, we were turning the variables on one at a time.So in the same fashion, our results are grouped by the variables turned on in the equation andby the phase velocity.4.1 Turning on Ay with Phase Velocity Greater Than CFigure 4.1: Phase diagram of Ay with a vp > c and a small input. Notice how circular it is whenthe input is small.214.1. Turning on Ay with Phase Velocity Greater Than CFigure 4.2: Phase diagram of Ay with a vp > c and a large input. Notice that when input islarge, it is an ellipse and it is more squared out.When the phase velocity is greater than c, we observe that the phase diagrams of A′y versusAy are stable and continuous (see figures 4.1 and 4.2). This indicates that although the waveis supported by nonlinearity and dispersion in the magnetosphere, the wave is periodic. Since asoliton is a single, solitary wave in space, it is not periodic. Therefore, this type of solution willnot reveal a soliton.What is interesting about these phase diagrams is that we can see the affects of nonlinearity.Previously, we were accounting for the nonlinearity in the magnetosphere through the vacuumdielectric and inverse magnetic permeability tensors, ij and µ−1ij . As shown in the equations3.22 and 3.23, the nonlinearity of these tensors are dependent on the electric and magnetic fields.Also, recall that the electric and magnetic fields are dependent on A′y. Thus, the larger the initialvalue of A′y, the stronger the nonlinearity.In the figure 4.1, a small initial A′y value was put in and as such, the affects of nonlinearityis muted. The phase diagram is very round, which is similar to the phase diagrams of a linearsystem, such as a linear harmonic oscillator.In the figure 4.2, a large initial A′y value was put in, making the affects of nonlinearity apparent.This is seen by the squaring of the phase diagram.224.2. Turning on Ay with Phase Velocity Less Than C4.2 Turning on Ay with Phase Velocity Less Than CFigure 4.3: Phase diagram of Ay with a vp < c and a small input.Figure 4.4: Phase diagram of Ay with a vp < c and a large input.As discussed in the previous section, phase diagrams can help to infer characteristics of thewave. Previously when the phase velocity was greater than c, the phase diagram was stable andcontinuous. Here, we see that the phase diagram is discontinous and unstable, indicating thatthe wave is aperiodic. Hence there is possibility of finding a soliton in this type of equation.Notice that the discontinuity of the phase diagram always occurs at A′y = ±1. From equation234.2. Turning on Ay with Phase Velocity Less Than C3.44, we see that this is because at A′y = ±1, we are dividing zero and the slope becomes infinite.These discontinuous phase diagrams 4.3 and 4.4 are not very useful because they are inaccurate.The reason is that we are numerically solving for these values using Runge Kutta 4 (RK 4), amethod that depends on the slope. When the slope heads towards infinity, the values that RK 4solves for afterwards will be thrown off.Hence, in order to get a useful phase diagram, we will not use RK 4 here. Instead, we willuse the analytically solved equation 3.45. We will plot the phase diagram, where Ay is a functionof A′y. That is, we will be plugging in values for A′y and then solving for Ay. In this way, we canplot out our phase diagram, Figure 4.5.As seen in Figure 4.5, the smaller A′y inputs have a rounder shape than the larger A′y inputs.The larger values seem to be more squared off. Recall from the previous section that this is dueto the larger A′y, the stronger the nonlinearity effects.Figure 4.5: Soliton is outlined in red.244.2. Turning on Ay with Phase Velocity Less Than C4.2.1 SolitonAs seen in Figure 4.5, the outline in red is a soliton. We know this because there is an abruptchange at coordinates (0.7, 1), (0.7, -1), (-0.7, 1), and (-0.7, -1) in Figure 4.5. This can beobserved if we start from the origin of the figure and travel outwards along the red outlinetowards coordinate (0.7, 1). Along this path, A′y and Ay are both positive, which makes sense.However, continue travelling along this red outline, we see that while the derivative of Ay, thatis A′y, is positive, the Ay is heading towards the negative direction, which does not make sense.The only sensible solution is that upon reaching the coordinate (0.7, 1), the solution jumps downtowards (0.7, -1) and continues inwards towards the origin.Recall from equation 3.7, that the magnetic field is proportional to A′y, hence if there is a jumpin A′y, then there is also a jump in the magnetic field. Based on the path that we just described, wesee that A′y is steadily increasing before it makes an abrupt change towards the negative direction.In the same way, we would expect the the magnetic field over space to continuously increase beforeit suddenly switches into the other direction (see Figure 4.7). Based on the magnetic field shownin Figure 4.7, we conclude that the soliton must be an infinite current sheet that sits right wherethe magnetic field switches direction.Figure 4.6: The discontinuity of A′y.254.3. Turning on Ay and Az Values, CircumpolarFigure 4.7: Discontinuity of the magnetic field.4.3 Turning on Ay and Az Values, CircumpolarHere we are applying the equations when both Ay and Az are turned on. We will only look atthe case when the phase velocity is greater than c. (The case when the phase velocity is less thanc does not seem to produce a coherent solution.)In the figures 4.8 and 4.9, we are observing circularly polarized waves. The numerical differencebetween the two figures is that figure 4.8 started with a smaller initial input than figure 4.9.Similar to the previous sections, the solution with the larger initial input was more affected bynonlinearity. Here we see that the phase diagram that is experiencing more nonlinearity has anangular precession in a Az versus Ay plot. This precession indicates that there is energy shiftingbetween the axes. The upper and lower boundary of the phase diagram “radius” is maintained,suggesting some sort of energy conservation.264.3. Turning on Ay and Az Values, CircumpolarFigure 4.8: Circularly polarized wave, where vp > c. Small initial input.Figure 4.9: Circularly polarized wave, where vp > c. Large initial input27Chapter 5ConclusionWe discussed about a nonperturbative method of finding stabilized waves travelling in the mag-netosphere of a magnetar. A nonpertubative method was necessary as the medium the wave waspropagating through was nonlinear. This method was achieved by making an ansantz that thewave was a plane wave with the parameter S = x− vt. In addition, our method of storing valuesin terms of magnetic vector potential, electric potential and plasma density has simplified theprocess of numerically solving for stabilized waves.While travelling through the magnetosphere, the wave will maintain its shape through thecombined effects of nonlinearity and dispersion caused by the QED vacuum in a strongly mag-netized field and by the strongly magnetized plasma, respectively. We were interested in caseswhere the wave steepening effect from the nonlinearity will balance with the effect from disper-sion, allowing us to oberve a soliton. We were able to find a soliton in the form of an infinitecurrent sheet.In this study we only looked at several variables and used a weak field approximation. In thefuture, we hope to extend our search by adding more variables and using a more general form ofthe field equations.28Bibliography[1] Heyl J S and Hernquist L 1997 J. Phys. A: Math. Gen. 30 6485[2] Klein J J and Nigam B P 1964 Phys. Rev. B 135 1279[3] Mazur D and Heyl J S 2010 arXiv:1002.2915v2 [astro-ph.HE][4] Mitra D and Rankin J M 2007 arXiv:0712.1338 [astro-ph][5] Ostlie D A and Carroll B W 2007. An Introduction to Modern Stellar Astrophysics 2nd Ed.Addison-Wesley San Franscisco[6] Pedrotti F L and Pedrotti L S 1993. Introduction to Optics Prentice-Hall New Jersey29Appendix AAppendixA sample C program is that solves nonlinear ODEs and plots out a 2D plot. When initialized,it will ask the user for initial input values for the equation. It will also ask the user what twovariables the user wishes to plot./** Nonlinear ODE solver and plotter.*/#include <stdio.h>#include <math.h>#include <stdlib.h>#include <assert.h>void initialize(double *, double *, double *, double *, double *, double *, double *,double *,double *, double *, int *, int *);void derivatives(double, double *, double *, double *, double *, double, double,double, double);void rungeKutta4(double *, double *, double *, double *, int, double, double, double*, double *, double, double, double, double, void (*derivatives)(double, double *,double *, double *, double *, double,double, double, double));void output( FILE *, double, double, double *, double *, double);int main(){FILE *outputFile;30Appendix A. AppendixFILE *pipe;//***Declaring Variables***double initialY, initialDy, finalt, initialZ,initialDz,C, By,Gamma,a;double h;double plotScale;int xAxis, yAxis;double y[2], dydt[2], yOut[2];double z[2], dzdt[2], zOut[2];double t;const char *plotLabel[7];outputFile = fopen("Data.dat", "w");pipe = popen("gnuplot -persist", "w");//***Read User Input From Screen And Feed To Variables ***initialize(&initialY, &initialDy, &initialZ, &initialDz,&By, &Gamma, &a,&finalt, &h, &plotScale,&xAxis, &yAxis);y[0] = initialY;y[1] = initialDy; // dx/dt = yz[0] = initialZ;z[1] = initialDz;C=1; // (v^2-c^2)/|v^2-c^2|t=0.;//****Computing and Writing Data Values to "outputFile"****output(outputFile, h, t, y, z, plotScale);while(t<=finalt){31Appendix A. Appendixt+=h;derivatives(t,y,dydt,z,dzdt,C,By,Gamma,a);rungeKutta4(y,z,dydt,dzdt,2,t,h,yOut,zOut,C,By,Gamma,a,derivatives);y[0]=yOut[0];y[1]=yOut[1];z[0]=zOut[0];z[1]=zOut[1];output(outputFile,h,t,y,z,plotScale);}fclose(outputFile);//******************PLOTTING SECTION ********************//***Plot Labels***plotLabel[0]="h";plotLabel[1]="t";plotLabel[2]="plotScale*(y[1]+z[1])";plotLabel[3]="y[0]";plotLabel[4]="y[1]";plotLabel[5]="z[0]";plotLabel[6]="z[1]";//***Actual Plotting Using An Imbedded Gnuplot***fprintf(pipe, "plot ’Data.dat’ using %i:%i title ’inital Ay=%4.2lf, dAy=%4.2lf,Az=%4.2lf, dAz=%4.2lf’ withdots\n",xAxis,yAxis,initialY,initialDy,initialZ,initialDz);fprintf(pipe, "set xlabel ’%20s’\n",plotLabel[xAxis-1]);fprintf(pipe, "set ylabel ’%20s’\n",plotLabel[yAxis-1]);fprintf(pipe, "replot\n");close(pipe);return 0;32Appendix A. Appendix}//****************METHOD CODE ********************void initialize(double *initialY, double *initialDy, double *initialZ, double*initialDz, double *By, double *Gamma, double *a, double *finalt, double *h,double *plotScale, int *xAxis, int *yAxis){double numOfPeriod;printf("\n\nAy, Az and By. \n");printf("We have set Bk=1, c=1, w=1 \n");//***Prompting User For Input Values and Plotting Axes***printf("Please enter the initial ay, ay’, az, az’,By, Gamma,alpha, final t,step size, and plot scaler\n");scanf("%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", initialY, initialDy, initialZ,initialDz, By, Gamma,a, finalt, h, plotScale);printf("PLOTTING: Enter x and y axes of plot. Choose from the following:\n 1.h\n 2. t\n 3. plotScale*(y[1]+z[1])\n 4. y[0]\n 5. y[1]\n 6. z[0]\n 7.z[1]\n");scanf("%i %i", xAxis, yAxis);}//******Taking Derivatives*********// Probably most important section of this script. The dydt[1] and dzdt[1] are wherethe// second order derivatives (of either A_y, A_z, \phi) that we calculated in Chapter 3are placed.void derivatives(double t, double *y, double *dydt, double *z, double *dzdt, double C,double By, double Gamma,double a){double w=1;33Appendix A. Appendixdouble Bk=1;double c=1;dydt[0] = y[1]; // slope of position (aka velocity)dydt[1] = **PUT SECOND DERIVATIVE EQUATION HERE**//This is where we put the equations that we found in Chapter 3Mathematical Details.//So if y = A_y, then you would put the equation for A_y’’ here.dzdt[0] = z[1]; // slope of position (aka velocity)dzdt[1] = **PUT SECOND DERIVATIVE EQUATION HERE**//This is where we put the equations that we found in Chapter 3Mathematical Details.//So if z = A_z, then you would put the equation for A_z’’ here.}//****************RK 4******************/void rungeKutta4(double *y, double *z, double *dydx, double *dzdx, int n, double x,double h, double *yOut, double *zOut, double C, double By, double Gamma, doublea,void (*derivatives)(double, double *, double *, double *, double *, double,double, double,double)){int i;double xh, hh, h6;double *y_otherSlope, *y_newSlope, *y_next;double *z_otherSlope, *z_newSlope, *z_next;y_otherSlope = (double *) malloc(n*sizeof(double));y_newSlope = (double *) malloc(n*sizeof(double));y_next = (double *) malloc(n*sizeof(double));z_otherSlope = (double *) malloc(n*sizeof(double));z_newSlope = (double *) malloc(n*sizeof(double));z_next = (double *) malloc(n*sizeof(double));hh = h*0.5;34Appendix A. Appendixh6 = h/6.;xh = x+hh;for(i=0; i<n; i++){y_next[i] = y[i]+hh*dydx[i];z_next[i] = z[i]+hh*dzdx[i];}(*derivatives)(xh,y_next,y_newSlope,z_next,z_newSlope,C,By,Gamma,a);for(i=0; i<n; i++){y_next[i] = y[i]+hh*y_newSlope[i]; //y_newSlope has been updated inderivatives(). So y_next will changez_next[i] = z[i]+hh*z_newSlope[i]; //y_newSlope has been updated inderivatives(). So y_next will change}(*derivatives)(xh,y_next,y_otherSlope,z_next,z_otherSlope,C,By,Gamma,a);for(i=0; i<n; i++){y_next[i] = y[i]+h*y_otherSlope[i];z_next[i] = z[i]+h*z_otherSlope[i];y_otherSlope[i] += y_newSlope[i]; //save a malloc call. ie. other slopeis now otherslope + y_newSlopez_otherSlope[i] += z_newSlope[i];}(*derivatives)(x+h,y_next,y_newSlope,z_next,z_newSlope,C,By,Gamma,a);for(i=0; i<n; i++){yOut[i] = y[i]+h6*(dydx[i]+y_newSlope[i]+2.0*y_otherSlope[i]);zOut[i] = z[i]+h6*(dzdx[i]+z_newSlope[i]+2.0*z_otherSlope[i]);}free(y_otherSlope);free(y_newSlope);free(y_next);free(z_otherSlope);free(z_newSlope);free(z_next);35Appendix A. Appendix}void output(FILE *outputFile, double h, double t, double *y, double *z, doubleplotScale) {fprintf(outputFile, "%lf %lf %lf %lf %lf %lf %lf\n", h, t,plotScale*(y[1]+z[1]), y[0], y[1], z[0], z[1]);}36
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Undergraduate Research /
- Searching for solitons in the magnetosphere of a magnetar
Open Collections
UBC Undergraduate Research
Searching for solitons in the magnetosphere of a magnetar Wong, Melody Candace Apr 30, 2013
pdf
Page Metadata
Item Metadata
Title | Searching for solitons in the magnetosphere of a magnetar |
Creator |
Wong, Melody Candace |
Date Issued | 2013-04 |
Description | We describe a non-pertrubative method of searching for solitary waves travelling through the magnetosphere of a magnetar. The wave will be supported by the combined effects of nonlinearity and dispersion caused by the presence of a quantum electrodynamic (QED) vacuum in a strongly magnetized field and by a strongly magnetized plasma of the magnetosphere. Using this method, we have found a soliton in the form of an infinite current sheet. The method and the results of this paper could be conducive to research studying the emission of strongly magnetized stars. |
Type |
Text |
Language | eng |
Series | University of British Columbia. PHYS 449 |
Date Available | 2014-04-10 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0085971 |
URI | http://hdl.handle.net/2429/46354 |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Peer Review Status | Unreviewed |
Scholarly Level | Undergraduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 52966-Wong_Melody_PHYS_449_2013.pdf [ 810.68kB ]
- Metadata
- JSON: 52966-1.0085971.json
- JSON-LD: 52966-1.0085971-ld.json
- RDF/XML (Pretty): 52966-1.0085971-rdf.xml
- RDF/JSON: 52966-1.0085971-rdf.json
- Turtle: 52966-1.0085971-turtle.txt
- N-Triples: 52966-1.0085971-rdf-ntriples.txt
- Original Record: 52966-1.0085971-source.json
- Full Text
- 52966-1.0085971-fulltext.txt
- Citation
- 52966-1.0085971.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.52966.1-0085971/manifest