Numerical Prediction of Changes in Atmospheric Chemical Compositions During a Solar Energetic Particle Event on Mars

Precipitation of solar energetic particles (SEPs) into planetary atmospheres causes changes in atmospheric chemical composition through ionization, dissociation, and excitation of atmospheric molecules. In contrast to the terrestrial atmosphere, where depletion of ozone in the polar mesosphere has been studied by observations and models for decades, there have been no studies on the effects of SEPs on the neutral chemical composition of Mars' present‐day atmosphere. This study provided the first estimate of the impacts of SEPs on neutral chemical composition in the present‐day Martian atmosphere coupling a Monte Carlo model and a one‐dimensional photochemical model. Our results showed that ozone density in the Martian atmosphere might decrease in the altitude range of 20–60 km with a factor 10 maximum enhancement occurring at 40 km during a Halloween‐class SEP event due to an enhanced concentration of HOx. The depletion of ozone occurred in the altitude range of 20–60 km, corresponding to the penetration of protons with an energy range of 4.6–46 MeV. Such ozone depletion should be detected by Trace Gas Orbiter/Nadir and Occultation for MArs Discovery during intense SEP event. A 75% depletion of the ozone density at 40 km can be expected during SEP events occurring once in every 1 year. Therefore, ozone concentration in the Martian atmosphere might sufficiently decrease during a SEP event as on Earth, but through different chemical pathways driven by CO2 ionization and CO recombination catalytic cycle.


Introduction
Solar energetic particles (SEPs) ionize, dissociate, and excite the atmospheric molecules when they precipitate into a planetary atmosphere, leading to changes in ion and neutral chemical compositions.The effects of SEPs on the chemical composition in the terrestrial atmosphere have been intensively studied for decades.For instance, the SEP events of late October and early November 2003, known as the Halloween event, induced a depletion of the ozone density by up to 40% and an enhancement of the odd hydrogen (HOx) and the odd nitrogen (NOx) in the polar mesosphere and stratosphere (e.g., Jackman et al., 2005;Randall et al., 2005;Seppälä et al., 2004).Precipitating energetic particles dissociate nitrogen molecules producing N and N( 2 D) in the atmosphere, the latter of which reacts with O 2 becoming NO (Crutzen et al., 1975;Rusch et al., 1981).Precipitating energetic particles ionize O 2 producing O 2 + in the atmosphere, which reacts with ambient water vapor to produce water cluster ions, and thus H and OH via recombination of the cluster ions with electrons (Solomon et al., 1981).The produced HOx and NOx catalytically destroy ozone in the polar mesosphere during SEP events.Such an effect was confirmed by the observed anti-correlation between HOx/NOx and ozone concentrations during SEP events (e.g., Crutzen et al., 1975;Jackman et al., 2005;Seppälä et al., 2004;Verkhoglyadova et al., 2015).
The impacts of SEPs on the Martian atmosphere are not local as in the Earth's atmosphere but global due to the absence of a strong global magnetic field on Mars.The imaging ultraviolet spectrograph (IUVS) instrument on board the Mars Atmosphere and Volatile EvolutioN (MAVEN) spacecraft has observed ultraviolet diffuse auroral emissions spanned across the whole night side of Mars during SEP events (Schneider et al., 2015(Schneider et al., , 2018)).Recent model studies revealed that the diffuse aurora is caused by the precipitation of SEP electrons (Gérard et al., 2017;Haider & Masoom, 2019;Schneider et al., 2015) and SEP protons (Nakamura, Terada, Leblanc, et al., 2022).Nakamura, Terada, Leblanc, et al. (2022) succeeded in reproducing the observed low altitude peak of the diffuse Supporting Information: Supporting Information may be found in the online version of this article.
10.1029/2022JA031250 2 of 30 auroral emission profiles by taking into account the contribution of SEP protons, indicating that SEP protons have a significant impact on the lower atmosphere of Mars.
In contrast to Earth, there have been no studies of the effects of SEPs on the neutral chemical composition of the Mars' present-day atmosphere.In this study, we used a Monte Carlo model Particle TRansport In Planetary atmospheres (PTRIP) developed by Nakamura, Terada, Leblanc, et al. (2022) to calculate the vertical profiles of the ionization and dissociation rates of atmospheric molecules during a SEP event.We also used a one-dimensional photochemical model developed by Nakamura, Terada, Tao, et al. (2022) and Nakamura et al. (2023) to investigate the changes in atmospheric chemical composition during the SEP event.Our target species are ozone, HOx, and NOx, since the high-resolution spectroscopy Nadir and Occultation for MArs Discovery (NOMAD) instrument on board the Trace Gas Orbiter (TGO) spacecraft has a high sensitivity to these three species.We also described the formaldehyde (H 2 CO) and nitrous oxide (N 2 O), which are of astrobiological interest because they are important for prebiotic synthesis of amino acid and nucleic acid (e.g., Airapetian et al., 2016;Lingam et al., 2018).

