# Deep Bayesian ICP Covariance Estimation

Andrea De Maio1 and Simon Lacroix1 1LAAS-CNRS, Université de Toulouse, 7
Avenue du Colonel Roche, 31031, Toulouse, France, E-mail:
andrea[dot]de[dash]maio[at]laas[dot]fr

###### Abstract

Covariance estimation for the Iterative Closest Point (ICP) point cloud
registration algorithm is essential for state estimation and sensor fusion
purposes. We argue that a major source of error for ICP is in the input data
itself, from the sensor noise to the scene geometry. Benefiting from recent
developments in deep learning for point clouds, we propose a data-driven
approach to learn an error model for ICP. We estimate covariances modeling
data-dependent heteroscedastic aleatoric uncertainty, and epistemic
uncertainty using a variational Bayesian approach. The system evaluation is
performed on LiDAR odometry on different datasets, highlighting good results
in comparison to the state of the art.

##  I Introduction

With the steady increase in 3D sensors availability, point clouds have seen a
rapid diffusion and adoption. Their use benefits from the 3D information that
they readily express, which allows to easily retrieve geometric properties and
yield the possibility to extract complex geometric features to classify object
shapes, perform scene segmentation, or assess complex situations. Point clouds
are also used to solve motion estimation and self-localization problems in
robotics. Scan registration algorithms provide an estimate of the 3D
transformation between the positions at which two point clouds are acquired.
Such processes are at the heart of LiDAR-based odometry [1] and SLAM
frameworks [2]. Originally introduced in [3], the Iterative Closest Point
(ICP) algorithm, in its many variants, gradually became the standard approach
to the point-cloud registration problem [4].

In order to qualify the result of ICP, and especially to integrate it within a
localization framework, it is essential to estimate its uncertainty. ICP
errors stem from various error sources. Among those, the structure of the
scene plays an important role: underconstrained situations, such as corridors
or mostly flat environments, yield larger errors than geometrically more
complex scenes (Fig. 1).

Building on the advances in deep learning techniques dealing with point
clouds, we present an approach that estimates data-dependent error models for
the ICP process in form of a covariance matrix. Our method learns
heteroscedastic aleatoric uncertainty from ICP input data and uses Bayesian
posterior approximation to capture the epistemic uncertainty.

Figure 1: Distribution of errors on a series of ICP estimates computed in
three different kind of environments: corridor-like (left), wide plain
(center), structured (right).

The following section introduces the ICP process, reviews the various sources
of errors and works relevant to ours: methods that estimate ICP covariance,
data-driven approaches that have proven successful to estimate covariances for
other estimation processes, and point-cloud based deep learning techniques.
Section III formalizes the problem and depicts the learning process. Section
IV evaluates our results against the state-of-the-art, showing an accurate
correlation between the real error of ICP and our predicted covariance.

##  II Background and related work

###  II-A Iterative Closest Point process

ICP tackles the problem of aligning a _reading_ point cloud
$\mathbfcal{P}\in\mathbb{R}^{3\times n}$ to a _reference_ point cloud
$\mathbfcal{Q}\in\mathbb{R}^{3\times m}$, estimating the true rigid
transformation between them. Using an initial estimation $\check{{\bm{T}}}\in$
SE(3) of the transformation (_e.g_ as provided by a prior on the motion), the
rigid transformation $\hat{{\bm{T}}}\in$ SE(3) that estimates the true
transformation is the one that minimizes the error function $e$:

$\hat{{\bm{T}}}=\operatorname*{arg\,min}_{{\bm{T}}}(e({\bm{T}}(\mathbfcal{P}),\mathbfcal{Q}))$
(1)

where ${\bm{T}}(\mathbfcal{C})$ is the transformation ${\bm{T}}$ applied to
the point cloud $\mathbfcal{C}$. The error function is the result of a
weighted average distance between appropriately matched point pairs or point
to plane associations. The optimization is an iterative process, and there is
no guarantee to reach the optimal solution.

###  II-B Sources of error for ICP

There are five main sources of error for ICP described in the literature:
wrong convergence, wrong initialization, sensor noise and bias,
underconstrained situations, and intrinsic randomness of the process [5, 6].

#### Convergence and initialization

Wrong convergence comes from the iterative nature of ICP, which does not
guarantee global convergence and can get stuck in local minima. This is
closely related to wrong initialization, which can be a predominant source of
error when the initial guess is far from the true solution [6].

#### Sensor noise and biases

Each sensor naturally suffers changes in its accuracy that depend on several
environmental conditions (_e.g_ temperature) [7]. Also, sensor calibration can
introduce distortion [8]. While these phenomena introduce a bias in the
estimate [9], there is an additional source of sensor noise added to each
point independently, which generally is a function of the distance,
orientation and physical nature or texture of the perceived material.

#### Underconstrained situations

These occur when the environment does not offer enough information to fully
estimate the transformation between the point clouds, mainly by precluding the
establishment of good matches between points. While in the planar 2D case such
situations can be exhaustively identified, namely the corridor and the circle,
in the 3D world this becomes a more challenging task.

#### Intrinsic ICP randomness

Depending on the configuration and implementation, a certain degree of
randomness is introduced by subprocesses of ICP, albeit the process itself is
deterministic. For example, filtering clouds with outlier rejection and sub-
sampling can yield different estimates, even with the same input data.

###  II-C Estimating ICP covariance

Various methods for estimating ICP covariance have been proposed. Monte-Carlo
algorithms [10] generate the covariance from a large number of samples of scan
registration results. Different clouds and initializations are randomly
subsampled to obtain a dispersion of the transformation that is used to
compute the covariance. Closed-form methods [5, 6, 11, 12] provide algebraic
solutions to the ICP covariance problem by linearizing the error functions
defined in the ICP algorithm. These methods are potentially ill-founded if
point reassociations happen at a scale smaller than the actual error [13].
Experiments on synthetic data have shown that this is often the case [14].
Early closed-form works, as the one in [5], have been demonstrated over
optimistic [2]. Recently, [14] proposed a learning approach to the ICP
covariance estimation problem. The work estimates a 3D covariance based on
preprocessed feature descriptors that are used in the learning scheme. While
such descriptors are generic, it is desirable to have a more flexible
representation of features impacting on the error of ICP.

It is worth to notice that most of these work tackle different sources of
error. The method presented in [6] shows the need for a proper initialization
of ICP and proposes a novel approach to 3D uncertainty of ICP accounting for
initialization errors. Censi’s pioneering work [5] studies the error due to
sensor noise and computes the covariance of this error as a function of the
error metric minimized by any ICP implementation. The work in [14], closer to
our approach, explores a data-driven approach to estimate the covariance
focusing on sensor noise and aiming at detecting underconstrained situations.

###  II-D Learning motion estimation error models

The literature in deep learning features a large amount of work in the
estimation of uncertainty. The approaches are often tied to the process for
which they predict error models. Interesting works in domains similar to ours
have been proposed, mainly tackling the Visual Odometry (VO) case. As ICP, VO
is a well-known family of techniques which estimate a rigid transformation
between two consecutive frames for the monocular or stereovision cases. End-
to-end methods perform frame to frame motion estimation [15] and covariance
prediction [16]. Both estimates are in this case delegated to neural networks,
and do not rely on any geometrical pipeline. The works in [17, 18]
respectively show covariance and joint bias-covariance estimation to the
geometric implementation of VO.

###  II-E Deep learning for point clouds

Recently, deep learning methods applied on point clouds have proven successful
for several tasks [19, 20, 21, 22]. In particular, PointNet [19] paved the way
for neural networks that directly use 3D point clouds for object
classification and segmentation purposes, and has been extended towards
learning point cloud features in metric space [20]. This extension increased
the effectiveness of PointNet to learn features on different scales, in the
manner of traditional convolutional networks. Other approaches rely on the
spherical projection of point clouds, usually coined as range (or depth) image
– which is, in most cases, the original structure of the sensor data. This is
used in the LiDAR odometry techniques [23, 24, 25] which makes use of the
entire scan to estimate the motion between two frames in an end-to-end manner.

##  III Data driven learning of ICP uncertainty

Excluding the process initialization, all error sources of ICP are inherited
from the input data. Sensor noise is directly reflected in the data, as well
as ICP subprocesses modifying its size, density and structure (_i.e_
randomness due to filtering). Additionally, the scene structure, with its
inherent geometry and features, largely impacts the error of ICP (as
illustrated in Fig. 1), as it defines the observability of the motion along
the different dimensions.

Addressing all these sources in a unified framework is a challenging task.
Leveraging recent progresses in deep learning for point clouds, we address the
problem of estimating the covariance of ICP that stems from the last three
error sources presented in Sec. II-B. Our hypothesis is that a learning scheme
that processes the raw input data can capture these error sources and properly
estimate ICP covariances.

In this section we describe our approach to estimate ICP covariances and
discuss the nature of the learned uncertainties. The optimization problem is
contextualized in a probabilistic framework and a description of the learning
architecture incorporating different types of uncertainty is provided.

###  III-A Minimization problem

Modeling errors for an estimation process can be tackled via a supervised
approach. Consider an estimator, ICP in our case, using sensory data
$\mathbfcal{C}_{i}=\\{\mathbfcal{P}_{i},\mathbfcal{Q}_{i+1}\\}$ with
$\mathbfcal{P}_{i}\in\mathbb{R}^{n_{i}}$ and
$\mathbfcal{Q}_{i+1}\in\mathbb{R}^{n_{i+1}}$. Its output
${}^{i}\hat{{\bm{T}}}_{i+1}$ is an estimate of the ground truth transformation
${}^{i}{\bm{T}}_{i+1}$. In practice, the point clouds used in the algorithm
are the result of several filters [4]. The actual input data is a pair of
decimated point clouds
$\bar{\mathbfcal{C}}_{i}=\\{\bar{\mathbfcal{P}}_{i},\bar{\mathbfcal{Q}}_{i+1}\\}$
with $\bar{\mathbfcal{P}}_{i}\in\mathbb{R}^{m_{i}},m_{i}\ll n_{i}$ and
$\bar{\mathbfcal{Q}}_{i+1}\in\mathbb{R}^{m_{i+1}},m_{i+1}\ll n_{i+1}$. The
error of ICP is computed as

${\bm{e}}_{i}\ =\ ^{i}\hat{{\bm{T}}}^{-1}_{i+1}\ \cdot\ ^{i}{\bm{T}}_{i+1}$
(2)

We represent an error ${\bm{e}}\in\mathrm{SE(3)}$ as a vector
${\bm{\xi}}\in\mathfrak{se}(3)$ member of the Lie algebra [26, 27]

