## 1 Introduction  
A prominent feature of brain signals recorded extracellularly at high temporal resolutions is the existence of rhythmic oscillations in several frequency bands that span from approximately 0.05 Hz to 500 Hz [1]. Extensive research in animals and humans has provided mounting evidence of the involvement of these rhythms in various brain functions, from sensorimotor to cognitive ones [1, 2], as well as their alterations in neurological and psychiatric disorders [3, 4]. Brain oscillations exhibit a number of spatiotemporal relationships, including correlations of rhythms at different frequencies (cross-frequency coupling, CFC). One type of CFC is phase-amplitude coupling (PAC), i.e., the modulation of the amplitude of a high-frequency rhythm by the phase of a low-frequency oscillation. The earliest observations of PAC were in the rat hippocampus [5] between the phase of theta waves (5-10 Hz) and the amplitude of gamma activity (40-150 Hz), termed theta-gamma coupling (TGC). PAC was later reported between other frequency bands, in other animals and humans, in recordings by both invasive [6, 7] and noninvasive [8] techniques. PAC was found to be correlated with cognitive states and performance [9, 10] and altered in a number of brain disorders [11, 12]. This suggests that it may play an important role in the functioning of the brain. A prominent hypothesis about the role of PAC in working memory suggested in [13] is that each remembered item is encoded by the spatial pattern of a group of cells, and every such group synchronously fires in a different high-frequency subcycle of the low-frequency oscillation. This was supported by studies of grid cells in the entorhinal cortex [14; 15]. The model suggested an explanation of the classic working memory capacity limit (span) of \(7\pm 2\)[16] as being determined by the number of gamma cycles within a theta cycle, although the evidence about this relationship has been mixed [17; 18]. The implication of PAC in brain function motivates its use for informing neuromodulatory interventions, including non-invasive transcranial electrical and magnetic stimulation (TES and TMS) [1; 19]. Recent studies have used individual frequencies of maximal TGC as parameters in protocols for transcranial alternating current stimulation (tACS) [20] or TMS [21]. This use implies that the determined frequencies of maximal coupling, at least to a certain approximation, constitute stable characteristics of each individual subject. However, there is insufficient data on the degree of session-to-session reliability of these parameters. One study [18] reported that the gamma frequencies of maximal coupling in two sessions had a (relatively low) Spearman correlation of 0.32, \(p=0.03\) (n = 44), while the theta frequencies were not correlated. This raises concerns about the possibility of reliable TGC peak frequency estimation. However, the data in the mentioned study were recorded in the resting state, and the situation may be different for measurements during cognitive tasks. Moreover, reliability may depend on the selection of data (electrodes and time periods) for analysis and the specific data processing algorithm employed. Indeed, the methodology of PAC calculation from EEG signals in the literature is highly heterogeneous. The general procedure comprises several steps, including preprocessing, data selection, filtering the signals in the frequency ranges of interest, calculating a measure of PAC strength, and its normalization. In each of the steps, several options have been suggested and used in different studies. A number of methodological papers have compared some of the existing PAC calculation algorithms [22; 23; 24; 25; 26; 27; 28; 29]. They mainly compared different modulation indices and, less often, the available approaches to filtering (wavelet transform with different numbers of cycles, band-pass filters, and the selection of their bandwidths).