{'fraction_non_alphanumeric': 0.034641906933811153, 'fraction_numerical': 0.02367377966009665, 'mean_word_length': 4.468527315914489, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sensitivity of helioseismic travel-times to the imposition of a Lorentz force limiter in computational helioseismology

The rapid exponential increase in the Alfv\'en wave speed with height above the solar surface presents a serious challenge to physical modelling of the effects of magnetic fields on solar oscillations, as it introduces a significant CFL time-step constraint for explicit numerical codes. A common approach adopted in computational helioseismology, where long simulations in excess of 10 hours (hundreds of wave periods) are often required, is to cap the Alfv\'en wave speed by artificially modifying the momentum equation when the ratio between Lorentz and hydrodynamic forces becomes too large. However, recent studies have demonstrated that the Alfv\'en wave speed plays a critical role in the MHD mode conversion process, particularly in determining the reflection height of the upward propagating helioseismic fast wave. Using numerical simulations of helioseismic wave propagation in constant inclined (relative to the vertical) magnetic fields we demonstrate that the imposition of such artificial limiters significantly affects time-distance travel times unless the Alfv\'en wave-speed cap is chosen comfortably in excess of the horizontal phase speeds under investigation.


Introduction
Sunspots are large cool regions on the solar surface associated with strong magnetic flux concentrations. They represent an important connection of the internal magnetic field of the Sun with the solar atmosphere. Their subsurface structure is a challenge to both theory and observation, and in particular to local helioseismology, which cannot yet fully deal with the dominant magnetic fields near the surface (for recent reviews see Moradi 2012;Gizon et al. 2010;Moradi et al. 2010). The ultimate aim is to better understand solar activity by seismically probing flows, magnetic fields, and thermal structures in sunspots and the wider active regions that host them.
A major issue in local helioseismology is the complex behaviour of trapped acoustic waves (p-modes) when they encounter an active re-gion (Cally & Bogdan 1993;Cally et al. 1994;Schunker & Cally 2006;Cally 2007;Cally & Goossens 2008;Schunker et al. 2013). The magnetic field turns them into a complex mixture of fast, slow and Alfvén waves through MHD mode conversion processes. The trapped modes can then partially escape into the solar atmosphere and partially reflect there to rejoin the internal wave field, with significant consequences for the local seismology.
Physical modelling of the effects of magnetic fields on solar oscillations can help explain helioseismic measurements in regions of strong magnetic field. Through forward modelling -the art of constructing computational models that mimic wave propagation through sunspots and matching the resultant wave statistics with observations -we have been able to gain much valuable insight into the interaction of helioseismic waves with magnetic fields (e.g., Cameron et al. 2008Cameron et al. , 2011Hanasoge 2008;Khomenko et al. 2009;Shelyag et al. 2009;Moradi et al. 2009;Schunker et al. 2013).
Simulating 3-D MHD wave propagation in the Sun is a complex and time-consuming task however, with one of the most pressing issues being how to treat the excessively large Alfvén wave speed (c a = B 0 / √ 4πρ 0 ; where B 0 and ρ 0 represent the magnetic field strength and density respectively) above the surface, which is brought about by the exponential drop in ρ 0 with height. With c a reaching several thousands of km s −1 , this leads to an extremely stiff numerical problem for the explicit numerical solvers, as the timestep (∆t ≈ ∆z/c a , where ∆z denotes the vertical grid resolution) is constrained by the Courant-Friedrichs-Lewy (CFL) condition, resulting in the need for very small ∆t for simulating even moderate magnetic field strengths.
Consequently, in order to simulate helioseismic data sets in a feasible amount of time, the most common approach in computational helioseismology has been to introduce a Lorentz force (F L ) scaling factor that limits it when the ratio between Lorentz and hydrodynamic forces (or in other words c a /c s , where c s is the sound speed) becomes too large. Typically F L is scaled by αc 2 s /(αc 2 s + c 2 a ) (where α is a free parameter that controls the amplitude of the limiter), resulting in c a being capped above the surface, commonly in the range of ∼ 20−60 km s −1 (Rempel et al. 2009;Cameron et al. 2011;Braun et al. 2012). Another similar approach has been to scale the magnetic field by a pre-factor such that c a never exceeds a certain predefined value (Hanasoge et al. 2012). The physical implications of artificially limiting c a in such manners have not been fully explored to date, with the general assumption being that the overlying atmosphere does not play a significant role in the seismology of the subphotosphere.
However, a number of recent studies have cast doubt on this assumption by demonstrating the critical role played by the exponentially increasing c a in the fast-to-Alfvén mode conversion process which takes place in the lower solar atmosphere (Cally & Goossens 2008;Cally & Hansen 2011;Hansen & Cally 2012;Khomenko & Cally 2011Felipe 2012). Using idealised MHD simulations, these studies show how the upward propagating helioseismic (fast) wave can reflect off the c a gradient back to the surface at the 'fast-wave reflection height' where the horizontal phase speed (v ph = ω/k h ; where ω denotes angular frequency and k h the horizontal wave number) roughly coincides with c a . These 3-D calculations also show that, depending on the local relative inclinations and orientations of the background magnetic field and the wavevector, the fast wave may undergo partial mode-conversion to either an upward or downward propagating Alfvén wave (see Figure 1 in Khomenko & Cally 2012) around the reflection height where they are nearresonant.
While the fast-to-Alfvén mode conversion process generally takes place above the surface, it has potentially serious consequences for helioseismology. This is because after they reflect off the c a gradient, the fast waves re-enter the solar interior wave field, meaning that their journey through the atmosphere and the phase changes they suffer in the conversion process must have some effect on the seismology. Inversions of observed timedistance travel times (e.g., Kosovichev et al. 2000;Couvidat et al. 2005) would normally but mistakenly interpret such phase changes as 'travel time shifts' due to subsurface inhomogeneities.
In a recent follow-up study,  quantified the implications of the returning fast and Alfvén waves for the seismology of the photosphere. They found substantial wave travel time shifts that acutely depend on magnetic field inclination and wave propagation orientation, in direct correspondence with the escaping acoustic and Alfvénic wave fluxes above the surface. However in another related study, it was observed that the imposition of a F L limiter actually suppresses the fast-to-Alfvén mode conversion process by artificially diminishing the Lorentz force above the surface . But nonetheless, reality (and indeed, finite computational resources) dictates that simulating artificial data sets on the time scales required for local helioseismic analysis will ultimately necessitate the use of a limiter to ensure a manageable CFL condition. Whether this has a flow-on effect to the helioseismic travel times is something that we wish to explore here, by modelling wave propagation through constant inclined magnetic fields and comparing the artificial helioseismic travel times derived from simulations with and without a F L limiter.

Numerical Model
We use the Seismic Propagation through Active Regions and Convection (SPARC; Hanasoge 2007) code to conduct the forward modelling component of our analysis. SPARC solves the 3-D linearised MHD equations in Cartesian geometry to investigate wave interactions with local perturbations. We employ a similar numerical setup to , with a 3-D computational box spanning 26.53 Mm in height z (covering −25 z 1.53 Mm) using 265 vertical grid points, with the grid spacing ranging from several hundred km deep in the interior to tens of kms in the near-surface layers. The box length is chosen to be 140 Mm in the horizontal directions x and y. We use 128 evenly spaced grid points in x and y, resulting in a horizontal resolution of ∆x = ∆y = 1.09 Mm/pixel. The vertical and horizontal boundaries of the box are absorbent, with perfectly matched layer (PML) boundary layers at the top and bottom and absorbing sponges lining the sides.
The background model consists of a convectively stabilised solar model (CSM B; Schunker et al. 2011) threaded by uniform magnetic field of B 0 = 500 G, inclined at angle θ (inclination from vertical, with 0 • θ 90 • ). In  we employed random stochastic wave sources to generate acoustic waves, but in this study we employ a perturbation source in vertical velocity (v z ) similar to Shelyag et al. (2009): (1) where t 0 = 300s, t 1 = 300s, σ t = 100s, σ r = 4∆x, and r 0 is the source position, located 160 km below the surface at (x 0 , y 0 ) = (0, 0). This source generates a broad spectrum of acoustic waves in the 3.33 mHz range (see Figure 1 a), mimicking wave excitation in the Sun.
We initially use the acoustic source to simulate wave propagation for ten different field inclinations (θ = 0 • , 10 • , 20 • ..., 90 • ), without invoking a limiter. The time-step required for each of these simulations in SPARC is ∆t = 0.1 s. We then repeat the simulations using the form of the limiter adopted by Hanasoge et al. (2012) to cap c a at a where ρ 0 denotes density (the subscript 0 indicates a time-stationary background quantity, whereas unsubscripted terms fluctuate), p the pressure, B magnetic field, B * 0 = √ κB 0 (where κ = αc 2 s /(αc 2 s + c 2 a )), g gravity (with direction −ê z ), and S the source term.
The first cap is placed at 20 km s −1 (α = 8, resulting in ∆t = 2.0 s), followed by 40 km s −1 (α = 30, ∆t = 1.0 s), 80 km s −1 (α = 125, ∆t = 0.5 s), and finally 160 km s −1 (α = 540, ∆t = 0.25 s). Figure 1 b) plots the various c a profiles of each simulation. Hence, in total, we conduct 50 unique simulations, each with a temporal duration of 1 hour. In addition to the magnetic simulations, we also complete a quiet-Sun reference simulation, using only the vertical stratification of the thermodynamic parameters in the CSM B model.
For the proceeding helioseismic analysis we use one hour simulated v z data cubes extracted at a constant geometrical height of 300 km above the surface.

