Towards an objective evaluation of EEG/MEG source estimation methods – The linear approach

The spatial resolution of EEG/MEG source estimates, often described in terms of source leakage in the context of the inverse problem, poses constraints on the inferences that can be drawn from EEG/MEG source estimation results. Software packages for EEG/MEG data analysis offer a large choice of source estimation methods but few tools to experimental researchers for methods evaluation and comparison. Here, we describe a framework and tools for objective and intuitive resolution analysis of EEG/MEG source estimation based on linear systems analysis, and apply those to the most widely used distributed source estimation methods such as L2-minimum-norm estimation (L2-MNE) and linearly constrained minimum variance (LCMV) beamformers. Within this framework it is possible to define resolution metrics that define meaningful aspects of source estimation results (such as localization accuracy in terms of peak localization error, PLE, and spatial extent in terms of spatial deviation, SD) that are relevant to the task at hand and can easily be visualized. At the core of this framework is the resolution matrix, which describes the potential leakage from and into point sources (point-spread and cross-talk functions, or PSFs and CTFs, respectively). Importantly, for linear methods these functions allow generalizations to multiple sources or complex source distributions. This paper provides a tutorial-style introduction into linear EEG/MEG source estimation and resolution analysis aimed at experimental (rather than methods-oriented) researchers. We used this framework to demonstrate how L2-MNE-type as well as LCMV beamforming methods can be evaluated in practice using software tools that have only recently become available for routine use. Our novel methods comparison includes PLE and SD for a larger number of methods than in similar previous studies, such as unweighted, depth-weighted and normalized L2-MNE methods (including dSPM, sLORETA, eLORETA) and two LCMV beamformers. The results demonstrate that some methods can achieve low and even zero PLE for PSFs. However, their SD as well as both PLE and SD for CTFs are far less optimal for all methods, in particular for deep cortical areas. We hope that our paper will encourage EEG/MEG researchers to apply this approach to their own tasks at hand.


Introduction
Electro-and magnetoencephalography (EEG and MEG) record electrical brain activity with high temporal resolution (in the millisecond range), i.e., at a time scale required to track fast perceptual and cognitive processes in real-time. However, it is well known that this strength of EEG and MEG is accompanied by limitations with respect to their spatial resolution Hämäläinen et al., 1993 ;Hauk et al., 2011 ). Unfortunately, exactly what these limitations are and especially what they mean for experimental researchers is still unclear. The question "What is the spatial resolution of EEG/MEG? " can only be answered with many ifs and buts, because the answer depends on a large number of parameters such as the measurement configuration (e.g., com-researchers without a methodological background to establish which methods are most suitable and which parameter settings most optimal for their task at hand. Commercial and freeware software packages for EEG/MEG analysis offer a large selection of source estimation methods, but rarely tools to objectively evaluate them Dalal et al., 2011 ;Delorme and Makeig, 2004 ;Gramfort et al., 2013 ;Litvak et al., 2011 ;Oostenveld et al., 2011 ;Tadel et al., 2011 ). Importantly, it can be hard for researchers to decide which aspects of a methodological study are relevant to their particular research question, e.g., whether assumptions about the source scenarios (such as single, multiple and distributed sources) are valid for their particular dataset. A method that works well for a few focal sources and early evoked potentials may fail for resting state activity or complex distributed sources underlying cognitive processes. Ideally, researchers should have at least an intuitive understanding of the EEG/MEG inverse problem and the distinguishing characteristics of different source estimation approaches, as well as the software tools to evaluate their performance on their own datasets. Here, in the tutorial part of this paper aimed at experimental researchers, we will describe a conceptual framework that intuitively describes the assumptions behind the most common distributed source estimation methods, namely L2-minimum-norm estimation (L2-MNE) and linearly-constrained minimum variance (LCMV) beamforming. We will use core elements of this framework, namely the resolution matrix and the associated point-spread and cross-talk functions (PSFs and CTFs) ( Hauk et al., 2011 ;Liu et al., 2002 ;Menke, 1989 ), to derive resolution metrics that objectively describe aspects of source estimates that researchers need to interpret their results (e.g., localization error and spatial extent of source leakage). We will illustrate this approach using openly available datasets, confirm results from previous methods comparisons ( Hauk et al., 2011 and extend them to a larger range of source estimation methods including weighted and normalized L2-MNE and two LCMV beamformers 1 . The reason why spatial resolution is of particular concern for EEG and MEG is the "ill-posed " inverse problem: Even if it were possible to perfectly record full three-dimensional electric and magnetic field distributions around the head, it would still be impossible to unambiguously determine the three-dimensional source current distribution inside the brain that gave rise to those measurements ( Sarvas, 1987 ;von Helmholtz, 1853 ). Several methods have been developed to find possible solutions to this inverse problem, using different constraints Friston et al., 2008 ;Fuchs et al., 1999 ;Michel et al., 2004 ;Wipf and Nagarajan, 2009 ). In cases where it is certain that the signal is generated by one dominant focal source the location, orientation and amplitude of the source current can be estimated with reasonable accuracy using a single equivalent dipole fit ( Lütkenhöner, 1998b ;Scherg and Berg, 1991 ). Unfortunately, in many EEG/MEG experiments, especially those with low SNR or involving complex brain processes, it is impossible to specify prior assumptions on source activity in a precise mathematical manner. In those cases, it is common to employ distributed source models, which attempt to explain the measured signal as source-current distributions within the whole brain volume or across the cortical surface. It is generally accepted that linear distributed-source solutions do not accurately reconstruct the true source distribution, but rather produce a "blurred " estimate of it. The degree of this blurring, i.e., the spatial resolution of these methods, and how it is affected by parameters such as sensor configuration and source depth, is still largely unknown. The evaluation of spatial resolution is more complicated for EEG/MEG than it is for functional magnetic resonance imaging (fMRI), as it depends on more parameters such as those described in the first paragraph above. Researchers therefore need tools for a comprehensive evaluation of the spatial resolution of their experimental setup at hand.
The most popular methods for EEG/MEG distributed source estimation are arguably L2-minimum-norm type methods (e.g., classical, optionally weighted, L2-minimum-norm estimation, L2-MNE; dynamic statistical parametric mapping, dSPM; standardized and exact low resolution electromagnetic tomography, s/eLORETA) and linearly constrained minimum variance (LCMV) beamformers. These methods all have in common that they are based on linear transformations of the data and can therefore be described and evaluated in the framework presented in this paper (with some caveats discussed later).
Linear transformations have some convenient properties that can be exploited for the evaluation of their spatial resolution ( Backus and Gilbert, 1968 ;Dale and Sereno, 1993 ;Hauk et al., 2019 ;Menke, 1989 ). Most importantly, the superposition principle holds, i.e., the source estimate of a signal produced by multiple sources is the sum of the source estimates produced by each individual source. Thus, knowing the source estimates for all possible point sources (i.e., their point-spread functions, PSFs) enables us, in principle, to construct the source estimate for any arbitrary configuration of sources. This framework is therefore well-suited to provide an objective and generalizable evaluation of spatial resolution for linear methods. Some important caveats with respect to the applicability of the superposition principle to (non-linear) intensity distributions and beamformers will be discussed below (see also Hauk, Stenroos et al., 2019 ).
In addition to point spread, another key concept for the evaluation of source estimation methods is that of cross-talk: How much is an estimate of activity for one source at a particular location affected by activity from other locations? This is described by cross-talk functions (CTFs), which for a particular source describe how all other possible sources with unit strength would affect the source estimate for this source of interest. This concept is especially relevant for the application of spatial filters (or "software lenses ") that are supposed to filter out activity from a particular brain location while suppressing activity from other locations ( Grave de Peralta Menendez et al., 1997 ;Hauk et al., 2019 ;Van Veen et al., 1997 ). While the concepts of PSFs and CTFs are related, there are important differences that will be discussed in a separate section below.
For practical purposes, linear methods can be described by means of linear algebra, i.e., using vector and matrix operations. This is a relatively intuitive approach to data analysis because it offers geometrical interpretations of mathematical concepts and is straightforward to use with programming languages such as Matlab and Python. However, we cannot assume that researchers from a non-methods background are proficient in linear algebra ( Hauk, 2020 ). While a tutorial in the relevant linear algebra is beyond the scope of this paper, we attempted to describe the basics of linear EEG/MEG source estimation in plain English as well as at the simplest adequate level using linear algebra. We also provide a brief primer on the relevant linear algebra in Appendix C. More information can of course be found in online resources or books (e.g., Cohen, 2021 ).
As we intend to keep matters as simple as possible, we open ourselves up to criticism that we simplify too much. Reality is always more complex than our models or approximations, but the question is whether those are still useful. Obvious concerns in the case of EEG/MEG source estimation are our approximations for head modeling (e.g. Stenroos et al., 2014 ), the way we deal with noise (e.g. Samuelsson et al., 2021 ), and the different theoretical frameworks that can be used to approach the EEG/MEG inverse problem (e.g. Henson et al., 2011 ;Sekihara and Nagarajan, 2015 ;Tarantola, 1994 ). In our view, the linear framework described in our paper provides a basis to understand and deal with more complex or more specialized issues of source estimation, also in different theoretical frameworks.
For an objective methods comparison, it is necessary to define resolution metrics that summarize important features of these distributions, which can then be visualized as whole-brain distributions, and compared across methods and across the parameter space. Researchers commonly interpret specific features of their source estimates, such as their peaks or their distribution within and across different brain regions. Thus, we need metrics that are location-specific and reflect meaningful features of the distributions that are relevant for their interpretation. Two aspects of spatial resolution are particularly important for the most common applications of EEG/MEG source estimation. First, the peaks of PSFs and CTFs should be close to the true source location. This has often been evaluated by means of the peak localization error, i.e., the Euclidean distance between the estimated peak and the true source location. Second, to be able to evaluate to which degree two distributions overlap with each other, we need a measure of spatial extent. This has for example been quantified using the standard deviation around the peak of the distribution, also called spatial deviation. Some of these metrics have been used previously to evaluate the benefit of combining EEG and MEG measurements ( Molins et al., 2008 ), to compare noisenormalized MNE methods ( Hauk et al., 2011 ), to compare different head models ( Stenroos and Hauk, 2013 ), and to compare different MEG sensor types ( Iivanainen et al., 2017 ).
In the present study, we will provide a tutorial to linear source estimation and demonstrate how the corresponding concepts can be used in practice to obtain intuitive as well as objective estimates of spatial resolution for different source estimation methods. We will compute resolution metrics for PSFs and CTFs using a combined EEG and MEG measurement configuration to compare some of the most popular source estimation methods in a best-case scenario. We extend the results of previous studies by including a larger number of methods in our comparison, by computing resolution metrics not just for PSFs but also for CTFs, and by providing a more detailed analysis of the relationship between spatial resolution and source depth. We would like to highlight that this approach to resolution analysis can easily be applied to any new data set, e.g., with different sensor configurations, head models and noise levels. It could therefore become a standard tool to guide the interpretation of EEG/MEG source estimation results.

