The Effect of Resistivity on the Periodicity of Oscillatory Reconnection

The oscillatory reconnection mechanism is investigated for a parameter study of eight orders of magnitude of resistivity, with a particular interest in the evolution of the oscillating current density at the null point and its associated periodicity. The resistive, nonlinear MHD simulations are solved in 2.5D for different levels of resistivity. Three methods (wavelet analysis, Fourier transform, and ANOVA) are used to investigate the effect of resistivity versus resultant period. It is found that there is an independence between the level of background resistivity and the period of the oscillatory reconnection mechanism. Conversely, it is found that resistivity has a significant effect on the maximum amplitude of the current density and the nature of its decay rate, as well as the magnitude of ohmic heating at the null.


INTRODUCTION
The solar magnetic field is a highly dynamic system, with motions constantly creating waves and flows within the field.Furthermore, a large number of regions with different magnetic polarities can lead to the formation of null points.Magnetic null points are singularities within the magnetic field where the field strength, and thus the Alfvén speed, are zero.Studies have shown, through the use of magnetic field extrapolations from magnetogram observations, that magnetic null points are ubiquitous within the solar atmosphere (e.g.Longcope 2005;Régnier et al. 2008).Due to this ubiquity, a large number of studies have been made on the role of null points, in particular how they interact with the motion of magnetohydrodynamic (MHD) waves through wave propagation (e.g.McLaughlin & Hood 2004, 2005, 2006a;McLaughlin et al. 2011;Sabri et al. 2022Sabri et al. , 2023) ) and mode conversion (e.g.McLaughlin & Hood 2006b;Thurgood & McLaughlin 2012;Tarr & Linton 2019;Pennicott & Cally 2021;Yadav et al. 2022) but also as sources of slow, fast and Alfvén waves (e.g.Lynch et al. 2014;Karpen et al. 2017;Thurgood et al. 2017).
Magnetic reconnection is a fundamental process in plasma physics where the topology of the magnetic field is rearranged.Strong currents are produced during reconnection events, which, due to electrical resistivity, allow magnetic fields lines to diffuse and change their connectivity (e.g.Parker 1957;Sweet 1958;Petschek 1964).In 2D, reconnection is understood to occur mainly at null points where it can be driven by a number of processes such as variations in the current density (Priest & Forbes 2000).Reconnection events can lead to heating of the surrounding plasma through ohmic heating and shock heating, as well as causing particle acceleration and localized motion.3D reconnection is a fundamentally different system and readers are directed to Priest et al. (2003), Pontin (2011) and Pontin & Priest (2022) for comprehensive reviews.
Generally, reconnection is studied using steady state models, either the Sweet-Parker model (Parker 1957;Sweet 1958) or the Petschek model (Petschek 1964).However, magnetic reconnection is a leading theory behind the origin of dynamic solar events such as solar flares (e.g.Shibata 1999;Su et al. 2013) and thus a time-dependent form of reconnection is more applicable to such dynamic events.One such time-dependent form is that of oscillatory reconnection (OR).
OR was first studied by Craig & McClymont (1991) who found that a null point, when perturbed from equilibrium and subsequently allowed to relax, could produce a reconnecting current sheet which oscillated about its original equilibrium orientation.It was then found by McLaughlin et al. (2009) that OR could be driven by a nonlinear fast magnetoacoustic wave which shocks and then collapses the null.OR has been studied in 3D by Thurgood et al. (2017) and more recently it has been found that it can be driven through the merging of magnetic flux ropes (Stewart et al. 2022).Karampelas et al. (2022a) simulated OR within 1 MK coronal plasma.
In particular, OR has become one of the leading contenders for the origin of quasi-periodic pulsations from solar flares (McLaughlin et al. 2018;Zimovets et al. 2021;Kou et al. 2022).OR has also become a potential theoretical explanation for the origin of other periodic solar phenomena, including periodicities within solar jets Figure 1.Contours of v ⊥ at times t = 0, 2.0, 2.7, 6.9, 12.5 and 60. t = 0 shows the initial pulse of v ⊥ (blue-to-red contour) and the equilibrium magnetic field, defined by Equation ( 7), where the black lines denote the magnetic field lines and the orange lines denote the seperatrices.The arrows indicate the direction of the magnetic field.t = 2.0 shows the formation of shocks on the wave pulse, which bend the magnetic field lines away from the normal.t = 2.7 shows the formation of the first current sheet.t = 6.9 and t = 12.5 show the change in the orientation of the current sheet.The black lines denote the evolving separatrices.t = 60.0 is the final snapshot from the simulation where the green line denotes the original separatrices at t = 0. (Hong et al. 2019), quasi-periodic flows (McLaughlin et al. 2012), periodicities within Type III radio bursts (Cattell et al. 2021), filament formation (Sun et al. 2023), flux rope formation and disappearance (Xue et al. 2019), stellar flare oscillations (Doyle et al. 2018), and a possible explanation for spikes and switchbacks observed by Parker Solar Probe (He et al. 2021).
In recent years, a number of studies have investigated the OR mechanism further to understand how the process is affected by various changes in parameters.In particular, Karampelas et al. (2022a,b) carried out a comprehensive study into how the period of the OR mechanism was affected by a number of environmental variables, such as temperature, density, background magnetic field strength, as well as the properties of the initial wave pulse.These studies led to the derivation of an empirical formula for the period of the OR oscillations in the electric current density and also led to the possibility of using OR as a plasma diagnostic in the corona (Karampelas et al. 2023).
In this paper, we will use a similar setup to that of McLaughlin et al. (2009) and study how variations in the levels of background resistivity affect the OR mech-anism.We shall focus on how a change in resistivity, across multiple orders of magnitude, affects the current density produced at the location of the magnetic null point.In Section 2 we describe our numerical approach to this study, outlining our physical domain, numerical setup and the code used.The results of this parameter study are presented in Section 3, with a focus on the effect on the amplitude of the oscillations of the current density (Section 3.2), the effect on the period in the system (Section 3.3), the effect on the decay rate of the current density oscillations (Section 3.4) and the corresponding ohmic heating produced (Section 3.5).In Section 4 we present our conclusions and general discussion resulting from this study.Appendix A reports on an additional study on the period of the system using Fourier Spectra, with Appendix B showing the results of simulations complementing those in Section 3.