${\bm{\xi}}^{\wedge}=\begin{bmatrix}{\bm{\rho}}\\\
{\bm{\phi}}\end{bmatrix}^{\wedge}=\
\begin{bmatrix}{\bm{\phi}}^{\wedge}&{\bm{\rho}}\\\
\mathbf{0}^{\mathsf{T}}&0\end{bmatrix}$ (3)

where the _hat_ operator $\wedge$ turns the vector
${\bm{\xi}}\in\mathbb{R}^{6}$ into a matrix
${\bm{\xi}}^{\wedge}\in\mathbb{R}^{4\times 4}$. The opposite linear map is
represented as the _vee_ operator $\vee$. The exponential map
$\mathrm{exp}({\bm{\xi}}^{\wedge})$ allows to retract an element of the Lie
algebra back to the group. Its inverse is the logarithmic map
$\mathrm{log}({\bm{T}})$. Once the error is defined as
${\bm{\xi}}=\mathrm{log}({\bm{e}})^{\vee}$, the associated uncertainty
${\bm{\Sigma}}\in\mathbb{R}^{6\times 6}$ can be expressed as

${\bm{\Sigma}}=\begin{bmatrix}{\bm{\Sigma}}_{{\bm{\rho}}{\bm{\rho}}}&{\bm{\Sigma}}_{{\bm{\rho}}{\bm{\phi}}}\\\
{\bm{\Sigma}}_{{\bm{\phi}}{\bm{\rho}}}&{\bm{\Sigma}}_{{\bm{\phi}}{\bm{\phi}}}\end{bmatrix}$
(4)

This uncertainty representation is generally valid for small perturbations
that are added to a given pose
$\hat{{\bm{T}}}={\bm{T}}\mathrm{exp}({\bm{\xi}})$ where
${\bm{\xi}}\sim\mathcal{N}(\mathbf{0},{\bm{\Sigma}})$[27]. The equation above
matches the error definition in Eq. 2, that is ${\bm{\xi}}\in\mathbb{R}^{6}$
can be viewed as the error between the real transformation and its estimate.

We define a dataset
$\mathcal{D}=\\{{\bar{\mathbfcal{C}}}_{i},{\bm{\xi}}_{i}|\forall
i\in[1,d]\\}$, where $d$ is the size of the dataset. The goal is to predict
the uncertainty of the error vector ${\bm{\xi}}_{i}$ along its dimensions. In
the data, we noticed that the unbiased error assumption is valid for the vast
majority of ICP estimates and the bias is very small otherwise. To keep the
uncertainty assumption valid, we encode a bias term for the error vector in
the distribution that ensures the error remains close to a zero-mean Gaussian.
Therefore, it is possible to estimate the parameters of the Gaussian
$\mathcal{N}({\bm{\mu}},{\bm{\Sigma}})$ as

$\operatorname*{arg\,max}_{{\bm{\mu}}_{1:d},{\bm{\Sigma}}_{1:d}}\sum_{i=1}^{d}p({\bm{\xi}}_{i}|{\bm{\mu}}_{i},{\bm{\Sigma}}_{i})$
(5)

To maximize the probability of drawing the correct error given the estimated
Gaussian parameters is equivalent to minimize the negative full log-likelihood
function [18]. To ensure a positive definite covariance matrix we use the LDL
decomposition as in [17].

###  III-B Uncertainty estimation

ICP covariance can be decomposed into two categories: epistemic and aleatoric
[28]. Both uncertainties are formulated as probability distributions over
different entities.

Aleatoric uncertainty models the uncertainty over the input data to the
process. It assumes an observation noise that can be either constant or vary
with the input. As discussed in Sec.II-B, we observed that noise levels change
with the scene (Fig. 1), hence we are interested in modeling the variation of
the noise in correlation with the input data. This is referred to as
heteroscedastic aleatoric uncertainty (as opposed to homoscedastic uncertainty
for the constant noise case).

Additionally, it is possible to model the variance intrinsic to the learned
model, that is how certain of a given prediction the network is. Epistemic
uncertainty can, in principle, be reduced if more data is used for training.
In our case, it accounts for the network’s poor predictions on out-of-data
samples such as scene features that were rarely encountered during training.
Epistemic uncertainty can be modeled by placing a distribution over the
network’s weights. Unfortunately, given a prior distribution
$\mathcal{N}_{{\bm{W}}}(0,\sigma)$ over the weights, it is computationally
difficult to analytically infer the posterior $p({\bm{W}}|{\bm{X}},{\bm{Y}})$,
that is the probability of drawing a set of weights given the network input
and output. There is a number of proposed methods to approximate this
posterior [29, 30]. The work in [31] produces an approximation using dropout
variational inference. This has proven to be an effective yet easily
obtainable Bayesian approximation when dealing with deep complex models.

In our context, From Eq.5, we consider the mean as a noisy error estimate, and
the covariance matrix as the aleatoric uncertainty due to the observation
noise. To approximate the posterior for a Bayesian neural network estimating
the ICP errors we resort to Monte-Carlo dropout. The epistemic uncertainty can
be captured by keeping dropout activated at test time. Sampling $N$ times the
predicted ICP error vector and ICP covariance yields a distribution
$[\hat{{\bm{\mu}}}_{1:N},\hat{{\bm{\Sigma}}}_{1:N}]=f_{{\bm{W}}}(\bar{\mathbfcal{C}})$
where $f_{{\bm{W}}}(\bar{\mathbfcal{C}})$ is the output of a Bayesian neural
network, function of the input point clouds.

The total covariance matrix, integrating aleatoric and epistemic uncertainty
is computed as

$\frac{1}{N}\sum_{n=1}^{N}(\hat{{\bm{\mu}}}_{n}-\bar{{\bm{\mu}}})(\hat{{\bm{\mu}}}_{n}-\bar{{\bm{\mu}}})^{\mathsf{T}}+\frac{1}{N}\sum_{n=1}^{N}\hat{{\bm{\Sigma}}}_{n}$
(6)

where $\bar{{\bm{\mu}}}=\frac{1}{N}\sum_{n=1}^{N}\hat{{\bm{\mu}}}_{n}$ is the
mean error vector over the sampled outputs. The first term in Eq. 6 computes
the uncertainty over the predicted error vector, which presents a variance
proportional to the network confidence over the estimate. The mean vector is
not used to correct ICP estimates as the original error is generally noisy. It
is rather evaluated in terms of dispersion after having approximated the model
posterior.

###  III-C Learning architecture

The improvements in point clouds based network involved a number of estimation
processes. Recently, FlowNet3D [32] proposed an approach to tackle 3D scene
flow expanding on the building blocks introduced in [20]. It provides the 3D
flow for each point in a reference point cloud deducing its motion from the
comparison with a second, successive, point cloud. Estimating the scene flow
is a problem closely related to the scan registration. Intuitively, a network
able to approximate such task results suitable for ICP related inferences. In
fact, the main difference between the two tasks lies in the generalization
needed to shift from a high dimensional flow to the estimate of a rigid
transformation in the ICP case. We leverage the main structure of FlowNet3D to
tackle the estimation of ICP errors, modifying few key aspects to adapt it to
our problem. To complete the architecture we add a final MLP layer reducing
the dimensionality to the number of parameters needed to predict the Gaussian
parameters. First we follow the original regression layer from FlowNet3D
reducing the feature space to $\mathbb{R}^{3}\times n$. The MLP is extended to
take as input the feature produced by the regression layer and to return the
parameters generating the covariance matrix and predicted error vector. All
the MLP (originally Linear-BatchNorm-ReLU) in FlowNet3D are stacked with
Dropout modules to allow the estimation of the epistemic covariance matrix.
Additionally, this proved to reduce overfitting as expected.

We initially assumed that including the set upconv layer would have not been
needed. Intuitively, as in a classical convolutional scheme, the regression
task we try to solve can be solved with direct dimensionality reduction, from
sensory data to parameters describing uncertainty [18, 17]. Nonetheless,
retaining upsampling layer produced better results in terms of predicted
uncertainty and mean. We associate this behavior with the architecture using
skip connections to use the same features extracted from the same point sets
originally used ICP. This creates a 1:1 matching between the data used by the
geometric estimator and the network, leaving less room for outliers. Our
network parameterization, with sampling rates and radius, can be found in
Table I.

Layer type | Radius | Sampling | MLP size
---|---|---|---
set conv | 2.0 |  1.0$\times$ | [32, 32, 64]
set conv | 4.0 |  0.25$\times$ | [64, 64, 128]
flow embedding | knn |  1.0$\times$ | [128, 128, 128]
set conv | 8.0 |  0.25$\times$ | [128, 128, 256]
set conv | 16.0 |  0.125$\times$ | [256, 256, 512]
TABLE I: Network architecture parameterization. Set upconv layers are the same
as in the updated version of [32].

##  IV Evaluation

###  IV-A ICP process and dataset

To learn errors produced by a geometric ICP algorithm we selected
_libpointmatcher_ [33]. It is a widely used open-source library offering a
chain that filters point clouds, establishes points associations, removes
outliers and estimates rigid transformations using an error minimizer. For our
setup we parameterized the process similarly to [14]. Cloud random sampling is
fixed at 2048 random points. We generate normals, used for the point-to-plane
minimizer, using the 10 nearest points before subsampling.

To select the initial transform for ICP, we add an error to the ground truth
${\bm{T}}$. In order not to draw errors from an overoptimistic distribution we
sample guesses for the initial transforms from a non-zero mean normal
distribution
$\mathcal{N}(\mathcal{N}({\bm{T}},a\cdot{\bm{T}}),b\cdot{\bm{T}})$, with
$a=0.25,b=0.2$. These figures yield more pessimistic estimates than commonly
used aiding sensor (_e.g_ odometry or inertial measurements) posing a
challenge to our system.

We select the KITTI odometry dataset [34], using LiDAR scans, which is made of
urban and rural road scenes. We train and test using all ground truth tagged
sequences but sequence 01, on which ICP fails to register most of the scans
due to scenes lacking features. To compare against the state of the art in
covariance estimation for ICP, we also present results on the challenging data
sets of [35]. Since for the LiDAR odometry case this dataset is rather small,
we refine our KITTI trained model with a few epochs of training, approximately
30 or less for early convergence, on all but one outdoor sequences and
validate on the one left out. As expected, features learned on KITTI are
ineffective to predict covariances on the _apartment_ and _stairs_ sequences,
which are highly structured small indoor environments. Our system learns to
model ICP errors with numerically small epistemic uncertainty. For reasons of
space, we only show results containing the combined aleatoric+epistemic
covariance when moving to unseen datasets, _i.e_ applying models learned on
KITTI to the ETH dataset.

###  IV-B Single-pair validation

