# Group-wise Correlation Stereo Network

Xiaoyang Guo
Institution1
Institution1 address
<EMAIL_ADDRESS>Second Author
Institution2
First line of institution2 address
<EMAIL_ADDRESS>Xiaoyang Guo1 Kai Yang2 Wukui Yang2 Xiaogang Wang1
Hongsheng Li1
1 The Chinese University of Hong Kong 2 SenseTime Research
{xyguo, xgwang<EMAIL_ADDRESS>{yangkai<EMAIL_ADDRESS>
###### Abstract

Stereo matching estimates the disparity between a rectified image pair, which
is of great importance to depth sensing, autonomous driving, and other related
tasks. Previous works built cost volumes with cross-correlation or
concatenation of left and right features across all disparity levels, and then
a 2D or 3D convolutional neural network is utilized to regress the disparity
maps. In this paper, we propose to construct the cost volume by group-wise
correlation. The left features and the right features are divided into groups
along the channel dimension, and correlation maps are computed among each
group to obtain multiple matching cost proposals, which are then packed into a
cost volume. Group-wise correlation provides efficient representations for
measuring feature similarities and will not lose too much information like
full correlation. It also preserves better performance when reducing
parameters compared with previous methods. The 3D stacked hourglass network
proposed in previous works is improved to boost the performance and decrease
the inference computational cost. Experiment results show that our method
outperforms previous methods on Scene Flow, KITTI 2012, and KITTI 2015
datasets. The code is available at https://github.com/xy-guo/GwcNet

##  1 Introduction

Accurate depth sensing is the core of many computer vision applications, such
as autonomous driving, robot navigation, and shallow depth-of-field image
synthesis. Stereo matching belongs to passive depth sensing techniques, which
estimates depth by matching pixels from rectified image pairs captured by two
cameras. The disparity $d$ of a pixel can be converted into depth by $Fl/d$,
where $F$ denotes the focal length of the camera lens, and $l$ is the distance
between two camera centers. Therefore, the depth precision improves with the
precision of disparity prediction.

Traditional stereo pipelines usually consist of all or portion of the
following four steps, matching cost computation, cost aggregation, disparity
optimization, and post-processing [18]. Matching cost computation provides
initial similarity measures for left image patches and possible corresponding
right image patches, which is a crucial step of stereo matching. Some common
matching costs include absolute difference (SAD), sum of squared difference
(SSD), and normalized cross-correlation (NCC). The cost aggregation and
optimization steps incorporate contextual matching costs and regularization
terms to obtain more robust disparity predictions.

Learning-based methods explore different feature representations and
aggregation algorithms for matching costs. DispNetC [14] computes a
correlation volume from the left and right image features and utilizes a CNN
to directly regress disparity maps. GC-Net [6] and PSMNet [2] construct
concatenation-based feature volume and incorporate a 3D CNN to aggregate
contextual features. There are also works [1, 19] trying to aggregate evidence
from multiple hand-crafted matching cost proposals. However, the above methods
have several drawbacks. The full correlation [14] provides an efficient way
for measuring feature similarities, but it loses much information because it
produces only a single-channel correlation map for each disparity level. The
concatenation volume [6, 2] requires more parameters in the following
aggregation network to learn the similarity measurement function from scratch.
[1, 19] stills utilizes traditional matching costs and cannot be optimized
end-to-end.

In this paper, we propose a simple yet efficient operation called group-wise
correlation to tackle the above drawbacks. Multi-level unary features are
extracted and concatenated to form high-dimensional feature representations
$\mathbf{f}_{l},\mathbf{f}_{r}$ for a left-right image pair. Then, the
features are split into multiple groups along the channel dimension, and the
$i$th left feature group is cross-correlated with the corresponding $i$th
right feature group over all disparity levels to obtain group-wise correlation
maps. At last, all the correlation maps are packed to form the final 4D cost
volume. The unary features can be treated as groups of structured vectors
[25], so the correlation maps for a certain group can be seen as a matching
cost proposal. In this way, we can leverage the power of traditional cross-
correlation matching cost and provide better similarity measures for the
following 3D aggregation network compared with [6, 2]. The multiple matching
cost proposals also avoid the information loss like full correlation [14].

The 3D stacked hourglass aggregation network proposed in PSMNet [2] is
modified to further improve the performance and decrease the inference
computational cost. $1{\times}1{\times}1$ 3D convolutions are employed in the
shortcut connections within each hourglass module without increasing too much
computational cost.

Our main contributions can be summarized as follows. 1) We propose group-wise
correlation to construct cost volumes to provide better similarity measures.
2) The stacked 3D hourglass refinement network is modified to improve the
performance without increasing the inference time. 3) Our method achieves
better performance than previous methods on Scene Flow, KITTI 2012, and KITTI
2015 datasets. 4) Experiment results show that when limiting the computational
cost of the 3D aggregation network, the performance reduction of our proposed
network is much smaller than previous PSMNet, which makes group-wise
correlation a valuable way to be implemented in real-time stereo networks.

##  2 Related Work

###  2.1 Traditional methods

Generally, traditional stereo matching consists of all or portion of the
following four steps: matching cost computation, cost aggregation, disparity
optimization, and some post-processing steps [18]. In the first step, the
matching costs of all pixels are computed for all possible disparities. Common
matching costs include sum of absolute difference (SAD), sum of squared
difference (SSD), normalized cross-correlation (NCC), and so on. Local methods
[30, 27, 15] explore different strategies to aggregate matching costs with
neighbor pixels and usually utilize the winner-take-all (WTA) strategy to
choose the disparity with minimum matching cost. In contrast, global methods
minimize a target function to solve the optimal disparity map, which usually
takes both matching costs and smoothness priors into consideration, such as
belief propagation [23, 9] and graph cut [11]. Semi-global matching (SGM) [4]
approximates the global optimization with dynamic programming. Local and
global methods can be combined to obtain better performance and robustness.

###  2.2 Learning based methods

Besides hand-crafted methods, researchers also proposed many learned matching
costs [29, 13, 21] and cost aggregation algorithms [1, 19]. Zbontar and Lecun
[29] first proposed to compute matching costs using neural networks. The
predicted matching costs are then processed with traditional cross-based cost
aggregation and semi-global matching to predict the disparity map. The
matching cost computation was accelerated in [13] by correlating unary
features. Batsos _et al_. proposed CBMV [1] to combine evidence from multiple
basic matching costs. Schonberger _et al_. [19] proposed to classify scanline
matching cost candidates with a random forest classifier. Seki _et al_.
proposed SGM-Nets [20] to provide learned penalties for SGM. Knobelreiter _et
al_. [10] proposed to combine CNN-predicted correlation matching costs and CRF
to integrate long-range interactions.

Following DispNetC (Mayer _et al_. [14]), there are a lot of works directly
regressing disparity maps from correlation cost volumes [17, 12, 22, 26].
Given the left and the right feature maps $\mathbf{f}_{l}$ and
$\mathbf{f}_{r}$, the correlation cost volume is computed for each disparity
level $d$,

$\textbf{C}_{corr}(d,x,y)=\frac{1}{N_{c}}\langle\mathbf{f}_{l}(x,y),\mathbf{f}_{r}(x-d,y)\rangle,$
(1)

where $\langle\cdot,\cdot\rangle$ is the inner product of two feature vectors
and $N_{c}$ denotes the number of channels. CRL [17] and iResNet [12] followed
the idea of DispNetC with stack refinement sub-networks to further improve the
performance. There are also works integrating additional information such as
edge features [22] and semantic features [26].

Recent works employed concatenation-based feature volume and 3D aggregation
networks for better context aggregation [6, 2, 28]. Kendall _et al_. proposed
GC-Net [6] and was the first to use 3D convolution networks to aggregate
context for cost volumes. Instead of directly giving a cost volume, the left
and the right feature $\mathbf{f}_{l}$, $\mathbf{f}_{r}$ are concatenated to
form a 4D feature volume,

$\textbf{C}_{concat}(d,x,y,\cdot)=\text{Concat}\left\\{\mathbf{f}_{l}(x,y),\mathbf{f}_{r}(x-d,y)\right\\}.$
(2)

Context features are aggregated from neighbour pixels and disparities with 3D
convolution networks to predict a disparity probability volume. Following GC-
Net, Chang _et al_. [2] proposed the pyramid stereo matching network (PSMNet)
with a spatial pyramid pooling module and stacked 3D hourglass networks for
cost volume refinement. Yu _et al_. [28] proposed to generate and select
multiple cost aggregation proposals. Zhong _et al_. [31] proposed a self-
adaptive recurrent stereo model to tackle open-world video data.

LRCR [5] utilized left-right consistency check and recurrent model to
aggregate cost volumes predicted from [21] and refined unreliable disparity
predictions. There are also other works focusing on real-time stereo matching
[7] and application friendly stereo [24].

Figure 1: The pipeline of the proposed group-wise correlation network. The
whole network consists of four parts, unary feature extraction, cost volume
construction, 3D convolution aggregation, and disparity prediction. The cost
volume is divided into two parts, concatenation volume (Cat) and group-wise
correlation volume (Gwc). Concatenation volume is built by concatenating the
compressed left and right features. Group-wise correlation volume is described
in Section 3.2. Figure 2: The structure of our proposed 3D aggregation
network. The network consists of a pre-hourglass module (four convolutions at
the beginning) and three stacked 3D hourglass networks. Compared with PSMNet
[2], we remove the shortcut connections between different hourglass modules
and output modules, thus output modules 0,1,2 can be removed during inference
to save time. $1{\times}1{\times}1$ 3D convolutions are added to the shortcut
connections within hourglass modules.

##  3 Group-wise Correlation Network

We propose group-wise correlation stereo network (GwcNet), which extends
PSMNet [2] with group-wise correlation cost volume and improved 3D stacked
hourglass networks. In PSMNet, the matching costs for the concatenated
features have to be learned from scratch by the 3D aggregation network, which
usually requires more parameters and computational cost. In contrast, full
correlation [14] provides an efficient way for measuring feature similarities
via dot product, but it loses much information. Our proposed group-wise
correlation overcomes both drawbacks and provides good features for similarity
measures.

###  3.1 Network architecture

The structure of the proposed group-wise correlation network is shown in
Figure 1. The network consists of four parts, unary feature extraction, cost
volume construction, 3D aggregation, and disparity prediction. The detailed
structures of the cost volume, stacked hourglass, and output modules are
listed in Table 1.

For feature extraction, we adopt the ResNet-like network used in PSMNet [2]
with the half dilation settings and without its spatial pyramid pooling
module. The last feature maps of conv2, conv3, and conv4 are concatenated to
form 320-channel unary feature maps.

The cost volume is composed of two parts, a concatenation volume and a group-
wise correlation volume. The concatenation volume is the same as PSMNet [2]
but with fewer channels, before which the unary features are compressed into
12 channels with two convolutions. The proposed group-wise correlation volume
will be described in details in Section 3.2. The two volumes are then
concatenated as the input to the 3D aggregation network.

The 3D aggregation network is used to aggregate features from neighboring
disparities and pixels to predict refined cost volumes. It consists of a pre-
hourglass module and three stacked 3D hourglass networks to regularize the
feature volumes. As shown in Figure 2, the pre-hourglass module consists of
four 3D convolutions with batch normalization and ReLU. Three stacked 3D
hourglass networks are followed to refine low-texture ambiguities and
occlusion parts by encoder-decoder structures. Compared with 3D aggregation
network of [2], we have several important modifications to improve the
performance and increase the inference speed, and details are described in
Section 3.3.

The pre-hourglass module and three stacked 3D hourglass networks are connected
to output modules. Each output module predicts a disparity map. The structure
of the output module and the loss function are described in Section 3.4.

###  3.2 Group-wise correlation volume