Particle TRansport In Planetary Atmospheres (PTRIP)
PTRIP is a Monte Carlo model that solves the transport of energetic electrons, protons and hydrogen atoms in planetary atmospheres and calculate the ionization, dissociation and excitation rates of atmospheric molecules (Nakamura, Terada, Leblanc, et al., 2022).PTRIP solves the equation of motion for each incident particle and takes into account all possible collisions between each incident particle and the atmospheric molecules.Collision type, energy loss, scattering angle, and energy of the secondary electrons are randomly determined at each timestep using cross sections, scattering angle distribution, and secondary electron energy spectrum.In this study, only SEP protons are injected into the Martian atmosphere because SEP protons have been shown to be the dominant ionization source at low altitudes (Nakamura, Terada, Leblanc, et al., 2022).One thousand particles were injected for each incident energy with incident angles isotropically distributed over the downward hemisphere.Electric fields and magnetic fields were ignored because the motion of protons with an energy greater than MeV is not affected by the electromagnetic environment of Mars.
The atmospheric density profiles used in this study are shown in Figure 1.The vertical profile of the neutral temperature was taken from Chaffin et al. (2017) standard case with an exobase temperature of 240 K and a surface temperature of 210 K.The number density of CO 2 on the surface corresponds to a surface pressure of 6.7 mbar.The volume mixing ratio of N 2 at the surface was set to 1.9% (Mahaffy et al., 2013).The vertical profiles of the CO 2 and N 2 densities were calculated by considering vertical eddy diffusion and binary diffusion using a photochemical model PROTEUS to be explained in Section 2.2 with boundary conditions on the surface number densities.As for the vertical profile of water vapor, the relative humidity below the tropopause was fixed at 22% and the same volume mixing ratio of water vapor was used above the tropopause as did in Koyama et al. (2021).Using the profiles of CO 2 , N 2 , and H 2 O as the initial ones, the steady state profiles of CO, O 2 , and As the incident proton energy increases especially above 10 MeV, the number of secondary electrons produced in the atmosphere via ionizing collisions becomes massive.To reduce the computational cost when taking into account the contribution of the secondary electrons, we have developed a method to calculate the normalized secondary electron flux in energy-altitude (E-z) grids, which is then multiplied by the number of secondary electrons produced within each E-z grid cell.This method precisely obtains the same ionization profiles with much less computational cost than the simulation tracing all the secondary electrons (see Figure S1).
The accuracy of the proton-impact ionization and elastic cross sections determine the accuracy of the simulation.All the collisional cross sections for protons, hydrogen atoms, and secondary electrons in this study are the same as in Nakamura, Terada, Leblanc, et al. (2022).The analytic equation of proton-impact ionization cross sections of CO 2 and N 2 fitted by Rudd et al. (1983) based on the experimental data within the energy range from 5 keV to 4 MeV was applied up to 220 MeV, and the theoretical screened Rutherford elastic cross section, as in Nakamura, Terada, Leblanc, et al. (2022).Figure 3a shows the production rate of CO 2 + at each incident proton energy using the incident proton flux displayed in Figure 2. PTRIP predicted that protons with energy below 100 MeV lose all their energy before reaching the surface, whereas protons with energy larger than 150 MeV could reach the surface.Such result is consistent with previous simulations of SEP proton penetration in the Martian atmosphere (Guo et al., 2018) performed with the GEANT4-based Monte Carlo model PLANETCOSMICS (Desorgher et al., 2006).Comparison with PLANETCOSMICS showed that the cross sections chosen in this paper are valid at high energy so that PTRIP is able to solve the transport of protons up to 220 MeV in the Martian atmosphere.The calculated vertical profiles of the production rates of CO 2 + , CO + , O + , C + , O 2 + , N 2 + , N + , N( 4 S), and N( 2 D) are shown in Figure 3b.It is noted that the bumpy structure seen in the production rate profiles stems from the selection of discrete energy bins.Contrary to the terrestrial atmosphere, most of the proton energy leads to the ionization of the Martian CO 2 molecules and the production rate of N 2 + is two orders of magnitudes smaller than that of CO 2 + .The production rates of CO 2 + , CO + , O + , and C + by the proton-CO 2 collisions were estimated by multiplying the ionization rate of CO 2 with the ratios 0.75, 0.07, 0.13, and 0.05 based on the measurements of proton-impact dissociative ionization cross sections, respectively (Knudsen et al., 1995; Molina-Cuberos  et al., 2001).It is noted that the production rates of CO + , O + , and O 2 + by the collisions between protons and CO, O, and O 2 were also calculated, but the contributions to the production rates of the former two ion species are negligibly small.The production rate of atomic nitrogen N( 4 S) and N( 2 D), which are important precursors for NOx chemistry, were estimated by multiplying the production rate of N 2 + with the following relative abundance, N 2 + :N + :N( 4 S):N( 2 D) = 1:0.22:0.73:0.95based on the 200 eV electron impact cross sections which was used as a proxy for the production by high energy particle in Krasnopolsky (2009).