We first assess the quality of the predicted covariances in a single-pair
fashion, and then check their consistence on a full trajectory. The metrics we
use are based on the relation between the observed errors and the predicted
covariance. This is due to the absence of the true distribution which would
allow for more classic figures of merit, such as the Kullback-Leibler
divergence. While it would be possible to obtain a pseudo-true distribution by
random sampling over the point cloud, this would affect the estimated
covariance as well, because the network prediction is tied to the ICP
configuration, and most importantly to the input data used by ICP.

As in [5, 6], we compute the _Normalized Norm Error (NNE)_ to evaluate the
scale of the predicted covariance w.r.t the actual error:

$\mathrm{NNE}=\frac{1}{N}\sum_{n=1}^{N}\sqrt{\frac{||{\bm{\xi}}_{n}||_{2}^{2}}{tr(\hat{{\bm{\Sigma}}}_{n})}}$
(7)

where $tr(\hat{{\bm{\Sigma}}}_{n})$ is the trace of the predicted covariance
matrix. The optimal value for NNE is one. Values below and above one
respectively represent pessimistic and optimistic uncertainty estimates. Table
II shows NNE values computed over test KITTI sequences. The scale of the
covariance is well predicted, and is only slightly pessimistic for two
validation sequences (05 and 07).

NNE (KITTI) | 05 | 06 | 07 | 09 | 10
---|---|---|---|---|---
Translation | 0.70 | 0.61 | 0.58 | 0.94 | 0.90
Rotation | 0.65 | 0.81 | 0.55 | 0.90 | 1.04
TABLE II: NNE translation and rotation values for the KITTI dataset.

The NNE is an immediate measure of the adequacy of the predicted variances,
but it does not account for off-diagonal covariances. For this purpose, the
Mahalanobis distance between ICP estimates and ground truth
${\bm{\xi}}_{n}=\mathrm{log}(^{n}\hat{{\bm{T}}}^{-1}_{n+1}\
{}^{n}{\bm{T}}_{n+1})^{\vee}$, weighted by the predicted covariance
${\bm{\Sigma}}_{n}$ averaged over each consecutive pair in a sequence, is a
better metric:

$\mathrm{D_{M}}=\frac{1}{N}\sum_{n=1}^{N}\sqrt{\frac{{\bm{\xi}}_{n}^{\mathsf{T}}{\bm{\Sigma}}_{n}^{-1}{\bm{\xi}}_{n}}{dim({\bm{\xi}})}}$
(8)

As in the NNE case, the optimal value for a normal distribution is one, with
distances below and above one respectively denoting pessimistic and optimistic
estimates.

_Single-Pair_ Mah. Dist (KITTI) | 05 | 06 | 07 | 09 | 10
---|---|---|---|---|---
Translation | 0.80 | 0.83 | 0.74 | 0.88 | 1.03
Rotation | 0.77 | 0.63 | 0.70 | 0.94 | 0.96
TABLE III: Single-pair Mahalanobis distance in terms of translation and
rotation for the KITTI dataset.

Table III shows the average Mahalanobis distance for the validation sequences.
All sequences present values positively close to 1, indicating that the system
is able to correctly capture the uncertainty under the Gaussian assumption.
For these sequences, the epistemic uncertainty is small and constant, which
denotes confidence in the system prediction. Minor variations are encountered
for seldom cases: an example is shown in Fig. 2 displaying the epistemic
covariance before, during, and after an area corresponding to _wide plain_
scenes shown in Fig. 1. Such spatial features are less frequently encountered
during training.

Figure 2: Epistemic and aleatoric $3\sigma$ intervals in a rarely seen area
lacking sharp geometrical features. The real error computed from the ground
truth is represented in black. The zoomed-in area shows epistemic uncertainty
increase, which corresponds to a wide plain environment.

###  IV-C Trajectory validation

In order to assess the quality of our predicted covariances on trajectories we
take a multi-step approach. First, we compose ICP estimates to generate the
final pose as ${}^{0}\hat{{\bm{T}}}_{f}\ =\ ^{0}\hat{{\bm{T}}}_{1}\
^{1}\hat{{\bm{T}}}_{2}\cdots\ ^{f-1}\hat{{\bm{T}}}_{f}$ for a sequence with
$f$ poses. The ground truth final pose ${}^{0}{\bm{T}}_{f}$ is computed
similarly. We propagate the covariance associated to each estimate using the
4th order approximation, as in [27], to obtain an uncertainty estimate for the
final ICP pose. This way, it is possible to evaluate the Mahalanobis distance
(Eq. 8) in a spatially consistent manner. The results are shown in Table IV.
They are only slightly optimistic, which is a remarkable outcome given the
length of KITTI trajectories (up to more than 2km).

Figure 3: Four propagated covariances in the 2$\sigma$ interval (95%) on considered ICP trajectories compared to ground truth. _Trajectory_ Mah. Dist (KITTI) | 05 | 06 | 07 | 09 | 10
---|---|---|---|---|---
Translation | 2.31 | 1.51 | 1.26 | 1.36 | 4.86
Rotation | 1.49 | 1.95 | 2.19 | 1.89 | 2.95
TABLE IV: Trajectory Mahalanobis distance in terms of translation and rotation
for the KITTI dataset.

###  IV-D Fine tuning on other datasets

We also evaluate the results on the smaller dataset of [35] and compare to the
state of the art. As our model requires a significant amount of data we
initially learn on a large dataset. KITTI makes available a considerable
number of LiDAR scans with different scenarios, but all correspond to outdoor
road environments. To apply our approach to smaller datasets acquired in very
different environments, a fine-tuning process is applied using the pre-trained
model on the KITTI dataset. This allows to obtain much better results compared
to re-training on such a small dataset and shows the importance of epistemic
uncertainty as an indicator of the performances of a model on unseen data.

| Single-Pair
---
Mahalanobis Distance (ETH)
| Gazebo
---
winter
| Gazebo
---
summer
| Mountain
---
| Haupt-
---
-gebaude
| Wood
---
autumn
| Wood
---
summer
| Trans (A) | 2.05 | 2.32 | 4.15 | 2.68 | 4.01 | 4.03
Ours (Direct) | Trans (E+A) | 2.04 | 2.12 | 3.88 | 2.49 | 4.19 | 3.78
| Trans (A) | 0.64 | 1.12 | 1.27 | 1.04 | 1.86 | 1.03
Ours (Finetuned) | Trans (E+A) | 0.59 | 1.06 | 1.21 | 0.98 | 1.70 | 0.90
| Rot (A) | 1.33 | 1.93 | 3.06 | 2.61 | 2.10 | 2.33
Ours (Direct) | Rot (E+A) | 1.30 | 1.89 | 3.03 | 2.47 | 2.07 | 2.19
| Rot (A) | 0.48 | 0.53 | 0.65 | 1.02 | 0.54 | 0.82
Ours (Finetuned) | Rot (E+A) | 0.44 | 0.49 | 0.65 | 0.96 | 0.51 | 0.71
TABLE V: Single-pair Mahalanobis distance in terms of translation and rotation for the ETH dataset using aleatoric only (A) and combined epistemic-aleatoric (E+A) uncertainty. A comparison between the pre-trained model on KITTI and after finetuning is presented. After finetuning, the epistemic uncertainty is less prominent and its (minor) impact on performances depends on whether the system is being pessimistic or optimistic. | Trajectory
---
Mahalanobis Distance (ETH)
| Gazebo
---
winter
| Gazebo
---
summer
| Mountain
---
| Haupt-
---
-gebaude
| Wood
---
autumn
| Wood
---
summer
CELLO-3D [14] | Trans | 0.1 | 0.2 | - | 0.3 | 0.1 | 0.1
_Brossard et al._ [6] | Trans | 1.8 | 1.0 | 1.2 | 1.8 | 1.2 | 1.5
Ours (Finetuned) | Trans (A) | 0.83 | 1.28 | 0.87 | 1.83 | 0.93 | 0.69
CELLO-3D [14] | Rot | 0.2 | 0.2 | - | 0.2 | 0.3 | 0.3
_Brossard et al._ [6] | Rot | 3.7 | 2.3 | 1.2 | 2.9 | 4.2 | 4.7
Ours (Finetuned) | Rot (A) | 0.71 | 1.03 | 0.81 | 1.48 | 1.69 | 1.65
TABLE VI: Trajectory Mahalanobis distance in terms of translation and rotation
for the ETH dataset using aleatoric only (A) uncertainty, compared to two
state-of-the-art approaches. Epistemic uncertainty is not showed after
finetuning.

Table V shows the result of our approach on single-pair estimates with the
dataset of [35]. While the environments are very different, the model performs
remarkably well in outdoor scenes. Even architectural features, such as in
_hauptgebaude_ can be leveraged to produce faithful covariance predictions as,
intuitively, the scene presents a structure (flat ground plane, lack of indoor
man-made elements) that is geometrically similar to the some scenes
encountered in KITTI.

Figure 4: Epistemic uncertainties normalized by the mean in the direct and
finetuned case.

More importantly, Table V highlights the impact of the finetuning, not only in
the expected improved results, but also in the magnitude of the epistemic
uncertainty. The improvement brought by considering the model uncertainty is
evident when testing the model trained on KITTI directly on an unseen dataset,
where, as expected, the network has little confidence in its prediction for
new scenarios. In fact, while some outdoor features may recur in the ETH
dataset, motion type and scene size are remarkably different when compared to
the KITTI dataset. After finetuning, not only the weight of the epistemic
covariance is greatly reduced as shown in Fig. 4, but the overall results are
closer to the optimal value.

Finally, Table VI compares our results against CELLO-3D [14] and Brossard et
al. [6] on aggregated transforms on the same dataset. It is worth noticing
that given the strong focus on ICP initialization from Brossard et al., it is
more relevant to compare our approach to CELLO-3D, as our work solely focus on
the errors introduced directly and indirectly by the data on a purely data-
driven approach. CELLO-3D is often over pessimistic while _Brossard et al._ is
generally optimistic. Especially on rotational uncertainty, our method
outperforms the two approaches by a significant margin.

##  V Conclusion

We presented an approach to estimate faithful covariances for the ICP scan
registration process. Using a data-driven paradigm we have shown that it is
possible to learn uncertainty for the selected algorithm, and to generalize to
a certain degree on multiple datasets. Leveraging a neural network
architecture conceived to directly process point clouds, we estimate
covariances with a close connection to the data used in the ICP algorithm.
Standard metrics show the reliability of our approach, and that it generally
outperforms the state-of-the-art. Future improvements will concern the network
architecture, possibly tailoring it to the scan registration problem. A
possible approach is to have the flow embedding layer (and the successive
layers) doubled to account for both point clouds motions, whereas in the scene
flow the problem is cast on the reference cloud only. Additionally, the role
of initialization and the use of the initial guess in the learning process can
also be considered.