Governing Equations
In this study the Lare2d code (Arber et al. 2001) is used to solve the resistive magnetohydrodynamic (MHD) equations in 2.5D, which are defined as: where ρ is density, P is plasma pressure, v is velocity, B is the magnetic field, E is the electric field, j is the electric current density, µ 0 is the magnetic permeability, ϵ is the specific internal energy density and η is resistivity.In these equations D Dt is the advective derivative.These equations can then be non-dimensionalised by letting: v = v o v * , B = B 0 B * , x = Lx * , y = Ly * , j = j 0 j * , ρ = ρ 0 ρ * , P = P 0 P * , E = E 0 E * , η = η 0 η * and t = t 0 t * .Here * denotes the dimensionless parameter and v 0 , B 0 , L, j 0 , ρ 0 , P 0 , E 0 , η 0 and t 0 correspond to constants with the dimensions of the variable they are non-dimensionalising.Due to the scale free nature of our equilibrium magnetic field (Equation 7), there is freedom in the choice of these constant values.For the remainder of this paper, the * is dropped from notation, with all variables being non-dimensionalised.Variables, such as velocity and time, can be re-dimensionalised through multiplication with their respective dimensional constants, such as v 0 = B 0 / √ µ 0 ρ 0 and t 0 = L/v 0 .
The Lare2d code utilises artificial shock viscosities (ν 1 and ν 2 ), which introduce a level of dissipation at steep gradients, allowing the accurate capture of MHD shocks (see Caramana et al. 1998, for details).The values that are chosen in this study for these viscosities are ν 1 = 0.5 and ν 2 = 0.1.

