Modeling of very low frequency ( VLF ) radio wave signal profile due to solar flares using the GEANT 4 Monte Carlo simulation coupled with ionospheric chemistry

X-ray photons emitted during solar flares cause ionization in the lower ionosphere ( ∼60 to 100 km) in excess of what is expected to occur due to a quiet sun. Very low frequency (VLF) radio wave signals reflected from the D-region of the ionosphere are affected by this excess ionization. In this paper, we reproduce the deviation in VLF signal strength during solar flares by numerical modeling. We use GEANT4 Monte Carlo simulation code to compute the rate of ionization due to a M-class flare and a X-class flare. The output of the simulation is then used in a simplified ionospheric chemistry model to calculate the time variation of electron density at different altitudes in the D-region of the ionosphere. The resulting electron density variation profile is then self-consistently used in the LWPC code to obtain the time variation of the change in VLF signal. We did the modeling of the VLF signal along the NWC (Australia) to IERC/ICSP (India) propagation path and compared the results with observations. The agreement is found to be very satisfactory.


Introduction
The most prominent causes of daytime lower ionospheric modulation (other than the diurnal variation) are the solar flares.The lowest part of the ionosphere, namely, the Dregion, formed in the daytime due mainly to the solar UV radiation, is significantly affected by these flares.The Xray and γ -ray from the solar flares penetrate down to the lower ionosphere causing enhancements of the electron and the ion densities (Mitra, 1974), though the part of the photon spectrum ranging from ∼2 to 12 keV is responsible for the modulation at the reflection heights of very low frequency (VLF) radio waves.The VLF signal (3-30 kHz) reflected mainly from the D-region ionosphere carries information of such enhancements.By analysing these signals between a transmitter and a receiver one could, for example, obtain enhancements in the electron-ion densities during such events (McRae et al., 2004).
In general, the entire radiation range of EUV, solar Lymanα, β, X-ray and γ -rays of solar origin (Madronich and Flocke, 1999) and from compact and transient cosmic objects as well as protons from cosmic rays and solar events act as perturbative sources and lead to the evolution of the entire atmospheric chemistry.The amount of energy deposited and the altitudes of deposition depend on the incoming radiation energy.Thus the chemical changes in different layers vary.Many workers have studied these changes theoretically as well as using laboratory experiments.In F2 region, the chemistry is simpler with N + , H + and O + with O + as the dominant ions (Mitra, 1974).Abundance of O + depends on solar UV rays (Wayne, 2000).Similarly O + and N + are dominant in F1 region and their densities are governed by EUV of 10-100 nm range (Mitra, 1974).E-region is mostly filled up by O 2 , N 2 and NO molecules.During higher energetic X-ray emission (1-10 nm), O 2 + and N 2 + are produced.Here, the molecular NO and O are formed mainly through dissociative recombination process.Mariska and Oran (1981) measured the evolution of E and F regions during solar flares using EUV data from the SOLRAD satellite system.The response of the D-region is dominated by hard X-ray and γ Ray injection from flares and Gamma Ray Bursts (Mitra, 1974;Wayne, 2000;Thomson et al., 2001;Höller et al., 2009;Tripathi et al., 2011;Mondal et al., 2012).It has nearly similar ion constituents as in the E-region along with hydrated ions like O 2 + (H 2 O) n , H + (H 2 O) n , etc. Several authors such as Mitra (1974), Mcrae (2004), Pal et al. (2012b), studied D-region electron density in detail.Beside these radiation events, regular solar proton events (SPE) having energies 1-500 MeV affect the atmospheric chemistry drastically as one such proton can ionize millions of molecules.Joeckel et al. (2003) reported interesting observations on 14 CO during SPE.Holloway and Wayne (2010), Jackman et al. (2001), Verronen et al. (2005), etc. have reported the effects on the concentrations of 14 CO, NO, NO 2 and O 3 in mesospheric and stratospheric regions.
Because the earth's ionosphere (and atmosphere in general) is a gigantic, noisy, yet free of cost detector of cosmic radiation, study of the evolution of its chemical composition in presence of high energy phenomena (both photons and charged particles) is of utmost importance not only to interpret the results of VLF and other radio waves through it, but also to understand long-term stability of this sensitive system.In the present paper, we apply the GEANT4 simulation (applicable to any high energy radiation detector) only to understand the changes in the D-region electron density (N e ) during flares as the ion-chemistry leaves negligible impression on VLF propagation characteristics.For that reason knowledge of only the dominant chemical processes in the D-region ionosphere would be enough to study the modulation due to events like solar flares (Rowe et al., 1974;Mitra, 1981;Chamberlain, 1987).One such simplified model, namely, the Glukhov-Pasko-Inan (GPI) model (Glukhov et al., 1992) can be applied for the events which cause excess ionization in the D-region.Inan et al. (2006) used this model (though modified for lower altitudes) to find the effects of Gamma Ray Bursts on the state of ionization.The model has also been adopted successfully by Haldoupis (Haldoupis et al., 2009) for his work on the early VLF perturbations associated with transient luminous events.
In the present paper, we model the VLF signal variation due to the D-region ionospheric modulation during solar flares by combining a Monte Carlo simulation (GEANT4) for electron-ion production with the GPI model in the Dregion to find the rate of free electron production at different altitudes by X-rays and γ -rays emitted from flares.Finally, the Long-Wave Propagation Capability (LWPC) code (Ferguson, 1998) has been used to simulate changes in VLF amplitude due to resulting variations of the electron number density.The LWPC code computes the amplitude and phase of the VLF signal for any arbitrary propagation path and ionospheric conditions along the path, including the effects of the earth's magnetic field.It uses the exponential D-region ionosphere defined by Wait and Spies (1964) and wave guide mode formulations developed by Budden (1951).
This program has been used by many workers and is quiet successful for modeling of long-range propagation of VLF signals even in presence of ionospheric anomalies (Grubor et al., 2008;Rodger et al., 1999;Chakrabarti et al., 2012;Pal et al., 2012a, b).Instead of Wait's model for the D-region ionosphere, we used our own result for that part of the ionosphere, by incorporating the results from the GEANT4 Monte Carlo simulation and the GPI scheme in the LWPC program.
We compared the simulated VLF signal amplitude deviations due to the effect of solar flares with those of the observational data detected by the ground based VLF receiver of the Ionospheric and Earthquake Research Centre (IERC) of Indian Centre for Space Physics (ICSP).To illustrate the success of our method, we use one X2.2 type and one M3.5 type solar flare data which occurred on 15 and 24 February 2011, respectively, and compare with the predicted VLF output from our model.Our interest is to reproduce the change in the VLF signal during the solar event and examine the decay pattern of the signal due to the long recovery time.
The organization of the paper is as follows: in Sect.2, we present the observational data.In Sect.3, we present detailed procedures to model the data.This includes successive usage of the GEANT4 simulation, GPI scheme and the LWPC code.In Sect.4, we present the results of our modeling which we compared with our observed data.Finally, in Sect. 5 we make concluding remarks.

