Neurovascular coupling (NVC) provides important insights into the intricate activity of brain functioning and may aid in the early diagnosis of brain diseases. Emerging evidences have shown that NVC could be assessed by the coupling between electroencephalography (EEG) and functional near-infrared spectroscopy (fNIRS). However, this endeavor presents significant challenges due to the absence of standardized methodologies and reliable techniques for coupling analysis of these two modalities. In this study, we introduced a novel method, i.e., the collaborative multi-output variational Gaussian process convergent cross-mapping (CMVGP-CCM) approach to advance coupling analysis of EEG and fNIRS. To validate the robustness and reliability of the CMVGP-CCM method, we conducted extensive experiments using chaotic time series models with varying noise levels, sequence lengths, and causal driving strengths. In addition, we employed the CMVGP-CCM method to explore the NVC between EEG and fNIRS signals collected from 26 healthy participants using a working memory (WM) task. Results revealed a significant causal effect of EEG signals, particularly the delta, theta, and alpha frequency bands, on the fNIRS signals during WM. This influence was notably observed in the frontal lobe, and its strength exhibited a decline as cognitive demands increased. This study illuminates the complex connections between brain electrical activity and cerebral blood flow, offering new insights into the underlying NVC mechanisms of WM.