We analyzed eyes-closed resting-state EEG data from two independent, publicly available datasets (Table 1). The PhysioNet EEGBCI dataset (Schalk et al., 2004) contains 109 subjects with eyes-closed baseline runs; no subjects were excluded at the subject level (artifact handling was performed at the epoch/channel level); we included all subjects with eyes-closed resting recordings meeting artifact criteria (±100 μV epoch rejection, bad-channel interpolation) (N = 109). The LEMON Mind-Brain-Body dataset (Babayan et al., 2019) provides extensively validated resting-state EEG from 227 participants aged 20-77 years (mean age 38.9 ± 19.2 years, 107 female); after excluding 16 subjects with poor signal quality, the final sample was N = 211. Total sample size was N = 320 (pooled; Fs/channels vary by dataset). Both datasets were processed using identical pipelines and resampled to common frequency resolution (1 Hz bins) before analysis. Both datasets employed eyes-closed resting-state paradigms with durations ranging from 1-5 min. Despite heterogeneity in sampling rates (160-2500 Hz) and channel counts (62-64), all datasets were preprocessed using identical pipelines and resampled to common frequency resolution for spectral analysis. Cross-dataset replication serves as a stringent test of generalizability across acquisition parameters.

EEG data were preprocessed using MNE-Python (Gramfort et al., 2013). Signals were bandpass filtered (1-45 Hz, FIR filter), and bad channels were identified using variance thresholds (>3 SD from mean channel variance) and interpolated using spherical splines. For multi-channel recordings, signals were re-referenced to the average reference. Data were segmented into 4-second epochs, and epochs with amplitude exceeding ±100 μV were rejected.

Power spectral density (PSD) was computed using Welch's method (4-s Hann windows, 50% overlap). Spectral centroids for theta (4-8 Hz) and alpha (8-13 Hz) bands were calculated as: f_centroid=Σ(f×P(f))/Σ(P(f)) where, f represents frequency and P(f) represents power at that frequency, summed within the band of interest. We chose spectral centroids over peak-based frequency estimates (e.g., Individual Alpha Frequency) for several methodological reasons. While centroids are sensitive to spectral tails and band overlap, this sensitivity is mitigated by: (1) strict band definitions (theta 4-8 Hz, alpha 8-13 Hz), (2) FOOOF aperiodic correction, and (3) our frontal theta validation demonstrating effects are not driven by posterior alpha leakage. First, centroids provide stable estimates even when oscillatory peaks are weak or absent, which is particularly relevant for the theta band where peaks are often less sharply defined than alpha (Corcoran et al., 2018). Second, centroids capture the full power distribution within a band rather than a single maximum, making them less susceptible to noise-induced spurious peaks (Haegens et al., 2014). Third, recent work has demonstrated that spectral centroids exhibit superior test-retest reliability compared to peak frequencies for characterizing transient spectral events (Vidaurre et al., 2016). We acknowledge that centroids may be influenced by broadband spectral shape; however, our aperiodic correction analysis (Section 2.7) addresses this concern. For posterior analysis, PSDs were computed per channel and averaged across posterior channels using standard 10-20 system positions (O1, O2, Oz, P3, P4, Pz, P7, P8) where alpha is typically maximal. These specific channels were selected because they are present in both datasets despite differences in total channel count (62-64 channels), ensuring methodological consistency across the pooled sample. For frontal theta analysis, frontal channels (Fz, F3, F4) were used.

For each subject, we computed the ratio of alpha to theta spectral centroids (R = f_alpha/f_theta). The Phi Coupling Index was defined as: PCI=log((|R-2.0|+ε)/(|R-ϕ|+ ε)) where, ϕ = 1.618034 and ε = 0.01 for regularization (chosen to balance numerical stability without distorting natural frequency ratio distributions; see sensitivity analysis). Positive PCI indicates the ratio is closer to ϕ (ϕ-organized), while negative PCI indicates proximity to the harmonic 2:1 ratio. We note that PCI quantifies spectral proximity to ϕ, not phase-amplitude coupling in the classical sense. This terminological distinction is important: classical cross-frequency coupling measures dynamical interactions between oscillations, whereas PCI measures a structural property of the frequency spectrum—the degree to which bands are organized according to ϕ (Aru et al., 2015).