Photochemical and RadiatiOn Transport Model for Extensive USe (PROTEUS)
We used a one-dimensional photochemical model named Photochemical and RadiatiOn Transport model for Extensive USe (PROTEUS), which has been developed by Nakamura, Terada, Tao, et al. (2022) and Nakamura et al. (2023), which was designed for adaptability to many planetary atmospheres, for example, Jovian ionosphere (Nakamura, Terada, Tao, et al., 2022) and Martian atmosphere (T.Yoshida et al., 2023).PROTEUS solves a system of continuity equations: where n i is the number density of i th species, P i is the production rate of i th species, L i is the loss rate of i th species, z is the altitude, and Φ i is the vertical flux of i th species.The vertical flux Φ i can be described as follows.
where D i and K are the binary and eddy diffusion coefficient, respectively, H i is the scale height of i th species, H is the mean scale height of the atmosphere, q i is the charge of i th species, e is the elementary charge, T e and T i are the temperatures of the electron and of the i th species, respectively, T is the neutral temperature, P e = n e kT e is the electron pressure, with k is the Boltzmann constant, and α i is the thermal diffusion coefficient.The temperature profiles of the electrons and all of the other species were assumed to be the same as the neutral temperature.
The vertical profiles of the binary diffusion coefficient were taken from Hunten (1973) and the vertical profile of the eddy diffusion was taken from Krasnopolsky (1993).The thermal diffusion coefficient for H and H 2 were taken from Hunten (1973) and those for other species were set to zero.
We have implemented 503 chemical reactions for 34 neutral and 48 charged species into PROTEUS; ionization rate by galactic cosmic ray was taken from Haider et al. (2009), the C-, H-, and O-bearing neutral chemistry was taken from Chaffin et al. (2017), chemical reactions related to formaldehyde were taken from Pinto et al. (1980), the N-bearing neutral chemistry was taken from Nair et al. (1994), the ionospheric chemistry was taken from Fox and Sung (2001), the CO 2 -bearing cluster ion chemistry was taken from Molina-Cuberos et al. (2001), and the water cluster ion chemistry was taken from Verronen et al. (2016).The reaction rate coefficient for the reaction N( 4 S) + CO 2 → NO + CO is still unknown and the reaction speed is known to be very slow since this reaction does not conserve spin angular momentum (Rawlins & Kaufman, 1976).Fox and Sung (2001) used rate coefficient for this reaction equal to 1.7 × 10 −16 cm 3 s −1 as a standard model in the Venusian atmosphere, which was estimated as an upper limit by Brown and Winkler (1970) and Herron and Huie (1968).There are several experimental studies for the upper limit of this reaction rate, however, they are controversial.Avramenko and Krasnen'kov (1967) obtained a temperature-dependent rate coefficient for this reaction as 3.2 × 10 −13 exp(−1,711/T) cm 3 s −1 in the temperature range 291-523 K, which gives a rate coefficient of 1.0 × 10 −15 cm 3 s −1 at a temperature of 300 K. Rawlins and Kaufman (1976) estimated an upper limit of 1.0 × 10 −19 cm 3 s −1 at temperature 300 K. Fernandez et al. (1998) tried to measure the rate coefficient of this reaction, however, they did not observe any reactions and they estimated an upper limit of 1.1 × 10 −17 cm 3 s −1 at 285 K.In this study, the rate coefficient for this reaction was set to 1.0 × 10 −19 cm 3 s −1 estimated as upper limit by Rawlins and Kaufman (1976).The chemical reactions used in this study are listed in Table A1 in Appendix A. The upper boundary of the photochemical model was set to 200 km altitude and the lower boundary was the surface.At the upper boundary, the escape fluxes of H and H 2 were calculated as Jeans escape and that of O was fixed at 1.2 × 10 8 cm −2 s −1 (Chaffin et al., 2017).
The solar extreme ultraviolet (EUV) flux was calculated by the EUVAC model (Richards et al., 1994) with a F10.7 value of 140 for a moderate solar condition.The solar flux between 0.5 and 1,000 nm used to calculate the dissociation rates of atmospheric molecules was taken from Woods et al. (2009).The model takes into account the radiative transfer and the absorption by atmospheric species.The detailed information about the absorption and dissociation cross sections are given in Nakamura et al. (2023).All the simulations in this study were performed at the subsolar point.
In this study, we first ran the model for only neutral species without SEP inputs in order to obtain a converged solution for the neutral species.Since the chemical and transport time scales of ions are around several hours (e.g., Cravens et al., 2017), we then ran the model for both neutral and ionized species without SEP inputs for 10 simulated days in order to reach a quasi-steady state density profiles of ions.The quasi-steady state density profiles for all the species were used as initial profiles to investigate the SEP effects.The production rates of CO 2 + , CO + , O + , C + , O 2 + , N 2 + , N + , N( 4 S), and N( 2 D) were calculated by PTRIP and used as inputs into the photochemical model.We ran the model with SEP inputs for 1 simulated day (24 hr) to investigate the change of density profiles during the SEP event.The SEP event was assumed to last 24 hr and the proton flux intensity and its spectral shape were assumed to be constant during that duration.After that duration, we ran the photochemical model for 10 simulated days (240 hr) assuming no SEP inputs in order to examine the recovery phase.