Methods
In this section, we will describe the basics of the EEG/MEG inverse problem, linear systems analysis and common source estimation methods in basic matrix notation. A brief primer on some of the relevant linear algebra is provided in Appendix C.

Linearity and the superposition principle
Many procedures in signal or image processing can be described as linear, at least in reasonable approximation, including Hifi systems, microscopes, telescopes or EEG/MEG source estimation. Linear systems, by definition, have the property that the output of any sum of inputs is just the sum of the outputs for the individual inputs ( "superposition principle "). This means that if we can approximate a complex input (e.g., multiple sources or a distributed source) as a sum of point sources, then we can also approximate the output by a corresponding sum of the outputs of the individual point sources, or "point-spread functions " (PSFs). This is illustrated for the cases of a hypothetical microscope ( https://en.wikipedia.org/wiki/Point_spread_function ), as well as for EEG/MEG sensor and source space data in Fig. 1 . Fig. 1 A shows example point-spread functions for a microscope, i.e., the degree of blurring for two circular objects. The ideal lens is supposed to provide an accurate image of the underlying objects, magnified by a certain factor. In reality, the glass may contain impurities, the grinding of the lens may not have been perfect, and there are also fundamental physical limits with respect to light diffraction depending on the light's wavelength. Thus, a point source and the sharp borders of the discs get blurred. The degree of blurring may depend on the location of the object relative to the lens, e.g., whether it is at the center or at the border, or whether the lens contains impurities in some parts but not others (not shown in this Figure).
The same principle applies to EEG/MEG topographies in signal space ( Fig. 1 B). These examples illustrate how the "blurring " of signals from point sources can result in partial overlap between PSFs, which can make them hard or impossible to separate when the two sources are present simultaneously. This is an important aspect of spatial resolution that can be described by metrics such as localization accuracy and spatial extent (which will be discussed in more detail later). Fig. 1 C demonstrates how the superposition principle applies to linear source estimation (here L2minimum-norm estimates; further examples for MNE-type methods and beamformers are presented in Figs. 3 and 6 , respectively).
More formally, a transformation T is linear if for any arbitrary input signals A and B and scalar values a and b the following is true: In the case of EEG/MEG source estimation, A and B can be signal patterns of dipolar current sources with unit strength, a and b source strengths, and the transformation T the "inverse operator " (usually a matrix) that turns the measured EEG/MEG data into current distributions in the brain. This means that we can characterize the properties of the transformation on the basis of the transformations of individual point sources with unit strength, i.e., their PSFs.

Resolution matrix
If we imagine our EEG/MEG measurements followed by source estimation as a "lens " with which we are looking through the skull into the brain, then the resolution matrix describes the blurring or distortion caused by this lens. In the case of EEG/MEG, the resolution matrix describes how the estimate for the sources is blurred and distorted by the combination of the measurement itself and the inverse operator applied to the data. Thus, it depends on the physical characteristics of the inverse problem, but also on parameter choices with respect to head model and inverse operator.
Before we can start thinking about estimating neural sources from our EEG/MEG data, we have to parametrize these sources such that we can include them in an explicit model (e.g., the general linear model in Eq. (2 ) below) and describe any additional constraints on them mathematically. Thus, we have to define a source space. The dominant sources of EEG/MEG signals are commonly assumed to be located in grey matter along the cortical surface ( Hämäläinen et al., 1993 ). In order to make the problem numerically tractable, this surface is usually discretized into a larger number ( ∼> 10,000) of vertices. This number is much larger than the number of measurement sensors. Each vertex is represented by an equivalent current dipole with a particular strength and orientation of current flow. The orientation of dipoles at vertices can be fixed a priori to be perpendicular to the cortical surface, in which case each vertex is associated with one value for the corresponding dipole strength in the model. If the orientation of the dipole is unknown, each vertex can be represented by 2 (MEG only) or 3 (EEG or EEG + MEG) perpendicular dipoles, in which case each vertex is associated with 2 or 3 values in the model, respectively. Because these dipoles reflect the activity of a vertex condensed into an infinitesimally small point, they can be considered as point sources as discussed in the previous section. In the following, the values of dipole strengths to be estimated across all vertices will be represented in a vector s ( n x1, n : number of vertices, or dipole strengths to be estimated).
The resolution matrix is an established tool for the characterization of spatial resolution for linear estimators ( Backus and Gilbert, 1970 ;Dale and Sereno, 1993 ;Grave de Peralta Menendez et al., 1997 ;Hauk, 2004 ;Krishnaswamy et al., 2017 ;Liu et al., 2002 ;Menke, 1989 ). More formally, if s is the true source vector and ̂ its estimate (in both cases each element represents a dipole in source space), then the resolution matrix describes the relationship between and ̂ in the form ̂ = . More details will be provided below, but the important message here is that R specifies what stands between our estimate and the truth. Thus, while we will not be able to know the true sources , the resolu- For two real objects (top left), the output image (right) is the real image convolved (or "blurred ") with the PSF. From https://en.wikipedia.org/wiki/Point-spread_function . (B) Illustration of the superposition principle for the EEG forward problem.A point source (current "dipole ") produces a characteristic voltage distribution, or topography, on the scalp (top left). This topography is represented by a column vector . The topographies of two separate sources add up at the scalp. If they are sufficiently distant from each other, their topographies hardly overlap (top right). As the sources move closer together, their topographies begin to overlap (bottom left). For close-by sources with different depths, the topography of the superficial source may overshadow the topography of the deep one (bottom right). The Gaussian-shaped EEG topographies for radial dipoles in this example were chosen for simplicity, and are only a very coarse approximation of reality. (C) The superposition principle for linear distributed source estimates. The source distributions for two point sources (PSFs) are shown in panels 1 and 2, respectively (L2-MNE, MEG only). The locations of the point sources are indicated by two black balls on the inflated cortical surface. The source distribution for both sources together is just the sum of their individual PSFs (panel 3). Red and blue colors indicate out-and in-going source currents with respect to the cortical surface, respectively. tion matrix can help us estimate how far we are away from the truth. More specifically, it can tell us for which sources we may at least get close, and for which sources we may not stand a chance.
The logic behind the construction of the resolution matrix is as follows: We know that our measured data are the result of the forward solution for the source distribution, and our linear source estimate is the result of the linear estimator applied to our measured data. If we combine the forward solution and the inverse estimator together into one transformation, we directly transform the unknown true source distribution into its estimate. Ideally, this transformation should be the identity, i.e., our estimate should be exactly the true source distribution. We know that we cannot achieve this for realistic models due to fundamental physical and mathematical limitations of the inverse problem, but an inspection of this transformation can tell us how close we are. This transformation is given by the resolution matrix, which will now be derived more formally.
The relationship between the source distribution ( n x1, n : number of point sources) and the data in n channels ( m x1, m : number of measurement channels) can be written in matrix form where ( m x n ) (is the so-called leadfield matrix which contains information about head geometry (shape and conductivities), measurement configuration (sensor types and position), and the physics of signal generation (quasi-static approximation of Maxwell's equation Sarvas, 1987 ), and ( m x 1 ) is the noise term. This forward problem is linear (illustrated in Fig. 1 B). Source estimation methods do not have to be linear, but we can try to find an optimal linear matrix ( n x m ) (as it is sometimes referred to as an inverse k ernel; optimal in a way to be defined), such that we obtain an estimate ̂ for the source distribution ̂ = (3) We can predict the signal produced by this estimate from Eq. (2) , i.e., we combine Eqs. (2) and (3) : The resolution matrix ( n x n ) provides the desired relationship between true and estimated sources, except for the last term that reflects how noise is transformed from signal space into source space. This will be relevant for regularization procedures later.

Point -spread and cross -talk functions (PSFs and CTFs)
The resolution matrix is still a rather abstract concept. How can we get some meaningful information out of it that can help us in practice? As illustrated in Fig. 1 , linear transformations can be meaningfully evaluated for point sources. PSFs and CTFs address two aspects of linear source estimation methods ( Fig. 2 shows examples for L2-MNE, Grave de Peralta Menendez et al., 1997 ;Hauk et al., 2011 ;Liu et al., 2002 ): (1) How does the method spread the activity of a point source with unit strength to other locations (PSF), and (2) how is the result of the method for one particular source affected by other sources (CTF)? In order to derive PSFs and CTFs more formally, it is useful to break down Eq. (3) for individual elements of and ̂ , i.e., corresponding to individual point sources.
PSFs are more widely used than CTFs in methodological comparisons of source estimation methods, but as we will argue below both are important to fully characterize resolution for a particular source. As the name suggests, PSFs represent how activity of a point source (such as one dipole at a particular location with a fixed orientation and unit activation) is spread across other sources in the source space by the source estimation procedure. More specifically, in order to obtain the PSF for point source i , , we first compute the topography in sensor space produced by source i with unit strength, and then subject this topography to our inverse operator. This topography is the i-th column of the leadfield matrix, i.e., . i . Therefore: In other words, the PSF for source i is the i th column of the resolution matrix R . This also shows that is a weighted sum of the columns of the inverse matrix , with the weights given by . i .
While PSFs show how point sources can potentially affect the estimates for other sources, we usually also face another (though related) question: If we are interested in activity in one particular location i , then how does activity from other locations affect the estimate at the location of interest? In other words, we are interested in the resolution of the specific estimator for location i , which is . (i.e., the i-th row of the inversion matrix ). The CTF for source i , , is then obtained by applying the linear estimator for source i to the topographies of all possible point sources in source space with unit strength. These topographies are the columns of the leadfield matrix . Thus In other words, the CTF for source i is the i-th row of the resolution matrix. It also shows that is the weighted sum of the rows of the leadfield matrix, with the weights given by i . .
In summary, here are the formulae for PSFs and CTFs next to each other: As stated above, PSFs are the columns of the resolution matrix R and linear combination of columns of the inverse matrix , and CTFs are the rows of the resolution matrix R and linear combinations of the rows of the leadfield matrix L . In general, the PSF and CTF for a source i will have different shapes, which is important for the comparison of different source estimation methods. However, the resolution matrix for the unweighted L2-minimum-norm estimate is symmetric (see below) and for this method PSFs and CTFs can be represented by the same distributions.
While the leadfield matrix L is determined by our forward model, the inversion matrix K can in principle be chosen arbitrarily to optimize our PSFs and CTFs. But whatever we do, linear algebra poses some fundamental limits ( Menke, 1989 ): (1) A distribution that cannot be approximated by any weighted combination of the rows of the leadfield matrix cannot be a CTF.
(2) The number of columns of the inverse estimator is the number of recording channels m , and a PSF is a linear combination of these columns. We can therefore get at most m well-shaped PSFs, but this does not guarantee that the remaining n -m PSFs are well-shaped (or even close) as well.
(3) Because appears in both formulae, we cannot optimize PSFs and CTFs independently of each other.
This implies that every optimization of PSFs and CTFs will require a trade-off between different resolution criteria. We will address this issue using quantitative resolution metrics below.

Noise and regularization
Source estimation attempts to determine the brain sources that generated the measured signals. It is highly desirable to suppress the influence of noise as much as possible, i.e., to minimize the term in Eq. (3) . This is addressed by regularization of the source estimates. In short, we accept some loss of spatial resolution and/or goodness-of-fit in return for a solution that is less affected by noise, which usually results in smoother source distributions ( Backus and Gilbert, 1968 ;Bertero et al., 1988 ;Menke, 1989 ). The degree of regularization is controlled by the regularization parameter, often denoted as λ (lambda). The larger this parameter, the lower the spatial resolution and/or goodness-of-fit, and the higher the smoothness of the source distribution.
The details of regularization procedures are not within the scope of this paper. However, for a fair methods comparison it is important to choose the regularization parameter in similar ways for all methods. This can for example be accomplished via the (assumed) SNR of the data. Furthermore, the regularization parameter determines the effect of the noise covariance matrix on the inversion matrix. If λ = 0 , the noise covariance matrix does not contribute at all. The larger λ, the more it matters. Thus, in order to compare methods and to generalize results to other data sets, the structure and possible regularization of the noise covariance matrix has to be comparable.

Linear source estimation methods
In this section we will describe some of the most common source estimation methods in the linear framework introduced above, namely L2-Minimum-Norm-type methods and beamformers. To keep the notation simple, we will describe these methods for the case of fixed source orientations (i.e., one scalar value per source grid point modeling the amplitude of the current dipole at this location). The use of the resolution matrix and resolution metrics to evaluate spatial resolution is applicable to any other linear distributed source estimation method as well.

L2 -minimum -norm -type methods
The name of L2-MNE reflects the fact that its solution has the minimal L2-norm among all possible solutions to Eq. (2) . This property does not seem to fit any specific physiological or physical constraints about brain activity. However, L2-MNE also produces the resolution matrix which is closest to the ideal case, i.e., the identity matrix, in the leastsquares sense ( Menke, 1989 ) (We provide a derivation in the Appendix). This property suggests that L2-MNE provides an optimal trade-off of different spatial resolution criteria. It has often been pointed out that L2-MNE shows a strong tendency to localize sources too close to the sensors. Two strategies have been proposed to alleviate this: depth-weighting  and (noise-)normalization ( Dale et al., 2000 ;Pascual-Marqui, 2002 ).
A common expression for the unweighted L2-MNE is where is the regularization parameter and ( m x m ) the noise covariance matrix. The resolution matrix in this case is which is a symmetric matrix, and therefore PSFs and CTFs for elements i are the same. Depth-weighting can be implemented by introducing a diagonal weighting matrix D into Eq. (8) : The values in D can be chosen to "boost " the impact of deep sources, e.g., based on the vector norms of the leadfield columns (topographies of deeper sources tend to have smaller norms). Note that the same weighting strategy can be chosen to include a priori information about source locations, e.g., from fMRI data: the more we boost activation of certain sources via D , the more we bias our source distribution towards these sources, possibly even when these sources are not active ( Ahlfors and Simpson, 2004 ). Thus, depth-weighting can also be interpreted as using a priori information to bias the source estimate towards deeper sources. For eLORETA, a diagonal weighting matrix D is computed in an iterative process to optimize the localization error of PSFs for the estimator ( Pascual-Marqui et al., 2011 ).
Another way to change the localization properties of the L2-MNE estimator is (noise-)normalization ( Dale et al., 2000 ;Pascual-Marqui, 2002 ). In this case, the L2-MNE matrix is multiplied by a diagonal weighting matrix : The resolution matrix then changes accordingly: Thus, we may find a W that improves some desirable properties of the resolution matrix. Because the matrices for the methods considered here are diagonal, each row i of the MNE resolution matrix is scaled by the factor . As a consequence, the shapes of the CTFs (rows of R ) do not change ( Hauk et al., 2011 ). Only the shape of PSFs (columns of R ), and therefore potentially locations of peaks and their spatial extent, is affected by this normalization procedure.
The two most popular methods of this type are dynamic statistical parametric mapping (dSPM) and standardized low-resolution electromagnetic tomography (sLORETA): For dSPM, the normalization matrix contains the minimum norm estimates of the noise at each source ( Dale et al., 2000 ), derived from the noise covariance matrix, i.e., For sLORETA ( Pascual-Marqui, 2002 ), the normalization uses the diagonal of the MNE resolution matrix It has been shown that dSPM and sLORETA have better peak localization performance for PSFs than L2-MNE ( Hauk et al., 2011 ), and even that sLORETA has zero dipole localization error ( Hauk et al., 2011 ;Pascual-Marqui, 2002 ) under ideal circumstances. However, as we will see localization performance alone does not allow similar conclusions about other aspects of PSFs (e.g., spatial extent, local extrema etc.), or about the shape of CTFs ( Hauk et al., 2011 ).

Beamforming and Spatial Filtering
Linear inverse methods can be interpreted as "spatial filters " or "software lenses " that attempt to estimate activity from sources of interest, while suppressing interference from all other possible sources. These spatial filters can be arranged as rows of an inverse matrix to provide the whole source distribution. In this logic, the ideal inverse matrix should yield a resolution matrix as close as possible to the identity matrix. Thus, an individual spatial filter's CTF should be as close as possible to the corresponding row of the identity matrix. If closeness to the identity matrix is required in the least-squares sense, this yields the L2-MNE solution ( Backus and Gilbert, 1968 ;Hauk, 2004 , see also Appendix A; Menke, 1989 ).
Beamformers approach this problem differently. They start from the requirement that a spatial filter for source i should be sensitive to activity from this source, i.e., it should yield the value 1 if this source is active with unit strength (unit gain beamformer) Van Veen et al., 1997 ): At the same time, the beamformer attempts to suppress noise and activity from other sources. In contrast to minimum-norm type methods, it does not do this by modeling other sources via the leadfield matrix, but by assuming that they are captured by the data covariance matrix D . Note that this matrix contains the noise covariance as well as possible signal covariance from other sources of interest. The requirement to suppress the effect of other activity on the estimate for source i can thus be written as This results in the unit-gain linearly constrained minimum variance (LCMV) beamformer: Variations exist that optimize beamformers for pairs of dipolar sources ( Brookes et al., 2007 ) and for time-and frequency-dependent beamforming ( Dalal et al., 2008 ;Woolrich et al., 2013 ).
Because the estimator is a vector that is multiplied to the data matrix ( . ) , we can use the concept of CTFs to describe its spatial resolution. In principle, the superposition principle holds, since the contributions of different sources add up. However, the situation is more complicated in this case: Because we are using the data (rather than noise) covariance matrix, the estimator depends on the sources of interest. Whenever the data change, so does the estimator. This is why beamformers are often called "adaptive spatial filters " ( Sekihara and Nagarajan, 2008 ) -they adapt to the data. In contrast, minimum-norm-type estimators are "static ", since they are based only on the leadfield and noise covariance matrices. This affects the generalizability of any resolution analysis: While the result for linear minimum-norm type methods apply to other data sets with similar measurement configurations, head models and noise statistics, the results for beamformers are only valid under those conditions that are represented in the data covariance matrix. Thus, while beamformers yield a linear transformation of the data, their adaptivity to data makes their PSFs and CTFs hard to generalize to different data sets, or even different analyses of the same dataset.
We evaluated two types of beamformers in this study: One using a covariance matrix based on pre-stimulus baseline activity, and one using a covariance matrix from a post-stimulus interval (referred to as "pre " and "post " in the following sections, e.g., in Fig. 7 )). For the former, we used the noise covariance matrix instead of the data covariance matrix in Eq. (18) : As we show in the Appendix, this resembles the maximum likelihood estimator (MLE) for a single dipole strength ( Lütkenhöner, 1998a ). This can still be interpreted as a special case of an LCMV beamformer. The data covariance matrix in any given experiment is the noise covariance matrix plus the signal covariance (assuming additivity of noise and signal, which is a common assumption). If we were going to estimate the sources in our baseline interval, then this filter could be interpreted as a beamformer, since in this case N would represent the activity due to the sources of interest. Baseline activity can also be considered as a type of resting state activity. We provide a more detailed description of this issue in Appendix B.

Resolution metrics
A meaningful resolution analysis of EEG/MEG source estimation methods should include metrics for the resolution categories localization accuracy and spatial extent. Table 1 lists several possible metrics, most of which have already been used in previous studies (but not necessarily together). We will use a subset of those in our Results section, as highlighted in the table.

Simulation set -up
In order to illustrate the use of PSFs and CTFs we computed individual PSFs and CTFs in Fig. 3 for the EEG/MEG sample dataset 2 provided by the MNE-Python software package. This data set contains MEG and EEG data (Elekta Vectorview) for simple auditory and visual stimuli as well as structural MRI data for one participant.
For a more extensive and quantitative methods comparison, we computed resolution metrics, histograms and correlations in the following figures using the open EEG/MEG data provided by Wakeman and Henson (2015) . This dataset consists of data for 16 participants from which simultaneous EEG and MEG were recorded in an Elekta Neuromag Vectorview scanner (70 electrodes, 102 magnetometers, 204 planar gradiometers), in addition to MRI, fMRI and DTI data (the MRI data were used to create head models for source estimation; the fMRI and DTI data were not used for the present study). Details of the data acquisition parameters are available in the previous publication ( Wakeman and Henson, 2015 ). The study reports that 0 to 14 MEG and 0-4 EEG channels per dataset were interpolated. Stimuli in this experiment were pictures of faces and scrambled faces, and participants were young and healthy adults. In the following, we will only summarize the most important parameters for our purposes.
For our simulations, we only required the sensor configurations, structural MRI images, as well as noise covariance matrices (details below). Pre-processing of EEG/MEG and MRI data, including source estimation, was carried out in MNE-Python software. Resolution analysis was performed in MNE-Python (version 0.23) ( Gramfort et al., 2013 ). The software used for our resolution analysis is available online ( https://github.com/olafhauk/EEGMEG_ResolutionAtlas ). The scripts are based on the corresponding example scripts on the MNE-Python web-site 3 .
The noise covariance matrices for each data set were computed concatenating pre-stimulus baseline intervals of 200 ms duration before stimulus onset. For the data covariance, this was done for a latency interval 50 to 250 ms post-stimulus presentation. The covariance matrix was estimated using the 'shrunk' regularization method ( Engemann and Gramfort, 2015 ). For regularization of the inverse estimators, the default signal-to-noise ratio in the MNE software was used (SNR = 3). This value is representative for many event-related EEG/MEG studies ( Farahibozorg et al., 2018 ;Molins et al., 2008 ;Rahimi et al., 2022 ).
High-resolution structural T1-weighted MRI images from the open dataset were acquired in a 3T Siemens Tim Trio (Siemens, Erlangen, Germany) scanner at the MRC Cognition and Brain Sciences Unit, University of Cambridge, UK, using a MPRAGE sequence, and processed using MNE-Python software. This included the segmentation of the cortical surface as source space, and the creation of a three-compartment boundary element model (BEM) for the forward solution. MEG sensor configurations and MRI images were co-registered based on the matching of about 50-100 digitized locations on the scalp surface with the reconstructed scalp surface from the structural MRI.
The source space consisted of the cortical surface with 4098 vertices per hemisphere (octahedron spacing 6). BEMs were created using linear collocation with 5120 triangles per surface. The surfaces represented the inner skull, outer skull and outer skin, respectively. Source orientations for the forward solution were constrained to be perpendicular to the cortical surface ( "fixed orientation constraint "). Inverse operators and forward solutions were used to compute the resolution matrix and the resolution metrics for each individual subject. Those results were then morphed based on a spherical representation of the cortex computed using the spherical registration of Freesurfer software ( Greve et al., 2013 ) and with ( ) ≥ * max ( ( ) ) ( ) : PSF or CTF depending on location r 0 < x < 1: fraction of maximum Relative area/volume ∑ i with ( s ) ≥ * max ( ( ) ) : surface area/volume represented by source i A: overall surface area/volume 0 < x < 1: fraction of maximum Fig. 3. PSFs and CTFs for two MNE-type source estimation methods in three brain areas. First Row: CTF and PSF for MNE in anterior temporal lobe (ATL), pars triangularis (PT) and pars opercularis (PO), respectively. Because MNE's resolution matrix is symmetric, its PSFs and CTFs can be visualized in the same image. Second and Third Row: As in the first row, but for eLORETA's PSFs and CTFs in separate rows. The displayed distributions are the first principal vectors for all PSFs/CTFs within the corresponding regions-of-interest (ROIs). The ROI borders are shown on the inflated cortical surface. PSFs/CTFs are displayed as signed distributions, i.e. red and blue colors indicate current flow in opposite directions with respect to the cortical surface (out-and in-flowing, respectively). MNE: minimum-norm estimator; eLOR: exact low resolution electromagnetic tomography. PSF: point-spread function; CTF: cross-talk function. averaged across subjects. Results will be shown on the average brain. We do not expect large qualitative differences between left and right brain hemispheres, and therefore only the left hemispheres will be shown for simplicity. Fig. 3 illustrates how individual PSFs and CTFs can be displayed for MNE and eLORETA in different regions-of-interest (ROIs) in order to address questions with respect to spatial resolution. We assume a scenario where we are mainly interested in activity in three ROIs, namely anterior temporal lobe (ATL), pars triangularis (PT) and pars opercularis (PO) (ROIs displayed as borders and indicated by black arrows). These ROIs could be of interest, for example, in a language experiment that attempts to disentangle semantic representation (involving ATL) and semantic control (involving IFG) processes (e.g. Jackson, 2021 ). Our specific ROIs were taken from the Desikan-Killiany Atlas ( Desikan et al., 2006 ) (with ATL as the anterior quarter of the middle temporal parcel).

Illustration of PSFs and CTFs for three ROIs
Using MNE we do not assume that activity occurs only in these ROIs (or that there is a certain (maximum) number of sources), and therefore we would like to control leakage from or into all other brain regions. The results in Fig. 3 are not meant to be quantitative or generalizable results, but shall illustrate how PSFs and CTFs can help with the interpretation of source estimation results for individual datasets. More generalizable results will be presented in the following figures.
Because ROIs typically comprise many point sources, we here used principal component analysis (PCA) to find the most representative PSFs and CTFs across all vertices for each ROI. For this purpose, we computed a PCA across all PSFs and CTFs in each ROI and chose the first principal vector, i.e., the distribution that explained most variance across all PSFs or CTFs that entered the PCA. These results were computed for one individual measurement configuration provided as a sample dataset in MNE-Python software package. PSFs/CTFs are displayed as signed distributions, i.e., red and blue colors indicate current flow in opposite directions with respect to the cortical surface (out-and in-flowing, respectively).
The first row of Fig. 3 shows the PSFs and CTFs for unweighted and unnormalized MNE. Because MNE's resolution matrix is symmetric, each distribution represents a PSF as well as CTF for the corresponding ROI. The aim is to find an estimator whose PSF and CTF are mostly confined to the target ROI, in the left-most column this is the ATL. This is mostly the case for ATL, although some leakage occurs in areas posterior to the ROI (in middle temporal regions) but also in inferior frontal regions (e.g., around the border between pars triangularis and opercularis). In the PSF interpretation this means that true activity in ATL may also affect source estimates posterior to ATL and to a lesser degree in inferior frontal areas. Conversely, in the CTF interpretation this means that activity posterior to ATL may leak into source estimates for ATL, and that weaker leakage can also occur from inferior frontal regions. In general, we can be fairly confident that we can separate activity from ATL and PT and PO, but we should be aware that there could be leakage from around ATL.
Turning our eye to PT, we also note that the major peaks of its PSF and CTF lie within its border. Most leakage occurs just around its borders, which also includes PO and to a lower degree ATL. As for ATL, this gives us some confidence that PT activity can be separated from PT and ATL, but that some leakage from around its borders is possible.
The situation is different for region PO. While the peaks of its PSF and CTF are within its borders, major peaks also occur outside, especially in PT but also in ATL and superior temporal lobe. We can conclude that activity from PO and ATL can be separated fairly well, but for PO and PT leakage is a serious concern.
We can now ask whether the picture changes significantly if we use eLORETA instead of MNE. PSFs and CTFs are different for eLORETA, and therefore presented in separate rows. eLORETA has the property of zero localization error for point sources (( Pascual-Marqui et al., 2011 ), but note that this is only the case for PSFs but not CTFs as shown above). As for MNE, the main peaks of PSFs and CTFs occur in the corresponding target ROIs, and the general pattern looks similar for eLORETA and MNE. However, eLORETA's PSFs are more spread out than MNE's, resulting in more leakage from the target ROIs to neighboring areas, including deeper areas in the Sylvian fissure. The CTFs for ATL and PT are very similar to MNE's. The CTF for PO receives less leakage from PT, albeit at the expense of more leakage from areas dorsal to PO and also from more distant locations in inferior parietal cortex. This suggests that the zero localization error of PSFs comes at the expense of more widespread distributions and therefore more leakage, and that the shape of CTFs can be changed but not arbitrarily optimized.
The choice of a best method in this scenario ultimately depends on the researcher's goals. If the goal is to minimize leakage only taking the three target ROIs into account, then eLORETA could be the better choice because of its CTF for PO. However, if we also want to minimize leakage into or from areas outside the three target ROIs, then MNE could be the safer choice.

Whole -cortex distributions of resolution metrics
The previous figure illustrated how plotting individual PSFs and CTFs for specific regions is a practical way to evaluate spatial resolution for specific research questions. However, it is impractical to plot these functions for all vertices, voxels or ROIs across the brain. Ideally, we would like to define some meaningful metrics that turn the most relevant features of PSFs and CTFs into a single number which we can then plot across the whole brain (see Table 1 for examples). These metric distributions can then also be analyzed and compared across different methods, measurement configurations, etc., as demonstrated in several previous studies (e.g. Hauk et al., 2011 ;Iivanainen et al., 2017 ;Molins et al., 2008 ). Fig. 4 presents the distributions of the resolution metrics peak localization error (PLE, left half) and spatial deviation (SD, right half) for PSFs and CTFs for five MNE-type methods (different rows). The corresponding histograms are presented in Fig. 5 . The PLE for MNE's PSF (top left) shows the expected pattern: It is close to zero in locations close to the sensors (e.g., around the gyri) and increases with increasing source depth into the sulci. It is largest in deep areas of the Sylvian fissure, with errors exceeding 5 cm. The corresponding histogram in Fig. 5 shows a peak at low PLE values below 2 cm, but also a long tail exceeding 5 cm. Mean and median values of the distribution are 1.7 cm (S.D. 1.3 cm) and 1.3 cm, respectively. Thus, on the basis of the well-known bias of unweighted MNE toward brain locations close to the sensors, we can conclude that peaks of PSFs occur in approximately correct locations for superficial sources, but will systematically be shifted towards the sensors for deep sources.
The other four methods have all been proposed to reduce the overall localization error. Indeed, depth-weighted MNE has generally lower PLE, but still some high values in deep brain areas. The corresponding histogram confirms that the range of PLE values is lower than for unweighted MNE, and especially the distribution's tail is shorter (Mean(SD)|Median: 0.8(0.6)|0.6 cm). dSPM has low PLEs in deep brain areas, but at the expense of moderate values in more superficial regions where MNE's PLE is lower. The corresponding histogram shows an approximately symmetrical bell-shaped distribution (Mean(SD)|Median: 1.2(0.5)|1.2 cm). Both sLORETA and eLORETA show the promised zero peak localization error for PSFs.
The comparison of PLE distributions for CTFs across methods in Fig. 4 is straightforward: They are all highly similar or in fact the same, and also highly similar to MNE's PSF just described. They all Fig. 4. Comparison of MNE-type methods with respect to spatial resolution metrics peak localization error and spatial deviation. The two left-most columns show peak localization error and the two right-most columns spatial deviation (see Table 1 ) in a lateral view of the left hemisphere. PSFs and CTFs are presented in different columns. The rows represent results for different source estimation methods. MNE: minimum-norm estimator; dSPM: dynamic statistical parametric mapping; e/sLOR: exact/standardized low resolution electromagnetic tomography. PSF: point-spread function; CTF: cross-talk function.
show low values in superficial brain areas and high peak localization errors in deeper brain areas. Their histograms are also highly similar with the same central tendencies (Mean(SD)|Median: 1.7(1.3)|1.3 cm). Thus, irrespective of the method used, source estimates at deep brain locations will always receive significantly more leakage from sources that are closer to the sensors ( Grave de Peralta Menendez et al., 1997 ;Hauk et al., 2011 ;Krishnaswamy et al., 2017 ;Liu et al., 1998 ).
The differences among methods with respect to spatial deviation (SD) of PSFs in Fig. 4 are more subtle than for PLE, but still clearly visible. MNE shows the lowest values in superficial areas across the cor-

Fig. 5.
Comparison of MNE-type methods with respect to histograms of spatial resolution metrics peak localization error and spatial deviation. Spatial arrangement of panels and abbreviations as in Fig. 4 . The vertical lines indicate the mean of the distribution. The numbers next to the vertical line are the mean (standard deviation) and median of the distribution. Histogram counts ( y -axis) were divided by 1000 for display. tical surface, but large value in deeper areas such as the Sylvian Fissure. This is confirmed by the corresponding histogram in Fig. 5: MNE has the lowest median (2.6 cm) and one of the lowest mean SDs (2.9 cm), but the distribution also has the longest tail, reflected in a high standard deviation (0.8 cm). The differences among the other methods are more subtle. All of them show highest values in deeper areas of the Sylvian Fissure and along the inferior temporal lobe, and this is particularly pronounced for dSPM and sLORETA. This is confirmed by the histograms, which show comparatively high values for dSPM (Mean(SD)|Median: 3.2(0.5)|3.1 cm) and sLORETA (Mean(SD)|Median: 3.0(0.4)|3.0 cm), lower values for dwMNE (Mean(SD)|Median: 2.9(0.3)|2.9 cm) and still lower for eLORETA (Mean(SD)|Median: 2.8(0.4)|2.8 cm).
As for PLE, the SD distribution for CTFs is highly similar across methods and also to those of MNE's PSF. The distribution has a long tail with largest values in deep brain areas such as the Sylvian Fissure. Means and medians are 2.9 cm and 2.6 cm, while standard deviations vary slightly across methods (MNE|dSPM|sLORETA: 0.8 cm; dwMNE|eLORETA: 0.7 cm).
Beamformers seemingly take a very different approach to solving the EEG/MEG inverse problem and can therefore be expected to show different PSFs and CTFs. Fig. 6 shows example PSFs and CTFs for two LCMV beamformers computed for the MNE-Python sample data set (similar to Fig. 3 ). The beamformers were computed with data covariance matrices from epochs corresponding to auditory and visual events, respectively. PSFs and CTFs are shown for a left occipital (Occ) and posterior superior temporal (PST) ROI. At the bottom of the figure we also show PSFs for a beamformer based on a data covariance matrix from the pre-stimulus baseline period as well as for unweighted L2-MNE. Fig. 6. PSFs and CTFs for two LCMV beamformers and MNE in two brain areas. PSFs and CTFs were computed with data covariance matrices for auditory (Aud) and visual (Vis) events, as well as for a pre-stimulus baseline interval (Base). MNE is shown at the bottom for comparison. The displayed distributions are the first principal vectors for all PSFs/CTFs within the corresponding regions-ofinterest (ROIs). The ROI borders are shown on the inflated cortical surface. PSFs/CTFs are displayed as signed distributions, i.e. red and blue colors indicate current flow in opposite directions with respect to the cortical surface (out-and in-flowing, respectively).
In general, PSFs and CTFs for all beamformers within the same ROI look similar, suggesting that the latency range for data covariance computation only has a subtle effect on the estimator. The result for the beamformer obtained from baseline data indicate that this beamformer is a good approximation of beamformers obtained for event-related data.
PSFs and CTFs for Occ are mostly confined within the ROI, with some leakage around its border. PSFs and CTFs for PST are more widespread, also covering areas outside the ROI such as in inferior parietal and middle temporal regions. However, there is hardly any overlap between PSFs and CTFs for the two regions, which mean that their activities can be well-separated from each other. The PSFs for PST also spread into deeper areas of the Sylvian Fissure, while this is not the case for the cor-responding CTFs. This demonstrates again that CTFs are fundamentally constrained by the leadfields.
For comparison we present PSFs/CTFs for MNE at the bottom of Fig. 6 . They resemble the corresponding results for the beamformers, but are less spread out.
The distributions of resolution metrics for LCMV beamformers are shown in Fig. 7 , as well as their histograms (similar to Figs. 4 and 5 for MNE-type methods). The results for the pre-and post-stimulus beamformers are similar, indicating that the exact choice of the latency interval for the data covariance does not significantly affect a beamformer's spatial resolution properties. The PLE of these beamformers is low, below 1 cm (for pre-stimulus beamformer PSF exactly zero; post-stimulus  Fig. 4 , but note different scaling). "Pre " and "Post " refer to beamformers computed for pre-(baseline) and post-stimulus intervals, respectively. Bottom two rows: Histograms corresponding to the distributions above (similar to Fig. 5 ). The vertical lines indicate the mean of the distribution. The numbers next to the vertical line are the mean (standard deviation) and median of the distribution. Histogram counts ( y -axis) were divided by 1000 for display. beamformer 0.6(0.2), 0.6; CTF pre and post: 4.3(1.1) cm, 4.1 cm). However, this is accompanied by high PLEs for CTFs as well as high spatial deviation values for both PSFs and CTFs across the whole cortex, with mean and median values higher than for any MNE-type methods in Figs. 4 and 5 (PSF pre: 5.0(0.4) cm, 4.9 cm; post: 5.0(0.4) cm, 5.0 cm; CTF pre: 5.3(0.7) cm, 5.2 cm; post: 5.3(0.7) cm, 5.3 cm). This reflects the fact that the derivation of LCMV beamformers does not explicitly include constraints on the spatial extent of either PSFs or CTFs.
We already noted several times that source depth strongly affects our resolution metrics. We report the correlations between resolution metrics and source depth for all six methods in Fig. 8 . For PSFs, MNE shows the highest correlations both for PLE (mean 0.92(standard deviation 0.01)) and for SD (0.84 (0.02)). dwMNE shows a moderate correlation for PLE (0.32(0.03)) but low and even slightly negative correlation for SD (-0.1(0.08)). dSPM has a low and negative correlation for PLE (-0.09(0.05)) but positive correlation for SD (0.19(0.07)). sLORETA, eLORETA and the pre-stimulus beamformer have zero localization error and therefore no meaningful correlation with source depth. However, sLORETA and eLORETA show moderately high correlations with SD (0.45(0.04), 0.29(0.05), respectively), and the beamformers show low negative correlations with high standard deviations (-0.23(0.26), -0.24(0.21), respectively). The correlations for the two beamformers are very similar, as expected from the previous results with respect to their resolution metrics.
For CTFs, the correlations with source depth are similarly high for all MNE-type methods both for PLE (all above or equal 0.92) and SD (all above or equal 0.8). Interestingly, correlations for beamformers are moderate for PLE (for both around 0.49) and SD ( ∼0.29), with high standard deviations compared to MNE-type methods (around 0.4 for both).
This high standard deviation may reflect the adaptive nature of LCMV beamformers, as they crucially depend on the data covariance matrices of the individual datasets. It is noteworthy that among MNE-type methods dwMNE and dSPM also show low negative correlations for PSFs (dSPM for PLE and dwMNE for SD, respectively). Depth-weighting for dwMNE and noise-normalization for dSPM have both been suggested to compensate for MNE's localization bias towards superficial sources. It is possible that such a compensation leads to an improvement of resolution metrics for deep relative to superficial sources in the case of PSFs. However, for CTFs all correlations were clearly positive. Note that the smaller depth-correlations for beamformers in the case of CTFs must be interpreted in the light of their generally larger SD values with respect to MNE-type methods.

Discussion
We presented a tutorial-style introduction to linear EEG/MEG source estimation and resolution analysis, and described several widely used distributed source estimation methods in this framework. We highlighted the importance of the resolution matrix, the need to evaluate spatial resolution on the basis of both point-spread and cross-talk functions (PSFs/CTFs) as well as multiple resolution metrics. We evaluated spatial resolution for several MNE-type methods and LCMV beamformers on the metrics peak localization error (PLE) and spatial deviation (SD), and particularly focused on their relationship with source depth.
Our results confirmed our assertion from the tutorial that source estimation methods should be evaluated based on CTFs as well as PSFs. PSFs tell us how one source would leak into others, and CTFs show how other sources would leak into the estimate for the source of interest Fig. 8. Correlations between spatial resolution metrics and source depth for seven source estimation methods. The bar graphs depict the correlations between peak localization error (top row) and spatial deviation (bottom row) for PSFs (left column) and CTFs (right column). Different bars represent means across datasets for different MNE-type methods. Error bars are standard deviations across participants. dMNE: depth-weighted MNE; sLOR/eLOR: s/eLORETA; Lpre/Lpost: LCMV beamformers computed with covariance matrices for pre-/post-stimulus interval. ( Fig. 2 ). Thus, we may be able to construct an estimator that produces a peak around a source of interest given that this source is truly active, i.e., it may produce a good PSF. However, it may receive a lot of leakage from other sources that may or may not be active at the same time, i.e., it may have a bad CTF. This leakage from other sources may even be larger than the desired "leakage " from the source into itself, reflected in large peaks of the CTF at locations distant from the source of interest. This is what we observed for all methods in deeper brain locations, based on the PLEs for CTFs ( Figs. 4 and 7 ). Based on the fact that the sensitivity of EEG/MEG recordings falls off rapidly with distance from the sensors ( Goldenholz et al., 2009 ;Hauk et al., 2011 ;Hillebrand and Barnes, 2002 ), and confirmed by our individual examples of PSFs and CTFs ( Figs. 3 and 7 ), this likely reflects strong leakage from sources in superficial brain regions into estimates for sources in deep brain regions (such as the Sylvian Fissure). Thus, a good localization error for PSFs is not enough to conclude that a source estimation method has good spatial resolution. Spatial resolution has to be evaluated based on the shapes of both PSFs and CTFs.
This still leaves us with the question whether our methods comparison yielded a method with best spatial resolution. Localization accuracy measured by PLE for PSFs can clearly be improved using weighted and normalized MNE-type methods ( Figs. 4 and 5 ) as well as beamformers ( Fig. 7 ), and the pre-stimulus baseline beamformer, sLORETA and eLORETA even have zero PLE. However, this is not reflected in a comparable improvement of spatial extent for PSFs, or generally for CTFs. In our comparison, eLORETA shows the lowest mean SD for CTFs, but MNE has the lowest median SD. Together with its zero PLE eLORETA may therefore be considered the optimal method in this comparison. We would like to highlight again that researchers can apply these methods to their own datasets in order to find the best method for their task at hand. For example, we computed our metrics across the whole cortical surface using fixed source orientations. Some researchers may prefer loose orientations or a volumetric source space, others may only be interested in a subset of regions (e.g., only visual brain areas).
Researchers may also choose different resolution metrics. We provided some examples in Table 1 . Previous studies have for example used the center-of-gravity rather than the peak to compute localization error, or spatial extent has been evaluated using measures of area ( Iivanainen et al., 2017 ). In previous studies we also used a measure of relative amplitude for PSFs/CTFs Hauk et al., 2011 ), which we omitted from the present study for simplicity. However, as our Fig. 1 B illustrates, a large difference in amplitudes between PSFs and CTFs can also have significant impact on their distinguishability.
We highlighted some special properties of LCMV beamformers in the framework of linear resolution analysis. While these estimators result in linear transformations of the data, they are adaptive and therefore data-dependent, and it is not straightforward to generalize their results across different data sets and applications. However, in our PSF/CTF examples ( Fig. 6 ) and for the resolution metrics ( Fig. 7 ) the differences between beamformers based on different data covariance matrices were subtle. This indicates that the data covariance matrices are dominated by background brain activity rather than relatively subtle evoked responses. Thus, it is possible that the beamformer based on pre-stimulus baseline activity is a good approximation for beamformers in general, and therefore suitable for standardized methods comparisons including beamformers. This should be tested in more detail in future studies.
The leakage problem has received particular attention in the domain of brain connectivity, as it may give rise to spurious connectivity or "ghost interactions " between spatially distinct sources ( Colclough et al., 2015 ;Farahibozorg et al., 2018 ;Hipp et al., 2012 ;Palva et al., 2018 ). LCMV beamformers are especially popular for the analysis of oscillatory brain activity and functional connectivity such as for resting state data ( Brookes et al., 2011 ;Colclough et al., 2016 ). The sensitivity of an estimator to source activity in a certain region to activity from other regions, as well as its ability to leak activity into other regions, depends on the spatial extent of its CTF and PSF, respectively. We assessed spatial extent using the metric spatial deviation (SD). Our two beamformers produced larger SD values than the MNE-type methods for both PSFs and CTFs. PLE for PSFs was low, but larger for CTFs. This is likely due to the fact that beamformer estimators are not constrained by the whole leadfield to optimally suppress all sources of no interest in the source space, but only those represented in the data covariance matrix. It has been argued previously that "this is the reason why adaptive spatial filters have a strange-looking beam response " (another term for cross-talk function ( Sekihara and Nagarajan, 2008 , p. 42). This raises the question whether LCMV beamformers are optimal for brain connectivity analysis with respect to their spatial resolution. Again, the concepts described in this paper will allow researchers to evaluate this in their own datasets.
We argued above that it is hard to determine a "best " method from our comparison, and that this depends on the specific goals and measurement set-ups of the researcher. However, we were able to determine limits that none of these methods can surpass, especially with respect to deep sources. While some methods achieve better localization accuracy in PSFs for deep sources, none of them achieved this for CTFs. All methods showed larger spatial extent for deep sources. The correlations between resolution metrics and source depth was positive for CTFs for all methods. We showed that CTFs are weighted sums of the rows of the leadfield matrix. Thus, if there is no linear combination of those rows that can achieve peaks in depth, then no method can achieve low PLE in depth. As a consequence, all estimates for deep sources may be contaminated by leakage from more superficial sources. This is supported by previous findings that topographies of deeper and superficial sources can be highly correlated, and therefore given realistic noise levels and their unknown spatial extent their sources are impossible to separate unless one uses a hierarchical and sparse modeling approach based on specific prior knowledge and modeling assumptions ( Krishnaswamy et al., 2017 ;Liu et al., 1998 ). We conclude that deep sources are hard to distinguish from superficial sources without additional well-justified constraints.
We performed our resolution analysis using combined EEG and MEG, which with the currently available technology represents the best-case scenario. We regularized our estimators in a comparable manner using standard methods and assuming a typical signal-to-noise ratio. However, as already pointed out by Fuchs et al. (1999) , the resolution matrix approach does not explicitly include noise in the computations of pointspread and cross-talk functions. Thus, our results are only perfectly accurate if our regularization procedure has completely suppressed the effect of noise on the source estimate, which is unlikely in realistic scenarios. Nevertheless, we hold the view that at low to moderate noise levels as in our study, these results are still meaningful for methods comparison and to determine an upper bound for the resolution of EEG/MEG source estimates. Recent studies have proposed ways to incorporate noise explicitly into the resolution matrix, e.g., by using an "empirical resolution matrix " ( Krishnaswamy et al., 2017 ;Samuelsson et al., 2021 ), which is computationally more expensive. We are planning to use this approach in the future, especially to investigate the effect of different types of noise on source estimation.
The empirical resolution matrix mentioned above can also be computed for non-linear source estimation methods, e.g., those that incorporate constraints on sparsity. Sparsity (i.e., the assumption of a small number of focal sources) is sometimes used as a nonlinear constraint to improve localization accuracy, e.g., in L1-minimum-norm-type methods ( Huang et al., 2014 ;Owen et al., 2012 ;Uutela et al., 1999 ), or multiple sparse priors ( Friston et al., 2008 ). However, it is important to note that for nonlinear methods the superposition principle does not hold, and such a resolution matrix would only reflect the specific source scenario used for its computation and cannot be generalized to more complex source configurations. Nonlinear methods are therefore hard to evaluate in a generalizable way, and their results strongly depend on the validity of the underlying assumptions. The simulation set-ups for nonlinear methods have to be chosen and evaluated carefully in order to represent source scenarios that can be realistically assumed to reflect brain activity in particular types of experiments.
Recent EEG/MEG research has increasingly focused on patterns of brain activations ( Kietzmann et al., 2019 ;Stokes et al., 2015 ). While any pattern can be modelled as the sum of point sources, it is not obvious to predict how an activity pattern in one ROI will affect the decoding of patterns in other ROIs. This will depend on multiple features of these patterns, e.g., their spatial frequency, homogeneity, etc. This requires a statistical approach, i.e., the simulation of a large number of representative patterns in one ROI and their effects in other ROIs, in order to compute a "multivariate pattern resolution matrix ". Our present study lays the ground for future work in this direction.
Finally, resolution metrics can be used to compare different sensor configurations with respect to their spatial resolution in different brain areas. This has already been exploited in previous studies comparing MEG-only with MEG + EEG setups Molins et al., 2008 ), as well as to evaluate novel non-cryogenic on-scalp magnetometers systems ( Iivanainen et al., 2017 ). The concepts presented in this paper can therefore help researchers to decide whether their research questions can be addressed using widely available and relatively cheap EEG, requires more costly MEG equipment (possibly in combination with EEG), or even requires novel on-scalp MEG systems. In the latter case, it can inform researchers about the most cost-and time-efficient sensor arrangement of on-scalp sensors (MEG as well as EEG) in order to target particular brain areas.
We conclude that we described a comprehensive, objective and computationally efficient approach to evaluate the spatial resolution of EEG/MEG source estimation methods that can easily be applied by experimenters to their task at hand. We hope that our tools will become widely available in all major EEG/MEG software packages. This will help EEG/MEG researchers in the interpretation of source estimation results and to determine the optimal measurement configurations and analysis strategies for their individual purposes.

Data and code availability statement
The EEG/MEG data used in this study are openly available datasets, and are also available from the MRC Cognition And Brain Sciences' data repository on request.

Credit authorship contribution statement
where N is the number of sources and the delta operator (1 if i = j , 0 otherwise).
Since the leadfield L is fixed and only the inversion matrix K is variable, and every element of the sum is greater than or equal to zero, the minimization of A1 can be achieved for every row of R independently: The solution can be found by taking the derivative of the minimization function and finding its zero point: The solution of this equation is Since ( ) −1 is independent of the source index i , and . is the i th column of the leadfield matrix, we can write the solution for all sources as which is the unweighted and unnormalized L2-minimum-norm estimator.

B) Maximum likelihood estimator for dipole strength as a special case of LCMV beamforming
Here, we demonstrate a few properties of LCMV beamformers, in particular their relationship to maximum-likelihood dipole strength estimation.
As in Eq. (19) , the spatial filter for an LCMV beamformer at location i is LCMV .

)
When data are pre-whitened, e.g., for source estimation, then the leadfield and data (including data covariance matrix) are multiplied by the square root of the symmetric noise covariance matrix, i.e.,

)
We exploited that −1∕2 N is symmetric and that the inverse of the product of invertible matrices is the product of the matrix inverses.
In the special case where the data covariance matrix is the same as the noise covariance matrix, where ̃ is the whitened leadfield . This is the i th column of the whitened leadfield matrix normalized by the square of its norm.
As an inverse matrix applied to all sources, this can be written as with 2 a diagonal matrix with the inverse squared norms of leadfield columns ( 1∕ ̃ T . ̃ . ) as its diagonal elements. Thus, when the prewhitening matrix is the same as the data covariance matrix, the LCMV beamformer is the transposed column-normalized whitened leadfield matrix. This estimator has zero dipole localization error, and has been used to demonstrate that this property alone is not a sufficient criteria for good spatial resolution ( Grave de Peralta et al., 2009 ). This is for example the case when the whitening matrix is estimated from a whole resting state dataset and the beamformer then applied to the same dataset, or when the whitening matrix is estimated from pre-stimulus intervals and the beamformer then applied to those pre-stimulus intervals. If whitening matrix and data covariance matrix differ, the difference of the LCMV beamformer to this expression depends on the difference between 1∕2 N −1 D 1∕2 N and the identity matrix. Thus, the more similar N and D , the more similar the LCMV beamformer is to the leadfield expression above.
We usually assume additivity of noise and signal, i.e. D = N + S where S is the covariance matrix reflecting only the sources of interest.
If we assume that the noise covariance matrix estimated from prestimulus intervals (comparable to resting state data) captures dominant ongoing brain activity, and that evoked responses are relatively small with respect to this activity (i.e., S is small compared to N ) , then this version of the LCMV beamformer is a good model to characterize LCMV beamformers in a generalizable and reproducible manner. This was confirmed by the results in this study.
The formulation of the LCMV beamformer above assumes that data are pre-whitened using the noise covariance matrix. If we include this pre-whitening of the data into the above expression and rearrange it, exploiting that N is symmetric, we obtain which is the maximum likelihood estimator for the estimation of the amplitude of a single dipole i . This demonstrates the relationship between LCMV beamforming and dipole scanning. For the maximum likelihood estimation of a single dipole strength i , we want to minimize the function where is the source strength of dipole i , d is the data vector, the noise covariance matrix and . i the topography of source i . The solution to this problem is (e.g., Lütkenhöner, 1998a ) = . −1 i.e., this estimator is the same as the LCMV beamformer with a noise covariance matrix above. As before the closer a data covariance matrix is to the noise covariance matrix, the closer the LCMV beamformer filters are to the maximum likelihood dipole strength estimator. Finally, we consider the case of a data covariance matrix under the assumption that all sources in the model are randomly activated, with uncorrelated time courses in the noise source vector N ( ) , following a uniform Gaussian distribution. In that case, the data covariance matrix would be with <> representing the expectation value over time. If we insert this into the previous equation for the LCMV beamformer, we obtain the L2-MNE estimator except for a different scaling. While the scaling will affect the PSFs, it will not affect the shape of the CTFs, and the CTFs results for L2-MNE will also apply to this type of beamformer.

