11institutetext: Centre for Digital Music, Queen Mary University of London
11email: {max.graf<EMAIL_ADDRESS>
# Combining Vision and EMG-Based Hand Tracking for Extended Reality Musical
Instruments

Max Graf 11 Mathieu Barthet 11

###### Abstract

Hand tracking is a critical component of natural user interactions in extended
reality (XR) environments, including extended reality musical instruments
(XRMIs). However, self-occlusion remains a significant challenge for vision-
based hand tracking systems, leading to inaccurate results and degraded user
experiences. In this paper, we propose a multimodal hand tracking system that
combines vision-based hand tracking with surface electromyography (sEMG) data
for finger joint angle estimation. We validate the effectiveness of our system
through a series of hand pose tasks designed to cover a wide range of
gestures, including those prone to self-occlusion. By comparing the
performance of our multimodal system to a baseline vision-based tracking
method, we demonstrate that our multimodal approach significantly improves
tracking accuracy for several finger joints prone to self-occlusion. These
findings suggest that our system has the potential to enhance XR experiences
by providing more accurate and robust hand tracking, even in the presence of
self-occlusion.

###### keywords:

Extended reality, extended reality musical instruments, hand tracking, surface
electromyography, deep learning

##  1 Introduction

Extended reality (XR) is an umbrella term encompassing virtual, augmented and
mixed reality (VR/AR/MR). In recent years, the increased popularity of XR
technology has seen the establishment of extended reality musical instruments
(XRMIs) as a research field [22]. Milgram et al. described the reality-
virtuality continuum [20], along which digital applications can be placed. It
stretches from real-world environments to fully virtual environments. Head-
mounted XR devices bridge this continuum. They are capable of rendering three-
dimensional imagery onto screens, removing the necessity for separate monitors
or mobile displays, blending the real and virtual worlds together. The rapid
development of XR technologies has opened up new possibilities for musical
creation, performance, and interaction, with the emergence of various XR-based
musical instruments and applications. Many XRMIs follow an embodied
interaction paradigm. These instruments offer novel opportunities for artists
to experiment with embodied interaction techniques, spatial sound design, and
immersive performances, thus expanding the boundaries of traditional music
making. XRMIs fall within the larger category of digital musical instruments
(DMIs). [26] suggest that the control of DMIs can be made intimate (personal
and familiar) by using appropriate control metaphors, low latency action-to-
sound, and continuous gesture recognition. This study is part of a larger
project that aims to support control intimacy in XRMIs.

Based on the current state of XR technology and prior work in XRMIs [3, 8], we
highlight that gesture sensing errors on XR devices are a bottlenek for
intimate musical control. Head-mounted XR devices (HMDs) rely on a set of
sensors to record data and provide embodied control interfaces for users,
e.g., head-tracking, hand tracking, and body pose detection. The transduction
of these real-world sensor data to digital representations depends on
computational methods. In this work we focus on the problem of hand tracking,
more specifically, accurate tracking of finger joints. Hand-tracking
algorithms often use visual information from camera sensors in conjunction
with machine learning techniques, for example, in the Oculus Quest 2 device
[10]. The accuracy of vision-based hand-tracking algorithms may be high [25],
but current recognition rates do not reach 100%. Self-occlusion - the
occlusion of finger joints by other parts of the hand - as well as challenging
lighting situations lead to failure cases in vision-based tracking systems.
Such error cases may produce instances of jitter, tracking loss, or glitches
in the virtual representation of the hands, which can have detrimental effects
on the usability and user experience in XRMIs, as shown in a previous study
[8].

This work aims to address such sensing-related issues through the use of
surface electromyography (sEMG) sensors. EMG sensors measure the electrical
potential produced during muscle contractions in the body. Surface
electromyograms can be obtained through electrodes that are positioned on the
surface of the skin, above muscle tissue regions. We present an investigation
into the potential of sEMG sensors and deep learning models to enhance hand-
tracking accuracy in XRMIs. Our approach combines sEMG data and vision-based
tracking methods to address sensing-related issues commonly encountered in
XRMI performance. Thereby, we aim to improve the tracking accuracy and
responsiveness of XR musical instruments, especially in situations where
vision-based tracking falls short.

The scope of this paper is limited to the exploration of sEMG and deep
learning techniques for hand-tracking in XR musical instruments. While our
findings may have broader applications in other areas of XR interaction, the
primary focus is on the improvement of XRMI design and user experience.
Through our work, we aim to contribute to the ongoing development of more
accurate, intuitive, and expressive extended reality musical instruments.

##  2 Background

The development of XRMIs has attracted growing interest as VR, AR, and MR
technologies continue to advance. Early studies in this domain focused on the
design and evaluation of virtual interfaces for musical performance and
interaction [18, 22, 7, 21]. Several works investigated user experience [6,
5], interaction techniques [2] and collaborative music making [19, 9]. More
recent studies have explored the creation of novel instruments and control
schemes [5, 3, 4, 8]. While there is no gold standard for XRMI design, many
XRMIs rely on hand-tracking to facilitate embodied interaction with the
instrument [21, 7, 4, 8]. Various vision-based tracking methods are employed,
including depth-sensing cameras [21, 7, 4], and machine learning-based
approaches [10]. Despite the progress in hand tracking research, limitations
such as occlusion, lighting issues, and computational complexity continue to
pose challenges for hand-controlled XRMI applications.

Surface electromyography (sEMG) has emerged as a promising alternative to
vision-based tracking methods for capturing user input in various
applications, including XR. Several studies have explored the use of sEMG data
and deep learning architectures for hand gesture recognition [17, 16, 1, 15].
These works have reported promising results, highlighting the potential of
employing sEMG and deep learning models for precise finger movement
estimation. However, some of these works depend on complex tracking setups [1]
or leverage low-resolution tracking data for training [15].

A notable limitation of these studies is the lack of shared training data and
code, hindering the reproducibility and comparability of the results across
different research efforts. Several datasets on the topic of finger joint
angle estimation through sEMG data have been published. However, they either
require specialised sEMG measuring equipment [12, 14, 13], making the
reproduction of results an expensive endeavour, or introduce temporal biases
into the dataset due to lack of synchronisation during the recording procedure
[11]. The absence of implementation details ([17, 16, 1]) makes it difficult
for other researchers to build upon these works, potentially slowing down the
progress in the field.

The potential benefits of integrating sEMG data and deep learning models with
vision-based tracking methods has not been thoroughly investigated in the
context of XRMIs. While machine learning has found its way into the NIME
community [23], the use of machine learning approaches to improve XRMI control
remains an underexplored area. This study aims to develop a multimodal hand
tracking approach that leverages sEMG data, deep learning models, and vision-
based tracking techniques. We share the training data and
code111https://github.com/maxgraf96/sEMG-myo-unity, fostering further research
and innovation in the domain of sEMG-based neural interfaces.

##  3 Deep Learning Model for Finger Joint Angle Estimation

We have developed a software pipeline for data collection, feature extraction
and modelling of sEMG data. We focus on eight finger joints that are prone to
self-occlusion: the metacarpophalangeal and proximal interphalangeal joints of
the index, middle, ring and pinky fingers. Specifically, we want to model the
rotations of the finger bones connected to these joints relative to the hand.
The thumb is excluded from our investigation. Modelling thumb rotations with
sEMG data is a hard problem, since the majority of muscles related to thumb
movements are located in the hand, rather than the forearm. With that in mind,
the goal of the model is to estimate the eight finger joint angles from a
window of sEMG data.

###  3.1 Data Collection

We collect surface EMG signal measurements using the Thalmic Labs Myo
armband222https://xinreality.com/wiki/Myo and vision-based hand tracking data
from the Oculus Quest 2 XR headset. Both devices are employed simultaneously
to capture the muscle activity and finger joint rotations, respectively. This
approach allows users to capture data without the need for external tracking
devices.

The Myo armband is a non-invasive wearable device that features eight sEMG
sensors. The armband is worn on the forearm, with the sensors evenly
distributed around the circumference of the arm, allowing it to capture the
activity of the forearm muscles during finger movements. We obtain sEMG data
from the Myo armband using the PyoMyo Python framework [24], extracting
rectified and smoothed signals at a sampling frequency of 50Hz. The Oculus
Quest 2 XR headset is equipped with four monochrome cameras that provide a
wide field of view, enabling it to capture hand positions and movements. The
built-in hand tracking algorithm [10] processes the camera data and estimates
the 3D rotations of the user’s hand joints in real time. In our system we
sample the hand joint rotations from the XR device at 50 Hz and synchronize
them with the sEMG data from the Myo armband.

One researcher recorded hand gestures and movements in a controlled
environment, with a focus on gestures relevant to XRMI interaction. This study
should be seen as a proof-of-concept for our methodology. Hence, for this
study, we focused on data from the right hand only. The gestures included
various finger flexions and extensions, as well as combinations of multiple
finger movements. They were performed at different speeds, forces, and
orientations. We conducted three data collection sessions across three days to
account for the natural variability in sEMG readings, ensuring a more robust
dataset. During the data collection process, the XR headset was strategically
positioned in diverse locations and orientations to minimise self-occlusion of
the hand. Data collection sessions lasted between ten and fifteen minutes,
resulting in a substantial amount of synchronized sEMG and hand tracking data.

###  3.2 Feature Extraction

Figure 1 shows the data flow in our pipeline. The selection of features was
informed by both the literature and a series of experiments.We use Python to
extract both time domain and frequency domain features from the sEMG data. The
pipeline takes 2D windows of sEMG samples, with $N$ number of samples and $C$
channels. We then compute the following features per channel: in the time
domain, mean absolute value (MAV), root mean square (RMS), and variance (VAR);
in the frequency domain, median frequency bin (MDF), mean frequency bin (MNF),
and peak frequency bin (PF). Additionally, wavelet coefficients at the fourth
level are extracted to provide further information about the signal’s
characteristics, as reported in [1]. Wavelet analysis provides a multi-
resolution representation of the sEMG signal, capturing both the time and
frequency characteristics of the data. This comprehensive feature
representation aims to capture essential information from the sEMG signals
both locally and globally, enabling holistic representation of the data.

Figure 1: Data preprocessing, feature extraction, and model pipeline

###  3.3 Model Architecture

The model architecture is a combination of a Long Short-Term Memory (LSTM)
network and multiple separate sets of fully connected (FC) layers, aiming to
capture both general trends in the sEMG signal data provided by the features
and high-frequency characteristics of the signal. Formally, our model learns a
mapping

$F:\mathbb{R}^{N\times C}\rightarrow\mathbb{R}^{M}$ (1)

where $N$ denotes the number of sEMG samples, $C$ represents the number of
sEMG channels and $M$ denotes the number of predicted joint angle values at
every time step. More accurately, we learn a mapping

$F(\mathbf{s})=\phi_{\text{final}}\Bigl{(}\phi_{\text{lstm}}(\mathbf{s}_{\hat{N}})\oplus\phi_{\text{feat}}(\psi_{\text{time-
freq}}(\mathbf{s}_{N}))\oplus\phi_{\text{wav}}(\psi_{\text{wavelet}}(\mathbf{s}_{N}))\oplus\phi_{\text{filt}}(\mathbf{s}_{\hat{N}})\Bigr{)}$
(2)

where

  * •

$F(\mathbf{s})$ represents the mapping function that takes $N\times C$ sEMG
samples ($\mathbf{s}$) and outputs $M$ finger joint angles.

  * •

$\mathbf{s}_{N}$ denotes all $N\times C$ sEMG samples.

  * •

$\mathbf{s}_{\hat{N}}$ denotes the last $\hat{N}\times C$ sEMG samples in the
data point.

  * •

$\phi_{\text{lstm}}$ represent the LSTM layers, and $\phi_{\text{feat}}$,
$\phi_{\text{wav}}$ and $\phi_{\text{filt}}$ represent the fully connected
layers processing time/frequency domain features, wavelet features and
filtered EMG data respectively.

  * •

$\phi_{\text{final}}$ denotes the final set of fully connected layers.

  * •

$\psi_{\text{time-freq}}$ and $\psi_{\text{wavelet}}$ denote the feature
extraction functions for time-frequency domain features and wavelet features,
respectively.

  * •

$\oplus$ represents the tensor concatenation operation.

The LSTM network, a type of recurrent neural network, operates on the subset
$\mathbf{s}_{\hat{N}}$, which contains the last $\hat{N}$ samples of the EMG
data, capturing the temporal dependencies within the most recent portion of
the sEMG signal. LSTMs can effectively learn medium-to-long-range dependencies
in time-series data and retain information across multiple time steps.
Additionally, $\mathbf{s}_{\hat{N}}$ is fed into a separate fully connected
layer, which was empirically found to improve the model’s performance. The
time and frequency domain features are computed over all $N$ samples and
processed by a set of fully connected layers. These layers are designed to
extract higher-level representations from the sEMG features, capturing general
patterns and trends in the data. The wavelet features are fed into a separate
set of fully connected layers, allowing the model to learn distinct patterns
associated with the wavelet coefficients. This additional information can help
the model to better discriminate between different types of hand movements and
gestures.

The outputs of the LSTM layers and the three sets of fully connected layers
(time-frequency domain features, wavelet features, and ${\hat{N}}$ sEMG
samples) are concatenated and passed to a final set of fully connected layers.
This combination of network components aims to capture a comprehensive
representation of the sEMG signal, taking into account both general trends and
high-frequency changes. The final output of the model is an estimation of the
eight finger joint angles described above.

###  3.4 Model Implementation and Training

We selected $N=150$, which gives a sampling window size of three seconds.
Related works use window sizes of five seconds [17, 16]. During our
experiments with the model architecture, we found that a smaller window size
of 150 samples did not reduce the quality of the predicted finger joint
angles, while yielding a performance gain in the data processing pipeline.
$\hat{N}=50$ was selected as a trade-off between LSTM accuracy and
performance. Higher values for $\hat{N}$ produced slightly better results, but
incurred a performance degradation, slowing down the model at inference time.

The collected data comprises approximately 80000 sEMG samples with
corresponding finger joint angle measurements. The sEMG samples were separated
into training and validation sets using a 90/10 ratio. Our deep learning model
is built using the PyTorch framework. The model was trained on a single NVidia
RTX 2080 Ti GPU for approximately 500000 steps with a batch size of 256. The
mean squared error function was employed as a loss metric. We applied an
exponentially decreasing learning rate, starting at 0.003 and reducing to
0.0003 over the first 10000 steps. During training, we tracked the mean
average difference between the predicted joint angles and the ground truth
angles for both training and validation sets. Our criterion for stopping the
training procedure was the moment of obtaining a mean joint angle difference
value of less than 1° across the validation set.

##  4 Multimodal XR Hand Tracking with sEMG and Vision-Based Tracking

In this section, we present our approach to multimodal XR hand tracking by
combining sEMG and vision-based tracking techniques. In our system, the
vision-based tracking data provides information about the overall hand
position and orientation in 3D space, while the sEMG-based model produces
granular information about individual finger joint rotations. The sEMG data
are continuously sampled, preprocessed, and passed to the trained deep
learning model to estimate the eight finger joint angles. The hand position
and orientation data are combined with the estimated finger joint angles to
generate a complete hand pose representation.

The deep learning model is optimized for real-time performance, ensuring that
the sEMG data can be processed with minimal latency. Our system operates at
50Hz, which incurs a latency of 20ms. For this work, data processing and model
inference took place on a consumer notebook. The notebook concurrently runs a
python server, responsible for sEMG data aggregation, preprocessing and model
inference, and a 3D XR environment on the Unity platform. At every time step,
the estimated finger joint angles are transferred from Python to Unity through
the low-latency ZeroMQ framework333https://zeromq.org/. A video demonstration
of our system is available
online444https://www.youtube.com/watch?v=ivl2g2t2oaI. It shows a side-by-side
comparison of the vision-based tracking system and our multimodal approach.

###  4.1 Evaluation

To validate the effectiveness of our multimodal hand tracking system, we
conducted an experimental evaluation, which compared the multimodal tracking
to the baseline vision-based hand tracking system provided by the XR device.
We simultaneously collected tracking data from the vision-based tracking
system and the multimodal tracking system. A Leap Motion sensor was used to
acquire ground truth labels for finger joint angles. The Leap Motion is a
high-precision vision-based hand tracking device that captures finger joint
angles and hand position in 3D space. The ground truth labels were then
compared to the results from both the vision-based and the multimodal tracking
system.

The experimental setup included a series of hand pose tasks, designed to cover
a wide range of hand movements, including gestures prone to occlusion. The
tasks were selected with regard to their utility in playing a keyboard-
inspired XRMI. The tasks were performed while wearing the Oculus Quest 2
headset and the Myo armband. The Leap Motion sensor was placed on a table to
record ground truth data. We recorded tasks under two conditions: 1) Full view
of the hand - here, the XR headset was positioned at a 50cm distance, 45°
above the hand to ensure optimal visual tracking conditions. 2) Self-occlusion
of the hand - in this condition, the distance was kept identical, but the
angle of the XR headset was lowered, such that the back of the hand occluded
the fingers. Figure 2 illustrates the six hand pose tasks devised: (i)
extending all fingers and making a fist; individual flexion and extension of
the (ii) index, (iii) middle, (iv) ring and (v) pinky fingers; (vi) sequential
flexion and extension of pinky, ring, middle and index fingers (similar to the
gesture of drumming on a table while waiting for something). Each task
involves the execution of the gesture at three different speeds: slow, over
approximately two seconds, moderate (one second), and fast (half a second). To
account for variability in the sEMG measurements, all tasks were executed
three times, over two days, under identical lighting conditions.

