Coupling of pupil- and neuronal population dynamics reveals diverse influences of arousal on cortical processing

Fluctuations in arousal, controlled by subcortical neuromodulatory systems, continuously shape cortical state, with profound consequences for information processing. Yet, how arousal signals influence cortical population activity in detail has so far only been characterized for a few selected brain regions. Traditional accounts conceptualize arousal as a homogeneous modulator of neural population activity across the cerebral cortex. Recent insights, however, point to a higher specificity of arousal effects on different components of neural activity and across cortical regions. Here, we provide a comprehensive account of the relationships between fluctuations in arousal and neuronal population activity across the human brain. Exploiting the established link between pupil size and central arousal systems, we performed concurrent magnetoencephalographic (MEG) and pupillographic recordings in a large number of participants, pooled across three laboratories. We found a cascade of effects relative to the peak timing of spontaneous pupil dilations: Decreases in low-frequency (2–8 Hz) activity in temporal and lateral frontal cortex, followed by increased high-frequency (>64 Hz) activity in mid-frontal regions, followed by monotonic and inverted U relationships with intermediate frequency-range activity (8–32 Hz) in occipito-parietal regions. Pupil-linked arousal also coincided with widespread changes in the structure of the aperiodic component of cortical population activity, indicative of changes in the excitation-inhibition balance in underlying microcircuits. Our results provide a novel basis for studying the arousal modulation of cognitive computations in cortical circuits.

Fluctuations in arousal, controlled by subcortical neuromodulatory systems, continuously shape cortical state, with profound consequences for information processing. Yet, how arousal signals influence cortical population activity in detail has only been characterized for a few selected brain regions so far. Traditional accounts conceptualize arousal as a 20 homogeneous modulator of neural population activity across the cerebral cortex. Recent insights, however, point to a higher specificity of arousal effects on different components of neural activity and across cortical regions. Here, we provide a comprehensive account of the relationships between fluctuations in arousal and neuronal population activity across the human brain. Exploiting the established link between pupil size and central arousal systems, 25 we performed concurrent magnetoencephalographic (MEG) and pupillographic recordings in a large number of participants, pooled across three laboratories. We found a cascade of effects relative to the peak timing of spontaneous pupil dilations: Decreases in low-frequency (2-8 Hz) activity in temporal and lateral frontal cortex, followed by increased high-frequency (>64 Hz) activity in mid-frontal regions, followed by linear and non-linear relationships with 30 intermediate frequency-range activity  in occipito-parietal regions. The non-linearity resembled an inverted U-shape whereby intermediate pupil sizes coincided with maximum 8-32 Hz activity. Pupil-linked arousal also coincided with widespread changes in the structure of the aperiodic component of cortical population activity, indicative of changes in the excitation-inhibition balance in underlying microcircuits. Our results provide a novel basis for 35 studying the arousal modulation of cognitive computations in cortical circuits.