Initial Conditions
The computational domain is defined as [x, y] ∈ [−97.5, 97.5] with a numerical resolution of 2560 × 2560 grid cells.A stretched grid is implemented within this domain in both the x and y directions, allowing for a uniform numerical resolution of 1280 × 1280 within the domain [x, y] ∈ [−5, 5].The implementation of this stretched grid allows for the uniform grid to be around the area of interest, the null point, but allows for the boundaries to be a significant distance away.
The system is initially defined as a cold plasma with a uniform density where the equilibrium magnetic field is defined similar to that of McLaughlin et al. (2009) but with a π/4 rotation: which is defined so that the solenoidal condition ∇ • B = 0 holds.It is worth noting here that this magnetic field becomes unphysical far from the null point due to the magnetic field strength tending towards unphysically large values.It has been shown in McLaughlin & Hood (2006a) that for equilibrium magnetic fields where this unphysical nature at large distances is removed, the key underlying physics near the null still hold.From this magnetic field, the vector potential A can be calculated, where B = ∇ × A. In the 2D coordinate system chosen, this equation becomes A = A z ẑ, with the equilibrium vector potential defined as: This initial magnetic field can be seen in Figure 1 for t = 0.0 along with its subsequent evolution.
In this system two main parameters related to the velocity are of interest, namely the parallel, v ∥ = v • B, and the perpendicular, The initial wave pulse is defined in the perpendicular component, v ⊥ , as a circular pulse around the null, as in McLaughlin et al. (2009) and Karampelas et al. (2022b).This initial condition, coupled with setting the parallel component v ∥ to zero, results in the initial wave pulse being a purely fast magnetoacoustic wave, defined by: where r = x 2 + y 2 .This initial wave pulse can be seen in Figure 1 for t = 0. Due to the definition of the initial wave pulse in this way, it will result in the wave splitting into two wavesone of these waves will move towards the null, with the other moving away -each with an amplitude of C.This study is solely focused on the incoming wave, with the outgoing wave being numerically removed via a damping region which is discussed in Section 2.3.
The choice for the value of constant C dictates at which distance from the null the wave pulse will begin to shock.In this system, a value of C = 1 is chosen, as per McLaughlin et al. (2009), which allows for the initial wave pulse to develop a shock before it reaches the null point.

Boundary Conditions and Damping Region
Zero gradient boundary conditions are implemented for velocity, magnetic field, density, temperature and energy.Due to the Alfvén speed increasing linearly with distance away from the null, waves are accelerated as they approach the boundary before reflecting.Since re- flected waves would interfere with the long-term evolution at the null, the damping of these oscillations and reflections is crucial for this study.Therefore a damping region is implemented to prevent these reflections from interfering with the region of interest around the null point.
As a result of the stretching of the grid in the numerical domain, a large damping region can be implemented due to the boundaries being positioned a large distance from the region of interest.This damping region is implemented for r ≥ 10 where, using a combination of viscosity and artificial kinetic energy removal, oscillations and reflections that enter this region are slowly damped away, never to return to the central region of interest.
The use of viscosity as a way of damping away the outgoing oscillations (for r ≥ 10) is implemented by introducing a viscosity term, ν∇ 2 v, into the momentum equation ( 2), where ν is the viscosity term.This term is highly spatially dependent (e.g. it is zero for r < 10) and is defined such that its value increases as it approaches the boundary.
The artificial removal of kinetic energy is implemented by reducing the amplitude of the oscillations within the damping region for each iteration of the Lare2d code.This is implemented using the equation: v = v/(1 + a), where v is velocity and a is a small amplitude defined by a number of hyperbolic tangent functions which are radially dependent.Due to the large size of the damping region, this process successfully damps away outgoing oscillations and any reflections.
It is important to note here that these damping techniques are not implemented within the uniform grid re- gion around the null.They only occur within the damping region (r ≥ 10).

Fast Magnetoacoustic Wave Driven Oscillatory Reconnection
Let us now consider a simulation where the value of the resistivity is set to η = 10 −4 , which we will refer to as the "baseline" simulation.It is noticed that due to the definition of the initial condition, Equations ( 9) and ( 10), the wave pulse splits into two oppositelypropagating waves, both with amplitude C. As discussed previously, the outer-propagating wave will move into the damping region away from the null where it is removed from the simulation.The inwardly-propagating wave pulse moves towards the null and develops shocks on either the leading or trailing edge, depending on the region of connectivity the wave is propagating through.As these shock fronts steepen, they begin to bend the magnetic field away from the normal and eventually collapse the null point into a current sheet (Figure 1 for t = 2.0 and t = 2.7 respectively).As this current sheet is formed, jets are formed at the ends of the sheet, which can be seen in Figure 1 for t = 2.7.We shall refer to this current sheet orientation, at approximately 135°, as orientation 1.
Subsequently, the length of the current sheet begins to shrink before changing orientation by 90°(which we will refer to as orientation 2), as seen in Figure 1 for t = 6.9.This process repeats, leading to an oscillating current sheet, which is the main characteristic of OR.McLaughlin et al. (2009) studied the mechanism for this oscillation in the current sheet orientation including how shocks form in the jets and affect the heating of the surrounding plasma.
This oscillation of the current sheet can be seen by calculating the j z component of the current density at the null point using Equation (6).From Figure 2 it can be seen that the current density remains at zero until a time of t ≈ 2.7, at which point the null point collapses as the first current sheet is formed (Figure 1 for t = 2.7).The value of j z drops considerably at this point, with the negative value corresponding to a current sheet of orientation 1.The value of the current then changes to become positive at a time of t ≈ 5, which corresponds to the formation of the current sheet in orientation 2 (seen in Figure 1 for t = 6.9).
It can also be seen that the oscillations in j z decay over time.This is due to the gradual shortening of the current sheets after each oscillation due to the magnetic field slowly relaxing back to the equilibrium state.It is noted here that the oscillations of j z do not oscillate around a zero value, but instead oscillate around a value of j z ≈ 0.87, as is shown by the blue dashed line in Figure 2.This is due to residual current den- Resistivity Level Order of Fit 0 3rd 10 −5 4th 10 −4.5 2nd 10 −4 2nd 10 −3.5 4th 10 −3 4th 10 −2.5 4th 10 −2 4th 10 −1.5 3rd sity that remains in the system, which results from the magnetic field being very slightly non-potential at the end of the simulation (due to localised heating from the periodic, orientation-changing reconnection jets).This can be seen in Figure 1 for t = 60, where there is a slight difference between the magnetic field line (black) and the t = 0 separatrices (green).