Changes in the Ion Compositions
The ion density profiles in the quasi-steady state calculated by the photochemical model without SEP input are shown in Figure 4a.There are two ionospheric peaks, one formed by the irradiation of solar EUV flux above 75 km altitude and a second ionospheric peak formed by galactic cosmic rays below 75 km altitude.The dominant ion is O 2 + above 75 km while water cluster ions dominate below 75 km.Below 50 km altitude, the densities of Figure 4b displays the calculated ion density profiles 1 day after the onset of the SEP event.Ion and electron densities are significantly increased below 100 km altitude, while there is no change above 100 km.It should be noted that the contribution of the SEP electrons was ignored in this study, which is the dominant ionization source at higher altitudes during SEP events (Nakamura, Terada, Leblanc, et al., 2022).Below 100 km altitude, O 2 + is the dominant ion in the altitude range from 50 to 100 km and water cluster ions are the dominant ions below 50 km.The electron density below 75 km is enhanced by 2-3 orders of magnitudes and it reaches 10 5 cm −3 in the altitude range from 50 to 75 km.Increase in electron density at low altitudes during SEP events could lead to absorption of radio emissions from spacecraft and cause the disappearance of surface echo of radar instrument (e.g., Espley et al., 2007;Harada et al., 2018;Lester et al., 2022;Morgan et al., 2006;Sánchez-Cano et al., 2019;Sheel et al., 2012).Our results showed that the contribution of SEP protons to the ion profiles is greater than that of galactic cosmic rays during a large SEP event, which is consistent with the previous model conclusion regarding the SEP effects on the ionosphere of Mars (Sheel et al., 2012).
The temporal variation of the electron density profile during the SEP event is shown in Figure 5a and after the SEP event in Figure 5b.The time scales of increase and decrease of the electron density during and after the SEP event are both about 10-100 s and the electron density profile recovers to the pre-SEP state in 1,000 s after the end of the SEP event.Chemical loss time scale of electrons can explain this temporal variation.The reaction rate coefficient for the recombination of O 2 + with electron is 7 × 10 −8 cm 3 s −1 and that of water cluster ions is about 5 × 10 −6 cm 3 s −1 assuming that the electron temperature is the same as the neutral temperature of about 130 K. Assuming that electron density is 10 5 cm −3 , the loss time scale of electron is about 140 s above 50 km at which altitude O 2 + dominates the ion density profiles and is 2 s below 50 km at which altitude water cluster ions dominates the ion density profiles.Our results show that the electron density is sensitive to the temporal variation of the SEP proton flux due to the short lifetime of ions and electrons at such low altitudes.

Changes in the Neutral Compositions
Figure 6a shows the vertical profiles of ozone, HOx (OH + HO 2 ), and NOx (NO + NO 2 ), before the SEP event, 1 day after the onset of the SEP event, 1 day after the end of the SEP event, and 10 days after the SEP event, calculated by the photochemical model using all the 503 chemical reactions.An enhancement of the HOx density is observed in the altitude range 20-60 km with a factor 10 maximum enhancement occurring at 40 km.Depletion of ozone density was confirmed in the altitude range 20-60 km with a factor 10 maximum decrease occurring at 40 km.The altitude range at which a HOx enhancement and an ozone depletion are observed corresponds to the maximum penetration depth of 4.6-46 MeV protons in the Martian atmosphere as shown in Figure 3.An enhancement of NOx density is observed in the altitude range 20-100 km with a factor 100 maximum change occurring at 50 km.During the recovery phase, the effect of the SEP on the HOx and ozone densities is still obvious 1 day after the end of the SEP event, and almost vanishes 10 days after the SEP event.The effect of the SEP on the NOx density still remains even 10 days after the SEP event due to the long chemical time scales of N and NO (∼10 4 -10 5 s at high altitudes).NO molecules move to lower altitudes, maintaining or increasing the NOx density over time, as seen in Figure 7 below 40 km altitude.
In order to clarify whether HOx or NOx contribute to the depletion of the ozone density, we perform another simulation using reactions excluding all the nitrogen-related ones.Figure 6b shows the vertical profiles of ozone and HOx simulated using 238 reactions without all the nitrogen-related ones at the four timings as in Figure 6a.The contribution of NOx to the variation of ozone is negligibly small.The depletion of ozone can be attributed to the enhancement of HOx during the SEP event, not to NOx.
The temporal variation of the O, O 3 , H, OH, and HO 2 densities at 40 km altitude where the variation of ozone is the largest are shown in Figure 7.The density of these species converged in 5 hr due to short chemical timescales, indicating that the duration of the SEP event does not affect the amount of change in density of these species.The temporal variation of NO at 50 km altitude is also shown in Figure 7.The NO density linearly increased during Profiles calculated using all the 503 reactions and (b) profiles calculated using 238 reactions without all the nitrogen-related reactions.Different lines represent the profiles at different timings for each species, before the solar energetic particle (SEP) event (denoted as "before the SEP"), 1 day after the onset of the SEP event (denoted as "during the SEP"), 1 day after the end of the SEP event (denoted as "1 day after the end of the SEP"), and 10 days after the end of the SEP event (denoted as "10 days after the end of the SEP").