The left unary features and the right unary features are denoted by
$\mathbf{f}_{l}$ and $\mathbf{f}_{r}$ with $N_{c}$ channels and in $1{/}4$ the
size of original images. In previous works [14, 6, 2], the left and right
features are correlated or concatenated at different disparity levels to form
the cost volume. However, both correlation volume and concatenation volume
have drawbacks. The full correlation provides an efficient way for measuring
feature similarities, but it loses much information because it produces only a
single-channel correlation map for each disparity level. The concatenation
volume contains no information about the feature similarities, so more
parameters are required in the following aggregation network to learn the
similarity measurement function from scratch. To solve the above issues, we
propose group-wise correlation by combining advantages of the concatenation
volume and the correlation volume.

The basic idea behind group-wise correlation is splitting the features into
groups and computing correlation maps group by group. We denote the channels
of unary features as $N_{c}$. All the channels are evenly divided into $N_{g}$
groups along the channel dimension, and each feature group therefore has
$N_{c}/N_{g}$ channels. The $g$th feature group
$\mathbf{f}_{l}^{g},\mathbf{f}_{r}^{g}$ consists of the
$g\frac{N_{c}}{N_{g}},g\frac{N_{c}}{N_{g}}+1,\dots,g\frac{N_{c}}{N_{g}}+(\frac{N_{c}}{N_{g}}-1)$th
channels of the original feature $\mathbf{f}_{l},\mathbf{f}_{r}$. The group-
wise correlation is then computed as

$\textbf{C}_{gwc}(d,x,y,g)=\frac{1}{N_{c}/N_{g}}\langle\mathbf{f}_{l}^{g}(x,y),\mathbf{f}_{r}^{g}(x-d,y)\rangle,$
(3)

where $\langle\cdot,\cdot\rangle$ is the inner product. Note that the
correlation is computed for all feature groups $g$ and at all disparity levels
$d$. Then, all the correlation maps are packed into a matching cost volume of
the shape $[D_{max}/4,H/4,W/4,N_{g}]$, where $D_{max}$ denotes the maximum
disparity and $D_{max}/4$ corresponds to the maximum disparity for the
feature. When $N_{g}{=}1$, the group-wise correlation becomes full
correlation.

Group-wise correlation volume $\textbf{C}_{gwc}$ can be treated as $N_{g}$
cost volume proposals, and each proposal is computed from the corresponding
feature group. The following 3D aggregation network aggregates multiple
candidates to regress disparity maps. The group-wise correlation successfully
leverages the power of traditional correlation matching costs and provides
rich similarity-measure features for the 3D aggregation network, which
alleviates the parameter demand. We will show in Section 4.5 that we explore
to reduce the channels of the 3D aggregation network, and the performance
reduction of our proposed network is much smaller than [2]. Our proposed
group-wise correlation volume requires less 3D aggregation parameters to
achieve favorable results.

To further improve the performance, the group correlation cost volume can be
combined with the concatenation volume. Experiment results show that the
group-wise correlation volume and the concatenation volume are complementary
to each other.

Name |  Layer properties | Output size
---|---|---
Cost Volume
unary_l/r |  N/A, S2 | $\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}320$
volume_g |  group-wise cost volume | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}40$
volume_c |  concatenation cost volume | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}24$
volume |  volume_g,volume_c: Concat | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}64$
Pre-hourglass
conv1 |  [$32{\times}32$, $3{\times}3{\times}3$, S1] ${\times}2$ | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}32$
conv2 |  [$32{\times}32$, $3{\times}3{\times}3$, S1] ${\times}2$ | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}32$
output |  conv1,conv2: Add | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}32$
Hourglass Module 1, 2, 3
input |  N/A | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}32$
conv1a |  $32{\times}64$, $3{\times}3{\times}3$, S2 | $\nicefrac{{D}}{{8}}{\times}\nicefrac{{H}}{{8}}{\times}\nicefrac{{W}}{{8}}{\times}64$
conv1b |  $64{\times}64$, $3{\times}3{\times}3$, S1 | $\nicefrac{{D}}{{8}}{\times}\nicefrac{{H}}{{8}}{\times}\nicefrac{{W}}{{8}}{\times}64$
conv2a |  $64{\times}128$, $3{\times}3{\times}3$, S2 | $\nicefrac{{D}}{{16}}{\times}\nicefrac{{H}}{{16}}{\times}\nicefrac{{W}}{{16}}{\times}128$
conv2b |  $128{\times}128$, $3{\times}3{\times}3$, S1 | $\nicefrac{{D}}{{16}}{\times}\nicefrac{{H}}{{16}}{\times}\nicefrac{{W}}{{16}}{\times}128$
deconv1* |  $128{\times}64$, $3{\times}3{\times}3$, S2, deconv | $\nicefrac{{D}}{{8}}{\times}\nicefrac{{H}}{{8}}{\times}\nicefrac{{W}}{{8}}{\times}64$
shortcut1* |  conv1b: $64{\times}64$, $1{\times}1{\times}1$, S1 | $\nicefrac{{D}}{{8}}{\times}\nicefrac{{H}}{{8}}{\times}\nicefrac{{W}}{{8}}{\times}64$
plus1 |  deconv1,shortcut1: Add&ReLU | $\nicefrac{{D}}{{8}}{\times}\nicefrac{{H}}{{8}}{\times}\nicefrac{{W}}{{8}}{\times}64$
deconv0* |  $64{\times}32$, $3{\times}3{\times}3$, S2, deconv | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}32$
shortcut0* |  input: $32{\times}32$, $1{\times}1{\times}1$, S1 | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}32$
output |  deconv0,shortcut0: Add&ReLU | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}32$
Output Module 0, 1, 2, 3
input |  N/A | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}32$
conv1 |  $32{\times}32$, $3{\times}3{\times}3$, S1 | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}32$
conv2** |  $32{\times}1$, $3{\times}3{\times}3$, S1 | $\nicefrac{{D}}{{4}}{\times}\nicefrac{{H}}{{4}}{\times}\nicefrac{{W}}{{4}}{\times}1$
score |  Upsample | $D{\times}H{\times}W{\times}1$
prob |  Softmax (at disparity dimension) | $D{\times}H{\times}W{\times}1$
disparity |  Soft Argmin (Equ. 4) | $H{\times}W{\times}1$
Table 1: Structure details of the modules. $H,W$ represents the height and the
width of the input image. S1/2 denotes the convolution stride. If not
specified, each 3D convolution is with a batch normalization and ReLU. *
denotes the ReLU is not included. ** denotes convolution only.

###  3.3 Improved 3D aggregation module

In PSMNet [2], a stacked hourglass architecture was proposed to learn better
context features. Based on the network, we apply several important
modifications to make it suitable for our proposed group-wise correlation and
improve the inference speed. The structure of the proposed 3D aggregation is
shown in Figure 2 and Table 1.

First, we add one more auxiliary output module (output module 0, see Figure 2)
for the features of the pre-hourglass module. The extra auxiliary loss makes
the network learn better features at lower layers, which also benefits the
final prediction.

Second, the residual connections between different output modules are removed,
thus auxiliary output modules (output module 0, 1, 2) can be removed during
inference to save computational cost.

Third, $1{\times}1{\times}1$ 3D convolutions are added to the shortcut
connections within each hourglass module (see dashed lines in Figure 2) to
improve the performance without increasing much computational cost. Since the
$1{\times}1{\times}1$ 3D convolution only has 1/27 multiplication operations
compared with $3{\times}3{\times}3$ convolutions, it runs very fast and the
time can be neglected.

Model |  Concat Volume |  Group Corr Volume |  Stack Hourglass |  Groups ${\times}$ Channels |  Init Volume Channel |  >1px (%) |  >2px (%) |  >3px (%) |  EPE (px) |  Time (ms)
---|---|---|---|---|---|---|---|---|---|---
Cat64-Base |  ✓ |  |  |  - |  64 |  12.78 |  8.05 |  6.33 |  1.308 |  117.1
Gwc1-Base |  |  ✓ |  |  1${\times}$320 |  1 |  13.32 |  8.37 |  6.62 |  1.369 |  104.0
Gwc10-Base |  |  ✓ |  |  10${\times}$32 |  10 |  11.82 |  7.31 |  5.70 |  1.230 |  112.8
Gwc20-Base |  |  ✓ |  |  20${\times}$16 |  20 |  11.84 |  7.29 |  5.67 |  1.216 |  116.3
Gwc40-Base |  |  ✓ |  |  40${\times}$8 |  40 |  11.68 |  7.18 |  5.58 |  1.212 |  122.2
Gwc80-Base |  |  ✓ |  |  80${\times}$4 |  80 |  11.69 |  7.17 |  5.57 |  1.214 |  133.3
Gwc160-Base |  |  ✓ |  |  160${\times}$2 |  160 |  11.58 |  7.08 |  5.49 |  1.188 |  157.3
Gwc40-Cat24-Base |  ✓ |  ✓ |  |  40${\times}$8 |  40+24 |  11.26 |  6.87 |  5.31 |  1.127 |  135.1
PSMNet [2] |  ✓ |  |  [2] |  - |  64 |  9.46 |  5.19 |  3.80 |  0.887 |  246.1
Cat64-original-hg |  ✓ |  |  [2] |  - |  64 |  9.47 |  5.13 |  3.74 |  0.876 |  241.0
Cat64 |  ✓ |  |  Ours |  - |  64 |  8.41 |  4.63 |  3.41 |  0.808 |  198.3
Gwc40 (GwcNet-g) |  |  ✓ |  Ours |  40${\times}$8 |  40 |  8.18 |  4.57 |  3.39 |  0.792 |  200.3
Gwc40-Cat24 (GwcNet-gc) |  ✓ |  ✓ |  Ours |  40${\times}$8 |  40+24 |  8.03 |  4.47 |  3.30 |  0.765 |  210.7
Table 2: Ablation study results of proposed networks on the Finalpass of Scene Flow datasets [14]. Cat, Gwc, Gwc-Cat represent only concatenation volume, only group-wise correlation volume, or the both. Base denotes the network variants without stacked hourglass networks. The time is the inference time for 480${\times}$640 inputs on a single Nvidia TITAN Xp GPU. The result of PSMNet [2] is trained with published code with our batch size, evaluation settings for fair comparison. Model |  KITTI 12 EPE (px) |  KITTI 12 D1-all(%) |  KITTI 15 EPE (px) |  KITTI 15 D1-all (%)
---|---|---|---|---
PSMNet [2] |  0.713 |  2.53 |  0.639 |  1.50
Cat64-original-hg |  0.740 |  2.72 |  0.652 |  1.76
Cat64 |  0.691 |  2.41 |  0.615 |  1.55
Gwc40 |  0.662 |  2.30 |  0.602 |  1.41
Gwc40-Cat24 |  0.659 |  2.10 |  0.613 |  1.49
Table 3: Ablation study results of our networks on KITTI 2012 validation and
KITTI 2015 validation sets.

###  3.4 Output module and loss function

For each output module, two 3D convolutions are employed to output a 1-channel
4D volume, and then the volume is upsampled and converted into a probability
volume with softmax function along the disparity dimension. Detailed
structures are shown in Table 1. For each pixel, we have a $D_{max}$-length
vector which contains the probability $p$ for all disparity levels. Then, the
disparity estimation $\widetilde{d}$ is given by the soft argmin function [6],

$\widetilde{d}=\sum_{k=0}^{D_{max}-1}k\cdot p_{k},$ (4)

where $k$ and $p_{k}$ denote a possible disparity level and the corresponding
probability. The predicted disparity maps from the four output modules are
denoted as
$\mathbf{\widetilde{d}}_{0},\mathbf{\widetilde{d}}_{1},\mathbf{\widetilde{d}}_{2},\mathbf{\widetilde{d}}_{3}$.
The final loss is given by,

$L=\sum_{i=0}^{i=3}\lambda_{i}\cdot\text{Smooth}_{L_{1}}(\mathbf{\widetilde{d}}_{i}-\mathbf{d}^{*}),$
(5)

where $\lambda_{i}$ denotes the coefficients for the $i$th disparity
prediction and $\mathbf{d}^{*}$ represents the ground-truth disparity map. The
smooth L1 loss is computed as follows,