Effect on Amplitude of j z Oscillation
Firstly, the focus will be on the amplitude of the current density produced at the null, in particular the current density produced at the formation of the first current sheet, which corresponds to the first local minimum in the current density.By comparing the values of j z (0, 0, t) at the formation of the first current sheet for each of the simulations, it can be understood how resistivity affects the magnitude of the current density produced through the OR mechanism.A comparison of these values can be seen in Figure 3.
From Figure 3, it can be seen that there are three different regions of interest.The first is the region consisting of points to the left-hand-side of the plot (points for η = 10 −8 to 10 −5 ).In this region, hereby referred to as region I, it is found that a change in the background resistivity has little, if any, effect on the amplitude of the maximum absolute value of the current density.We interpret this is due to the simulations not being of high enough resolution to resolve these low levels of resistivity, and we conclude that they instead are dominated by numerical resistivity.This is confirmed when comparing these values with the value from the ideal simulation (η = 0).Due to this, no conclusions will be made based solely on these simulations.
Secondly, there is a region in the center of the plot for levels of resistivity (points for η = 10 −4.5 to 10 −2 ).In this region, hereby referred to as region II, we can see a clear linear trend between resistivity and the amplitude of the |j z (0, 0, t)| values.
The third region is that for high levels of resistivity (points for η = 10 −1.5 to 10 −1 ).In this region, hereby referred to as region III, one can see the levels of current density flatten out as they approach zero.This trend is a trivial one as the value of |j z (0, 0, t)| cannot drop below zero due to us taking the absolute value of the current.
A comparison of the first local minima can be seen in Figure 4, where the effect of resistivity on the amplitude of the j z (0, 0, t) can be seen.

Effect of Resistivity on Periodicity
The period of the OR system can be found by focusing on the oscillating j z (0, 0, t) signal seen in Figure 2. In the next subsections, two methods for analyzing the period of the system are discussed, along with the comparison for both the values and evolution of the period for different values of resistivity.The first method that will be used to analyze the period for the OR system is that of wavelet analysis, where for this analysis the Morlet wavelet is employed.

Wavelet Analysis
The results of this analysis is a contour of power against time and period (Figure 5: middle).In this contour, the dark red [blue] signifies a high [low] power.Notice also on this plot the regions for the cone of influence (COI), defined by the white slanted lines, and the 60% confidence interval (CI), defined by the white contour encompassing the main peak.
The period of the system can be found by analyzing the results of the wavelet power contour within the region between the COI and the 60% CI.
This wavelet analysis begins by taking cuts of the wavelet power contour at every time step.This gives a profile of power against period for each time step.To find the predominant period at each cut a Gaussian is fitted to the profiles, with the center µ being the predominant period, and its standard deviation σ.Plotting these predominant periods (blue points) and their respective standard deviations (orange bars) against time gives the plot seen in Figure 5: bottom for η = 10 −4 .(The wavelet plots of the other simulations of varying levels of η can be found in Appendix B, Figure 14) From this plot, the period of the system can be calculated by taking the average value of the period within the region between the COI and the CI.The method of Propagation of Error (PoE) is used to incorporate all of the standard deviations (between the COI and CI) in Figure 5: bottom to calculate the error on this period.This process is then repeated across each of the simulations for different levels of resistivity resulting in Figure 6 which shows the respective periods (blue points) and attributed errors (red error bars).The black dashed lines split this plot into three regions as in Figure 3.
From Figure 6, it can be seen that any change in period when resistivity is changed is negligible within the error bars.Therefore, we can fit a zeroth order trend to this plot, p = 8.50, shown by the blue line in Figure 6: left.This fit is calculated only using points from region II, η = 10 −4.5 to η = 10 −2 , due to the size of the errors for the points where η = 10 −1.5 and η = 10 −1 .
To calculate a resultant uncertainty of this constant period, PoE is again utilized to incorporate all the individual errors on the period values for each simulation.This uncertainty can be seen via the blue dashed lines in Figure 6.These calculations yield a value for the period of p wavelet = 8.50 ± 0.89, irrespective of the value for the resistivity.
To corroborate these results, Fourier transforms of the j z (0, 0, t) profiles were calculated with the results found in Appendix A.