Introduction
Variations in cortical state profoundly shape information processing and, thus, cognition (Busse et al., 2017;Fu et al., 2014;Zagha et al., 2013). These variations occur continuously due to subtle fluctuations in the level of arousal, and even in the absence of changes in overt behavior (Harris & Thiele, 2011;McGinley, Vinck, et al., 2015). Two key mediators of arousal-45 dependent variations are the brainstem nucleus locus coeruleus (LC), which supplies noradrenaline (NE), and the basal forebrain (BF), which supplies acetylcholine (ACh) (Harris & Thiele, 2011;Hasselmo, 1995;Lee & Dan, 2012;Steriade, 1996). It has long been thought that these nuclei are organized homogeneously and innervate cortical target regions diffusely, and indiscriminately across cortical regions. This has led to the idea that the arousal system 50 acts as a global "broadcast signal" and uniform controller of cortical state (Aston-Jones & Cohen, 2005;Harris & Thiele, 2011;Leopold et al., 2003;Turchi et al., 2018).
More recently, however, this view has been challenged: Several neuronal subpopulations with distinct projection targets have been found in both the LC and BF (Chandler et al., 2013(Chandler et al., , 2019Sarter et al., 2009;Schwarz & Luo, 2015;Totah et al., 2018;Zaborszky et al., 2015;55 Záborszky et al., 2018). In addition, neuromodulator receptors exhibit a rich and diverse distribution across the cortex (Burt et al., 2018;van den Brink et al., 2019). In this report, we comprehensively show how this structural and molecular heterogeneity translates into a multitude of arousal effects on neuronal population activity across the human cortex.
Rodent research has mostly been limited to circumscribed sensory cortical areas (Joshi & Gold, 2020;McGinley, David, et al., 2015;Reimer et al., 2014; but see Shimaoka et al., 2018). 75 Insights from non-invasive magnetic resonance imaging in humans, on the other hand, have revealed spatially rich correlations of cortical activity with pupil-linked arousal (Groot et al., 2021;Schneider et al., 2016;Yellin et al., 2015). However, the low temporal resolution of these measurements complicates direct comparisons with results from electrophysiological recordings in rodents and non-human primates. Therefore, we still lack a comprehensive 80 understanding of which features of cortical population activity are shaped by arousal and how such effects are distributed across cortex.
To close these gaps, we assessed the co-variation of spontaneous fluctuations in pupillinked arousal and whole-head magnetoencephalography (MEG) recordings in resting human participants across brain regions (Figure 1), components of cortical activity (oscillatory and 85 aperiodic), and temporal lags. Our findings provide a comprehensive picture of linear and nonlinear relationships between intrinsic fluctuations in pupil-linked arousal and neural population activity in different brain regions. First, we demonstrate opposite effects of pupillinked arousal on low-and high-frequency components of cortical population activity across large parts of the brain, similar to previously observed effects in rodent sensory cortices. 90 These distinct effects occurred in a systematic temporal sequence. Second, we also uncovered hitherto unknown, non-linear (inverted U-shape) arousal effects in occipitoparietal cortex. Third, we identified spatially widespread correlations between pupil-linked arousal and the structure of aperiodic activity, suggesting underlying changes in cortical excitation-inhibition balance (Gao et al., 2017;Waschke et al., 2021). 95

105
We analysed data from concurrent whole-head MEG and pupil recordings in 81 human participants, collected at three MEG laboratories (Glasgow, UK; Hamburg & Münster, Germany). All recordings were carried out under constant illumination and in soundattenuated environments, thus removing fluctuations in pupil size due to changes in sensory input. Participants were instructed to keep their eyes open, while fixating and resting 110 otherwise. The length of individual recordings varied between 5 -10 min (see Methods for further details).
For each of the three MEG laboratories, we obtained time courses of MEG activity and of pupil diameter fluctuations with highly similar spectral characteristics (Supplementary Figures  S1 and S2). In what follows, we refer to the time course of pupil diameter as 'pupil'. Previous 115 work in rodents suggests that this time course and its temporal derivative indicate distinct (cholinergic and noradrenergic, respectively) neuromodulatory inputs to cortex (Reimer et al., 2016). Thus, we also tested for relationships of MEG activity with the first temporal derivative Robust pupil-brain couplings across the spectrum of cortical population activity We evaluated the strength of pupil-brain coupling across a wide range of frequency bands in MEG sensor space by means of mutual information (MI) which measures all (linear and 125 nonlinear) associations between two signals (Ince et al., 2017). We calculated the MI between pupil traces and power time courses of band-limited cortical activity across 25 log-spaced center frequencies from 2 to 128  For all three MEG laboratories, standardised-MI spectra showed similar patterns of widespread pupil-brain associations across frequencies with prominent peaks in the 8-32 Hz range (Figure 2A). MI scalp distributions within this frequency range further underpinned the 135 commonalities between recording sites despite different MEG systems and setups ( Figure  2B). Significant relationships with the pupil were also evident for the low-frequency (4)(5)(6)(7)(8) and high-frequency (64-128 Hz) components of MEG-power ( Figure 2A).

Fluctuations in pupil diameter lag behind changes in cortical population activity in a frequency-dependent manner
Having established robust pupil-brain couplings, we quantified the temporal relationship between both signals. Changes of neural activity in brainstem arousal centers such as the LC generate changes in neural activity that exhibit considerable temporal lag (~500 ms in 155 A Pupil-size covariation with cortical activity across canonical frequency bands -by laboratory macaques) and temporal low-pass characteristics (Breton-Provencher & Sur, 2019;de Gee et al., 2017;Hoeks & Levelt, 1993;Joshi et al., 2016;Korn & Bach, 2016). This temporal smoothing is produced by the peripheral pathways (nerves and muscles) controlling pupil motility (Hoeks & Levelt, 1993;Korn & Bach, 2016). Thus, one might expect the effect of neuromodulatory arousal signals to occur earlier in cortical population activity (a single 160 synapse from LC to cortex, unless mediated by a third structure, such as the thalamus) than in pupil diameter.
We computed cross-correlations between fluctuations in frequency-resolved sensor-level power and pupil (Methods). We found the maximum correlation magnitude between pupil diameter and power fluctuations at a frequency-dependent lag ( Figure 3A, B): For activity in 165 the low (around 2-4 Hz) and high frequency bands (around 64-128 Hz), correlation magnitudes peaked at lags of 832 ms (negative peak) and 562 ms (positive peak), respectively. This result was highly consistent across the three MEG sites ( Figure 3A). These latencies indicated power fluctuations that preceded corresponding variations in pupil diameter. Correlations with power fluctuations in the 8 -16 Hz range, on the other hand, 170 peaked at a lag closer to zero (210 ms; negative peak), indicating a temporally closer relationship ( Figure 3A, B; see Discussion).
Fluctuations in pupil, and pupil-derivative, may primarily reflect cholinergic and noradrenergic activity, respectively (Reimer et al., 2016). Consistent for all three data sets, the cross-correlation analysis for pupil-derivative revealed peaks closer to zero lag for low 175 and high frequencies (2 -4 Hz negative peak at 277 ms, 64 -128 Hz positive peak at -90 ms; Supplementary Figure S3). In contrast, correlations with activity in the 8 -16 Hz range peaked somewhat later (700 ms; negative peak; Supplementary Figure S3). The shorter lags for the pupil derivative are consistent with previous findings in rodents (Reimer et al., 2016).
The cortical innervation profiles of subcortical arousal centers and the cortical receptor 180 composition for the neuromodulators released by these centers exhibit substantial heterogeneity across cortex (Ramos & Arnsten, 2007;van den Brink et al., 2018van den Brink et al., , 2019. We, therefore, reasoned that the nature (frequency-dependence and/or temporal lag) of pupilbrain couplings may differ considerably across cortical areas. We next tested the topographical profile of the cross-correlations. To account for the different channel layouts 185 of the three MEG systems, we averaged cross-correlations, computed for individual MEG channels, within 39 non-overlapping channel groups, consisting of channels that were grouped based on their location along the anterior-to-posterior axis (see Methods). In keeping with the previous analyses, we collapsed the data into three non-overlapping frequency bands of interest: 2 - 4 Hz,[8][9][10][11][12][13][14][15][16]  Consistent with the results obtained from sensor averages, a negative correlation between pupil diameter fluctuations and power in the low frequencies (2-4 Hz) peaked at a lag of 210 around 900 ms, with a spatial profile that was relatively uniform across brain areas (Figure 3C,left). For the high frequencies (64-128 Hz), on the other hand, the positive peak was confined to posterior regions, with a lag similar to the low frequency band (Figure 3C,right).
In the 8 to 16 Hz-range, the picture was more nuanced (Figure 3C,center): The close-to zerolag positive correlation dominated posterior sensors whereas anterior sensors registered a 215 negative correlation with a peak around 1000 ms. This suggests that the correlation between pupil diameter and power in the 8 -16 Hz range may be driven by independent mechanisms that operate on distinct time scales (see Discussion).
The spatio-temporal correlation patterns between the pupil-derivative and band-limited power is shown in Supplementary Figures S3 and S4. 220 In summary, we identified lags between spontaneous fluctuations in pupil diameter and low-and high-frequency activity that are consistent with previously reported lags between noradrenergic and cholinergic activity, respectively, and spontaneous pupil diameter fluctuations in rodents. In addition, we found that pupil-brain coupling in the intermediate 8- 16 Hz range peaked with a time lag closer to maximum pupil dilation for posterior sensors.

225
Our findings thus suggest that the correlations in different frequency bands, and for different cortices within the intermediate frequency range, may be mediated by different neuromodulatory systems (see Discussion).

Spatial and spectral dissociations of pupil-power correlations
Having established the distinct temporal coupling profiles for different frequency bands, 230 we now tested which cortical areas received the strongest modulation. Sensor space analyses suggested that effects of pupil-brain coupling were not uniformly distributed along the anterior-posterior axis (Figure 3). To map out cortical loci of frequency-specific pupilbrain coupling in detail we used spectrally resolved source-reconstructed MEG data. Given the high consistency of the sensor-space results, we report data pooled across the three 235 recording sites for the following analyses.
Based on the results obtained from the cross-correlation analysis and in line with previous reports (Hoeks & Levelt, 1993;Joshi et al., 2016;Reimer et al., 2016), we shifted the pupil signal 930 ms forward (with respect to the MEG signal). In other words, we shifted the pupil signal to compensate for the lag that had been observed between external manipulations of 240 arousal (Hoeks & Levelt, 1993) as well as spontaneous noradrenergic activity (Reimer et al., 2016) and changes in pupil diameter. Using the forward-shifted pupil signal, we computed correlations with fluctuations in spectral power (Figure 4), separately for 8799 source locations covering the 246 regions of the Brainnetome atlas (Fan et al., 2016;see Methods). For each region, pupil-power associations were quantified by means of Spearman's rank 245 correlation. All steps were repeated for correlations between pupil-derivative and power fluctuations, however without accounting for the lag as previous research showed nearinstantaneous correlations between neuromodulatory activity and fluctuations in pupil derivative (Reimer et al., 2016).
Averaged across all source locations, correlations between intrinsic pupil and power 250 fluctuations were frequency-dependent ( Figure 4A): Consistent with findings from invasive recordings in primary visual cortex of rodents (Reimer et al., 2014;Vinck et al., 2015), we found a robust negative correlation between pupil diameter and power in the low frequencies (2-8 Hz) and a numerically small, albeit statistically robust, positive correlation in the high frequencies (between 50 Hz and 100 Hz). In order to investigate the spatial distribution of the 255 correlations, we sorted the regions of the Brainnetome atlas along the anterior-posterior axis. Whereas the correlations in the low frequencies covered large parts of the brain (except for a few posterior regions, Figure 4B), pupil-brain coupling in other frequency ranges exhibited a spatially richer structure, with positive (8 -32 Hz) and negative (> 50 Hz) correlations in more occipital regions and correlations of opposite sign in more anterior sources ( Figure 4B). 260 Previous studies have reported largely similar effects of pupil-linked arousal on population activity in the different sensory cortices (McGinley, David, et al., 2015;Reimer et al., 2014;Vinck et al., 2015; but see Shimaoka et al., 2018). Capitalising on our source reconstructions, we tested if the correlations between pupil diameter and cortical activity varied across cortical regions of interest (see inset maps in Figure 4C), including three sensory areas (visual cortex, 265 auditory cortex and somatosensory cortex) as well as anterior cingulate cortex, a region receiving strong innervations from locus coeruleus (Chandler et al., 2013). Patterns of pupilpower correlations were highly region-specific: Correlations in the primary visual cortex were positive within the 8-32 Hz range (Figure 4C,top). Correlations in the auditory cortex, on the other hand, were dominated by negative correlations in the low frequency range (< 8 Hz; 270 Figure 4C, second to top). Somatosensory cortex exhibited negative correlations in the low and intermediate frequencies (2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16) and, at the same time, positive correlations in the high frequency range (~ 64 Hz; Figure 4C, second to bottom). Correlations with activity in anterior cingulate were significant across most frequencies, with a spectral pattern similar to auditory cortex ( Figure 4C, bottom). 275

285
between pupil and power (in clockwise direction, starting from the top left) at 4 Hz,11.3 Hz,19.0 Hz, and 64 Hz. All maps are masked at P = 0.05 (two-tailed paired t-test; FDR-adjusted, q=0.1).
For further insights into the spatial distribution of the observed correlations, we mapped correlations of pupil and pupil-derivative time courses with band-limited power fluctuations 290 across all 8799 source locations, again estimated at a lag of 930 ms for pupil and a lag of 0 ms for pupil-derivative. To this end, we focused on four frequency bands of interest, with the following center frequencies ( Figure 4D): 4 Hz,11.3 Hz,19 hz, and 64 Hz, corresponding to the classical theta, alpha, beta, and gamma frequency ranges. In the low frequencies (centered at 4 Hz), correlations were negative in the majority of source locations, with bilateral 295 peaks in the anterior sections of the hippocampus ( Figure 4D, top left). Significant positive correlations between pupil and activity in the alpha-/beta range (11.3 Hz and 19 Hz) were located in the visual cortices of both hemispheres ( Figure 4D, top right and bottom left). Significant negative correlations were found in left and right motor and somatosensory cortices ( Figure 4D, top right), although limited to activity in the alpha-range (centered at 11.3 300 Hz). Similarly, correlations with high frequency activity (centered at 64 Hz) were non-uniform across space, with negative correlations located in visual cortices bilaterally, and positive correlations across large parts of the prefrontal cortex ( Figure 4D, bottom right).
Repeating the analyses for the pupil derivative, we found largely similar spatial patterns for the correlations, albeit with reduced magnitudes (Supplementary Figure  Briefly summarised, source localisation of pupil-brain coupling revealed frequency-band specific patterns of arousal influences on cortical processing. While coupling with low-310 frequency activity was relatively uniform, coupling with intermediate-and high-frequency activity showed effects of opposite signs for different cortices, even within frequency bands, indicating that arousal has discernible effects on different sensory cortices, as well as higherorder associative regions of the human brain at rest.

Nonlinear relations between pupil-linked arousal and band-limited cortical activity 315
Above described analyses focused on quantifying linear relationships between intrinsic fluctuations in pupil diameter, as well as pupil-derivative, and band-limited cortical population activity. However, arousal may exert more complex influences on cortical activity that roughly follow a nonlinear, inverted-U shape, known as the Yerkes-Dodson law (Yerkes & Dodson, 1908;Aston-Jones & Cohen, 2005;McGinley, David, et al., 2015). Evidence suggests that 320 gain of neuronal input-output functions as well as behavioural performance peak at intermediate arousal levels (He & Zempel, 2013;McGinley, Vinck, et al., 2015;Waschke et al., 2019), flanked by lower gain and performance at lower (e.g. disengagement and drowsiness) and higher arousal states (e.g. stress). We, therefore, reasoned that local pupilpower couplings may exhibit similar nonlinear relationships. To test this, we divided the data 325 into non-overlapping temporal segments of 2s length and computed power spectra as well as mean pupil size (or mean pupil-derivative) for each segment. Next, we averaged the power spectra within each of 14 evenly sized pupil bins, ranging from smallest to largest pupil diameter on each block.
Focusing on three predefined frequency bands, the relation between pupil diameter and 330 pupil-derivative with low-frequency (2-4 Hz) and high-frequency activity (64-128 Hz) was wellcharacterized by a linear relationship ( Figure 5A, left and right) when averaged across all source locations. In contrast, activity in the 8-16 Hz range exhibited a high degree of nonlinearity that followed an inverted-U-relationship ( Figure 5A, center).
To test for potential nonlinear relationships exhaustively across frequencies and space, we 335 fitted a quadratic polynomial to the power-by-pupil bin functions, separately for each source region and frequency (see Methods). Next, we mapped the coefficient of the quadratic term of the polynomial (β2) across space and frequencies. In order to account for differences in absolute power (offset), we normalized (i.e., z-scored) the power-by-pupil functions prior to model fitting. The analysis revealed that inverted-U relationships were exclusive to the range 340 of 8-32 Hz (i.e., β2 significantly smaller than zero; two-sided paired t-test; Figure 5B). Focusing on the quadratic nonlinearity in four regions of interest (as used in the previous 355 analyses) pointed towards locally emphasised intermediate-frequency inverted-U shaped relationships (auditory and somatosensory cortices, Figure 5C). A spatially more fine-grained analysis of these effects revealed that the inverted-U relationship in the alpha band (~11 Hz) peaked in the left superior parietal lobule ( Figure 5D, center left). The inverted-U relationship between activity in the low beta band (~19 Hz) and pupil diameter, on the other hand, 360 additionally comprised bilateral parietal and temporal regions ( Figure 5D, right). We also found non-linear relationships when analysing the pupil derivative instead of the pupil diameter, exhibiting different spatial patterns (Supplementary Figure S7).
Taken together, we found evidence for a non-linear inverted-U relationship between pupillinked arousal and band-limited cortical activity that was confined to the classical alpha/beta 365 frequency band (8-32 Hz) and parieto-occipital cortices. Put differently, cortical activity in this frequency range, and in the respective cortical regions, was maximal at medium pupil diameters that indicate intermediate optimal levels of arousal. Activity in surrounding lowand high-frequency ranges showed exclusively linear effects.

Pupil-linked arousal also predicts changes in aperiodic components of cortical 370 population activity
So far, we have shown how pupil-linked arousal co-varies with band-limited cortical activity. However, cortical electrophysiological processes can also be characterised by a broadband, or aperiodic component. This component can be understood as the decay rate of power with increasing frequency, and is therefore also known as the 1/f component. The 375 slope of the 1/f component reflects fluctuations in attentional state (Waschke et al., 2021) and has been suggested to track the ratio between excitation and inhibition (in short: E/I) in the underlying neuronal circuits (Gao et al., 2017). As neuromodulators linked with the regulation of arousal have been shown to change cortical E/I (Froemke, 2015;Pfeffer et al., 2018Pfeffer et al., , 2020, we next investigated the relation between intrinsic fluctuations in pupil diameter (as well as 380 its temporal derivative) and the aperiodic component in the power spectrum. Our goal was to characterize correlations between intrinsic fluctuations in pupil-linked arousal and the slope of the aperiodic component, as a potential marker of cortical E/I.
To this end, we parameterized the power spectra of consecutive data epochs of the source-projected MEG data in a frequency range from 3 to 40 Hz, separately for all source 385 locations (Donoghue et al., 2020). The resulting parameter estimates described the periodic and the aperiodic components per epoch (data segmentation as in the previous section; see Methods).
We correlated fluctuations in the estimated spectral exponent of the aperiodic components, a measure that expresses the steepness of the 1/f slope, with fluctuations in 390 pupil diameter. The spectral exponent decreased linearly as a function of pupil diameter ( Figure 6A) -note that an decrease in the spectral exponent means a shallower decay of power as a function of frequency (see Methods for details). This resulted in a significant negative correlation when averaged across all source locations ( Figure 6B, grey dots), with little variability when zooming in on four regions of interest (Visual cortex: mean r = -0.0266; 395 P = 0.0013; Auditory cortex: mean r = -0.0484 , P = 8.12*10 -8 ; Somatosensory cortex: mean r = -0.0412, P = 3.75*10 -5 ; Anterior cingulate cortex: r = -0.0325, P = 1.13*10 -5 ; two-tailed paired t-tests against zero). In fact, the negative correlation between the spectral exponent and pupil diameter was significant in the vast majority of brain regions (8197/8799 of source locations, or 93.2%, after adjusting for a false discovery rate of 0.1; Figure 6C). A correlation 400 was notably absent in several extrastriate occipital cortices.
We found a similarly widespread negative correlation between the spectral exponent and pupil-derivative. Albeit weaker in magnitude, this effect additionally encompassed the entire occipital cortex and showed stronger local maxima in parietal cortices, when compared with the correlation with pupil (Supplementary Figure S8A-C). 405 In sum, these results speak to a global effect of pupil-linked arousal on the excitationinhibition ratio across the brain, with few potential exceptions (occipital cortices). Minor differences between pupil-and pupil-derivative-based spatial patterns of these effects may again point towards the involvement of different neuromodulators.

Couplings between pupil-linked arousal and periodic components of cortical population activity are not induced by couplings with aperiodic components
Given the pupil-linked changes in the slope of the (aperiodic component of) the cortical power spectrum, one concern is that these introduced spurious changes in band-limited 425 power (Donoghue et al., 2020), thus potentially producing spurious correlations between fluctuations in pupil diameter and band-limited power. To test whether the correlations between pupil and spectral power reported above were merely a result of changes in the slope of the power spectra, we subtracted the fitted aperiodic component from the empirical power spectra (in the range from 3 Hz to 40 Hz), separately for each time segment and source 430 location. Subsequently, we correlated the residual power estimates with pupil diameter, across all segments.
The pattern of correlations between pupil diameter and power remained largely unchanged ( Figure 6D): Fluctuations in pupil diameter were negatively correlated with activity in the low frequencies, around 4 Hz, as well as with activity in a higher frequency range around 32 Hz. 435 Band-limited activity in the range between 8 and 16 Hz exhibited robust positive correlations.
Correlations between residual power and pupil-derivative, on the other hand, were strongly reduced in magnitude and not significant for most frequencies (Supplementary Figure S8D).
These control analyses confirm that the couplings between pupil (and pupil-derivative) and band-limited cortical activity reported in the preceding sections were genuine, reflecting 440 modulations of specific cortical oscillations by arousal.

Discussion
Pupil-linked arousal exerts pivotal influences on the states of sensory cortices (McGinley, David, et al., 2015;Reimer et al., 2014;Vinck et al., 2015) and associative cortical regions (Joshi & Gold, 2020). However, subcortical arousal systems also project to a wealth of other 445 regions across the cerebral cortex. This warranted a detailed investigation of how these projections shape cortical population activity. To achieve this, we exhaustively mapped the relationship between intrinsic fluctuations in pupil-linked arousal and resting-state human MEG-recorded cortical activity. In a large sample (N=81), our mapping revealed diverse and consistent linear and non-linear arousal effects on band-limited, as well as aperiodic cortical 450 activity.

Arousal modulates cortical activity across space, time and frequencies
Our findings are consistent with previously reported links between both locomotionrelated and pupil-related changes in arousal and LFP recordings in sensory cortices of rodents (McGinley, David, et al., 2015;Reimer et al., 2014Reimer et al., , 2016. Briefly put, periods of 455 elevated arousal were accompanied by a decrease in low-and an increase in high-frequency activity, similar to cortical state changes typically observed during increased attention (Harris & Thiele, 2011;Schwalm & Jubal, 2017). The observed temporal relationships also align with the lags between manipulations of arousal state and pupil dilations of around 500-1000 msec in humans (Hoeks & Levelt, 1993;Wierda et al., 2012) and a ~500 msec delay between LC 460 activation and peak pupil dilation in non-human primates (Joshi et al., 2016). In rodents, Reimer et al. (2016) and Breton-Provencher & Sur (2019) reported peak correlations between noradrenergic axonal activity and pupil diameter at around 900 ms, whereas the correlation between basal forebrain-regulated cholinergic activity and pupil diameter peaked around 500 ms (Reimer et al., 2016). Here, we report comparable peak lags for low (832 ms) and high 465 frequency activity (562 ms), suggesting that different neuromodulatory systems may drive correlations in different frequency ranges in the human brain. Future research using brainwide maps of the distribution of neuromodulatory receptors (van den Brink et al., 2019) may help to further disentangle the differential contributions of noradrenaline and acetylcholine as well as their various receptor subtypes. 470 In addition to the bipartite division into low-and high frequency effects known from rodent sensory cortices, we also report linear pupil-brain coupling in a third intermediate alpha/betarange frequency band. Specifically, a negative correlation, localised to somatomotor cortices and specific to alpha-band (8-12 Hz) activity, peaked earlier than a positive correlation that localised to occipital visual cortices and encompassed alpha and beta bands (see Figures  475 3C & 4D). Of all effects found, the effect in visual cortices had the shortest lag with respect to pupil dilation (~200msec). This short-lag positive correlation may stem from interactions between the locus coeruleus and regions of the thalamus (McCormick, 1989;McCormick et al., 1991;Stitt et al., 2018), in particular those implicated in the generation of occipital alphaband activity via thalamo-cortical loops (e.g. the Pulvinar; Saalmann et al., 2012). Put 480 differently, a more indirect multi-synaptic route to neuromodulation of occipital alpha activity may delay the lag relative to other effects and put it closer to the time of maximum pupil dilation. Another, not mutually exclusive, possibility is that the arousal-dependent low-lag alpha-/beta-band activity is linked to the activity of a distinct subpopulation of GABAergic neurons in the locus coeruleus, that have been shown to exhibit a temporally closer 485 relationship with pupil diameter compared to the noradrenergic neurons (Breton-Provencher & Sur, 2019).
Note that reported changes in alpha activity in the visual cortex cannot be trivially explained by pupil dilation altering retinal illumination. In that case, one would expect a negative relationship, i.e. decreasing alpha power with increasing pupil size (e.g., Donner & 490 Siegel, 2011). Moreover, changes in alpha activity precede changes in pupil dilation, albeit with a shorter lag compared to activity in other frequency ranges (see Figure 3).
Correlations between cortical activity and pupil derivative largely mirrored the results from pupil diameter, although most links were found to be weaker in magnitude. This is likely due to the fact that both signals (pupil diameter as well as its temporal derivative) correlate to 495 some extent with noradrenergic and cholinergic activity (Reimer et al., 2016).

An arousal-triggered cascade of activity in the resting human brain
The temporal profile of the diverse relationships between pupil-linked arousal and local cortical activity described above points towards a cascade of interrelated effects that the resting brain may undergo routinely. Several scenarios may produce such a cascade: First, 500 episodes of arousal could arise seemingly at random, from underlying complex non-linear neural dynamics (Robinson et al., 2003;Kringelbach et al., 2020). Second, they could be coupled to other physiological processes, such as respiration (Kluger & Gross, 2020;Zelano et al., 2016). Similarly to pupil fluctuations (Bouma & Baghuis, 1971;Pomè et al., 2020;Turnbull et al., 2017), respiration has an average base rate of ~0.2 Hz (Del Negro et al., 2018;505 Fleming et al., 2011). Finally, as a third scenario, the observed cascade may imply a specific cognitive sequence: Most participants will likely have 'mind-wandered', or engaged in 'taskunrelated thoughts' (Groot et al., 2021) while resting. The sequence of effects we observe may therefore be a transient re-orientation towards the external sensory environment -a 'monitoring sweep'. This is consistent with the increase in pupil-linked arousal, and its inverse 510 relationship with an early dip in low-frequency activity, i.e. the cessation of a globally synchronised cortical state, and agrees with a concurrent widespread increase in highfrequency cortical activity that signifies increased cortical excitability to process external sensory input (Hanslmayr et al., 2011;Harris & Thiele, 2011).
One aspect of mind-wandering is the retrieval of information from memory that heavily 515 relies on theta rhythms in the hippocampus (Siegle & Wilson, 2014). Here, we localised the strongest pupil (and pupil-derivative) correlations with theta-range power fluctuations in bilateral hippocampal regions (see Figure 4D), suggesting an arousal modulation of memory processes. In support, hippocampal locus-coeruleus-noradrenaline modulation has long been established in the rodent brain (Sara, 2009;Segal & Bloom, 1976b, 1976a, and 520 McGinley et al. (2015) reported negative correlations between pupil-linked arousal and the occurrence of theta-range hippocampal ripple events. Moreover, a recent human neuroimaging study found that pupil-linked arousal modulated blood oxygenation in (para-) hippocampal regions during memory processes (Clewett et al., 2018).
Further on, the observed sequence of effects concludes with a peak in alpha/beta activity 525 that is triggered roughly at the time of maximum correlations of arousal with low-and highfrequency activity. The increase in synchronised alpha/beta activity marks the return to a strong regime where external sensory (visual) input is attenuated, potentially indicating a return to an introspective focus. A strong alpha rhythm thereby imposes widespread bouts of "pulsed inhibition" on visual cortex that suppress local gamma band activity (Jensen & 530 Mazaheri, 2010; Haegens et al., 2011; for a review see Van Diepen et al., 2019). This relationship between alpha and gamma activity has been well described in non-human primate- (Spaak et al., 2012) and human visual cortices, with a role for cholinergic neuromodulation (Bauer et al., 2012). Corroborating this scenario, we also observe a decrease in local gamma-band activity in visual cortices. 535 The three scenarios described here are not mutually exclusive and may explain one and the same phenomenon from different perspectives. A pivotal manipulation to test these assumptions will be to contrast this cascade with couplings of pupil-linked arousal and cortical activity between different behavioural states.

Nonlinear relationship between periodic cortical activity and pupil-linked arousal 540
Nonlinear relations between arousal and other factors are well-characterized, most prominently in the form of the Yerkes-Dodson law, which posits an inverted-U-shaped relation between arousal and cognitive performance (Yerkes & Dodson, 1908). More recently, an inverted-U relation has also been identified for pupil-linked arousal fluctuations and detection performance and low-frequency (2-10 Hz) membrane potential fluctuations in 545 primary auditory cortex (McGinley, David, et al., 2015). Here, we found a pronounced and comparable inverted-U-shaped relationship between power in the alpha-and beta-bands, and pupil and pupil-derivative alike: The inverted-U alpha-band nonlinearity was most strongly expressed in inferior parietal cortex, extending towards higher visual as well as temporal regions in the beta band. In contrast to McGinley et al. (2015), we found no evidence 550 for an U-shaped relationship between pupil and high-frequency (50-100 Hz) activity.
The inverted-U shaped relationship between arousal and behavior may result from the differential activation of neuromodulatory receptors in particular noradrenergic adrenoreceptors with varying affinity (Berridge & Spencer, 2016) or nicotinic and muscarinic acetylcholine receptors (Bentley et al., 2011). Similar mechanisms may therefore underlie the 555 quadratic relationship observed here. More specifically, in the case of noradrenaline, the activation of high-affinity α2 receptors at medium arousal level may (directly or indirectly) lead to increased alpha-and beta-band activity. At elevated arousal levels, on the other hand, correspondingly high levels of noradrenaline lead to the activation of low-affinity α1 receptors, possibly resulting in a decrease of alpha-and beta-band activity. 560 From a functional perspective, an inverted-U relationship seems paradoxical. One would expect alpha activity to approach a minimum at intermediate states of pupil-linked arousal, i.e., show a U-shaped relationship instead, to optimally facilitate the processing of external sensory input (Jensen & Mazaheri, 2010). However, our participants were resting, while neither receiving nor expecting visual input. In this situation an inward focus of attention, 565 suppressing external visual information may have reversed the expected relationship (also see Hong et al., 2014).
Of note, we find evidence for both, linear positive and non-linear inverted-U shape relationships between alpha/beta activity and pupil-linked arousal. Cortical sources of the two effects coincide spatially (cp. Figures 4D & 5C). Previous research has suggested that 570 parieto-occipital cortices harbor at least two distinct alpha rhythms (Barzegaran et al., 2017;Keitel & Gross, 2016) with separable characteristics (Benwell et al., 2019) that may serve different functions (Capilla et al., 2014;Sokoliuk et al., 2019). Therefore, spatially overlapping linear and non-linear effects may reflect distinct influences of arousal on the generative processes underlying separate alpha rhythms. 575

Arousal modulation of cortical excitation-inhibition ratio
In addition to interactions between pupil-linked arousal and band-limited activity, we show that pupil diameter correlates systematically with aperiodic brain activity. Recent theoretical and empirical evidence suggest that the slope of the aperiodic component of neural power spectra, the spectral exponent, reflects the ratio between excitatory and inhibitory neuronal 580 processes in the underlying neuronal circuits (Colombo et al., 2019;Gao et al., 2017).
Here, we report a negative linear relation between pupil diameter and the spectral exponent. In other words, periods of elevated arousal were accompanied by a flattening of the power spectrum. In simulated neural networks, a shallower slope of the aperiodic component is related to increased excitation relative to inhibition (Gao et al., 2017;Trakoshis 585 et al., 2020). Thus, increased arousal (i.e., dilated pupils) may co-occur with increased E/I in distributed neuronal circuits, whereas periods of low arousal (i.e., constricted pupils) would indicate states of relatively decreased E/I. This interpretation is also consistent with previous findings arguing that noradrenaline increases cortical E/I (Pfeffer et al., , 2020, possibly through a decrease in inhibition (Froemke, 2015;Martins & Froemke, 2015). We found this 590 relationship of E/I and pupil-linked arousal in the majority of regions, thus suggesting a global and uniform effect. In sum, pupil-linked arousal may provide a window into yet another cortical circuit property, cortical E/I, relevant for cognitive computation (Cavanagh et al., 2020;Kosciessa et al., 2021;Lam et al., 2017;B. K. Murphy & Miller, 2003;Pettine et al., 2021;Pfeffer et al., 2020). 595

Conclusion
We exhaustively mapped how fluctuations in pupil-linked arousal and cortical activity covary in the human brain at rest. The diverse linear and non-linear relationships we describe in our data suggest profound influences of arousal on cortical states and challenge a global and unspecific role for arousal neuromodulation in cortical function. Instead, our results support 600 the view of specific neuromodulatory influences that differ depending on the targeted cortical region, can express in varying frequencies of band-limited activity -or as broadband effects, and vary in their timing. The present data provide the basis for studying how these influences change under cognitive task demands and when engaging in different behaviours. In fact, neuromodulatory effects are often neglected when studying the cortical population activity 605 underlying cognitive function -given a central role for arousal, our results argue for including it as a covariate in future studies.

Participants and Data Acquisition (Hamburg) 610
30 healthy human participants (16 women, mean age 26.7, range 20-36) participated in the study after informed consent. All included participants were healthy, with no current or previous diagnosis of psychiatric or neurological disorder (full list of exclusion criteria can be found in Pfeffer et al., 2020). The study was approved by the Ethics Committee of the Medical Association Hamburg. Two participants were excluded from analyses, one due to excessive 615 MEG artifacts, the other due to not completing all recording sessions (see below). Thus, we report results from N=28 participants of which 15 were women. The present dataset is a subset of a larger data set that entailed selective pharmacological manipulations (Pfeffer et al., , 2020. For the present article, only data from the placebo / resting state condition was analyzed. 620 Each participant completed three experimental sessions (scheduled at least 2 weeks apart), consisting of drug or placebo intake at two time points, a waiting period of 3 hours, and an MEG recording session. During the recordings, participants were seated on a chair inside a magnetically shielded chamber. Each recording session consisted of two resting state measurements as well as four blocks consisting of two variants of a behavioral task. 625 Each block was 10 minutes long and followed by a short break of variable duration.
Brain activity was recorded using a whole-head CTF 275 MEG system (CTF Systems, Inc., Canada) at a sampling rate of 1200 Hz. Eye movements and pupil diameter were recorded with an MEG-compatible EyeLink 1000 Long Range Mount system (SR Research, Osgoode, ON, Canada) and electrocardiogram (ECG) as well as vertical, horizontal and radial EOG was 630 acquired using Ag/AgCl electrodes. During all recordings, the participants were instructed to keep their eyes open and fixate a green fixation dot in the center of a gray background of constant luminance, which was projected onto a screen from outside the magnetically shielded recording chamber at a refresh rate of 60 Hz.

Participants and Data Acquisition (Münster)
Forty right-handed volunteers (21 women, age 25.1 ± 2.7 y (Mean ± SD), range 21-32 y) participated in the study. All participants reported having no respiratory or neurological disease and gave written informed consent prior to all experimental procedures. The study was approved by the local ethics committee of the University of Muenster. 640 Participants were seated upright in a magnetically shielded room while we simultaneously recorded 5 minutes of MEG and eye tracking data. MEG data was acquired using a 275 channel whole-head system (OMEGA 275, VSM Medtech Ltd., Vancouver, Canada) at a sampling frequency of 600 Hz. Pupil area of the right eye was recorded at a sampling rate of 1000 Hz using an EyeLink 1000 plus eye tracker (SR Research). During recording, participants 645 were to keep their eyes on a fixation cross centred on a projector screen placed in front of them. To minimise head movement, participants' heads were stabilised with cotton pads inside the MEG helmet.

Participants and Data Acquisition (Glasgow) 650
MEG resting-state data were acquired for 24 healthy, right-handed participants (9 women; age 23.5 ± 5.5 y (Mean ± SD), range 18-39 y). The study was approved by the local ethics committee (University of Glasgow, College of Science and Engineering; approval number 300140078). Participants gave written informed consent prior to testing.
Recordings were obtained with a 248-magnetometer whole-head MEG system (MAGNES 655 3600 WH, 4-D Neuroimaging), set to a sampling rate of 1017 Hz, and an EyeLink 1000 eye tracker (SR Research), sampling the pupil area of the left eye at 1000 Hz, both situated in a magnetically shielded room (Vacuumschmelze). Head position was measured at the start and end of the recording. To this end, a set of five coils were placed on the head of the participant. Coil positions and head shape were digitised using a FASTRAK stylus (Polhemus Inc., VT, 660 US). During recordings, participants sat upright and fixated on a green point (measure) projected centrally on a screen with a DLP projector (60 Hz refresh rate) for a minimum of 7 minutes (420 s).
Resting state recordings were part of a study into audio-visual speech processing that was recorded in two sessions. Audio-visual data have been published elsewhere (Keitel et 665 al., 2018(Keitel et 665 al., , 2020. In this study, the resting state block was recorded as the last block of the second session after 7 -9 experimental blocks, or more than 65 min into the session. Note that the resting state data have not been reported in any of the previous publications.
Data of two participants were excluded from analysis after screening recordings, one due to strong intermittent environmental artifacts in their MEG recordings and another due to a 670 technical fault in tracking the pupil.

MEG Preprocessing (Hamburg)
The sensor-level MEG signal was cleaned of extra-cranial artifacts using semi-automatic artifact detection routines implemented in the FieldTrip toolbox . 675 Broken and noisy recordings channels were identified through visual inspection and were marked as artifactual across all channels (+/-500 ms). Subsequently, the data were downsampled to 400 Hz split into low ([0.5-2]-40 Hz; the lower cutoff was variable across (but identical within) subjects at 0.5, 1 or 2 Hz) and high (>40 Hz) frequency components, using a 4th order Butterworth filter. The separate data segments were independently 680 submitted to Independent Component Analysis (ICA; Hyvarinen, 1999). The split into lowand high-frequency data aimed to facilitate the detection of sustained muscle artifacts especially present in the high-frequency range. Artifactual independent components were identified through visual inspection of their topography and power spectra. On average, 23 ± 14 components (Mean ± SD) were subtracted from the raw data (placebo condition only). 685 After artifactual components were removed, the low-and high-frequency segments were combined into one single data set.

MEG Preprocessing (Münster)
MEG preprocessing was done in Fieldtrip  for MATLAB (The 690 MathWorks Inc). First, noisy recording channels were identified by visual inspection and rejected from further analyses. Next, we adapted the synthetic gradiometer order to the third order for better MEG noise balancing (ft_denoise_synthetic). Power line artefacts were removed using a discrete Fourier transform (DFT) filter on the line frequency of 50 Hz and all its harmonics (including spectrum interpolation; dftfilter argument in ft_preprocessing). Next, 695 we applied independent component analysis (ICA) on the filtered data to capture eye blinks and cardiac artefacts (ft_componentanalysis with 32 extracted components). On average, artefacts were identified in 2.35 ± 0.83 components (M ± SD) per participant and removed from the data. Finally, cleaned MEG data were downsampled to 400 Hz.

MEG Preprocessing (Glasgow)
For MEG and eye tracking data preprocessing and analysis, we used Matlab (The MathWorks Inc), involving in-house MATLAB routines and the FieldTrip toolbox . The MEG signal was resampled to the sampling rate of the eye tracker (1000 Hz) and denoised using in-built MEG reference sensors. We rejected between 2-27 noisy 705 MEG channels (median = 7.5 ± 4 absolute deviation) by visual inspection using FieldTrip's ft_rejectvisual. Continuous MEG recordings were then screened for noisy segments. Highpass filtered (2 Hz cutoff, 4th-order Butterworth, forward-reverse two-pass) MEG time series were subjected to an independent component analysis (ICA; Hyvarinen, 1999). Prior to the ICA, a principal component analysis (PCA) projected MEG data into a 32-dimensional 710 subspace. For each participant, we visually identified 2-4 components (median = 3) capturing eye (blinks, movements) and heart-beat artefacts. Finally, resampled continuous MEG recordings were projected through the subspace spanned by the remaining components.
MEG and eye tracker time axes were realigned by cross-correlating traces of horizontal and vertical eye movements that were simultaneously recorded by the eye tracker and the 715 MEG. The lag corresponding to the maximum of the cross-correlation functions of horizontal and vertical eye movement traces was used as offset for the realignment.

Pupil preprocessing (all sites)
Pupil area traces were converted to pupil diameter to linearise our measure of pupil size. 720 Using additional Matlab routines (available from https://github.com/anne-urai/ pupil_preprocessing_tutorial, as used in (Urai et al., 2017), blinks were identified by an automatic (and visually validated) procedure and linearly interpolated. In a second pass, a similar procedure was used on smaller eye tracker artifacts (spikes, jumps). Next, canonical responses to blinks and saccades were estimated and then removed from pupil time series 725 (Hoeks & Levelt, 1993;Knapen et al., 2016;Wierda et al., 2012). To that end, pupil time series were band-pass filtered (pass band: 0.005 -2 Hz, 2nd-order Butterworth, forward-reverse two-pass), then downsampled to a sampling rate of 400 Hz. Finally, we computed the firstorder derivative of each band-pass filtered and down-sampled pupil trace.

Spectral decomposition of pupil time series (all sites)
The spectral content of pupil traces was evaluated using the multi-taper approach as implemented in FieldTrip. We used the default settings, except for a spectral smoothing of 735 0.035 Hz, and zero-padding time series to 2 19 (minimum power of 2 that exceeded the length of any pupil trace) to unify frequency resolutions across centres and traces of different lengths. Power spectra were calculated on standardized (z-scored) pupil traces for the 0.005-2 Hz range, sampled in logarithmic steps.
Resulting power spectra were also subjected to the spectral parameterization toolbox 740 ('fooof'; Donoghue et al., 2020) to extract Hippus peak frequency as well as the exponent describing the spectral slope as a means to characterise the aperiodic activity. We used the 'fooof' MATLAB wrapper with the default settings, except setting the maximum number of peaks ('max_n_peaks') to 5, and the peak detection threshold ('peak_threshold') to 0.5.

MEG spectral analysis (all sites) 745
In order to derive spectral estimates from the sensor-level data, we followed an approach outlined previously (Hipp et al., 2012), based on Morelet's wavelets (Tallon-Baudry & Bertrand, 1999): Following Hipp et al. (2012), we derived spectral estimates for 25 logarithmically spaced center frequencies, covering a broad frequency range from 2 Hz to 128 Hz. Spectral estimates were derived for successive temporal windows with an overlap of 80%. Windows that contained periods marked as artifactual (see Preprocessing) were discarded from all further analyses. From the spectral estimates, sensor-level power envelopes were obtained 755 through: ( , ) = | 89:8 ( , ) | ; (Eq. 2) where 89:8 denotes the complex spectral estimates for segment and frequency . 760

MEG source reconstruction (all sites)
Accurate source models were generated using the FieldTrip interface to SPM8 (Litvak et al., 2011), and the Freesurfer toolbox (Fischl, 2012). Source models were based on individual T1-weighted structural magnetic resonance images (MRIs), acquired on a 3T Siemens Tim Trio system with a 12-channel head coil (Glasgow), a Siemens 3T Prisma with a 20-channel 765 head coil (Münster), and a Siemens Magnetom Trio (Hamburg). Anatomical scans were coregistered with the MEG coordinate system using a semiautomatic procedure (Gross et al., 2013), then segmented and linearly normalised to a template brain (MNI space). We used a single shell as the volume conduction model (Nolte, 2003).
We projected sensor-level data into source space using frequency-specific DICS (Dynamic 770 Imaging of Coherent Sources) beamformers (Gross et al., 2001) with a regularisation parameter of 5% and optimal dipole orientation (singular value decomposition method). As the source model we used a volumetric grid. Grid points had a spacing of 5 mm, resulting in 8799 dipoles covering the whole brain. Source level power envelopes were computed from the source-projected complex signal: 775 where A denotes the spatial filter, 8<= the source-level power envelopes and 89:8 the sensor-level analytic signal. The 5 mm spacing was chosen so as to allow for an integer spatial downsampling of the Brainnetome atlas (Fan et al., 2016) which we used to map 780 cortical (and subcortical) regions into 246 distinct areas.

Mutual information between pupil and spectral activity in sensor space
To quantify global pupil-power coupling in sensor space we calculated mutual information (MI) between frequency-specific power envelopes and pupil time traces. To this end, pupil time series were downsampled to the temporal resolution of the respective power envelope, 785 and segments containing artifacts in power envelopes were also removed from the pupil traces. After normalising both time series by means of a Gaussian-copula based approach (Ince et al., 2017), mutual information was computed per frequency and MEG sensor.
Observed MI was evaluated against permutation distributions of surrogate MIs based on randomly time-shifting pupil traces 200 times, while avoiding a window of +/-10 sec around 790 the match of both original time series. Using the log-transformed observed MI, and the mean and standard deviation of equally log-transformed permutation distributions we computed zvalues for each frequency and sensor. In a final step, these z-values were pooled across sensors and their mean across participants tested against zero, for each recording site separately. Test results are reported for uncorrected and FDR-corrected thresholds in Figure  795 2.

Cross-correlation analysis (sensor space)
In order to allow for the comparison across MEG recording sites (with different sensor layouts), we first ordered the recording channels into 39 equally sized, non-overlapping spatial bins, sorted from anterior to posterior channels. We next averaged the cross 800 correlograms per bin and additionally averaged the cross correlations across three distinct frequency ranges: 2- 4 Hz,[8][9][10][11][12][13][14][15][16] Due to the varying width of the wavelets for each center frequency, the temporal resolution, and, therefore, the number of samples, varied across frequencies. Hence, prior to averaging, we linearly interpolated the frequencyspecific cross-correlograms such that the number of samples was equal to the number of 805 samples of the highest frequency of each band (i.e., 4 Hz,16 Hz and 128 Hz).

Source-level power spectra and spectral parametrization
In order to correlate fluctuations in pupil diameter (and its temporal derivative) with fluctuations in the slope of the aperiodic component of the power spectra (henceforth referred 810 to as the 'spectral exponent') we first computed source-level power spectra. To this end, we source-projected the broad-band signal using linear beamforming (LCMV; Van Veen et al., 1997) instead of DICS (Gross et al., 2001). The procedure was similar to DICS, with the difference that the spatial filters were obtained using the cross-spectral density matrix averaged across all 25 frequency-bands of interest. In keeping with the previous analysis, we 815 used a regularization parameter of 5%.
Next, we divided the resulting source-level broad-band signal into non-overlapping temporal segments of 2s length. For each segment, power spectra were computed using MATLAB's 'pwelch' function, ranging from 2 Hz to 128 Hz, with a frequency resolution of 0.5 Hz. In addition, for each temporal segment, the mean pupil diameter as well as the mean of 820 the pupil-derivative were computed.
To separate aperiodic and periodic components of the power spectra, we utilized the spectral parameterization toolbox 'fooof' for Python 3.7 (Donoghue et al., 2020). The algorithm iteratively fits neuronal power spectra as the sum of a Lorentzian function and a variable number of Gaussians of varying width, height and mean (corresponding to the center 825 frequency). To constrain the fitting procedure, we limited the number of periodic components to 6 and defined a minimum peak height of 0.05, with the minimum and maximum bandwidth of each peak set to 1 Hz and 8 Hz. The knee of the Lorentzian was set to zero.
The algorithm was applied to the power spectra of each individual temporal segment (see above), which provided time-resolved estimates of periodic and aperiodic components. While 830 the slope of the aperiodic component is relatively stable for lower frequencies up until ~40-50 Hz, the rate at which power changes can differ substantially for higher frequencies. This means that fitting power spectra uniformly across the entire frequency range may result in poorer fits and, consequently, misleading parameter estimates. Thus, and largely consistent with previous work (Pfeffer et al., 2020;Waschke et al., 2021), neuronal power spectra were 835 fitted in the frequency range from 3 Hz to 40 Hz.

Polynomial model
We quantified non-linear components of the relation between pupil diameter (and pupil derivative) and band-limited activity fluctuations in spectral power. To this end, we sorted the aforementioned temporal segments (see section above: 'Source-level power spectra and 840 spectral parametrization') based on pupil diameter (or pupil derivative) into 14 nonoverlapping and equidistant bins. Next, we computed the average power for each of the bins and normalized (z-scored) the obtained values. Quadratic relationships were quantified by fitting the power vs. pupil (or pupil-derivative) function with a polynomial of degree 2, separately for each of the 25 MEG center frequencies: 845 where denotes the source location (with being 1 to 246 regions of the Brainnetome atlas), is the MEG center frequency (from 2 to 128 Hz, in 25 logarithmic steps) and is the mean 850 pupil diameter for each pupil bin. The coefficients for the quadratic and the linear component are denoted with 2 and 1 , respectively, and 0 describes the offset. The polynomial was fitted by minimizing the sum of the squares of the difference between the polynomial model and the power vs. pupil (pupil-derivative) functions (as implemented in MATLAB's 'polyfit' function). 855

Acknowledgments
TP and THD thank Christiane Reissmann and Karin Deazle for assistance with MEG recordings and the acquisition of structural brain scans. CK, AK, GT and JG thank Gavin 1175 Paterson and Frances Crabbe for expert assistance with MEG recordings and structural brain scans. CK thanks Linbi Hong, Josef Faller, and Paul Sajda for comments on the research. DSK and JG thank Hildegard Deitermann, Ute Trompeter, and Karin Wilken for assistance with MEG recordings and structural brain scans. All authors thank Ruud L. van den Brink for helpful comments on the manuscript.   dots denote P<0.05 after adjusting for false discoveries (q = 0.1; two-sided paired t-test). The gray shaded area depicts the standard error of the mean across subjects. Positive and negative values are indicative of a U-shaped and inverted-U-shaped relationship, respectively (C) Normalized β2 coefficient, averaged across all subjects, separately for four regions of interest: visual cortex, auditory cortex, somatosensory cortex and anterior cingulate cortex (from left to right). Insets show the approximate extent of the regions of interest. (D) Spatial distribution of 1280 the coefficient of the quadratic term (see Methods for details), at center frequencies 11 Hz (left) and 19 Hz (right). Masked at P = 0.05 (two-sided paired t-test, adjusted for a false discovery rate of 0.1).