To measure the degree of finger occlusion, we integrated a ray-casting system
with the XR application. We cast rays from the XR headset’s 3D position to the
eight finger bones whose rotations we measured every time a sample was taken.
Rays intersecting other parts of the hand, e.g., the back of the hand, were
used to mark the respective finger bones as occluded. This allowed us to
quantify the level of occlusion per finger for every recording.

Figure 2: Finger movements in the six hand pose tasks

###  4.2 Analysis Methods & Results

To evaluate the performance of the vision-based and multimodal tracking
systems across every task, we obtained matrices of difference values between
the estimated joint angles and the ground truth angles for both systems at
every time step. The matrices were aggregated across the three sessions. We
assessed the normality of the difference matrices using the Shapiro-Wilk test
for each of the six hand pose tasks. We then applied the Wilcoxon signed-rank
test to see whether there was a significant difference between the results
produced by the vision-based and multimodal tracking systems.

Across tasks, the results of the Shapiro-Wilk test showed p-values $<.001$,
indicating that the data in the difference value matrices did not fit a normal
distribution. Therefore, we proceeded with the non-parametric Wilcoxon signed-
rank test for further analysis. Figure 3 shows the results obtained from the
mean joint angle differences across all finger joints per task, for both
occlusion conditions (full view and occluded). Under the occluded condition,
our model produces significantly lower deviations from the ground truth data
across all tasks, compared to the vision-based tracking system. On average, it
improves the finger joint angle tracking accuracy by five to 15 degrees across
all fingers.

Table 1 shows the results obtained from the Wilcoxon signed-rank tests,
aggregated across all eight tracked joints per task for both conditions. The
p-values indicate significant differences between the difference value
matrices. Additionally, the table lists the average finger occlusion results
obtained through the raycasting occlusion measure described above. The
occlusion results describe the mean portion per task, in which the fingers
were occluded by another part of the hand.

Figure 3: Average deviation in degrees between the finger joint angles generated by the vision-based (V) and multimodal (M) tracking systems and the ground truth data for each task. Error bars show standard deviation; three asterisks indicate a significant difference between the V and M values ($p<.001$) Table 1: P-values and average occlusion measurement results across tasks under both conditions | Full view | Occluded
---|---|---
Tasks | P-value | Occlusion (%) | P-value | Occlusion (%)
(i) | 1.0 | 16.78 | <.001 | 93.00
(ii) | 1.0 | 2.65 | <.001 | 70.58
(iii) | 1.0 | 15.25 | <.001 | 63.11
(iv) | 1.0 | 16.78 | <.001 | 53.02
(v) | 1.0 | 27.94 | <.001 | 78.82
(vi) | 1.0 | 10.36 | <.001 | 59.48

##  5 Discussion

The results of the evaluation showed that the multimodal hand tracking system
outperformed the pure vision-based hand tracking system across all tasks under
the occluded condition. These findings support our hypothesis that the
integration of sEMG-based finger joint angle estimation can help overcome
occlusion-related limitations in vision-based hand tracking, resulting in more
accurate and reliable XRMI interactions. Under the full view condition, the
vision-based hand tracking produced fewer errors in all tasks. This was
expected, as the vision-based hand tracking system operates optimally under
full view of the hand.