Analysis of Variance and Prediction Intervals
Section 3.3.1 yields values for the period with the assumption that this period is constant throughout the duration of the simulations.However, in Figure 5: bottom, there appears to be a change in the period with respect to time.To test the significance of this change, Analysis of Variance (ANOVA) tests (Scheffe 1999) can be carried out by testing the significance of higher-order terms in a polynomial fit.This results in the lowestorder polynomial which accurately fits to the data to be found.A visual representation of this ANOVA test on the baseline simulation can be seen in Figure 7.
This technique is carried out on the plots of period against time in Figure 5: bottom for each of our simulations, with the resulting order of polynomial fit for each simulation given in Table 1.Due to the high levels of decay for the η = 10 −1 simulation, we do not include this simulation in this analysis as the errors from the wavelet power contour are too large.Due to the dominance of numerical resistivity for simulations within region I, these simulations are also omitted from this analysis.
It can be noticed from Table 1 that despite changing the level of resistivity across multiple orders of magnitude, the evolution of the periods of each simulation can be described by similarly ordered polynomials.As a further test of these fits, the prediction intervals of the fitted lines against the data are calculated.To incorporate the error bars calculated in Figure 5: bottom, the Von Neumann rejection method is used to produce a Gaussian distribution of 1000 random points at each time step.These Gaussian distributions take the values for parameters µ and σ, directly from the values in the bottom panel of Figure 5 (i.e.Section 3.3.1).The prediction intervals are then calculated using these distributions.
The prediction interval was calculated using the inbuilt function "Predict" in R, with the underlying statistics calculated by the equation: where y p± is the prediction interval, ŷh is the fitted data, n is the number of data points, M SE is the mean square error, x k is the predictor value, x is the mean value and t a/2,n−2 is the critical t-value where a is the confidence level with n − 2 degrees of freedom.In this study we take the 95% prediction interval.
Figure 8 shows the prediction intervals for each simulation.In this plot, the range of time is fixed to that of the shortest range between the COI and 60% CI (Figure 5: middle) across the simulations.This is done so that a direct comparison of the evolution for the period of each simulation can be made.It can be seen from Figure 8 that the fits for each simulation (solid lines) lie within the prediction interval of each simulation (dashed lines).This suggests that the variation in the fits across the simulations is negligible in relation with the attributed errors.Due to this it can be stated that none of the models are significantly different from one another.
3.4.Effect on Decay Rate of j z (0, 0, t) Oscillation A further study can be made into how resistivity affects the decay rate of the j z oscillation at the null.This is done using a similar method to that used in Karampelas et al. (2022b), namely considering the local maxima of the oscillations in j z (0, 0, t) and finding the rate at which these maxima decay throughout the oscillation.
To find the decay rate of these maxima, the natural logarithm (ln) is taken of these local maxima and then a linear trend fitted, where the linear trend is simply defined as y lin = at + b, where a is a decay constant and b is a constant.This fit can be seen in Figure 9: bottom row, where the values of the decay constants can be found in the legend with their corresponding errors.
This fitting is repeated for each of the simulations, where for each simulation, the number of local maxima to which the fit is made is varied, to investigate the effect that this has on the decay rate values calculated.The plots of the decay rates for each simulation can be found in Appendix B, Figure 16.The decay rates for each sim- ulation can be seen in Figure 10.In this plot, decay rate values calculated from 5, 6 and 7 peaks are indicated by the triangle, cross and circle markers respectively, with the same separating lines (black dashed) as in Figure 3 used, with the resistivity regions labeled in the same fashion.