C) A brief primer on linear algebra for EEG/MEG source estimation
This appendix is intended as a brief primer for the linear algebra used to describe EEG/MEG source estimation in this manuscript, especially with respect to resolution matrix, point-spread and cross-talk functions. An extensive introduction to linear algebra can for example be found in Cohen (2021) .

Vectors, matrices, matrix multiplication
In everyday life, we commonly encounter "real " or "scalar " numbers (e.g., 1, 2.3, -10.3 etc.) and use them in basic algebraic operations (addition, subtraction, multiplication, division, etc.). A vector is an ordered list of such numbers In this case the vector x has n elements, i.e., it is of dimension n. Vector addition and subtraction work on an element-by-element basis, i.e., and correspondingly for subtraction. However, vector multiplication usually refers to the scalar product of two vectors as follows: The result is called scalar product because it is just a single number. If x and y have zero mean across their elements, their scalar product * is proportional to their covariance. This illustrates that this definition is convenient to describe useful relationships between vectors. The L2norm of a vector is defined based on its product with itself, i.e.
A matrix is a rectangular array of numbers: In this case, X is a matrix of dimension m x n , i.e., it has m rows and n columns. It can be considered as an ordered list of vectors, either along columns . j = ( X 1j … X ij ⋯ X mj ) (where j can be between 1 and n ) or rows . = ( X 1 … X ij … X in ) (where i can be between 1 and m ). In the context of matrices, different vectors can have the same number of elements but still different dimension (i.e., they can be row of column vectors, i.e. they can potentially fill a row or a column of a matrix). For column vectors we indicate their dimensions by ( n × 1), while for row vectors it is (1 × n ). Thus, in this way vectors can be considered as matrices with one singleton dimension. The transformation of turning a row vector into a column vector or vice versa is called transposition. For matrices, a transposition "flips over" rows and columns, such that the i -th column of a matrix becomes its transpose's i -th row, and vice versa.
As for vectors, for matrices addition and subtraction are defined element-by-element. The definition of matrix multiplication is less intuitive at first glance. For matrices X ( m x n ) and Y ( n x p ): This only works if the number of columns of X and the number of rows of Y are equal (in this case equal to n ). The resulting matrix P of this matrix multiplication has m rows (the number of rows of X ) and p columns (the number of columns of Y ), i.e. P has dimension ( m x p ).
This expression becomes more intuitive if we consider every element of * as the scalar product (as defined above) between two vectors, namely rows of X and columns of Y . More precisely, every element ( i, j ) of the matrix product = * is the product of the i th row of X and the j th column of Y : As in the example above, if the appropriate rows and columns of X and Y are zero mean then the elements of this matrix product are proportional to the covariances of all combinations of these rows and columns. This is a compact way to describe the relationships among patterns in different matrices.
From the previous expression we can also see that elements of a particular column of P are the result of the multiplication of the same column of Y with successive rows of X . Another interpretation is that every column of P is a linear combination (i.e. weighted sum) of columns of X , with the weightings from one column of Y . It takes a moment to get one's head around this. Similarly, every row of P is constructed with one row of X and successive columns of Y . Thus, the rows of P are linear combinations of the rows of Y . In matrix and vector multiplication the symbol " * " is often omitted. In the following, we will use the " * " for clarity except where we refer to equations in the main part of this manuscript.
For example, the EEG/MEG forward solution is given by = , and since s is a vector (i.e., a matrix with only one dimension), this results in ( Eq. (2 ), with m : number of sensors, n : number of sources) Thus, the predicted data for a given source distribution s are a linear combination of the columns of the leadfield matrix L , weighted by the source strengths associated with each column. This is not surprising, as the columns of L are the topographies for individual sources with unit strength.
The estimate of the sources ̂ = ( Eq. (3 )) with an inverse operator K ( n x m ) can be interpreted in a similar manner: ̂ is a linear combination of the columns of K , with the weightings given by the data.
A central concept of our manuscript is the resolution matrix R , which is the product of an inverse operator matrix K ( n x m ) and the leadfield L ( m x n ). Thus: R has the dimension ( n x n ), i.e., it is a square matrix with both dimensions equal to the number of sources. It provides a fundamental relationship between true and estimated sources ( Eq. (4 ), but we ignore noise here for simplicity):