Despite the promising results, our study has several limitations. Due to the
nature of sEMG data, the tracking performance of the multimodal approach is
unlikely to extend to other users without fine-tuning the deep learning model.
Surface EMG signals differ substantially between individuals and can be
influenced by factors such as muscle fatigue, electrode placement, and
individual anatomical differences. It will be valuable to investigate the
system’s performance across different users and under varying conditions. The
identification of sEMG data representations that allow for generalisation
under consideration of these factors without requiring extensive amounts of
data is still an ongoing research topic. However, our work allows XR users
with access to sEMG devices to train their own models using our pipeline and
code.
The performance of our multimodal hand tracking system was evaluated using a
single type of XR headset and sEMG armband. Future research should explore
more complex occlusion scenarios, as well as test the system’s performance
across different hardware setups and sEMG devices, to better understand the
generalisability of our findings.
With that in mind, we see numerous avenues for further research. The
integration of additional tracking modalities, such as depth sensing or
inertial measurement units (IMUs), could further enhance the robustness and
accuracy of the multimodal hand tracking system by enabling stronger
representations of the underlying data. A future study will explore the impact
of our multimodal hand tracking system on usability, user experience and task
performance in XRMI interactions, and provide insights into the practical
implications of our findings. By conducting user studies with tasks that
require precise hand movements and are susceptible to occlusion, the benefits
of our system for real-world applications could be better understood.
Our study provides evidence that the combination of vision-based tracking and
sEMG-based finger joint angle estimation can effectively address occlusion
issues in hand tracking for XRMI interactions. The findings suggest that the
multimodal hand tracking system has the potential to enhance user experiences
and enable more immersive and natural interactions in virtual environments.

##  6 Conclusion

In this paper, we introduced a multimodal hand tracking system designed to
address occlusion issues in XRMI interactions by combining vision-based
tracking with sEMG-based finger joint angle estimation. The goal of this study
was to demonstrate the potential of our proposed system to improve hand
tracking accuracy and robustness, even when the hand is partially occluded.

While our results show promise, the experimental setup was relatively simple,
and further research should explore more complex scenarios and investigate the
system’s performance across different hardware and user conditions. Future
work could also integrate additional tracking modalities and machine learning
techniques to enhance the robustness and accuracy of the system.

Our multimodal hand tracking system demonstrates the potential to improve XRMI
interactions by addressing occlusion issues in vision-based hand tracking. As
XR technologies continue to evolve, the integration of complementary tracking
modalities, such as sEMG and vision-based tracking, will likely play a crucial
role in enhancing user experiences and enabling more immersive and natural
interactions in virtual environments.

## References

  * [1] Cries Avian, Setya Widyawan Prakosa, Muhamad Faisal and Jenq-Shiou Leu  “Estimating Finger Joint Angles on Surface EMG Using Manifold Learning and Long Short-Term Memory with Attention Mechanism”  In _Biomedical Signal Processing and Control_ 71, 2022, pp. 103099  DOI: 10.1016/j.bspc.2021.103099
  * [2] Florent Berthaut  “3D Interaction Techniques for Musical Expression”  In _Journal of New Music Research_ 49.1, 2020, pp. 60–72  DOI: 10.1080/09298215.2019.1706584
  * [3] Sam Bilbow  “Developing Multisensory Augmented Reality As A Medium For Computational Artists”  In _Proceedings of the Fifteenth International Conference on Tangible, Embedded, and Embodied Interaction_ New York, NY, USA: Association for Computing Machinery, 2021, pp. 1–7  URL: https://doi.org/10.1145/3430524.3443690
  * [4] Sam Bilbow  “Evaluating Polaris~ - An Audiovisual Augmented Reality Experience Built on Open-Source Hardware and Software”  In _NIME 2022_ , 2022
  * [5] Anıl Çamcı, Matias Vilaplana and Ruth Wang  “Exploring the Affordances of VR for Musical Interaction Design with VIMEs”  In _Proceedings of the International Conference on New Interfaces for Musical Expression_ , 2020, pp. 121–126  URL: https://zenodo.org/record/4813268
  * [6] Thomas Deacon, Tony Stockman and Mathieu Barthet  “User Experience in an Interactive Music Virtual Reality System: An Exploratory Study”  In _Bridging People and Sound_ 10525 Cham: Springer International Publishing, 2017, pp. 192–216  DOI: 10.1007/978-3-319-67738-5˙12
  * [7] J. Fillwalk  “ChromaChord: A Virtual Musical Instrument”  In _2015 IEEE Symposium on 3D User Interfaces (3DUI)_ , 2015, pp. 201–202  DOI: 10.1109/3DUI.2015.7131770
  * [8] Max Graf and Mathieu Barthet  “Mixed Reality Musical Interface: Exploring Ergonomics and Adaptive Hand Pose Recognition for Gestural Control”  In _International Conference on New Interfaces for Musical Expression_ , 2022  DOI: 10.21428/92fbeb44.56ba9b93
  * [9] Rob Hamilton and Chris Platz  “Gesture-Based Collaborative Virtual Reality Performance in Carillon”  In _Proceedings of the 2016 International Computer Music Conference_ , 2016, pp. 5
  * [10] Shangchen Han et al.  “MEgATrack: Monochrome Egocentric Articulated Hand-Tracking for Virtual Reality”  In _ACM Transactions on Graphics_ 39.4, 2020, pp. 87:87:1–87:87:13  DOI: 10.1145/3386569.3392452
  * [11] Xuhui Hu et al.  “Finger Movement Recognition via High-Density Electromyography of Intrinsic and Extrinsic Hand Muscles”  In _Scientific Data_ 9.1 Nature Publishing Group, 2022, pp. 373  DOI: 10.1038/s41597-022-01484-2
  * [12] Néstor J. Jarque-Bou et al.  “A Calibrated Database of Kinematics and EMG of the Forearm and Hand during Activities of Daily Living”  In _Scientific Data_ 6.1 Nature Publishing Group, 2019, pp. 270  DOI: 10.1038/s41597-019-0285-1
  * [13] Ning Jiang  “Gesture Recognition and Biometrics Electromyography (GRABMyo) Dataset”  IEEE, 2022  URL: https://ieee-dataport.org/documents/gesture-recognition-and-biometrics-electromyography-grabmyo-dataset
  * [14] Piotr Kaczmarek, Tomasz Mańkowski and Jakub Tomczyński  “putEMG—A Surface Electromyography Hand Gesture Recognition Dataset”  In _Sensors_ 19.16 Multidisciplinary Digital Publishing Institute, 2019, pp. 3548  DOI: 10.3390/s19163548
  * [15] Hyunin Lee, Dongwook Kim and Yong-Lae Park  “Explainable Deep Learning Model for EMG-Based Finger Angle Estimation Using Attention”  In _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ 30, 2022, pp. 1877–1886  DOI: 10.1109/TNSRE.2022.3188275
  * [16] Yang Liu, Chengdong Lin and Zhenjiang Li  “WR-Hand: Wearable Armband Can Track User’s Hand”  In _Proceedings of the ACM on Interactive, Mobile, Wearable and Ubiquitous Technologies_ 5.3, 2021, pp. 118:1–118:27  DOI: 10.1145/3478112
  * [17] Yilin Liu, Shijia Zhang and Mahanth Gowda  “NeuroPose: 3D Hand Pose Tracking Using EMG Wearables”  In _Proceedings of the Web Conference 2021_ , WWW ’21  New York, NY, USA: Association for Computing Machinery, 2021, pp. 1471–1482  DOI: 10.1145/3442381.3449890
  * [18] Teemu Mäki-Patola, Juha Laitinen, Aki Kanerva and Tapio Takala  “Experiments with Virtual Reality Instruments”  In _Proceedings of the 2005 Conference on New Interfaces for Musical Expression_ , NIME ’05  SGP: National University of Singapore, 2005, pp. 11–16
  * [19] Liang Men and Nick Bryan-Kinns  “LeMo: Supporting Collaborative Music Making in Virtual Reality”  In _IEEE 4th VR Workshop on Sonic Interactions for Virtual Environments (SIVE)_ , 2018  DOI: 10.13140/RG.2.2.27033.93289
  * [20] Paul Milgram, Haruo Takemura, Akira Utsumi and Fumio Kishino  “Augmented Reality: A Class of Displays on the Reality-Virtuality Continuum”  In _Telemanipulator and Telepresence Technologies_ 2351, 1994  DOI: 10.1117/12.197321
  * [21] Alec G. Moore et al.  “Wedge: A Musical Interface for Building and Playing Composition-Appropriate Immersive Environments”  In _2015 IEEE Symposium on 3D User Interfaces (3DUI)_ Arles, France: IEEE, 2015, pp. 205–206  DOI: 10.1109/3DUI.2015.7131772
  * [22] Stefania Serafin et al.  “Virtual Reality Musical Instruments: Guidelines for Multisensory Interaction Design”  In _Proceedings of the Audio Mostly 2016_ Norrköping Sweden: ACM, 2016, pp. 266–271  DOI: 10.1145/2986416.2986431
  * [23] Théo Jourdan and Baptiste Caramiaux  “Machine Learning for Musical Expression: A Systematic Literature Review”, 2023
  * [24] Peter Walkington  “PyoMyo”, 2021  URL: https://github.com/PerlinWarp/pyomyo
  * [25] Frank Weichert, Daniel Bachmann, Bartholomäus Rudak and Denis Fisseler  “Analysis of the Accuracy and Robustness of the Leap Motion Controller”  In _Sensors_ 13.5 Multidisciplinary Digital Publishing Institute, 2013, pp. 6380–6393  DOI: 10.3390/s130506380
  * [26] David Wessel and Matthew Wright  “Problems and Prospects for Intimate Musical Control of Computers”, 2001  DOI: 10.5281/zenodo.1176382


# Non-unitary three-neutrino mixing in the early Universe

Stefano Gariazzo  Pablo Martínez-Miravé  Olga Mena  Sergio Pastor  and Mariam
Tórtola

###### Abstract

Deviations from unitarity in the three-neutrino mixing canonical picture are
expected in many physics scenarios beyond the Standard Model. The mixing of
new heavy neutral leptons with the three light neutrinos would in principle
modify the strength and flavour structure of charged-current and neutral-
current interactions with matter. Non-unitarity effects would therefore have
an impact on the neutrino decoupling processes in the early Universe and on
the value of the effective number of neutrinos, $N_{\rm eff}$. We calculate
the cosmological energy density in the form of radiation with a non-unitary
neutrino mixing matrix, addressing the possible interplay between parameters.
Highly accurate measurements of $N_{\rm eff}$ from forthcoming cosmological
observations can provide independent and complementary limits on the
departures from unitarity. For completeness, we relate the scenario of small
deviations from unitarity to non-standard neutrino interactions and compare
the forecasted constraints to other existing limits in the literature.

##  1 Introduction