## References

  * [1] J. Zhang and S. Singh, “LOAM: Lidar Odometry and Mapping in Real-time.,” in Robotics: Science and Systems, 2014.
  * [2] E. Mendes, P. Koch, and S. Lacroix, “ICP-based pose-graph SLAM,” in 2016 IEEE International Symposium on Safety, Security, and Rescue Robotics (SSRR), pp. 195–200, IEEE, 2016.
  * [3] P. Besl and N. McKay, “A Method for Registration of 3-D Shapes,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 14, no. 2, pp. 129–256, 1992.
  * [4] F. Pomerleau, F. Colas, and R. Siegwart, “A review of point cloud registration algorithms for mobile robotics,” Foundations and Trends® in Robotics, vol. 4, no. 1, pp. 1–104, 2015.
  * [5] A. Censi, “An accurate closed-form estimate of ICP’s covariance,” in Proceedings of IEEE International Conference on Robotics and Automation (ICRA), pp. 3167–3172, IEEE, 2007.
  * [6] M. Brossard, S. Bonnabel, and A. Barrau, “A new approach to 3D ICP covariance estimation,” IEEE Robotics and Automation Letters, vol. 5, no. 2, pp. 744–751, 2020.
  * [7] F. Pomerleau, A. Breitenmoser, M. Liu, F. Colas, and R. Siegwart, “Noise characterization of depth sensors for surface inspections,” in 2012 2nd International Conference on Applied Robotics for the Power Industry (CARPI), pp. 16–21, IEEE, 2012.
  * [8] J.-E. Deschaud, “IMLS-SLAM: scan-to-model matching based on 3D data,” in 2018 IEEE International Conference on Robotics and Automation (ICRA), pp. 2480–2485, IEEE, 2018.
  * [9] J. Laconte, S.-P. Deschênes, M. Labussière, and F. Pomerleau, “Lidar measurement bias estimation via return waveform modelling in a context of 3D mapping,” in 2019 International Conference on Robotics and Automation (ICRA), pp. 8100–8106, IEEE, 2019.
  * [10] A. G. Buch, D. Kraft, et al., “Prediction of ICP pose uncertainties using Monte Carlo simulation with synthetic depth images,” in 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 4640–4647, IEEE, 2017.
  * [11] S. M. Prakhya, L. Bingbing, Y. Rui, and W. Lin, “A closed-form estimate of 3D ICP covariance,” in 2015 14th IAPR International Conference on Machine Vision Applications (MVA), pp. 526–529, IEEE, 2015.
  * [12] N. Gelfand, L. Ikemoto, S. Rusinkiewicz, and M. Levoy, “Geometrically stable sampling for the icp algorithm,” in Fourth International Conference on 3-D Digital Imaging and Modeling, 2003. 3DIM 2003. Proceedings., pp. 260–267, IEEE, 2003.
  * [13] S. Bonnabel, M. Barczyk, and F. Goulette, “On the covariance of ICP-based scan-matching techniques,” in 2016 American Control Conference (ACC), pp. 5498–5503, IEEE, 2016.
  * [14] D. Landry, F. Pomerleau, and P. Giguere, “CELLO-3D: Estimating the Covariance of ICP in the Real World,” in 2019 International Conference on Robotics and Automation (ICRA), pp. 8190–8196, IEEE, 2019.
  * [15] S. Wang, R. Clark, H. Wen, and N. Trigoni, “DeepVO: Towards end-to-end visual odometry with deep recurrent convolutional neural networks,” in Robotics and Automation (ICRA), 2017 IEEE International Conference on, pp. 2043–2050, IEEE, 2017.
  * [16] S. Wang, R. Clark, H. Wen, and N. Trigoni, “End-to-end, sequence-to-sequence probabilistic visual odometry through deep neural networks,” The International Journal of Robotics Research, vol. 37, no. 4-5, pp. 513–542, 2018\.
  * [17] K. Liu, K. Ok, W. Vega-Brown, and N. Roy, “Deep inference for covariance estimation: Learning Gaussian noise models for state estimation,” in 2018 IEEE International Conference on Robotics and Automation (ICRA), pp. 1436–1443, IEEE, 2018.
  * [18] A. De Maio and S. Lacroix, “Simultaneously Learning Corrections and Error Models for Geometry-based Visual Odometry Methods,” IEEE Robotics and Automation Letters, vol. 5, no. 4, pp. 6536–6543, 2020.
  * [19] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “Pointnet: Deep learning on point sets for 3d classification and segmentation,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 652–660, 2017.
  * [20] C. R. Qi, L. Yi, H. Su, and L. J. Guibas, “Pointnet++: Deep hierarchical feature learning on point sets in a metric space,” Advances in Neural Information Processing Systems (NeurIPS), vol. 30, 2017.
  * [21] H. Thomas, C. R. Qi, J.-E. Deschaud, B. Marcotegui, F. Goulette, and L. J. Guibas, “Kpconv: Flexible and deformable convolution for point clouds,” in Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), pp. 6411–6420, 2019.
  * [22] C. R. Qi, O. Litany, K. He, and L. J. Guibas, “Deep Hough voting for 3D object detection in point clouds,” in Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 9277–9286, 2019.
  * [23] Y. Cho, G. Kim, and A. Kim, “Deeplo: Geometry-aware deep lidar odometry,” arXiv preprint arXiv:1902.10562, 2019.
  * [24] Z. Li and N. Wang, “Dmlo: Deep matching lidar odometry,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 6010–6017, IEEE, 2020.
  * [25] Q. Li, S. Chen, C. Wang, X. Li, C. Wen, M. Cheng, and J. Li, “Lo-net: Deep real-time lidar odometry,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 8473–8482, 2019.
  * [26] J. Sola, J. Deray, and D. Atchuthan, “A micro Lie theory for state estimation in robotics,” arXiv preprint arXiv:1812.01537, 2018.
  * [27] T. D. Barfoot and P. T. Furgale, “Associating uncertainty with three-dimensional poses for use in estimation problems,” IEEE Transactions on Robotics, vol. 30, no. 3, pp. 679–693, 2014.
  * [28] A. Kendall and Y. Gal, “What uncertainties do we need in bayesian deep learning for computer vision?,” in Advances in neural information processing systems, pp. 5574–5584, 2017.
  * [29] M. Stephan, M. D. Hoffman, D. M. Blei, et al., “Stochastic Gradient Descent as Approximate Bayesian Inference,” Journal of Machine Learning Research, vol. 18, no. 134, pp. 1–35, 2017.
  * [30] J. Postels, F. Ferroni, H. Coskun, N. Navab, and F. Tombari, “Sampling-free epistemic uncertainty estimation using approximated variance propagation,” in Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), pp. 2931–2940, 2019.
  * [31] Y. Gal and Z. Ghahramani, “Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning,” in Proceedings of the 33rd International Conference on Machine Learning (ICML), 2016.
  * [32] X. Liu, C. R. Qi, and L. J. Guibas, “FlowNet3D: Learning scene flow in 3D point clouds,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pp. 529–537, 2019.
  * [33] F. Pomerleau, S. Magnenat, F. Colas, M. Liu, and R. Siegwart, “Tracking a Depth Camera: Parameter Exploration for Fast ICP,” in Proc. of the IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 3824–3829, IEEE Press, 2011.
  * [34] A. Geiger, P. Lenz, and R. Urtasun, “Are we ready for Autonomous Driving? The KITTI Vision Benchmark Suite,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pp. 3354–3361, IEEE, 2012.
  * [35] F. Pomerleau, M. Liu, F. Colas, and R. Siegwart, “Challenging data sets for point cloud registration algorithms,” The International Journal of Robotics Research, vol. 31, no. 14, pp. 1705–1711, 2012.


# Ballistic-to-hydrodynamic transition and collective modes
for two-dimensional electron systems in magnetic field

Kirill Kapralov1,2<EMAIL_ADDRESS>Dmitry Svintsov1 1Center for
Photonics and 2d Materials, Moscow Institute of Physics and Technology
(National Research University), Dolgoprudny 141700, Russia.
2Kotelnikov Institute of Radioengineering and Electronics of Russian Academy
of Sciences, Mokhovaya 11-7, Moscow 125009, Russia

###### Abstract

The recent demonstrations of viscous hydrodynamic electron flow in two-
dimensional electron systems poses serious questions to the validity of
existing transport theories, including the ballistic model, the collision-
induced and collisionless hydrodynamics. While the theories of transport at
hydrodynamic-to-ballistic crossover for free 2d electrons are well
established, the same is not true for electrons in magnetic fields. In this
work, we develop an analytically solvable model describing the transition from
ballistic to hydrodynamic transport with changing the strength of electron-
electron collisions in magnetic fields. Within this model, we find an
expression for the high-frequency non-local conductivity tensor of 2d
electrons. It is valid at arbitrary relation between frequency of external
field $\omega$, the cyclotron frequency $\omega_{c}$, and the frequency of e-e
collisions $\tau^{-1}_{ee}$. We use the obtained expression to study the
transformation of 2d magnetoplasmon modes at hydrodynamic-to-ballistic
crossover. In the true hydrodynamic regime, $\omega\tau_{ee}\ll 1$, the 2DES
supports a single magnetoplasmon mode that is not split at cyclotron
harmonics. In the ballistic regime, $\omega\tau_{ee}\gg 1$, the plasmon
dispersion develops splittings at cyclotron harmonics, forming the Bernstein
modes. A formal long-wavelength expansion of kinetic equations (”collisionless
hydrodynamics”) predicts the first splitting of plasmon dispersion at
$\omega\approx 2\omega_{c}$. Still, such expansion fails to predict the zero
and negative group velocity sections of true magnetoplasmon dispersion, for
which the full kinetic model is required.

††preprint: APS/123-QED

##  I  Introduction

In common semiconductor materials, scattering of electrons by phonons and
impurities leads to the diffusive Ohmic transport. At the same time, in pure
materials, the characteristic momentum relaxation length $l_{mr}$ for
scattering by phonons and impurities can exceed the size of the device $W$.
This leads to ballistic dc transport with predominant scattering of electrons
at the sample boundaries Mayorov  _et al._ (2011a); Heremans  _et al._
(1992a). A very interesting situation appears when momentum-conserving
electron-electron collisions are intense, so that the respective free path
$l_{ee}\ll l_{mr}$ and $l_{ee}\ll W$. In this case, electrons form a viscous
fluid with its dynamics described by the Navier-Stokes equations. This
corresponds to the hydrodynamic dc transport, where current whirlpools and
negative nonlocal resistance Gurzhi (1968); Levitov and Falkovich (2016);
Torre  _et al._ (2015); Pellegrino  _et al._ (2016); Chandra  _et al._ (2019a)
can be observed. This hydrodynamic regime in electronic systems has attracted
significant interest and has been experimentally demonstrated in graphene
Bandurin  _et al._ (2016); Kumar  _et al._ (2017); Bandurin  _et al._ (2018);
Berdyugin  _et al._ (2019); Ku  _et al._ (2020), (Al,Ga)As heterostructures
Braem  _et al._ (2018), GaAs quantum wells Gusev  _et al._ (2018, 2021) and
Weyl semimetals Gooth  _et al._ (2017, 2018); Osterhoudt  _et al._ (2021).

