Harmonic Amplitude Summation for Frequency-tagging Analysis

Abstract In the approach of frequency tagging, stimuli that are presented periodically generate periodic responses of the brain. Following a transformation into the frequency domain, the brain's response is often evident at the frequency of stimulation, F, and its higher harmonics (2F, 3F, etc.). This approach is increasingly used in neuroscience, as it affords objective measures to characterize brain function. However, whether these specific harmonic frequency responses should be combined for analysis—and if so, how—remains an outstanding issue. In most studies, higher harmonic responses have not been described or were described only individually; in other studies, harmonics have been combined with various approaches, for example, averaging and root-mean-square summation. A rationale for these approaches in the context of frequency-based analysis principles and an understanding of how they relate to the brain's response amplitudes in the time domain have been missing. Here, with these elements addressed, the summation of (baseline-corrected) harmonic amplitude is recommended.


INTRODUCTION Frequency Tagging
It has long been known that a stimulus presented at a periodic rate elicits a response from an observer's brain at exactly that rate. For example, a light flickering on and off at a periodic rate, 14 times a second, elicits a measurable response in the EEG of a human observer 14 times a second (Adrian & Mathews, 1934). In the time domain, a response is evident as periodic changes in the brain's response amplitude across time. Following Fourier transformation into a frequency domain representation (Danielson & Lanczos, 1942;Fourier, 1822), the response is evident as a high amplitude "peak" at exactly the fundamental stimulus presentation rate (frequency = F ) and/or its higher harmonics, that is, at frequencies that are integer multiples of F (2F, 3F, etc.; Regan, 1966Regan, , 1989. The approach of presenting stimuli and analyzing neural responses at the frequency of stimulation is referred to by many names: "Frequency tagging" (Srinivasan, Russell, Edelman, & Tononi, 1999;Tononi, Srinivasan, Russell, & Edelman, 1998) is the one that will be used here. Other names for this approach differ mainly on their point of reference: to the responses that appear consistently periodic to stimuli presented at high rates, that is, "steady-state" responses, for example, "steady-state visual-evoked potentials" (Norcia, Appelbaum, Ales, Cottereau, & Rossion, 2015;Heinrich, 2010;Di Russo, Teder-Salajarvi, & Hillyard, 2002;Regan, 1966Regan, , 1989 and "auditory steady state potentials/responses" (Watkin, 2008;Geisler, 1960) or "traveling wave" responses (Engel, Glover, & Wandell, 1997); to the "stimulation" mode itself ("fast periodic visual stimulation"; Rossion, 2014;Rossion, Retter, & Liu-Shuang, 2020); or the "analysis" occurring in the frequency domain ("Fourier analysis/synthesis"; Zhou, Melloni, Poeppel, & Ding, 2016;Bach & Meigen, 1999;Movshon, Thompson, & Tolhurst, 1978; or simply "frequency (domain) analysis," for example, as in McKeefry, Russell, Murray, & Kulikowski, 1996). Despite the varying terminologies, the principles of the approach are the same. In a similar vein, various types of stimulation modalities (visual, auditory, somatosensory, cross-modal) and recording methods (EEG, electroretinogram, fMRI, single-cell recordings, etc.) may be applied with various participant groups (human adults, children, infants, nonhuman primates, cats, rodents, frogs, insects, etc.), resulting in some practical differences, but the same fundamentals, of the approach.
In (cognitive) neuroscience research, the frequencytagging approach is associated with undeniable advantages. As noted early on, this approach is well suited for specifically relating brain processes to external events: "This gives a method of tracing the visual messages in the brain, for by means of the flicker rhythm they can be made easy to recognize" (Adrian, 1944, p. 361). More recently, its objectivity and sensitivity (i.e., high signalto-noise ratio) have been highlighted, and the use of the paradigm is undoubtedly on the rise, having been extended from the study of basic sensory processes and their modulation by spatial/selective attention to the direct measurement of higher levels of cognition in recent years (see Norcia et al., 2015, for a review). However, frequency tagging is still fundamentally limited by outstanding conceptual and methodological ambiguities in dealing with responses occurring across harmonics.

Higher Harmonics
Another way of describing frequency tagging is the following: Given a periodic stimulus, responses of the brain periodic to that stimulus are investigated. In this formulation, it is evident that the brain's responses may occur at the rate of stimulation, F, but also at the other rates periodic to the stimulation: the higher harmonics (2F, 3F, etc.). For example, a stimulus modulated 8 times a second, at 8 Hz, may generate responses that are evident as amplitude peaks in the frequency domain representation of the brain recording at 8 Hz (F, the first harmonic, corresponding to the fundamental frequency 1 ), but also at 16 Hz (2F, the second harmonic) and 24 Hz (3F, the third harmonic). Because only responses at higher harmonics are periodic to the fundamental frequency, it is uniquely at the higher harmonics, rather than at a diffuse band, that higher frequency constituents of frequencytagged brain responses are present.
Although responses of the brain are not always generated at the higher harmonics, they often do occur Zhou et al., 2016;Norcia et al., 2015;Rossion, 2014;Heinrich, 2010;Vialatte, Maurice, Dauwels, & Cichocki, 2010;Bach & Meigen, 1999;Regan, 1966). Note that responses are not always generated at F either; for a classic example, in the case of alternating symmetrical stimulus inputs (e.g., pattern-reversing checkerboards), the brain responds only at 2F and higher even harmonics (Cobb, Morton, & Ettlinger, 1967;reviewed in Norcia et al., 2015; for different examples, Zhou et al., 2016;Heinrich, 2010;Movshon et al., 1978). Furthermore, note that, throughout this article, only harmonics that are specific to their fundamental frequency are addressed, which is always the case when a single stimulus presentation frequency is used (but for an extension to other cases, with more complex stimulation paradigms, please see the Which Harmonics to Consider? section).
How often is there considerable amplitude at the higher harmonics? Can the cost of omitting higher harmonics in published studies be evaluated? Unfortunately, most studies do not report whether or not there were responses at higher harmonics, as mentioned above. Moreover, when higher harmonic responses were reported to be present, they were often not described, for example, "Peaks were also present at the harmonics of the stimulus frequency but were not analyzed in this study" (Srinivasan et al., 1999, p. 5438); "Higher harmonics may play a role, especially at lower temporal frequencies (see, for example, the double peaks in the 12-Hz data in Figure 1), but these are not considered here" (Kremers, Rodrigues, Silveira, & da Silva Filho, 2010, p. 579); "Note however that 2 Hz is a harmonic of 1 Hz and may actually be a relevant spectral region to consider (albeit outside the scope of Figure 1. A demonstration of the importance of considering higher harmonics in frequency-tagged response analyses. (A) Two synthetic periodic signals, each composed of five harmonic frequencies. (B) In the frequency domain, a consideration of only one harmonic (at the fundamental frequency) describes Signal 2 as larger than Signal 1. This description is not in agreement with typical time domain response analyses (e.g., peak amplitudes) to compare these signals. Note that this figure will be revisited (expanded) in the Interpreting Harmonics section. this report)" (Kosem, Gramfort, & Van Wassenhove, 2014, Figure S2; evident in figures but not discussed or included in analyses: Schettino, Porcu, Gundlach, Keitel, & Müller, 2020;Eidelman-Rothman et al., 2019;Chadnova et al., 2018;Winawer et al., 2013;Hönegger et al., 2011;Katzner et al., 2009;Pastor, Artieda, Arbizu, Valencia, & Masdeu, 2003;Pastor et al., 2002;Kaspar, Hassler, Martens, Trujillo-Barreto, & Gruber, 2001;Müller et al., 1998), preventing a wide-scale review.
In some cases, the conclusions of studies considering and not considering higher harmonics can be compared. For example, in frequency-tuning studies, without considering higher harmonics, maximal visual responses were reported to stimuli modulated at about 10-15 Hz with EEG (e.g., Ding, Sperling, & Srinivasan, 2006;Pastor et al., 2003;Regan, 1966; see also Vialatte et al., 2009). However, when higher harmonics were considered, the lowest stimulation frequency tested (3 Hz) yielded the maximal visual EEG responses, being over 3 times higher than the responses to 12 Hz stimulation with natural images   Figure S3). Similarly, maximal auditory responses were reported to stimuli modulated at about 40 Hz with EEG (e.g., Pastor et al., 2002;Ross et al., 2000;Galambos, Makeig, & Talmachoff, 1981), but a consideration of higher harmonics produced the maximal auditory EEG responses at the lowest stimulation frequency tested (0.75 Hz; Tlumak et al., 2011, Figure 4). Overall, although it is thus impossible to ascertain what the impact of unreported or uncharacterized higher harmonic responses in most studies may have been, it is likely that it was often considerable.