Effect on Ohmic Heating
One effect that the generation of current density can have on a system is heating, in particular ohmic heating, which is directly related to the currents produced and the levels of resistivity in the system, with the relation defined by ηj 2 .
Figure 11 shows the evolution of ohmic heating at the null for each simulation.It can be seen from this plot, due to direct correlation with the oscillations in j z (0, 0, t), that the level of ohmic heating also oscillates with the formation of each subsequent current sheet.The lines for η = 10 −8 , 10 −7 and 10 −6 are not included in this figure due to them lying on top of the line of η = 10 −5 , however their points are included in Figure 12.
For a direct comparison, the maximum value of the ohmic heating for each simulation is taken and compared.Just like for the amplitude of j z (0, 0, t) study previously, the sampling of the maximum value corresponds to the formation of the first current sheet.The comparison for this maximum value can be seen in Figure 12.

CONCLUSIONS
In this paper a study into the OR mechanism in a 2D, initially cold plasma is made, developing on previous works by investigating the effect of resistivity on the system.
By solving the resistive MHD equations, using the Lare2d code, the level of background resistivity is varied across eight orders of magnitude (η = 10 −8 to η = 10 −1 , with also an ideal MHD case where η = 0), allowing a study of how this variation affects key characteristics of the OR mechanism such as the current density produced, ohmic heating and period of the oscillating system.The choice of initial condition in this paper consists of a circular wave pulse surrounding a null point, such as in previous works by McLaughlin et al. (2009) and Karampelas et al. (2022aKarampelas et al. ( , 2023)), which initiates OR.
For each of these simulations, the current density profile at the location of the null point, j z (0, 0, t), were calculated with the maximum absolute value of the amplitudes of the current oscillations compared.It was found that for "low" levels of resistivity (region I of Figure 3), there was little or no change in the amplitude of |j z (0, 0, t)|.We explain this by concluding that this is due to the numerical resistivity of the simulation dominating the system in this region.This can be seen by the straight line plotted through these points.This interpretation is also corroborated by the good correspondence between this line and the value for η = 0, where the resistivity in the ideal simulation is of course solely numerical.Due to this, no conclusions can be made from these simulations.
It is found that for 'high' levels of resistivity, region III in Figure 3, the level of the current flattens out as the absolute value approaches zero.This is shown by the green, decaying exponential fit in Figure 3.
In region II, the behavior is not affected by the numerical resistivity, nor is it affected by the current levels nearing zero.It is found that there is a direct relationship between resistivity and the amplitude of the currents produced in the system.This is indicated by the linear fit to the points in region II (light green line in Figure 3).
From the current density profiles, the period of the oscillations within the profiles were calculated using wavelet analysis.It was found that there is no significant change in the period when resistivity is increased across multiple orders of magnitude, suggesting an independence between resistivity and the periodicity of OR for our choice of parameters and assumptions.It is found that there is a good correspondence between the constant period values calculated from the wavelet analysis and Fourier transform, i.e. p wavelet = 8.50±0.89and p fourier = 8.59 ± 1.50.
These values for the period are in computational units due to the system being non-dimensionalized.However, one can re-dimensionalize these values by multiplying them by the dimensional constant t 0 = L/v 0 , where G and ρ 0 = 10 −12 kg m −3 for the characteristic length, magnetic field strength and density respectively, we find t 0 = 7.78 s.This value gives a period for our simulations of p wavelet = 66.13 ± 6.92 s and p fourier = 66.83 ± 11.67 s, calculated using wavelet analysis and Fourier transforms respectively.However, it is important to note that Equation ( 7) is scale free.Therefore there is freedom in the choice of values when re-dimensionalizing.
In addition, a statistical approach was also used through ANOVA and a prediction interval calculation to study the evolution of the period in the current density profiles from the wavelet analysis.It was found that the evolution of the period for each simulation can be described by similarly ordered fits.When the prediction intervals for each simulation were calculated, based on these fits, it was found that each evolution lies within the prediction interval for all other fits.Therefore we conclude that these fits are not significantly different from the other.This again suggests that there is an independence between resistivity and periodicity.
We conclude that there is an independence between resistivity and the periodicity of OR.
This result of an independence between resistivity and periodicity of OR has to be discussed in the context of previous studies, specifically that of Craig & McClymont (1991).This result seems to contradict the result of their paper.In their study they found a direct link between resistivity and period of the oscillating current sheet.They found the frequency of the oscillation for the fundamental mode to be defined by ≈ 2 ln S where S = 1/η is the Lundquist number.However, it is important to note that in Craig & McClymont (1991), a linearized system within a closed cylindrical domain was solved, whereas in the study presented in this paper the OR mechanism is investigated using the nonlinear MHD equations with effective open boundaries.Thus we conclude that the study presented here and that of Craig & McClymont (1991) are investigating fundamentally different systems.
From the current density profiles for each of the simulations, an investigation was made into how resistivity affects the decay rate for the oscillations of j z (0, 0, t).By finding the local maxima within the profiles, the decay rate was calculated by taking the natural logarithm (ln) of the maxima and making a linear fit to those points.In Figure 10, it was found that the decay rate increases as the level of resistivity increases, with the exact nature of this increase being defined by an exponential function (solid black line).In the region where numerical resistivity dominates (region I), no fit was made.When calculating the relationship between resistivity and decay rates, only the values in region II were fitted.It was found that the decay rates calculated for 'high' levels of resistivity in region III flatten out suggesting a possible maximum decay rate for the oscillations.
The ohmic heating produced at the null point was calculated, with the oscillations in the levels of ohmic heating for each of the simulations being directly linked to the oscillations in the current density (Figure 11).For a more direct comparison, the maximum value of this ohmic heating for each of our simulations was taken.
Similar to previous sections, Figure 12 is split into regions I, II and III.Region I, where numerical resistivity dominates, demonstrates a linear increase in ohmic heating as resistivity increases.This linear trend is due to fact that, as seen in the previous sections, we see little change in the levels of current density in the system for different choices of η.Therefore the increase that we see in the ohmic heating is solely due to the increase in the level of the background resistivity.
In region II, a distinct quadratic trend between ohmic heating and resistivity is found.This trends suggests that there may be an "optimal" value of resistivity in terms of maximizing ohmic heating.
This quadratic trend also extends to region III, and follows from the study on the amplitude of the current (Section 3.2), where for high resistivity levels we see a high level of damping on j z (0, 0, t) (Figure 3).This strong effect on the current in turn will have an effect on the ohmic heating, which can be seen in Figure 12.
To summarise, this study shows that there is an independence between background resistivity and the periodicity of the oscillatory reconnection mechanism for the parameter values chosen.However, it is also found that for values above the level of numerical resistivity, background resistivity has a distinct effect on the amplitude and decay rate of the current density oscillations at the null, as well as the levels of ohmic heating.  5.For each panel, the top panel shows the evolution of jz(0, 0, t).Middle panel shows the wavelet power contour, where dark red denotes a high power and dark blue denotes low power.The white lines spanning from (t = 0, p = 0) to (t = 20, p = 20) and (t = 40, p = 20) to (t = 60, p = 0) denote the cone of influence and the white contour encompassing the main peak in the spectrum denotes the 60% confidence interval.The bottom panel shows the predominant period at each time-step (dark blue) with their attributed standard deviation (orange).
Similar to Section 3.3.1, the uncertainty of this fit is calculated using PoE to incorporate the error bars in Figure 13: right, and are shown by the blue dashed lines.This results in a value for the period of p fourier = 8.59 ± 1.50, which shows a good agreement with the period calculated from the wavelet analysis.Figures 14,15,and 16 show the plots for the wavelet analysis, Fourier transform, and decay rates respectively for each simulation of the parameter study.Similarly, η = 10 −4 is not shown due to being available in Figure 13.The mean, μ, and standard deviation, σ, are given in the legend for each panel.