In recent years, an unprecedented accuracy has been reached in the extraction
of the neutrino mixing parameters [1, 2, 3]: solar, atmospheric, reactor and
accelerator neutrino oscillation experiments have provided conclusive evidence
for neutrino masses in nature. However, the underlying mechanism for
generating the neutrino masses remains unknown. Although there is a plethora
of possible theoretical frameworks aiming to settle this issue, generally, a
common mechanism consists on the generation of the neutrino masses through the
mediation of heavy neutral leptons (HNLs). Then, the lepton mixing matrix has
to account for the possible admixture between light neutrinos and the HNLs. In
seesaw models [4, 5, 6, 7, 8], the smallness of neutrino masses would indicate
a very high scale for the conventional realisations of this mechanism, leading
to a negligible mixing between light neutrinos and HNLs [9], rendering a
mixing matrix very close to unitary. Nonetheless, in ‘low scale’ realisations
of the seesaw mechanism [10, 11, 12, 13, 14, 15], much more relevant
departures from unitarity can be expected. In this context, two general
scenarios can be envisioned: (i) not-so-heavy neutral leptons, with masses
such that they can be produced in oscillation experiments, or (ii) canonical
HNLs, which are kinematically inaccessible in oscillation experiments. The
former is being scrutinised due to its potential to explain short-baseline
neutrino oscillation anomalies (see e.g. [16, 17, 18, 19, 20]), and its
cosmological effects have been studied extensively in the past [21, 22, 23,
24]. However, when discussing non-unitarity (NU) of the three-neutrino mixing
matrix, we shall consider exclusively the latter. Direct searches for HNLs and
precision flavour observables [25] allow to constrain the light-heavy mixing
for a wide range of masses. These measurements rely on the fact that the
charged-current leptonic interactions are modified. For neutrino experiments
operating at energies below the mass of the HNLs, neutrino neutral-current
interactions are also altered. As a consequence, complementary limits from
both the non-unitarity mixing and the changes in the structure of the
interactions can also be obtained within this context.

In the early Universe, electroweak processes kept neutrinos in thermal contact
with electrons, positrons and photons. Once the collisions became ineffective
at temperatures around the MeV, neutrinos decoupled from the cosmic plasma.
Consequently, any modification of neutrino interactions would change the
picture of neutrino decoupling and alter the resulting background of relic
neutrinos. The role of neutrinos is encoded in the effective number of
relativistic species $N_{\rm eff}$, which quantifies relativistic degrees of
freedom different from photons that contribute to the cosmological radiation
density. Within the canonical picture of neutrino decoupling, $N_{\rm
eff}=3.0440$ [26, 27, 28], in excellent agreement with the most recent
determination from the combination of several cosmological probes, $N_{\rm
eff}=2.99^{+0.34}_{-0.33}$ at 95% confidence level (CL) from Planck 2018 (TT,
TE, EE + lowE + lensing + BAO) [29]. This measurement is expected to improve
significantly with forthcoming data from Cosmic Microwave Background (CMB)
measurements from the Simons Observatory ($\sigma(N_{\rm eff})\simeq
0.05-0.07$) [30] and CMB-S4 ($\sigma(N_{\rm eff})\simeq 0.02-0.03$) [31].

In this work, we address the impact of a non-unitary three-neutrino mixing on
the process of neutrino decoupling and determine the expected complementary
bounds on departures from unitarity that will be obtained by near future CMB
surveys via their measurement of $N_{\rm eff}$. The structure of this
manuscript is as follows. In Section 2, we introduce the formalism and
notation adopted for the study of non-unitarity three-neutrino mixing and in
Section 3 we discuss the main implications for neutrino cosmology. The
obtained results are presented in Sections 4 and 5, where we also discuss a
mapping to the non-standard interactions scenario in case of small departures
from unitarity, as well as the complementarity with other searches. We
conclude in Section 6.

##  2 Non-unitarity in the three-neutrino mixing scenario

Let us denote the $n$-vector representing the neutral lepton states with
definite mass, both light and heavy, by $\upnu_{L}$. In this basis, the
charged- and neutral-current (CC and NC, respectively) Lagrangians read

$\displaystyle\mathcal{L}_{\rm
CC}=\frac{g}{\sqrt{2}}(W^{-})_{\mu}\bar{e}_{L}\gamma^{\mu}{\rm
K}\upnu_{L}+{\rm H.c.}\quad{\rm and}\quad\mathcal{L}_{\rm
NC}=\frac{g}{2\cos\theta_{W}}Z_{\mu}\bar{\upnu}_{L}({\rm K}^{\dagger}{\rm
K})\upnu_{L}+{\rm H.c.}\,,$ (2.1)

where we have introduced the matrix K, which relates the three active flavour
eigenstates and the $n$ mass eigenstates.111Note that the $3\times n$ matrix K
also includes the rotation required to define the charged-leptons mass basis.
At low energies, it is possible to express the CC and NC interactions in terms
of four-fermion interactions, i.e.

$\displaystyle\mathcal{L}_{\rm CC}=-2\sqrt{2}G_{\rm F}\sum_{i,j}({\rm
K}^{\dagger})_{ie}{\rm
K}_{ej}\left(\bar{\nu}_{i}\gamma^{\mu}P_{L}\nu_{j}\right)\left(\bar{e}\gamma_{\mu}P_{L}e\right)~{},$
(2.2)

and

$\displaystyle\mathcal{L}_{\rm NC}=-2\sqrt{2}G_{\rm
F}\sum_{X=L,R}g_{X}\sum_{i,j}\left({\rm K}^{\dagger}{\rm
K}\right)_{ij}\left(\bar{\nu}_{i}\gamma^{\mu}P_{L}\nu_{j}\right)\left(\bar{e}\gamma_{\mu}P_{X}e\right)\,,$
(2.3)

where in the case of Eq. (2.2) we have performed a Fierz transformation. Here
$G_{\rm F}$ is the Fermi constant and the index $X=\\{L,R\\}$, so that $P_{X}$
denotes the chiral projectors $P_{R,L}=(1\pm\gamma_{5})/2$.

Note that, when computing the amplitude of the decay and scattering processes,
the sum over mass eigenstates should be limited to the largest mass eigenstate
kinematically accessible. Thus, it is useful to separate the matrix K in two
blocks,

$\displaystyle{\rm K}=\begin{pmatrix}N&S\end{pmatrix}\,,$ (2.4)

where $N$ describes the mixing among the three lightest states and $S$
accounts for the mixing between the three lightest states and the $n-3$ heavy
ones.

In what follows, we shall assume that only three mass eigenstates lie below
the energy scale of interest here. In that case, the sums in Eqs. (2.2)-(2.3)
only include the three lightest mass eigenstates, and only the $3\times 3$
submatrix $N$ is involved. Notice that the unitarity of the full $n\times n$
lepton mixing matrix does not require the submatrix $N$ appearing in the
effective CC and NC Lagrangians to be unitary. In order to account for
departures of non-unitarity in $N$, we adopt the following parametrisation
[32]

$N=\begin{pmatrix}\alpha_{11}&0&0\\\ \alpha_{21}&\alpha_{22}&0\\\
\alpha_{31}&\alpha_{32}&\alpha_{33}\end{pmatrix}U\,,$ (2.5)

where $U$ is the usual unitary three-neutrino mixing matrix, $\alpha_{ii}$ are
real parameters, and $\alpha_{ij}~{}(i\neq j)$ are complex and provide
additional sources of CP violation 222The relation between the $\alpha_{ij}$
parameters and the parametrisation of the full $n\times n$ mixing matrix is
detailed in Appendix A.. The off-diagonal parameters are related to the
diagonal ones through the triangular inequality [33],

$|\alpha_{ij}|\leq\sqrt{(1-\alpha^{2}_{ii})(1-\alpha^{2}_{jj})}\,.$ (2.6)

Additionally, from the unitarity of the full $n\times n$ mixing matrix [34],
the following conditions apply

$\displaystyle\alpha^{2}_{11}|\alpha_{21}|^{2}$
$\displaystyle\leq(1-\alpha^{2}_{11})(1-\alpha^{2}_{22}-|\alpha_{21}|^{2})\,,$
(2.7) $\displaystyle\alpha^{2}_{11}|\alpha_{31}|^{2}$
$\displaystyle\leq(1-\alpha^{2}_{11})(1-\alpha^{2}_{33}-|\alpha_{31}|^{2}-|\alpha_{32}|^{2})\,,$
(2.8) $\displaystyle|\alpha_{22}\alpha_{32}+\alpha^{*}_{21}\alpha_{31}|^{2}$
$\displaystyle\leq(1-\alpha^{2}_{22}-|\alpha_{21}|^{2})(1-\alpha^{2}_{33}-|\alpha_{31}|^{2}-|\alpha_{32}|^{2})\,.$
(2.9)

Non-unitary three-neutrino mixing would not only alter CC and NC interactions,
but also the picture of neutrino oscillations, whose observations at both
short- and long-baseline experiments can be used to constrain the non-
unitarity parameters. Current bounds are summarized in Table 1, from [34].

NU parameters | $|\alpha_{21}|$ | $|\alpha_{31}|$ | $|\alpha_{32}|$ | $\alpha_{11}$ | $\alpha_{22}$ | $\alpha_{33}$
---|---|---|---|---|---|---
3$\sigma$ bounds | $<0.025$ | $<0.075$ | $<0.02$ | $>0.93$ | $>0.98$ | $>0.72$
Table 1: Current constraints on the non-unitary parameters from Ref. [34].

##  3 The impact of non-unitary three-neutrino mixing in neutrino decoupling

In the following, we shall concentrate in a scenario in which HNLs can not be
produced at temperatures $\sim$MeV and, therefore, their distribution is
Boltzmann-suppressed. Consequently, only the evolution of the three lightest
mass eigenstates needs to be studied. In previous studies of neutrino
decoupling, the evolution of the neutrino distributions was mostly addressed
in terms of flavour states [27, 28], but an adequate treatment of non-
unitarity requires to work in the mass basis (see e.g. [26]). The reason is
that the truncation in the sum over accessible states is only well defined for
mass eigenstates.

We adopt the usual definitions for the comoving variables $x=m_{e}\,a$,
$y=p\,a$, $z=T_{\gamma}\,a$, in terms of the electron mass $m_{e}$, the
cosmological scale factor $a$, the neutrino momentum $p$, and the photon
temperature $T_{\gamma}$ [23, 28]. The evolution of the density matrix
$\varrho$ for tree neutrinos in the mass basis reads [26]

$\frac{{\rm d}\varrho(y)}{{\rm d}x}=\sqrt{\frac{3m^{2}_{\rm
Pl}}{8\pi\rho}}\left\\{-i\frac{x^{2}}{m_{e}^{3}}\left[\frac{\mathbb{M}}{2y}-\frac{2\sqrt{2}G_{\rm
F}ym_{e}^{6}}{x^{6}}\left(\frac{\mathbb{E}_{\ell}+\mathbb{P}_{\ell}}{m_{W}^{2}}\right),\varrho\right]+\frac{m_{e}^{3}}{x^{4}}\mathcal{I(\varrho)}\right\\}\,,$
(3.1)

where $m_{\rm Pl}$ denotes the Planck mass, $\rho$ is the energy density of
the Universe, $m_{W}$ is the mass of the $W$ gauge bosons and we have defined
the diagonal neutrino mass matrix $\mathbb{M}={\rm diag}(0,\Delta
m^{2}_{21},\Delta m^{2}_{31})$. We have also introduced the energy density of
charged leptons, $\mathbb{E}_{\ell}$, and the pressure of charged leptons,
$\mathbb{P}_{\ell}$. In this scenario, it is useful to define the following
matrices