Chemical Pathways of Ozone Depletion
First, we describe the chemical pathways that occur in the Martian atmosphere during the SEP event, focusing on the reason for the depletion of ozone.The production of HOx during the SEP event starts from the ionization of CO 2 due to proton impact (p*) and secondary electron impact (e*).
The produced CO 2 + then react with the ambient CO 2 and O 2 to become O 2 + .Reactions (R3) and (R4) are essentially dominant at low altitudes below 80 km because three-body reactions become significant at low altitudes, while those reactions are not efficient in the ionosphere altitudes.
Water cluster ions H 3 O + (OH) and H + (H 2 O) n are finally broken by the recombination with electrons or negative ions (R -) to produce H and OH in the atmosphere.
Depletion of ozone cannot easily be attributed to the reaction with the enhanced H and OH.If either H or OH is the direct reactant to destroy ozone, temporal variation of ozone should coincide with those of H and/or OH considering that the loss time scale of ozone is around 100 s.However, as seen in Figure 7, the temporal variation of ozone does not coincide with those of H and OH; ozone does not follow the peak shape seen in the H and OH variation at 3 hr after the onset of SEP event.On the other hand, temporal variations of O and O 3 are similar to that of HO 2 .According to our simulation, the contribution of the reaction HO 2 + O 3 is essentially negligible for the ozone loss and HO 2 + O is the dominant reaction for the loss of O during the SEP event.Therefore, the following scenario is a plausible explanation for the ozone depletion during the SEP event.The increase of the H and OH densities induce the production of HO 2 via the following catalytic cycle of CO recombination.
During this catalytic cycle, O is lost by the reaction with HO 2 (R17), which results in the decrease of the production rate of ozone by the following main path of ozone production.
The H and OH densities are regulated by the reaction with HO 2 , which results in the decrease of H and OH 3 hr after the onset of the SEP event.On the other hand, the HO 2 density does not decrease because HO 2 is mainly removed by the reaction of O + HO 2 and this reaction rate decreases due to the decrease of the O density.This scenario suggests that the depletion of ozone during the SEP event in the Martian atmosphere is different from that in the terrestrial atmosphere, starting by the ionization of CO 2 , the depletion of ozone is induced by the decrease of the O density due to an enhanced HO 2 density produced by the catalytic cycle of the CO recombination (R16-R18).
Duration of the SEP event could affect the variation of the density profiles.As seen in Figure 7, the NO density increases in time due to the long chemical timescales of N and NO (∼10 4 -10 5 s in the altitude range 60-80 km) and the downward motion of N and NO, so that the enhancement of the NO density is sensitive to the duration of the SEP event.Since the enhancement of HOx density and the decrease of ozone density reached their converged values in 5 hr, the duration of the SEP event does not affect the amount of variation for those species during the SEP event except if the SEP event lasts less than 5 hr.As seen in Figures 8a and 8b, since the H density increases in time above 45 km altitude, the duration of the SEP event could affect the amount of the decrease of the ozone density because it could affect the downward flux of H after the SEP event.
After the end of the SEP event, the recovery of the ozone density occurs through two consecutive phases.In the first phase, the ozone density rapidly increases during 5 hr as displayed in Figure 7.This recovery is purely due to the short lifetime of H, OH, and HO 2 , and the short production timescales of O and O 3 .Indeed the lifetime of a chemical species can roughly be estimated by n i /L i , those of H, OH, and HO 2 are about 100, 4, and 50 s at 40 km altitude, respectively.The timescale of production of a chemical species can be estimated by n i /P i , and those of O and O 3 are 1,000 and 500 s at 40 km altitude, respectively.In the second phase, ozone density slowly recovers to a pre-SEP level on a week timescale.Such slow recovery phase is due to the downward motion of H atoms from the upper altitudes at which the lifetime of H atom is longer.For instance, the lifetimes of H atom at 60, 70, and 80 km altitudes are 1 hr, 14 hr, and 20 days.Lifetime of H atom rapidly increases with altitude because H atom is mainly lost by a three-body reaction with O 2 (R16) and a two-body reaction with O 3 , both of which are less efficient at high altitudes where atmospheric density and ozone density rapidly decrease.Figures 8c and 8d displays the temporal variation of the H density and downward flux after the SEP event.The H density rapidly decreases below 45 km altitude due to its short lifetime, while it slowly decreases above 45 km altitude because of the long lifetime and descent of H from the altitudes above.The downward flux of H increases below 70 km after the end of the SEP event, which supplies enough H to the lower altitudes to maintain its density.The downward velocity of H is around 1 cm s −1 (obtained by dividing the flux by the number density), which means that it takes several days for H atoms to move downward from 80 to 40 km altitudes.