Theta-alpha convergence was quantified using a bounded metric to avoid singularities: Convergence = 1/(|f_alpha-f_theta| + 0.5). This bounded formulation avoids singularities and maps the physiological range of frequency differences (minimum ~0 Hz at 8 Hz boundary, maximum ~9 Hz) to convergence values between 0.1 and 2.0 (maximum convergence = 2.0 occurs when fα ≈ fθ near 8 Hz). Higher values indicate frequencies that approach the ~8 Hz boundary. For sensitivity analysis, we also computed the unbounded form (1/|Δf|) and an alternative 8 Hz symmetry metric.

Because both PCI and convergence are derived from the same frequency variables, their correlation could potentially arise from mathematical coupling. We employed six validation procedures:

Null model simulation: We generated 100,000 synthetic datasets by sampling f_theta and f_alpha from their empirical marginal distributions (preserving univariate distributions while destroying joint dependence). The observed correlation was compared to this null distribution.

Per-dataset replication: PCI-convergence correlations were computed separately for each of the two datasets to assess consistency across independent samples with different acquisition parameters.

Robust statistics: Spearman rank correlation was computed alongside Pearson to assess robustness to outliers and distributional assumptions.

ϕ-specificity parameter sweep: We computed generalized PCI using reference constants from 1.3 to 2.2 (step = 0.05) and correlated each with convergence. If ϕ-organization is meaningful, the PCI-convergence association should remain robust and significant at ϕ = 1.618.

Epsilon sensitivity: PCI was recomputed across ε = (0.001, 0.01, 0.1, 0.5, 1.0) to assess stability to regularization parameter choice.

Frontal theta validation: To address concerns that posterior theta may reflect alpha leakage rather than true theta oscillations, we computed PCI using frontal theta (Fz, F3, F4) with posterior alpha in the LEMON dataset (N=211), where both frontal and posterior channels were available.

To verify that PCI reflects oscillatory rather than aperiodic (1/f) organization, we performed sensitivity analysis using the FOOOF algorithm (Donoghue et al., 2020) with parameters: frequency range 1-40 Hz, peak width limits 1-12 Hz, maximum peaks unlimited, minimum peak height 0.1, aperiodic mode ‘fixed'. The aperiodic component was parameterized and subtracted in log-power space; negative residuals were set to zero before recomputing spectral centroids and PCI.

Pearson and Spearman correlations assessed relationships between PCI and convergence. Bootstrap resampling (10,000 iterations) provided 95% confidence intervals. For null model comparison, z-scores were computed as (r_observed - mean(r_null))/SD(r_null). All analyses used two-tailed tests with α = 0.05.

Across 320 subjects, the mean theta-alpha frequency ratio was 1.677 (SD = 0.142), deviating only 3.6% from ϕ = 1.618. This is substantially closer to ϕ than to the harmonic 2:1 ratio (16% deviation). A large majority of subjects (80.0%, N = 320) showed ϕ-organization (PCI > 0), indicating their theta-alpha ratios were closer to the golden ratio than to the harmonic 2:1 ratio (Figure 1).

PCI showed a strong positive association with theta-alpha convergence [r = 0.54, p < 10⁻²⁵, 95% CI (0.46, 0.62); Figure 2]. This relationship was highly robust to outliers (Spearman ρ = 0.82 (higher rank correlation reflects monotonic association with bounded transform), p = 1.50 × 10⁻⁶⁸), with the rank correlation exceeding the linear correlation.