Time-Distance Analysis
In time-distance helioseismology wave travel times are generally calculated from fits to the temporal cross correlation function between two points (source and receiver) at the solar surface (Duvall et al. 1993). However since wave excitation in our computational box is generated by a single source, there is no need to compute the velocity correlations between different surface points, meaning that time-distance diagrams can be constructed by plotting v z as functions of time for all horizontal points (see Figures 1 c, d). Moreover, since the simulations are in 3-D, by simply selecting a receiver point away from the central axis in the xy-plane, we have the ability to choose the magnetic field orientation with respect to the vertical plane of wave propagation, which we refer to as the "azimuthal" field angle (φ, where Prior to calculating the travel times, we first filter the data cubes in two frequency ranges: 3 and 5 mHz by employing a Gaussian filter with a dispersion of 0.5 mHz. We then measure the phase travel time perturbations δτ (i.e., the differences in the phase travel times between the magnetic and nonmagnetic simulations) using Gabor wavelet fits (Kosovichev & Duvall 1997) to the time-distance diagram at particular wave travel distances (∆) away from the source, as functions of field inclination (θ) and azimuthal direction (φ). A rectangular window of width 14 minutes centred on the first-bounce ridge selects the fitting interval in time lag. The fits are done by minimising the misfit between the Gabor wavelet and the wave form. Figure 1 e) shows an example of a wave form and the associated best-fit Gabor wavelet from one of our calculations.
We measured δτ for a range of distances, but for the sake of brevity we present the results for ∆ = 11.6, 24.35 and 42.95 Mm below. The horizontal phase speeds v ph associated with these distances are 16.3, 34.8 and 46.8 km s −1 respectively 2 .