Detectability of Changes in Chemical Composition by TGO/NOMAD
Vertical volume mixing ratio profiles for several species before and at the end of the SEP event are summarized in Figure 9.The detection lower limits in the solar occultation (SO) mode of TGO/NOMAD are 1 ppbv for HO 2 , 0.1 ppbv for NO 2 , 0.001 ppbv for N 2 O, and 0.03 ppbv for H 2 CO, and the detection limit of TGO/ NOMAD ultraviolet and visible spectrometer (UVIS) channel in SO is 0.05 ppbv for O 3 (Vandaele et al., 2015).
However, these estimations were performed for clear sky condition (i.e., no aerosols in the atmosphere), which is not representative for most of the cases.The actual detection limits are estimated to be 10 times worse than those values with the presence of moderate abundance of aerosols (Vandaele et al., 2018).In fact, the detection limit of HCl for clear sky condition was estimated to be 0.03 ppb (Vandaele et al., 2018), however, it is ∼0.3 ppb for the actual observations (Aoki et al., 2021).As shown in Figure 9, the simulated volume mixing ratios are 10-100 ppbv for O 3 , 0.01-1 ppbv for HO 2 , 10 −5 -10 −1 ppbv for NO 2 , 10 −4 -10 −3 ppbv for N 2 O, and 10 −15 ppbv for H 2 CO within the altitude range where changes induced by the SEP event are significant.The simulated change of the concentration of O 3 is sufficiently above the detection limit of TGO/NOMAD even considering the effects of aerosols.Currently more than 10 vertical profiles can be obtained in a day by the instrument, and the results for the observation of the ozone density by TGO/NOMAD have been published by several papers (Khayat et al., 2021;Patel et al., 2021).Therefore, the depletion of ozone during a Halloween-class SEP event should be detected by TGO/ NOMAD.For other species, the simulated changes in the concentration of HO 2 and NO 2 are just at the detection limit of TGO/NOMAD for clear sky condition, so it would be challenging to detect the enhancement of those  species during the SEP event due to the small signal-to-noise ratio.As for N 2 O and H 2 CO, which are important for prebiotic synthesis of amino acid and nucleic acid, the simulated concentrations are far below the detection limit of TGO/NOMAD.
The dayglow and nightglow emission intensities of NO γ and δ bands are expected to change during SEP events due to the enhanced NO density.
Figure 10 shows the variation in the production rate of NO(C 2 Π) due to the recombination of N( 4 S) and O( 3 P) before, during, and after the SEP event.We found that the production rate of NO(C 2 Π) at the peak altitude decreases during the SEP event, indicating the suppression of the NO nightglow emissions (γ and δ bands).That is because the enhanced NO density around 80 km altitudes suppresses the concentration of N( 4 S) by the principal loss reaction for odd nitrogen N( 4 S) + NO → N 2 + O.The production rate of NO(C 2 Π) increases below 60 km altitudes where the N( 4 S) density increases during the SEP event due to the downward transport of N( 4 S) from higher altitudes, but their contribution to the total emission is negligible.The nightglow of NO, produced primarily by the recombination of N( 4 S) and O( 3 P), decreases during SEP events, however, the dayglow of NO is mainly produced by the solar fluorescence of NO, which is directly associated with the density of NO, are expected to enhance during SEP events.Note that our simulations are performed at the subsolar point, and do not represent the variation on the nightside.
There are several limitations to the comparison between our one-dimensional model and future observations.Our simulations were performed at the subsolar point, while the SO observations are performed near the terminator.Since the ozone density rapidly increases across the terminator of Mars from the dayside to the nightside (Piccialli et al., 2021(Piccialli et al., , 2023)), three-dimensional simulation is needed for the direct comparison with SO observations.Furthermore, due to the orbital geometry of TGO spacecraft, the SO observations are performed frequently at high latitudes, where dynamics play an important role in the vertical profile of ozone (Lefèvre et al., 2021).The lack of dynamics in our one-dimensional photochemical model is a major limit in particular at high latitudes.Previous studies have shown that the ozone concentration is known to be strongly anti-correlated with the concentration of water vapor (Lefèvre et al., 2021).Our simulation results should therefore depend on the initial vertical profile of water vapor.Since the goal of our present work is to provide a first prediction of changes in chemical composition during a SEP event in the Martian atmosphere, the effects of dynamics, seasonal variation of water vapor content, and other uncertainties such as seasonal and spatial variation of eddy diffusion coefficient (N.Yoshida et al., 2022) are out of scope.
We would like to note that our simulation suggests that a SEP event also induces enhancement of H 2 O 2 volume mixing ratio at 20-60 km (e.g., from 0.1 to 10 ppbv at 40 km altitude, see Figure 9).This is not detectable by current available instruments, however, H 2 O 2 has strong absorption lines in the sub-mm range and future spectroscopic instruments in that spectral range with limb-viewing (e.g., Kasai et al., 2012) may capture the simulated enhancement of H 2 O 2 .

