## 2 Materials and Methods
### Synchronization between Participants' HRV and Earth's Magnetic Field  
In order to assess the degree of synchronization between groups' HRV and local magnetic field spectral power (computed as described in previous subsection), techniques proposed and validated in [11; 44] were employed. The key concept of the applied approach is the existence of such mapping that allows to reduce data signal to one integer characterizing the dynamical properties of the analyzed series from the geometrical perspective. Namely, the integer corresponds to the time lag value that maximizes the area of the attractor reconstructed from the analyzed signal.  
The procedure denoted in [11] as Algorithm C was slightly modified in order to yield a single scalar value corresponding to the synchronization between analyzed data signals. The short outline of the modified algorithm is presented below.  
Consider two synchronously sampled data series \(RR=(RR_{1},\ldots,RR_{n})\) and \(M=(M_{1},\ldots,M_{n})\) corresponding to participant's HRV (RR) and local magnetic field (M) spectral power, respectively.  
* Data signals RR and M are divided into \(T\) 5-min-long observation windows (according to HRV analysis standards) [36] and optimal time lag values are computed for each observation window (using the procedure presented in [11] and summarized above) for HRV time series \(\left(\tau_{*}^{(RR)}=\left(\tau_{*1}^{(RR)},\ldots,\tau_{*T}^{(RR)}\right)\right)\) and magnetic field spectral power data signal \(\left(\tau_{*}^{(M)}=\left(\tau_{*1}^{(M)},\ldots,\tau_{*T}^{(M)}\right)\right)\). This procedure maps RR and M data series to integer vectors, describing dynamical and geometrical properties of those series.
* Obtained optimal time lag vectors are smoothed via the computation of moving average in order to identify averaged changes in time lags for HRV \(\left(\overline{\tau}_{*}^{(RR)}\right)\) and  
Figure 1: The original \(D_{i}\) (blue line) and smoothed \(\widehat{D_{i}}\) (red line) discrete derivative of the difference between standard deviations of filtered and non-filtered data. The degree of tolerance used in this study was set to \(\varepsilon=0.01\) (black line). Thus, the optimal size of the filtering data window was selected as the maximal data window satisfying \(\widehat{D_{i}}>\varepsilon\). It can be observed from this figure that the optimal length filtering data window corresponding to such criteria is 36,000 elements, i.e., 10 h of local magnetic field spectral power data.  
magnetic field data \(\left(\overline{\tau}_{*}^{(M)}\right)\). The size of the observation window used for the moving average is 20 data points.
* Pearson correlation coefficient \(C^{(RR,M)}=\rho\left(\overline{\tau}_{*}^{(RR)},\overline{\tau}_{*}^{(M)}\right)\) is calculated between averaged optimal time lag vectors. Resulting value \(C^{(RR,M)}\) is considered as the degree of synchronization between analyzed HRV and magnetic field spectral power signals. Naturally, higher value of \(C^{(RR,M)}\) corresponds to a higher degree of synchronization and vice versa (\(C^{(RR,M)}\in[-1,1]\)). Note that the assumption of the normality of data series used in further analysis was verified.  
Let \(RR^{(k)}=\left(RR_{1}^{(k)},\ldots,RR_{n}^{(k)}\right),k=1,\ldots,K\) is a set of HRV data signals, where \(RR^{(k)}\) corresponds to \(k\)-th participant from one of the considered countries. The mean synchronization between groups' HRV and respective magnetic field spectral power signal is calculated as the mean correlation coefficient \(C^{(RR^{(1,K)},M)}=\frac{1}{K}\sum_{k=1}^{K}C^{(RR^{(k)},M)}\), where \(C^{(RR^{(k)},M)}\) is the synchronization between \(k\)-th participant's HRV and magnetic field data signals computed via the procedure described above.