Astronomy & Astrophysics manuscript no. paper c ESO 2012 (cid:13) January 24, 2012 Solar Particle Acceleration at Reconnecting 3D Null Points A. Stanier1, P. Browning1, and S. Dalla2 1 Jodrell Bank CentreforAstrophysics,School ofPhysicsand Astronomy,UniversityofManchester, ManchesterM13 9PL, UK e-mail: [email protected] 2 Jeremiah Horrocks Institute,University of Central Lancashire, Preston PR1 2HE, UK Preprint online version: January 24, 2012 2 ABSTRACT 1 0 Context. The strong electric fields associated with magnetic reconnection in solar flares are a plausible mechanism to 2 accelerate populations of high energy, non-thermal particles. One such reconnection scenario, in a fully 3D geometry, occurs at a magnetic null point. Here, global plasma motion can give rise to strong currents in the spine axis or fan n plane. a J Aims. To understandthemechanism of charged particle energy gain in both theexternaldrift region and thediffusion region associated with 3D magnetic reconnection. In doing so we aim to evaluate the efficiency of resistive spine and 3 fan models for particle acceleration, and findpossible observables for each. 2 Methods. We use a full orbit test particle approach to study proton trajectories within electromagnetic fields that are exact solutions to the steady and incompressible magnetohydrodynamic equations. We study the acceleration physics ] R of single particle trajectories and find energy spectra from many particle simulations. The scaling properties of the accelerated particles with respect to field and plasma parameters is investigated. S Results.Forfan reconnection,strongnon-uniformelectricdriftstreamlinescanaccelerate thebulkofthetestparticles. h. Thehighestenergygainisforparticlesthatenterthecurrentsheet,whereanincreasing“guidefield”stabilisesparticles p againstejection.Theenergyisonlylimitedbythetotalelectricpotentialenergydifferenceacrossthefancurrentsheet. - The spine model has both slow external electric drift speed and weak energy gain for particles reaching the current o sheet. r Conclusions.Theelectromagneticfieldsoffanreconnectioncanaccelerateprotonstothehighenergiesobservedinsolar st flares,gainingupto0.1GeVforanomalousvaluesofresistivity.However,thespinemodel,whichgaveaharderenergy a spectrum in theideal case, is not an efficient accelerator after pressure constraints in theresistive model are included. [ Key words. Sun: corona - Sun: flares - Sun: particle emission - Sun: X-rays, gamma rays - Magnetic reconnection - 1 Acceleration of particles v 6 4 1. Introduction electrons is a significant fraction of the emission site den- 8 4 sity, setting tough efficiency constraints on any proposed . ObservationsofHardX-ray(HXR)andγ-rayemissionfrom acceleration mechanism. 1 0 solarflaresbytheRHESSIspacetelescope(Lin et al.2002) It is well accepted that magnetic reconnection plays 2 suggest that a large proportion of magnetic energy is con- the key role in the dissipation of magnetic energy during 1 verted into kinetic energy of non-thermal accelerated par- a flare and there is a growing body of observational sig- : ticles. The dominant HXR sources are chromosphericfoot- natures for the process (see McKenzie 2011). The site of v points of the flaring loops, at which there is continuum reconnection in the standard (CSHKP) flare model (eg. i X free-freeemissionfromabeamofenergeticelectronsincol- Priest 2000) is above the thermal loops, not in disagree- r lision with ambient plasma (Brown 1971). This continuum mentwiththesiteofthenon-thermalcoronalHXRsource. a spectrumgivesbeamelectronenergiesfromaround10keV Super-Dreicer (Dreicer 1959) electric fields associatedwith up to almost 100 MeV (Lin 2006). There is line emission reconnection are one plausible mechanism for particle ac- at the γ-ray end of the spectrum from processes involv- celeration and much theoretical work has been done to in- ing accelerated ions such as neutron-capture and nuclear vestigate the efficiency of this mechanism (for review, see de-excitation(seeeg.Vilmer et al.2011,forareview),sug- Zharkova et al. 2011). gesting ions with energies up to 100 MeV/nucleon. ∼ Earlywork onchargedparticle trajectorieswithin a re- When the emission from the foot-points is weak, or connecting current sheet concentrated on single particle when they are occulted by the solar limb, a weaker HXR motion and energy gain in idealised field configurations. emissionsourceis sometimesobservedabovethe topofthe Speiser (1965) found that chargedparticles in the simplest flareloops(Masuda et al.1994).Recentobservationsoftwo currentsheet,ofoppositelydirectedmagneticfieldandcon- such flares indicate that this emission is non-thermal and stantelectricfield,aretrappedsotheenergygainislimited that the source is actually the acceleration site for a sig- onlybythesheetlength.Withanadditionalsmallandcon- nificantnumber ofenergeticelectrons(Krucker et al.2010; stantmagneticfieldcomponentnormaltothecurrentsheet Ishikawa et al. 2011). The estimated number of energetic plane,the particlesare turned by 90◦ and ejected fromthe 1 A. Stanier et al.: Solar Particle Acceleration at Reconnecting3D Null Points current sheet into the external drift region. Zhu & Parks electric field. For ideal fan reconnection the singular elec- (1993), Litvinenko & Somov (1993) and Litvinenko (1996) tric field occurs in the fan plane, driven by a shearing flow also included a third component of the magnetic field, a across the spine axis. Craig et al. (1995), Craig & Fabling guide field parallel to the electric field. Above a critical (1996) and Craig et al. (1997) found exact solutions to the guide field the trajectory is stabilised against ejection and steady and incompressible resistive MHD equations at 3D the energy gain is once again only bounded by the sheet null points by considering a flux-pileup disturbance field length (Litvinenko 1996). superposed with the background potential magnetic field. This early analytical work was extended with 2D (or The disturbance field induces a current sheet in the spine 2.5D where the fields are invariant in the third dimen- axisandinthefanplanefortheresistivespine andresistive sion) test particle simulations, many of which have been fan modelsrespectively.Craig & Fabling (1998)foundcor- carried out using simple prescribed magnetic and electric responding time-dependent solutions to these steady mod- fields that would be expected in a reconnection solution els,andthe numericalsimulationsofHeerikhuisen & Craig to the MHD equations. These simulations consider the ef- (2004) found reconnectionscalings in agreementwith both fect of the guide field on trajectories and energy spectra steady state and time dependent models at peak recon- within the Harris current sheet (Zharkova & Gordovskyy nection rate. The 3D MHD simulations of Pontin et al. 2004, 2005; Wood & Neukirch 2005) and magnetic X- (2007a,b) also found a hybrid of the spine and fan mod- points (Vekstein & Browning 1997; Hannah & Fletcher els, named spine-fan reconnection, when compressibility is 2006; Hamilton et al. 2005). included.Recentnumericalandanalyticalstudygivesaddi- Heerikhuisen et al. (2002) and Craig & Litvinenko tionalmodelsfornullreconnectionwhenthe globalplasma (2002) used magnetic and electric fields from the exact an- motionisrotationalratherthanashearflow(seeforreview alytical solutions of Craig & Henton (1995) to the 2D in- Priest & Pontin 2009; Pontin et al. 2011). compressible, resistive MHD equations. Also, an approach Itisnotyetknownifthesenullpointsareeffectiveparti- combining numerical MHD simulations with a test parti- cleaccelerators.PreviousworkbyDalla & Browning(2005, cle code has also been used to study 2D forced reconnec- 2006,2008) andBrowning et al.(2010) usedthe idealelec- tion (Gordovskyy et al. 2010a,b). These simulations can tromagnetic fields of Priest & Titov (1996) in a test parti- include compressibility and time evolution, making them cle code, finding the ideal spine reconnection model was morerealisticforcoronalplasma.However,analyticalsolu- effective to accelerate protons and electrons to high en- tionsareessentialtostudyaccelerationduetoreconnection ergies (max 107 eV) for solar coronal parameters. The ∼ in very large Lundquist number plasma at present. ideal fan reconnection model was less effective for protons, The complexity of the coronal magnetic field in a flar- partly as the geometry of the electric drift streamlines was ing Active Region motivates the study of test particle less efficient at delivering particles to regions of high elec- motion in fully 3D reconnection geometries. Reconnection tric field. Guo et al. (2010) used null point magnetic and models in 3D are comparatively new, but it is clear that electric field configurationsfromMHD simulations,finding there are significant qualitative differences from the famil- that strong electric fields due to convective plasma mo- iar 2D models (see eg. Pontin 2011). Reconnection in 3D tion can be efficient at accelerating protons but less so can occur both with and without magnetic null points. for electrons. Litvinenko (2006) used the WKB method of However, the simplest 3D magnetic configuration is based Bulanov & Cap(1988)toshowthatsingleprotonsandelec- on the potential magnetic field about such a null point. tronsclosetothenullinthereconnectingfancurrentsheet Here, the solenoidal condition defines the magnetic topol- can achieve the high energies observed in flares. However, ogy of a 1D spine line and a 2D fan surface (called γ- this energy is limited as the particles become unstable in line and Σ-surface by Lau & Finn 1990) that separates the sheet, due to the potential background field, and are different magnetic flux domains (a linear description of ejected. magnetic configurations at nulls is given by Parnell et al. In this paper we examine test particle trajectories and 1996). Although these null points cannot be measured in energy spectra of protons in electromagnetic fields which the corona at present there is some indirect evidence for are exact solutions to the 3D, steady-state, incompress- their existence. Nulls are common features of magnetic ible and resistive MHD equations at magnetic null points. topology models that reconstruct the magnetic field from These are the resistive spine and resistive fan solutions photospheric magnetogramsinto the corona(see Longcope given in Craig et al. (1995); Craig & Fabling (1996) and 2005, for a review). The application of this method at Craig et al.(1997).Weconsidertrajectoriesthatstartboth several flare sites suggests the importance of these nulls intheouteridealregion,forcomparisonwithparticleaccel- in certain flares (Des Jardins et al. 2009; Aulanier et al. eration results in the ideal models (Browning et al. 2010), 2000; Fletcher et al. 2001). Reconnection at magnetic null andthosethatstartdirectlyinsidetheresistivefancurrent points is also thought to be important for the emer- sheet, to compare with the analytical work of Litvinenko gence of new flux from beneath the photosphere into the (2006). corona (T¨or¨ok et al. 2009; Maclean et al. 2009; Liu et al. The paper is organised as follows. In Section 2 we de- 2011). scribe the model fields used, along with parameters cho- The type of reconnection that occurs at a 3D null de- sen considering the pressure constraints and optimisations pends upon the magnetic configuration and global plasma of Craig et al. (1997) and Craig & Watson(2000). We also flow. Priest & Titov (1996) proposed two models of recon- derive the electric fields and potentials used in the code, nection using a potential magnetic field and prescribed and give approximate external drift velocity scalings. In global flows that satisfy the ideal MHD equations. In ideal Section3wegivethe resultsoftestparticlesimulationsfor spine reconnection a shear flow acrossthe fan plane causes thermaldistributions ofprotonsstartinginthe driftregion frozen-in flux inflow that converges on the spine axis. At for each model. We choose typical high energy particles the spine the field reconnects in the presence of singular from these simulations and follow the single trajectories to 2 A. Stanier et al.: Solar Particle Acceleration at Reconnecting3D Null Points understand the energy gain mechanism. We see how drift proper radial null (Priest & Titov 1996) where the back- timesandenergyspectrascalewiththedifferentparameter groundmagneticfieldlinesinthe fanplanelieintheradial choices in the fan model. In Section 4 we give a summary direction. This backgroundfield is then written as and conclusionson the efficiency of eachmodel for acceler- α ating protons. P = (xxˆ+yyˆ 2zzˆ), (8) 2 − with α giving the sign and strength of the field. For the 2. Model Fields and the Test Particle Code spine model, the displacement field distorts the fan plane We consider two of the reconnection solutions at 3D in the z-direction QS = Z(x,y)zˆ. For the fan model, it null points found by Craig et al. (1995) and developed distorts the spine in the x-direction QF = X(z)xˆ (the in Craig & Fabling (1996) and Craig et al. (1997). The more generalfan case given in Craig et al. (1997) of QF = fields satisfy the resistive,steady-state,andincompressible X(z)xˆ+Y(z)yˆ has not been covered here). MHD equations.These are normalisedin the usual way by the characteristic magnetic field strength B , at a typical 0 2.1. Spine Analytic Fields lengthscaleL fromthenullpoint,andbydensityρ .This 0 0 choice leads to the natural units for velocities in terms of The disturbance field for the spine model in cylindricalco- the external Alfv´en speed, v0 = vAe = B0/√ρ0µ0. The ordinates (r,φ,z) is thermal pressure is normalised by the dynamic pressure at Alfv´en speed, p0 = ρ0vA2e, so that the dimensionless pres- Q =Z(r,φ)zˆ= Bsr M 3,2, r2 sin(φ)zˆ, (9) sureontheL0 boundaryishalftheplasmabetape =βe/2. S rη (cid:18)2 −rη2(cid:19) The dimensionless resistivity, η, is given by (Craig et al. 1997) in terms of the confluent hypergeomet- ηd ric (Kummer) function M(a,b,ζ) (Abramowitz & Stegun η = S−1, (1) L0vAeµ0 ≡ 1972). The flux pile-up factor Bs gives the approximate strength of the magnetic field at a dimensionless distance where ηd is the dimensional resistivity (Spitzer resistivity r from the spine axis, where r is defined as η η in the case of purely collisional plasma), µ is the mag- 0 netic permeability, and S is the Lunquist number which is 4η typiFcoarllycovmerpyleltaerngeessin, tthhee smoalainr cporroopneartSie∼s o1f01t2he−s1o0lu14t.ion rη ≡ 4η¯≡s|α|(1−λ2). (10) p aregivenhere;see Craig & Fabling (1996) and Craig et al. It is the radius of a cylindrical region centred on the spine (1997) for more detail. After normalisation the govern- axis where resistive effects become significant (Craig et al. ing equations consist of the momentum equation, which in 1997). curled form is Theformofthedisplacementfieldinequation(9)isonly (u ∇)ω (ω ∇)u=(B ∇)J (J ∇)B, (2) a solutionto the governingequations providedα<0.This · − · · − · givesfrozen-inplasmainflowalongthefanplaneconverging and the induction equation, on the spine and outflow in the z directions away from ± thenullpoint.Themagneticfieldintheouter(ideal)region (u ∇)B (B ∇)u=η 2B, (3) is also directed inwards along the fan plane and outwards · − · ∇ along the spine axis. Some representative magnetic field with the solenoidal and incompressibility conditions, linesareshowninFigure1(a),thedisplacementtermshears ∇ B =0, ∇ u=0. (4) thefanplaneatφ= π/2towardsthespine axiswhilethe ± · · fieldlines in φ=0,π of the fan plane remain perpendicular Here J is the current density and ω is the vorticity in to the spine. termsofthebulkplasmavelocityu.Inthisnormalisedform To integrate particle trajectories using a test particle they are model we require the electric field. We calculate this from J =∇ B, ω =∇ u. (5) the uncurled form of equation (3) as × × (199T5h),eCthrareige&dimFaenbsliinogna(1l9a9n6a)lyatnicdsCorluaitgioentsaolf.(C1r9a9i7g)ehtaavle. E(r,φ)= η∂Zrˆ+ (1 λ2)PrZ η∂Z φˆ, (11) r ∂φ − − ∂r magnetic and flow fields of the form (cid:20) (cid:21) where P =αr/2 is the radial part of the potential field. B =λP +Q, (6) r This electric field is curl-free (as required for steady- u=P +λQ, (7) state) and so we can calculate the electric potential V to use as a check of energy conservation. This can be found where the scalar 0 λ < 1 gives the shear between the by integrating E = ∇V to get B and u fields. Th≤e vector field P(x,y,z) is a potential − backgroundfieldofstrengthα,andQisadisturbancefield α(1 λ2)r2f(r) V(r,φ)=cosφ − ηrf′(r) , (12) of strength Bs which gives rise to current in the models. 2 − For comparisonwith the particle accelerationresults at (cid:20) (cid:21) 3DnullsinidealMHD(Dalla & Browning2005)wechoose wheref(r)istheradialpartofthedisplacementfieldin(9), the z-axis to be aligned with the spine, with z = 0 as the Z(r,φ)=f(r)sinφ. fan plane. It must be noted that this choice of axis differs We can study the behaviour of these fields at small from that used by Craig et al. (1997). We study only the and large distances using the truncated power series and 3 A. Stanier et al.: Solar Particle Acceleration at Reconnecting3D Null Points asymptotic formulae for the Kummer function respectively (Abramowitz & Stegun 1972). For all our cases the third argumentinthe Kummerfunction isnegative.We havefor 0 ξ 1, zz//LL 00 ≤ ≪ M(a,b, ξ) 1 aξ/b, (13) − ≈ − and for ξ 1, ≫ 1 Γ(b) M(a,b, ξ) ξ−a, (14) 0.5 − ≈ Γ(b a) − 0 in terms of the Gamma function Γ(b). Near the spine axis -0.5 -0.5 r 0, the only contribution to the electric field is from -1 -0.25 ≈ current in the x-direction and 0 x/L E(r r ) ηJ(0)= ηBsxˆ. (15) -0.5 -0.25 0 0.25 0.5 0.25 0 ≪ η ≈ rη y/L0 0.5 The full current distribution is plotted in Figure 1(b), it (a) forms two cylindrical vortex structures that are localised with respect to the resistive region and invariant in the z- direction. At largedistances fromthe spine, the electric fieldgoes as E(r r ) −2ηBssinφφˆ. (16) η ≫ ≈ √π r This has the same functional form as the ideal spine solu- tion ofPriest & Titov (1996), the subject ofprevious work on3Dnull-pointparticleaccelerationbyDalla & Browning (2005,2006)andBrowning et al.(2010).Indeed,bysimple choiceofparameterswecouldsetthemagneticandelectric fieldstoasymptoticallymatchthoseoftheidealcase.Some care must be takenhere as Dalla & Browning(2005) stud- iedpositivenulls,whereα>0,andwithE(0<φ<π)>0. We haveopposite signfor bothelectric andmagneticfields givingthesameelectricdriftinflowquadrantsbutdifferent sign for the convective electric field. Particles that become non-adiabatic in the external region r rη, and gain en- (b) ≫ ergyparalleltheelectricfield,willrotateaboutthespinein the oppositedirectiontothoseinDalla & Browning(2005, Fig.1: a) Representative magnetic field lines for the spine 2006).Inthispaperwewillonlyqualitativelycomparepar- modelwithparametersλ=0.75,B =3.4,α= 2,η =3 s ticle trajectories in the ideal and resistive spine models as 10−3.Thefieldlinesareseededfromthetopand−baseofth×e anasymptoticmatchwillgiverisetounphysicalhydromag- spineaxis.b)Showingthedirectionandrelativestrengthof netic pressures on the edge of the resistive region r rη thecurrentinaplaneofconstantzforthesameparameters. ≈ that were absent in the simplified ideal model (see below). Here, r = √4η¯ 0.12 is the size of the resistive region η ≈ The thermalpressureprofileforthe spine modelcanbe centred on the spine axis. found from integrating the uncurled form of equation (2). It is given in Craig et al. (1997) to be 1 p=p P2+Z2 +λαzZ, (17) boundary (or normalising by suitable solar coronal values, 0 − 2 v = 6.5 106 ms−1, B = 0.01 T, gives E 1/40). Ae 0 wherep isthegaspressu(cid:0)reatthen(cid:1)ullpoint,thefirstterm This value×can be equated with equation (16). C≈rucially, 0 inside the brackets is due to dynamic pressure from the to match the external electric field in the resistive model backgroundflow and the other two terms are from balance to the fixed amplitude electric field in the ideal spine re- with magnetic pressure. All terms except for p0 are nega- connection model requires the scaling Bs ∼η−1 as η is re- tive,asα<0,soconstraintsmustbeputonthevaluesofα duced to suitable solar coronal values (Craig et al. (1997) and B in order to avoid unphysical negative pressures as showed that if we require displacement field at the bound- s discussed in Litvinenko et al. (1996); Litvinenko & Craig ary Z(1,π/2) 1 this also gives Bs η−1). However, this ∼ ∼ (1999);Craig et al.(1997)andCraig & Watson(2000).We scaling gives rise to large magnetic pressure on the sheet give some of the arguments here for the sake of complete- edge. The maximum of the displacement field occurs at ness (see above references for more detail). r rη where Z(rη) Bs giving magnetic pressure there ≈ ∼ The strong electric field (fast electric drift) simulations from equation (17) Z2 B2 η−2. To avoid negative ∼ s ∼ for the ideal spine model studied by Dalla & Browning thermal pressure in the model this requires the null point (2005) were characterised by the dimensional value of the pressure p > (Z(r ))2 η−2 which is unphysically large 0 η ∼ electric field E = 1500 V/m on the r = 1, φ = π/2 for the values of η considered. 0 4 A. Stanier et al.: Solar Particle Acceleration at Reconnecting3D Null Points Craig et al. (1997) showed that B must be limited to s a saturation value on r = r , giving weak electric fields η and small amplitude displacement field on the boundary zz//LL 00 Z(1) 1. Also, at r =1, z 1 we have dynamic pressure due t≪o bulk fluid inflow p ≪ p P2 where P(1) α. ≈ 0 − ∼ 3 We must constrain α B or this dynamic pressure will s ≤ require the gas pressure at the null to be even larger. The 1.5 maximum value we can take for p is the largest possible 0 external hydromagnetic pressure p available to drive e,max 0 the reconnection. We follow Craig et al. (1997) and take pe,max =Be2,max/2 where the maximum external magnetic -1.5 5 field is that of a sunspot at the photosphere, B =0.3 e,max 2.5 T. This gives a normalised saturation value Bs,max =30. -3-5 0 So far we do not know the value α should take, but -2.5 0 -2.5 y/L0 expect that the bulk fluid exhaust from the reconnection x/L 2.5 regionisoftheorderofthelocalAlfv´enspeed.Theexhaust 0 5-5 on the edge of the current sheet at a global distance from the null, r =r , φ=π/2, z =1, is given by Fig.2: Representative magnetic field lines for the fan solu- η tion with parameters λ = 0.75, B = α = 10, η = 10−6. v(r ,π,1) λB α, s | η 2 |≈ s− The lines are again seeded from the top (solid lines) and base (dashed lines) of the spine. where the local Alfv´en speed is v (r ,π,1) = B(r ,π,1) B λα | A η 2 | | η 2 |≈ s− regioncentredonthefanplane,z =0.Thisformofsolution for our choice of normalisation. As we are not interested is only valid for α > 0 which gives a positive null point, in the case where λ = 1 (where there is no shear between the field is washed in from the global boundaries at z = the velocity and magnetic fields) we have α B for 1 and it exits the simulation box radially along the fan Alfv´enic exhaust. This is the largest magnitud≈e −of αs we ±plane. Some representative magnetic field lines are shown cantakewithouthavingproblemsduetodynamicpressure. in Figure 2; the displacement field shears the spine axis as It also leads to the thinnest current sheet and thus max- it approaches the fan plane giving rise to strong current imises the current density in the resistive region. However, inside the resistive region. as Craig & Watson (2000) show, the dissipation rate is The electric field is E =yˆ ηX′(z) (1 λ2)P X(z) +zˆ (1 λ2)P X(z) , z y Wη =η J2dV πηBs, (18) − − − (22) ≈ (cid:2) (cid:3) (cid:2) (cid:3) Z and the electric potential is which has no α dependence as the increase in current den- V(y,z)= αy(1 λ2)[η¯X′(z)+zX(z)], (23) sity due to resistive region thinning is cancelled by the rη2 − − dependence of the total dissipation volume. whereJ(z =0)=X′(0)yˆiscurrentdensityatthecentreof Theelectricdriftvelocityintheexternalregionisgiven the sheet. The current only has z-dependence; it is infinite by in extent inthe x and y directions.This is clearly unrealis- tic, although resistive MHD simulations by Pontin et al. ηB sinφ 2zrˆ rzˆ v (r r ) s − − [v ] (19) (2007b) find that spine-fan reconnecting current sheets E ≫ η ≈ λα√π r(r2/4+z2) Ae formed due to shear flows around a null point spread out | | (cid:18) (cid:19) along the fan plane in the incompressible limit. Note that which is very slow when α = B . It is thus necessary to | | s analyticmultiplenullsolutionsfoundbyCraig et al.(1999) limit the magnitude of α so that results can be obtained havefinitecurrentsheets,avoidingthisproblem.Inoursim- with reasonable integration times. For the simulations in ulationsbelowweconsiderparticleaccelerationonlywithin Section 3 we use B = 10, α = 0.1, this limits the re- connection exhaust cslose to the sp−ine current sheet to sub- arestrictedrangeof5L0,effectivelylimitingthesizeofthe current sheet. Alfv´enic speeds. The thermal pressure profile for the fan model is (Craig et al. 1997), 2.2. Fan Analytic Fields p=p (P2+X2)/2 αλxX/2. (24) 0 − − The displacement field for the fan model is However, in this case a displacement field of order unity B z 3 3 z2 on the z = 1 boundary, X(1) 1, gives the scaling QF =X(z)xˆ= η¯1s/2M(cid:18)4,2,−2η¯ (cid:19)xˆ, (20) Bdrso∼maηg−n1e/t4ic(Cprraesigsuerteaoln.1t9h9e7)c.uTrrhei∼nstgsivheesetmeudcghewceoamkeprahreyd- (Craig et al. 1997). We define z as to the spine model but it is still too large for the values of η η considered. Again we saturate B =30 and we have s,max 2η α B to avoid problems from dynamic pressure. s z 2η¯ , (21) ≤ η ≡ ≡sα(1 λ2) Craig & Watson (2000) show that the Ohmic dissipa- − tion rate from the fan model is p the approximate height at which X takes the maximum W =η J2dV ηB2/z (25) value,Xmax Bs.Itisameasureoftheheightofaresistive η ∼ s η ≈ Z 5 A. Stanier et al.: Solar Particle Acceleration at Reconnecting3D Null Points and so in this case, for fixed (saturated) B , the maximum s dissipation occurs with the thinnest current sheet (the so called optimised solution). The thinnest sheet we can have subjectto the dynamic pressureconstraintis whenα=B s (for any fixed value of λ). Also, as this choice gives the largest current density, it maximises the resistive electric field within the sheet which is interesting for particle ac- celeration. As above, this choice of α sets the bulk fluid exhaust at x2+y2 =1, z =z to the local Alfv´en speed. η Using the asymptotic approximation (14) we find that the z-component of the electric drift that brings the parti- cles to the fan plane is, for x, z η1/2, ≫ (1 λ2)P P X(z) v − x z B3/4η1/4, (26) Ez ≈ λP2 ∼ s for the optimised solution α=B . This gives electric drift s inflowforpositivex,z (asP <0),andoutflowforpositive z z and negative x. It is much faster than the spine case due Fig.3: Angular distribution of protons from the null point to the more favourable scaling with resistivity. There are for spine model at t = 1.6 106T , at which the en- also fast drift streamlines in the x-y plane that that can ω,p × ergy spectrum is steady state, for parameters λ = 0.75, be found from the numerical (or approximate analytical) B = 10, α = 0.1, η = 10−6. The initial distribution was solution of s − dx dy Maxwellian at T =86 eV in the upper right inflow region. = , (27) v v Ex Ey we numerically plot these streamlines on top of the single All magnetic fields mentioned in Section 2 have dimen- particle trajectory results. sions of B = 0.01 T, typical for the solar corona. We set 0 v = v = 6.5 106 ms−1 (corresponding to a number 0 Ae 2.3. Test Particle Code and Parameter Choice density n = 1.1×26 1015 m−3). All dimensionless times 0 × quoted are in terms of the gyro-period, T =2πm/(qB ). ω 0 We modify the test particle code of Dalla & Browning To examine single particle trajectories in both models (2005, 2006, 2008) and Browning et al. (2010) to use the we choose values of B = 10, λ = 0.75, η = 10−6. This η electromagnetic fields given above (from the solutions of s value is rather large, towards the highest possible anoma- Craig et al.1997).AVariable-StepVariable-OrderAdam’s lous resistivities (with L = 104 m), but is useful to ob- method, where the step size is recalculated to properly 0 serveparticlesenteringthecurrentsheet.Wevaryallthree resolve gyro-motion, is used to integrate the relativistic of these parameters to produce scalings of energy gain and Lorentz equation. drift times, where we use values as low as those expected dp q p in purely collisional plasma i.e. η =10−12. = E+ B , (28) dt m γm × (cid:18) (cid:19) where p is the momentum of the particle, γ is the Lorentz 3. Results factor,q andm arethe chargeandrestmassandE andB 3.1. Spine Global Trajectories are the analytic expressions for the electric and magnetic fields for each model. Initially, we place a distribution of 5000 protons with We use the expressionsfor the electric potential V, cal- Maxwellian velocities of temperature T = 106 K (86 culated in equations (12) and (23) to calculate the electric eV) in the spine model fields. The protons have positions potentialenergyateachtimestep.Withthisweverifythat from a uniform random distribution at a global distance the total energy x2+y2+z2 =1 from the null point. We only discuss here W =ǫ +qV (29) k protonsthat startin the upper rightinflow regionoflongi- isconservedwhereǫ =(γ 1)mc2.Foreachsimulationwe tude 0<φ<180◦ and latitude 0<β <90◦ (here φ=0 is k − find that this is conservedup to 5 significant figures. Also, the x-axis and β =0 is the fan plane). to check the code handles non-adiabatic motion in strong Figure 3 shows the final spatial distribution and ener- magneticfieldgradientsofacurrentsheetwereproducethe gies of the particles at time t = 1.6 106T 10 s, at ω,p × ≈ resultsofSpeiser(1965),includingtheejectiontimeforthe which the energy spectrum in Figure 6 becomes steady- case with backgroundfield. state.Thoseprotonsstartinginthe lowerleftinflow region We choose L , the normalisinglength scale, to be L = have final distributions as in Figure 3 after reflections in 0 0 104mforglobalsimulationstokeepintegrationtimesshort. bothφ=0andβ =0apartfromstatisticaldifferences.The This size of simulationbox canbe consideredas a localre- parameters used here are B = 10, η = 10−6, α = 0.1. s − gion around the null at which the linear background field This value of α limits the bulk flow exhaust speed to be and flow in equation (8) are good approximations. We use sub-Alfv´enic but it increases the electric drift speed in the alargervalueofL =106mforsimulationswhereparticles externalregion(see equation(19)) due to weakermagnetic 0 are initially distributed within the current sheet, as veloc- field on r =1. This gives reasonable simulation times, but ities are typically much faster here. Note that a change in therearestillsomeparticlesintheupperrightinflowquad- L also changes the value of η as given in equation (1). rant at the end of the simulation. 0 6 A. Stanier et al.: Solar Particle Acceleration at Reconnecting3D Null Points Fig.4:Typicalprotonglobaltrajectoryfrompopulation’A’ for parameters λ = 0.75, B = 10, α = 0.1, η = 10−6. s − The particle is taken from the many particle simulation having initial position x = ( 0.52,0.80,0.29) and velocity − v =( 0.0044,0.0013, 0.0088)v .Themagneticfieldlines Ae − − (thindashed)areaprojectionofthefieldfromtheplaneof Fig.5: Typical proton trajectory from population ’B’ the trajectory φ 120◦. Inset shows the 3D trajectory of in the many particle simulation, with initial posi- ≈ the proton as it crosses the spine-axis, the solid line in the tion x = ( 0.54,0.78,0.31)L0 and velocity v = − centre is the line B (x,y,0.95)=0. ( 0.004, 0.006,0.002)v .Thedashedlinesshowthe pro- z − − Ae jection of the magnetic field from the plane of motion, φ 110◦, onto the y-z plane. Inset shows the motion in ≈ the x-y plane close to the spine axis. The purple arrows There are two main populations of accelerated parti- show the direction and relative magnitude of the gradient cles. The population labelled ’A’ in Figure 3 is close to the driftvelocityandthedash-dottedlinesshowcontoursofthe fan plane, β . 10◦, with energy ǫ & 1 keV, and with electricpotential,withtheintersectingtickmarkindicating k longitude |90|◦ . φ . 90◦ comprising of about 8% of the lower potential to the right. − totalprotonnumber.The maximumparticleenergyofthis population is about 15 keV. Note that the current in the spine axis is aligned with φ = 0 through the centre of this A typical proton trajectory from population ’A’ is population. There are also some high energy protons scat- showninFigure4.Theprotonwhichstartsat(x ,y ,z )= 0 0 0 tered at large positive latitudes for φ . 0, and at large ( 0.52,0.80,0.29)in the upper right hand inflow quadrant negative latitudes for φ & 0. To look more closely at what in−itially moves away from the null but mirror bounces and is happening herewe willchoosea typicalprotonfromthis travelsbacktowardsthespinealongthefanplane.Theelec- population and follow its trajectory below. tricdriftspeedincreasestowardsthespinecausingthepro- For those particles that have crossed the fan plane, tontoentertheresistiveregion,whichhasradiusr 0.01, η ≈ β = 0, into the lower right outflow quadrant, the spatial aboutthe spine axis.Itentersat(x,y,z) ( 0.01,0,0.95) andenergydistributionlookssimilartotheidealspinecase after t = 3 105T 2 s (inset). A≈t t−his point the ω,p × ≈ in Dalla & Browning (2006). The accelerated population proton becomes unmagnetised as the gyro-radius becomes which has ǫ & 1 keV and β . 85◦ is labelled ’B’. This comparable to the length-scale of magnetic field gradient k − population is about 6% of the total protons in the simula- ρ/L > 1 (we typically find that gyro-motion starts to ∇B tion and the maximum kinetic energy in this population is break down when ρ/L &10−2). ∇B ǫ 12 keV. The angular distribution differs slightly Theprotonisthendirectlyacceleratedinthex-direction k,max ≈ with the ideal case in that there are few particles found parallelto the currentatthe spine with du/dt qE /mas 0 between the latitudes 70◦ < β < 85◦; particles appear it crosses r =0. For small displacements in the≈y-direction − − to be closer to the negative spine axis in the resistive case. astrongLorentzforceduetotheB fieldreturnsittoy =0 z 7 A. Stanier et al.: Solar Particle Acceleration at Reconnecting3D Null Points line. These oscillations are Speiser-like (Speiser 1965) with frequency approximately ω t1/2. ∝ The Speiser-like motion finishes and the first gyrations start(notshown)whentheprotonreachesr 5r atwhich η ≈ ρ/L .1. However,the energygain of ǫ 11 keV is lo- ∇B k ≈ calisedtowithin x 2r ,duringwhichthe trajectorydoes η ≈ notdeviatemuchfromthex-direction(notethey-axisscale intheinsetofFigure4).Ineffect,theprotonhasleftthelo- calised current sheet while unmagnetised but before it can be ejectedbythe backgroundfieldcomponents,incontrast to2Dcurrentsheetconfigurationswithweakguidefield(eg. Speiser 1965; Litvinenko 1996). Figure 4 may give the im- pression that the particle is being ejected, however, this is just the centre of the Speiser-like oscillations following the B (x,y,z = const.) = 0 line (which here is not straight as z in the usual 2D configurations). This behaviour is evident considering the F = qv B force for the unmagnetised × proton if B is the dominant component of the magnetic z field. Fig.6:Energyspectrum fromthe many particle simulation for protons in the spine model, with parameters λ = 0.75, it hAasftewreathkeelpercottroicndbreifcto,mves re-mva.gnItetfioselldowast trhe≈fie5lrdη- Bs = 10, α = −0.1, η = 10−6. For particles leaving the E ≪ ω R=5 sphere, the energy at the time of crossing is used. lines closely and mirror bounces travelling back towards the spine: there the proton is takenup to high latitude be- fore it bounces again. This mirror bouncing is the reason for the ’scattered’ acceleratedprotons in Figure 3, some of 3.2. Fan Global Trajectories which are at large latitudes. The many particle simulation for the fan model is shown Atypicalparticletrajectorychosenfrompopulation’B’ in Figure 7 for the optimised solution B = α = 10, s isshowninFigure5.Theprotonstartsat( 0.54,0.78,0.31) with η = 10−6, λ = 0.75. The initial distribution has − and drifts towards the spine but bounces and crosses the thermal energy ǫ = 86 eV with uniform random posi- k fan plane instead. It exits the simulation box down the tion in the upper inflow quadrant 90◦ < φ < 90◦ and base of the spine axis, reaching an energy ǫk = 6.72 keV 0 < β < 90◦. The final angular dist−ribution is taken from as it crosses z = 5. As there is no electric field in the when the proton distribution reaches a steady state in en- − z-direction, the energy gain must occur due to motion in ergy at t = 4000T 0.025 s. This is more than two ω,p the x-y plane, which is also shown in Figure 5. The proton ≈ orders of magnitude faster than the spine model for simi- enterstheregionclosetothespineaxisparalleltoacontour lar parameters (even after the spine drift was increased by of the electric potential, but then drifts across the contour limitingα<B )astheexternalelectricdrift(equation26) s duetostronggradientdrift.Whiletheprotongainsenergy, scales more favourably with the resistivity. Protons that thegradientdriftislargerthantheelectricdriftbyafactor cross the R = 5L spherical boundary from the null point 0 of 2 with the latter directed inwards towards the current before this time are stopped and the energy and angular sheet. The proton is stopped as it reaches z = 5L which 0 position at time of crossing is used. The t=0 angular dis- − we do consistently throughout these simulations. At the tribution has some structure in terms of final energy gain time of stopping it is actually losing energy as it re-crosses as the initial random thermal velocities are dominated by the same electric potential contours. However, some other the strong electric drift. protonsfromthemany-particlesimulationinFigure3reach Within this structure there is some asymmetry in φ. the current sheet at low latitudes, gaining higher energy. Indeed,wedonotexpectsymmetrybetweenparticlesdrift- The energy spectrum for the spine simulation is shown ing clockwise and anti-clockwise about the null as in the in Figure 6. If protons cross the R=5L spherical bound- ideal case (Dalla & Browning 2006) now that there is a 0 ary we use the energy at the instant of crossing (if this current in the y direction. Those protons with ǫk 107 is not done some protons reach order 102L which be- eV at φ 20◦ (the yellow vertical band to the ≈left of comes unrealistic as the background fi∼eld incr0eases with- the green≈ve−rtical band in Figure 7(a)) do not enter the out bound away from the null, also causing the time-step current sheet, but gain high energy, as they are unmag- to decrease and simulation time to increase). The initial netised slightly, ρ/L∇B 10−3, due to very fast electric ∼ Maxwellian spectrum hardens to a broken power law with drifts close to the sheet. Here the first adiabatic invariant, maximum energy of about ǫ 15 keV. This maximum the constancy of µ, is also violated. k ≈ energy canbe understood as the difference in potential en- Typically, the high energy protons of Figure 7 start ei- ergy across the spine current sheet, ǫ qEx where ther close to the x-axis at low to mid latitudes (about 7% k acc E E ηB /r [v B ]andx 2r ∼[L ]isthe accel- of the total number at latitude β & 1◦ with final energy 0 s η Ae 0 acc η 0 ≈ ≈ ≈ eration distance (from r . x . r ), as the electric field ǫ & 10 MeV), or they start at very low latitude close η η k,fin dropsoffquicklyfor x −>r .Forthe parametersused,this to the fan plane (< 1% of total at β . 1◦ and ǫ & 10 η k,fin | | gives ǫ 13 keV. This approximateexpressionhas no de- MeV). At t = 4000T these energetic protons are found k ω,p pendence≈upon the parameter α, so the limiting of α<B at β 0 either side of φ = 90◦; the y-axis in the direction s ≈ should not have a large effect on this result. of the fan current. 8 A. Stanier et al.: Solar Particle Acceleration at Reconnecting3D Null Points 1 2 (a) (a) 2 1 (b) Fig.7: a) Angular distributions of protons in fan model at t = 0, with initial temperature T = 86 eV. The x-axis is φ =0 and the fan plane is β = 0. Protons are coloured by (b) thefinalenergyatt=4000.Parametersusedareλ=0.75, Fig.8: Two typical proton trajectories from the many par- B = α = 10, η = 10−6. b) Angular distributions at time s ticlesimulationsinthefanmodel.Proton’1’isrepresented t=4000T whentheenergyspectrumhasreachedsteady ω,p byathinlinewithinitialposition(0.86,0.41,0.30),andpro- state. ton ’2’ by a thick line with initial position (0.8,0.003,0.6). The parameters are λ=0.75, B =α=10, η =10−6. The s solid lines are representative magnetic field lines (seeded At t = 4000Tω,p there are a small number of high en- from the top of the spine axis and projected into the 2D ergy protons scattered at high latitudes (about 0.1% with planes)andthearrowsshowthedirectionandrelativemag- ǫk > 10 MeV). These enter the current sheet temporarily nitude of the electric drift velocity. a) In the x-y plane, atnegativelongitudefarfromthenullpoint,butexitagain wheretheelectricdriftarrowsarefromtheedgeofthecur- without any Speiser-like motion. They become slightly un- rent sheet z = z . b) In the x-z plane close to the current magnetised,withmaximumρ/L∇B 10−2,followingcom- sheet, where theηelectric drift arrows are plotted on y =0. ≈ plicated trajectories. As they are not typical they are not The initial positions are not shown in this plane. investigated further in the external region, but the be- haviourwithin the currentsheet is discussed below (shown in Figure 12). Figure 8 shows the trajectory of two typical protons takenfromthesimulation.Proton’1’startsat(x ,y ,z )= dominant as the particle exits the simulation box paral- 0 0 0 (0.86,0.41,0.30)anddriftsaroundthenullpointduetothe lel to the negative x-axis. The first adiabatic invariant is strong azimuthal electric drift. Although it drifts down to- not violated, µ=const. and the maximum ρ/L 10−4 ∇B ∼ wards the current sheet, it reaches a minimum height of at closest point of approach to the sheet. The proton is z 15z before it flows into the outflow quadrant, not strongly magnetised throughout. Despite not reaching the η ≈ entering the sheet. The main velocity contribution is elec- current sheet the energy gain is still considerable,reaching tricdriftasitmovesaroundthenullpoint,butv becomes 0.5 MeV as it crosses the R=5 sphere. k 9 A. Stanier et al.: Solar Particle Acceleration at Reconnecting3D Null Points a spherical surface of radius R = 10 from the null point, instead of R = 5 that has been used consistently through- out these simulations. Now the ’flat tail’ at 107.5−8 eV in Figure 9 becomes a ’bump on tail’ centred at 108 eV (not shown) disconnected from the main distribution. Here, the powerlawpartremainsmostlyunchanged.Thepopulation of protons that is trapped in the sheet as it crosses R = 5 due to the strong ’guide field’ remains trapped at R = 10 where B (y) has doubled in strength. y 3.3. Fan Current Sheet Trajectories The simulations considered thus far concern proton tra- jectories starting from the external region, at a distance R = 1 from the null point. However, most of the protons entering the current sheet do so far from the null point. In the following, protons are initially distributed within the fan current sheet close to the null, to study the transition from non-adiabatic to adiabatic motion. Fig.9: Energy spectrum for the many particle fan simula- Firstly, we place particles within the sheet so that they tion for protons with parameters λ = 0.75, B = α = 10, s η =10−6.ForparticlesleavingtheR=5sphere,theenergy are initially unmagnetised by the By(y) component of the background field. They are magnetised only by the strong at the time of crossing is used. B (x,z).Theprotonsareuniformlydistributedinthearea x x < 1; y = 0; z < z with initial thermal energy η | | | | T =86eV.Figure10showsthepositionof2000protonsat Particle 2 starts at (0.8,0.003,0.6)L . The azimuthal 0 t=2500T (a), and t=17500T (b), during this sim- electricdriftisweakclosetothex-axisandtheprotondrifts ω,p ω,p ulationfortheparametersλ=0.75,B =α=5,η =10−8. downtothefancurrentsheet.Itentersthesheetat(x,y)= s We increase the dimensional box length to L = 106 m as ( 0.02,0.15) and becomes unmagnetised: ρ/L > 1 and 0 − ∇B velocities inthe currentsheetaretypically fast,giving rea- µ is not conserved. We observe Speiser-like oscillations as sonableintegrationtimes.Thismakesourresultsmorecom- theprotonisacceleratedinthey-direction.Att=0.846ms parabletotheapproximateanalyticsolutionsofLitvinenko after enteringthe sheet,itpassesoutofthe simulationbox (2006). Note that η decreases due to the increase in L in at R = 5. Here, the particle is still within the sheet with 0 equation(1).Weagainartificiallystoptheparticlesasthey v =0.36c andǫ =67MeV. Using this time period in the k k cross the R=5 spherical surface. directaccelerationformula,y =qE t2/2m,withtheelectric 0 At t=2500T most of the protons are strongly mag- ωp field on z = 0, E0 = ηBs/η¯1/2 [vAeB0] from (22), gives netised by the Bx(x,z) magnetic field. Inside the current y 5.Thustheprotonisdirectlyacceleratedinthecurrent sheet, z <z wecanuseequation(13)togetapproximate sh≈eet for the entire length of the simulation box. However, express|io|ns foηr the electric and magnetic fields, thismotionisnotSpeiser-likethroughoutasρ/L <10−2 ∇B when the proton reaches y = 1.5L . The proton reaches a E E yˆ ηB /z yˆ [v B ], (30) 0 y s η Ae 0 ≈ ≈ global distance in the y-direction and becomes magnetised by the background B component of the magnetic field, λαx B z λαy y B + s , , λαz [B ], (31) which acts as a kind of guide field. When the simulation ≈ 2 z 2 − 0 is run without stopping the particle at R = 5, the proton (cid:18) η (cid:19) is not ejected from the currentsheet throughoutthe whole Ez is small except at global distance in y (see below). simulation time t=4000Tω,p. For a proton starting at x = 0, y = 0, z = zη, on the This particle enters the current sheet at a distance edge of the current sheet, the background components of R 0.15 from the null point; however, this distance is the magneticfieldarenegligible.Theprotondriftstowards not≈typical for the many particle simulation in Figure 7. the vertical centre of the sheet vEz (η/z)zˆ [vAe]. It ≈ − In the simulation 9.3% of the total particles reach the cur- becomes unmagnetisedatthe fanplane,z 0,closeto the ≈ rentsheet,afterameantimeofabout800T .Theaverage nullpointandisdirectlyacceleratedinthe y-direction.We ω,p distance from the null point of particles entering the sheet comparethistrajectorytotheanalyticalWKBsolutionsof is R 2.2; some remain magnetised by the background Litvinenko (2006). The ejection time for a non-relativistic magne≈tic field inside the sheet. protonthatisunmagnetisedclosethenullpoint,x 0, z ≈ ≈ The energyspectrum for the fansimulationis shownin 0, in the fan current sheet in our parameters is Figure9.Almostalltheprotonsareacceleratedintoanon- thermal power law distribution f(E) E−γ, with slope m2BsL0 1/3 γ 1.5. For most particles, this effic∝ient acceleration is tejec ≈ q2z B λ2α2E , (32) ≈ (cid:18) η 0 y(cid:19) due to the fast electric drift speed in the fan model being muchlargerthantheinitialthermalvelocity.Thespectrum (Litvinenko2006)providedthatthe protonremainswithin appears to have reached a steady state by t = 4000T ; the non-adiabatic region and the displacement magnetic ω,p however, it also depends upon the position at which pro- field gradient is much stronger than the gradient from the tons are stopped as they leave the simulation box. As a backgroundcomponent, B /z λα. The secondassump- s η ≫ test we repeat the simulation but stopping the protons at tionisvalidfor oursimulation;however,wedonotobserve 10