Dependence on the Intensity and Spectral Shape of the SEP Proton Flux
We discuss the role of the SEP proton flux intensity and of its spectral shape on our results.The SEP event occurring on the 28 October 2003 was one of the most intense SEP events at the Earth, which would occur once every 10 years (Birch & Bromage, 2022;Kataoka, 2020).The only instrument that can measure the absolute intensity of high-energy protons at Mars is the SEP instrument on board MAVEN, with an energy range limited to 6 MeV, which limits the accuracy of fitting parameters for the SEP proton flux.We will therefore use the measured proton flux at 1 MeV as a proxy instead of using normalization constant C in the expression of the proton flux spectrum in Equations 1 and 2 to expect future SEP events which could cause ozone depletion.Equations 1 and 2 were therefore re-written by replacing C by the proton flux at 1 MeV f 1MeV in unit cm −2 s −1 str −1 keV −1 as:   2022) by a factor of 6 to match the observed auroral emission intensity (Nakamura, Terada, Leblanc, et al., 2022) and convert the unit into cm −2 s −1 str −1 keV −1 .During the SEP event in December 2014 on Mars, the f 1MeV value reached 0.8 cm −2 s −1 str −1 keV −1 at peak timing.The range of f 1MeV value of 0.1-10 covers from calm to extreme SEP events.
The results with a slope at low energy γ a = 0.65 are shown in Figure 10, those with a slope at low energy γ a = 1.23 are shown in Figure 11, and those with a slope at low energy γ a = 1.81 are shown in Figure 12.As a whole, a hard-spectral slope at low energy, that is, small value of γ a , a large flux intensity at 1 MeV f 1MeV , and a large break energy E B resulted in the large amount of ozone variation.As already described, ozone depletion occurs in the altitude range 20-60 km corresponding to the penetration of protons with energy 4.6-46 MeV, small value of γ a with large E B lead to the largest intensity of proton flux at high energy, leading to a large production of H and OH in the atmosphere.The dependence of the power spectral slope at high energy γ b varies with E B .The amount of ozone variation is almost independent of γ b for large E B , since the proton flux intensity varies with energy gradually above several times E B as γ b changes, and the proton flux remained the same in the energy range 4.6-46 MeV even as γ b changed.The qualitative signature of γ b dependence did not vary with γ a .The ozone variation does not depend on γ b and is linearly dependent on E B if the break energy E B is larger than 6.18 MeV, near the upper limit of MAVEN/SEP instrument.If the change of slope in the proton flux spectrum does not appear in the MAVEN/ SEP flux data during a future SEP event, we can only obtain γ a and f 1MeV .In this case, we do not have to take care of the power law spectral slope γ b above the upper energy limit of the instrument.The uncertainty would only be related to the break energy E B , leading to an uncertainty in the ozone variation by several factors.If the break energy E B can also be determined using the MAVEN/SEP flux data, we then take care of the uncertainty in γ b , leading to the uncertainty in the ozone variation by several factors.As an example, during the September 2017 SEP event on Mars, γ a , f 1MeV , and E B were 0.8, 5 cm −2 s −1 str −1 keV −1 , and 3 MeV, respectively (Nakamura, Terada, Leblanc, et al., 2022).Those values were similar to the case of Figures 11a and 11b, allowing us to estimate the depletion of ozone density to ∼50-75% at 40 km altitude and 5%-10% in column, which could sufficiently be detected.
Finally, we discuss the frequency of SEP events that cause ozone depletion.Statistical analysis by Gopalswamy (2018) showed the relationship between the proton flux >10 MeV and the cumulative distribution of the number of SEP events for 263 SEP events in 1976-2016.Statistical analysis by Kataoka (2020) for 261 SEP events in 1976-2020 obtained almost the same relationship as Gopalswamy (2018).The proton flux >10 MeV was calculated by integrating the proton flux over energy >10 MeV for each pair of γ a , γ b , f 1MeV , and E B in Figures 11-13, which was then converted into the frequency of SEP events using the relationship by Gopalswamy (2018).The frequency of SEP events for each pair of those spectral parameters is shown in Figures 11-13 in red lines.In all cases, a 50% depletion of the ozone density at 40 km and a 4% depletion of the column ozone density would correspond to the frequency of SEP events 2 year −1 , a 75% depletion of the ozone density at 40 km and an 8%-10% depletion of the column ozone density would correspond to the frequency of SEP events 1 year −1 , a >90% depletion of the ozone density at 40 km and a 14%-16% depletion of the column ozone density would correspond to the frequency of SEP events 0.3 year −1 .Since the depletion of the ozone density occurs around 40 km altitude where SEP protons with energy greater than ∼10 MeV induce reactions (Figure 3a), the depletion of the ozone density is the same for the same proton flux >10 MeV for each pair of spectral parameters.Therefore, a significant depletion of Mars' atmospheric ozone can occur not only during extreme SEP events such as a Halloween-class SEP event, but also even during relatively frequent SEP events at Mars.

Conclusions
We provide the first prediction on the evolution of the atmospheric neutral chemical composition in the Martian atmosphere during a large SEP event by using a Monte Carlo model PTRIP and a one-dimensional photochemical The depletion of the ozone density is expected to be detectable by TGO/NOMAD, while the variation of other species is not.We perform a sensitivity test of ozone variation with respect to the intensity and spectral shape of the SEP proton flux spectrum.We find that a hard spectral slope at low energy and larger break energy results in a large amount of ozone depletion.A 75% depletion of the ozone density at 40 km altitude and an 8%-10% depletion of the column ozone density can be expected during SEP events occurring 1 year −1 on average based on a statistical analysis of SEP events in 1976-2016.Our model reveals, for the first time, that ozone concentration can decrease significantly during a large SEP event in the Martian atmosphere as on Earth, but via different chemical pathways driven by CO 2 ionization and CO recombination catalytic cycle.It is now a good opportunity for us to predict the effects

Table A1
Continued

•
The solar energetic particle (SEP) effects on the neutral chemical composition in the Martian atmosphere are investigated by numerical simulations • Ozone depletion due to enhanced HOx during a large SEP event could be detected by Trace Gas Orbiter/Nadir and Occultation for MArs Discovery • Chemical pathways leading to ozone depletion on Mars differ from those on Earth

Figure 1 .
Figure 1.Vertical profiles of (a) number densities of CO 2 , O 2 , CO, O, N 2 , and H 2 O, and (b) neutral temperature used in the simulation.The profiles of CO 2 , N 2 , and H 2 O were used as the initial ones in the photochemical model PROTEUS, and those of O 2 , CO, and O were the steady state profiles obtained by PROTEUS.The profiles except for H 2 O were used in the Monte Carlo model PTRIP.

Figure 2 .
Figure 2. Energy spectrum of the solar energetic particle (SEP) proton flux during the 28 October 2003 SEP event scaled to the Mars' position at 1.5 AU.