Results & Analyses
The contour plots in Figures 2 and 3 show the resulting 3 (left column) and 5 (right column) mHz δτ for ∆ = 11.6 and 42.95 Mm respectively, as functions of θ and φ. The bottom panels (g, h) of each figure represent the results derived from the simulations in which we did not impose a limiter. The general pattern of δτ for these cases are strikingly similar to those presented in , with the behaviour being strongly linked to mode conversion in the atmosphere (see their Figures 2 and 5). A detailed discussion is presented in , which we summarise here in short: 1977; Jefferies et al. 2006) the upward propagation of acoustic waves into the atmosphere is severely inhibited due the acoustic cutoff frequency (ω c , being just over 5 mHz in the atmosphere), resulting in δτ values of a few seconds being recorded.
• At intermediate φ, the fast wave loses more energy near its apex to the Alfvén wave, contributing a positive δτ that partially cancels the underlying negative δτ . Panels (a-f) in Figures 2 and 3 show the δτ derived from simulations which have F L limiters imposed to cap c a at 20 (a, b), 40 (c, d) and 160 (e, f) km s −1 . When comparing the overall δτ behaviour in these panels to those derived in the absence of a limiter (g, h), we notice distinctive differences in both the magnitude and general behaviour of δτ across both θ and φ, evident at both 3 and 5 mHz and being most prevalent around θ ≈ 60 • − 90 • and φ ≈ 90 • (a similar δτ pattern was also observed for ∆ = 24.35 Mm, but is not shown here). This can be seen more clearly in Figures 4 and 5, which plot δτ for θ = 70 • and 80 • across φ, for all three ∆ and F L limiters studied. Differences of up to 20 seconds in δτ can be observed around φ = 80 • with the c a cap at 20 km s −1 . The situation improves somewhat as the cap is lifted to 40 km s −1 , but differences of up to ∼ 10 seconds still persist for some distances. But as the c a cap is raised to 80 and then 160 km s −1 , the differences in δτ become progressively smaller and we observe Mm. Left column represents 3 mHz and right 5 mHz. Note that 5.2 cos 70 • = 1.8, so even 3 mHz is well above the ramp-reduced atmospheric cutoff at these high inclinations.
a convergence to the δτ values derived without a limiter. We did not observe complete convergence, even with the c a cap at 160 km s −1 . This behaviour is not completely unexpected. As discussed earlier, the fast wave reflection height is determined by where v ph ≈ c a . This implies that by artificially limiting c a above the surface, we are allowing fast waves with v ph above c a to reach the absorbing PML layer at the top of the computational domain and therefore never to return to the subsurface seismic field. As is evident  in Figures 2-5, this has profound consequences for travel-time measurements if c a is capped below the v ph associated with the chosen ∆. This explains why we observe the most significant differences in δτ associated with the c a cap at 20-40 km s −1 (remembering that the largest v ph sampled, associated with ∆ = 42.95 Mm, is 46.8 km s −1 ), and at the relative inclinations and azimuthal orientations typically associated with maximal fast-to-Alfvén mode conversion (Cally & Goossens 2008;Khomenko & Cally 2011Felipe 2012). Only once the c a is placed comfortably above the v ph being studied do we begin to see a convergence to the δτ derived without a limiter. Furthermore, since fast-to-Alfvén conversion is also spread over many scale heights for wavenumbers typical of local helioseismology (Cally & Hansen 2011), it explains why we don't observe a complete convergence in δτ , even with the c a cap at 160 km s −1 . It is also worth noting that even in some of the cases with the limiter set at 160 km s −1 , the result-ing δτ discrepancies of ∼ 1 − 2 seconds are comparable with the estimates from Schunker et al. (2013) for the sensitivity of travel times to changes in the subsurface structure of sunspots. In fact, even 160 km s −1 could be insufficient in some cases for practical helioseismology with large ∆. Finally, we note that we also conducted a number of test cases using the same simulation setup as described in Section 2, but with B 0 = 1 and 1.5 kG, and apart from larger amplitudes in δτ , the overall results were almost identical to those derived from the 500 G cases.

Summary & Discussion
Computational helioseismology generally entails high resolution, long (temporal) duration 3-D MHD simulations to study the interaction of helioseismic waves with magnetic fields. To alleviate the severe CFL time-step constraints introduced by the exceedingly high c a above the surface in magnetic regions, a number of explicit numerical codes limit the strength of F L in low plasma-β regions, essentially capping the c a gradient above the surface in the process. While this approach can increase the explicit time-step limit to any desired or practical value, it also severely impacts on the fast-wave reflection height (c a ≈ v ph ), and thus by extension the fast-to-Alfvén mode conversion process, which recent studies have shown to be problematic for helioseismology.
Using 3-D MHD simulations of waves in homogenous inclined magnetic fields we find that, in the absence of a F L limiter, time-distance δτ are sensitive to magnetic field inclination and wave propagation orientation, consistent with the recent results of . We also find that the imposition of an artificial F L limiter can have a significant impact on these δτ , unless the c a cap is placed well above the horizontal phase speed associated with the wave travel distance being studied -thus assuring that the reflection height of helioseismic fast wave in the lower solar atmosphere remains (relatively) unharmed.