The parameter space for diffusive, ballistic, and hydrodynamic regimes changes
at finite frequency of external field $\omega$. A number of interesting
effects in ac transport was predicted theoretically for ideal electron fluid
Moessner  _et al._ (2018); Semenyakin and Falkovich (2018); Cohen and
Goldstein (2018); Alekseev (2018) where e-e collisions are so intense that
$\omega\tau_{ee}\ll 1$. With increasing the frequency, the electrons’ motion
becomes effectively ballistic if no collisions appear during the field cycle,
i.e. at $\omega\tau_{ee}\gg 1$ and $\omega\tau_{mr}\gg 1$. In many
experimentally relevant situations, particularly, at THz or GHz frequencies,
$\omega$ and $\tau_{ee}^{-1}$ are comparable in magnitude. This leads to the
need to construct models of the transport intermediate between hydrodynamic
and ballistic regimes Müller  _et al._ (2009, 2008); Narozhny  _et al._
(2015); Foster and Aleiner (2009); Svintsov  _et al._ (2012); Lucas  _et al._
(2016); Mayorov  _et al._ (2011b); Svintsov (2018); Chandra  _et al._ (2019b);
Holder  _et al._ (2019); Gupta  _et al._ (2021); Ledwith  _et al._ (2019);
Heremans  _et al._ (1992b).

It’s important to mention that ideal electron fluid is formed at very strong
e-e collisions, contrary to ideal electron gas where electrons do not interact
at all. Indeed, the dc kinematic viscosity $\nu_{kin}=v_{0}^{2}\tau_{ee}/4$
($v_{0}$ is the Fermi velocity) tends to zero for very short free path time
$\tau_{ee}$. On the other hand, the transport with long $\tau_{ee}$ can be
equivalently termed as ’highly viscous’ or ’ballistic’. The true Navier-Stokes
hydrodynamics is inapplicable in this regime. Yet it is possible to expand the
transport equations in the long-wavelength limit, generating the nearly local
’collisionless hydrodynamics’ or ’hydrodynamics of highly viscous fluid’.

Such equations predict the existence of a novel collective mode, the
transverce zero sound, in the absence of a magnetic field Lucas and Sarma
(2018); Khoo and Villadiego (2019). In finite magnetic field $B$, the dynamic
viscosity $\nu(\omega,B)$ was predicted to have the viscoelastic resonance
Alekseev (2018); Alekseev and Alekseeva (2019); Alekseev (2019). The latter
should occur at double cyclotron frequency $\omega=2\omega_{c}$. Interaction
between 2d plasmons and viscoelastic resonance was predicted to lead to
emergence of novel collective mode, dubbed as ”transverse zero magnetosound”
Alekseev and Alekseeva (2019).

The predicted mode structure of highly non-ideal electron fluid contradicts to
that of ballistic electrons, though both regimes should be equivalent and
appear at $\omega\tau_{ee}\gg 1$. The ballistic theory states that plasmon
dispersion interacts with multiple cyclotron resonances forming the so called
Bernstein modes Chiu and Quinn (1974); Bernstein (1958); Sitenko and Stepanov
(1957); Batke  _et al._ (1985); Gudmundsson  _et al._ (1995); Batke  _et al._
(1986); Holland  _et al._ (2004); Roldan  _et al._ (2011); Volkov and
Zabolotnykh (2014); Bandurin  _et al._ (2022). The fundamental Bernstein mode,
like conventional magnetoplasmon, starts from $\omega_{c}$ at $q=0$. As it
approaches $2\omega_{c}$, it reaches a plateau, and then falls downward with
negative group velocity. The group velocity on the plateau is zero, and
density of states is singular. The plateau frequency is shifted from
$2\omega_{c}$ by minigap. The next branch of the Bershtein modes starts from
$2\omega_{c}$ at $q=0$ and forms a plateau near the triple cyclotron
frequency. Similar branches occur at each multiple of the cyclotron frequency.

In this article, we are aiming to resolve the abovementioned contradictions
between ’ballistic’ and ’highly-viscous hydrodynamic’ models for 2DES in
magnetic field. To this end, we develop a classical kinetic model that
describes magnetotransport in ballistic, hydrodynamic, and intermediate
transport regimes (with $\varepsilon_{F}\gg T$). The model is based on
Boltzmann kinetic equation with model momentum and particle conserving e-e
collision integral. We calculate the nonlocal high-frequency conductivity with
non-zero magnetic field, which is a necessary block for calculating many
light-matter interaction characteristics Ciracì  _et al._ (2012); Gonçalves
_et al._ (2021); Lundeberg  _et al._ (2017a) (Section II).

Using the expression for conductivity, we study the evolution of the
dispersion of magnetoplasmons during the transition from the hydrodynamic
regime to the ballistic regime (Section III). At a high frequency of electron-
electron collisions corresponding to the hydrodynamic regime, an conventional
magnetoplasmon with a single cyclotron gap is observed. As the collision
frequency decreases, more and more pronounced dispersion splittings at
multiple cyclotron frequencies appear, which in the ballistic regime take on
the explicit form of Berstein modes. We inspect the approximation of
”collisionless” hydrodynamics (Section IV) and find that is applicable only
for wave vectors smaller than the anticrossing point $q^{*}$ of the two lowest
Bernstein modes, coinciding with them in the region of applicability. This
means that transverse zero magnetosound is an artifact of long-wave expansion
of transport equations beyond the region of applicability.

Figure 1: Real part of longitudinal conductivity $\sigma^{\prime}_{xx}$ in the
hydrodynamic regime with $\tau_{ee}=0.1$ ps (red line), intermediate regime
with $\tau_{ee}=1$ ps (green line), and in the ballistic regime with
$\tau_{ee}=10$ ps (blue line). The conductivity is normalized by collisionless
Drude conductivity $\sigma_{D}=ine^{2}/m\omega$. All curves are obtained at
carrier density $n_{0}=2\times 10^{12}$ cm-2, $f_{c}=1$ THz, and wave vector
$q=1$ $\mu$m-1, $qR_{c}\approx 6.3$.

##  II  Non-local high frequency conductivity tensor

We start with evaluation of the non-local conductivity tensor of 2DES
$\hat{\sigma}({\bf q},\omega)$ in a classically strong perpendicular magnetic
field $B\parallel\mathbf{e}_{z}$. To this end, we find the current response in
a weak electric field $\delta{\bf
E}=\delta\mathbf{E}_{0}e^{i(\mathbf{qr}-\omega t)}$ with
$\mathbf{r},\mathbf{q}\perp\mathbf{e}_{z}$. It is found from the solution of
linearized kinetic equation

$-(i\omega-i\mathbf{qv_{p}})\delta f-e\delta\mathbf{E_{0}v_{p}}\dfrac{\partial
f_{0}}{\partial\mathbf{\epsilon}}+\dfrac{1}{c}\mathbf{v_{p}}\times\mathbf{B}\dfrac{\partial\delta
f}{\partial\mathbf{p}}={\rm St}_{ee}\\{\delta f\\},$ (1)

where distribution function $f=f_{0}+\delta fe^{i(\mathbf{qr}-\omega t)}$,
${\bf v_{p}}$ is the electron velocity, $c$ the speed of light and ${\rm
St}_{ee}\\{\delta f\\}$ the electron-electron collision integral. In a
practically important limit $\varepsilon_{F}\gg T$, only electrons with
velocities close to the Fermi velocity $v_{F}$ participate in high-frequency
kinetic processes, which makes this model applicable for describing both
classical 2DES and graphene. Contrary to electron-impurity and electron-phonon
collisions, e-e collisions conserve the net momentum of colliding particles.
To account for this fact, we adopt the e-e collision integral in generalized
relaxation-time approximation Svintsov (2018); Abrikosov and Khalatnikov
(1959); Bhatnagar  _et al._ (1954); Gross and Krook (1956); Mermin (1970);
Lucas and Hartnoll (2018); Guo  _et al._ (2017); Lucas (2017). The collision
integral pushes the distribution function $\delta f$ towards local equilibrium

${\rm St}_{ee}\\{\delta f\\}=\frac{\delta f-\delta f_{hd}}{\tau_{ee}},\ \ \
\delta f_{hd}=-\frac{\partial
f_{0}}{\partial\varepsilon}\left(\delta\mu+\mathbf{p}\delta\mathbf{v}\right),$
(2)

where $\delta\mu$ and $\delta\mathbf{v}$ are found from the conservation of
particle number and momentum:

$\displaystyle\sum_{\mathbf{p}}\left(\delta f-\delta
f_{hd}\right)=0,\qquad\sum_{\mathbf{p}}\mathbf{p}\left(\delta f-\delta
f_{hd}\right)=0.$ (3)

The kinetic equation with e-e collision integral (2) and conservation laws (3)
are sufficient to describe the behaviour of electron liquid at large wave
vectors $q\sim\omega/v_{0}$, and across the whole hydrodynamic-to-ballistic
crossover. The scheme of solution is described in Appendix A. First, one
passes to the polar coordinates in momentum space and expands the distribution
function in Fourier series. Physically, the coefficients of these series are
cyclotron harmonics. The kinetic equation in cyclotron harmonic representation
is purely algebraic. Its symbolic solution at the first stage contains two
unknown variables, $\delta\mu$ and $\delta{\bf v}$. On the second stage, these
’hydrodynamic quantities’ are found from the conservation laws (3). In a
practically important limit $\varepsilon_{F}\gg T$, this leads us to a system
of linear equations describing both ballistic and hydrodynamic behavior of
2DES:

$\left(\begin{array}[]{ccc}1-i\gamma Y^{(0)}_{00}&-2i\gamma
Y^{(1)}_{00}&-2\gamma Y^{(0)}_{10}\\\ -i\gamma Y^{(1)}_{00}&1-2i\gamma
Y^{(2)}_{00}&-2\gamma Y^{(1)}_{10}\\\ \gamma Y^{(0)}_{10}&2\gamma
Y^{(1)}_{10}&1-2i\gamma
Y^{(0)}_{11}\end{array}\right)\left(\begin{array}[]{ccc}\frac{\delta
n}{n_{0}}\\\ \frac{\delta v_{x}}{v_{0}}\\\ \frac{\delta
v_{y}}{v_{0}}\end{array}\right)=-\left(\begin{array}[]{ccc}0\\\
iY^{(2)}_{00}E_{x}+Y^{(1)}_{10}E_{y}\\\
iY^{(0)}_{11}E_{y}-Y^{(1)}_{10}E_{x}\end{array}\right)\dfrac{2e}{\omega m},$
(4)