$\text{Smooth}_{L_{1}}(x)=\left\\{\begin{aligned} &0.5x^{2},&\text{if
}|x|<1\\\ &|x|-0.5,&\text{otherwise}\\\ \end{aligned}\right.$ (6)

##  4 Experiment

In this section, we evaluate our proposed stereo models on Scene Flow datasets
[14] and the KITTI dataset [3, 16]. We present ablation studies to compare
different models and different parameter settings. Datasets and implementation
details are described in Section 4.1 and Section 4.2. The effectiveness and
the best settings of group-wise correlation are explored in Section 4.3. The
performance improvement of the new stacked hourglass module is discussed in
Section 4.4. We also explore the performance of group-wise correlation when
the computational cost is limited in Section 4.5.

###  4.1 Datasets and evaluation metrics

Scene Flow datasets are a dataset collection of synthetic stereo datasets,
consisting of Flyingthings3D, Driving, and Monkaa. The datasets provide 35,454
training and 4,370 testing images of size 960${\times}$540 with accurate
ground-truth disparity maps. We use the Finalpass of the Scene Flow datasets,
since it contains more motion blur and defocus and is more like real-world
images than the Cleanpass. KITTI 2012 and KITTI 2015 are driving scene
datasets. KITTI 2012 provides 194 training and 195 testing images pairs, and
KITTI 2015 provides 200 training and 200 testing image pairs. Both datasets
provide sparse LIDAR ground-truth disparity for the training images.

For Scene Flow datasets, the evaluation metrics is usually the end-point error
(EPE), which is the mean average disparity error in pixels. For KITTI 2012,
percentages of erroneous pixels and average end-point errors for both non-
occluded (Noc) and all (All) pixels are reported. For KITTI 2015, the
percentage of disparity outliers D1 is evaluated for background, foreground,
and all pixels. The outliers are defined as the pixels whose disparity errors
are larger than $\text{max}(3\text{px},0.05d^{*})$, where $d^{*}$ denotes the
ground-truth disparity.

###  4.2 Implementation details

Our network is implemented with PyTorch. We use Adam [8] optimizer, with
$\beta_{1}=0.9$, $\beta_{2}=0.999$. The batch size is fixed to 16, and we
train all the networks with 8 Nvidia TITAN Xp GPUs with 2 training samples on
each GPU. The coefficients of four outputs are set as $\lambda_{0}=0.5$,
$\lambda_{1}=0.5$, $\lambda_{2}=0.7$, $\lambda_{3}=1.0$.

For Scene Flow datasets, we train the stereo networks for 16 epochs in total.
The initial learning rate is set to 0.001. It is down-scaled by 2 after epoch
10, 12, 14 and ends at 0.000125. To test on Scene Flow datasets, the full
images of size 960${\times}$540 are input to the network for disparity
prediction. We set the maximum disparity value as $D_{max}=192$ following
PSMNet [2] for Scene Flow datasets. To evaluate our networks, we remove all
the images with less than 10% valid pixels ($0{\leq}d{<}D_{max}$) in the test
set. For each valid image, the evaluation metrics are computed with only valid
pixels.

For KITTI 2015 and KITTI 2012, we fine-tune the network pre-trained on Scene
Flow datasets for another 300 epochs. The initial learning rate is 0.001 and
is down-scaled by 10 after epoch 200. For testing on KITTI datasets, we first
pad zeros on the top and the right side of the images to make the inputs in
size 1248${\times}$384.

###  4.3 The effectiveness of Group-wise correlation

In this section, we explore the effectiveness and the best settings for the
group-wise correlation. In order to prove the effectiveness of the proposed
group-wise correlation volume, we conduct several experiments on the Base
model, which removes the stacked hourglass networks and only preserves the
pre-hourglass module and the output module 0. Cat-Base, Gwc-Base, and Gwc-Cat-
Base are the base models with only concatenation volume, only group-wise
correlation volume, or both volumes.

Experiment results in Table 2 show that the performance of the Gwc-Base
network increases as the group number increases. When the group number is
larger than 40, the performance improvement becomes minor and the end-point
error stays around 1.2px. Considering the memory usage and the computational
cost, we choose 40 groups with each group having 8 channels as our network
structure, which corresponds to the Gwc40-Base model in Table 2.

All the Gwc-Base models except Gwc1-Base outperform the Cat-Base model which
utilizes concatenation volume, which shows the effectiveness of the group-wise
correlation. The Gwc40 model reduces the end-point error by 0.1px and the
3-pixel error rate by 0.75%, and the time consumption is almost the same. The
performance can be further improved by combining group-wise correlation volume
with concatenation volume (see Gwc40-Cat24-Base model in Table 2). The group-
wise correlation could provide accurate matching features, and the
concatenation volume provides complementary semantic information.

Figure 3: Our model Gwc-Cat achieves much better performance than Cat when the
number of channels decreases. The models with 32 base channels correspond to
the Cat64 model (concatenation volume) and the Gwc40-Cat24 model (group-wise
correlation and concatenation volume). The channels of the cost volume and all
3D convolutions decrease by the same factor as the base channel.

|  |
---|---|---
|  |
|  |
(a) Visualization results on the Scene Flow datasets.

|  |
---|---|---
|  |
|  |
(b) Visualization results on the KITTI 2012 dataset.

|  |
---|---|---
|  |
|  |
(c) Visualization results on the KITTI 2015 dataset.

Figure 4: Depth visualization results on the test sets of Scene Flow [14], KITTI 2012 [3] and KITTI 2015 [16] datasets. From left to right, input left images, predicted disparity maps, and error maps. | All (%) | Noc (%) |  Time
---|---|---|---
|  D1-bg |  D1-fg |  D1-all |  D1-bg |  D1-fg |  D1-all |  (s)
DispNetC [14] |  4.32 |  4.41 |  4.34 |  4.11 |  3.72 |  4.05 |  0.06
GC-Net [6] |  2.21 |  6.16 |  2.87 |  2.02 |  5.58 |  2.61 |  0.9
CRL [17] |  2.48 |  3.59 |  2.67 |  2.32 |  3.12 |  2.45 |  0.47
iResNet-i2e2 [12] |  2.14 |  3.45 |  2.36 |  1.94 |  3.20 |  2.15 |  0.22
PSMNet [6] |  1.86 |  4.62 |  2.32 |  1.71 |  4.31 |  2.14 |  0.41
SegStereo [26] |  1.88 |  4.07 |  2.25 |  1.76 |  3.70 |  2.08 |  0.6
GwcNet-g (Gwc40) |  1.74 |  3.93 |  2.11 |  1.61 |  3.49 |  1.92 |  0.32
Table 4: KITTI 2015 test set results. The dataset contains 200 images for training and 200 images for testing.  | >2px (%) | >3px (%) | >5px (%) | Mean Error (px) |  Time
---|---|---|---|---|---
|  Noc |  All |  Noc |  All |  Noc |  All |  Noc |  All |  (s)
DispNetC [14] |  7.38 |  8.11 |  4.11 |  4.65 |  2.05 |  2.39 |  0.9 |  1.0 |  0.06
MC-CNN-acrt [29] |  3.90 |  5.45 |  2.43 |  3.63 |  1.64 |  2.39 |  0.7 |  0.9 |  67
GC-Net [6] |  2.71 |  3.46 |  1.77 |  2.30 |  1.12 |  1.46 |  0.6 |  0.7 |  0.9
iResNet-i2 [12] |  2.69 |  3.34 |  1.71 |  2.16 |  1.06 |  1.32 |  0.5 |  0.6 |  0.12
SegStereo [26] |  2.66 |  3.19 |  1.68 |  2.03 |  1.00 |  1.21 |  0.5 |  0.6 |  0.6
PSMNet [6] |  2.44 |  3.01 |  1.49 |  1.89 |  0.90 |  1.15 |  0.5 |  0.6 |  0.41
GwcNet-gc (Gwc40-Cat24) |  2.16 |  2.71 |  1.32 |  1.70 |  0.80 |  1.03 |  0.5 |  0.5 |  0.32
Table 5: KITTI 2012 test set results. The dataset contains 194 images for
training and 195 images for testing.

###  4.4 Improved stacked hourglass

In this paper, we applied several modifications to the stacked hourglass
networks proposed in [2] to improve the performance of cost volume
aggregation. From Table 2 and Table 3, we can see that the model with the
proposed hourglass networks (Cat64) increases EPE by 7.8% on Scene Flow
datasets and 5.8% on KITTI 2015 compared with the model Cat64-original-hg
(with the hourglass module in [2]). The inference time for $640{\times}480$
inputs on a single Nvidia TITAN Xp GPU also decreases by 42.7ms, because the
auxiliary output modules can be removed during inference to save time.

###  4.5 Limit the computational cost of 3D network

We explore to limit the computational cost by decreasing channels in the 3D
aggregation network to verify the effectiveness of the proposed group-wise
group correlation. The results are shown in Figure 3. The base number of
channels are modified from the original 32 to 2, and the channels of the cost
volume and all 3D convolutions are reduced with the same factor. As the number
of channels decreasing, our models with group-wise correlation volume (Gwc-
Cat) perform much better than the models with only concatenation volume (Cat).
The performance gain enlarges as more channels reduced. The reason for this is
that the group-wise correlation provides good matching cost representations
for the 3D aggregation network, while the aggregation network with only
concatenation volume as inputs needs to learn the matching similarity function
from scratch, which usually requires more parameters and computational cost.
As a result, the proposed group-wise correlation could be a valuable method to
be implemented in real-time stereo networks where the computational costs are
limited.

###  4.6 KITTI 2012 and KITTI 2015

For KITTI stereo 2015 [16], we split the training set into 180 training image
pairs and 20 validation image pairs. Since the results on the validation set
are not stable, we fine-tune the pretrained model for 3 times and choose the
model with the best validation performance. From Table 3, the performance of
both Gwc40-Cat24 and Gwc40 is better than the models without group-wise
correlation (Cat64, Cat64-original-hg). We submit the Gwc40 model (without
concatenation volume) with the lowest validation error to the evaluation
server, and the results on the test set are shown in Table 4. Our model
surpasses the PSMNet [2] by 0.21% and SegStereo [26] by 0.14% on D1-all.

For KITTI 2012 [3], we split the training set into 180 training images and 14
validation image pairs. The results on the validation set are shown in Table
3. We submit the best Gwc40-Cat24 model on the validation set to the
evaluation server. The evaluation results on the test set are shown in Table
5. Our method surpasses PSMNet [2] by 0.19% on 3-pixel-error and 0.1px on mean
disparity error.

##  5 Conclusion

In this paper, we proposed GwcNet to estimate disparity maps for stereo
matching, which incorporates group-wise correlation to build up the cost
volumes. The group-wise correlation volumes provide good matching features for
the 3D aggregation network, which improves the performance and reduces the
parameter requirements of the aggregation network. We showed that when the
computational cost is limited, our model achieves larger gain than previous
concatenation-volume based stereo networks. We also improved the stacked
hourglass networks to further improve the performance and reduce the inference
time. Experiments demonstrated the effectiveness of our proposed method on the
Scene Flow datasets and the KITTI dataset.

## References

  * [1] K. Batsos, C. Cai, and P. Mordohai.  Cbmv: A coalesced bidirectional matching volume for disparity estimation.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2060–2069, 2018.
  * [2] J.-R. Chang and Y.-S. Chen.  Pyramid stereo matching network.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5410–5418, 2018.
  * [3] A. Geiger, P. Lenz, and R. Urtasun.  Are we ready for autonomous driving? the kitti vision benchmark suite.  In Conference on Computer Vision and Pattern Recognition (CVPR), 2012.
  * [4] H. Hirschmuller.  Accurate and efficient stereo processing by semi-global matching and mutual information.  In Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on, volume 2, pages 807–814. IEEE, 2005.
  * [5] Z. Jie, P. Wang, Y. Ling, B. Zhao, Y. Wei, J. Feng, and W. Liu.  Left-right comparative recurrent model for stereo matching.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3838–3846, 2018.
  * [6] A. Kendall, H. Martirosyan, S. Dasgupta, P. Henry, R. Kennedy, A. Bachrach, and A. Bry.  End-to-end learning of geometry and context for deep stereo regression.  In Proceedings of the IEEE International Conference on Computer Vision, pages 66–75, 2017.
  * [7] S. Khamis, S. Fanello, C. Rhemann, A. Kowdle, J. Valentin, and S. Izadi.  Stereonet: Guided hierarchical refinement for real-time edge-aware depth prediction.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 573–590, 2018.
  * [8] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [9] A. Klaus, M. Sormann, and K. Karner.  Segment-based stereo matching using belief propagation and a self-adapting dissimilarity measure.  In 18th International Conference on Pattern Recognition (ICPR’06), volume 3, pages 15–18. IEEE, 2006.
  * [10] P. Knobelreiter, C. Reinbacher, A. Shekhovtsov, and T. Pock.  End-to-end training of hybrid cnn-crf models for stereo.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2339–2348, 2017.
  * [11] V. Kolmogorov and R. Zabih.  Computing visual correspondence with occlusions using graph cuts.  In Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001, volume 2, pages 508–515 vol.2, July 2001.
  * [12] Z. Liang, Y. Feng, Y. G. H. L. W. Chen, and L. Q. L. Z. J. Zhang.  Learning for disparity estimation through feature constancy.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2811–2820, 2018.
  * [13] W. Luo, A. G. Schwing, and R. Urtasun.  Efficient deep learning for stereo matching.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5695–5703, 2016.
  * [14] N. Mayer, E. Ilg, P. Hausser, P. Fischer, D. Cremers, A. Dosovitskiy, and T. Brox.  A large dataset to train convolutional networks for disparity, optical flow, and scene flow estimation.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 4040–4048, 2016.
  * [15] X. Mei, X. Sun, W. Dong, H. Wang, and X. Zhang.  Segment-tree based cost aggregation for stereo matching.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 313–320, 2013.
  * [16] M. Menze and A. Geiger.  Object scene flow for autonomous vehicles.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3061–3070, 2015.
  * [17] J. Pang, W. Sun, J. S. Ren, C. Yang, and Q. Yan.  Cascade residual learning: A two-stage convolutional neural network for stereo matching.  In ICCV Workshops, volume 7, 2017.
  * [18] D. Scharstein and R. Szeliski.  A taxonomy and evaluation of dense two-frame stereo correspondence algorithms.  International journal of computer vision, 47(1-3):7–42, 2002.
  * [19] J. L. Schonberger, S. N. Sinha, and M. Pollefeys.  Learning to fuse proposals from multiple scanline optimizations in semi-global matching.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 739–755, 2018.
  * [20] A. Seki and M. Pollefeys.  Sgm-nets: Semi-global matching with neural networks.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 231–240, 2017.
  * [21] A. Shaked and L. Wolf.  Improved stereo matching with constant highway networks and reflective confidence learning.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 4641–4650, 2017.
  * [22] X. Song, X. Zhao, H. Hu, and L. Fang.  Edgestereo: A context integrated residual pyramid network for stereo matching.  In Asian Conference on Computer Vision (ACCV), 2018.
  * [23] J. Sun, N.-N. Zheng, and H.-Y. Shum.  Stereo matching using belief propagation.  IEEE Transactions on pattern analysis and machine intelligence, 25(7):787–800, 2003.
  * [24] S. Tulyakov, A. Ivanov, and F. Fleuret.  Practical deep stereo (pds): Toward applications-friendly deep stereo matching.  In Advances in Neural Information Processing Systems, pages 5875–5885, 2018.
  * [25] Y. Wu and K. He.  Group normalization.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 3–19, 2018.
  * [26] G. Yang, H. Zhao, J. Shi, Z. Deng, and J. Jia.  Segstereo: Exploiting semantic information for disparity estimation.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 636–651, 2018.
  * [27] Q. Yang.  A non-local cost aggregation method for stereo matching.  In 2012 IEEE Conference on Computer Vision and Pattern Recognition, pages 1402–1409. IEEE, 2012.
  * [28] L. Yu, Y. Wang, Y. Wu, and Y. Jia.  Deep stereo matching with explicit cost aggregation sub-architecture.  In Thirty-Second AAAI Conference on Artificial Intelligence, 2018\.
  * [29] J. Zbontar and Y. LeCun.  Computing the stereo matching cost with a convolutional neural network.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1592–1599, 2015.
  * [30] K. Zhang, J. Lu, and G. Lafruit.  Cross-based local stereo matching using orthogonal integral images.  IEEE transactions on circuits and systems for video technology, 19(7):1073–1079, 2009.
  * [31] Y. Zhong, H. Li, and Y. Dai.  Open-world stereo video matching with deep rnn.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 101–116, 2018.


# Direct kinetic Alfvén wave energy cascade in the presence of imbalance

T. Passot  P.L. Sulem and D. Laveder  Université Côte d’Azur, Observatoire de
la Côte d’Azur, CNRS, Laboratoire J.L. Lagrange, Boulevard de l’Observatoire,
CS 34229, 06304 Nice Cedex 4, France

###### Abstract

A two-field Hamiltonian gyrofluid model for kinetic Alfvén waves retaining ion
finite Larmor radius corrections and parallel magnetic field fluctuations is
used to study direct turbulent cascades from the MHD to the sub-ion scales.
For moderate energy imbalance and weak enough magnetic fluctuations, the
spectrum of the transverse magnetic field and that of the most energetic wave
display a steep transition zone near the ion scale, while the parallel
transfer (and thus the parallel dissipation) remains weak. In this regime, the
perpendicular flux of generalized cross helicity displays a significant decay
past the ion scale, while the perpendicular energy flux remains almost
constant. A phenomenological model suggests that the interactions between co-
propagative waves present at the sub-ion scales can play a central role in the
development of a transition zone in the presence of a helicity barrier.

##  1 Introduction

Alfvén wave (AW) turbulence, that plays an important role in magnetized
collisionless plasmas such as the solar wind (Bruno & Carbone, 2013, 2016), is
often imbalanced in the sense that the energy carried by the waves propagating
in the forward and backward directions relatively to the ambient field have
significantly different magnitudes (Tu et al., 1989; Lucek & Balogh, 1998;
Wicks et al., 2013). The degree of imbalance is dependent on the type of wind
(Tu et al., 1990) and also on the distance from the Sun (Roberts et al., 1987;
Marsch & Tu, 1990; Chen et al., 2020). For example, in a classical model of
solar wind, outgoing AWs are emitted at the Sun’s surface, as a consequence of
reconnection processes in the chromospheric magnetic network, while in-going
waves result from reflection of these waves on density gradients (Perez &
Chandran, 2013; Chandran & Perez, 2019), or from parametric decay instability
(Malara et al., 2001). In the case of low-frequency highly transverse
perturbations, corresponding to the assumptions underlying the gyrofluid model
derived in Passot et al. (2018), imbalance is associated with a non-zero value
of an invariant referred to as generalized cross helicity (GCH), which
identifies with the cross-helicity at the magnetohydrodynamic (MHD) scales and
with the magnetic helicity at sub-ion scales. As cross-helicity cascades to
small scales and magnetic helicity to large scales (Miloshevich et al., 2021),
the questions arises of the behaviour of the GCH flux arriving from large
scales to the transition zone between MHD and sub-ion scales. Two publications
recently addressed this issue. One uses a gyrofluid model (in the limit of
small values of the electron beta parameter) which isolates the AW dynamics
(Meyrand et al., 2021). The other is based on a hybrid-kinetic PIC simulation
which enables the coupling of low-frequency quasi transverse AW with quasi-
parallel high-frequency ion-cyclotron waves (ICW) (Squire et al., 2021). The
role of ICWs and of the dissipation associated with the proton cyclotron
resonance at the spectral break separating the MHD and the sub-ion ranges was
indeed stressed, using the data from the Messenger and Wind spacecraft
missions (Telloni et al., 2015; Woodham et al., 2018).

In the regimes studied by Meyrand et al. (2021), the nonlinear parameter,
defined as the ratio of the inverse nonlinear time to the Alfvén wave
frequency is typically above unity. The dynamics then rapidly develops a
significant transfer of energy towards the large parallel wavenumbers,
resulting in the breaking of the asymptotics underlying gyrofluid models.
Despite this shortcoming, the gyrofluid simulation reveals the existence of a
helicity barrier which depletes the energy flux towards the small
perpendicular scales, leading to a transfer in the parallel direction at the
ion scale, together with a steepening of the energy (and magnetic) spectra.
Energy accumulates at large scales and imbalance increases until saturation
occurs under the action of a strong parallel dissipation. In contrast with
usual fluid turbulence, the energy level that is reached depends on the
magnitude of the parallel viscosity coefficient. In this framework,
dissipation, which plays a central role in the dynamics, is not associated
with any realistic physical phenomenon. The above scenario was then validated
by the kinetic simulation of Squire et al. (2021), where the dissipation
originates from ion-cyclotron resonance. In this simulation, performed at
$\beta_{e}=0.3$, the level of fluctuations remains moderate because cyclotron
resonance is rapidly triggered. In fact, an argument based on critical
balance, together with the requirement that $k_{\|}d_{i}=1$ (where $d_{i}$ is
the ion inertial length) at perpendicular scales such that
$k_{\perp}\rho_{i}=1$ ($\rho_{i}$ denoting the ion Larmor radius), leads to a
saturation amplitude scaling like $\beta_{i}^{1/2}$. As a consequence, when
$\beta_{e}$ is taken of order unity (for comparable ion and electron
temperatures), one may expect that high amplitude fluctuations should develop
before a stationary regime can establish. Recent observations using Parker
Solar Probe have comforted the results of the above simulation, with evidence
of the coincidence between the presence of ICWs and of magnetic spectra
displaying an ion transition zone associated with a polarity change, as
expected when Alfvén waves trigger the formation of ICWs (Bowen et al., 2021).

Transition zones are commonly observed in the solar wind at 1AU, where the
plasma beta parameter is of order unity or larger (see e.g. Huang et al.
(2021) and references therein). When the amplitude of the fluctuations is
moderate (such as in the fast wind), the above argument suggests that ICWs can
hardly be excited at the perpendicular ion scale. In this context, an
alternative mechanism should be at work, like for example the effect of co-
propagating Alfvén wave interactions that become relevant at ion scales
(Voitenko & De Keyser, 2016; Gogoberidze & Voitenko, 2020). In order to
address this question, it is appropriate to use a two-field gyrofluid model
(Passot et al., 2018), where driving and dissipation have been supplemented,
which only retains the Alfvén waves, but requires that the fluctuation
amplitudes remain small and the parallel dynamics sub-dominant (conditions
that can be checked a posteriori). We chose to prescribe the ”energy
imbalance” (referred to as imbalance in the following) by maintaining constant
the amplitudes of the largest modes retained in the simulations, a driving
process classically used in hydrodynamics turbulence (Siggia & Patterson,
1978). It turns out that, in the quasi-stationary phases of the simulation,
the GCH and energy injection rates are roughly constant (equal to the total
dissipation rates), the ”injection imbalance”, measured by their ratio,
resulting from the dynamics.

The paper is organized as follows. Section 2 describes the model and the
numerical setup. Section 3 discusses the gyrofluid simulations, showing that,
for moderate imbalance, the spectrum of the magnetic fluctuations (and also
that of the most energetic wave) displays a steep transition zone near the ion
scale, while the parallel transfer (and thus the parallel dissipation) remains
weak. The helicity barrier is mostly conspicuous on the perpendicular GCH
flux, which displays a significant decay past the ion scale, while the
transverse energy flux remains almost constant. In the situation where
parallel energy and cross-helicity dissipations are negligible, the effect of
co-propagating wave interactions seems to be the only alternative mechanism at
the origin of the transition zone within the context of the present model.
Section 4 discusses this mechanism using a Leith-type phenomenological
spectral model. Section 5 is the conclusion.

##  2 Equations and numerical setup

Let us consider a situation where the plasma dynamics results from weak
perturbations of a homogeneous equilibrium state characterized by a density
$n_{0}$, isotropic ion and electron temperatures $T_{0i}$ and $T_{0e}$
(characterized by their ratio $\tau=T_{0i}/T_{0e}$), and subject to an ambient
magnetic field of amplitude $B_{0}$ along the $z$-direction. The various
characteristic scales are conveniently measured in terms of the sonic Larmor
radius $\rho_{s}=c_{s}/\Omega_{i}$, where $c_{s}=\sqrt{T_{0e}/m_{i}}$ is the
sound speed and $\Omega_{i}=eB_{0}/(mc)$ the ion gyrofrequency. In particular,
the ion Larmor radius $\rho_{i}=v_{th\,i}/\Omega_{i}$ (where the ion thermal
velocities is given by $v_{th,i}=(2T_{0i}/m_{i})^{1/2}$) and the ion inertial
length $d_{i}=v_{A}/\Omega_{i}$ (where $v_{A}=B_{0}/(4\upi
n_{0}m_{i})^{1/2}=c_{s}\sqrt{2/\beta_{e}}$ is the Alfvén velocity), satisfy
$\rho_{i}=\sqrt{2\tau}\rho_{s}$ and $d_{i}=\sqrt{2/{\beta_{e}}}\rho_{s}$.

The two-field gyrofluid model is derived from the full gyrofluid model of
Brizard (1992) under specific scaling assumptions. Such a simple reduced fluid
model does not retain physical phenomena like kinetic dissipation, but also
misses some nonlinear couplings with other modes. Discarding these effects,
one here focuses on the nonlinear dynamics of AWs, or kinetic Alfvén waves
(KAWs). Fully kinetic simulations will ultimately be necessary to validate the
findings. Like the parent model, the model is not restricted to the small
$\beta_{e}$ regime, unless physics at the electron inertia scale $d_{e}$ is
retained (in the absence of electron FLR terms in the two-field equations, the
electron Larmor radius must indeed remain smaller than $d_{e}$). More
precisely, leaving aside issues related to kinetic damping, the model
reproduces Reduced MHD in the limit of large transverse scales, and Electron
Reduced MHD at sub-ion scales, both systems describing AWs or KAWs, uncoupled
to slow waves, with no restriction on $\beta_{e}$ (especially since the
parallel magnetic field component is retained). Only at the ion scales, where
the system includes a dispersive term associated with the electron pressure
gradient, can a coupling between Alfvén and slow waves become relevant (see
Appendix E of Schekochihin et al. (2009)). We here neglect the electron
inertia, so the model covers a spectral range extending from the MHD to the
sub-ion scales, within an anisotropic scaling where the transverse scales are
much smaller than the parallel ones.

###  2.1 The two-field gyrofluid model

The model (which retains no dissipation process) is written as equations for
the electron gyrocenter number density $N_{e}$ and the parallel component of
the magnetic potential $A_{\|}$. When electron inertia is neglected, it takes
the form

$\displaystyle\partial_{t}N_{e}+[\varphi,N_{e}]-[B_{z},N_{e}]+\frac{2}{\beta_{e}}\nabla_{\|}\Delta_{\perp}A_{\|}=0$
(1) $\displaystyle\partial_{t}A_{\|}+\nabla_{\|}(\varphi-N_{e}-B_{z})=0.$ (2)

Here, $\Delta_{\perp}=\partial_{xx}+\partial_{yy}$ is the Laplacian in the
plane transverse to the ambient field and
$[f,g]=\partial_{x}f\partial_{y}g-\partial_{y}f\partial_{x}g$ the canonical
bracket of two scalar functions $f$ and $g$. Furthermore, $\Gamma_{n}$ denotes
the (non-local) operator $\Gamma_{n}(-\tau\Delta_{\perp})$ which in Fourier
space reduces to the multiplication by the function $\Gamma_{n}(\tau
k_{\perp}^{2})$, defined by $\Gamma_{n}(x)=I_{n}(x)e^{-x}$ where $I_{n}$ is
the modified Bessel function of first type of order n. For a scalar function
$f$, the parallel gradient operator $\nabla_{\parallel}$ is defined by
$\nabla_{\parallel}f=-[A_{\parallel},f]+\partial_{z}f$.

The equations are written in a nondimensional form, using the following units:
$\Omega_{i}^{-1}$ for time, $\rho_{s}$ for the space coordinates (and thus
$\rho_{s}^{-1}$ for the wavenumber components), $B_{0}$ for the parallel
magnetic fluctuations $B_{z}$, $n_{0}$ for the electron gyrocenter density
$N_{e}$, $T_{e}/e$ for the electric potential $\varphi$ and $B_{0}\rho_{s}$
for the parallel magnetic potential $A_{\|}$.

Introducing the operators
$L_{1}={2}/{\beta_{e}}+(1+2\tau)(\Gamma_{0}-\Gamma_{1})$,
$L_{2}=1+(1-\Gamma_{0})/\tau-\Gamma_{0}+\Gamma_{1}$,
$L_{3}=(1-\Gamma_{0})/\tau$ and $L_{4}=1-\Gamma_{0}+\Gamma_{1}$, the parallel
magnetic fluctuations are given by $B_{z}=M_{1}\varphi$, with
$M_{1}=L_{1}^{-1}L_{2}$, and the electrostatic potential by
$\varphi=-M_{2}^{-1}N_{e}$, where $M_{2}=L_{3}+L_{4}L_{1}^{-1}L_{2}$ is
positive definite. Thus, $B_{z}$ and $\varphi$ can both be expressed in terms
of $N_{e}$.

At the linear level, the phase velocity $v_{ph}=\omega/k_{z}$ is given by

$v_{ph}^{2}\equiv\left(\frac{\omega}{k_{z}}\right)^{2}=\frac{2}{\beta_{e}}k_{\perp}^{2}\frac{1-{\widehat{M}}_{1}+{\widehat{M}}_{2}}{{\widehat{M}}_{2}},$
(3)

where the caret refers to the Fourier symbol of the operator. The associated
operator $V_{ph}$ is given by

$V_{ph}=s(-\Delta_{\perp})^{1/2}(1-M_{1}+M_{2})^{1/2}M_{2}^{-1/2},$ (4)

where $s=(2/\beta_{e})^{1/2}$ is the equilibrium Alfvén velocity in sound
speed units. In physical space, the eigenmodes, that can be viewed as
generalized Elsasser potentials, are given by $\mu^{\pm}=\Lambda\varphi\pm
sA_{\|}$, where
$\Lambda=(-\Delta_{\perp})^{-1}(1+M_{2}-M_{1})^{1/2}M_{2}^{1/2}$.

In the absence of dissipation processes (needed in the turbulent regime), the
system (1)-(2) preserves the energy ${\mathcal{E}}$ and the generalized cross-
helicity ${\mathcal{C}}$,

$\displaystyle{\mathcal{E}}=\frac{1}{2}\int\Big{(}\frac{2}{\beta_{e}}|\bnabla_{\perp}A_{\|}|^{2}+\frac{4\delta^{2}}{\beta_{e}^{2}}|\Delta_{\perp}A_{\|}|^{2}-N_{e}(\varphi-
N_{e}-B_{z})\Big{)}d^{3}{x},$ (5) $\displaystyle{\mathcal{C}}=-\int
N_{e}A_{\|}d^{3}{x}.$ (6)

In terms of $\mu^{\pm}$, the two invariants rewrite ($D=(-\Delta)^{1/2}$),

$\displaystyle{\mathcal{E}}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\int\left\\{(D\mu^{+})^{2}+(D\mu^{-})^{2}\right\\}d^{3}x$
(7) $\displaystyle{\mathcal{C}}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\int
V_{ph}^{-1}\\{D\mu^{+})^{2}-(D\mu^{-})^{2}\\}d^{3}x.$ (8)