$\displaystyle\left(Y_{L}\right)_{ij}$
$\displaystyle=g_{L}\left(N^{\dagger}N\right)_{ij}+(N^{\dagger})_{ie}N_{ej}\,,$
(3.2) $\displaystyle\left(Y_{R}\right)_{ij}$
$\displaystyle=g_{R}\left(N^{\dagger}N\right)_{ij}\,,$ (3.3)

with the $i$ and $j$ indices running from 1 to 3, such that the charged-lepton
energy density and pressure,

$\displaystyle\mathbb{E}_{\ell}=\rho_{e}\,\mathbb{E}_{\rm NU}~{},$
$\displaystyle\mathbb{P}_{\ell}=P_{e}\,\mathbb{E}_{\rm
NU}\simeq\rho_{e}\,\mathbb{E}_{\rm NU}/3~{},$ (3.4)

are expressed in terms of the electron energy density and the pressure of
charged leptons ($\rho_{e}$ and $P_{e}$, respectively) and also via a matrix
accounting for the contribution of charged leptons to the matter potential

$\mathbb{E}_{\rm NU}=Y_{L}-Y_{R}\,.$ (3.5)

Finally, the term $\mathcal{I(\varrho)}$ in (3.1) encodes the collision
integrals, arising from the neutrino-electron scattering and annihilation
processes, as well as from the neutrino-neutrino interactions. The terms
accounting for neutrino-electron interactions are proportional to $G_{\rm
F}^{2}$ and also depend on the matrices $Y_{L}$ and $Y_{R}$ defined in Eqs.
(3.2)-(3.3), see Appendix B for details. Neutrino-neutrino interactions are
also modified if the three-neutrino mixing is not unitary. The low-energy
Lagrangian accounting for neutrino self-interactions reads

$\mathcal{L}_{\nu SI}=-2\sqrt{2}G_{\rm
F}\sum_{i,j,k,m}\left(N^{\dagger}N\right)_{ij}\left(N^{\dagger}N\right)_{mk}(\bar{\nu}_{i}\gamma^{\mu}P_{L}\nu_{j})(\bar{\nu}_{k}\gamma_{\mu}P_{L}\nu_{m})\,.$
(3.6)

It has been shown that neutrino self-interactions lead to a refractive term in
the Hamiltonian, which does not modify the picture qualitatively and it only
gives a sub-leading contribution to $N_{\rm eff}$ 333The role of neutrino
self-interactions in cosmology has been studied extensively and was proposed
as a solution to the so-called Hubble tension, see e.g. Ref. [35] and
references therein. Note, however, that self-interactions invoked in those
scenarios are orders of magnitude larger than the one expected to arise from
the non-unitarity of the three-neutrino mixing., hence we ignore such terms.

Notice that we have treated the Fermi constant $G_{\rm F}$ as a fundamental
quantity that we can precisely measure experimentally. However, in the
presence of non-unitary three-neutrino mixing, what it is measured is an
effective parameter depending on the mixing between light neutrinos and HNLs
[32]. The effective parameter measured in the case of $\beta$-decay is

$G^{\beta}_{\rm F}=G_{\rm F}\sqrt{(NN^{\dagger})_{ee}}=G_{\rm
F}\alpha_{11}~{},$ (3.7)

whereas for muon decay it is

$G^{\mu}_{\rm F}=G_{\rm
F}\sqrt{(NN^{\dagger})_{ee}(NN^{\dagger})_{\mu\mu}}=G_{\rm{F}}\sqrt{\alpha^{2}_{11}(\alpha^{2}_{22}+|\alpha_{21}|^{2})}\,.$
(3.8)

The measurement from muon decay has achieved a higher precision and,
consequently, we take it as a reference, and consider $G^{\mu}_{\rm
F}=1.1663787(6)\times 10^{-5}\,{\rm GeV}^{-2}$ [36].

##  4 Results

We make use in what follows of the publicly available code FortEPiaNO (Fortran
Evolved Primordial Neutrino Oscillations) [23, 28], which has been modified in
order to work in the mass basis and to include the effects from a possible
non-unitary three-neutrino mixing. The numerical settings adopted here allow
to compute the value of $N_{\rm eff}$ with a precision better than $10^{-3}$,
which is safely smaller than the expected sensitivity from future surveys and,
therefore, accurate enough for the study of non-unitary three-neutrino mixing.

###  4.1 The role of diagonal parameters $\mathbf{\alpha_{ii}}$

According to the parametrisation in Eq. (2.5), non-unitarity could manifest as
a deviation from unity in the parameters $\alpha_{ii}$ (see Appendix A for the
relation with the parametrisation of the $n\times n$ unitary mixing matrix).
To understand their impact, we address the role of each one individually, as
well as the interplay between them.

Figure 1 depicts the dependence of $N_{\rm eff}$ on each of the individual
diagonal parameters, $\alpha_{ii}$. In the left panel, we show the exact
calculation of $N_{\rm eff}$ as a function of $\alpha_{ii}$, allowing only one
of them being different from zero simultaneously. It is known that, in other
non-standard scenarios, the main contribution of new physics to $N_{\rm eff}$
comes from the modification of the collision integrals. For illustrative
purposes, the top right panel shows the values of $N_{\rm eff}$ after
considering the dependence of $G_{\rm F}$ on each of the $\alpha_{ii}$
parameters only in the collision terms (and not in the Hamiltonian) whereas
the bottom right panel displays the opposite situation (note the different
scale for $N_{\rm eff}$).

Figure 1: Values of $N_{\rm eff}$ as a function of the diagonal non-unitarity
parameters $\alpha_{ii}$, when only one of them is allowed to be non-zero. For
comparison, we show the exact calculation (left panel), where the vertical
dotted lines indicate the existing limits (see Table 1), and the cases when
the dependence of $G_{\rm F}$ on the $\alpha_{ii}$ parameters is considered
only in the collision integrals (upper right panel) or in the Hamiltonian
(lower right panel).

Notice that, since the relation between the true value of $G_{\rm F}$ and the
one measured from muon decay does not depend on $\alpha_{33}$, the calculated
$N_{\rm eff}$ is almost insensitive to this parameter. Indeed, the main impact
of non-unitarity on neutrino decoupling arises from the mismatch between the
true and the observed value of $G_{\rm F}$. The second key point is that, as
it can be noticed from the right plots of Figure 1, the collision terms
dominate the dependence of $N_{\rm eff}$ on the non-unitarity parameters:
changes on the effective $G_{\rm F}$ alter both the scattering and
annihilation processes involving neutrinos, electrons and positrons, leading
to a delayed decoupling and hence to a larger value of $N_{\rm eff}$. Albeit
matter effects also depend on the effective $G_{\rm F}$, this is not as
relevant for neutrino decoupling as the interactions with the cosmic plasma.
From the current measurement of $N_{\rm eff}=2.99^{+0.34}_{-0.33}$ at 95$\%$
C.L. [29], it is possible to set weak limits on the non-unitarity of the
three-neutrino mixing matrix,

$\displaystyle\alpha_{11}>0.07\,,$ $\displaystyle\alpha_{22}>0.15\,,$ (4.1)

both at $95\%$ C.L. As an exercise, we predict the improvement on these limits
from future CMB observations. Let us assume a forecasted sensitivity to
$N_{\rm eff}$ with an uncertainty $\sigma=0.02$ [31], together with a Gaussian
posterior distribution for this parameter. Then, one can define the $\chi^{2}$
test:

$\displaystyle\chi^{2}=\frac{[(N_{\rm eff})_{0}-N_{\rm
eff}]^{2}}{\sigma^{2}}\,,$ (4.2)

where $(N_{\rm eff})_{0}$ corresponds to the measured value of $N_{\rm eff}$,
assumed to match the standard theoretical prediction, i.e. $(N_{\rm
eff})_{0}=3.044$. Under these assumptions, we find that future CMB
measurements will improve considerably the former $95\%$ C.L. upper limits:

$\displaystyle\alpha_{11}>0.29\,,$ $\displaystyle\alpha_{22}>0.50\,.$ (4.3)

In the following, we shall present an exhaustive study of the parameter
correlations. Given the $\chi^{2}$ test definition in Eq. (4.2), one can
extract the confidence regions expected when two parameters are allowed to
differ from unity simultaneously in the analyses.

Figure 2: The three panels show the 1$\sigma$, 2$\sigma$ and 3$\sigma$ allowed
regions expected from the measurement of $N_{\rm eff}$ when two of the
diagonal parameters are allowed to vary simultaneously. The top panel shows
the $\alpha_{11}-\alpha_{22}$ plane, whereas the bottom left and bottom right
panels correspond to the $\alpha_{11}-\alpha_{33}$ and
$\alpha_{22}-\alpha_{33}$ planes, respectively. The isocontour of $N_{\rm
eff}$=$(N_{\rm eff})_{0}$=3.044 is shown in white, together with a shaded band
indicating the numerical uncertainty, $(N_{\rm eff})_{0}\pm 0.001$.

Figure 2 depicts the expected sensitivity to the non-unitarity diagonal
parameters, projected in two-dimensional contours. The three coloured regions
correspond to the 1$\sigma$, 2$\sigma$ and 3$\sigma$ deviations from the
standard value for $N_{\rm eff}$, shown by the white line. The white shaded
band corresponds to an estimate of the numerical error in the standard
prediction, given the numerical settings chosen in this work.

In the top panel, we show how the expected limit on $N_{\rm eff}$ translates
into allowed regions in the $\alpha_{11}-\alpha_{22}$ two-dimensional plane.
Although these regions are larger than the ones derived from the comparison of
the measurements of $G_{\rm F}$ from $\beta$ and muon decays or from neutrino
oscillation experiments, they provide independent constraints on the non-
unitarity diagonal parameters. The lower panels correspond to the sensitivity
in the $\alpha_{11}-\alpha_{33}$ and $\alpha_{22}-\alpha_{33}$ two-dimensional
planes. As it was already shown in Figure 1, values of $\alpha_{33}$ which
deviate significantly from unity do not result in a relevant change in $N_{\rm
eff}$: in particular, for values of $\alpha_{33}<0.5$, $N_{\rm eff}$ becomes
almost independent of $\alpha_{33}$. The reason for the small dependence of
$N_{\rm eff}$ on $\alpha_{33}$ comes from the fact that the former only
affects interaction rates involving the $\tau$ neutrino, which have a minor
role in the neutrino decoupling process. On the contrary, $\alpha_{11}$ and
$\alpha_{22}$ tend to result in a larger value of $N_{\rm eff}$ because of
their effect on $G_{\rm F}$.

###  4.2 Including off-diagonal parameters