where $\gamma=(\omega\tau_{ee})^{-1}$ is the dimensionless collision
frequency, $v_{0}$ the Fermi velocity, $n_{0}$ equilibrium carrier density and
$\mathbf{q}\parallel\mathbf{e}_{x}$. The dimensionless factors

$Y_{ij}^{(k)}=\sum_{s=-\infty}^{\infty}\left(\dfrac{s}{qR_{c}}\right)^{k}\dfrac{J^{(i)}_{s}[qR_{c}]J^{(j)}_{s}[qR_{c}]}{1-s\omega_{c}/\omega+i\gamma},$
(5)

where $J^{(i)}_{s}[x]$ is the $i$-th order derivative of the Bessel function,
$R_{c}=v_{0}/\omega_{c}$ a cyclotron radius and $\omega_{c}$ a cyclotron
frequency.

These equations can be interpreted as generalized hydrodynamics; the term
’generalized’ implies applicability at arbitrarily large wave vectors. When
expanded to terms $\propto q^{2}$ (with $qR_{c}\ll 1$), we obtain
”collisionless” hydrodynamics ($\omega\tau_{ee}\gg 1$) or ”true” hydrodynamics
($\omega\tau_{ee}\ll 1$):

$\displaystyle-i\omega\delta n+i(\mathbf{q}\delta\mathbf{v})n_{0}=0,$ (6)
$\displaystyle-i\omega\delta\mathbf{v}=\omega_{c}[\delta\mathbf{v}\times\mathbf{e_{z}}]-i\mathbf{q}\dfrac{v_{0}^{2}}{2}\delta
n-\nu_{xx}(\omega)q^{2}\delta\mathbf{v}-$
$\displaystyle\qquad\qquad\qquad\qquad-\nu_{xy}(\omega)q^{2}\delta[\mathbf{v}\times\mathbf{e_{z}}]+\dfrac{e\delta\mathbf{E}_{0}}{m},$
(7)

where the ac shear viscosity coefficients $\nu_{xx}(\omega)$ and
$\nu_{xy}(\omega)$ depend on magnetic field and frequency as

$\displaystyle\nu_{xx}(\omega)=\nu_{0}\dfrac{1-i\omega\tau_{ee}}{(1-i\omega\tau_{ee})^{2}+4\omega_{c}^{2}},$
(8)
$\displaystyle\nu_{xy}(\omega)=\nu_{0}\dfrac{2\omega_{c}\tau_{ee}}{(1-i\omega\tau_{ee})^{2}+4\omega_{c}^{2}},$
(9)

where $\nu_{0}=v_{0}^{2}\tau_{ee}/4$ is the viscosity in the absence of
magnetic field. Whereas ”true” hydrodynamics, due to the small relaxation
time, is also valid for large wave vectors, collisionless hydrodynamics has a
limited range of applicability ($qv_{0}\ll\omega_{c}$) which requires special
care when applying this approximation for large wave vectors.

Finally, the components of the conductivity tensor $\hat{\sigma}(\omega,q)$
can be found from

${\bf j}=en_{0}{\delta\bf v}=\hat{\sigma}({\bf q},\omega)\delta{\bf E}_{0},$
(10)

where the relation between ${\delta\bf v}$ and $\delta{\bf E}_{0}$ is found
from equations (4) .

It is instructive to track the changes in longitudinal conductivity
$\sigma_{xx}({\bf q},\omega)$ while varying the strength of e-e collisions.
Fig.1 illustrates this evolution and shows the frequency dependence of ${\rm
Re}\sigma_{xx}$ at fixed $q=1\ \mu m^{-1}$ and $\omega_{c}/2\pi=1$ THz. In the
ballistic regime ($\tau_{ee}=10$ ps, blue curve), the conductivity has
resonances at multiple cyclotron frequencies which amplitude decreases with
increasing the resonance order. As the frequency of e-e collisions
$\tau_{ee}^{-1}$ becomes comparable to $\omega_{c}$, multiple resonances
rapidly fade away. Only the main cyclotron resonance persists in the
hydrodynamic regime, and its width shrinks again at very short values of
$\tau_{ee}$, i.e. at low viscosity.

##  III Magnetoplasmon dispersion

Figure 2: Magnetoplasmon dispersion curve visualized through loss function
${\rm Im}\varepsilon^{-1}(q,\omega,\omega_{c})$ in (a) the hydrodynamic
transport regime in graphene, $\tau_{ee}=0.1$ ps, (b) intermediate regime,
$\tau_{ee}=1$ ps, (c) in ballistic regime, $\tau_{ee}=10$ ps. All curves are
obtained at carrier density $n_{0}=2\times 10^{12}$ cm-2, $f_{c}=1$ THz with
gate distance $d=100$ nm, $\chi$ = 4

We proceed to inspect the evolution of collective excitations in magnetized
2DES across the hydrodynamic-ballistic crossover. To address this problem, we
consider the excitation of collective modes by a point dipole located above
the 2DES in the quasi-static regime ($\omega/q\ll c$). We find that total
electric potential in the 2DES plane can be represented as (Appendix B):

$\varphi({\bf q},\omega)=\frac{\varphi_{0}({\bf q},\omega)}{\varepsilon({\bf
q},\omega)},$ (11)

where $\varphi_{0}({\bf q},\omega)$ is the electric potential of the point
source in the absence of 2DES and $\varepsilon({\bf q},\omega)$ the 2DES
dielectric function. In the quasi-static regime, it is affected only by the
longitudinal component of the conductivity tensor $\sigma_{xx}({\bf
q},\omega)$ and is given by

$\varepsilon(\omega,q)=1+i\dfrac{2\pi
q}{\chi\omega}(1-e^{-2qd})\sigma_{xx}({\bf q},\omega),$ (12)

where $\chi$ is the effective dielectric constant, $d$ the distance to the
gate electrode located below the 2DES and the longitudinal conductivity
$\sigma_{xx}$ was obtained before in Sec. II.

It is convenient to visualize the dispersion of plasmons using the loss
function ${\rm Im}\varepsilon^{-1}$. The peaks in loss function correspond to
a resonant enhancement of external electric field, i.e. to the collective
electrostatic waves – (magneto)plasmons. Using the obtained generalized
equations of hydrodynamics, it is possible to trace the evolution of
collective modes with an increase in the frequency of electron-electron
collisions up to the transition to the hydrodynamic regime.

The emerging loss functions are displayed in Fig. 2. Our calculations show
that the splittings of plasmon dispersion at multiple cyclotron frequencies
are gradually blurred with increasing e-e collision frequency
$\tau_{ee}^{-1}$. In the hydrodynamic regime, the dispersion collapses into an
ordinary magnetoplasmon

$\omega^{2}=\omega_{c}^{2}+s(q)^{2}q^{2},$ (13)

where $s(q)=\sqrt{\dfrac{2\pi n_{0}e^{2}}{m\kappa q}(1-e^{-2qd})}$.

Figure 3: Magnetoplasmon dispersion curve in graphene visualized through loss
function ${\rm Im}\varepsilon^{-1}(q,\omega,\omega_{c})$ with $\tau_{ee}=20$
ps, carrier density $n_{0}=6.7\times 10^{9}$ cm-2, $f_{c}=0.05$ THz with gate
distance $d=50$ nm, $\chi$ = 4. Black dashed line shows the boundary of the
domain of Landau damping, yellow dashed line shows the conventional
magnetoplasmon dispersion curve.

Now it’s important to track the transition from the obtained results to the
zero magnetic field. In the absence of a magnetic field in the ballistic
regime, the phase velocity of plasmons does not fall below the Fermi velocity
in Dirac materials; thus, the dispersion of plasmons does not fall into the
Landau damping region, and at low values of the carrier density $n_{0}$ and a
low gate distance $d$, the dispersion asymptotically tends to the boundary
$\omega=qv_{0}$ Svintsov (2018); Lundeberg  _et al._ (2017b). However, in the
presence of a weak magnetic field, this boundary is destroyed due to the
splitting of the Bershtein modes (Fig. 3)

Figure 4: Magnetoplasmon dispersion curve in graphene visualized through loss
function ${\rm Im}\varepsilon^{-1}(q,\omega,\omega_{c})$ in (a) the
”collisionless” hydrodynamics model, dashed line shows the asymptotics of the
lowest mode corresponding to the transverse magnetosound mode
$\omega^{2}=4\omega_{c}^{2}+u_{0}^{2}q^{2}/4$; (b) in the full theory, dashed
lines correspond to modes obtained in the ”collisionless” hydrodynamics
approximation. All curves are obtained at carrier density $n_{0}=2\times
10^{12}$ cm-2, $f_{c}=1$ THz with gate distance $d=100$ nm, $\chi$ = 4 and
$\tau_{ee}=10$ ps

##  IV Collisionless hydrodynamic model

In this section, we study the range of applicability of collisionless
hydrodynamics for the description of magnetoplasmon modes. Using the
continuity equation (6) and the Navier-Stokes equations (7), we find the
conductivity of a two-dimensional system by definition (10). We plot the loss
function in the collisionless hydrodynamic model (6-7) and in the ballistic
transport model (4) (Fig.4). In the first one (Fig.4a) the zeros of the
dielectric function (12) in this case give two magnetoplasmon modes. The first
of them goes from the cyclotron frequency $\omega_{c}$, at small wave vectors
it coincides with the conventional magnetoplasmon (13), and at large wave
vectors (after the anti-crossing point $q^{*}$) it coincides with the
transverse magnetosound

$\omega^{2}=4\omega_{c}^{2}+\dfrac{v_{0}^{2}q^{2}}{4},$ (14)

which is related to perturbations of shear stress of a charged Fermi liquid in
a magnetic field. Alekseev and Alekseeva (2019) The second mode starts from
the doubled cyclotron frequency and after the anti-crossing point tends to
conventional magnetoplasmon. It is possible to approximately estimate the
anti-crossing point $(\omega^{*},q^{*})$ as the point at which the value of
doubled cyclotron frequency $2\omega_{c}$ is reached by the dispersion of
conventional magnetoplasmon (13):

$\omega^{*}\approx 2\omega_{c},\ \ \ \ \ \ \
q^{*}\approx\dfrac{\sqrt{3}\omega_{c}}{s}.$ (15)