Figure 2 .
Figure2.Plot of jz(0, 0, t) at the null point against time for baseline simulation (η = 10 −4 ).The green dashed line denotes jz = 0 with the blue dashed line denoting jz = 0.87, which is the value that the jz profile tends to.

Figure 3 .
Figure 3. Plot of resistivity vs maximum amplitude of |jz|, with the black dashed lines giving approximate locations for different trends and splitting the plot into three regions, labelled I, II and III respectively.The red line shows a flat trend for points where numerical resistivity dominates.The light green line gives a linear trend for intermediate resistivity values (in region II).The dark green line shows the trend for high resistivity values as the amplitude of the current approaches zero.

Figure 4 .
Figure 4. Overplot of jz(0, 0, t) oscillations, showing effect of resistivity on amplitude of current density produced at the formation of the first current sheet.The line for each value of η is given in the legend.

Figure 5 .
Figure 5. Top panel shows the jz(0, 0, t) oscillation for η = 10 −4 .Middle panel shows the wavelet power contour, where dark red denotes a high power and dark blue denotes low power.The slanted white lines denote the cone of influence and the white contour encompassing the main peak in the spectrum denotes the 60% confidence interval.The bottom panel shows the predominant period at each time-step (dark blue) with their standard deviation (orange).

Figure 6 .
Figure 6.Resistivity vs period from wavelet analysis.The blue dots indicate the mean period values with the red bars denoting their attributed errors.The blue line indicates a mean period calculated from the period values for η = 10 −4.5 to η = 10 −2 with the blue dashed lines indicating the uncertainty on this value calculated through Propagation of Error.The blue star indicates the results from the ideal simulation.The black dashed lines split the plot into three regions for low (I), intermediate (II) and high (III) levels of resistivity as in Figure 3.