Allowing for non-zero off-diagonal parameters, $\alpha_{ij}>0$ with $i<j$,
changes the prediction for the effective number of relativistic species, as
these parameters modify the structure of neutrino interactions with the plasma
as well. In addition, for non-zero $\alpha_{21}$, the relation between the
effective Fermi constant measure from muon decay, $G_{\rm F}^{\mu}$, and the
corresponding fundamental constant $G_{\rm F}$ is altered. In what follows, we
present the results of adding one extra non-zero off-diagonal parameter to the
three diagonal ones, limiting the analysis to real values of
$\alpha_{ij}$444This choice is motivated by the results in [37], where it was
shown that the CP phase does not alter the predicted $N_{\rm eff}$ in the case
of unitary mixing..

We show in Figure 3 how non-zero values of $\alpha_{21}$ satisfying Eqs.
(2.7)-(2.9) enlarge the parameter space allowed in the
$\alpha_{11}-\alpha_{22}$ two-dimensional plane after marginalising over the
$\alpha_{21}$ parameter, see the upper left panel of Figure 2 for comparison.
Notice that non-zero values of the $\alpha_{21}$ parameter can partially
compensate the larger value of $N_{\rm eff}$ expected when $\alpha_{11}$ and
$\alpha_{22}$ depart from one, thanks to its effect on $G_{\rm F}$.
Conversely, non-zero values of the off-diagonal parameters $\alpha_{31}$ and
$\alpha_{32}$ were found to introduce only sub-leading changes in the planes
$\alpha_{11}-\alpha_{33}$ and $\alpha_{22}-\alpha_{33}$, respectively. This
implies that $\alpha_{31}$ and $\alpha_{32}$ do not reduce the value of
$N_{\rm eff}$ but rather increase it, leaving therefore the two-dimensional
contours in the $\alpha_{11}-\alpha_{33}$ and $\alpha_{22}-\alpha_{33}$ planes
almost unchanged.

Figure 3: Expected allowed regions at 1$\sigma$, 2$\sigma$ and 3$\sigma$
significance from the measurement of $N_{\rm eff}$ in the
$\alpha_{11}-\alpha_{22}$ plane, when $\alpha_{21}$ is allowed to vary while
respecting the unitarity condition Eq. (2.7) and after marginalizing over
$\alpha_{21}$, see main text for details.

##  5 Discussion

###  5.1 From non-unitarity to non-standard interactions

In the case in which the deviations from unitarity of the three-neutrino
mixing matrix are small, one can write the effective Lagrangian from charged-
and neutral- current interactions in terms of the active neutrino states, as

$\displaystyle\mathcal{L}_{NU}=$ $\displaystyle-2\sqrt{2}G_{\rm
F}\left[\sum_{\alpha,\beta}\left((NN^{\dagger})_{\alpha
e}(NN^{\dagger})_{e\beta}+g_{L}\left[NN^{\dagger}\right]^{2}_{\alpha\beta}\right)\left(\bar{\nu}_{\alpha}\gamma^{\mu}P_{L}\nu_{\beta}\right)\left(\bar{e}\gamma_{\mu}P_{L}e\right)\right.$
$\displaystyle\left.+g_{R}\sum_{\alpha,\beta}\left[NN^{\dagger}\right]^{2}_{\alpha\beta}\left(\bar{\nu}_{\alpha}\gamma^{\mu}P_{L}\nu_{\beta}\right)\left(\bar{e}\gamma_{\mu}P_{R}e\right)\right]\,.$
(5.1)

This means that, if we take into account that the value of $G_{\rm F}$
measured from muon decay is given by Eq. (3.8), then one can map the scenario
of small deviations from unitarity to the one of neutrino-electron neutral
current non-standard interactions (NSI). If we define the corresponding
neutral current Lagrangian with NSI as

$\displaystyle\mathcal{L}_{SM+NSIe}=$ $\displaystyle-2\sqrt{2}\,G_{\rm
F}\left(\overline{\nu}_{e}\gamma^{\mu}P_{L}e\right)(\overline{e}\gamma_{\mu}P_{L}\nu_{e})$
$\displaystyle-2\sqrt{2}\,G_{\rm
F}\sum_{X,\alpha}\left[g_{X}\left(\overline{\nu}_{\alpha}\gamma^{\mu}P_{L}\nu_{\alpha}\right)(\overline{e}\gamma_{\mu}P_{X}e)+\sum_{\beta}\varepsilon_{\alpha\beta}^{eX}\left(\overline{\nu}_{\alpha}\gamma^{\mu}P_{L}\nu_{\beta}\right)(\overline{e}\gamma_{\mu}P_{X}e)\right],$
(5.2)

then, the NSI-like parameters read as

$\displaystyle\varepsilon^{eL}_{\alpha\beta}=-(\delta_{\beta e}\delta_{\alpha
e}+g_{L}\delta_{\alpha\beta})+\frac{(NN^{\dagger})_{\alpha
e}(NN^{\dagger})_{e\beta}+g_{L}\left[NN^{\dagger}\right]^{2}_{\alpha\beta}}{\sqrt{\alpha^{2}_{11}(\alpha^{2}_{22}+|\alpha_{21}|^{2})}}~{},$
(5.3)

and

$\displaystyle\varepsilon^{eR}_{\alpha\beta}=-g_{R}\delta_{\alpha\beta}+g_{R}\frac{\left[NN^{\dagger}\right]^{2}_{\alpha\beta}}{\sqrt{\alpha^{2}_{11}(\alpha^{2}_{22}+|\alpha_{21}|^{2})}}\,.$
(5.4)

The involved matrix elements, ($NN^{\dagger})_{\alpha\beta}$, can be computed
at leading order [33] and then, a deeper physical understanding of the results
becomes possible. The equivalent NSI parameters resulting from mapping
departures from unity in the diagonal parameters are shown in Table 2. If the
interactions between neutrinos and the cosmic plasma become weaker, one
expects an earlier decoupling, whereas stronger interactions would keep
neutrinos in thermal contact for a longer time, and, consequently, neutrinos
could receive more energy from the electron-positron pair annihilation.

NU Parameter | Effective NSI
---|---
$\alpha_{11}<$ 1 |  $\varepsilon^{eL}_{\alpha\beta}$ = | $(1+g_{L})(-1+\alpha_{11}^{3})\delta_{\beta e}\delta_{\alpha e}+g_{L}\left(\frac{1}{\alpha_{11}}-1\right)(\delta_{\alpha\mu}\delta_{\beta\mu}+\delta_{\alpha\tau}\delta_{\beta\tau})$
| $\varepsilon^{eR}_{\alpha\beta}=$ | $-g_{R}(1-\alpha^{3}_{11})\delta_{\alpha e}\delta_{\beta e}+g_{R}\left(\frac{1}{\alpha_{11}}-1\right)(\delta_{\alpha\mu}\delta_{\beta\mu}+\delta_{\alpha\tau}\delta_{\beta\tau})$
$\alpha_{22}<$ 1 |  $\varepsilon^{eL}_{\alpha\beta}$ = | $-(1+g_{L})\left(1-\frac{1}{\alpha_{22}}\right)\delta_{\alpha e}\delta_{\beta e}$
|  | $-g_{L}(1-\frac{1}{\alpha_{22}})\delta_{\alpha\tau}\delta_{\beta\tau}-g_{L}(1-{\alpha^{3}_{22}})\delta_{\alpha\mu}\delta_{\beta\mu}$
| $\varepsilon^{eR}_{\alpha\beta}=$ | $-g_{R}\left(1-\frac{1}{\alpha_{22}}\right)(\delta_{\alpha e}\delta_{\beta e}+\delta_{\alpha\tau}\delta_{\beta\tau})-g_{R}(1-\alpha^{3}_{22})\delta_{\alpha\mu}\delta_{\beta\mu}$
$\alpha_{33}<$ 1 |  $\varepsilon^{eL}_{\alpha\beta}$ = | $g_{L}(-1+\alpha_{33}^{4})\delta_{\alpha\tau}\delta_{\tau\beta}$
| $\varepsilon^{eR}_{\alpha\beta}=$ | $g_{R}(-1+\alpha_{33}^{4})\delta_{\alpha\tau}\delta_{\tau\beta}$
Table 2: Mapping between the diagonal non-unitarity (NU) parameters and
effective NSI with electrons.

From Table 2, it is possible to notice that the behaviour of the predictions
in Figure 1 result from the modification of the interactions, which, depending
on the values of $\alpha_{11}$ and $\alpha_{22}$, can shift to earlier or
later times the neutrino decoupling process. Conversely, the results obtained
for $\alpha_{33}$ are clearly understood from the fact that a value
$\alpha_{33}<1$ reduces the strength of the neutral-current interactions,
keeping tau neutrinos in thermal contact with the plasma, and leading to an
earlier decoupling.

When two or more parameters are needed to parametrise the non-unitarity of the
three-neutrino mixing, the picture becomes more complicated since the
different effects can compensate. Nevertheless, one can get an idea of the
results obtained when three parameters were allowed to vary simultaneously.
For instance, $\alpha_{31}>0$ and $\alpha_{32}>0$ result in the appearance of
flavour-changing neutral currents, which lead to a delayed decoupling [38].
The scenario of $\alpha_{21}$ is different since it modifies both the
structure of interactions and their strength.

Note that the treatment above is only valid for small departures from
unitarity. Otherwise, the calculations must be performed in the mass basis in
order to consider only the kinematically accessible states.

###  5.2 Comparison to other limits

At Earth, several precision and flavour observables can also be used to
constrain deviations from unitarity in the three-neutrino mixing matrix.
Whereas a detailed global analysis of all the existing constraints is beyond
the scope of this work, for completeness, we discuss some of the most
restrictive bounds from references [25, 39, 33, 40, 41].

At tree level, the weak mixing angle is related to the Fermi constant by

$\displaystyle\sin^{2}\theta_{W}\cos^{2}\theta_{W}=\frac{\alpha\pi}{\sqrt{2}G_{\rm
F}\,m^{2}_{Z}}=\frac{\alpha\pi}{\sqrt{2}G^{\mu}_{F}\,m^{2}_{Z}}\sqrt{(NN^{\dagger})_{ee}(NN^{\dagger})_{\mu\mu}}\,,$
(5.5)

where $\alpha$ is the fine structure constant. Requiring the agreement between
measurements of the weak mixing angle and its indirect inference using $G_{\rm
F}$ allows to set limits on departures from unitarity for the matrix $N$.

The prediction of weak decay widths of the Z boson also depends on the
assumptions about the unitarity of the lepton mixing matrix. At tree level,
the decay width to two fermions depends linearly on $G_{\rm F}$, whereas for
the invisible decay width,

$\displaystyle\Gamma_{\rm inv}\propto G_{\rm
F}\sum_{i,j=1,2,3}|(N^{\dagger}N)_{ij}|^{2}\,.$ (5.6)