In this case, up to the anti-crossing point $q^{*}$, there is a good agreement
with the two lowest Bernstein modes observed in the full theory (Fig.4b) (but
modes starting at multiple cyclotron frequencies are not observed). However,
for large wave vectors, serious discrepancies are observed. In the
collisionless hydrodynamics model, the lowest mode exceeds twice the cyclotron
frequency and coincides with the magnetosonic mode, while in the full theory,
the fundamental Bernstein mode never exceeds twice the cyclotron frequency and
forms a plateau with a singular density of states.

##  V Discussion and conclusion

In the limit $\omega\tau_{ee}\gg 1$, the constructed full kinetic theory
describes Bernstein modes. Magnetosound waves predicted in the framework of
collisionless hydrodynamics described by the equations (7) in Ref. Alekseev
and Alekseeva (2019) coincide with the full theory only in a narrow $q$-range.
Since the lowest mode in the collisionless hydrodynamics model tends to
magnetosound mode (14) only after the anti-crossing point (15), where this
model is inapplicable, it can be argued that transverse magnetosound modes are
an artifact of spectrum cutoff in generalized hydrodynamic equations (4) and
are not observed in the full theory.

While the physics behind these differences is clear, it’s interesting to see
and compare the emerging absorption spectra. The dielectric function
$\varepsilon(q,\omega,\omega_{c})$ corresponds to the screening by electrons
of the two-dimensional system of the external field (Appendix B). The loss
function Im$[\varepsilon^{-1}]$ is responsible for the magnetoplasmon-assisted
absorption. It peaks at the magnetoplasmon modes, and therefore should enhance
the near-field magnetoabsorption of an inhomogeneous field incident on the
2DES. Bandurin  _et al._ (2022)

Both transverse magnetosound and Bernstein modes are predicted in the loss
function (Fig.4). In both cases, for pure materials, an asymmetric (with
respect to the magnetic field) resonance enhancement of magnetoabsorption in
the vicinity of the doubled cyclotron frequency is predicted. However, in the
case of magnetosound, the sharp gradow side of the resonance is located in the
direction of high magnetic fields, and in the case of Bershtein modes, on the
contrary, in the direction of low magnetic fields, which makes it possible to
distinguish one from the other in measurements.

Measurements of the photoresistence demonstrate similar asymmetric resonances
at double, triple, and even quadruple the cyclotron frequency in pure graphene
in a magnetic field Bandurin  _et al._ (2022) that matches the Bernstein modes
seen in the full theory and are not predicted by the model of collisionless
hydrodynamics. In the full theory the sharp gradow side of the resonance is
located in the direction of low magnetic fields. In Ref. Alekseev and
Alekseeva (2019) it was shown for a thin strip of a 2DES in the case of
collisionless hydrodynamics that the sharp gradow side of the resonance is
located in the direction of high magnetic fields. However, the shape of the
resonance strongly depends on the width of the strip and becomes almost
symmetric in the case of an infinite 2DES. It can be assumed that in this case
the shape of the resonance should strongly depend on the shape of the sample.
However, measurements of various samples show an asymmetry consistent with the
full theory Dai  _et al._ (2010); Bandurin  _et al._ (2022), which is another
argument that these resonances are caused by Bernstein modes.

In summary, a model of high-frequency dynamics of a two-dimensional electronic
system has been constructed, which allows one to study the transition from the
ballistic regime to the hydrodynamic. It is shown that the conventional
magnetoplasmon in the hydrodynamic regime, with a decrease in the frequency of
electron-electron collisions, acquires characteristic splittings at multiple
cyclotron frequencies and becomes Bernstein modes in the ballistic limit.
According to calculations, the model of ”collisionless” hydrodynamics, which
predicts transverse magnetosound, is applicable to the description of
magnetoplasmon modes only for wave vectors not exceeding the anti-crossing
point (15) of the two lowest Berstein modes, and magnetosound mode are beyond
the range of applicability of this model.

###### Acknowledgements.

The work was supported by grants # 16-19-10557 and # 21-12-00287 of the
Russian Science Foundation. The authors thank I.A. Dmitriev for helpful
discussions.

### Appendix A: Solution of the kinetic equation

To solve the kinetic equation (1), we pass to polar coordinates (with
$\mathbf{q}\parallel\mathbf{e}_{x}$) and represent the equation as:

$\displaystyle-i(\alpha+i\gamma_{c})\delta f+i\beta\cos{\varphi}\delta
f+\dfrac{\partial\delta f}{\partial\varphi}=\dfrac{\partial f_{0}}{\partial
p}\left(\dfrac{eE_{x}}{\omega_{c}}\cos{\varphi}+\dfrac{eE_{y}}{\omega_{c}}\sin{\varphi}+\dfrac{\gamma_{c}\delta\mu}{v_{0}}+\dfrac{\gamma_{c}p_{F}\delta
v_{x}}{v_{0}}\cos{\varphi}+\dfrac{\gamma_{c}p_{F}\delta
v_{y}}{v_{0}}\sin{\varphi}\right),$ (16)

where $\alpha=\omega/\omega_{c}$, $\beta=qR_{c}$,
$\gamma_{c}=(\omega_{c}\tau_{ee})^{-1}$, and $\omega_{c}$ is the cyclotron
frequency. Then, we seek for the function $\delta f$ in the form $\delta
f=ge^{-i\beta\sin{\varphi}}$, and expand $g$ as a series of angular harmonics
$\varphi$:

$g=\sum_{s=-\infty}^{+\infty}g_{s}e^{is\varphi}.$ (17)

The solution for $g_{s}$ reads as:

$\displaystyle g_{s}=iJ_{s}(\beta)\dfrac{\partial f_{0}}{\partial
p}\dfrac{\left(\frac{eE_{x}}{\omega_{c}}\frac{s}{\beta}+\frac{\gamma_{c}\delta\mu}{v_{0}}+\frac{\gamma_{c}p_{F}\delta
v_{x}}{v_{0}}\frac{s}{\beta}\right)}{\alpha+i\gamma_{c}-s}+$
$\displaystyle+J_{s}^{\prime}(\beta)\dfrac{\left(\frac{eE_{y}}{\omega_{c}}+\frac{\gamma_{c}p_{F}\delta
v_{y}}{v_{0}}\frac{s}{\beta}\right)}{\alpha+i\gamma_{c}-s},$ (18)

where
$J_{s}(\beta)=\dfrac{1}{2\pi}\int_{0}^{2\pi}e^{i(\beta\sin{\varphi}-s\varphi)}d\varphi$
is the $s$-th order Bessel function and $J^{\prime}_{s}(\beta)$ its
derivative. Using the equations of conservation of particles and momentum (3),
we obtain the equations of high-frequency dynamics (4)

### Appendix B: Magnetoabsorption calculation

In this section, we first consider the absorption of an inhomogeneous wave by
an infinite two-dimensional electron system, then consider the case when an
inhomogeneous field is created by planar charge density $\rho(\mathbf{r})\sim
e^{-i\omega t}$ located a short distance $L$ from the plane of a 2DES ($qL\ll
1$) (Fig.5). The absorption given by the Joule’s law:

$\displaystyle
P=2\int\operatorname{Re}{[\mathbf{E}(\mathbf{q},\omega)\mathbf{j}^{*}]}d\mathbf{q},$
(19)

where $\mathbf{E}(\mathbf{q},\omega)$ is the total field in the plane of the
two-dimensional system. Further, we take into account the Ohm’s law
$\mathbf{j}=\hat{\sigma}\mathbf{E}$ with conductivity tensor

$\hat{\sigma}=\begin{pmatrix}\sigma_{xx}&-\sigma_{xy}\\\
\sigma_{xy}&\sigma_{xx}\end{pmatrix}.$ (20)

and obtain an expression for the absorbed power in the form

$\displaystyle
P=2\int\dfrac{d\mathbf{q}}{(2\pi)^{2}}\left(\sigma^{\prime}_{xx}|\mathbf{E}(\mathbf{q},\omega)|^{2}+\sigma^{\prime\prime}_{xy}\operatorname{Im}{[E^{*}_{x}(\mathbf{q},\omega)E_{y}({\mathbf{q},\omega})]}\right).$
(21)

In the plane of a two-dimensional system, the field is screened by electrons
of 2DES. In the electrostatic approximation ($qc\gg\omega$), suitable for
describing a sharp increase in absorption by a strongly inhomogeneous
Bernstein mode field, the total potential in the plane of the system is given
by the expression

$\mathbf{\varphi}|_{z=0}=\dfrac{2\pi\rho(\mathbf{q})e^{-qL}}{\chi
q\varepsilon({\bf q},\omega)},$ (22)

where $\varepsilon({\bf q},\omega)$ is the dielectric function and $\chi$
substrate dielectric constant.

Figure 5: Schematic representation of a two-dimensional system, above which an
inhomogeneous charge density is located near

To find the dielectric function $\varepsilon({\bf q},\omega)$ of a planar
symmetric non-gated two-dimensional system with a known conductivity tensor
$\hat{\sigma}$, we solved the field equation (with $qL\ll 1$)

$\left(\dfrac{\partial^{2}}{\partial
z^{2}}-q^{2}\right)\varphi(\mathbf{q},z)=\dfrac{4\pi}{\chi}(\rho_{ind}(\mathbf{q})+\rho(\mathbf{q}))\delta(z),$
(23)

where $\omega\rho_{ind}(\mathbf{q})=\mathbf{q}\mathbf{j}$ and
$\mathbf{j}=\hat{\sigma}\mathbf{E}$. Solving the obtained equations, we obtain
an expression for the dielectric function

$\varepsilon({\bf q},\omega)=1+i\dfrac{2\pi
q}{\chi\omega}\sigma_{xx}(\omega,q).$ (24)

Taking into account the screening by a gate located at a distance $d$, we
obtain the expression (12)

In this case, the expression for absorption takes the form

$P_{\rm near}=2\int{d{\bf
q}\frac{\omega}{2\pi\chi^{2}q}|\rho(\mathbf{q})|^{2}\operatorname{Im}\frac{1}{\varepsilon({\bf
q},\omega)}}.$ (25)

The specific type of charge density $\rho(\mathbf{q})$ depends on the problem.
Absorption caused by an inhomogeneous field of charges resulting from the
diffraction of an incident plane wave on a thin contact was considered in
detail in Ref. Bandurin  _et al._ (2022)