Figure 7 .
Figure 7.The analysis of variance study for the baseline simulation (η = 10 −4 ).The dark blue dots show the predominant period at each time step with the orange error bars taken from the bottom panel of Figure 5.The lime, cyan and magenta lines denote the 0th, 1st and 2nd order fits respectively.

Figure 8 .
Figure 8. Prediction interval fits for simulations where η = 10 −5 to η = 10 −1.5 as well as the ideal case (η = 0).The solid lines denote the predicted values with the dashed lines showing the attributed errors.The colors for each simulation can be found in the legend.

Figure 9 .
Figure9.Calculation of the decay rate of the jz(0, 0, t) oscillation when η = 10 −4 .In the top row, the blue points correspond to the location of the local maxima of the jz(0, 0, t) profile (solid black line).The bottom row shows the values of the local maxima but on a natural logarithmic scale (blue points) with a linear fit to these values for 5, 6, and 7 local maxima shown by the magenta, cyan and lime lines respectively.. The legend shows the fitted decay constant of the linear fit with the attributed error.

Figure 10 .
Figure 10.Resistivity vs. decay rate taken from linear fit to log-value plots (Fig. 9 bottom row).The plot is split into three regions (I, II and III), as in Figure 3, shown by the black dashed lines.Each simulation has multiple points depending on how many local maxima were used when fitting the data, indicated in the legend.The black solid line is an exponential fit, fitted to the points in the intermediate resistivity range (region II) for the decay rates calculated using 6 local maxima.

Figure 11 .
Figure11.Time evolution of ohmic heating at the null for different levels of resistivity, presented as log 10 (ηj 2 z ).The lines for each simulation can be found using the legend.

Figure 12 .
Figure 12.Resistivity vs. maximum ohmic heating.The black dashed lines split the plot into three regions for low (I), intermediate (II) and high (III) levels of resistivity as in Figure 3.The solid red line shows a linear trend, fitted to the values for low resistivity where numerical resistivity dominates.The green line shows a quadratic trend for the intermediate and higher values of resistivity.

Figure 14 .
Figure14.Panels showing wavelet analysis of current density signals, jz(0, 0, t) at the null for increasing levels of background resistivity, from η = 0 and η = 10 −7 to η = 10 −1 .Results for η = 10 −8 , 10 −7 and 10 −6 are not shown due to their similarity with η = 0 and 10 −5 .Similarly, η = 10 −4 is not shown due to being available in Figure5.For each panel, the top panel shows the evolution of jz(0, 0, t).Middle panel shows the wavelet power contour, where dark red denotes a high power and dark blue denotes low power.The white lines spanning from (t = 0, p = 0) to (t = 20, p = 20) and (t = 40, p = 20) to (t = 60, p = 0) denote the cone of influence and the white contour encompassing the main peak in the spectrum denotes the 60% confidence interval.The bottom panel shows the predominant period at each time-step (dark blue) with their attributed standard deviation (orange).

Figure 15 .
Figure15.Panels showing the Fourier spectra for each simulation (black line) that is normalized to the maximum power, with a fitted Gaussian (red line) to the main peak of the spectra, for increasing levels of background resistivity, from η = 0 and η = 10 −5 to η = 10 −1 .Results for η = 10 −8 , 10 −7 and 10 −6 are not shown due to their similarity with η = 0 and 10 −5 .Similarly, η = 10 −4 is not shown due to being available in Figure13.The mean, μ, and standard deviation, σ, are given in the legend for each panel.

Table 1 .
Order of fitted polynomials to period vs time plots from ANOVA results.