## 1 Introduction  
From the neuroscience perspective, each physiological or cognitive process will produce a particular pattern of electrical interactions linking neurons from different brain regions. Therefore, the electrical response related to the interactions between neuron assemblies allows studying the brain function. Under such a perspective, the neuroscientists have found brain electrical activity (BEA) patterns associated with motor functions, cognitive processes, and neuropathologies [1]. However, the wide range of possible mental conditions hampers the BEA analysis and poses very challenging tasks [2]. In particular, capturing BEA by placing a set of electrodes over the scalp, known as electroencephalography (EEG), gathers and amplifies currents reflected in the brain cortex from all the possible brain sources, yielding a mixture of latent activity sources at each channel. For discovering such a latent activity, the literature considers different types of EEG analyses ranging from time and spectral domain processing [3, 4], through connectivity measures between channels [5], to the complex network analysis [6]. Due to the non-stationarity behavior of EEG data, the classical temporal analyses cannot decode differences among several mental states or conditions. Besides, several studies demonstrate that the variations in the EEG oscillatory patterns play a fundamental role in the maintenance of brain functions and the identification of different neural conditions [7]. Hence, spectral decomposition approaches look for relevant information at the brain rhythms (namely alpha, beta, theta, delta, and gamma). For instance, performing working memory or creative tasks evokes discriminative oscillatory patterns [8; 9]. In addition, neuropathologies, such as Alzheimer's disease, cause abnormal cortical neural synchronization at resting-state rhythms [10]. The general spectral analysis procedure consists of a channel-wise frequency band splitting, followed by an identification of relevant time intervals, usually supervised by a specialist. The subsequent stages compute descriptors from the time-frequency representation. The power spectral density (PSD) and spectral entropy are among the most considered descriptors due to their straightforward interpretation [11; 12]. Although EEG frequency analysis has proven to extract useful information for brain function understanding, the channel-wise approach still lacks the interpretability of several cognitive processes because each channel only holds a reflected version of BEA from a neural assembly [13]. Besides, the low spatial resolution of EEG restricts the extracted information about the behavior of some brain regions involved in high-level cognitive or physiological processes [4]. In recent years, connectivity analysis techniques account for channel interdependencies to enhance EEG spatial resolution through functional relationships captured in the BEA. Metrics from several domains attempt to quantify different properties of the BEA. Specifically, the coherence (COH), phase value index (PVI), and phase-locking value (PLV) capture pair-wise channel dependencies in the frequency domain. Moreover, the latter two gained attention due to their non-linear capability for unraveling latent connectivity patterns, which have proven valuable for applications, such as motor imagery and emotion recognition [14; 15]. Nevertheless, the selection of the connectivity measure is not entirely straightforward. Their dependence on an estimated cross-spectrum and subject data variability yield to low generalization capability in a wide range of scenarios [16; 17]. The above issues demand reliable brain connectivity approaches to automatically identify the relevant spectral information from the inherent EEG uncertainty. In this regard, a probabilistic modeling framework, such as Gaussian Processes (GPs), along with an appropriate covariance function possesses the capability for characterizing the latent processes from BEA data [18; 19].