## 3 Data and Methodology  
For the statistical correlations and the distribution fitting investigations on the 10 sub-region seismicity levels, we selected an earthquake catalog comprising the \(M_{w}\geq 4.1\) events that occurred in the territory of Greece during 1975-2021, which we have taken from the regional catalog compiled by the Geophysics Department of the Aristotle University of Thessaloniki (Available online: [[http://geophysics.geo.auth.gr/ss](http://geophysics.geo.auth.gr/ss)]([http://geophysics.geo.auth.gr/ss](http://geophysics.geo.auth.gr/ss)) (accessed on 30 April 2022)). The earthquake magnitudes are expressed in moment magnitude scale (\(M_{w}\)), as obtained directly from waveform modeling or equivalent \(M_{w}\) based on scaling relations suggested by [PERSON] et al. [29]. We divided the catalog into 10 sub-catalogs corresponding to the 10 distinctive sub-areas. Figure 2 shows the epicentral distribution of the earthquakes comprised in each sub-catalog with different colors for each sub-area for the reader's ease. For each sub-catalog the samples of the earthquake intercept times, \(\Delta T\), data are then defined for earthquakes with magnitude thresholds of 0.1 bin increment, starting from \(M_{w}\)\(\geq\)4.1 (\(\Delta T_{1}\), \(\Delta T_{2}\),... for the magnitude thresholds of \(M_{w}\geq 4.1\), \(M_{w}\geq 4.2\),... and so on for each sub-area). These data sets are the inputs for applying the methodology, which is described below. We seek for the threshold magnitude, \(M_{thr}\), of a certain data set, above which the interevent times are proved to be independent. For this purpose, we investigated the correlation of \(\Delta T\) in each subset, with the subsets being created for certain magnitude bins. We examined the correlations by calculating the Autocorrelation (_ACF_) and Partial Autocorrelation (_PACF_) functions, a widely used approach for the same purpose (e.g., [30] for global earthquake data; [31] in induced seismicity of Poland). The Autocorrelation Function (_ACF_),  
\[\rho_{x}=\frac{\sum_{i=1}^{n-k}(x_{i}-\overline{x})(x_{i+k}-\overline{x})}{ \sum_{i=1}^{n}(x_{i}-\overline{x})^{2}} \tag{1}\]  
where \(n\) is the number of the observations, \(k\) is the number of the lags and \(\overline{x}\) is the mean value of the observations, and is used for the investigation of correlations between past and future values of a given time series by assuming a confidence level (the 95% confidence level in the current application). If a certain \(\rho_{x}\) value at lag \(k\) exceeds the \(\pm\)95% confidence level, the 95% confidence level is \(\pm\)95% confidence level. interval then the given time series could be considered as correlated for the specific lag, otherwise, it could be considered statistically independent. Then, the possible correlations detected by the _ACF_ can be confirmed in the same way based on the values of the Partial Autocorrelation Function that are given by  
\[r_{k,k}=\frac{\rho_{k}-\sum_{i=1}^{k-1}r_{k-1,j}\rho_{k-j}}{1-\sum_{i=1}^{k-1}r _{k-1,j}\rho_{j}} \tag{2}\]  
where \(r_{k,j}=r_{k-1,j}-r_{k,k}r_{k-1,j-i}\) for \(j=1,2,\ldots,k\). Additionally, the Ljung-Box Q-test [32], which is a modified version of the Box-Pierce Portmanteau test, is applied as an alternative method to investigate correlations among \(\Delta T\). This test is implemented under the null hypothesis that the time series exhibits no autocorrelation for a fixed number of lags, \(L\), against the alternative that some (statistically significant) autocorrelation exists. The test statistic is given by  
\[Q=n(n+2)\sum_{i=1}^{L}\left(\frac{\rho_{k}^{2}}{n-k}\right). \tag{3}\]  
The asymptotic distribution of \(Q\) is Chi-square (\(x^{2}\)) with \(L\) degrees of freedom. If the statistic of the test for the given lag \(L\) is less than the critical value of the Chi-square (\(Q_{stat}<Q_{critical}\)), the null hypothesis cannot be rejected.