Can Proton Beam Heating Flare Models Explain Sunquakes?

Solar Dynamics Observatory (SDO)/Helioseismic and Magnetic Imager (HMI) observations reveal a class of solar flares with substantial energy and momentum impacts in the photosphere, concurrent with white-light emission and helioseismic responses, known as sunquakes. Previous radiative hydrodynamic modeling has demonstrated the challenges of explaining sunquakes in the framework of the standard flare model of “electron beam” heating. One of the possibilities to explain the sunquakes and other signatures of the photospheric impact is to consider additional heating mechanisms involved in solar flares, for example via flare-accelerated protons. In this work, we analyze a set of single-loop Fokker–Planck and radiative hydrodynamics RADYN+FP simulations where the atmosphere is heated by nonthermal power-law-distributed proton beams which can penetrate deeper than the electron beams into the low atmospheric layers. Using the output of the RADYN models, we calculate synthetic Fe i 6173 Å line Stokes profiles and from those the line-of-sight observables of the SDO/HMI instrument, as well as the 3D helioseismic response, and compare them with the corresponding observational characteristics. These initial results show that the models with proton beam heating can produce the enhancement of the HMI continuum observable and explain qualitatively the generation of sunquakes. The continuum observable enhancement is evident in all models but is more prominent in ones with E c ≥ 500 keV. In contrast, the models with E c ≤ 100 keV provide a stronger sunquake-like helioseismic impact according to the 3D acoustic modeling, suggesting that low-energy (deka- and hecto-keV) protons have an important role in the generation of sunquakes.