Should Higher Harmonic Responses Be Combined, and If So, How?
At present, harmonics are surrounded by many questions: Why do they occur? What do they represent? Which, or how many, harmonics should be considered? Should they be taken into account for response identification and measurement, and if so, how? Indeed, the lack of understanding and standard practice regarding higher harmonics has limited the ease (i.e., objectivity) of frequency-tagged response identification and measurement. This is particularly significant because objectivity is given as a primary advantage of the frequency-tagging technique, contributing to its increasing application in (cognitive) neuroscience research and clinical applications (e.g., see Rossion et al., 2020;Norcia et al., 2015).
In the following, a validated methodology for combining (baseline-corrected) harmonic amplitudes through simple summation will be provided. This approach derives from a theoretical basis of how signals over time are represented through mathematical transformations into the frequency domain (see the Frequency Domain Representations section), extended to experimental responses in practice (see the Frequency-tagged Responses in the Frequency Domain section). This approach was indicated empirically by Retter and Rossion (2016a), and it has since been applied in a number of studies, however, primarily by those authors or associated research groups (e.g., van de Walle de Ghelcke, Rossion, Schiltz, & Lochy, 2021;Damon, Leleu, Rekow, Foncet, & Baudouin, 2020;Fisher, Towler, Rossion, & Eimer, 2020;Gwinn & Jiang, 2019;Dwyer, Xu, & Tanaka, 2019;Van der Donck et al., 2019;Beck, Rossion, & Samson, 2018;Chemin, Huang, Mulders, & Mouraux, 2018;De Keyser, Mouraux, Quek, Torta, & Legrain, 2018;Guillaume, Mejias, Rossion, Dzhelyova, & Schiltz, 2018;Gwinn, Matera, O'Neil, & Webster, 2018;Leleu et al., 2018;Xu, Liu-Shuang, Rossion, & Tanaka, 2017). To be of further use to the scientific community, the approach requires deeper methodological evaluation and, especially, evaluation in a theoretical context, which is the goal throughout this article. From this, some practical guidelines are offered (see the Combining Harmonic Responses section) and implications are drawn for the interpretation of harmonic responses more generally (see the Interpreting Harmonics section).

FREQUENCY DOMAIN REPRESENTATIONS Sine Waves
When a signal is transformed into the frequency domain (by means of a Fourier transform), it becomes represented through a combination of sine waves, which are the fundamental units of the frequency domain. Although there are many texts on the mathematics of frequency transformations and representations (e.g., Gonzalez & Woods, 2018;Forinash & Christian, 2016;Patel, 2012;Strang, 2007;Smith, 1997;Press, Falnnery, & Teukolsky, 1993), a basic understanding of sine waves and their combination is a sufficient foundation for the interpretation of multiharmonic responses of the brain (Regan, 1989).
Briefly, sine waves are trigonometric functions that describe periodic signals in terms of frequency, amplitude, and phase ( Figure 2). The frequency of a sine wave describes the number of cycles (of 360°; equivalent to 2π radians) per unit of time or space (time is typically given in units of cycles/sec = sec −1 = Hertz = Hz). Note that the cycles of sine waves are periodic and could repeat their pattern infinitely, as a circle could be endlessly traced (the sine wave, as in Figure 2B, derives from the y-axis values of a unit circle, as shown in Figure 2A). The amplitude of a sine wave is defined along the y-axis (e.g., in Figure 2B, the sine wave, spanning from −1 to 1, has an amplitude of 1; its unit in EEG recordings is typically microvolts). The phase of a sine wave is a measure of its starting angle (indicated by theta in Figure 2A), with an arbitrary beginning at zero (in units of degrees or radians), as is shown in Figure 2B. Changes in frequency, amplitude, and phase, which help demonstrate these properties, are illustrated in Figure 2C. For those whom it helps to see it mathematically, the expression of a sine wave, as a function of x, is: y(x) = asin(2πfx + ϕ), where a = amplitude (scaling on the y-axis), f = frequency (by cycles), and ϕ = phase (x-axis shifts).

A Lot of Sine Waves
A frequency domain representation of a signal is essentially a lot of sine waves. That is, when a signal is transformed into the frequency domain, the resultant x-axis describes the frequency of its constituent sine waves. The other descriptors of sine waves, amplitude, and phase are described in the transformed, complex-valued y-axis at each frequency, which is typically plotted as separate amplitude and/or phase frequency spectra. 2 The resolution (x-axis sampling) of the frequency domain spectrum is the inverse of the signal recording length and the range spans from zero to half of the signal sampling rate (note that these properties have practical implications for frequency-tagging experimental design, e.g., as addressed in Bach & Meigen, 1999). The combination, through summation, of these sine waves described in the frequency domain reconstructs the original signal in the time domain. Here, the focus will be on periodic signals over time, but note that frequency domain analyses can be applied in many settings (e.g., signals over space or over two dimensions).
In the simplest case, a periodic signal that is a perfect sine wave is represented in the frequency domain by a single frequency, representing a single sine wave ( Figure 3A). Another way to understand this is to observe that, in this case, the frequency, amplitude, and phase of a single sine wave in the frequency domain are sufficient to reconstruct the original signal in the time domain. In most cases, signals are more complex (i.e., nonsinusoidal), but this does not pose a problem: A combination of sine waves at different frequencies can sum to model any signal. In a classic example, a periodic squarewave signal is shown to be represented with a sum of sine waves specific to its periodicity ( Figure 3B). Nonperiodic signals, for example, ERPs to temporally jittered stimuli, can also be represented in the frequency domain, but because they are not specific to limited frequencies, their interpretation does not correspond to that of frequency-tagged signals ( Figure 3C). Although frequency-based analyses of nonperiodic signals may be applied (e.g., see Herrmann, Rach, Vosskuhl, & Strüber, 2014;Başar & Schürmann, 1994;Regan, 1989; see also Chemin et al., 2018), these are outside the present focus on frequency tagging.
It may be observed that a simple sine wave ranging from −1 to 1 in the time domain has an amplitude of 1 in the frequency domain, but that the relationship between the time domain and frequency domain amplitudes for multiharmonic signals is more complex (compare Figure 3A and B). However, there is a direct relationship between these dimensions: Because energy is conserved from the time to the frequency domain, the sum of the root-mean-square amplitudes of the time domain signal equals the sum of the squared root mean amplitudes of the frequency domain signal (Parseval's relation; Parseval des Chênes, 1806; e.g., see Smith, 1997). 3 For example, the sum of the squared amplitudes per cycle of the time domain signal in Figure 3A is equal to 0.5, and the sum of the squared root-mean-square amplitudes of its discrete frequency domain signal is equal to 0.5. Multiharmonic signals also preserve this relationship, although their time domain amplitude range does not directly relate to their frequency domain amplitude (being affected by phase; see the What About Phase? section).

FREQUENCY-TAGGED RESPONSES IN THE FREQUENCY DOMAIN One Harmonic; A Lot of Harmonics
Frequency tagging is an approach in which stimuli are presented periodically to generate periodic responses of the brain that can thus be identified in the frequency domain at specific frequencies harmonic to the stimulation (i.e., Figure 2. The sine wave. (A) A unit circle, with a radius of 1, illustrates the underlying trigonometry of sine (and cosine) functions. (B) A sine wave derives from the y-axis values of points on the unit circle as a function of angle (here, this is emphasized with corresponding angle colors). A complete cycle of a sine wave contains 360°. (Note that a cosine wave, which derives from the x-axis values, has the same shape but with a 90°phase shift.) (C) Sine waves are described in terms of frequency (cycles per unit), amplitude ( y-axis scale), and phase (x-axis shift). The given examples exhibit sequential changes in these properties (beginning in reference to B, then from left to right; the changed property is indicated in bold). Amp. = amplitude; Freq. = frequency. the fundamental and the higher harmonic frequencies). In the following, examples will be taken only for harmonics that are specific to a single tagged frequency: Again, for determining specific harmonics in the context of multiple tagged frequencies, see the Which Harmonics to Consider? section).
According to the principles of frequency analyses, a simple, sinusoidal brain response would be represented only at the fundamental frequency F, whereas more complex brain responses would be represented with a combination of F and its higher harmonics, 2F, 3F, and so on. Generally, this is evidenced with experimental brain responses. In the event that the brain responses are nearly sinusoidal, the response is dominated by amplitude at F (e.g., at a high stimulus presentation rate: Figure 4A); in the event that the responses of the brain are complex, a combination of sine waves at different frequencies (i.e., the higher harmonics) can sum to model any signal ( Figure 4C). In many studies, complex, nonsinusoidal responses of the brain evoked over time are represented in the frequency domain not only at F, but with considerable amplitude at its higher harmonics (as addressed in the Introduction, e.g., Rossion et al., 2020;Retter & Rossion, 2016a;Zhou et al., 2016;Norcia et al., 2015;Alonso-Prieto et al., 2013;Tlumak et al., 2011;Heinrich, 2010;Vialatte et al., 2009Vialatte et al., , 2010Kremers & Scholl, 2001;Ross et al., 2000;Bach & Meigen, 1999;Sieving, Arnold, Jamison, Liepa, & Coats, 1998;Regan, 1966Regan, , 1989Donker, 1975;Van der Tweel & Verduyn Lunel, 1965;Brazier, 1964).
Indeed, higher harmonic responses may be accounted for in relation to the complex (i.e., nonsinusoidal) responses of the brain, in accordance with the principles of frequency domain analysis of periodic signals (as in Rossion et al., 2020;Zhou et al., 2016;Norcia et al., 2015;Heinrich, 2010;Regan, 1989). This account explains that higher harmonics are present when complex brain responses are present, but does not implicate a specific source of complex brain responses (see Heinrich, 2010). However, it is important to note that complex brain responses are not a product of frequency tagging and may equivalently occur with nonperiodic (event-related) stimulus presentation modes. 4 The harmonics do not represent new information, specific to the frequency domain: They are merely highlighted in an alternative, frequency domain representation of the original time domain signal (certain variations may be represented more or less clearly in each domain). These domains are interchangeable: As time domain data can be transformed into the frequency domain, frequency domain data can also be inversely transformed back into the time domain; indeed, time domain responses can be reconstructed from frequency domain harmonic amplitudes and phases (e.g., Ruhnau, Keitel, Lithari, is represented with a combination of many specific harmonic frequencies (lines colored correspondingly across top and bottom panels). Literally, the sum of these (and higher, not illustrated) colored lines' amplitude at each time point reconstructs the original signal. (C) A nonperiodic, ERP signal is represented with a combination of many nonspecific frequencies (data from Retter et al., 2020). Note several properties of the frequency domain signal: (1) the 0 frequency bin reflects the mean amplitude (DC offset) of the signal; (2) the x-axis resolution is the inverse of the signal recording duration; (3) although the frequency domain is plotted only until 10 Hz here, its range spans further (up to half of the signal sampling rate); and (4) although only the phase of tagged frequencies is indicated on the lower row here, each frequency in the spectrum has a corresponding phase value. Weisz, & Neuling, 2016;Sieving et al., 1998). The next section focuses on the combination of harmonic responses; however, individual harmonics will be further addressed in the Interpreting Harmonics section.

Imperfect Signals: Accounting for Baseline Noise
In theoretical examples of signal transformation into the frequency domain (as in the Frequency Domain Representations section), the signal is pure signal. In frequency tagging, as in all brain recordings, the signal (here, i.e., the responses of the brain at the tagged frequencies) also carries "noise," a term that refers to both non-event-related brain activity and artifacts (e.g., Luck, 2005;Regan, 1989). Note that there are many methods for correcting for noise (i.e., isolating signal) in frequency-tagging research, although a discussion of these is more general than the scope of this article (see instead, e.g., Meigen & Bach, 2000; Appendix 2 of Norcia et al., 2015). In the examples given here, a simple correction for noise will be applied that subtracts a local baseline from the amplitude of the frequencies of interest (e.g., Retter & Rossion, 2016a). 5 The baseline is defined as the mean amplitude of a symmetrical range of neighboring (i.e., continuously adjacent 6 ) frequency bins (for theoretical justification, see Norcia et al., 2015;Regan, 1989; e.g., Boremanse, Norcia, & Rossion, 2013;Meigen & Bach, 2000;Peterzell & Norcia, 1997; with power rather than amplitude: Mouraux et al., 2011;Vialatte et al., 2009;Srinivasan et al., 1999). This method is used to provide a measure of signal amplitude in the frequency domain that is relatable (i.e., both in the same unit) to amplitude in the time domain and with a noise level at zero, while compensating for local variations of noise inherent to brain recordings across the frequency spectrum.

COMBINING HARMONIC RESPONSES Combining Harmonic Response Amplitude
The combination of sine waves is simple: Sine waves sum linearly to reconstruct a signal. However, with the goal of identifying and measuring overall response amplitude in the frequency domain, because sine waves carry both amplitude and phase information, their sum is not intuitive to interpret in terms of amplitude only (or phase only). Note that there are alternative approaches for combining harmonics that incorporate both amplitude and phase; however, these approaches make use of phase as an indicator of reliability (coherence), typically across short stimulation durations (e.g., Delorme & Makeig, 2004;Strasburger, 1987;Jervis, Nichols, Johnson, Allen, & Hudson, 1983). Furthermore, although the amplitude of a single sine wave in the frequency domain directly relates to its time domain peak amplitude, that is, half its positive to negative peak range (e.g., Figure 3A; similarly, see Figure 4A), the frequency domain amplitude of complex time domain signals (summed sine waves) is not as easily visualized from the time domain. Perhaps for these reasons, various approaches have been taken for the combination of frequency-tagged multiharmonic brain response amplitude in the frequency domain (e.g., as mentioned previously, averaging or root-mean-square summing, respectively, Liu-Shuang et al., 2014;Hou et al., 2003).
The summation of harmonic amplitudes is recommended here for identifying and measuring the overall brain response (based on Retter & Rossion, 2016a; see also Heinrich, 2009). In the study of Retter and Rossion (2016a), this approach was validated empirically by qualitative comparison of time and frequency domain responses, in the situation where several equivalent time domain EEG responses were produced by several slow target stimulus presentation frequencies. There, it was observed that despite different distributions of harmonic amplitudes stemming from the different fundamental target stimulus presentation frequencies (1.1-2.5 Hz), the summation of baseline-subtracted harmonic amplitude across a common frequency range led to equivalent overall amplitudes that related to approximately equivalent response amplitude peaks in the time domain by visual inspection (see Figure 5A and B, Row 1, here, for examples of reprocessing of that data in combination with the underlying harmonic distributions). Moreover, a faster stimulus presentation rate (4.2 Hz), which produced visually lower amplitude deflections in the time domain ( Figure 5C, Row 1), also produced a lower summed harmonic response amplitude.
Here, these data are revisited quantitatively, with typical time domain interpretations of response amplitudes: peak-to-peak amplitude of the largest deflections and the area under the curve of the response deflections ( Figure 5, Row 4). Note that an exact comparison of specific deflections, as is more commonly done in relating brain responses, is possible for the conditions at 1.1 and 1.4 Hz, but that a different response pattern is observed for the condition at 4.2 Hz, preventing such a direct comparison.
The summed harmonic response amplitude is shown to be congruent with these measures ( Figure 5D). Critically, other approaches for harmonic combination would not have led to these conclusions when comparing conditions. Here, in the frequency domain, a large fundamental harmonic amplitude relates to a smaller number of harmonic responses (with an amplitude above 0.1 μV). Therefore, for example, averaging the harmonic responses would have generated lower amplitude responses the slower the stimulus presentation frequency (1.1 Hz < 1.4 Hz < 4.2 Hz; Figure 5D). For another example, the root-mean-square harmonic amplitude would have generated the highest amplitude for the highest stimulus presentation rate (4.2 Hz; Figure 5D). For a last example, using non-baseline-corrected amplitudes would have produced a larger response at 1.1 Hz than 1.4 Hz, because "noise" would have been included at more and lower frequency (noisier) harmonics at 1.1 Hz.
Thus, summing baseline-subtracted harmonic amplitudes is advantageous for a correspondence with interpretations of time domain brain responses. This approach has been used to quantify and compare overall response amplitude in a number of studies following Retter and Rossion (2016a), as mentioned previously (e.g., including time domain correspondences: De Keyser et al., 2018;Leleu et al., 2018; frequency domain analyses only: Damon et al., 2020;Dwyer et al., 2019;Gwinn & Jiang, 2019;Beck et al., 2018;Chemin et al., 2018;Guillaume et al., 2018;Gwinn et al., 2018;Xu et al., 2017). 7 If measures other than amplitude are desired (e.g., signal-to-noise ratio, z scores, or another statistic), the harmonic amplitudes can be extracted with an inclusion of a baseline frequency range (i.e., as a "chunk" of X Hz, centered around each frequency of interest), and then summed before these baseline-relative computations (Retter & Rossion, 2016a; see also Appendix 2 of Norcia et al., 2015; Box 2 of Rossion et al., 2020). In this way, a single statistical measure can be applied to the combined harmonic amplitude relative to its combined baseline amplitude (i.e., "noise"). Note that different approaches for combining harmonics may serve different ends, in that they describe different aspects of the signal, for example, the root-mean-square amplitude relates to the equivalent power of a flat (nonsinusoidal) signal; however, these aspects must be justified in relation to their physiological meaning.

What About Phase?
As addressed previously, there is a direct relationship between signal amplitude in the time and frequency domains. As a reminder, this relationship is given by Parseval's relation, which states that energy is conserved across the time domain (where energy equals the sum of the squared amplitudes) and frequency domain (where energy equals the sum of the squared root-mean-square amplitudes). In light of this, the amplitude across the harmonics relates to the overall amplitude of the signal in the time domain, regardless of phase.
However, to fully relate signals across the time and frequency domains, both the amplitude and phase of the representative frequency domain sine waves need to be taken into account. Without phase information, the fluctuation of amplitude across time (e.g., affecting local amplitude peaks) cannot be determined. Therefore, there is a cost toward relating time and frequency domain signals when excluding phase information. However, this cost is reasonably minor, for example, as relative phase changes across harmonics, it is possible that the latency of signal peaks varies, but that their amplitude does not ( Figure 6A and B). When relative phase does affect peak amplitudes, this influence is limited (e.g., compare Figure 6B and C). Moreover, despite relative phase changes, the area under the curve of the time domain signal may remain approximately constant ( Figure 6A-C; see also Heinrich, 2010). 8 Figure 5. The combination of harmonic amplitude. For A-C: Row 1: example brain responses, recorded with EEG (channel PO10 displayed here), were elicited from periodic visual stimulation of natural face (vs. object) images at various frequencies (data from Retter & Rossion, 2016a). Row 2: Frequency domain representations of these responses. The amplitude of harmonic responses above 0.1 μV are plotted in color, and these harmonic sine waves are superimposed in the corresponding color in Row 1. Row 3: The colored harmonic responses above are summed (set at Bin 0), following a baseline subtraction of "noise," defined as the average amplitude of the two adjacent frequency bins. Row 4: Similar response amplitudes are demonstrated in the time domain in Panels A and B, consistent with Row 3. Time outside one cycle duration is shadowed in gray, and the response amplitude range is emphasized between the horizontal red and blue lines. Moreover, in frequency tagging, it is worth remembering that the phase is not arbitrary: The phase of each relevant harmonic is determined relative to the time domain signal. In other words, the aligning positive and negative peaks of the sine waves across harmonic frequencies correspond to the time of the positive and negative peaks of the signal in the time domain (see again Figure 3B, for an example of deconstructive [when the signal is 0] and constructive [when the signal is 1], phase-locked harmonic sine wave superpositioning). This leads to phase differences across harmonics that are similar to describe time domain signals with similar temporal dynamics, despite the use of different stimulation frequencies ( Figure 6D; see also Strasburger, 1987). Thus, the influence of phase on combined harmonics is largely invariant of the stimulus presentation frequency, given consistent temporal dynamics of the response (as hinted at empirically, e.g., Retter & Rossion, 2016a;Appelbaum et al., 2006). Finally, it is worth noting a couple of helpful restrictions in the context of frequency tagging for combining harmonic responses: Only one sine wave is represented at each frequency bin, and nonharmonic frequencies are not considered, such that the response is fully periodic at the cycle duration of the fundamental frequency, F. 9

Which Harmonics to Consider?
Before combining harmonics, a decision of which harmonics to consider is required. To this extent, harmonics of interest (similarly to a ROI) must be defined. The first criterion for determining harmonics of interest is whether a higher harmonic is specific to its fundamental frequency. As mentioned previously, in frequency-tagging paradigms using a single-stimulation frequency, the higher harmonics are always specific to the fundamental frequency. However, in paradigms deploying multiple stimulation frequencies, unspecific harmonics may occur, which are often excluded from the analyses (for a paradigm-focused review, see Norcia et al., 2015). For example, two stimuli may be simultaneously presented at different spatial locations, one at 8 Hz (F 1 ) and the other at 6 Hz (F 2 ). If a response occurred at 24 Hz, it would not be specific to either stimulus, being the third harmonic of 8 Hz (3F 1 ) and the fourth harmonic of 6 Hz (4F 2 ), and would therefore be excluded from the analyses of responses to each stimulus (for further examples, see Table 1).
The use of a limiting frequency range is recommended here, either as determined a priori or from an assessment of the highest harmonic meeting a threshold (in terms of amplitude, signal-to-noise ratio, or significance). This is recommended because the upper frequency limit of harmonic responses, although affected by the overall strength of the signal, generally relates to the highest frequency that is strongly represented in the signal (see the Interpreting Harmonics section). The upper frequency limit of harmonic responses is thus often conserved across fundamental stimulation frequencies (see Figure 2 of Retter & Rossion, 2016a).
Note that the highest harmonic of interest can be determined either at the group level across conditions or as presented by any participant for any condition, but that typically a common range of frequencies of interest should be used across participants and conditions (e.g., Jacques et al., 2016). In this approach, there may be harmonic frequencies included for consideration at which there is no signal (e.g., in some participants, conditions, or regions of interest); however, including a small number of such frequencies is likely less detrimental (given that an appropriate baseline noise correction is applied, e.g., so that approximately zero amplitude values are added) than missing some frequencies containing a weak signal. Similarly, although responses are typically expected to occur consecutively across harmonic frequencies, in the event that a small number of within-range harmonic frequencies do not contain signal (above threshold), including them is typically tolerable (e.g., Rossion et al., 2020;Liu-Shuang et al., 2014).
Finally, in some cases, harmonic responses appear to be qualitatively different from one another. This may occasionally be related to physiological sources: for example, different harmonic response patterns are generated from the recordings of frequency-tagged responses from singleversus double-opponent cortical cells (Movshon et al., 1978). However, more often, physiological sources may only be tentatively inferred, for example, when different EEG scalp topographies are observed at different harmonic frequency ranges (e.g., Rossion, 2014; see the What Do Higher Harmonic Responses Represent? section). In this case, is it appropriate to select subranges of qualitatively homogeneous harmonics to consider and/or combine? Perhaps, although it should be remembered that harmonic responses are not independent of one another (e.g., Retter & Rossion, 2016a;Zhou et al., 2016) and therefore should also be described individually and/or summed all together (see the Should Higher Harmonic Responses Be Combined, and If So, How? (Reprise) section).
It is not advised to select or subgroup harmonics a priori in accordance to only their number, unless this is explicitly derived from the stimulation paradigm (see Table 1). For example, there is a persistent history of considering the first versus second harmonic response (Saupe, Schröger, Andersen, & Müller, 2009;Pastor, Valencia, Artieda, Alegre, & Masdeu, 2007;Kremers & Scholl, 2001;Falsini et al., 1999;Burns, Elsner, & Kreitz, 1992;Baker & Hess, 1984; the first vs. second harmonic, rather than odd vs. even harmonics: Kim, Grabowecky, Paller, & Suzuki, 2011;Kim, Grabowecky, Paller, Muthu, & Suzuki, 2007). This relates to early interpretations of the first harmonic reflecting asymmetries in responses following on and off stimulation cycles (e.g., Clynes, Kohn, & Lifshitz, 1964) and the second harmonic being typically dominant with pattern reversal stimulation (see the following section). However, the presence of third, fourth, and further higher harmonics-and their interdependence-is indicative of the limits of such an oversimplification in most cases.

INTERPRETING HARMONICS
Why Are There Higher Harmonics?

Nonsinusoidal Brain Responses
Higher harmonic responses represent complex neural responses in the time domain. At a fundamental level, these harmonic responses are like any other frequency domain representations: They are sine waves described by frequency, amplitude, and phase ( Figure 2). Although only one sine wave is required to describe a sinusoidal signal in the time domain, a combination of (a lot of ) sine waves is required to describe complex signals in the time domain (Figure 3). Frequency-tagged brain responses are periodic in the time domain, and thus, only sine waves periodic to their fundamental frequency (i.e., the harmonics) are mathematically available to describe them. Simple brain responses require few harmonics, whereas

Single frequency A-A-A-A-A-A-A-…
F --F and its harmonics Adrian and Matthews (1934), Regan (1989)

Multiple frequencies (F 1 : A-A-A-A-A-A…)
F 1 (e.g., 8 Hz) F 2 (e.g., 6 Hz) Any coinciding harmonics of F 1 and F 2 (e.g., 24 Hz) Exclude overlapping harmonics from analysis of both F 1 and F 2 Regan and Heron (1969), Regan (1989) ( F harmonics coinciding with 2F and its harmonics Exclude 2F (even) harmonics from the analysis of the F (odd) harmonics Tyler and Kaitz (1977), Victor and Zemon (1985) Oddball In the paradigm example sequence illustrations: A = one stimulus or stimulus type; B = another stimulus or stimulus type. Special cases: In the case that A and B stimuli in a symmetry/asymmetry paradigm lead to symmetrical brain responses (e.g., if representing pattern reversals), only even harmonics are observed (Norcia et al., 2015;Hou et al., 2003;Cobb et al., 1967); in a combined symmetry/asymmetry and oddball design (Braddick et al., 2005), the odd harmonic analysis is unaffected. In the case that multiple frequencies lead to intermodulation (i.e., additive and subtractive interaction frequencies and their harmonics), the analysis of the intermodulation harmonics should exclude the overlapping harmonics of F 1 and F 2 (e.g., Gordon et al., 2019;Boremanse et al., 2013;Applebaum et al., 2009;Hou et al., 2003;Burns et al., 1992;Zemon & Ratliff, 1984). In the case that a stepwise sweep design is applied to a symmetry/asymmetry paradigm, this does not affect the harmonic analysis (see Norcia et al., 2015).
complex responses require more harmonics (Figure 4). For example, lower stimulation frequency responses often have more harmonics, because there are relatively more harmonic frequencies available within a relevant frequency range ceiling ( Figure 5).
Limitations of a nonlinearity account. Higher harmonic responses have often been interpreted as being caused by nonlinearities in the stimulus presentation and/or brain responses (e.g., Gordon, Hohwy, Davidson, van Boxtel, & Tsuchiya, 2019;Norcia et al., 2015;Burns et al., 1992;Regan, 1966Regan, , 1989Troelstra, 1971;Spekreijse, 1969;Van der Tweel & Verduyn Lunel, 1964, 1965; see also Shapley, 2009). This relates to early attempts to present stimuli perfectly sinusoidally (e.g., a uniform field being modulated sinusoidally in luminance; since van der Tweel et al., 1958). The rationale was that if the brain's response to a sinusoidal stimulus was linear at the level of recording, it too would be perfectly sinusoidal in following this stimulus and would be represented in the frequency domain by a response only at the fundamental frequency, that is, without higher harmonics (see the Frequency Domain Representations section). Contrarily, higher harmonic brain responses were produced in most cases and were attributed to nonlinearities in the brain's responses themselves (e.g., nonlinear action potential firing, neural population response dynamics; Shapley, 2009;Skottun et al., 1991;Movshon et al., 1978). 10 However, more recent studies suggest that the amount of nonlinearity, or complex temporal frequency content, in stimulus presentation may not correspond with the amount of nonlinearity in the brain's response at the population level: There was little to no difference in the higher harmonic amplitude distributions in response to (imperfect) sinusoidal versus squarewave (i.e., abrupt on/off ) complex stimulus presentation (Dzhelyova, Jacques, & Rossion, 2017;Retter, 2016;Fawcett, Barnes, Hillebrand, & Singh, 2004;Burns et al., 1992). Moreover, although the inherent nonlinearity of the brain's responses could account for higher harmonics, in practice, the amplitude of the higher harmonics is not always above noise level (as addressed in the Introduction) or is very low, suggesting only a modest contribution of this factor (see the Frequency-tagged Responses In The Frequency Domain section). That being said, one source of complexity in the brain's responses (even those underlying the first harmonic) may be these nonlinearities.

What Do Higher Harmonic Responses Represent?
Higher harmonic responses represent the relevant frequency characteristics of the response in the time domain (e.g., Rossion et al., 2020;Zemon & Gordon, 2018;Retter & Rossion, 2016a;Zhou et al., 2016;Gaume et al., 2014;Heinrich, 2010;Norcia, Sato, Shinn, & Mertus, 1986;Baker & Hess, 1984;Galambos et al., 1981). That is, dynamics of the time domain response best represented at different frequency ranges will produce more amplitude in those frequency ranges in the frequency domain ( Figure 7; compare with the frequency domain representation in Figure 3C; for another example of frequency representations over time, see Makeig et al., 2002). To visualize the impact of (a range of ) individual harmonics, partial harmonic time domain reconstructions have been plotted (Sieving et al., 1998;Baker & Hess, 1984). Note that the amplitude distribution of harmonics across frequencies is, however, affected by the fundamental stimulation frequency and the overall amplitude of the signal. Moreover, individual harmonics do not represent independent or temporally separated aspects of a time domain response (e.g., Retter & Rossion, 2016a;Zhou et al., 2016).
In line with the above interpretation, harmonic responses may be (gradually) influenced, quantitatively and/or qualitatively, by the frequency at which they fall. For example, Retter and Rossion (2016a, Figure 2C) described harmonic EEG responses that were gradually characterized by frequency, in terms of amplitude and scalp lateralization, commonly across conditions with different, Figure 7. A time domain signal can be fit with segments of sine waves of different frequencies.
Although this is not analogous to a frequency transformation, it hints at the range of frequencies that may be optimal for representing this signal (over time).
low stimulation frequencies. Harmonic responses may also appear to group into somewhat distinct frequency ranges, for example, visually evoked EEG responses above about 10 Hz having a more medial-occipital ("low-level") scalp topography Rossion, 2014; see also Zemon & Gordon, 2018). In general, higher harmonics may be more associated with earlier response dynamics, as the onset slopes of (ERP) responses are typically steeper than the offsets (Norcia et al., 1986;as in Figure 7). In some cases, different experimental effects may be pronounced at some harmonic frequencies and not others (e.g., an increased response because of attention: Saupe et al., 2009;Pei, Pettet, & Norcia, 2002; opposing effects following current stimulation: Ruhnau et al., 2016). Indeed, limited frequency ranges may be appropriate for measuring certain properties (e.g., a lower frequency range for chromatic than luminance signals; see Parry et al., 2012;Burns et al., 1992). Thus, in addition to combining harmonics, investigating individual (ranges of ) harmonics may provide insight into the functional dynamics of the neural processes that occur at their respective frequencies.
It is not recommended to display harmonics independently by their sequential number (e.g., across different fundamental stimulation frequencies: Herrmann, 2001;Ross et al., 2000;Troelstra, 1971). Indeed, because harmonic responses are characterized by their frequency, they are generally not well characterized by sequential number, irrespective of frequency, unless differentially tagged in the stimulation paradigm (see the Which Harmonics to Consider? section). However, one case in which harmonic frequency does not characterize responses well is at high stimulus presentation frequencies, for which the first harmonic is dominant. At high frequencies, overlapping responses interfere with one another, decreasing the measured response's complexity and amplitude and therefore harmonic contents (see the Frequency-tagged Responses in the Frequency Domain section). Thus, a dominant exponential decrease of harmonic amplitude may be observed across sequential harmonics (compare Figure 4B and C; Retter et al., 2020,  Figure S2). However, this must not be taken as evidence of distinct brain responses being evidenced at the first versus second (vs. third, etc.) harmonic (e.g., see the discussion of Saupe et al., 2009).
Finally, it may be noted that consecutive response interference is not only dependent upon the duration of the brain response being measured but is also influenced by the spatiotemporal dynamics of the interfering response components. In models of consecutive responses, destructive superpositioning of overlapping (ERP) response components may reduce or even eliminate (selective harmonic) responses (constructive superpositioning is also possible; see Heinrich, 2010, Figures 1 and 2; see also Footnote 4). In practice, overlapping brain responses do not occur independently of one another (see Retter et al., 2020;Keysers & Perrett, 2002). However, particularly when distinct response sources are implicated (e.g., light-and dark-preferring neural populations or across sensory modalities), differential interactions of stimulation frequency and selective harmonic amplitudes could also provide insights into the underlying dynamics of the brain responses.

Should Higher Harmonic Responses be Combined, and If So, How? (Reprise)
Limiting the analysis of frequency-tagged responses to a (nonpredominant) fundamental harmonic is not recommended: Taking harmonic responses into account leads to substantially improved response quantification and classification (Figure 8; Cetin et al., 2020;Zemon & Gordon, 2018;Retter & Rossion, 2016a;Chen et al., 2015;Tlumak et al., 2011;Cebulla et al., 2006;Muller-Putz et al., 2005;Davila et al., 1998). The combination of harmonic amplitudes through summation is justified through the principles of frequency-based analyses, and it leads to a combined response measurement that relates to typical time domain amplitude measurements. It is useful for comparing brain response amplitudes overall, especially those with different temporal dynamics or following different stimulus presentation rates. Note that it is does not preclude, but rather is complementary, to the description of harmonic responses individually.

CONCLUSION
Stimuli that are presented periodically generate periodic responses of the brain that are often complex (i.e., nonsinusoidal). To capture and describe these complex brain responses overall, (baseline-corrected) frequency-tagged harmonic response amplitude can be combined through simple summation.

Diversity in Citation Practices
A retrospective analysis of the citations in every article published in this journal from 2010 to 2020 has revealed a persistent pattern of gender imbalance: Although the proportions of authorship teams (categorized by estimated gender identification of first author/last author) publishing in the Journal of Cognitive Neuroscience ( JoCN) during this period were M(an)/M = .408, W(oman)/M = .335, M/ W = .108, and W/ W = .149, the comparable proportions for the articles that these authorship teams cited were M/M = .579, W/M = .243, M/ W = .102, and W/ W = .076 (Fulvio et al.,JoCN,33:1,. Consequently, JoCN encourages all authors to consider gender balance explicitly when selecting which articles to cite and gives them the opportunity to report their article's gender citation balance. Notes 1. A note on nomenclature: Here, the "first" harmonic is the fundamental stimulation frequency. In another existent convention, the "first" harmonic is the double of the fundamental stimulation frequency. 2. Note that in frequency domain transformations, technically, a complex-valued combination of sine waves and cosine waves is used to represent the signal. In many discrete, fast Fourier transforms, the sine component carries only the amplitude information, and the cosine component carries only the phase information. Here, "sine waves" are referred to as complex entities themselves, combining both amplitude and phase information. Additionally, note that phase spectra are rarely plotted in the same format as amplitude spectra, because (1) their values are circular (e.g., with equivalent distance from 0°to 359°and from 0°to 1°; see an alternative plotting example in Figure 6D) and (2) the phase at non-frequency-tagged signal frequencies is random (i.e., full-range noise). 3. In the time domain, energy is equal to power over time, which is equal to the sum of the squared amplitudes. In the frequency domain, energy is equal to the sum of the squared root-meansquare amplitudes (root-mean-square amplitude of a sine wave = amplitude/√2; see Smith, 1997). 4. Although largely beyond the scope here, note that it is extensively debated whether frequency-tagged EEG responses reflect the (linear) superposition of ERPs or whether they reflect an interaction with endogenous oscillations in the brain (see Zoefel, Oever, & Sack, 2018;Retter & Rossion, 2016a;Heinrich et al., 2015;Norcia et al., 2015;Keitel, Quigley, & Ruhnau, 2014;Gruss, Wieser, Schweinberger, & Keil, 2012;Capilla et al., 2011;Heinrich, 2010;Makeig et al., 2002;Herrmann, 2001;Galambos, Makeig & Talmachoff, 1981;Donker, 1975). 5. Note that, further, the signal and noise may interact nonlinearly (supra-additively), but when the signal is much (e.g., 3 or 4 times) greater than noise, this modest contribution is negligible (Norcia, Tyler, Hamer, & Wesemann, 1989;Peli, McCormack, & Sokol, 1988;Strasburger, 1987;see Bach & Meigen, 1999). Particularly for weaker signals, a more conservative baseline correction could be computed as the square root of the signalpower-minus-noise-power (see Appendix 1 of Norcia et al., 1989). 6. Please note an exception: Sometimes, the first adjacent bin, on either side of the frequency bin of interest, is excluded, so as to avoid overspill, that is, signal leakage, in the noise estimate at high-frequency resolutions (e.g., Rossion, Alonso-Prieto, Boremanse, Kuefner, & Van Belle, 2012). 7. See also publications motivated by the as yet unpublished findings of Retter and Rossion (2016a) from the Face Categorization Lab (https://face-categorization-lab.webnode.com/; Jacques, Retter & Rossion, 2016;Jonas et al., 2016;Liu-Shuang, Torfts & Rossion, 2016;Lochy, Van Reybroeck, & Rossion, 2016;Retter & Rossion, 2016b;Dzhelyova & Rossion, 2014a, 2014b. 8. In any case, it remains unclear at present how phase could be combined meaningfully across harmonics, beyond time domain latency, to relate to functional, physiological processes (see Strasburger, 1987). 9. In frequency tagging with a single-stimulation frequency, there are not multiple sine waves of the same frequency that interfere with each other, for example, sine waves with a 180°phase difference that produce complete interference (zero amplitude sum). Furthermore, when a complete F cycle is considered, the amplitude is representative of the complete signal (which is not the case for multiple [nonharmonic] frequencies that may sum differentially across short time segments of evaluation). Note that these conditions are not always met in other contexts in which frequency-based analyses are performed. 10. Higher (intermodulation) harmonic responses have then been used to characterize the nonlinearities of brain responses, leading to the exclusion of some potential nonlinear models in specific cases (e.g., Baker & Wade, 2017;Hou et al., 2003;Victor & Conte, 2000;Regan & Regan, 1988b).