Null model comparison: Two complementary null models assessed whether the observed correlation could arise from mathematical coupling alone. Null A (Marginal-Resampling, 100,000 draws): Sampling theta and alpha from empirical marginal distributions while destroying joint dependence produced mean r = 0.00 (SD = 0.10). The observed r = 0.54 exceeded this null by z = 5.4 (0 of 100,000 samples exceeded observed). Null B (Band-Constrained Monte Carlo, 10,000 draws per dataset): Random frequencies within physiological bands (θ: 4-8 Hz, α: 8-13 Hz) matched to empirical means/SDs produced negative null mean correlations (PhysioNet: mean r = −0.175, SD = 0.093; LEMON: mean r = −0.293, SD = 0.061). Observed correlations exceeded Null B by z = 8.6 (PhysioNet) and z = 12.9 (LEMON), with 0 of 10,000 exceeding observed. This validation approach—comparing observed statistics to null distributions derived from the data—is the established method for assessing non-trivial connectivity in neurophysiological research (Toppi et al., 2012).

Per-dataset replication: The PCI-convergence association replicated across both datasets (Table 2), with correlations ranging from r = 0.497 to r = 0.628, all highly significant (p < 10⁻⁹).

ϕ-specificity parameter sweep: The correlation between generalized PCI and convergence peaked near the golden ratio. At ϕ = 1.618, the correlation was within 0.03 of the optimum, and correlations declined sharply for reference constants further from ϕ (Figure 3), indicating that ϕ is near-optimal rather than arbitrary.

Epsilon sensitivity: PCI-convergence correlations remained significant across all epsilon values tested: ε = 0.001 (r = 0.41), ε = 0.01 (r = 0.54), ε = 0.1 (r = 0.62), ε = 0.5 (r = 0.74), ε = 1.0 (r = 0.78). The direction and significance of the association were stable across three orders of magnitude. Correlations remained positive and significant across ε, while effect size increased with larger ε. Sensitivity analysis using unbounded convergence (1/|Δf|) in the combined sample yielded r = 0.62 in the combined sample, qualitatively consistent with the primary bounded result (same direction and significance) (r = 0.54) (p < 10¹⁴), and an alternative 8 Hz symmetry metric (defined as negative distance: –(|8 – fθ| + |fα – 8|), so higher values indicate greater convergence) yielded highly similar results (r = 0.65), supporting robustness to convergence metric definition (Figure 4).

To address concerns that posterior theta may reflect alpha leakage rather than true theta oscillations, we computed PCI using frontal theta (Fz, F3, F4) with posterior alpha in the LEMON dataset (N = 211). Frontal and posterior theta centroids showed moderate correlation (r = 0.514), indicating partially overlapping but distinct signals.

Critically, frontal theta PCI-convergence correlation was substantially stronger (r = 0.718, p ≈ 10⁻³⁵) than posterior theta (r = 0.497). This is opposite to what volume conduction would predict—if posterior theta were merely alpha leakage, frontal theta should show weaker or no effect. Instead, using the anatomically appropriate source for theta yields even stronger correlations (Figure 5), providing convergent validation that the ϕ-organization reflects genuine theta-alpha relationships rather than measurement artifacts.

After removing the aperiodic component using FOOOF, the PCI-convergence association remained virtually unchanged (r = 0.574, compared to r = 0.54 uncorrected). The association remained similar after aperiodic correction, suggesting it reflects oscillatory organization rather than aperiodic spectral slope.

The frontal theta analysis provides crucial validation, consistent with recent source localization studies demonstrating that posterior theta originates from distinct neural generators. Güth et al. (2025) used MEG and simultaneous EEG-fMRI to localize right posterior theta to the parahippocampal gyrus during spatial navigation, supporting its cortical origin rather than volume conduction from frontal sources. Similarly, van Driel et al. (2012) employed current source density analysis to isolate local posterior theta activity and distinguish it from frontal theta. A potential concern with posterior channel analysis is that theta activity measured posteriorly might reflect volume-conducted alpha rather than true theta oscillations. If this were the case, computing PCI from frontal theta—where theta is known to be anatomically maximal—should produce weaker or absent effects. Instead, we observed a substantially stronger correlation (r = 0.718 vs. r = 0.497). This dissociation provides strong evidence that the ϕ-organization is consistent with genuine cross-frequency relationships rather than measurement artifacts arising from spatial mixing of signals.