## References

  * Mayorov  _et al._ (2011a) A. S. Mayorov, R. V. Gorbachev, S. V. Morozov, L. Britnell, R. Jalil, L. A. Ponomarenko, P. Blake, K. S. Novoselov, K. Watanabe, T. Taniguchi, _et al._ , Nano letters 11, 2396 (2011a).
  * Heremans  _et al._ (1992a) J. Heremans, M. Santos, and M. Shayegan, Applied physics letters 61, 1652 (1992a).
  * Gurzhi (1968) R. Gurzhi, Soviet Physics Uspekhi 11, 255 (1968).
  * Levitov and Falkovich (2016) L. Levitov and G. Falkovich, Nature Physics 12, 672 (2016).
  * Torre  _et al._ (2015) I. Torre, A. Tomadin, A. K. Geim, and M. Polini, Physical Review B 92, 165433 (2015).
  * Pellegrino  _et al._ (2016) F. M. Pellegrino, I. Torre, A. K. Geim, and M. Polini, Physical Review B 94, 155414 (2016).
  * Chandra  _et al._ (2019a) M. Chandra, G. Kataria, D. Sahdev, and R. Sundararaman, Physical Review B 99, 165409 (2019a).
  * Bandurin  _et al._ (2016) D. Bandurin, I. Torre, R. K. Kumar, M. B. Shalom, A. Tomadin, A. Principi, G. Auton, E. Khestanova, K. Novoselov, I. Grigorieva, _et al._ , Science 351, 1055 (2016).
  * Kumar  _et al._ (2017) R. K. Kumar, D. Bandurin, F. Pellegrino, Y. Cao, A. Principi, H. Guo, G. Auton, M. B. Shalom, L. A. Ponomarenko, G. Falkovich, _et al._ , Nature Physics 13, 1182 (2017).
  * Bandurin  _et al._ (2018) D. A. Bandurin, A. V. Shytov, L. S. Levitov, R. K. Kumar, A. I. Berdyugin, M. B. Shalom, I. V. Grigorieva, A. K. Geim, and G. Falkovich, Nature communications 9, 1 (2018).
  * Berdyugin  _et al._ (2019) A. I. Berdyugin, S. Xu, F. Pellegrino, R. K. Kumar, A. Principi, I. Torre, M. B. Shalom, T. Taniguchi, K. Watanabe, I. Grigorieva, _et al._ , Science 364, 162 (2019).
  * Ku  _et al._ (2020) M. J. Ku, T. X. Zhou, Q. Li, Y. J. Shin, J. K. Shi, C. Burch, L. E. Anderson, A. T. Pierce, Y. Xie, A. Hamo, _et al._ , Nature 583, 537 (2020).
  * Braem  _et al._ (2018) B. A. Braem, F. Pellegrino, A. Principi, M. Röösli, C. Gold, S. Hennel, J. V. Koski, M. Berl, W. Dietsche, W. Wegscheider, _et al._ , Physical Review B 98, 241304 (2018).
  * Gusev  _et al._ (2018) G. Gusev, A. Levin, E. Levinson, and A. Bakarov, Physical Review B 98, 161303 (2018).
  * Gusev  _et al._ (2021) G. Gusev, A. Jaroshevich, A. Levin, Z. Kvon, and A. Bakarov, Physical Review B 103, 075303 (2021).
  * Gooth  _et al._ (2017) J. Gooth, A. C. Niemann, T. Meng, A. G. Grushin, K. Landsteiner, B. Gotsmann, F. Menges, M. Schmidt, C. Shekhar, V. Süß, _et al._ , Nature 547, 324 (2017).
  * Gooth  _et al._ (2018) J. Gooth, F. Menges, N. Kumar, V. Sü$\beta$, C. Shekhar, Y. Sun, U. Drechsler, R. Zierold, C. Felser, and B. Gotsmann, Nature communications 9, 1 (2018).
  * Osterhoudt  _et al._ (2021) G. B. Osterhoudt, Y. Wang, C. A. Garcia, V. M. Plisson, J. Gooth, C. Felser, P. Narang, and K. S. Burch, Physical Review X 11, 011017 (2021).
  * Moessner  _et al._ (2018) R. Moessner, P. Surówka, and P. Witkowski, Physical Review B 97, 161112 (2018).
  * Semenyakin and Falkovich (2018) M. Semenyakin and G. Falkovich, Physical Review B 97, 085127 (2018).
  * Cohen and Goldstein (2018) R. Cohen and M. Goldstein, Physical Review B 98, 235103 (2018).
  * Alekseev (2018) P. Alekseev, Physical Review B 98, 165440 (2018).
  * Müller  _et al._ (2009) M. Müller, J. Schmalian, and L. Fritz, Physical review letters 103, 025301 (2009).
  * Müller  _et al._ (2008) M. Müller, L. Fritz, and S. Sachdev, Physical Review B 78, 115406 (2008).
  * Narozhny  _et al._ (2015) B. Narozhny, I. Gornyi, M. Titov, M. Schütt, and A. Mirlin, Physical Review B 91, 035414 (2015).
  * Foster and Aleiner (2009) M. S. Foster and I. L. Aleiner, Physical Review B 79, 085415 (2009).
  * Svintsov  _et al._ (2012) D. Svintsov, V. Vyurkov, S. Yurchenko, T. Otsuji, and V. Ryzhii, Journal of Applied Physics 111, 083715 (2012).
  * Lucas  _et al._ (2016) A. Lucas, R. A. Davison, and S. Sachdev, Proceedings of the National Academy of Sciences 113, 9463 (2016).
  * Mayorov  _et al._ (2011b) A. S. Mayorov, R. V. Gorbachev, S. V. Morozov, L. Britnell, R. Jalil, L. A. Ponomarenko, P. Blake, K. S. Novoselov, K. Watanabe, T. Taniguchi, _et al._ , Nano letters 11, 2396 (2011b).
  * Svintsov (2018) D. Svintsov, Physical Review B 97, 121405 (2018).
  * Chandra  _et al._ (2019b) M. Chandra, G. Kataria, D. Sahdev, and R. Sundararaman, Physical Review B 99, 165409 (2019b).
  * Holder  _et al._ (2019) T. Holder, R. Queiroz, T. Scaffidi, N. Silberstein, A. Rozen, J. A. Sulpizio, L. Ella, S. Ilani, and A. Stern, Physical Review B 100, 245305 (2019).
  * Gupta  _et al._ (2021) A. Gupta, J. Heremans, G. Kataria, M. Chandra, S. Fallahi, G. Gardner, and M. Manfra, Physical Review Letters 126, 076803 (2021).
  * Ledwith  _et al._ (2019) P. Ledwith, H. Guo, A. Shytov, and L. Levitov, Phys. Rev. Lett. 123, 116601 (2019).
  * Heremans  _et al._ (1992b) J. Heremans, M. Santos, and M. Shayegan, Applied physics letters 61, 1652 (1992b).
  * Lucas and Sarma (2018) A. Lucas and S. D. Sarma, Physical Review B 97, 115449 (2018).
  * Khoo and Villadiego (2019) J. Y. Khoo and I. S. Villadiego, Physical Review B 99, 075434 (2019).
  * Alekseev and Alekseeva (2019) P. Alekseev and A. Alekseeva, Physical review letters 123, 236801 (2019).
  * Alekseev (2019) P. Alekseev, Semiconductors 53, 1367 (2019).
  * Chiu and Quinn (1974) K. Chiu and J. Quinn, Physical Review B 9, 4724 (1974).
  * Bernstein (1958) I. B. Bernstein, Physical Review 109, 10 (1958).
  * Sitenko and Stepanov (1957) A. Sitenko and K. Stepanov, Soviet Phys. JETP 4 (1957).
  * Batke  _et al._ (1985) E. Batke, D. Heitmann, J. Kotthaus, and K. Ploog, Physical review letters 54, 2367 (1985).
  * Gudmundsson  _et al._ (1995) V. Gudmundsson, A. Brataas, P. Grambow, B. Meurer, T. Kurth, and D. Heitmann, Physical Review B 51, 17744 (1995).
  * Batke  _et al._ (1986) E. Batke, D. Heitmann, and C. Tu, Physical Review B 34, 6951 (1986).
  * Holland  _et al._ (2004) S. Holland, C. Heyn, D. Heitmann, E. Batke, R. Hey, K. Friedland, and C.-M. Hu, Physical review letters 93, 186804 (2004).
  * Roldan  _et al._ (2011) R. Roldan, M. Goerbig, and J.-N. Fuchs, Physical Review B 83, 205406 (2011).
  * Volkov and Zabolotnykh (2014) V. Volkov and A. Zabolotnykh, Physical Review B 89, 121410 (2014).
  * Bandurin  _et al._ (2022) D. A. Bandurin, E. Mönch, K. Kapralov, I. Y. Phinney, K. Lindner, S. Liu, J. H. Edgar, I. A. Dmitriev, P. Jarillo-Herrero, D. Svintsov, and S. D. Ganichev, Nature Physics (2022), .
  * Ciracì  _et al._ (2012) C. Ciracì, R. Hill, J. Mock, Y. Urzhumov, A. Fernández-Domínguez, S. Maier, J. Pendry, A. Chilkoti, and D. Smith, Science 337, 1072 (2012).
  * Gonçalves  _et al._ (2021) P. Gonçalves, T. Christensen, N. M. Peres, A.-P. Jauho, I. Epstein, F. H. Koppens, M. Soljačić, and N. A. Mortensen, Nature communications 12, 1 (2021).
  * Lundeberg  _et al._ (2017a) M. B. Lundeberg, Y. Gao, R. Asgari, C. Tan, B. Van Duppen, M. Autore, P. Alonso-González, A. Woessner, K. Watanabe, T. Taniguchi, _et al._ , Science 357, 187 (2017a).
  * Abrikosov and Khalatnikov (1959) A. Abrikosov and I. Khalatnikov, Reports on Progress in Physics 22, 329 (1959).
  * Bhatnagar  _et al._ (1954) P. L. Bhatnagar, E. P. Gross, and M. Krook, Physical review 94, 511 (1954).
  * Gross and Krook (1956) E. P. Gross and M. Krook, Physical Review 102, 593 (1956).
  * Mermin (1970) N. D. Mermin, Physical Review B 1, 2362 (1970).
  * Lucas and Hartnoll (2018) A. Lucas and S. A. Hartnoll, Physical Review B 97, 045105 (2018).
  * Guo  _et al._ (2017) H. Guo, E. Ilseven, G. Falkovich, and L. S. Levitov, Proceedings of the National Academy of Sciences 114, 3068 (2017).
  * Lucas (2017) A. Lucas, Physical Review B 95, 115425 (2017).
  * Lundeberg  _et al._ (2017b) M. B. Lundeberg, Y. Gao, R. Asgari, C. Tan, B. Van Duppen, M. Autore, P. Alonso-González, A. Woessner, K. Watanabe, T. Taniguchi, _et al._ , Science 357, 187 (2017b).
  * Dai  _et al._ (2010) Y. Dai, R. Du, L. Pfeiffer, and K. West, Physical review letters 105, 246802 (2010).