3D simulations of realistic power halos in magneto-hydrostatic sunspot atmospheres: linking theory and observation
Carlos Rijs, S.P. Rajaguru, Damien Przybylski, Hamed Moradi, Paul S. Cally, Sergiy Shelyag
DD RAFT VERSION O CTOBER
13, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
3D SIMULATIONS OF REALISTIC POWER HALOS IN MAGNETO-HYDROSTATIC SUNSPOT ATMOSPHERES:LINKING THEORY AND OBSERVATION C ARLOS R IJS ,S.P. R
AJAGURU ,D AMIEN P RZYBYLSKI ,H AMED M ORADI ,P AUL
S. C
ALLY ,S ERGIY S HELYAG (Dated: October 13, 2018)
Draft version October 13, 2018
ABSTRACTThe well-observed acoustic halo is an enhancement in time-averaged Doppler velocity and intensity powerwith respect to quiet-sun values which is prominent for weak and highly inclined field around the penumbra ofsunspots and active regions. We perform 3D linear wave modelling with realistic distributed acoustic sources ina MHS sunspot atmosphere and compare the resultant simulation enhancements with multi-height SDO obser-vations of the phenomenon. We find that simulated halos are in good qualitative agreement with observations.We also provide further proof that the underlying process responsible for the halo is the refraction and return offast magnetic waves which have undergone mode conversion at the critical a = c atmospheric layer. In addition,we also find strong evidence that fast-Alfvén mode conversion plays a significant role in the structure of thehalo, taking energy away from photospheric and chromospheric heights in the form of field-aligned Alfvénwaves. This conversion process may explain the observed "dual-ring" halo structure at higher ( > Subject headings: sun: magnetic fields – sun: oscillations – sun: photosphere – sun: chromosphere – sunspots INTRODUCTION
A complete picture of the interaction between wave mo-tions and magnetic field in the solar photosphere and chromo-sphere is not yet available to solar phycisists.Significant uncertainties still exist in the computation of he-lioseismological inversions in active regions for instance, es-pecially given that the atmosphere above photospheric levelsundoubtedly plays a role in muddying the seismic observablesat the surface (Cally & Moradi 2013).The theory of mode conversion provides a framework as tohow active regions act as a gateway between the subsurfaceand the overlying atmosphere and modify the properties ofotherwise normal acoustic p - modes.The first and most important property of active regions to beexplained in a mode conversion context was the well knownabsorption of p - modes (Braun et al. 1987). Upon initialsuggestion by Spruit & Bogdan (1992), it was eventually de-termined that both conversion to the field-aligned slow mode(which travels downwards into the interior) and to the up-wards travelling acoustic mode (for non-trapped waves) werethe responsible mechanisms (Cally & Bogdan 1997; Callyet al. 2003).The acoustic halo was first noted in Dopplergrams alongsidethe aforementioned p - mode absorption as a peculiar en-hancement in 6 mHz power (with respect to average quiet sunvalues) which extended several Mm radially outwards fromthe umbra (Braun et al. 1992; Brown et al. 1992; Toner &Labonte 1993).Later, in studies utilising the Michelson Doppler Imager (MDI) (Scherrer et al. 1995) onboard the
Solar and Helio- [email protected] School of Mathematics, Monash University, Clayton, Victoria 3800,Australia Monash Centre for Astrophysics (MoCA), Monash University, Clay-ton, Victoria 3800, Australia Indian Institute of Astrophysics, Koramangala II Block, Bangalore560034, India Trinity College, The University of Melbourne, Royal Parade,Parkville, Victoria 3052, Australia spheric Observatory (SOHO), it was noted that the enhance-ment was not present in measurements of the continuum in-tensity (Hindman & Brown 1998; Jain & Haber 2002).This suggests that either there is a process at work affectingobserved power somewhere in the height range between theintensity continuum height and the Doppler velocity observa-tion height, or that the mechanism causing the enhancementis not a process that is detectable in measurements of intensi-ties.It turns out that the former case is much more likely, as in-tensity halos taken from spectral lines at greater heights havesince been observed and studied in detail (Moretti et al. 2007;Rajaguru et al. 2013).Also using MDI, Schunker & Braun (2011) examined 7 daysof observations of the active region AR 9787 and showed thathalos are manifested for relatively horizontally aligned, weak-to-moderate magnetic fields (150 G < | B | < k ) for a given frequency ( ν ) (compared tothe ridges from an area of the quiet sun) and that this effect ismore pronounced for larger k , which sugggests that shallowerwaves are being more strongly affected in the enhancementregion.The most comprehensive observational halo study to date, byRajaguru et al. (2013) utilised the Helioseismic and MagneticImager (HMI) and
Atmospheric Imaging Assembly (AIA) in-struments onboard the
Solar Dynamics Observatory .The authors conducted a multi height analysis of several ac-tive regions, measuring the time-averaged power from intensi-ties and velocities corresponding to 6 different heights. Fromthe intensity continuum at z = 0 (the base of the photosphere,where the optical depth is unity) to Doppler velocities of theFe I 6173.34 Å line at around z = 140 km to intensities mea-sured from the AIA 1600 Å and 1700 Å chromospheric spec-tral lines, halo properties were compared and analysed in de-tail. The findings can be summarised as follows:1. The halo is present for non-trapped frequencies, be- a r X i v : . [ a s t r o - ph . S R ] D ec ginning at 5.5 - 6 mHz (as observed by all referencesabove) and is present up to at least 9-10 mHz. The 6mHz halo is the strongest in measurements of the Fe I6173.34 Å Doppler velocity at z = 140 km.2. The halo magnitude is a clear function of height.There is no enhancement in the time-averaged intensitycontinuum ( z = 0) power or in the derived line-wingDoppler velocity ( z = 20 km). For weak-field regions atthese heights, there is also a uniform wave power abovethe acoustic cutoff, which is to be expected. However,at z = 140 km (the aforementioned HMI Dopplervelocity line) the situation is markedly different, andthe halo comes into effect.3. The halo is clearly present in observations of the chro-mosphere, as measured by AIA. The time-averagedpower of the 1600 Å and 1700 Å wavelength channels(corresponding approximately to z = 430 and 360 kmrespectively) shows a halo in the 7-10 mHz range, thatspreads radially with height, agreeing with the sugges-tions of Finsterle et al. (2004).4. The spatial extent and structure of the halo changesabove about 8 mHz. This higher frequency halo isseen in power maps to be thinner and more confinedspatially than the more diffuse structure seen at 6 mHz.Radially outwards from this higher ν field is a region ofslightly reduced power, which in turn is surrounded bya diffuse, weak halo region, extending radially manyMm into quiet regions (Rajaguru et al. 2013).In this study we are interested in providing a consistenttheoretical explanation for the acoustic halo. There are ofcourse a variety of existing theories as to the mechanismbehind the phenomenon.By conducting radiative simulations for instance, Jacoutotet al. (2008) determined that strong magnetic fields can alterthe scale size of granulation cells, which in turn can modifythe local excitation frequency of resultant photosphericwaves. They found that the stronger field also increases theamplitude of non-trapped waves at frequencies consistentwith halos.Kuridze et al. (2008) show semi-analytically that waveswith m > m is the azimuthal wave number) canbecome trapped under field free canopy regions, resulting inan enhancement of higher frequency wave power.Hanasoge (2009) suggests that the halo is a consequence ofthe equilibrium state of the solar surface, and that the localoscillation can be shifted to a lower mode mass (Bogdan et al.1996) due to scattering from the magnetic flux tube.We will discuss why these theories do not appear viable inlight of our simulation results in the discussion at the end ofthis paper. Mode conversion
In Rijs et al. (2015), we performed 3D simulations todetermine whether there was promise in the suggestion ofKhomenko & Collados (2009) that it is in fact the overlying atmosphere that is directly responsible for the halo. Specif-ically that the addition of energy from high frequency non-trapped waves which have travelled above the Alfvén-acoustic equipartition ( a = c ) layer and undergone mode conversionand refraction are responsible. In this case, the process ofmode conversion describes the intrinsic physics.At greater depths below the solar photosphere, the plasma β (where β = P g / P m , with P g and P m being the gas and magneticpressures respectively) increases. Several Mm below the sur-face the plasma is dominated by hydrodynamic physics andwaves are governed by the standard gas sound speed ( c ).Conversely (assuming one is in the proximity of an activeregion of some sort), well above the surface the gas density( ρ ) drops, the β becomes small and waves are governed morestrongly by the Alfvén speed ( a ), where a ∝ | B | / √ ρ , and | B |is the local magnetic field strength.There is therefore a layer of the atmosphere (roughly where β = 1) where a and c equate - the so called a = c layer. Atthis height, the phase speeds of the magnetoacoustic fast andslow waves are equal, allowing the two modes to interact. En-ergy can be channeled from the fast to the slow branch or viceversa.The fast wave is largely acoustic in nature when a < c andmagnetic when a > c , and it is this fast magnetoacoustic wavethat will refract and then reflect at the fast wave turning height(where ω = a k h , with ω = 2 πν and k h being the horizontalcomponent of the wavenumber, k ), returning downwards fromabove the a = c layer.Energy is preferentially converted from the fast-acousticmode to the fast-magnetic mode if there is a large attack anglebetween the wavevector of the incident wave and the orienta-tion of the magnetic field. If the attack angle is small thenenergy will be primarily channeled into the field aligned slowmode. This perhaps explains why halos are observed amongsthorizontal field; The line of sight component of the Dopplervelocity is largely vertical (when observing at disk center) andprovides a large attack angle with the horizontal field.The theory can also explain the spreading of the halo that isobserved with height (Rajaguru et al. 2013), given that the a = c layer is located at greater radial distances from the um-bra as a function of height.Waves with frequencies below the acoustic cut-off are gen-erally unable to reach the a = c height, as they have re-flected inwards, which is presumably why halos are only ob-served at non-trapped frequencies. The fast magneto-acousticwave provides the excess energy in observable regions, whichwould otherwise not be present in the quiet sun (see Cally2006; Schunker & Cally 2006 for further details on mode con-version or Cally 2007 for a succinct review of the theory).Khomenko & Collados (2009) have performed simulationswith both monochromatic and gaussian wave sources in amagneto-hydrostatic (MHS) sunspot atmosphere and show aclear correlation between the power halo and a suspiciousincrease in RMS velocities for non-trapped waves resultingfrom the interference pattern generated by downwards travel-ling fast waves.In Rijs et al. (2015) we extended this work in 3D. By per-forming forward modelling simulations with a spatially lo-calised gaussian (in space, time and frequency) wave pulse,the halo structure resulting from the vertical component ofvelocity ( v z ) was analysed as a function of radius, height andfrequency. A clear correlation between the position of the a = c layer and the halo was shown and the dependancy of thehalo on the overlying atmosphere was exhibited.In this work we perform simulations in similar MHS sunspotatmospheres to those of Rijs et al. (2015). However we nowuse a realistic distributed wave source, modelled as a slabof point sources at some depth below the photosphere. Thesources are tuned to mimic the observed photospheric powerspectrum, peaking at the 5 minute oscillation period ( ν = 3 . THE SIMULATION
In this section we present an overview of our simulations,including the details of the sunspot atmosphere used, a sum-mary of the distributed wave source and details regarding thecalculation of synthetic instensities, phase shifts and veloci-ties.
The MHS atmosphere
A detailed description of the sunspot model we are usingcan be found in Przybylski et al. (2015), where the modelof Khomenko & Collados (2008) is optimised in order toincrease spectropolarimetric accuracy and produce morerealistic line formation regions.In short, the MHS configuration combines the self-similarsub-photospheric model of Low (1980) with the potentialconfiguration of Pizzo (1986). Convective stability is en-forced by the method of Parchevsky & Kosovichev (2007).The model makes use of the Model S for the distribu-tion of quiet subphotospheric thermodynamic variables(Christensen-Dalsgaard et al. 1996) and the Avrett umbrafor the non-quiet variables (Avrett 1981). The VALIIICchromosphere (Vernazza et al. 1981) is smoothly joined ontothese distributions to complete the full model, yielding asunspot-like magnetic field configuration embedded into theatmosphere.The sunspots we use in this instance are similar to those usedin Rijs et al. (2015) and Moradi et al. (2015), except for someparameters, such as the peak field strength, the inclinationat the surface, and the simulation box size, which have beenmodified.The sunspot model not only provides the freedom to choosethe peak field strength at the surface of the photosphere butalso the depth of the Wilson depression (the height at whichthe atmosphere becomes optically thin is depressed in highfield regions such as the umbra). As such we make use oftwo model atmospheres in this study, one with a peak surfacefield strength of | B | = 1 . | B | = 2 . τ ) = 0 (where τ is the optical depth scale,as calculated from the known thermodynamic values at everypoint in the box) and follows the Wilson depression. Thesurface corresponds to the height of formation of the 5000 Åintensity continuum ( z = 0) in this atmosphere. Forward modelling
As in our previous work, we use the SPARC code for for-ward modelling (Hanasoge 2007; Hanasoge et al. 2007). Thecode has been used several times for wave-sunspot interaction studies (Moradi & Cally 2013, 2014; Moradi et al. 2015).The code solves the ideal linearised MHD equations in carte-sian geometry.As input, we define a background atmosphere and instigatewave propagation for the desired simulation length. The back-ground atmosphere can be any magnetic plasma such as thesunspot atmospheres mentioned above or any quiet-sun atmo-sphere, provided it is convectively stable.The output is the perturbations to the background states of thepressure ( p ), ρ , magnetic field ( B = [ B x , B y , B z ]) and velocity( v = [ v x , v y , v z ]).The computational domain in both cases is square in the hori-zontal with 256 points in each of the x and y directions (where L x = L y = 200 Mm) yielding a horizontal spatial resolution of δ x = 0 . a ) increases rapidly above thesurface, we use the Alfvén speed limiter described by Rem-pel et al. (2009), which was also used in Rijs et al. (2015).This allows us to sidestep the requirement of using a pro-hibitively small simulation time step, imposed by the Courant-Friedrichs-Lewy (CFL) condition. Work has been done to as-certain whether the use of an Alfvén speed limiter has a detri-mental effect on helioseismic travel times (Moradi & Cally2014), with the conclusion being that one must be certain thatthe artificial capping is occuring well above heights where anyrelevant physics is occuring (such as the fast wave reflectionheight or the a = c layer).We have set our limiter at a value of a lim = 90 km/s, yieldinga simulation time step of around 0.2 seconds.Figure 1 shows a cut through the centre of our 2.7 kG sunspotatmosphere (along the plane at y = 0). Overlaid are the a = c and a = 90 km/s contours, as well as the photospheric surface,with a Wilson depression of around 400 km.The vertical inclination contours show the rather rapid drop-off in field inclination, with the field reaching 30 degrees fromthe horizontal some 20 Mm from the umbra (at the surface).To reiterate, the mode conversion effects occur around the a = c layer, and so it is important that fast waves are givenspace to refract back downwards as they naturally would be-fore the limiter at a = 90 km/s takes effect. We have takencare to ensure that this is the case and that the modificationof the atmosphere will not affect these returning fast waves.In this regard, simulations have been run with Alfvén limitervalues up to 200 km/s, with no change to the halo propertiesobserved.Regarding our wave source function, we are attempting tomodel the uncorrelated stochastic wave field seen on the solarphotosphere. This wave field is, in reality, generated by sub-surface convective cells. We choose a depth of 150 km belowthe surface and initiate a source function, S , in the manner of Figure 1.
A cut through the sunspot center. Field inclination contours are shown for typical umbral/penumbral and penumbral/quiet sun values of 45 and 60degrees from the vertical respectively. The surface or photosphere layer, where log( τ ) = 0, is shown by the solid black curve. The dashed curve is the a = c equipartition layer for this atmosphere and the dash-dotted curve is the a = 90 km/s layer, where the Alfvén limiter is in effect. The background contour is log( a )in km/s as it would appear without any Alfvén limiter in application. In our simulations a is constant above the a = 90 km/s curve. Note that the aspect ratio ofthe figure is highly stretched, with the horizontal dimension spanning 200 Mm and the vertical spanning only around 2 Mm. Figure 2.
Panel a): The power spectrum of the wave source function used, tuned to provide a solar-like peak. b:) Arbitrarily normalised power ridges in (cid:96) - ν space for 6 hours of simulation time, calculated at the surface ( z = 0) from v z . Hanasoge et al. (2007), i.e. S ( x , y , z , t ) = ˆ S ( x , y , t ) f ( z ) (1)where the horizontal function ˆ S ( x , y , t ) is a plane of spatialdelta functions which are excited at a cadence of 30 seconds,and the function f ( z ) is a gaussian function in depth withFWHM of approximately 100 km centered at 150 km belowthe surface.The source power spectrum has been tuned such that it more-or-less fits the spectrum of power observed on the surface ofthe quiet sun. Panel a) of Figure 2 displays this spectrum,with a peak in power at around 3.3 mHz, and non zero powerpresent until above 10 mHz. Panel b) shows the power ridgesin (cid:96) - ν space calculated from 6 hours of v z output at the sur-face.In taking into account the fact that strong umbral fields in-hibit subsurface convection and wave propagation, we do notexcite waves in the umbra of the sunspot itself, smoothly sup-pressing the source amplitude as the magnetic field strengthincreases.Wave propagation is initiated and run for 6 hours of solar timein total using both the 1.4 kG and 2.7 kG sunspot atmospheres(as separate simulations).We analyse the power manifested in synthetic intensities cor-responding to the 5000 Å continuum intensity, the AIA 1700Å and 1600 Å intensity bands as well as both the vertical andhorizontal components of the velocity perturbation ( v z and v h respectively), which correspond observationally to the line-of-sight components of velocity when observing at disk cen-tre ( v z ) and at the limb ( v h ).In reality, the HMI Doppler camera (Scherrer et al. 2012)measures velocities from the Fe I 6173.34 Å line, which hasits peak of formation at a height of around 140 km (Flecket al. 2011; Rajaguru et al. 2013), while the AIA (Lemen et al.2012) 1700 Å and 1600 Å wavelength intensity channels areformed at approximate heights of 360 km and 430 km respec-tively (Fossum & Carlsson 2005; Rajaguru et al. 2013).Thus, in comparing the structure of power enhancementspresent in our simulations with the observed power behaviourfrom Rajaguru et al. (2013) we extract simulation velocity sig-nals from a height of z = 140 km. We then calculate the syn-thetic intensities corresponding to the two above AIA chan-nels as well as the 5000 Å continuum intensity for our 6 hourwave propagation simulations. The approximate 1600 Å and1700 Å intensities are calculated by interpolating the ATLAS9continuum and line opacity tables (Kurucz 1993) using theplasma parameters from the simulation and integrating themtogether with the corresponding LTE source function alongthe lines-of-sight for each column in the sunspot models. Theroutine used for the intensity calculation is similar to that ofJess et al. (2012). The filter bandwidths are set to 10 Å forboth simulated AIA channels. No line-of-sight velocity ormagnetic field information is used in this radiation intensitycalculation. COMPARISONS WITH OBSERVATIONS - VERTICAL VELOCITIESAND INTENSITIES
This section details the comparisons between the powerstructures present in our 6 hour simulations and those ob-served in the active region NOAA 11092.As shown in observations, the acoustic halo is a phenomenonespecially sensitive to | B | and to the local field inclination. We firstly demonstrate here some of the similarities and dif-ferences in these properties exhibited by the artificial sunspotsand the real active region.Figure 3 compares the topology of | B | and the unsigned fieldinclination from the vertical ( γ ) at the surface our 2.7 kGsunspot atmosphere and NOAA 11092.For the observations of NOAA 11092, | B | is calculatedfrom the disambiguated vector maps with components B x , B y and B z , where | B |= (cid:113) B x + B y + B z . γ in degrees is then sim-ply defined as γ = 90 − (180 /π ) | arctan( B z / B h ) | where B h = (cid:113) B x + B y .As can be seen in the figure, the field strength of NOAA 11092drops off in a similar manner to the artificial sunspot, how-ever the small scale features present in the real active regionintroduce many variations in both the field and its inclinationwhich are not modelled in our simulations. The behaviourof γ around NOAA 11092 with radius for example is not thesmooth monotonically increasing function yielded by the 2.7kG sunspot model. We can therefore expect some differencesbetween observed and simulated halo structure will result.Firstly, we compare the acoustic power for v z - from boththe weak (1.4 kG) and strong (2.7 kG) sunspot atmospheres -with the 14 hour time averaged Doppler velocity power fromNOAA 11092.Power maps are shown in Figure 4 for a range of frequenciesof interest. The power at every point has been divided by theaverage power of a quiet corner of the simulation domain, inorder to represent an enhancement over quiet values. In bothsimulations, the enhancement comes into effect at around 5mHz, when waves are in the non-trapped regime, just as inthe observations.The differences between the two sunspot simulations (rows 1and 2) are immediately evident, with the 2.7 kG sunspot ex-hibiting a larger umbra. A consequence of having a strongermagnetic field strength is also that the a = c height will belower in the atmosphere, resulting in a spreading of this con-tour for a particular observation height. It is clear that the haloappears correlated with the a = c contour in both cases.An intriguing feature of the simulated halos is the clear dual-ring structure present for higher frequencies. The inner ringappears to conform qualitatively well at a glance with the ob-servational halo. However the rings appear to be interruptedby a region of mild power deficit (with respect to the quietsun).Although not immediately visible in the power maps in thebottom row of Figure 4, observed halos do exhibit a simi-lar structural change when observed at increasingly high fre-quencies. This feature can clearly be seen in power maps ofobserved Doppler velocity in Rajaguru et al. (2013) and Han-son et al. (2015) at 8 and 9 mHz respectively.In section 5 of this paper we discuss how fast-Alfvén conver-sion likely leads to this dual-ring structure.Comparing power maps in this way is of only so much use. Tomore rigorously compare the structure of observed and simu-lated power halos we plot unfiltered power enhancements asfunctions of | B | and ν (i.e. no frequency filter is applied dur-ing the fourier transform.) In this way we may fully examinethe spectral structure of the halo (Figure 5).Also present in Figure 5 is the power calculated from the AIA1700 Å and AIA 1600 Å intensity bands, which we have syn-thetically calculated in our simulations in order to compare toobservations. Figure 3.
Panels a) and b) show | B | and the unsigned field inclination from vertical ( γ ) respectively for NOAA 11092. Panels c) and d) are the counterpart plotsfor the 2.7 kG simulated sunspot atmosphere. Figure 4.
Top row - 6 hr time-averaged v z power maps at the height of formation of the Fe 6173.34 Å line ( z = 140 km) for 4 illustrative frequency ranges for theweak sunspot case (1.4 kG). Middle row - The same power maps for the stronger field case (2.7 kG). Bottom row - 14 hr time averaged observational Dopplervelocity power maps of the active region NOAA 11092 for the same frequency ranges. The green contour in rows 1 and 2 is the a = c contour at z = 140 km. The left hand panels show the NOAA11092 power structurefor the Doppler velocity and intensities, whereas the righthand panels correspond to simulation output for the 2.7 kGsunspot atmosphere. To be clear, panel b) of Figure 5 corre-sponds directly to the power maps in the middle row of Figure4, it is simply unfiltered in frequency space and so is inclusiveof the entire spectral structure. The power at every point hasbeen binned according to the local value of | B | and then aver-aged so as to reveal not only the spectral structure of the halobut also how it behaves with respect to field strength. The first thing to notice in Figure 5 is that the simulated v z power structure (panel b) matches up reasonably well withthe observed Doppler power (panel a). The halo has formedover relatively weak field ( 50 G < | B | < 700 G) as expected.In the simulation, | B | decreases (and γ increases) smoothlyand uniformly as one moves away from the umbra. As suchthis field strength range corresponds to nearly horizontal in-clinations of 55 ◦ < γ < 75 ◦ .This seems to also agree with all other observational reportsof enhancements which place the halo amongst moderate to Figure 5.
Panel a): Unfiltered Doppler velocity power as a function of | B | and ν . b): 2.7 kG simulation unfiltered v z power. c) & d): observed and 2.7 kGsynthetic unfiltered AIA 1700 power respectively. e) & f): observed and 2.7 kG synthetic unfiltered AIA 1600 power respectively. weak and horizontally inclined field (Jain & Haber 2002;Schunker & Braun 2011; Rajaguru et al. 2013).The dual-ring structure can clearly be seen at higher frequen-cies in panel b), manifesting as the second lobe of enhance-ment for very weak field. Wedged between the two rings (ataround ( | B | , ν ) = (100 ,
6) is the clear region of power reduc-tion.Looking at greater heights in the form of the AIA 1700 Å and1600 Å intensities (corresponding to heights of 360 km and430 km above the base of the photosphere respectively) wealso see a general agreement in ν, | B | space. The spreading ofthe magnetic canopy at these heights has resulted in the inten-sity halos forming at much weaker field locations both in theobservations and the simulations.The magnitudes of the enhancements in the simulations areconsistently larger than the observed values, as evident fromthis figure. This is a feature that was also noted in Rijs et al.(2015) and can most likely be attributed to the fact that oursunspot is symmetric and its magnetic field inclination is asteep, monotonically decreasing function of radial distance.The MHS structure is such that horizontal field is enforcedat the side boundaries of the simulation domain and so thereis a large expanse of nearly horizontal field. As explainedpreviously, the fast-slow mode conversion mechanism for thegeneration of the halo relies on a large attack angle betweenwavevector and field and so, in analysing v z power enhance-ments, it is reasonable to expect that this horizontal field willbe very conducive to the conversion of energy into magneticfast waves and hence, a prominent halo.Power derived from the 5000 Å intensity continuum (at z = 0)was also calculated synthetically to compare with the observa-tional intensity continuum power. It is well known that halosdo not appear in measurements of intensity continuum powerand we also found this to be the case, with no enhancementpresent.Another intersting result, shown in Figure 6, is the compari-son between observed and simulated phase shifts. A net up-ward or downward propagation of waves in an atmosphere can be diagnosed by calculating the temporal cross-spectrumof any wave quantity sampled at two different heights.For example, for velocities v ( z , t ) and v ( z , t ) sampled at twodifferent heights z and z , the phase shift corresponding to aheight evolution of the wave is given by the argument or phaseof the complex cross-spectrum, φ , ( ν ) = arg[ V ( z , ν ) V ∗ ( z , ν )] , (2)where V is the Fourier transform of v . In the above conven-tion, a positive phase-shift would mean that the wave is prop-agating from height z to z , while the opposite holds for anegative phase-shift.The phase shift contour maps of Figure 6 describe the phaseshifts of waves at the AIA 1700 Å and 1600 Å intensity for-mation ( z = 360 km and 430 km respectively) with respect tothose at the height of formation of the intensity continuum.The simulation yields a clean band of positive phase shifts athalo frequencies with respect to those at the surface at weakfield regions. The same basic pattern is seen in the observa-tions, however there is some extended phase shift structureat higher field strengths in the AIA 1700 Å power (panel a)which is not replicated in the simulation.The simulation phase shifts are also of a greater magnitudethan observations - particularly in the case of the AIA 1600intensities.These variations in features are not too surprising. Consider-ing Figure 3 we see that NOAA 11092 exhibits a much morerapid horizontality of field away from the umbra than seen inthe MHS model. We show in section 5 how these bands ofpositive phase shifts at given observation heights may be in-trinsically related to the process of fast-Alfvén mode conver-sion. The physics of fast-Alfvén mode conversion are stronglytied to the local magnetic field inclination. Therefore the rea-son that NOAA 11092 exhibits such an extended phase shiftstructure into higher field regions (and our MHS sunspot doesnot) may be in part due to the more horizontal field at thoseradii for the active region. THE THEORETICAL UNDERPINNINGS OF HALOS
Figure 6.
Phase shifts at the heights of formation of the AIA 1700 Å and 1600 Å lines of all waves with respect to those at the surface. Panels a) and c)correspond to observations and b) and d) to the 2.7 kG simulation.
In order to prove that the halo is produced by the return ofreflected fast waves, we examine several intriguing featurespresent in our simulations. Firstly, in a similar manner toRijs et al. (2015) we perform several identical simulations tothe 2.7 kG case examined above, except with incrementallysmaller Alfvén limiter values.After undergoing mode conversion at around the a = c height,fast magnetic waves will begin to refract and then ultimatelyreflect at the point in the atmosphere where the horizontalphase speed equals the Alfvén speed (i.e. where ω/ k h = a ).By reducing the height of the artificial ‘cap’ on the atmo-sphere we are allowing less and less room for fast waves to re-fract and deposit extra energy onto observable heights. Wavesthat impinge on the altered region of constant a will simplytravel upwards and out of the local area. As the original sim-ulation had a value of a lim = 90 km/s, we run simulations with a lim = 40, 20 and 12 km/s and analyse the power in a similarmanner to Figure 5, i.e. as a function of | B | and ν . The com-parison is shown in Figure 7.Panel d) corresponds to the case where the limiter is onlybarely above the a = c height, enabling the mode conversionto take effect but yielding virtually no room for fast waves toreturn. Moreover in the intermediate cases of panels b) andc), the magnitude is reduced as the more vertically orientedwaves are escaping to the top of the box, yielding contribu-tions from only the more horizontally inclined waves.Clearly the halo is entirely dependent on the overlying atmo-sphere and by restricting the refraction and return of the fastwaves the enhancement is entirely absent.The second theoretical check we perform is to compare thestructure of the halo resulting from both the horizontal andvertical components of the velocity. A reasonable attack an-gle between the horizontal component of the wavevector, k h and B is still entirely likely in our simulations, as the fieldis never entirely horizontal. Also, as noted by Khomenko &Collados (2009) we can expect that it would at least be of asimilar strength to the v z halo, as waves are largely horizontalat around the refraction height. Figure 8 shows the comparison. A clear feature is that the v h enhancement occurs at preferentially higher field strengththan the v z enhancement. This feature also makes sense asthe field inclination is more vertical at these radii, providing alarger attack angle.It would be extremely useful if there were any center-to-limbobservational studies of halo features, so that we could com-pare the horizontal Doppler component with our v h . Zharkovet al. (2009) have performed an analysis of the umbral "bellybutton" as a function of observation angle, but as of yet, nosuch studies focusing on halo properties have been conducted.Finally, and perhaps most importantly, we explain the "dual-ring" power enhancement structure seen in the power mapsearlier and in observations.In Figure 9 we compare v z power (once again at the standardobservational height of 140 km) for both the 1.4 kG and the2.7 kG simulations with the phase shifts at the same height.The phase shifts in this case are those calculated at z = 140km height, with respect to waves at z = 0, so we are only look-ing at the phase shifts that the waves experience over a heightchange of 140 km in the simulation.The black curves have been added simply by eye to aid in thecomparisons here. In both simulations there is a similar phaseshift pattern to that observed in both the simulated and ob-served intensities at greater heights, however the magnitudeis less here as the waves have travelled a shorter vertical dis-tance.The key fact to note is that the strong branch of positive phaseshifts corresponds precisely to the region between the dualrings of power enhancement. This enhancement gap in | B | , ν space is the dark ‘moat’ seen between the two halo rings atvarious frequencies in the simulation power maps of Figure 4.The halo itself shows no real phase shift which most likelyindicates a mixture of upwards and downwards travellingwaves. This is to be expected at high, non-trapped frequenciesas waves rise upwards towards the a = c layer and are refractedback downwards. The halo structure itself does not appearto change too significantly with respect to the peak magnetic Figure 7.
Panel a): Unfiltered v z power halos in the case a lim = 90 km/s. b), c) and d) show the same quantity from simulations with progressively lower valuesof a lim . Figure 8.
Panel a) is the standard binned v z power for the 2.7 kG atmosphere. Panel b) is the binned power corresponding to the horizontal component ofvelocity, v h . The dashed vertical line is the position of the a = c for the observational height of z = 140 km. field strength of the model, apart from the noted correlationwith the a = c layer. We certainly do not see any noticeablechange in the peak halo frequency, as Khomenko & Collados(2009) suggested may be the case. This is most likely dueto the fact that, although the peak field strengths of the twomodels are considerably different in the umbra (1.4 kG and2.7 kG), at the halo radius (some 20 Mm out) the differencein the field strength will not be so significant.The pertinent question is: why are there only upwards travel-ling waves in the moat in between the concentric halos? FAST WAVE DAMPING AND ALFVÉN WAVES
The answer would appear to lie in the process of fast-Alfvénmode conversion, the basics of which are described in Cally (2011) and Cally & Hansen (2011).Fast-Alfvén mode conversion has been well studied in bothsunspot-like (Moradi & Cally 2014; Moradi et al. 2015) andsimple magnetic field geometries: Pascoe et al. (2011, 2012)have studied the damping of transverse kink waves in termsof the associated Alfvén resonance and Cally & Goossens(2008) and later Khomenko & Cally (2011) have conductedparameter studies with monchromatic wave sources and sim-ple inclined field magnetic structures. The finding of the lattertwo works was that fast wave energy is converted to the field-aligned Alfvén wave at favoured field inclinations ( θ ) andwavevector-to-field angles ( φ ). The process is also stronglydependent on both ν and k h .0 Figure 9.
The top row corresponds to the weak field simulation, with peak field strength 1.4 kG, and the bottom row is from the strong field case (2.7 kG peak).On the left are phase shifts calculated at z = 140 km in height. On the right are the standard binned and unfiltered v z power distributions. The black curves aredrawn by eye to denote where the phase shifts would be in the power plots. Once again, the dashed vertical line is the position of the a = c for the observationalheight. In the case of our distributed source simulations, waves ex-hibit a distribution of wavenumbers and frequencies in a sim-ilar manner to the quiet sun and so the picture is somewhatmuddied in comparison to such simulations. We can expecthowever that fast-Alfvén conversion will in some way act onfast waves as they reach the Alfvén resonance near their upperturning point (on the order of a few hundred kilometres abovethe a = c , depending on k h ).As the halo appears to be generated by downwards turningfast waves, we would anticipate that some of this returningenergy may be lost to the field aligned Alfvén wave, whichwill follow the local field lines until reaching the top (or theside) of the simulation domain.In Figure 4 we noted the strong concentric halos and the gapof power enhancement in between them. Figure 9 shows thismore comprehensively and associates this dark ring with astrong positive phase shift.We suggest that the reason that the halo is not one continuousregion is that for specific field inclinations, fast mode energyis lost to the Alfvén wave.Figure 10 suggests this to be the case. Each panel of the figurecorresponds to a specific frequency filtering. The top halves ofthe panels are the same as the panels in the middle row of Fig-ure 4, i.e. filtered power maps corresponding to the strongerfield 2.7 kG simulation at the Doppler velocity observationalheight of 140 km.The bottom halves of the panels show the magnetic energyassociated with the Alfvén wave in the form of the Poyntingvector, S , where S = 1 µ ( − v × B ) × b , (3)where v and b indicate the perturbations to the velocity andthe background field respectively. The bottom panels showthe vertical component of the vector, S z , corresponding to theupcoming Alfvén flux, and are calculated at the very top ofthe simulation domain, at a height z = 2 Mm, just before the PML comes into effect at the top of the box. In each case thevelocity has been pre-filtered around the associated frequencyrange prior to the calculation of S z to match the power maps.It is worth remembering that we have applied a cap to a above a = 90 km/s in the atmosphere and so any upwards travellingAlfvén waves will encounter our modified atmosphere andtravel at a constant speed to the top of the box, instead ofbeing subject to a rapidly increasing Alfvén speed.The correlation between the Alfvén flux and the positionof the dark ring is immediately noticeable, especially inthe 5.5 and 6.5 mHz cases. Note that upwards travellingAlfvén waves will follow the field and that there is some fieldspreading with height in this MHS atmosphere which is why S z is diffuse and does not align precisely with the dark ring atobservation heights.It seems clear that this Alfvén wave energy is responsible forthe strong band of positive phase shifts (and thus upwardstravelling waves) in the dark moat. There is no wave energyleft to return downwards at these radii and field inclinations.Furthermore this supports the fast-wave halo mechanismrather strongly as the two processes are critically interlinked. DISCUSSION AND CONCLUSION
Linear forward modelling in realistic MHS sunspot atmo-spheres has yielded acoustic halos that match up quite wellwith observations, both spatially and spectrally. Apart fromthe magnitudes of the enhancements themselves, most ob-served features seem to be reproduced in our simulations, notjust when comparing Doppler and vertical velocities, but alsointensity halos at multiple heights in the chromosphere. Asin the observations we see no power enhancement in calcula-tions of the time-averaged intensity continuum power.We have also presented convincing evidence that the mecha-nism responsible for halo formation is the refraction and re-turn of magneto-acoustic fast waves at non-trapped frequen-cies. The halo appears very sensitive to the position of the a = c layer in the atmosphere, which is the critical loca-1 Figure 10.
Power map - Poynting vector composites. Top halves are xy power maps at z = 140 km, filtered around the respective frequencies. Bottom halves are S z in units of ergs/ m s , calculated at z = 2 Mm. Note that the Poynting vector scaling is not consistent from plot to plot, as there is significalty less energy arrivingat the top of the box for each subsequently higher non-trapped frequency range. tion for fast-wave mode conversion. With our realistic dis-tributed wave source, we see a strong relationship betweenthe strength of the halo and the extent to which fast waves areallowed to return downwards. This suggests that the halo iscompletely governed by the overlying a > c atmosphere andthe extra energy injected to observable heights by these re-turning fast waves.The theory also predicts that an enhancement should bepresent in the power of v h , as this component will also in-teract with the field, and that this enhancement should be con-centrated toward more vertical field (as the horizontal compo-nent makes a larger attack angle with vertical field); this wasshown to be the case as well. Unfortunately center-to-limbobservational studies of the halo do not yet exist and so wecannot compare this horizontal velocity halo to the real thing.Our simulations are performed in a MHS atmosphere, solvingthe linear MHD equations and using a wave excitation mech-anism that approximates the wave bath of the solar photo-sphere. The fact that we see halos in such simulations (whichare of course, entirely non-radiative and do not in any way in-clude convective effects) suggests that the halo is not createdby any convective cell-magnetic field interaction as suggestedby Jacoutot et al. (2008).Similarly, the idea of Kuridze et al. (2008) that m > v z halosare greater than this by a factor of 2 or even 3, depending onfrequency.There are several possible explanations for this discrepancy.As we have postulated, the halo enhancement most likely oc-curs as a result of fast waves interacting with the sunspot mag-netic field at large attack angles. This yields a large conver-sion fraction to the fast magnetic wave which refracts anddeposits additional energy in the photopshere and chromo-sphere. The penumbral field structure of active regions dif-fers significantly from the simple MHS model used here how-ever. Our atmosphere does not explicitly include an umbra orpenumbra, but rather consists of a smoothly decreasing fieldstrength and vertical inclination component, yielding signifi-cant regions of smooth, nearly-horizontal field. Non-trappedwaves which reach the a = c equipartition layer will have alarge vertical component and so we would expect these waves2to interact strongly with primarily horizontal field. In na-ture, penumbrae contain fine structure, with bright and darkfilaments giving rise to the now well-observed combed mag-netic field configuration (Scharmer et al. 2002; Bellot Ru-bio et al. 2004). At the outer penumbral boundary, studieshave shown up to a 60 ◦