###  2.2 Numerical setup

Equations (1)–(2) with $\beta_{\rm e}=2$ and $\tau=1$ were solved in a
isotropic periodic box of size $L$, using a spectral method for the space
variables and a third-order Runge-Kutta scheme for the time stepping. An
eight-order hyper-dissipation operator
$\nu_{\perp}(\Delta_{\perp})^{4}+\nu_{z}\partial^{8}_{z}$ acting on $N_{e}$ or
on $A_{\|}$ is supplemented in the equations for the corresponding quantities.
For convenience, the parameters $\nu_{\perp}$ and $\nu_{\|}$ are referred to
as transverse (or perpendicular) and parallel viscosities. The system is
driven by freezing the amplitude of the modes whose transverse wavenumber
stands in the first spectral shell and the parallel one corresponds to
$k_{z}(L/2\upi)=\pm 1$, as in Siggia & Patterson (1978). A main interest of
this procedure is that it permits simulations with a high imbalance and a
moderate nonlinearity parameter. In such a regime, the injection rates of
energy and of GCH are not prescribed, in contrast with the situations where
the system is driven by an external random force or by a negative damping
(Meyrand et al., 2021).

| $N$ | $L/2\upi$ | I | $\chi_{0}$ | $\chi_{\rm final}$ | $\nu_{perp}$ | $\nu_{z}/\nu_{perp}$ | $\Delta t$ | Color code
---|---|---|---|---|---|---|---|---|---
$R_{1}$ |  $784^{3}$ |  11 |  100 |  0.47 |  0.61 |  $3\cdot 10^{-9}$ |  1 |  $2\cdot 10^{-3}$ |
$R_{2}$ |  $280^{3}$ |  5.5 |  100 |  0.50 |  0.68 |  $5\cdot 10^{-9}$ |  1 |  $3.125\cdot 10^{-3}$ | red
$R_{3}$ |  $280^{3}$ |  5.5 |  10 |  0.43 |  0.60 |  $2\cdot 10^{-8}$ |  1 |  $3.125\cdot 10^{-3}$ | green
$R_{4}$ |  $280^{3}$ |  5.5 |  1 |  0.36 |  0.56 |  $5\cdot 10^{-8}$ |  1 |  $3.125\cdot 10^{-3}$ | orange
$R_{5}$ |  $280^{3}$ |  5.5 |  100 |  0.87 |  1.2 |  $2\cdot 10^{-8}$ |  1 |  $3.125\cdot 10^{-3}$ | magenta
$R_{6}$ |  $280^{3}$ |  5.5 |  100 |  0.50 |  0.45 |  $5\cdot 10^{-9}$ |  0.025 |  $3.125\cdot 10^{-3}$ | blue
$R_{7}$ |  $280^{3}$ |  5.5 |  100 |  0.50 |  0.43 |  $5\cdot 10^{-9}$ |  0 |  $3.125\cdot 10^{-3}$ | cyan
$R_{8}$ |  $720^{3}$ |  5.5 |  1400 |  0.49 |  0.55 |  $2\cdot 10^{-10}$ |  1 |  $8\cdot 10^{-4}$ | violet
Table 1: Parameters of the runs, together with the values $\chi_{0}$ and
$\chi_{\rm final}$ of the nonlinear parameter at the beginning and at the end
of the simulation, estimated by the square root of the magnetic energy at this
time. The color code refers to the color of the lines in figures 7-11.