Figure 3 .
Figure 3. (a) Ionization rate of CO 2 calculated by PTRIP at each incident energy using the solar energetic particle (SEP) proton flux of a Halloween-class SEP event.(b) Vertical profiles of the production rates of CO 2 + , CO + , O + , C + , O 2 + , N 2 + , N + , N( 4 S), and N( 2 D) calculated by PTRIP.

Figure 4 .
Figure 4. Vertical profiles of the ion species calculated by PROTEUS.(a) Before the solar energetic particle (SEP) event and (b) 1 day after the onset of the SEP event.

Figure 5 .
Figure5.(a and b) Temporal variations of the electron density profile during the solar energetic particle (SEP) event and after the end of the SEP event, respectively.Horizontal axes are the time from the onset of the SEP event and the time from the end of the SEP event, respectively.

Figure 6 .
Figure 6.Vertical profiles of O 3 , HOx (OH + HO 2 ), and NOx (NO + NO 2 ) simulated by the photochemical model.(a)Profiles calculated using all the 503 reactions and (b) profiles calculated using 238 reactions without all the nitrogen-related reactions.Different lines represent the profiles at different timings for each species, before the solar energetic particle (SEP) event (denoted as "before the SEP"), 1 day after the onset of the SEP event (denoted as "during the SEP"), 1 day after the end of the SEP event (denoted as "1 day after the end of the SEP"), and 10 days after the end of the SEP event (denoted as "10 days after the end of the SEP").

Figure 7 .
Figure 7. (a) Temporal variation of the O, O 3 , H, OH, and HO 2 densities at 40 km altitude and the NO density at 50 km from the onset of the solar energetic particle (SEP) event to 1 day after the end of the SEP event.The horizontal axis is the time from the onset of the SEP event on the left panel and the time from the termination of the SEP event on the right panel.The timings of the onset and the end of the SEP event are shown above each panel.(b) Temporal variation of the O, O 3 , H, OH, and HO 2 densities at 40 km altitude and the NO density at 50 km from the end of the SEP event to 10 days after.

Figure 9 .
Figure9.Vertical volume mixing ratio profiles for several species before the solar energetic particle (SEP) event (dashed curves) and after the end of the SEP event (solid curves).The mixing ratio of H 2 CO is multiplied by a factor of 10 12 .

Figure 8 .
Figure 8. (a and b) Simulated vertical profiles of hydrogen atom (H) density and downward flux during the solar energetic particle (SEP) event, respectively.(c and d) Vertical profiles of H density and downward flux after the end of the SEP event, respectively.

Figure 10 .
Figure10.The reaction rate of N( 4 S) + O( 3 P) → NO(C 2 Π) before (black), 1 day after the onset of solar energetic particle (SEP) event (red), and 10 days after the end of SEP event.
the Band parameters for the SEP proton spectra was given byDesai et al. (2016).They analyzed 46 SEP events for about 16 years from 1998 to 2014 to characterize the variation of the Band parameters among SEP events.They found that the mean values for the low-and high-energy power law slopes γ a and γ b were 1.23 and 3.63, respectively.The standard deviation (σ) of γ a and γ b were 0.58 and 1.12, respectively.We re-analyzed the distribution of E B for 44 of 46 events inDesai et al. (2016) at which E B were determined and find that the 25th, 50th, and 75th percentile values of the break energy E B were 2.92 MeV, 6.18 MeV, and 12.4 MeV, respectively.We investigated the ozone variation as a function of γ b and f 1MeV for a value of γ a = 0.65, 1.23, and 1.81 (ranging from mean value minus 1σ to mean value plus 1σ), and E B of 2.92 MeV, 6.18 MeV, and 12.4 MeV.The f 1MeV value range was set to 0.1-10 cm −2 s −1 str −1 keV −1 according to previous SEP events.The f 1MeV value was equal to 5 cm −2 s −1 str −1 keV −1 during the SEP event on 28 October 2003.During the SEP event in September 2017 on Mars, which was the most intense SEP event observed since MAVEN spacecraft insertion around Mars, the f 1MeV value was equal to 5 cm −2 s −1 str −1 keV −1 if we multiply the SEP proton flux in Figure4ofNakamura,  Terada, Leblanc, et al. (

Figure 11 .
Figure 11.Dependence of the ozone density variation on parameters of solar energetic particle (SEP) proton flux spectra.Left panels represent the variation of the ozone density at 40 km altitude and right panels represent the variation of the ozone column density.Vertical axis is the power law spectral slope at high energy γ b and horizontal axis is the proton flux at 1 MeV.The power law spectral slope at low energy γ a is 0.65 in all panels.(a and b) The break energy E B = 2.92 MeV, (c and d) E B = 6.18MeV, and (e and f) E B = 12.4 MeV.The red solid, dashed, dash-dot, and doted lines represent the frequency of SEP event, once every 3 months, 6 months, 1 year, and 3 years, respectively.
(H 2 O) 2 via the following chain reactions.OH radical is produced through the production of H + (H 2 O), and H + (H 2 O) 2 .
R4) O 2 + then reacts with the ambient CO 2 and H 2 O molecules to produce water cluster ions H 3 O + (OH), H + (H 2 O), and H +