̂ =
The point-spread function for a source i is the response of our source estimate to a point source with unit strength. This corresponds to a source vector that contains only zeros except at the element corresponding to the source of interest i , i.e. = ( 0 … 1 … 0 ) . Thus: = * * = * According to the definition of matrix multiplication and its possible interpretations described above, this results in one column vector, which is the sum of columns of R weighted with the elements of . Because only the i-th element of is non-zero, the resulting PSF is the i th column of R .
The cross-talk function of a source i describes how all sources in s potentially leak into the estimate of source i . This estimate is given by the i-th element of ̂ . According to ̂ = , this estimate is calculated by multiplying the i-th row of R with s , i.e., ̂ = . * . The elements of this row . therefore tell us how much the different elements of s can potentially affect the estimate of the i-th source. . is thus the cross-talk function for source i .
In summary, the resolution matrix R contains point-spread and crosstalk functions as its columns and rows, respectively.

Matrix inversion
For scalar numbers, taking the inverse is straightforward: multiplying the original number and its inverse with each other yields 1 (e.g., 3 * 1/3 = 1). The equivalent of the scalar number 1 for matrices is the identity matrix, which consists of ones along its diagonal but is zero everywhere else, e.g.,

)
Multiplying this matrix to any other suitable matrix M will produce the same matrix M . In analogy to the inverse of scalar numbers, the inverse of a matrix M is defined as a matrix which -multiplied to the matrix M -yields the identity matrix. For example, for the diagonal matrix ( 2 0 0 3 ) , the inverse matrix is ) .
However, is there an inverse matrix K for the matrix = ( ) ∶= ( 1 0 0 1 ) Thus, the sums K 11 + K 12 and K 21 + K 22 must be both non-zero and zero -which is not possible. In this case, an inverse matrix does not exist. The reason is that the rows and columns of M are multiples of each other and do not pose independent constraints on the elements of the inverse matrix (they are "linearly dependent "). The maximum number of linearly independent rows and columns of a matrix is also called its "rank ". A matrix is only invertible when it is square and has "full rank ", i.e. when the rank is equal to its maximum dimension.
So far we have only considered square matrices. From the above it should be clear that for rectangular matrices (matrices with unequal numbers of rows and columns) an inverse as defined above cannot exist. For example, if there are fewer rows than columns, then the rows may all be linearly independent but the columns cannot be. The number of degrees of freedom that we have to create linearly independent vectors across rows and columns is the minimum of the two dimensions. Thus, if these numbers are not equal, one dimension will always loose out. This is the case for the EEG/MEG inverse problem, where the leadfield commonly has many more columns (sources) than rows (sensors). For this reason, we cannot invert Eq. (2 ) in a straightforward manner to get our true sources. Thus, we have to look for other principled ways to at least get close to our goal. For cases where a unique inverse matrix does not exist, one can still compute so called "pseudo-inverses ", which share at least some properties with the inverse as defined above. One can include additional constraints based on a priori knowledge about the solutions. Or one can reformulate the problem in a way that leads to an invertible matrix, e.g., by focusing only on certain features of the solution. A more detailed description of these approaches is beyond the scope of this appendix. The main manuscript deals with some of the most popular approaches to this problem in the context of EEG/MEG source estimation.