##  3 The GCH barrier

###  3.1 Spectral transition zone and arrest of the GCH cascade

In addition to the imbalance given by the ratio
$I={\mathcal{E}}^{+}/{\mathcal{E}}^{+}$ of the energies of the forward and
backward propagating waves, another quantity governing the dynamics is
provided by the nonlinearity parameter $\chi$, defined as the ratio of the
inverse nonlinear time to the Alfvén frequency. At early times, when the
maintained modes are dominant, a simple phenomenological argument (Miloshevich
et al., 2020) leads to $\chi\approx(k_{\perp 0}/k_{z0})|B_{\perp}(k_{\perp
0})|=|B_{\perp}(k_{\perp 0})|$, where $B_{\perp}$ is the transverse magnetic
fluctuations (measured in units of the ambient field), and the index zero
refers to the maintained modes for which $k_{\perp 0}=k_{z0}$ in the
(anisotropic) units of the gyrofluid model. At later times, in the isotropic
simulation box, the rms value of the magnetic fluctuation is used as an
estimation of $\chi$. In this paper, we concentrate on the case where the
nonlinearity parameter keeps moderate values (initially $\chi\leq 0.5$),
unless otherwise stated.

Figure 1: Transverse energy spectrum $E^{\pm}$ of the forward and backward
propagating waves (left panel) and energy spectrum of the transverse and
parallel magnetic fluctuations $E_{B_{\perp}}$ and $E_{B_{z}}$, averaged on
three Alfvén crossing times $L/v_{A}$, for run $R_{1}$.