INTRODUCTION
The energy release process in solar flares affects all layers of the solar atmosphere, from the photosphere to the corona.The standard 'thick-target' flare model assumes that a substantial part of the flare energy is released in the solar corona in the form of a high-energy (deka-keV) electron distribution traveling downward along magnetic field lines and heating the upper chromosphere (Hudson 1972).Radiative hydrodynamic modeling of the atmospheric response to the electron beam heating, initiated by Kostiuk & Pikelner (1975), revealed upflows of chromospheric plasma into the flare loops (commonly referred to as 'chromospheric evaporation'), accompanied by a dense downward propagating shock (chromospheric 'condensation') that is strongly compressed due to the radiative losses (Livshits et al. 1981;Fisher et al. 1985;Kosovichev 1986).These hydrodynamic models found that the downward-moving shock quickly decays, within 60 s or so, which is borne out by observations (see also Ashfield & Longcope 2021;Ashfield et al. 2022;Kerr 2022Kerr , 2023)).This has the implication that the downward-moving shock is insufficient for explaining deep perturbations of the solar photosphere.Livshits et al. (1981) suggested that the white-light emission can be produced by chromospheric condensations (see also Kowalski et al. 2015Kowalski et al. , 2017, who studied both solar and stellar models of continuum emission from chromospheric condensations).The radiative back-warming process can also play an important role in white-light emission generation (Machado et al. 1989).
Using observations from the Michelson Doppler Imager (MDI) on Solar and Heliospheric Observatory (Scherrer et al. 1995), Kosovichev & Zharkova (1998) showed that flares can produce a significant impact in the photosphere, which is sufficient for the generation of helioseismic waves ('sunquakes').The localized impulsive impacts and sunquakes are also observed in the photospheric observations of the Helioseismic and Magnetic Imager (HMI) onboard the Solar Dynamic Observatory (Scherrer et al. 2012).Such impacts are typically accompanied by enhancements of the continuum emission (for example, all 18 sunquake events studies by Buitrago-Casas et al. 2015, were accompanied by such enhancements) as well as strong variations in Doppler shift and magnetic field and are identified as sunquake sources (Sharykin & Kosovichev 2020).In the recent work, Wu et al. (2023) suggested that the high-energy tail (E e > 300 keV) of a non-thermal electron distribution is a preferred driver of the sunquakes.This result was obtained based on the analysis of 20 strong flares of the Solar Cycle 24 and fitting the hard X-ray emission spectra observed by The Reuven Ramaty High-Energy Solar Spectroscopic Imager (RHESSI, Lin et al. 2002).
Recently, the flare community has made extensive use of radiative hydrodynamics simulations performed using the RADYN code, developed by Carlsson & Stein (1992) and adopted for flare modeling by Abbett & Hawley (1999) and Allred et al. (2005Allred et al. ( , 2015)).In agreement with the earlier models, these RADYN simulations demonstrated that electron beams can only weakly affect the photospheric layers through direct heating.Sadykov et al. (2020) performed modeling of the Fe I 6173 Å Stokes profiles and corresponding SDO/HMI LOS observables for single-loop RADYN electron beam-driven simulations.Those simulations were available as a part of the F-CHROMA project (Carlsson et al. 2023, http://www.F-CHROMA.org/).The highest HMI continuum intensity observable 1 enhancement of about 3%, accompanied by HMI observable Doppler velocities of ∼0.4 km s −1 , were found for the model with the total energy of E total =10 12 erg cm −2 , low cut-off energy of E c = 25 keV, and a power-law spectral index δ = 3.The electrons were injected for 20 s into the solar atmosphere in a triangular profile that peaked at t = 10 s for F-CHROMA models, so that correspondingly, the average injected flux in that model was equal to F d = 5 × 10 10 erg cm −2 s −1 .While the perturbations of the SDO/HMI observables in that grid of simulations could not explain the SDO/HMI derived continuum intensity enhancements observed during white light flares (10-100% depending on the flare's soft X-ray class, Song & Tian 2018), models with higher electron beam fluxes (Kowalski 2022) are expected to show a significant white light emission.The F-CHROMA models also clearly do not result in the velocity signals of several km s −1 , which are needed for the initiation of sunquakes (Stefan & Kosovichev 2020).
One of the possibilities to explain the deep perturbations of the photosphere, including sunquakes, is to consider additional heating mechanisms involved in solar flares, for example, Alfvén wave heating (Reep & Russell 2016;Kerr et al. 2016) and heating by non-thermal proton beams (Procházka et al. 2018).It is likely that non-thermal protons are present in flares, and may even carry energy equivalent to that of the non-thermal electron distribution (Emslie et al. 2012).However, largely owing to poor constraints on the properties of the distribution they are often ignored in flare modelings, and numerical studies of their role in the Sun's atmospheric response to flares are rare.For further discussion see the introduction to Kerr et al. (2023).
In this work, we analyze single-loop RADYN proton beam simulations in a wide set of beam parameters and impose the pressure perturbations from these models into the 3D helioseismic model of sunquakes (Stefan & Kosovichev 2020).Our goal is to answer the question: can single-loop RADYN proton heating simulations, coupled with the 3D acoustic model of the Sun, cause a photospheric impact and explain the initiation and propagation of helioseismic signals detected in sunquakes?The paper is structured as follows.Section 2 describes RADYN proton beam heating simulations, synthesis of the Fe I 6173 Å line profiles and SDO/HMI line-of-sight (LOS) observables for these simulations, and the related results.Section 3 describes the development of the acoustic model simulations driven by the perturbations imposed from the RADYN models and the related results.The summary of the paper results is presented in Section 4 and is followed by a discussion in Section 5. 1 Hereafter when we refer to HMI's continuum intensity observable we drop the 'HMI' We employ the unified computational model for solar and stellar flares, RADYN (Allred et al. 2015), which selfconsistently combines the equations of radiation transport, non-equilibrium atomic level populations, charge conservation, and hydrodynamics following the injection of flare energy into one leg of a semi-circular loop, assumed to be symmetrical.RADYN was recently coupled with the FP code (Allred et al. 2020) which models the transport of high-energy particles from their injection in the solar corona to thermalization in the low atmosphere by solving the Fokker-Planck equation.The FP code models the evolution of the distribution function taking into account Coulomb collisions, nonuniform ionization, magnetic mirroring, the return electric currents, and synchrotron emission reactions, and offers an improvement upon prior implementations of non-thermal particles in RADYN.FP implements Coulomb collisions of non-thermal particles with an ambient plasma of arbitrary temperature and is applicable throughout the cold-target to warm-target regimes.Importantly for this study, for non-thermal protons with energies in the deka-to hecto-keV range, the ambient plasma is typically a warm-target, which is much less effective at slowing the protons than cold-target collisions (Tamres et al. 1986).
All RADYN proton beam runs analyzed in this work have approximately the same energy flux rates, F d = 10 11 erg cm −2 s −1 , with two values for the energy power-law spectral indexes: δ = 3 and 5, and eight values for the low-energy cutoff, E c = 25 keV, 50 keV, 100 keV, 150 keV, 250 keV, 500 keV, 1000 keV, and 3000 keV.The proton beam was injected into the atmosphere during the first 20 s of each simulation, after which the atmosphere evolved for an additional 80 s.Of these 16 total runs 15 are considered here.The run, parameterized by F d = 10 11 erg cm −2 s −1 , δ = 5, and E c = 25 keV, required a very small time step due to large perturbations of the atmosphere driven by beam heating, and was performed for a shorter duration.The overall properties of the completed runs are summarized in Tables 1 and 2. Note that the pre-flare atmosphere used in the F-CHROMA grid of simulations is rather different than that used here.Also, we inject energy flux at a faster rate than the F-CHROMA grid, which is injected more gradually over a 20 s triangular pulse, and the energy flux is twice larger on average.

Modeling of Fe I 6173 Å Stokes Profiles and SDO/HMI Observables
The synthesis of the Fe I 6173 Å Stokes Profiles for the RADYN simulations is similar to the procedure described in Sadykov et al. (2020).In brief, the NLTE Fe atomic level populations and resulting emission from bound-bound transitions are calculated under the assumption of statistical equilibrium using the RH1.5D radiative transfer code (Pereira & Uitenbroek 2015;Rybicki & Hummer 1991, 1992;Uitenbroek 2001).Since RH15D solves the statistical equilibrium equations, non-equilibrium effects are not included.The dynamic hydrogen populations and the free electron densities are imported from the RADYN model and not recalculated, somewhat mitigating the requirement to assume statistical equilibrium in the solution of Fe I.The full Stokes profiles and, correspondingly, the Right-Circular and Left-Circular Polarizations (RCP and LCP) are solved for the Fe I 6173 Å transition, with no effects of the background polarization taken into account.The models are augmented with various settings of the imposed vertical magnetic field unchanged during the simulations.Although we consider 500 G and 1000 G magnetic field settings (selected for illustration purposes only), the choice of the magnetic field does not impact the hydrodynamic simulation results, and only affects the calculations of the Stokes profiles and HMI observables.The Stokes profiles are computed from RADYN run output that is typically at a cadence of 0.1 s.
The properties of the Fe I 6173 Å line are computed (1) directly from the fully-resolved Stokes profiles, (2) by applying the SDO/HMI Line-of-Sight (LOS) pipeline to the Stokes profiles and assuming that the polarization signals are obtained instantly for every wavelength/filter, and (3) by applying the SDO/HMI Line-of-Sight (LOS) pipeline to the Stokes profiles following the proper temporal sequence of polarization measurements (Schou et al. 2012).For the properties directly estimated from the fully-resolved Stokes profiles, we compute the continuum intensities (calculated as the average intensity at ±1 Å from the reference wavelength of the line, λ ref =6173.34Å), line depths (defined as the strongest deviations from the continuum intensities), and line Doppler shifts (calculated using the center of gravity approach; e.g., Sadykov et al. 2019).The same observables are obtained using a simplified version of the HMI data analysis algorithm following the original method described by Couvidat et al. (2012Couvidat et al. ( , 2016)).The simplified pipeline is described in detail in Sadykov et al. (2020) and its application to RADYN electron beam heating runs is tested therein.Here, we provide a brief summary of the pipeline.
The LCP and RCP profiles of the Fe I 6173 Å are sampled in six wavelength points by assuming a Gaussian-like transmission profile in the wavelength space (with FWHM∼76 m Å).The wavelength points are centered at the rest wavelength of Fe I 6173 Å (λ ref =6173.34Å) and are sampled at 68.8 m Å apart.From these six measurements, the first and second Fourier components are computed and used for the estimation of the HMI observables (continuum intensity, line depth, Doppler velocity, and LOS magnetic field; see Couvidat et al. 2012) with the corresponding correction factors (Couvidat et al. 2016).The Fe I 6173 Å line profile is assumed to be Gaussian.In contrast with the previous work (Sadykov et al. 2020), the line width is not assumed to be fixed but recalculated for every application of the SDO/HMI pipeline.Keeping the line width fixed and equal to the unperturbed line profile width does not change the results and conclusions qualitatively.We sample the line profiles at each time snapshot and compute 'instantaneous observables' (i.e., observables where the temporal sequence of measurements is not taken into account).
The HMI LOS pipeline measurement sequence takes 45 s to be completed.Given that the heating phase lasts only 20 s, the observables will depend on the time of the heating with respect to the measurement sequence timing.In this work, for every time moment of the simulation, we assume that the SDO/HMI pipeline is centered temporally at this time moment, and compute the observables.Therefore the dynamics of the observable presented further needs to be interpreted as what 'can possibly be' observed by HMI during the proton beam heating event rather than what 'is' observed.That is, the best-case scenario is HMI happens to catch the start of the impulsive heating of a single pixel.The temporal sequence of polarizations and wavelengths are assumed following Schou et al. (2012).The polarization profiles for t < 0 s are assumed to be the same as the profiles of an unperturbed atmosphere at t = 0 s, and the same as for t = 100 s time moment for any t > 100 s.This allows one to compute the observables during the heating phase and for the last 20 s of the run.

Results from RADYN
The summary of the strongest perturbations (both the physical properties of the atmosphere and the properties of the Fe I 6173 Å line formation and appearance) for all considered RADYN runs is presented in Tables 1 and 2. The perturbations of the atmospheric parameters are computed for the heights of 100 km and 300 km which correspond to the heights where the Fe I 6173 Å line is typically formed for the quiet Sun conditions (Norton et al. 2006;Kitiashvili et al. 2015).We note here that during the flare process, the Fe I 6173 Å may experience a significant chromospheric contribution (Monson et al. 2021) which we do not investigate in this work.The height of 100 km is chosen for analysis instead of the 0 km because the latter is close to the bottom boundary of the modeling domain (which is just ∼60-65 km below the 0 km height) and affected by the boundary conditions.The selection of the upper height of h = 300 km is consistent with maximum values of τ λ = 1 modeled for the continuum near the Fe I 6173 Å line presented in the tables (although, as noted earlier, this does not preclude the possibility of the contribution of higher levels of the atmosphere to the Fe I 6173 Å formation).The atmospheric perturbations at h = 300 km for the presented runs range from several percent for the considered parameters (for moderate E c ∼ 250 keV values) to several tens of percent for the extremely high or extremely low values of the E c .The perturbations are mostly several percent at the height of h = 100 km for these runs and rarely reach tens of percent.Vertical velocities rarely reach as much as several hundreds of m s −1 at h = 100 km.At that height, the largest downward velocity2 is in Model 9 with v z ∼ − 0.794 km s −1 .However, vertical velocities are significantly faster at a height of h = 300 km, where, in Model 9, the peak downward velocity is v z ∼ − 2.78 km s −1 .The electron beam-induced velocities for the F-CHROMA grid were significantly lower (Monson et al. 2021;Sadykov et al. 2020).
For the atmospheric parameters (temperatures, pressures, and vertical velocities at the heights of h = 100 km and h = 300 km), the strongest perturbations occur for the lowest and highest values of the low-energy cutoff parameter, E c , corresponding to the cases of the highest number of protons or the highest energy per proton.
Tables 1 and 2 also present the spectroscopic line parameters as computed from the fully resolved Stokes profiles (i.e. with no application of the SDO/HMI LOS algorithm at this point).The strongest perturbations of the Fe I line profiles (in terms of their continuum intensities, line depths, and Doppler shifts) and changes to the τ = 1 heights occur in the cases of very high or very low E c .In particular, the continuum intensity enhancement reaches more than 40% and the strongest blueshift (representing an upflow) reaches almost v D ∼ 0.97 km s −1 for Model 8 (E c = 3000 keV, δ = 3), which has the highest average energy per proton considered.As mentioned previously, Model 9 (E c = 50 keV, δ = 5) has the strongest redshift of v D ∼ − 2.78 km s −1 at h = 300 km.
Figure 1 illustrates the maximum values of the relative continuum intensity enhancements, I c /I 0 c , and the strongest redshift, min v D (normally interpreted as a downflow speed, though opacity effects make a one-to-one relation difficult), as functions of the E c and δ parameters of the proton beam energy spectrum.While both the continuum intensity enhancement and strongest redshifts tend to increase with the change of the E c for both δ setups, the trends for I c /I 0 c and v D are slightly different.The enhancement of the continuum near the Fe I 6173 Å spectral line is in the range of 5-10% for most of the runs and increases consistently with the increase of E c .The situation is the opposite for Doppler shifts, where the strongest redshifts tend to increase with the decrease of E c .Models 8 and 9 have the strongest atmospheric and line profile perturbations so are chosen for more detailed analysis.Figure 2 displays the evolution of the temperature and pressure relative to their unperturbed values (at t = 0 s) and the vertical velocities for these two models.One can see that both models resulted in significant perturbations of the lower atmosphere.Figure 3 visualizes the behavior of temperature, vertical velocity, gas pressure, and the τ = 1 height of the Fe I 6173 Å line core and nearby continuum as a function of time.Model 8 has a stronger temperature response and a larger change in the τ = 1 heights than does Model 9.In Model 9, during the time interval of ≈30-90 s, there is extensive enhanced pressure in the region of line formation that is not evident in Model 8.In contrast to Figure 2e, Figure 2b shows the downward propagation of the pressure enhancement peak in the lower atmosphere, below 750 km.In contrast, for Model 8 the increase of the gas pressure happens during the first 20 s of the run during the heating phase.

Analysis of Stokes profiles
Figures 4 and 5 illustrate the evolution of the Fe I 6173 Å left-and right-circular polarization profiles (LCP and RCP) throughout the simulation.Overall, the line profiles presented in Figure 4 (corresponding to Model 9) do not experience notable changes until after the initial heating phase, when they become asymmetric: the blue wing of the profiles deepens between t ≈20-30 s and then enhances, with the strongest enhancement happening at around t∼60 s. Figure 3b demonstrates that the strongest enhancement of the blue wing follows the strongest downward motions of the atmosphere at h = 300 km by just several seconds.In contrast to that behavior, the dynamics in Model 8 result in more enhanced, complex shapes of the line profiles during the heating phase and right after it, as seen in Figure 5a-c.In particular, the Fe I line LCP and RCP profiles demonstrate an emission feature superimposed with the absorption and have an enhanced blue wing after 10 s of the start of the run.The enhancement of the Fe I polarization profiles is also reflected in τ = 1 heights: the heights lie in a more shallow region, within h ∼50-200 km (see Figure 3h).
In addition to properties derived directly from the LCP and RCP profiles, we have applied the simplified SDO/HMI pipeline (Sadykov et al. 2020, also described in Section 2.2) and computed the instantaneous and time-dependent lineof-sight (LOS) observables.The results are presented in Figure 6 for the case of a 500 G imposed vertical magnetic field, and do not change qualitatively for other values of the field.The instantaneous observables demonstrate the same patterns as the spectral line properties (continuum intensities near the line, line depths, and Doppler shifts) derived from the full-resolution line profiles for Models 8 and 9.The systematic offset of the synthetic observables with respect to the quantities derived from their full-resolution counterparts is typically a result of the non-Gaussian shape of the Fe I line (see Figures 4 and 5) while the SDO/HMI LOS pipeline relies on the Gaussian-shape assumption (Couvidat et al. 2016).The behavior of the time-dependent observables (taking into account the timing of the polarization measurements at different wavelengths) is close to the instantaneous observables for Model 9 (Figure 6a-d).However, for Model 8 the behavior of the line depth (Figure 6f) and Doppler velocity (Figure 6g) observables differ dramatically.In particular, while the redshifts measured from the fully resolved LCP and RCP profiles reach only ≈-0.7 km s −1 , the time-dependent SDO/HMI observables demonstrate Doppler shifts greater than -2 km s −1 .The vertical magnetic field observable fluctuates for this model approximately in the range of 300-900 G, while its true value was always unchanged and equal to 500 G.Such dynamics are due to the strong impulsive heating whose duration (20 s) is less than the duration of the SDO/HMI LOS pipeline procedure (45 s).Model 9 (Fig. 6a-d) shows a much better agreement between the properties of the line profile derived directly from the spectrum and the SDO/HMI observables with respect to Model 8 (Fig. 6e-h), most likely because of the more gradual evolution of the Fe I polarization profiles evident in Figures 4 and 5. Overall, it once again confirms that SDO/HMI LOS observables, including a continuum intensity observable, have to be interpreted with caution during solar flares ( Švanda et al. 2018).

3D ACOUSTIC MODELS DRIVEN BY RADYN SIMULATIONS
In this section, we use time-dependent accelerations derived from the proton beam simulations to excite acoustic waves in a whole-Sun acoustic model with solar background stratification.The total duration of all beam inputs to the acoustic model is 20 seconds.An example of the input acceleration profiles, as described in the next section, for low-energy cut-offs E c = 50 keV and E c = 3000 keV (Models 9 and 15) interpolated onto the hydrodynamic model's background mesh is shown in Figure 7,

Model Description
The 3D acoustic model treats acoustic oscillations as linear, adiabatic perturbations to pressure, density, and velocity (Stefan & Kosovichev 2020).The background stratification for the solar interior is derived from the Standard Solar Model (Christensen-Dalsgaard et al. 1996) which smoothly transitions to the atmosphere used in the RADYN simulations.The transition location in the Standard Solar Model, at R = 695.707Mm relative to the solar center, is chosen where the mass density is equal to that at z = 0 Mm in the RADYN mesh.The difference between the original photosphere (R = 695.996Mm) and the new photosphere corresponds to a Wilson depression of z W ≈ 300 km, consistent with the expected photospheric depression in the presence of active region magnetic fields.
The acoustic model itself is semi-spectral, with the radial derivatives evaluated numerically and polar-and azimuthalangle derivatives evaluated spectrally using spherical harmonics.While the acoustic model does not include radiative damping, the generated oscillations are damped according to the horizontal wavenumber, with damping parameters derived from quiet-Sun pressure wave (p-mode) data reported by E J Rhodes et al. (2011).We consider the choice of quiet-Sun damping parameters, as opposed to active region-like, appropriate here as the majority of the sunquake wavefront propagates outside the generating active region, where the magnetic field is moderate or weak.We take advantage of the advanced treatment of radiation in RADYN by deriving input accelerations for the acoustic model using the simulated perturbations to gas pressure from the applied proton beam heating.These gas pressure perturbations are smaller in RADYN than would be generated in the acoustic model as RADYN accounts for the energy lost in optically thin and NLTE optically thick radiation.The input accelerations are computed from the radial gas pressure gradients in the RADYN simulations, vr = (1/ρ 0 ) • ∂P /∂r, where ρ 0 is the initial state's mass density.
We assume a constant cross-sectional area, and the horizontal profile of the input accelerations is considered to be Gaussian with an FWHM of 1.5 Mm.The FWHM is based on HMI observations of sunquake kernels that range in size from one to several pixels; this corresponds to an impact site between 0.75 Mm and 2 Mm across.Where the lower end of the supplied accelerations ends, slightly below R = 695.707Mm, the functions are appended by a Gaussian with a drop-off closely matching the unappended input.A similar drop-off is applied to the upper end of the supplied accelerations, which extend to the top of the modeled corona, to avoid boundary effects that may be caused by providing input close to the upper boundary.The Gaussian upper drop-off begins 350 km from the upper boundary with an FWHM of 125 km.An additional description of the acoustic model, including numerics, is provided in Stefan & Kosovichev (2020).

Results from the Acoustic Model
We are primarily interested in the behavior of the photospheric radial velocity, as this is generally the largest component of the line-of-sight velocity in SDO HMI Dopplergrams for observations close to the disk center.We then examine the resulting photospheric p-mode wavefront, with absolute maximum radial velocity for each case shown as a function of horizontal distance from the beam target in Figure 8.Note that acoustic-gravity waves-not typically observed in actual sunquake events-are generated in addition to the usual p-mode wavefront.In the cases where the cut-off energy is greater than 250 keV, these acoustic-gravity waves have amplitudes that exceed the p-mode wavefront.We, therefore, examine the absolute maximum amplitude at each distance within 5 minutes of the wavefront travel time predicted by ray theory; however, it is not possible to disentangle the two wavefronts for short (<5-7 Mm) distances.
We observe core velocities in the acoustic model which are significantly greater than those in the RADYN simulations, in particular for the simulations with low cut-off energies.For example, the acoustic model predicts a magnitude of the radial velocity for Model 9 (E c = 50 keV, δ = 5) of 68.7 km s −1 ; in such cases, we do not expect the predictions that are close to the source to be reliable because of both the linear nature of the model and the lack of radiative damping.Conversely, the magnitude of the radial velocity for Model 15 (E c = 3000 keV, δ = 5) at the beam core is only 751 m s −1 which securely falls in the linear regime.Thus we consider the entire range of the acoustic model to be reliable for the higher cutoff-energy cases.
The relationship between the absolute maximum radial velocities of each case remains largely the same over distance as compared to the beam core.In general, the greatest absolute velocities are associated with small low-energy cut-offs with the wavefront amplitudes decreasing with increasing low-energy cut-off.There is relatively little deviation in this relationship when increasing the spectral index from δ = 3 to δ = 5, though we note that the radial velocities for the E c = 500 keV cases do change appreciably from the δ = 3 case (Model 6) to δ = 5 case (Model 13).In the associated acceleration profile from each model, we observe significantly stronger evaporation in the δ = 5 case as well as a downward-propagating acceleration front that penetrates slightly more deeply in Model 13 than in Model 6.While the initial evaporation in each case is similarly impulsive-that is, the evaporation fronts propagate upwards with similar speeds-more of the beam energy is released in the initial evaporation of Model 13.For models with lower energy cut-offs, there is a similar increase in acceleration magnitude when moving from δ = 3 to δ = 5.The downward-propagating acceleration front is weak for low-energy cut-off E c = 1000 keV (Models 7 and 14) and nonexistent for low-energy-cutoff E c = 3000 keV (Models 8 and 15), and the radial velocity in these two cases changes the least with spectral index.The lack of the downward-propagating front in Model 15 is clearly seen in the comparison with Model 9 in Figure 7.
We now look more closely at the maximum radial velocity at a horizontal distance of X = 18 Mm indicated by the vertical dashed line in Figure 8, where the sunquake wavefront is expected to reach its maximum amplitude (aside from the beam core).Explicitly plotting these velocities on a log-log scale, as in Figure 9, we find that the generated sunquakes fall into two separate regimes.There is a low-energy cut-off regime extending up to E c = 250 keV where the sunquake amplitudes are similar to observations (on the order of 100s of m s −1 ), and a high-energy cut-off regime beginning at E c = 1000 keV where the sunquake amplitude is significantly lower than in observations.

SUMMARY OF THE RESULTS
To conclude the analysis of RADYN models and synthetic Fe I 6173 Å emissions and SDO/HMI LOS observables, we claim that there are two regimes found in which the perturbations of the line profiles and the atmosphere were significant and resulting in a potential helioseismic response or a white-light flare (see Figure 1).We have selected the two most extreme models.Model 9 (E c = 50 keV, δ = 5, F d = 10 11 erg cm −2 s −1 ) has the highest rate of the momentum deposition into the atmosphere (except the model E c = 25 keV, δ = 3, F d = 10 11 erg cm −2 s −1 which, interestingly, led to weaker perturbations).The model causes gradual but strong changes in the atmospheric parameters, resulting in ≈12% increase in flare continuum during the heating phase.The Doppler velocity derived from well-resolved polarization profiles has the strongest redshift of v D ∼ − 0.8 km s −1 about 40 s after the heating.The SDO/HMI observables, in general, closely follow the properties of the full-resolution line profile.The other model, Model 8 (E c = 3000 keV, δ = 3, F d = 10 11 erg cm −2 s −1 ), has the highest average energy per proton.It causes a much stronger enhancement of the continuum near the Fe I 6173 Å line (of the order of 42%) with respect to Model 9 and also results in the Fe I 6173 Å redshifts of v D ∼ − 0.7 km s −1 at the end of the heating phase.The SDO/HMI observables in this model differ significantly from the counterparts derived from the full-resolution line profiles.Overall, both models demonstrate the impact on the deep layers of the solar atmosphere, and the spectral analysis of the Fe I 6173 Å line does not yet allow us to claim whether the proton beams of high energy per proton or low energy per proton are the preferential candidates for causing sunquakes.
To provide more insights, one can utilize the responses of the atmosphere to the proton beam heating and impose this response into the 3D acoustic models (Section 3).While the Doppler velocities in the RADYN models do not differ much at either extreme of the low-energy cutoff, the physical velocities in the acoustic model are significantly stronger for Model 9 (E c = 50 keV, δ = 5, F d = 10 11 erg cm −2 s −1 ) than for Model 8.The continuum intensity enhancement derived from the RADYN simulation is, on the opposite, stronger for the δ = 3 cases, with Model 8 producing the strongest enhancement.There appear to be two separate regimes in the low-energy cutoff spectrum: small low-energy cutoff beams which produce a strong seismic signal and large low-energy cutoff beams which produce comparatively more intense continuum emission.Observational evidence for such a duality can be found in, for example, the analysis of Pedram & Matthews (2012) where white light enhancement for a sample of nine flares was found to be noticeably weaker in sunquake-generating events.However, both the RADYN model and acoustic model show an increase in the Doppler velocity and physical velocity, respectively, in the E c = 3000 keV beam when the spectral index increases from δ = 3 to δ = 5.It is possible that such large low-energy cutoffs can explain sunquake-generating events with significant white-light enhancement, though an extension of this analysis to higher low-energy cutoffs and spectral indices is necessary to verify the trend.Overall, we show not only that sunquakes can be generated from a high injected flux of protons with relatively low cutoff energy, but that these sunquakes have much higher amplitudes than those predicted from simulations with large cutoff energies.This is particularly significant as large low-energy cutoffs are more difficult to physically justify.

DISCUSSION
The role of the proton beams in the solar flare energy deposition has been discussed for more than five decades ( Švestka 1970;Simnett 1986), yet without a clear understanding of the energy fraction carried by these beams.Some works (Emslie et al. 2012) suggest that the energy transported by the proton beams in solar flares is comparable to that of the electron beams.Proton beams are an attractive candidate as a mechanism to explain sunquake excitation as they deposit energy significantly deeper in the solar atmosphere than electron beams.They also carry more momentum with respect to the electron beams of similar energy due to the difference in the particle masses.Given that the energy spectrum of the proton beam is determined by the power law, with the constant A found by normalizing to the known deposited energy rate, the total momentum flux of the beam (in the non-relativistic limit) and the average energy per particle are then and Although we do not analyze the momentum transport in detail in this work, it is helpful to provide some estimates.The injected momentum flux for the proton beam in Model 9 (with F d = 10 11 erg cm −2 s −1 , δ = 5, and E c = 50 keV) is M d ∼ 5.5 × 10 2 g cm −1 s −2 .The total momentum deposited per unit area during 20 s therefore is M tot d = M d ∆t ∼ 1.1 × 10 4 g cm −1 s −1 .Assuming the photospheric density ρ ≈ 10 −7 − 10 −8 g cm −3 (considering the quiet Sun and sunspot model atmospheres, correspondingly; Allred et al. 2015) and the characteristic scale height of H ≈ 100 km, the bulk velocity of the plasma is v ≈ M tot d /(ρH) ≈ 0.1 − 1 km s −1 .The latter velocity estimate is comparable to what is expected for sunquakes, and we can therefore conclude that the momentum transport by the proton beams may play an important role in sunquake initiation.
We have explored whether proton beams of the deposited energy flux (F d = 10 11 erg cm −2 s −1 ) are capable of exciting sunquakes with amplitudes similar to observations.For fixed cut-off energy, we find that the beam with spectral index δ = 5 produces a higher amplitude wavefront than the corresponding δ = 3 beam in nearly every case; the only outlier here is E c = 3000 keV, though the difference is at most only 0.5 m s −1 .This is generally consistent with Eqn. 3 which shows that, given the same E c and F d , larger values of δ result in larger momenta.Also, the momentum flux for fixed cut-off energy varies weakly with the spectral index, M d ∝ (δ − 2)/(δ − 3 2 ), and is only significant for smaller values of δ.For example, a proton beam with δ = 7 is expected to deposit only 6% more momentum than the same beam with δ = 5.However, for the considered δ = 3 and δ = 5, the difference in M d for the beams with the same E c is nearly 29%.The momentum is also inversely proportional to the √ E c which explains the increasing trend of the Doppler velocities in Fig. 1b.Decreasing the cut-off energy has a similar effect on the energy spectrum of the proton beam as increasing the spectral index, though the momentum flux depends more strongly on the cut-off energy, M ∝ 1/ √ E c .This dependence is only noticeable for E c ≥ 250 keV; for cut-off energies less than this, the decrease in sunquake amplitude with increasing cut-off energy is much milder.
The F-CHROMA models considered in Sadykov et al. (2020) have a peak energy deposition rate of F d = 10 11 erg cm −2 s −1 .However, the average energy deposition rate was twice lower than that value, which may impact the close comparison of F-CHROMA models and the RADYN proton beam heating models considered in this work.Figure 10 illustrates a comparison of the response of the atmospheric parameters at h = 300 km for the proton beam model F d = 10 11 erg cm −2 s −1 , E c = 50 keV, δ = 5 (Model 9) and the electron beam heating model with F d = 10 11 erg cm −2 s −1 , E c = 15 keV, δ = 5 from Graham et al. (2020).The RADYN electron beam heating model had a duration of heating of 20 s similar to the one in proton models, and a total duration of the run of 60 s.One can see in Figure 10 that, while experiencing weaker but yet comparable temperature enhancement, the electron beam heating model does not result in strong downward velocities and pressure enhancements at h = 300 km.
The correlation of the white light emission (including SDO/HMI observable continuum) with the hard X-ray sources in solar flares (Watanabe & Imada 2020;Battaglia & Kontar 2012) suggests that the electron beam heating and radiative back-warming can enhance the white light emission.Fig. 1a demonstrates that the proton beams can contribute to the enhancements of the continuum near the Fe I 6173 Å line as well.The Figure also shows a tendency of the white-light enhancements to increase with the E c increase and preferring lower δ values.This is in qualitative agreement with Eqn. 4 which demonstrates that (1) for the same E c the δ = 3 beam particles have, on average, 50% more energy per particle, and (2) the < E p > is directly proportional to E c .The higher values of < E p > would indicate that the proton beam particles, on average, have to encounter thicker media to lose their energy through thermalization and, therefore, have to penetrate deeper into the solar atmosphere.All of the considered proton beam models exceeded this enhancement: 3.9% was the weakest enhancement value observed for Models 11 and 12, Models 6-9 and 15 had an enhancement of more than 10%, and the remaining Model enhancements were concentrated within the 6%-10% range.Taking into account these numbers, the continuum enhancements of more than 10% presented in several works (Procházka et al. 2018;Song & Tian 2018) can be contributed by proton beams.It is also important to notice that the enhancements reported in this work are derived directly from the computed continua near the Fe I 6173 Å line and not after applying the SDO/HMI observable algorithm, which can potentially generate artificial enhancements of the continua (Mravcová & Švanda 2017;Švanda et al. 2018).
A population of the very high-energy (>30 MeV) protons can be diagnosed by producing the 2.223 MeV neutroncapture γ-ray line (Shih et al. 2009).However, the lower-energy proton distribution is currently almost 'invisible' to the observer.Recently, Kerr et al. (2023) studied the Orrall-Zirker effect (Orrall & Zirker 1976) by performing proton beam driven RADYN simulations, with F d =10 9 -10 11 erg cm −2 s −1 , E c = 150 keV, δ = 5.Though their models predicted a much weaker, and very transient, signal than earlier experiments suggested, Kerr et al. (2023) did find a detectable non-thermal enhancement of Lyman lines produced via a charge exchange between the protons in the beam and the ambient plasma and indicated the potential possibility to diagnose the injected beam via Lyβ observations from the Spectral Imaging of the Coronal Environment (SPICE) Instrument onboard the Solar Orbiter. Figure 4 of (Orrall & Zirker 1976) also indicates that the enhancements of the Lyα wings are most notable for lower δ values and for the ≈30 keV protons (also noticed in Simnett 1995).This provides an opportunity to observationally test the modeling-based selection rule that lower-energy proton beams are responsible for sunquakes, especially given that many strong flares have helioseismic counterparts (Sharykin & Kosovichev 2020).

2.
MODELING FE I 6173 Å STOKES PROFILES AND SDO/HMI OBSERVABLES FROM RADYN PROTON BEAM SIMULATIONS 2.1.Description of RADYN Proton Beam Heating Runs 6173 Å line formation for RADYN proton beam simulations with δ = 3, F d = 10 11 erg cm −2 s −1 .The subscript '0' denotes the parameter value at t=0 s.The 'min' and 'max' operators return the minimum and the maximum values of the parameters during the run.τ Ic denotes the optical depth for the continuum near the Fe I

Figure 1 .
Figure 1.The maximum intensity enhancement (a) and the strongest Doppler redshift (b) as a function of the low-energy cutoff, Ec, for the considered models.The visualization is based on the RADYN model summaries in Table1 and 2.

Figure 6 .
Figure 6.The continuum intensity (a), line depth (b), Doppler velocity (c), and magnetic field (d) inferred from the synthetic Fe I 6173 Å line profiles from Model 9 with an imposed vertical uniform 500 G magnetic field.The black dashed curves correspond to measurements inferred from the RH-predicted line profiles.The red solid curves show "instantaneous" observables obtained with the HMI algorithm applied to instantaneous line profiles.The blue curves show the observables obtained with the HMI algorithm applied with the actual observing sequence timing centered at the reference time.The black dashed horizontal lines mark the zero level of the observables.Panels e-h illustrate the same quantities for Model 8.

Figure 7 .
Figure 7. Time-dependent acceleration profiles derived from the RADYN proton beam simulations for Models 9 (left) and 15 (right).

Figure 8 .
Figure 8. Log of the respective sunquake wavefront amplitudes as a function of distance from the excitation source.The typical p-mode wavefront is isolated from other waves generated in the acoustic model by considering the maximum amplitude only within five minutes of the predicted wavefront arrival time.The sample of velocities presented in Figure 9 is marked by the vertical dashed line.

Figure 10 .
Figure 10.Illustration of the temperature (a), vertical velocity (b), and gas pressure (c) as a function of time at the height of 300 km for Model 9 (black line) and for a RADYN simulation of non-thermal electron beam heating characterized by a power-law with Ec = 15 keV, δ = 5, F d = 10 11 erg cm −2 s −1 (considered in Graham et al. 2020, red line).

Table 1 .
Properties of the atmosphere and the Fe I

Table 2 .
Same as Table 1 but for the models with δ = 5, F d = 10 11