Also at tree level, the mass of the W boson can be expressed in terms of the
Fermi constant, $G_{\rm F}$, as

$\displaystyle m^{2}_{W}=\frac{\alpha\pi}{\sqrt{2}G_{\rm
F}\sin^{2}\theta_{W}}\,.$ (5.7)

Again, bounds on departures from unitarity arise from requiring the agreement
between direct measurements and the value of $m_{W}$ derived using
$G^{\mu}_{\rm F}$.

As mentioned above, the values of $G_{\rm F}$ derived from $\beta$ and muon
decays would receive different corrections. As a consequence, limits on the
unitarity of the Cabibbo-Kobayashi-Maskawa matrix can be recasted, since

$\displaystyle\sum_{i=1}^{3}|V_{ui}|^{2}=\left(\frac{G^{\beta}_{\rm
F}}{G^{\mu}_{\rm
F}}\right)^{2}=\frac{1}{\alpha_{22}^{2}+|\alpha_{21}|^{2}}\,.$ (5.8)

In the presence of non-unitarity, couplings between leptons and gauge bosons
are no longer flavour independent, implying that semileptonic decay rates can
be flavour-dependent. Thus, measurements of pion and kaon decays also restrict
deviations from unitarity. One can also consider universality tests in the
leptonic sector and W-boson decays. For instance,

$\displaystyle R^{W}_{\alpha\beta}=\frac{\Gamma(W\rightarrow
l_{\alpha}\bar{\nu}_{\alpha})}{\Gamma(W\rightarrow
l_{\beta}\bar{\nu}_{\beta})}=\sqrt{\frac{(NN^{\dagger})_{\alpha\alpha}}{(NN^{\dagger})_{\beta\beta}}}\,.$
(5.9)

These and other probes have been studied extensively in the literature and
provide stringent limits to possible deviations from unitarity of the three-
neutrino mixing matrix [39, 41, 33, 40]. Table 3 reports the constraints on
the $\alpha_{ij}$ parameters obtained in [33] when considering constraints
from neutrinos and charged leptons. Constraints from neutrinos alone [34] are
complementary albeit not completely competitive.

$\alpha_{11}$ | $\alpha_{22}$ | $\alpha_{33}$ | $|\alpha_{21}|$ | $|\alpha_{31}|$ | $|\alpha_{31}|$
---|---|---|---|---|---
$>$ 0.9961 |  $>$ 0.9990 |  $>$ 0.9973 |  $<$ 2.6$\cdot 10^{-3}$ |  $<$ 5.0$\cdot 10^{-3}$ |  $<$ 2.4$\cdot 10^{-3}$
Table 3: Combined limits on the non-unitarity parameters $\alpha$ from charged
lepton and neutrino measurements at 90% C.L. (considering 6 degrees of
freedom), from [33].

##  6 Concluding remarks

In this work, we have considered the impact of non-unitarity in the three-
neutrino mixing matrix on the process of neutrino decoupling. The presence of
non-unitarity, a smoking-gun low-energy effect for the presence of heavy
neutral leptons, alters the structure of both charged and neutral current
interactions, thus affecting how neutrinos decouple from the thermal plasma in
the early Universe.

In light of the forecasted accuracy expected from forthcoming cosmological
surveys, we have derived the deviations on $N_{\rm eff}$ induced by a non-
unitarity mixing in the three-neutrino sector. In this regard, future
cosmological measurements can serve as a novel high-redshift tool to provide
model-independent and laboratory-complementary constraints on the non-
unitarity parameters arising from new physics. Furthermore, they can be used
as a cross-check to either falsify (or confirm) neutrino mass-generation
scenarios.

We employ the commonly exploited parameterisation for the non-unitarity
parameters in Eq. (2.5). If the diagonal parameters $\alpha_{11}$ and
$\alpha_{22}$ are different from unity, they will delay neutrino decoupling,
resulting in a larger value of $N_{\rm eff}$. The lower limits we obtained
here when considering the future CMB-S4 mission combined with large-scale
structure observations are $\alpha_{11}>0.29$ and $\alpha_{22}>0.50$, both at
$95\%$ C.L. Concerning the third non-unitarity diagonal parameter
$\alpha_{33}$, it would barely alter the overall picture. We have also
explored the possible interplay between the diagonal and off-diagonal non-
unitarity parameters when either two or three of them are allowed to vary
simultaneously. Interestingly, a non-zero value $\alpha_{21}$ relaxes
considerably the constraints in the $\alpha_{11}-\alpha_{22}$ plane.

We conclude here emphasizing the fact that, even if the constraints on the
non-unitarity parameters expected from cosmological observations are among the
weakest when compared to other terrestrial probes, they provide an independent
limit and a consistency test of the neutrino decoupling picture. Furthermore,
they allow to set constraints on the putative non-unitary three-neutrino
mixing at a completely different epoch in the universe evolution. Therefore,
they should not be regarded as competitive limits but as supportive searches.

## Acknowledgments

SG and OM acknowledge financial support from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 754496 (project FELLINI) and the ITN grant H2020-MSCA-
ITN-2019/860881-HIDDeN, respectively. In Valencia, this work was supported by
the Spanish grants PID2020-113775GB-I00 and PID2020-113644GB-I00
(MCIN/AEI/10.13039/501100011033), as well as projects PROMETEO/2019/083 and
CIPROM/2021/054 (Generalitat Valenciana). PMM is also supported by the grant
FPU18/04571 (MICIU). He is grateful for the hospitality of the Particle and
Astroparticle Physics Division of the Max-Planck-Institut für Kernphysik
(Heidelberg), the Centro de Física Teórica de Partículas (Instituto Superior
Técnico, Lisbon) and the AstroNu group at Niels Bohr Institute (Copenhagen)
during the development of this project.

##  Appendix A Parametrisation of the non-unitary three-neutrino mixing matrix

The non-unitarity of the three-neutrino mixing matrix can be parametrised by a
triangular matrix. The parameters introduced can be related to the full
$n\times n$ unitary lepton mixing matrix, $U_{n\times n}$. If we define the
following complex rotation matrix,

$\displaystyle\omega_{ij}=\begin{pmatrix}1&0&&\cdots&0&\cdots&&&0\\\
0&1&&&&&&&\vdots\\\ \vdots&&c_{ij}&\cdots&0&\cdots&\eta_{ij}&&\\\
&&\vdots&\ddots&&&\vdots&\\\ &&0&&1&&0&&\\\ &&\vdots&&\ddots&&\vdots&&\\\
&&\bar{\eta}_{ij}&\cdots&0&\cdots&c_{ij}&&\vdots\\\ \vdots&&&&&&&1&0\\\
0&&&\cdots&0&\cdots&&0&1\end{pmatrix}\,,$ (A.1)

where we use the notation $c_{ij}=\cos\theta_{ij}$,
$\eta_{ij}=\sin\theta_{ij}\,e^{-\phi_{ij}}$ and
$\bar{\eta}_{ij}=-\sin\theta_{ij}\,e^{i\phi_{ij}}$, the full lepton mixing
matrix can be parametrised as

$\displaystyle U_{n\times
n}=\omega_{n-1\,n}\,\omega_{n-2\,n}...\,\omega_{1\,n}\omega_{n-2\,n-1}\,\omega_{n-3\,n-1}...\,\omega_{1\,n-1}...\,\omega_{23}\,\omega_{13}\,\omega_{12}\,.$
(A.2)

The diagonal parameters $\alpha_{ii}$ in Eq. (2.5) are real and expressed as
[32, 34]

$\displaystyle\alpha_{ii}=c_{i\,n}\,c_{i\,n-1}...\,c_{i\,4}\,,$ (A.3)

and the off-diagonal ones are

$\displaystyle\alpha_{21}=\,$ $\displaystyle
c_{2\,n}\,c_{2\,n-1}...\,c_{25}\,\eta_{24}\,\bar{\eta}_{14}+c_{2\,n}...\,c_{26}\,\eta_{25}\,\bar{\eta}_{15}\,c_{14}+...+\eta_{2\,n}\,\bar{\eta}_{1n}\,c_{1\,n-1}...\,c_{14}\,,$
(A.4) $\displaystyle\alpha_{31}=\,$ $\displaystyle
c_{3\,n\,}c_{3\,n-1}...\,c_{35}\,\eta_{34}\,c_{24}\,\bar{\eta}_{14}+c_{3\,n}...\,c_{36}\,\eta_{35}\,c_{25}\,\bar{\eta}_{15}\,c_{14}+...+$
$\displaystyle\hskip
207.7052pt+\eta_{3\,n}\,c_{2\,n}\,\bar{\eta}_{1\,n}\,c_{1\,n-1}...\,c_{14}\,,$
(A.5) $\displaystyle\alpha_{32}=\,$ $\displaystyle
c_{3\,n}\,c_{3\,n-1}...\,c_{35}\,\eta_{34}\,\bar{\eta}_{24}+c_{3\,n}...\,c_{36}\,\eta_{35}\,\bar{\eta}_{25}\,c_{24}+...+\eta_{3\,n}\,\bar{\eta}_{2\,n}\,c_{2\,n-1}...\,c_{24}\,.$
(A.6)

##  Appendix B Collision integrals in the mass basis

Following the notation and procedure for integral reduction in [28], the
neutrino-electron collision integral is separated in the terms accounting for
scattering and annihilations,

$\displaystyle\mathcal{I}_{\nu
e}[\varrho(x,y)]=\frac{G^{2}_{F}}{(2\pi)^{3}y^{2}}\\{I_{\nu e}^{\rm
scatt}[\varrho(x,y)]+I_{\nu e}^{ann}[\varrho(x,y)]\\}\,,$ (B.1)

using the comoving variables defined in Section 3. The term accounting for
scattering reads

$\displaystyle I_{\nu e}^{\rm scatt}[\varrho(x,y)]=\int{\rm d}y_{2}{\rm
d}y_{3}\sum_{a=L,R}\sum_{b=L,R}A_{ab}(x,y,y_{2},y_{3})F_{\rm
sc}^{ab}\left(\varrho^{(1)},f_{e}^{(2)},\varrho^{(3)},f_{e}^{(4)}\right)\,,$
(B.2)

whereas the one related to annihilations is

$\displaystyle I_{\nu e}^{\rm ann}[\varrho(x,y)]=\int{\rm d}y_{2}{\rm
d}y_{3}\sum_{a=L,R}\sum_{b=L,R}B_{ab}(x,y,y_{2},y_{3})F^{\rm
ann}_{ab}\left(\varrho^{(1)},\varrho^{(2)},f_{e}^{(3)},f_{e}^{(4)}\right)\,.$
(B.3)

The exact dependence on the comoving variables of the scattering kernels,
$A_{ab}$ and $B_{ab}$, can be found in [28]. Their expressions do not depend
on the basis in which the evolution of the density matrix is computed (mass
basis here, flavour basis in [28]) nor on the strength and structure of the
interactions.