Figure 2: Energy and GCH fluxes in the transverse direction for run $R_{1}$
averaged in the same time interval as in Fig. 1.

Figure 3: Parallel energy (left panel) and GCH (right panel) dissipations
${\mathcal{D}}_{z}({\mathcal{E}})$ and ${\mathcal{D}}_{z}({\mathcal{C}})$
respectively, for run $R_{1}$. The run was started using equal perpendicular
and parallel viscosity coefficients $\nu=5\times 10^{-8}$ using a resolution
of $420^{3}$ collocation points (red line); a first restart was then performed
with $\nu=5\times 10^{-9}$ and $560^{3}$ collocation points (green line), and
another one with $\nu=3\times 10^{-10}$ and $784^{3}$ collocation points (blue
line). The inserts show the corresponding perpendicular energy and GCH
dissipations ${\mathcal{D}}_{\perp}({\mathcal{E}})$ and
${\mathcal{D}}_{\perp}({\mathcal{C}})$. Due to the lack of space, the time
labels were not specified for the inserts that cover the same time interval as
the main graphs.

We shall first illustrate our main results on a large-resolution simulation
(run $R_{1}$) that displays extended inertial ranges. At the end of this
simulation, the system did not yet reach a stationary state. Nevertheless, as
mentioned in section 3.2, the observed dynamics is similar to the one found in
a lower-resolution simulation (run $R_{2}$), carried out until saturation is
reached. For the $R_{1}$ simulation characterized by an initial imbalance
$I=100$ and an initial value $\chi_{0}=0.47$ of the nonlinearity parameter,
performed in a computational box of size $L=11\times 2\upi$, figure 1 displays
the transverse spectra (integrated over $k_{z}$) $E^{\pm}(k_{\perp})$ of the
forward and backward propagating waves (left panel), and
$E_{B_{\perp}}(k_{\perp})$, $E_{B_{z}}(k_{\perp})$ of the transverse and
parallel magnetic fluctuations (right panel), averaged over the time interval
[14656, 14760] corresponding to $1.5$ Alfvén crossing times (see figure 3 for
a characterization of this interval relatively to the full simulation). The
spectra $E^{+}(k_{\perp})$ and $E_{B_{\perp}}(k_{\perp})$ display a transition
zone at scales slightly larger than the sonic Larmor radius (corresponding to
$k_{\perp}=1$), with a spectral index close to $-3$, located between a sub-ion
range with a spectral index $-2.7$ and a MHD range whose extension is strongly
limited by the numerical constraints which make difficult the simulation in a
box that would ideally cover a more extended range at large scale, while also
retaining the small-scale dynamics. The transition zone is seen to extend up
to the wavenumber where the spectrum $E_{B_{z}}(k_{\perp})$ changes from a
quasi-flat range to a decay range displaying the same spectral exponent as
$E_{B_{\perp}}(k_{\perp})$. The large-scale $B_{z}$ fluctuations are much
smaller than those of $B_{\perp}$ as expected when only Alfvén waves are
excited at large scale (see e.g. Kobayashi et al. (2017)).

Figure 2 shows the perpendicular fluxes of energy and GCH, averaged over the
same time interval as the spectra displayed in figure 1. One observes that the
energy flux is quasi-constant, even if the temporal fluctuations are
comparable to the mean value (not shown) at scales larger than the transition
one. Differently, the GCH flux undergoes a conspicuous decay at the transition
range, a consequence of the helicity barrier whose existence was pointed out
by Meyrand et al. (2021) in similar simulations (although in the asymptotic
regime $\beta_{e}\to 0$), performed in the large $\chi$ regime where the
energy flux to small transverse scales is also significantly inhibited by the
barrier.

Because of numerical constraints, run $R_{1}$ was performed in three steps,
associated with increased resolutions ($N=420^{3}$, $560^{3}$ and $784^{3}$
mesh points) and decreased (equal) parallel and perpendicular viscosities.
This procedure permitted us to analyze the sensitivity of the parallel and
perpendicular dissipations of energy (figure 3 left) and GCH (figure 3 right).
We observe that the perpendicular dissipations are essentially not affected by
the transition from one step to the following one, as expected in a usual
turbulence cascade: reduction of the viscosity does not affect the dissipation
and only shifts the dissipation range to smaller scales. In contrast, we
observe that this is not the case for the parallel dissipation of energy and
even more of GCH, which both decrease with the viscosity, confirming the
absence of a standard turbulent cascade in the direction of the ambient field,
at least in the regime of small or moderate nonlinear parameters considered in
the present study.

Figure 4: Parallel (left) and perpendicular (right) dissipation of GCH in the
$(k_{\perp},k_{z})$-plane (logarithmic coordinates) for run $R_{1}$.

Figure 5: Parallel dissipation of the energies of the $+$ (left) and $-$
(right) waves in the $(k_{\perp},k_{z})$-plane for run $R_{1}$.

Figure 6: Left panel: parallel (red) and perpendicular (blue) dissipations of
the energy of the $+$ wave, and parallel (green) and perpendicular (magenta)
dissipations of the energy of the $-$ wave, integrated on $k_{z}$, versus
$k_{\perp}$ for run $R_{1}$. Right panel: parallel (red) and perpendicular
(blue) dissipations of GCH, integrated on $k_{z}$, versus $k_{\perp}$, for the
same run.

Figure 4 displays two-dimensional color plots of the time-averaged parallel
(left panel) and perpendicular (right panel) GCH dissipations, defined by
$D_{z}^{C}(k_{\perp},k_{z})=\nu_{z}k_{z}^{8}E_{C}(k_{\perp},k_{z})$ and
$D_{\perp}^{C}(k_{\perp},k_{z})=\nu_{\perp}k_{\perp}^{8}E_{C}(k_{\perp},k_{z})$
respectively, where both the viscosities $\nu_{\perp}$ and $\nu_{z}$ are here
equal to $3\times 10^{-10}$ (and where $E_{C}(k_{\perp},k_{z})$ stands for the
GCH power spectrum, after angle averaging in the transverse plane). In
addition to the usual small-scale dissipation in both parallel and
perpendicular directions, one clearly observes the existence of a significant
parallel dissipation of GCH, which is larger near a transverse wavenumber
close to $k_{\perp}=1$ than at the smallest scales.

These observations can be understood by considering the contributions of the
$+$ and $-$ waves, separately. Figure 5 shows the two-dimensional color plots
of the time-averaged parallel dissipations
$D_{z}^{+}(k_{\perp},k_{z})=\nu_{z}k_{z}^{8}E^{+}(k_{\perp},k_{z})$ and
$D_{z}^{-}(k_{\perp},k_{z})=\nu_{z}k_{z}^{8}E^{-}(k_{\perp},k_{z})$ of the
energies associated with the forward and backward propagating waves. A
dissipation spot is visible near $k_{\perp}=1$ in the case of the $+$waves
only. In contrast, comparable dissipations of the forward and backward
propagating waves take place at small scales, both in the parallel and
perpendicular directions. Indeed, the (weak) perpendicular transfer
$\epsilon^{-}$ of $E^{-}$ is almost not affected by the GCH barrier, while
$\epsilon^{+}$ is reduced past the barrier to a value close to $\epsilon^{-}$
since at these scales imbalance has almost disappeared (see also Meyrand et
al. (2021), in a case where the helicity barrier affects the energy transfer
more significantly than in the present simulation). In order to reach a
stationary state, the transfer in the $z$-direction and the resulting parallel
dissipation thus have to become sizeable near the ion scale for the $+$ waves.
At small scale, the dissipation of the two types of waves are comparable and
thus almost cancel out in the GCH dissipation, making the parallel dissipation
of GCH to be dominant at transverse wavenumbers close to $k_{\perp}=1$, as
observed in figure 4.

A synthetic view of the dissipations in terms of the transverse wavenumber is
obtained when integrating on $k_{z}$ the dissipations $D_{z}^{+}$,
$D_{z}^{-}$, $D_{\perp}^{+}$, $D_{\perp}^{-}$ associated with the energy of
the backward and forward propagating waves, and the analogous dissipations
$D_{z}^{C}$, $D_{\perp}^{C}$ for the GCH. The resulting functions
$d_{z}^{+}(k_{\perp})$, $d_{z}^{-}(k_{\perp})$, $d_{\perp}^{+}(k_{\perp})$,
$d_{\perp}^{-}(k_{\perp})$, $d_{z}^{C}(k_{\perp})$, $d_{\perp}^{C}(k_{\perp})$
are displayed in figure 6 (left panel: dissipation of energies; right panel:
dissipations of GCH), under the same time averaging as in figures 4 and 5. The
perpendicular energy dissipation is concentrated at small scales, with
comparable magnitudes for both directions of propagation. Differently, for the
parallel dissipation, while the energy of the backward propagating waves is
dissipated at small scales only, the dissipation of the forward propagating
waves, and thus the GCH, also displays a maximum near the transition zone.

It turns out that the effect of the helicity barrier is enhanced when the
large-scale imbalance is increased, resulting in a steeper transition zone, in
agreement with solar wind observations reported in Huang et al. (2021). This
is exemplified in figure 7 (left) which displays the time-averaged spectrum
$E^{+}(k_{\perp})$ for three runs with different initial imbalance : $I=10$
(run $R_{3}$), $I=100$ (run $R_{2}$), $I=1400$ (run $R_{8}$). Time averages
are performed in the intervals [20000, 25000], [30000, 40000] and [7064,
7096], respectively. Compared to run $R_{1}$, they have a similar value of
$\chi$ but are integrated in a box of size $L=5.5\times 2\upi$, with
resolutions of $280^{3}$ (for $I$ = 10, 100) and $720^{3}$ (for $I$ = 1400)
mesh points. When considering the ratio $E^{+}(k_{\perp})/E^{-}(k_{\perp})$
(figure 7, right), we observe that it becomes close to unity for
$k_{\perp}\gtrsim 2$, whatever the value of the imbalance at large scales.

Figure 7: Influence of the imbalance on the perpendicular $E^{+}(k_{\perp})$
spectrum (left panel) and on the ratio $E^{+}(k_{\perp})/E^{-}(k_{\perp})$
(right panel): $I=10$ (run $R_{3}$), $I=100$ (run $R_{2}$), $I=1400$ (run
$R_{8}$).

###  3.2 Temporal evolution of the dissipations and of the ideal invariants

While at zero imbalance the energy quickly reaches its saturation value (with
nevertheless significant temporal oscillations), when the imbalance is
increased, the time needed to reach saturation can become very large. This is
shown in figure 8 which also indicates that, for a small initial value of
$\chi$, the typical saturation time increases with $I$ (runs $R_{2}$, $R_{3}$.
$R_{4}$).111In order to perform simulations significantly longer than run
$R_{1}$, we resorted to reduce the numerical resolution to $280^{3}$ mesh
points and to use a smaller computational box, of size $L=5.5\times 2\upi$.

Interestingly, for run $R_{2}$ with an energy imbalance $I=100$, the injection
imbalance, given by the GCH to energy injection ratio, is relatively small,
with a value close to $0.1$ at saturation, as seen in figure 9 (left). Here,
the injection rates are given by the slope of the cumulative injections,
estimated as the sums of the energy or the GCH and of their respective
cumulative dissipations. Figure 9 (right) displays the energy and GCH fluxes
for run $R_{2}$ in the steady state. They look qualitatively similar to those
of run $R_{1}$, and in particular confirm the persistence of the GCH barrier
in the saturated regime, in spite of the small value of the injection
imbalance. In contrast, an energy barrier is hardly visible in this regime.