Observational data
The VLF receivers at the Indian Centre for Space Physics (ICSP) are continuously monitoring several worldwide VLF transmitter signals including NWC (19.8 kHz), VTX (18.2 kHz) and JJI (22.2 kHz).We use the VLF data for the NWC to IERC/ICSP propagation path (distance 5691 km) to compare with our simulation results.In Fig. 1a we present an example of variations of the NWC signal at 19.8 kHz for a typical solar active day (solid curve), with many flares and a normal "quiet" day (dashed curve).There is a sudden rise of wave amplitude followed by a slow decay associated with each of the flares.SoftPAL instrument was used in obtaining the data.In the present paper, we analyze stronger flares (one X2.2 type and one M3.5 type) to begin with.In future, we shall target weaker flares as they require more careful subtraction of the background.While analysing flares in question, we use the GOES and RHESSI satellite data (Sui, et al., 2002) to obtain the light curve and spectrum of those flares.

GEANT4 Monte Carlo simulation
GEANT4 is a well-known object-oriented detector simulation program (Agostinelli et al., 2003)

GEANT4 Monte Carlo simulation
GEANT4 is a well known object-oriented detector simulation program (Agostinelli et al., 2003), used widely for the study of the detector and instrument response to high energy radiations.The earth's ionosphere is a giant detector and as such, the software used to analyse detectors in high energy physics could be used here as well.Most of the ionizations occur in the ionosphere due to the collision of molecules with secondary electrons produced by initial photo-ionization.The program includes sufficiently adequate physical processes for the production of electron ion pairs in the atmosphere by energetic photon interactions.GEANT4 (GEANT is the abbreviation of Geometry and Tracking) is a openly available and 4 energy radiations.The earth's ionosphere is a giant detector and as such, the software used to analyze detectors in high energy physics could be used here as well.Most of the ionizations occur in the ionosphere due to the collision of molecules with secondary electrons produced by initial photo-ionization.The program includes sufficiently adequate physical processes for the production of electron ion pairs in the atmosphere by energetic photon interactions.
GEANT4 (GEANT is the abbreviation of Geometry and Tracking) is a openly available and widely tested toolkit (freely available from the website http://geant4.cern.ch/).The Monte Carlo processes in this application deal with the transport and interaction of particles with matter in a predefined geometry and store all the information of track, nature and number of secondaries efficiently.With this application we can find accurately the ionization processes and tracks without any prior knowledge and assumption of chemical or physical model of the region.The included electromagnetic processes related with X-ray incidence are well developed and suitable for application in the context of the ionosphere.For simulation with UV, the application at this stage is not quite suitable, but with some development in the toolkit area, it is possible to use it for Monte Carlo simulation with UV and EUV photons in the ionosphere.We simulated only the rate of ionization at different altitudes due to solar X-ray photons during flares using GEANT4.For the chemistry part related to the interaction of the produced electron, ions and the neutrals in the D-region ionosphere, we adopt the GPI scheme.
The earth's atmosphere, constructed in the detector construction class consists of concentric spherical shells, so that the atmosphere is divided into several layers.The distribution of layers is as follows, 100 layers above the earth surface, 1 km each followed by 20 layer above this, 5 km each and 12 more layers, 25 km each.They are formed with average molecular densities and other parameters at those heights.
The data is obtained from NASA-MSISE-90 atmospheric model (Hedin, 1991) of the atmosphere.
Initially, the Monte Carlo process is triggered by primary particle generation class of GEANT4 with photons of energy ranging from 1-100 keV.We are not interested in higher energy gamma ray photons as their energy deposition height is much lower than that of the D-region ionosphere.First a uniform spectrum (equal number of photons per keV bin) is chosen for input in primary generation class, such that the incident photons in a specific bin have random energies within the range of the bin.This gives a distribution of electron-ion production rate in a matrix form, the ij th element of which gives the electron production rate (cm −3 s −1 ) in i th km above sea level due to the incident photons in the j th keV bin.We produce such altitude dependent distribution by simulation with different incident angles of the photons.The real altitude distributions of electron-ion production rate during a flare are obtained by multiplication of the matrix of corresponding time (hence angle of incidence) with the actual flux spectrum (see Fig. 2) at that time.For this we collect the back-ground subtracted solar spectra during the flares from RHESSI satellite data (http://hesperia.gsfc.nasa.gov/rhessi2/).This method helps us to save a significant amount of computation time.
In some previous work (e.g., Glukhov et al., 1992;Haldoupis et al., 2009) the rate of ionization was approximated from the total energy deposition by a photon, by dividing it with a pre-assigned value of the average energy for production of an electron-ion pair.We did not require such an approximation, instead, we follow the electron ionization interactions down to the lowest possible level in the Monte Carlo process.For this, the lower energy threshold for electron and photon processes is extended down to 10 eV in the electromagnetic processes in GEANT4 toolkit.From simulations, we found the average energy required per electron-ion pair production to be ∼ 31 eV.

S. Palit et al.: Modeling VLF signals due to Solar flares with GEANT4 simulation
processes is extended down to 10 eV in the electromagnetic processes in GEANT4 toolkit.From simulations, we found the average energy required per electron-ion pair production to be ∼ 31 eV.
In Fig. 3 we show the altitude distribution of the rate of electron density produced by the spectra shown in Fig. 2(a-b) at the peak flare time assuming the Sun to be at the zenith.From Fig. 3(a-b), it is evident that the ionization effect due to the solar X-rays is dominant in the D-region (∼ 60-100 km) only.The simulation also shows that the altitude of maximum ionization is somewhat lower for X-class flare (∼ 81 km) than that of (∼ 88 km) M-class.This is consistent with 6 In Fig. 3 we show the altitude distribution of the rate of electron density produced by the spectra shown in Fig. 2(a-b) at the peak flare time assuming the Sun to be at the zenith.From Fig. 3(a-b), it is evident that the ionization effect due to the solar X-rays is dominant in the D-region (∼ 60-100 km) only.The simulation also shows that the altitude of maximum ionization is somewhat lower for X-class flare (∼ 81 km) than that of (∼ 88 km) M-class.This is consistent with 6 In Fig. 3 we show the altitude distribution of the rate of electron density produced by the spectra shown in Fig. 2a and b at the peak flare time assuming the sun to be at the zenith.
From Fig. 3a and b, it is evident that the ionization effect due to the solar X-rays is dominant in the D-region (∼60-100 km) only.The simulation also shows that the altitude of maximum ionization is somewhat lower for X-class flare (∼81 km) than that of (∼88 km) M-class.This is consistent with the fact that the relative abundance of high energy photons is larger in X-class spectra.
The effects of the zenithal variation of the sun is included in the simulation to account for the variation of the flux in long duration flares.Fig. 4 shows the simulated electron production rate (cm −3 s −1 ) at various altitudes for the M-class spectrum (Fig. 2a) at the flare peak if it had occurred at different times of the day.In finding the rate of production of electron density over the whole period of the flares, we convolve the normalized (by corresponding spectrum) electron density distribution profiles with the light curves of the flares, taking into account the dependency of the electron density on the angular positions.
We simulate the condition at the middle of the VLF propagation path (Fig. 1b) and assume that the resulting electron density distribution apply for the whole path.Along the propagation path in question, the solar zenith angle variation is at the most ∼ 30 • .Figure 4 shows that the approximation is justifiable, especially if the solar zenithal angle is close to 0 at the mid-propagation path.
Using this method, we compute the free electron production rates per unit volume for both of the flares we wish to model.Figure 5a and b show the GOES light curve of the Mclass and the X-class flares respectively.The corresponding free electron production rates per unit volume, at different altitudes as time proceeds during the flares, obtained from the Monte Carlo simulation are shown in Fig. 5c and d.The latter curves generally follow the variation of X-rays.
ing electron density distribution apply for the whole path.Along the propagation path in the solar zenith angle variation is at the most ∼ 30 • .Figure 4 shows that the approximation le, especially if the solar zenithal angle is close to 0 at the mid-propagation path.
iation of electron production rates with the angular position of the Sun at different heights assuming spectrum.
his method, we compute the free electron production rates per unit volume for both of the wish to model.Figs.5a and 5b show the GOES light curve of the M-class and the X-class ectively.The corresponding free electron production rates per unit volume, at different s time proceeds during the flares, obtained from the Monte Carlo simulation are shown in nd 5d respectively.The latter curves generally follow the variation of X-rays.

GPI model
goal is to compute the changes in atmospheric chemistry due to the flare.These could be using a simplified model, such as the GPI model (Glukhov et al., 1992).In this model, the ensity evolution of four types of ion species are taken into account for the estimation of the ty changes in the D-region.These are: electrons, positive ions, negative ions and positive 7 Fig. 4. Variation of electron production rates with the angular position of the sun at different heights assuming an M-class spectrum.

The GPI model
Our next goal is to compute the changes in atmospheric chemistry due to the flare.These could be obtained using a simplified model, such as the GPI model (Glukhov et al., 1992).In this model, the number density evolution of four types of ion species are taken into account for the estimation of the ion density changes in the D-region.These are: electrons, positive ions, negative ions and positive cluster ions.The positive ions N + comprise of mainly O + 2 and NO + 2 .The negative ions N − include O − 2 , CO − 3 , NO − 2 , NO − 3 etc.The positive cluster ions N + x are usually of the form H + (H 2 O) n .The GPI model is meant for the nighttime lower ionosphere.Modified set of equations exists (Lehtinen and Inan, 2007), which are applicable at even lower heights (∼50 km), where the inclusion of negative ion clusters N − x is important and the scheme has been used for the daytime lower ionosphere also (Inan et al., 2007).However, for the altitudes (above 60 km) of our interest, where the VLF wave is reflected, the simplified GPI model is sufficient.While considering the solar photons, we concentrate only on the ionization produced by the X-rays from solar flares, not from the regular UV and Lyman alpha lines etc. Though it is a simplification, this will produce good results as long as the electron-ion production rates due to X-ray dominates over the regular production rates.The model is suitable for the strong flares (at least higher C-class and above) and is valid up to the time of recovery when the above approximation holds.
According to the model, the time evolution of the ion densities follow ordinary differential equations, the right hand sides of which are the production (+ve) and loss (-ve) terms.
The neutrality of the plasma requires, Here, β is the electron attachment rate (Rowe et al., 1974), the value of which is given by, N O 2 and N N 2 represent the number densities of molecular Oxygen and Nitrogen respectively, and T is the electron temperature.The neutral atom concentrations at different altitudes, required for the calculations of the coefficients, were obtained from the NASA-MSISE atmospheric model.The detachment coefficient, i.e., the coefficient of detachment of electrons from negative ions has contributions from both photodetachment (at daytime) and also processes other than photodetachment.The coefficient for the photodetachment process is ∼ 10 −2 to 5 × 10 −1 and for the other process it is ∼ 10 −1 to 10 −2 at 80 km (Bragin, 1973).According to Lehtinen and Inan (2007) and Pasko and Inan (1994), the value of the detachment coefficient γ varies widely from 10 −23 N s −1 to 10 −16 N s −1 .Following Glukhov et al. (1992) the value of γ is taken to be 3×10 −18 N s −1 , where N is the total number density of neutrals.In our calculations the value of γ comes out to be 2.4 × 10 −2 at 75 km to 6.6 × 10 −3 at 84 km, which is reasonable.Previous studies (Lehtinen and Inan, 2007;Rowe et al., 1974) predict that the effective coefficient of dissociative recombination α d may have values from 10 −7 to 3 × 10 −7 cm −3 s −1 .We have used 3 × 10 −7 cm −3 s −1 for α d .
The effective recombination coefficient of electrons with positive cluster ions α c d has the value ∼ 10 −6 -10 −5 cm −3 s −1 .Here the value of 10 −5 cm −3 s −1 is adopted as suggested by Glukhov et al. (1992).The value of the effective coefficient of ion-ion recombination processes for all types of positive ions with negative ions is taken to be 10 −7 cm −3 s −1 (Mitra, 1968;Rowe et al., 1974).B is the effective rate of conversion from the positive ion (N + ) to the positive cluster ions (N + x ) and has the value 10 −30 N 2 s −1 (Rowe et al., 1974;Mitra, 1968), in agreement with Glukhov et al. (1992).
The term I , which corresponds to the electron-ion production rate, can be written as I = I 0 + I x .Here I 0 is from the constant source of ionization which are responsible for the  (Lehtinen and Inan, 2007), that are applicable at even lower heights (∼ 50 km), where the inclusion of negative ion clusters N − x is important and the scheme has been used for the day time lower ionosphere also (Inan et al., 2007).However, for the altitudes (above 60 km) of our interest, where the VLF wave is reflected, the simplified GPI model is sufficient.While considering the solar photons, we concentrate only on the ionization produced by the X-rays from solar flares, not from the regular UV and Lyman alpha lines etc. Though it is a simplification, this will produce good results as long as the electronion production rates due to X-ray dominates over the regular production rates.The model is suitable for the strong flares (at least higher C-class and above) and is valid up to the time of recovery when the above approximation holds.
According to the model, the time evolution of the ion densities follow ordinary differential equations, the right hand sides of which are the production (+ve) and loss (-ve) terms.from the GPI scheme.We show this for the X-class flares only.It is to be noted that they are the outcome of simplified chemical model devoid of physical processes such as the transport of ions ambient distributions of ions.In our case, I 0 corresponds grossly to the part produced by the background ultraviolet ray, Lyman alpha lines, etc. from the sun, producing the D region.The term I x is the modulated component of the ionization rate.Here I x is generated by the X-ray photons (in the energy range from 2 to 12 keV, important for the D-region) coming from the solar coronal region during solar flares only.
To account for the time variation of the number density of the ion species, we divide each term into a constant and a time dependent parts, i.e., N e is replaced with N 0e + N e (t), N + with N + 0 + N + (t), N − with N − 0 + N − (t) and N + x with N + 0x + N + x (t).A set of equations governing the equilibrium conditions can be obtained from Eqs. (1-4), by replacing the ion density terms with their ambient values, I by I 0 and equating the rates to zero (Glukhov et al., 1992).
Using the conditions of equilibrium and the above replacements in Eqs.(1-4), we get the following set of equations from the GPI scheme.We show this for the X-class flares only.It is to be noted that they are the outcome of simplified chemical model devoid of physical processes such as the transport of ions through the layers of Ionosphere etc.These processes are being studied and would be reported elsewhere.

Modeling VLF response with LWPC
Now that we have the electron densities (N e ) at various altitudes, we can use them to obtain the time variation of the deviation of the VLF signal amplitude during a flare event.We use the Long Wave Propagation Capability (LWPC) code (Ferguson, J.A., 1998) to accomplish this.LWPC is a very versatile and well-known code and the current LWPC version has the scope to include arbitrary ionospheric variations as the input parameters for the upper waveguide.The electron density profiles (Fig. 6) represent the modulated D-region ionosphere due to the solar flares and can be used as 11 The ambient values of the electron and ion densities are obtained from the IRI model (Rawer et al., 1978).Unperturbed negative ion density N − 0 is obtained from the relationship N − 0 = λN 0e , where λ is the relative composition of the negative ions.The values of positive ion densities are calculated from the charge neutrality condition, i.e., Using the ionization rate per unit volume (I x ) obtained from the GEANT4 simulation and the coefficients of the ODEs calculated from MSIS data, we solve the ODEs (Eqs.7-10) by Runge-Kutta method to find the residual electron density at various altitudes during the course of the flare.Figure 6a and b show the electron density variations at different altitudes throughout the flares.The enhancement in electron density from ambient values can be clearly seen.
Figure 7a and b shows the altitude variations of the electron densities as obtained from the GPI model.These are plotted at three different times including the peak time.The times are shown by arrows in Fig. 6a and b, The ambient electron density during the X-class flare was higher than that of the M-class flare.This is because the background X-ray flux during the X-class flare was stronger.
Though we do not expect drastic changes in the VLF radio signal due to the changes in the chemical composition at the altitudes of our interests, for the sake of completeness, we show in Fig. 8, the altitude profile of three types of ion densities at three different times evaluated numerically from the GPI scheme.We show this for the X-class flares only.It is to be noted that they are the outcome of simplified chemical model devoid of physical processes such as the transport of ions through the layers of Ionosphere etc.These processes are being studied and would be reported elsewhere.

Modeling VLF response with LWPC
Now that we have the electron densities (N e ) at various altitudes, we can use them to obtain the time variation of the deviation of the VLF signal amplitude during a flare event.We use the Long-Wave Propagation Capability (LWPC) code (Ferguson, 1998) to accomplish this.LWPC is a very versatile and well-known code and the current LWPC version has the scope to include arbitrary ionospheric variations as the input parameters for the upper waveguide.The electron density profiles (Fig. 6) represent the modulated D-region ionosphere due to the solar flares and can be used as the inputs of the LWPC code to simulate the VLF signal behavior.The electron-neutral collision frequency in the lower ionosphere plays an important role for the propagation of VLF signals.We used the collision frequency profiles between electrons and neutrals as described by Kelley (2009).This is given by the following equation, ν e (h) = 5.4 × 10 −10 n n T 1/2 , where n n and T are the neutral density (in cc) and electron temperature (in K), respectively.Assuming the temperature of the electrons and the ions are same and using the ideal gas equation, the above equation can be expressed as (Schmitter, 2011), ν e (h) = 3.95 × 10 12 T −1/2 exp(−0.145h).

Results
So far, we used various ways to obtain altitude variations of the electron densities before and during flares.We use the LWPC model to compute the resulting VLF amplitude.To find the observational deviation in the VLF signal due to flares, we subtract the flare day data from the VLF data of the nearest quiet day.For the signal generated by the LWPC code the inputs of the LWPC code to simulate the VLF signal behaviour.The electron-neutral collision frequency in the lower ionosphere plays an important role for the propagation of VLF signals.We used the collision frequency profiles between electrons and neutrals as described by Kelley (2009).This is given by the following equation, where, n n and T are the neutral density (in cc) and electron temperature (in K) respectively.Assuming the temperature of the electrons and the ions are same and using the ideal gas equation, the above equation can be expressed as (Schmitter, 2011), ν e (h) = 3.95 × 10 12 T −1/2 exp(−0.145h).
12   also the subtraction is made, but with a flat VLF profile obtained from the ambient electron density shown in Fig. 6a and  b.The variation of ambient ion densities during the course of the day affects the recombination process of the modulated D-region ionosphere when the flare duration is long.This is evident from Eqs. ( 7)-( 10) that the rate of change of the electron density depends on the ambient concentration.If it was not the case, the subtraction would have made no difference.
The response of the VLF signal to this electron distribution as obtained from the LWPC code after appropriate subtraction of the background is shown in Fig. 9a, b.We see that the variation in VLF signal strengths resemble with those recorded for NWC to IERC/ICSP propagation path by the IERC monitors.For the M-class flare, the resulting VLF signals from the model matches with the observation up to a certain time beyond which modelled slope becomes gradually lower than that of the observation.For the X-class flare, the modelled slope becomes gradually higher than that of the observation.The M-class flare occurred at the local afternoon when the diurnal ambient electron densities at various altitudes are decreasing.During the decay of the flare when the effect of ambient ion densities start to dominate, observed slope of decay become higher in value (than that predicted from the model, where this effect is ignored) while still following the trend of variation of observed VLF amplitude at that time.The X-class flare, on the other hand, occurred in the morning when ambient electron densities are increasing with time due to the diurnal variation.In this case, the slope of observed decay is a lower.
It can also be seen (comparing Fig. 9 with Fig. 6) that the differences in slopes of decay between the models and the observations start to occur when the electron densities fall below a certain level for both the flares (∼2000 cm −3 at 80 km).Furthermore, since the M-class flare spectrum is relatively softer than that of the X-class one, the ambient electron density starts to dominate the VLF amplitude a bit earlier during decay of the M-class flare.We can see that the deviation of the modelled slope from the observed one starts earlier for M-class flare than the other.

Concluding remarks
In this paper, we have used a three step approach to compute the deviation of the VLF signal amplitude during a solar flare.The steps are (i) GEANT4 Monte Carlo simulation for ionization rate estimation, (ii) GPI-scheme of ionospheric chemistry analysis and finally (iii) LWPC modeling of subionospheric VLF signal propagation.This combined model is self-consistent, realistic and uses a few justifiable assumptions.Our method is capable of estimating the VLF signal behavior directly using the corresponding incident solar flux variation on D-region ionosphere during a solar flare.
In our model computation, we have ignored the effects of the variation of ambient concentrations, i.e., the diurnal variations (ambient) of N 0e , N + 0 , etc. are not taken into account during the flares while solving the ODEs.So as long as the electron densities produced by the X-rays from the flares dominate over the ambient electron densities due to regular ionization sources, such as the ultraviolet, Lyman alpha photons, etc., the model should follow observed variations in the VLF data.We can see in Fig. 9 that it is in fact the case.In future, we shall study the effects of the occurrence times of the flare on the modeling and also the effects of radiation on chemical changes in other altitudes using our method.
We are interested in further modeling of the ionization effect of EUV and UV photons so that we can remove the discrepancies between the modelled and observed decay of the perturbation and extend the modeling throughout the duration of the flare.This will require a knowledge of some more details of the ionospheric chemistry and fluid dynamics.Also, the inclusion of inhomogeneous distribution of ionospheric dynamics during flares and adoption of more realistic values of ionospheric parameters will lead to results in better agreement with observations.In future, these modifications will be made in the time analysis of the flare decay and the sluggishness of ionosphere.

Fig. 2 .
Fig. 2. Spectra at three phases (dotted: ascent phase; solid: peak; dash-dotted: descent phase of the (a) M-class flare and the (b) X-class flare which are used in this paper.The spectra were obtained from the RHESSI satellite data and the OSPEX software.

Fig. 3 .
Fig. 3. Electron production rates per unit volume at peak times for the (a) M-class flare and the (b) X-class flare assuming the Sun to be at zenith.

Fig. 2 .
Fig. 2. Spectra at three phases (dotted: ascent phase; solid: peak; dash-dotted: descent phase) of the (a) M-class flare and the (b) X-class flare which are used in this paper.The spectra were obtained from the RHESSI satellite data and the OSPEX software.

Fig. 2 .
Fig. 2. Spectra at three phases (dotted: ascent phase; solid: peak; dash-dotted: descent phase of the (a) M-class flare and the (b) X-class flare which are used in this paper.The spectra were obtained from the RHESSI satellite data and the OSPEX software.

Fig. 3 .
Fig. 3. Electron production rates per unit volume at peak times for the (a) M-class flare and the (b) X-class flare assuming the Sun to be at zenith.

Fig. 3 .
Fig. 3. Electron production rates per unit volume at peak times for the (a) M-class flare and the (b) X-class flare assuming the sun to be at zenith.

Fig. 5 .
Fig. 5. Light curves of the flares from GOES satellite of an (a) M-class and (b) X-class flare.The time variation of the electron production rates at different heights as obtained from GEANT4 simulation for these flare respectively are shown in (c) and (d).

8Fig. 5 .Fig. 6 .Fig. 7 .
Fig. 5. Light curves of the flares from GOES satellite of an (a) M-class and (b) X-class flare.The time variation of the electron production rates at different heights as obtained from GEANT4 simulation for these flare respectively are shown in (c) and (d).

Fig. 6 .
Fig. 6.Electron density at different layers during the (a) M-class and (b) X-class flares.

Fig. 7 .
Fig. 7. Altitude profiles of electron density as obtained from GPI model at three different times for an (a) M-class and (b) X-class flare.Solid line shows the ambient profile before the beginning of the flare while the dashed and the dot-dashed lines represent the profile during the peak and at a time during recovery.

Fig. 8 .
Fig. 8. Altitude profiles of three types of ion densities as obtained from the GPI scheme: (a)Negative ions, (b) Positive ions, and (c) Positive Cluster ions.The solid, dashed and dot-dashed lines represent respectively the ambient profile, the profile at the peak and the profile at a time during recovery of the X-type flares as discussed in the text.

Fig. 8 .
Fig. 8. Altitude profiles of three types of ion densities as obtained from the GPI scheme: (a) negative ions, (b) positive ions, and (c) positive cluster ions.The solid, dashed and dot-dashed lines represent respectively the ambient profile, the profile at the peak and the profile at a time during recovery of the X-type flares as discussed in the text.

Fig. 9 .
Fig. 9. Simulation results of VLF amplitude modulations (dotted line) and corresponding observed VLF data (solid line) are plotted with UT for an (a) M-class flare and an (b) X-class flare respectively.

Fig. 9 .
Fig. 9. Simulation results of VLF amplitude modulations (dotted line) and corresponding observed VLF data (solid line) are plotted with UT for an (a) M-class flare and an (b) X-class flare, respectively.