The fact that we are working in the mass basis only manifests in the phase-
space factors $F^{\rm ann}_{ab}$ and $F^{\rm scatt}_{ab}$, which depend on the
matrices $Y_{L}$ and $Y_{R}$ defined in Eqs. (3.2)-(3.3). Although not stated
explicitly, the density matrix $\varrho$ is also in the mass basis. Then, the
expressions of the factors accounting for the phase-space of scattering and
annihilation processes read

$\displaystyle F^{\rm
scatt}_{ab}\left(\varrho^{(1)},f_{e}^{(2)},\varrho^{(3)},f_{e}^{(4)}\right)$
$\displaystyle=$ $\displaystyle
f_{e}^{(4)}(1-f_{e}^{(2)})\left[Y_{a}\varrho^{(3)}Y_{b}(1-\varrho^{(1)})+(1-\varrho^{(1)})Y_{b}\varrho^{(3)}Y_{a}\right]$
(B.4) $\displaystyle-$ $\displaystyle
f_{e}^{(2)}(1-f_{e}^{(4)})\left[\varrho^{(1)}Y_{b}(1-\varrho^{(3)})Y_{a}+Y_{a}(1-\varrho^{(3)})Y_{b}\varrho^{(1)}\right]\,,$
$\displaystyle F^{\rm
ann}_{ab}\left(\varrho^{(1)},\varrho^{(2)},f_{e}^{(3)},f_{e}^{(4)}\right)$
$\displaystyle=$ $\displaystyle
f_{e}^{(3)}f_{e}^{(4)}\left[Y_{a}(1-\varrho^{(2)})Y_{b}(1-\varrho^{(1)})+(1-\varrho^{(1)})Y_{b}(1-\varrho^{(2)})Y_{a}\right]$
(B.5) $\displaystyle-$
$\displaystyle(1-f_{e}^{(3)})(1-f_{e}^{(4)})\left[Y_{a}\varrho^{(2)}Y_{b}\varrho^{(1)}+\varrho^{(1)}Y_{b}\varrho^{(2)}Y_{a}\right],$

where $a,b=L,R$. These expressions depend both on the electron momentum
distribution function, $f_{e}$, and on the density matrix
$\varrho^{(i)}=\varrho(y_{i})$, where $y_{i}$ is the comoving momentum of
particle $i$.

## References

  * [1] P.F. de Salas et al., _2020 global reassessment of the neutrino oscillation picture_ , _JHEP_ 02 (2021) 071 [2006.11237].
  * [2] I. Esteban, M.C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, _The fate of hints: updated global analysis of three-flavor neutrino oscillations_ , _JHEP_ 09 (2020) 178 [2007.14792].
  * [3] F. Capozzi, E. Di Valentino, E. Lisi, A. Marrone, A. Melchiorri and A. Palazzo, _Unfinished fabric of the three neutrino paradigm_ , _Phys. Rev. D_ 104 (2021) 083031 [2107.00532].
  * [4] P. Minkowski, _$\mu\to e\gamma$ at a Rate of One Out of $10^{9}$ Muon Decays?_, _Phys. Lett. B_ 67 (1977) 421.
  * [5] T. Yanagida, _Horizontal gauge symmetry and masses of neutrinos_ , _Conf. Proc. C_ 7902131 (1979) 95.
  * [6] J. Schechter and J.W.F. Valle, _Neutrino Masses in SU(2) x U(1) Theories_ , _Phys. Rev. D_ 22 (1980) 2227.
  * [7] J. Schechter and J.W.F. Valle, _Neutrino Decay and Spontaneous Violation of Lepton Number_ , _Phys. Rev. D_ 25 (1982) 774.
  * [8] R.N. Mohapatra and G. Senjanovic, _Neutrino Masses and Mixings in Gauge Models with Spontaneous Parity Violation_ , _Phys. Rev. D_ 23 (1981) 165.
  * [9] W. Grimus and L. Lavoura, _The Seesaw mechanism at arbitrary order: Disentangling the small scale from the large scale_ , _JHEP_ 11 (2000) 042 [hep-ph/0008179].
  * [10] R.N. Mohapatra and J.W.F. Valle, _Neutrino Mass and Baryon Number Nonconservation in Superstring Models_ , _Phys. Rev. D_ 34 (1986) 1642.
  * [11] E.K. Akhmedov, M. Lindner, E. Schnapka and J.W.F. Valle, _Left-right symmetry breaking in NJL approach_ , _Phys. Lett. B_ 368 (1996) 270 [hep-ph/9507275].
  * [12] E.K. Akhmedov, M. Lindner, E. Schnapka and J.W.F. Valle, _Dynamical left-right symmetry breaking_ , _Phys. Rev. D_ 53 (1996) 2752 [hep-ph/9509255].
  * [13] M. Malinsky, J.C. Romao and J.W.F. Valle, _Novel supersymmetric SO(10) seesaw mechanism_ , _Phys. Rev. Lett._ 95 (2005) 161801 [hep-ph/0506296].
  * [14] M. Malinsky, T. Ohlsson, Z.-z. Xing and H. Zhang, _Non-unitary neutrino mixing and CP violation in the minimal inverse seesaw model_ , _Phys. Lett. B_ 679 (2009) 242 [0905.2889].
  * [15] M. Malinsky, T. Ohlsson and H. Zhang, _Non-unitarity effects in a realistic low-scale seesaw model_ , _Phys. Rev. D_ 79 (2009) 073009 [0903.1961].
  * [16] S. Gariazzo, C. Giunti, M. Laveder, Y.F. Li and E.M. Zavanin, _Light sterile neutrinos_ , _J. Phys. G_ 43 (2016) 033001 [1507.08204].
  * [17] C. Giunti and T. Lasserre, _eV-scale Sterile Neutrinos_ , _Ann. Rev. Nucl. Part. Sci._ 69 (2019) 163 [1901.08330].
  * [18] A. Palazzo, _Exploring Light Sterile Neutrinos at Long Baseline Experiments: A Review_ , _Universe_ 6 (2020) 41.
  * [19] S. Schoppmann, _Status of Anomalies and Sterile Neutrino Searches at Nuclear Reactors_ , _Universe_ 7 (2021) 360 [2109.13541].
  * [20] M. Archidiacono and S. Gariazzo, _Two Sides of the Same Coin: Sterile Neutrinos and Dark Radiation, Status and Perspectives_ , _Universe_ 8 (2022) 175 [2201.10319].
  * [21] S. Hannestad, I. Tamborra and T. Tram, _Thermalisation of light sterile neutrinos in the early universe_ , _JCAP_ 07 (2012) 025 [1204.5861].
  * [22] N. Saviano, A. Mirizzi, O. Pisanti, P.D. Serpico, G. Mangano and G. Miele, _Multi-momentum and multi-flavour active-sterile neutrino oscillations in the early universe: role of neutrino asymmetries and effects on nucleosynthesis_ , _Phys. Rev. D_ 87 (2013) 073006 [1302.1200].
  * [23] S. Gariazzo, P.F. de Salas and S. Pastor, _Thermalisation of sterile neutrinos in the early Universe in the 3+1 scheme with full mixing matrix_ , _JCAP_ 07 (2019) 014 [1905.11290].
  * [24] S. Hagstotz et al., _Bounds on light sterile neutrino mass and mixing from cosmology and laboratory searches_ , _Phys. Rev. D_ 104 (2021) 123524 [2003.02289].
  * [25] S. Antusch and O. Fischer, _Non-unitarity of the leptonic mixing matrix: Present bounds and future sensitivities_ , _JHEP_ 10 (2014) 094 [1407.6607].
  * [26] K. Akita and M. Yamaguchi, _A precision calculation of relic neutrino decoupling_ , _JCAP_ 08 (2020) 012 [2005.07047].
  * [27] J. Froustey, C. Pitrou and M.C. Volpe, _Neutrino decoupling including flavour oscillations and primordial nucleosynthesis_ , _JCAP_ 12 (2020) 015 [2008.01074].
  * [28] J.J. Bennett et al., _Towards a precision calculation of $N_{\rm eff}$ in the Standard Model II: Neutrino decoupling in the presence of flavour oscillations and finite-temperature QED_, _JCAP_ 04 (2021) 073 [2012.02726].
  * [29] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [30] Simons Observatory collaboration, _The Simons Observatory: Science goals and forecasts_ , _JCAP_ 02 (2019) 056 [1808.07445].
  * [31] CMB-S4 collaboration, _CMB-S4 Science Book, First Edition_ , 1610.02743.
  * [32] F.J. Escrihuela, D.V. Forero, O.G. Miranda, M. Tortola and J.W.F. Valle, _On the description of nonunitary neutrino mixing_ , _Phys. Rev. D_ 92 (2015) 053009 [1503.08879].
  * [33] F.J. Escrihuela, D.V. Forero, O.G. Miranda, M. Tórtola and J.W.F. Valle, _Probing CP violation with non-unitary mixing in long-baseline neutrino oscillation experiments: DUNE as a case study_ , _New J. Phys._ 19 (2017) 093005 [1612.07377].
  * [34] D.V. Forero, C. Giunti, C.A. Ternes and M. Tortola, _Nonunitary neutrino mixing in short and long-baseline experiments_ , _Phys. Rev. D_ 104 (2021) 075030 [2103.01998].
  * [35] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri et al., _In the realm of the Hubble tension—a review of solutions_ , _Class. Quant. Grav._ 38 (2021) 153001 [2103.01183].
  * [36] E. Tiesinga, P.J. Mohr, D.B. Newell and B.N. Taylor, _Codata recommended values of the fundamental physical constants: 2018_ , _Rev. Mod. Phys._ 93 (2021) 025010.
  * [37] J. Froustey and C. Pitrou, _Primordial neutrino asymmetry evolution with full mean-field effects and collisions_ , _JCAP_ 03 (2022) 065 [2110.11889].
  * [38] P.F. de Salas, S. Gariazzo, P. Martínez-Miravé, S. Pastor and M. Tórtola, _Cosmological radiation density with non-standard neutrino-electron interactions_ , _Phys. Lett. B_ 820 (2021) 136508 [2105.08168].
  * [39] S. Antusch and O. Fischer, _Probing the nonunitarity of the leptonic mixing matrix at the CEPC_ , _Int. J. Mod. Phys. A_ 31 (2016) 1644006 [1604.00208].
  * [40] F.J. Escrihuela, L.J. Flores and O.G. Miranda, _Neutrino counting experiments and non-unitarity from LEP and future experiments_ , _Phys. Lett. B_ 802 (2020) 135241 [1907.12675].
  * [41] E. Fernandez-Martinez, J. Hernandez-Garcia and J. Lopez-Pavon, _Global constraints on heavy neutrino mixing_ , _JHEP_ 08 (2016) 033 [1605.08774].