Furthermore, as illustrated in figure 10, the ratio of perpendicular to
parallel energy and GCH dissipations at saturation decreases as $I$ increases.
When saturation is reached with a small value of $\chi$, perpendicular
transfer and dissipation still strongly dominate the parallel ones. This is no
longer the case when $I$ is large, as illustrated with the run $R_{5}$. Note
that in this simulation, the initial value of $\chi$ was increased in order to
shorten the time needed to reach saturation. For such a large value of $\chi$,
the parallel dissipation rapidly overcomes the perpendicular one, as described
in Meyrand et al. (2021). It is to be noted that this situation is not
specific of imbalanced turbulence but also occurs for balanced simulations at
large $\chi$ (not shown), indicating that the breakdown of the asymptotics is
not limited to imbalanced runs. Moreover, it is conspicuous that the ratio of
the perpendicular to parallel dissipation is always smaller for the GCH than
for the energy, pointing out that the barrier acts dominantly on the GCH.

As the nonlinear parameter $\chi$ is increased, the transition zone is less
extended and displays a shallower slope (not shown), except possibly when the
imbalance is also increased, as in the simulations presented in Meyrand et al.
(2021).

Figure 8: Time evolution (averaged over 3 Alfvén crossing times) of the energy
(left panel) and GCH (right panel), normalized to their initial values, for
runs $R_{2}$, $R_{3}$ and $R_{4}$, with similar values of the nonlinear
parameter $\chi$ (whose time evolution is displayed in the insert within the
right panel) and different imbalances: I=100 (red), I=10 (green), I=1
(orange), respectively; a fourth run ($R_{5}$) with I=100 and a larger $\chi$
is also displayed in magenta.

Figure 9: Left: Cumulative injections of energy (solid line) and GCH (dashed
line). Right: energy (solid line) and GCH (dashed line) fluxes of Run 2,
averaged on the time interval $[30000,40000]$.

Figure 10: Time evolution (averaged over 3 Alfvén crossing times) of the ratio
between perpendicular and parallel dissipations of energy (left) and GCH
(right), for the simulations shown in figure 8 (same color code).

###  3.3 Transition zone with negligible parallel dissipation

The transition zone not only exists in regimes where $\nu_{z}=\nu_{\perp}$ (as
in all the runs presented in the previous sections), but persists when the
parallel viscosity $\nu_{z}$ is decreased while keeping $\nu_{\perp}$
unchanged. Figure 11 (left) shows the transverse energy spectra
$E^{\pm}(k_{\perp})$ when $\nu_{z}$ is reduced by a factor $40$ relatively to
its value in the run $R_{2}$. Parallel viscosity can in fact be taken
arbitrary small and even zero without significantly affecting the accuracy of
the numerical simulation during a significant integration time. As seen on
figure 11 (right), when reducing $\nu_{z}$ the perpendicular dissipation
remains unchanged, while the parallel one (which is always subdominant)
decreases further and the contribution of the spectral range close to
$k_{\perp}=1$ still dominates that of the small transverse scales. Regimes
with small parallel dissipation are in fact relevant for the solar wind where
the dynamics being quasi-transverse (Sahraoui et al., 2010), the parallel
transfer, and consequently the parallel dissipation at the smallest scales,
are much weaker than the transverse ones. The above results suggest that the
existence of the transition zone is not a consequence of the parallel
dissipation, at least for small values of the nonlinear parameter. A possible
origin of the transition in this case could be the interactions between co-
propagating waves that take place at the sub-ion scales, as suggested by
Voitenko & De Keyser (2016) and Gogoberidze & Voitenko (2020). This issue is
addressed in the next section, using a semi-phenomenological model.

Figure 11: Left panel: $E^{\pm}$ energy spectra, averaged over the time
interval [15000, 20000], for run $R_{6}$ with a parallel viscosity smaller by
a factor 40 than in run $R_{2}$. Right panel: perpendicular (solid lines) and
parallel (dotted lines) energy dissipations (averaged over 3 Alfvén crossing
times) for simulations performed with different parallel viscosities, all the
other parameters being kept unchanged (red: run $R_{2}$; blue: run $R_{6}$;
cyan: run $R_{7}$).

##  4 Modeling of the transition zone

###  4.1 A phenomenological model including co-propagating wave interactions