Our findings extend theoretical proposals by Pletzer et al. (2010) that ϕ-based frequency ratios provide optimal neural decoupling. The observation that 80% of subjects cluster near ϕ supports their hypothesis that the golden ratio represents “the highest physiologically possible desynchronized state.” Furthermore, our results align with Rodriguez-Larios and Alaerts (2019) and Rodriguez-Larios et al. (2020), who found that resting states favor non-harmonic ~1.6:1 coupling ratios while active cognition shifts toward harmonic 2:1 arrangements.

The consistent replication across two datasets with markedly different acquisition parameters (sampling rates: 160-2500 Hz; channel counts: 62-64; different laboratories and recording protocols) strengthens confidence that the finding is generalizable. The range of correlations (r = 0.50-0.63) across datasets is consistent with a robust effect subject to typical sampling variability.

Several methodological aspects merit discussion. First, PCI and convergence share mathematical dependency on the same frequency variables (f_theta, f_alpha), which could inflate correlations. However, our null model testing demonstrated that the observed correlation substantially exceeds what would be expected from this mathematical coupling alone (>5 SD). Additionally, the cross-dataset replication and frontal theta validation provide converging evidence that the effect is not merely a statistical artifact.

Second, the epsilon sensitivity analysis demonstrates that our results are robust to regularization parameter choice, with correlations remaining significant across three orders of magnitude (ε = 0.001 to 1.0). This addresses reviewer concerns about parameter arbitrariness.

Third, with our bounded convergence metric (1/(|Δf|+0.5)), values typically range from ~0.15 (for Δf≈6 Hz) to ~0.67 (for Δf≈1 Hz), reflecting the physiological range of theta-alpha frequency differences. Fourth, while alternative convergence formulations (e.g., the 8 Hz symmetry metric) yielded slightly different correlation magnitudes (r = 0.63-0.65), all showed the same pattern of results, indicating that the singularity in the original metric does not drive our conclusions.

Several limitations warrant consideration. First, while we demonstrate association, the functional significance of ϕ-organization requires investigation with behavioral paradigms. Does higher ϕ-organization predict cognitive flexibility, meditation experience, or other meaningful outcomes? Second, this analysis focused on resting-state data; whether ϕ-organization shifts systematically during cognitive tasks—as suggested by Rodriguez-Larios et al.—remains to be directly tested with within-subject designs. Third, the spectral centroid approach, while robust to noise, may be affected by individual differences in spectral shape; future work could examine whether peak-based frequency estimates yield comparable results. Fourth, we did not assess test-retest reliability; comparing centroid-based results to peak-based frequency estimates (e.g., Individual Alpha Frequency) would further validate our approach; establishing the temporal stability of individual differences in ϕ-organization would strengthen claims about its trait-like nature. Additionally, testing whether similar ϕ-organization patterns emerge in other frequency band pairs (e.g., alpha-beta during attentional tasks) would clarify the generality of the “golden rhythms” framework. Fifth, while all recordings were nominally eyes-closed resting state, we did not control for vigilance fluctuations, drowsiness, or micro-sleep episodes using objective markers (e.g., EOG, spectral slope dynamics), which are known to affect theta and alpha frequencies. Without objective vigilance markers (e.g., EOG monitoring, spectral slope dynamics), we cannot definitively rule out that some variance in ϕ-organization reflects transient state differences rather than stable trait-like properties. Future studies should incorporate vigilance monitoring to address this limitation.

Golden ratio organization in human EEG theta-alpha frequencies is a robust association, replicating across 320 subjects from two independent datasets with diverse acquisition parameters. The 80% prevalence of ϕ-organized spectral architecture, the consistency across datasets (r = 0.50-0.63), and the striking frontal theta validation (r = 0.718) suggest this may reflect a consistent pattern of theta-alpha organization in resting-state EEG. These results provide converging evidence against a simple spatial-mixing explanation and are consistent with ϕ-related organization of theta-alpha frequencies. The findings have potential implications for understanding neural computation, developing biomarkers of brain state, and potentially explaining why ϕ appears so frequently across biological systems.