Figure 12: Spectra $E^{\pm}((k_{\perp})$ and local slopes predicted by the
phenomenological model for $d=3$ (top) and $d=0$ (bottom), with $C=15$,
$\epsilon=1$, $\eta(k_{\perp})=0.82/v_{ph}^{1.1}(k_{\perp})$.

We here concentrate on a diffusion model in spectral space which was first
derived in Passot & Sulem (2019), and then reproduced in Miloshevich et al.
(2020) using a heuristic approach based on the conservation of the energy and
and of the GCH, together with a phenomenological estimate of the transfer
times. This model does not capture the interactions between co-propagating
waves, not only in the MHD range, where they do not exist, but also in the
sub-ion range where co-propagative KAWs can actually interact. In Passot &
Sulem (2019), following Voitenko & De Keyser (2016), we enriched the model by
retaining interactions between triad wavenumbers which are comparable but not
necessarily asymptotically close. We introduce the nonlinear characteristic
frequencies for the mode ${\boldsymbol{k}}$ propagating in the $\pm$
direction,
$\gamma_{\boldsymbol{k}}^{\pm(\uparrow\uparrow)}=V_{\boldsymbol{k}}^{(\uparrow\uparrow)}\sqrt{k_{\perp}^{3}{E}^{\pm}(k_{\perp})}$
and
$\gamma_{\boldsymbol{k}}^{\pm(\uparrow\downarrow)}=V_{\boldsymbol{k}}^{(\uparrow\downarrow)}\sqrt{k_{\perp}^{3}{E}^{\mp}(k_{\perp})}$
(where $V_{\boldsymbol{k}}^{(\uparrow\uparrow)}$ and
$V_{\boldsymbol{k}}^{(\uparrow\downarrow)}$ estimate the strength of the
corresponding interactions), associated with the interactions between co-
propagating $(\uparrow\uparrow)$ or counter-propagating $(\uparrow\downarrow)$
waves, respectively. Restricting the discussion to the case where both waves
undergo strong nonlinear interactions, we can define global inverse transfer
times as
$(\tau_{tr,G}^{\pm})^{-1}=\gamma_{\boldsymbol{k}}^{\pm(\uparrow\uparrow)}+\gamma_{\boldsymbol{k}}^{\pm(\uparrow\downarrow)}$.
The ratio
$\alpha(k_{\perp})=V_{\boldsymbol{k}}^{(\uparrow\uparrow)}/V_{\boldsymbol{k}}^{(\uparrow\downarrow)}$
should be an increasing function of $k_{\perp}$, that is essentially zero in
the MHD range and saturates to a finite value in the dispersive range.
Identifying the inverse transfer times
$\gamma_{\boldsymbol{k}}^{\pm(\uparrow\downarrow)}$ associated with the
contra-propagating waves with the nonlinear frequencies
$\sqrt{k_{\perp}^{3}v_{ph}^{2}{E}^{\mp}(k_{\perp})}$, the new inverse transfer
times then rewrite

$(\tau_{tr,G}^{\pm})^{-1}\approx
k^{3/2}_{\perp}v_{ph}\left(({E}^{\mp})^{1/2}+\alpha(k_{\perp})({E}^{\pm})^{1/2}\right).$
(9)

Defining $u^{\pm}={E}^{\pm}/k_{\perp}$, the equations governing the evolution
of ${E}^{\pm}(k_{\perp})$ is rewritten, in the stationary regime,

$\frac{d}{dk_{\perp}}u^{\pm}(k_{\perp})=-C\frac{\varepsilon\pm\eta
v_{ph}}{k^{5}_{\perp}v_{ph}}\frac{1}{u^{\mp}(k_{\perp})^{1/2}+\alpha(k_{\perp})u^{\pm}(k_{\perp})^{1/2}},\\\
$ (10)

where the energy and GCH fluxes $\epsilon$ and $\eta$ are a priori functions
of $k_{\perp}$ and $C$ is a numerical constant. This equation is similar to
Eq.(7.34) of Passot & Sulem (2019) when $\nu=0$, except that the weighting
factor of the co-propagating wave interactions is here taken to be $\alpha(k)$
instead of its square.

The function $\alpha$, which is taken as
$\alpha(k_{\perp})=3((1+ck_{\perp}^{2})^{1/2}-1)/((1+ck_{\perp}^{2})^{1/2}+1)$
in Voitenko & De Keyser (2016), with $c$ denoting a slowly varying function of
$k_{\perp}$, can also be written

$\alpha(k_{\perp})=d\left(\frac{v_{ph}(k_{\perp})-s}{v_{ph}(k_{\perp})+s}\right),$
(11)

where $s=\sqrt{2/\beta_{e}}$ is the Alfvén velocity (in sound speed units) and
$d$ a free constant parameter. An important remark, related to the existence
of the helicity barrier, is that Eq. (10) does not have a satisfactory
solution when $\eta\neq 0$. The linear growth of the phase velocity $v_{ph}$
at small scales implies that for $k_{\perp}$ large enough, depending on the
sign of $\eta$, one of the spectra ${E}^{+}$ or ${E}^{-}$ increases with
$k_{\perp}$.

###  4.2 Numerical integration of the phenomenological model

As discussed in Passot & Sulem (2019), Eq. (10) can be considered either as an
initial value problem when the energy and generalized cross-helicity fluxes
$\varepsilon$ and $\eta$ are given and the spectra ${E}^{\pm}$ prescribed at a
wavenumber $k_{\perp}=k_{0}$ or, alternatively, as a nonlinear eigenvalue
problem for $\varepsilon$ and $\eta$, when the spectra are specified at
$k_{0}$ and prescribed to decrease to zero at infinity. Numerical integration
of these equations is best performed as an initial value problem where $k_{0}$
is taken close to the largest wavenumber of the spectral domain. Indeed, the
eigenvalue problem (where ideally $k_{0}$ is chosen at the lower end of the
spectral domain) requires a shooting method and turns out to be extremely
difficult to solve numerically, due to its instability. Indeed, under the
effect of any small perturbation, the small-scale spectrum takes the form of
an absolute equilibrium or goes to zero at a finite value of $k_{\perp}$,
depending on the boundary conditions.

In the integrations presented below, we choose $\beta=2$, $\tau=1$,
$k_{0}=200$ and small enough initial conditions such as to obtain a solution
where the spectra tend to zero for a value of $k_{\perp}$ slightly in excess
of $k_{0}$. As suggested by the numerical simulations of the gyrofluid
equations, in the conditions where $\chi$ is not too large, $\varepsilon$ can
be taken constant. We used $\varepsilon=1$ (close to the almost constant
energy flux observed in the gyrofluid simulations) and $C=15$. Differently,
the gyrofluid simulations show that the flux of GCH undergoes a significant
decrease at the barrier 222As discussed in Passot & Sulem (2019) and
Miloshevich et al. (2020), a constant GCH flux would lead to unrealistically
large amplitudes and large imbalance at large scale, a consequence of the fact
that transfer in the z-direction is not retained by the diffusion model. To
comply with this observation, we were led to prescribe that $\eta$ decreases
faster than $1/v_{ph}(k_{\perp})$. For simplicity, we here chose to take
$\eta(k_{\perp})=\eta_{0}/v_{ph}^{\gamma}(k_{\perp})$. The numerical factor
$\eta_{0}=0.82$ ensures that the imbalance at large scales tends to a value
close to $100$, as in most of the numerical simulations presented in Section
3, while the exponent $\gamma=1.1$ leads to a small-scale imbalance which is
also compatible with the simulations. A larger value of $\gamma$ would produce
a smaller imbalance at small scales. Note that decaying transfer rates also
arise in nonlinear diffusion models of turbulence retaining Landau damping
(Howes et al., 2011).

In figure 12, we display the ${E}^{\pm}$ spectra (left) and their local
spectral exponents (right) for the cases $d=3$ (top) and $d=0$ (bottom), the
other parameters being taken as mentioned above. In the absence of co-
propagating wave interactions ($d=0$), the ${E}^{+}$ spectrum exponent
undergoes a transition, at $k_{\perp}$ close to unity, between the Kolmogorov
exponent $-5/3$ and a value slightly steeper than the theoretical value
$-7/3$, as expected in the presence of imbalance (Passot & Sulem, 2019). When
$d=3$ (a value compatible with kinetic theory, used by Voitenko & De Keyser
(2016)), a clear overshoot of the spectral index is observed, associated with
the existence of a transition zone which starts for $k_{\perp}$ slightly
smaller than unity, the steeper slope being equal to -3.28 at
$k_{\perp}=1.18$. Near this wavenumber, the ${E}^{-}$ spectrum displays a weak
flattening.

Several remarks can be made concerning the influence of the various parameters
of the model, which, at this stage, remains only qualitative. First, the slope
of the transition zone is steeper (shallower) as $d$ in increased
(respectively decreased). This further indicates the role of the co-
propagating wave interactions in the generation of this transition zone.

At larger imbalance (e.g. when taking $\eta_{0}=0.92$, which leads to a large-
scale imbalance close to $I=1000$), the transition zone is slightly steeper,
with a slope reaching $3.52$ at $k_{\perp}=1.11$.

In the absence of co-propagating wave interactions, and when the variation of
$\eta$ is such that $\eta v_{ph}(k_{\perp})$ is constant, the imbalance
remains independent of $k_{\perp}$. This is no longer the case when $\eta$
decreases faster than $1/v_{ph}(k_{\perp})$. As $\gamma$ is increased from the
value $1.1$ chosen in the figure, the imbalance at small scale is strongly
reduced but the transition zone takes a different form. While the slope of the
${E}^{+}$ spectrum is not as steep as for smaller values of $\gamma$, the
${E}^{-}$ spectrum becomes shallower. The case where $\eta$ decreases faster
can be associated with a stronger transfer and/or dissipation in the
$z$-direction and thus a larger value of $\chi$. This effect is seen in some
simulations but a systematic study has not been performed.

##  5 Conclusion

Using a two-field gyrofluid model retaining AWs only, we have observed that,
for $\beta_{e}$ of order unity or larger and a weak or moderate energy
imbalance, a transition zone can exist in the absence of a significant
parallel energy dissipation, in contrast with the mechanism based on the
dissipation induced by the ion-cyclotron resonance, recently proposed by
Squire et al. (2021). This transition zone is found to be steeper at larger
imbalance, in line with Parker Solar Probe observations (Huang et al., 2021).
In these simulations, an helicity barrier also forms, associated with a decay
of the perpendicular CGH flux at the ion scales, while the corresponding
energy flux is almost scale independent. Such simulations are consistent with
the possible effect of the interactions between co-propagating AWs. The
efficiency of this mechanism is supported by the analysis of a
phenomenological model based on diffusion equations in the spectral space, in
the spirit of the Leith’s model (Leith, 1967) for hydrodynamic turbulence. It
turns out that, in the presence of co-propagative waves, and only in this
case, the energy barrier, modeled by a decay of the GCH flux with the
transverse wavenumber, leads to a spectral transition zone. As a consequence
of weak or negligible parallel dissipation, the injected energy accumulates
and a relatively long time is needed for turbulence to reach a steady state.
As a consequence, in the context of a more realistic model, one may suspect
that before ion-cyclotron wave can be excited, other dissipative effects can
come into play, and among them, Landau damping or stochastic heating. Kinetic
simulations for values of the beta parameters of order unity or larger are
necessary to address these issues.

## Acknowledgments

This work was granted access to the HPC resources of CINES/IDRIS under the
allocation A0090407042. Part of the computations have also been done on the
“Mesocentre SIGAMM” machine, hosted by Observatoire de la Côte d’Azur.

## Declaration of Interests

The authors report no conflict of interest.

## References

  * Bowen et al. (2021) Bowen, T. A., Squire, J., Bale, S. D., Chandran, B., Duan, D., Klein, K. G., Larson, D., Mallet, A., McManus, M. D., Meyrand, R., Verniero, J. L. & Woodham, L. D. 2021 The in situ signature of cyclotron resonant heating , arXiv: 2111.05400v1.
  * Brizard (1992) Brizard, A. 1992 Nonlinear gyrofluid description of turbulent magnetized plasmas. Phys. Fluids B 4, 1213–1228.
  * Bruno & Carbone (2013) Bruno, R. & Carbone, V. 2013 The solar wind as a turbulence laboratory. Living Rev. Solar Phys. 10, 2.
  * Bruno & Carbone (2016) Bruno, R. & Carbone, V. 2016 Turbulence in the Solar Wind, Lectures Notes in Physics, vol. 928. Springer.
  * Chandran & Perez (2019) Chandran, B. D. G. & Perez, J. C. 2019 Reflection-driven magnetohydrodynamic turbulence in the solar atmosphere and solar wind. J. Plasma Phys. 85, 905850409.
  * Chen et al. (2020) Chen, C. H. K., Bale, S. D., Bonnell, J. W., Borovikov, D., Bowen, T. A., Burgess, D., Case, A. W., Chandran, B. D. G., de Wit, T. Dudok, Goetz, K., Harvey, P. R., Kasper, J. C., Klein, K. G., Korreck, K. E., Larson, D., Livi, R., MacDowall, R. J., Malaspina, D. M., Mallet, A., McManus, M. D., Moncuquet, M., Pulupa, M., Stevens, M. L. & Whittlesey, P. 2020 The evolution and role of solar wind turbulence in the inner heliosphere. Astrophys. J. Supp. 246 (2), 53.
  * Gogoberidze & Voitenko (2020) Gogoberidze, G. & Voitenko, Y. M. 2020 Model of imbalanced kinetic Alfvén turbulence with energy exchange between dominant and subdominant components. MNRAS 497, 3472–3476.
  * Howes et al. (2011) Howes, G. G., Tenbarge, J. M. & Dorland, W. 2011 A weakened cascade model for turbulence in astrophysical plasmas. Phys. Plasmas 18 (10), 102305–102305.
  * Huang et al. (2021) Huang, S. Y., Sahraoui, F., Andrés, N., Hadid, L. Z., Yuan, Z. G., He, J. S., Zhao, J. S., Galtier, S., Zhang, J., Deng, X. H., Jiang, K., Yu, L., Xu, S. B., Xiong, Q. Y., Wei, Y. Y., Dudok de Wit, T., Bale, S. D. & Kasper, J. C. 2021 The ion transition range of solar wind turbulence in the inner heliosphere: Parker solar probe observations. Astrophys. J. Lett. 909 (1), L7.
  * Kobayashi et al. (2017) Kobayashi, S., Sahraoui, F., Passot, T., Laveder, D., Sulem, P. L., Huang, S. Y., Henri, P. & Smets, R. 2017 Three-dimensional simulations and spacecraft observations of sub-ion scale turbulence in the solar wind: Influence of Landau damping. Astrophys. J. 839 (2), 122.
  * Leith (1967) Leith, C. E. 1967 Diffusion approximation to inertial energy transfer in isotropic turbulence. Phys. Fluids 10, 1409–1416.
  * Lucek & Balogh (1998) Lucek, E. A. & Balogh, A 1998 The identification and characterization of Alfvénic fluctuations in Ulysses data at midlatitudes. Astrophys. J. 507, 984–990.
  * Malara et al. (2001) Malara, F., Primavera, L. & Veltri, P. 2001 Nonlinear evolution of the parametric instability: numerical predictions versus observations in the heliosphere. Nonlinear Proc. Geophys. 8 (3), 159–166.
  * Marsch & Tu (1990) Marsch, E. & Tu, C.-Y. 1990 On the radial evolution of MHD turbulence in the inner heliosphere. J. Geophys. Res. 95 (A9), 8211–8229.
  * Meyrand et al. (2021) Meyrand, R., Squire, J., Schekochihin, A. A. & Dorland, W. 2021 On the violation of the zeroth law of turbulence in space plasmas. J. Plasma Phys. 87, 535870301\.
  * Miloshevich et al. (2020) Miloshevich, G., Laveder, D., Passot, T. & Sulem, P. L. 2020 Modeling imbalanced collisionless Alfvén wave turbulence with nonlinear diffusion equations. Astrophys. J. Lett. 888, L7.
  * Miloshevich et al. (2021) Miloshevich, G., Laveder, D., Passot, T. & Sulem, P. L. 2021 Inverse cascade and magnetic vortices in kinetic Alfvén-wave turbulence. J. Plasma Phys. 87, 905870201\.
  * Passot et al. (2018) Passot, T., Sulem, P.L. & Tassi, E. 2018 Gyrofluid modeling and phenomenology of low-$\beta_{e}$ Alfvén wave turbulence. Phys. Plasmas 25, 042107.
  * Passot & Sulem (2019) Passot, T. & Sulem, P. L. 2019 Imbalanced kinetic Alfvén wave turbulence: from weak turbulence theory to nonlinear diffusion models for the strong regime. J. Plasma Phys. 85, 905850301\.
  * Perez & Chandran (2013) Perez, J. C. & Chandran, B. D. G. 2013 Direct numerical simulations of reflection-driven, reduced magnetohydrodynamic turbulence from the sun to the Alfvén critical point. Astrophys. J. 776 (2), 124.
  * Roberts et al. (1987) Roberts, D. A., Goldstein, M. L., Klein, L. W. & Matthaeus, W. H. 1987 Origin and evolution of fluctuations in the solar wind: Helios observations and Helios-Voyager comparisons. J. Geophys. Res.: Space Physics 92 (A11), 12023–12035.
  * Sahraoui et al. (2010) Sahraoui, F., Goldstein, M. L., Belmont, G., Canu, P. & Rezeau, L. 2010 Three dimensional anisotropic $k$ spectra of turbulence at subproton scales in the solar wind. Phys. Rev. Lett. 105, 131101.
  * Schekochihin et al. (2009) Schekochihin, A. A., Cowley, S. C., Dorland, W., Hammett, G. W., Howes, G. G., Quataert, E. & Tatsuno, T. 2009 Astrophysical gyrokinetics: kinetic and fluid turbulent cascades in magnetized weakly collisional plasmas. Astrophys. J. Suppl. 182, 310–377.
  * Siggia & Patterson (1978) Siggia, E. & Patterson, G.S. 1978 Intermittency effects in a numerical simulation of stationary three-dimensional turbulence. J. FLuid Mech. 86, 567–592.
  * Squire et al. (2021) Squire, Jonathan, Meyrand, Romain, Kunz, Matthew W., Arzamasskiy, Lev, Schekochihin, Alexander A. & Quataert, Eliot 2021 The helicity barrier: how low-frequency turbulence triggers high-frequency heating of the solar wind. arXiv e-prints , arXiv: 2109.03255.
  * Telloni et al. (2015) Telloni, D., Bruno, R. & Trenchi, L. 2015 Radial Evolution of Spectral Characteristics of Magnetic Field Fluctuations at Proton Scales. Astrophys. J. 805 (1), 46\.
  * Tu et al. (1990) Tu, C. Y., Marsch, E. & Rausenbauer, H. 1990 The dependence of MHD turbulence spectra on the inner solar wind stream structure near solar minimum. Geophys. Res. Lett. 17, 283–286.
  * Tu et al. (1989) Tu, C.-Y., Marsch, E. & Thieme, K. M. 1989 Basic properties of solar wind MHD turbulence near 0.3 AU analyzed by means of elsasser variables. J. Geophys. Res. 94 (A9), 11739–11759.
  * Voitenko & De Keyser (2016) Voitenko, Y. & De Keyser, J. 2016 MHD-kinetic transition in imbalanced Alfvénic turbulence. Astrophys. J. Lett. 832, L20.
  * Wicks et al. (2013) Wicks, R. T., Roberts, D. A., Mallet, A., Schekochihin, A. A., Horbury, T. S. & Chen, C. H. K. 2013 Correlations at large scales and the onset of turbulence in the fast solar wind. Astrophys. J. 778, 177.
  * Woodham et al. (2018) Woodham, Lloyd D., Wicks, Robert T., Verscharen, Daniel & Owen, Christopher J. 2018 The role of proton cyclotron resonance as a dissipation mechanism in solar wind turbulence: a statistical study at ion-kinetic scales. Astrophys. J. 856 (1), 49.