11institutetext: Solar-Terrestrial Centre of Excellence — SIDC, Royal
Observatory of Belgium, Avenue Circulaire 3, B-1180, Brussels, Belgium;
11email<EMAIL_ADDRESS>22institutetext: Solar-Terrestrial Centre of
Excellence, Space Physics Division, Royal Belgian Institute for Space
Aeronomy, Brussels, Belgium 33institutetext: Skobeltsyn Institute of Nuclear
Physics, Moscow State University, Leninskie gory, 119991, Moscow, Russia

# Initiation of Alfvénic turbulence by Alfven wave collisions: a numerical
study

S. V. Shestov 11 Y. M. Voitenko 22 A. N. Zhukov, 1133

(28 September 2021)

In the framework of compressional magnetohydrodynamics (MHD), we study
numerically the commonly accepted presumption that the Alfvénic turbulence is
generated by the collisions between counter-propagating Alfvén waves (AWs). In
the conditions typical for the low-beta solar corona and inner solar wind, we
launch in the three-dimensional simulation box two counter-propagating AWs and
analyze polarization and spectral properties of perturbations generated before
and after AW collisions. The observed post-collisional perturbations have
different polarization and smaller cross-field scales than the original waves,
which supports theoretical scenarios with direct turbulent cascades. However,
contrary to theoretical expectations, the spectral transport is strongly
suppressed at the scales satisfying the classic critical balance of
incompressional MHD. Instead, a modified critical balance can be established
by colliding AWs with significantly shorter perpendicular scales. We discuss
consequences of these effects for the turbulence dynamics and turbulent
heating of compressional plasmas. In particular, solar coronal loops can be
heated by the strong turbulent cascade if the characteristic widths of the
loop sub-structures are more than 10 times smaller than the loop width. The
revealed new properties of AW collisions have to be incorporated in the
theoretical models of AW turbulence and related applications.

###### Key Words.:

Magnetohydrodynamics (MHD) - Turbulence - Plasmas - Methods: numerical

##  1 Introduction

Recent studies have revealed that the turbulence in magnetized plasmas is
greatly affected by the Alfvén wave effects. The well-documented example is
the solar-wind turbulence whose nature is essentially Alfvénic and turbulent
fluctuations can be approximately described as Alfvén waves (AWs) (Belcher &
Davis, 1971; Bruno & Carbone, 2013). The standard magnetohydrodynamic (MHD)
description of Alfvénic turbulence in astrophysical and laboratory plasmas is
based on the interaction of oppositely propagating incompressible wave packets
(Iroshnikov, 1963; Kraichnan, 1965).

Following significant previous work on the weak turbulence in incompressible
MHD (Sridhar & Goldreich, 1994; Montgomery & Matthaeus, 1995; Ng &
Bhattacharjee, 1996; Galtier et al., 2000), the more recent work (Howes &
Nielson, 2013) has described the mechanism of turbulent energy transfer via AW
collisions in more detail. The authors showed analytically that two colliding
counter-propagating AWs with wavevectors
$\mathbf{k}_{0}^{-}=k_{\perp}^{-}\mathbf{\hat{y}}+k_{\parallel}\mathbf{\hat{z}}$
and
$\mathbf{k}_{0}^{+}=k_{\perp}^{+}\mathbf{\hat{x}}-k_{\parallel}\mathbf{\hat{z}}$
first produce a specific intermediate wave with
$\mathbf{k}_{2}=k_{\perp}^{+}\mathbf{\hat{x}}+k_{\perp}^{-}\mathbf{\hat{y}}$,
and then its interaction with the initial waves produces the tertiary waves
with wavevectors
$\mathbf{k}_{3}^{-}=k_{\perp}^{+}\mathbf{\hat{x}}+2k_{\perp}^{-}\mathbf{\hat{y}}+k_{\parallel}\mathbf{\hat{z}}$
and
$\mathbf{k}_{3}^{+}=2k_{\perp}^{+}\mathbf{\hat{x}}+k_{\perp}^{+}\mathbf{\hat{y}}-k_{\parallel}\mathbf{\hat{z}}$.
Here $\mathbf{\hat{x}}$, $\mathbf{\hat{y}}$ and $\mathbf{\hat{z}}$ are the
unit Cartesian vectors such that $\mathbf{\hat{z}}$ is parallel to the
background magnetic field $\mathbf{B}_{0}$. These analytical results have been
confirmed by both gyrokinetic simulations in the MHD limit (Nielson et al.,
2013) and experimentally in the laboratory (Drake et al., 2013, 2014, 2016).
Since the energy is transferred to AWs with higher perpendicular wavenumbers,
this process represents an elementary step of the direct turbulent cascade in
which energy is transferred from larger to smaller scales.

Goldreich & Sridhar (1995) introduced the critical balance conjecture and
developed their famous model of strong anisotropic MHD turbulence. The
critical balance assumes that the linear (wave-crossing) and nonlinear (eddy
turnover) times are equal at each scale. Whereas the critical balance remains
a physically reliable hypothesis not strictly derived from basic principles,
it allows for a phenomenological prediction of turbulence properties, in
particular the energy spectrum $\sim k_{\perp}^{-5/3}$ and anisotropy of
turbulent fluctuations. The Goldreich & Sridhar model gave rise to many
important insights in the turbulence nature and resulted in many theoretical,
numerical, and experimental studies (see e.g. Verniero & Howes, 2018; Verniero
et al., 2018; Mallet et al., 2015, and references therein). It is worth noting
that the critical balance conjecture is essentially a statement implying
persistence of linear wave physics in the strongly turbulent plasma.

Despite extended investigations of the critically balanced turbulence, many
actual problems remain open, such as the non-zero cross-helicity effects in
the presence of shear plasma flows (Gogoberidze & Voitenko, 2016), or non-
local effects in AW collisions (Beresnyak & Lazarian, 2008). Also, the plasma
compressibility can introduce surprising effects in the behavior of MHD waves
(Magyar et al., 2019).

Numerical simulations of turbulence are usually done either via numerical
codes for reduced MHD or using analytical frameworks (Beresnyak, 2014, 2015;
Mallet et al., 2015; Perez et al., 2020), pseudo-spectral (Chandran & Perez,
2019) and gyrokinetic (Verniero et al., 2018). Pezzi et al. (2017a, b)
performed simulations using compressible MHD, Hall MHD, and hybrid Vlasov-
Maxwell codes; the 2.5D geometry used in these works did not allow to take
into account nonlinear terms
$\sim\left(\mathbf{v}^{\pm}\cdot\mathbf{\nabla}\right)\mathbf{v}^{\mp}$ and
$\sim\left(\mathbf{b}^{\pm}\cdot\mathbf{\nabla}\right)\mathbf{b}^{\mp}$
($v^{\pm}$ and $b^{\pm}$ are velocity and magnetic fluctuations in $\pm$
waves) for AWs with $\mathbf{k}_{\perp}^{+}\times\mathbf{k}_{\perp}^{-}\neq
0$.

Using compressible MHD model in 3D, we study numerically the commonly accepted
presumption that the AW turbulence is generated by the collisions between
counter-propagating AWs, particularly the wavenumber dependence of the
amplitudes of induced waves. Our simulations reveal that the AW collisions can
occur in two regimes, the first one corresponding to the case of strong
turbulence which follows theoretical explanation, and the second one
corresponding to larger scales which obviously is governed by a different
mechanism.

##  2 Physical and Numerical setup

The simulations were performed in 3D using the numerical code MPI-AMRVAC
(Porth et al., 2014). The code applies the Eulerian approach for solving the
compressible resistive MHD equations:

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})=0,$ (1)
$\frac{\partial(\rho\mathbf{v})}{\partial
t}+\nabla\cdot\left(\mathbf{v}\rho\mathbf{v}-\mathbf{BB}\right)+\nabla
p_{\mathrm{tot}}=0,$ (2) $\frac{\partial\mathbf{B}}{\partial
t}+\nabla\cdot\left(\mathbf{vB}-\mathbf{Bv}\right)=-\nabla\times(\eta\mathbf{J}),$
(3) $\frac{\partial e}{\partial
t}+\nabla\cdot\left(\mathbf{v}e-\mathbf{BB}\cdot\mathbf{v}+\mathbf{v}p_{\mathrm{tot}}\right)=\nabla\cdot\left(\mathbf{B}\times\eta\mathbf{J}\right),$
(4)

where $e$, $\rho$, $\mathbf{v}$, $\mathbf{B}$ are the total energy density,
mass density, velocity, and magnetic field,
$p=(\gamma-1)(e-\rho\mathbf{v}^{2}/2-{B}^{2}/2)$ is the thermal pressure,
$p_{\mathrm{tot}}=p+{B}^{2}/2$ is the total pressure,
$\mathbf{J}=\nabla\times\mathbf{B}$ is the electric current density, $\eta$ is
the electrical resistivity, and $\gamma$ is the ratio of specific heats. The
magnetic field is measured in units for which the magnetic permeability is 1.
Since in this study we are not interested in dissipative processes, we take
$\gamma=5/3$, and $\eta=0$. We used three following normalization constants:
the length $L_{\mathrm{N}}=1$ Mm, the magnetic field $B_{\mathrm{N}}=20$ G,
and the density $\rho_{\mathrm{N}}=1.67\times 10^{-15}$ g cm-3. This
determined normalization for other physical quantities: electron concentration
$n_{\mathrm{N}}=10^{9}$ cm-3, speed
$v_{\mathrm{N}}={B_{\mathrm{N}}}/\sqrt{4\pi\rho_{N}}=1\,380$ km s-1, and time
$t_{\mathrm{N}}=L_{\mathrm{N}}/v_{\mathrm{N}}=0.7246$ s.

The simulations are performed in 3D in Cartesian geometry with a rectangular
numerical box. The background magnetic field $B_{0}=20$ G is directed along
$z$-axis. Equilibrium plasma parameters are taken typical for the solar
coronal base: $n_{e}=10^{9}$ cm-3 ($\rho_{0}=1.67\times 10^{-15}$ g cm-3) and
temperature $T=1$ MK, which determines the plasma beta parameter
$\beta=0.017$. The Alfvén speed in equilibrium plasma is
$v_{A}=B_{0}/\sqrt{4\pi\rho_{0}}=1\,380$ km s-1 or $v_{A}=1$ in normalized
units, and the sound speed is $C_{S}=\sqrt{\gamma\beta/2}v_{A}=0.11v_{A}$.

In order to induce counter-propagating Alfven waves, we set the components of
magnetic field and velocity at the $z$-boundaries of the simulation volume.
The forward wave propagating in $+z$ direction along $\mathbf{B}_{0}$ is
initiated at $z=0$ by the following forcing:

$\displaystyle b_{x}$ $\displaystyle=$ $\displaystyle b\sin\left(\omega
t\right)\sin(k_{\perp}^{-}y);$ (5) $\displaystyle v_{x}$ $\displaystyle=$
$\displaystyle-u\sin\left(\omega t\right)\sin(k_{\perp}^{-}y);$ (6)
$\displaystyle v_{z}$ $\displaystyle=$ $\displaystyle
A_{p}\left[1-\sin\left(2\omega t\right)\right]\sin(k_{\perp}^{-}y),$ (7)

and the backward wave propagating in $-z$ direction is initiated at
$z=z_{\mathrm{max}}$:

$\displaystyle b_{y}$ $\displaystyle=$ $\displaystyle b\sin\left(\omega
t\right)\sin(k_{\perp}^{+}x);$ (8) $\displaystyle v_{y}$ $\displaystyle=$
$\displaystyle u\sin\left(\omega t\right)\sin(k_{\perp}^{+}x);$ (9)
$\displaystyle v_{z}$ $\displaystyle=$ $\displaystyle
A_{p}\left[-1+\sin\left(2\omega t\right)\right]\sin(k_{\perp}^{+}x),$ (10)

where $\omega=k_{z}v_{A}=2\pi\left(\lambda_{z}\right)^{-1}v_{A}$ is the
angular frequency, the parallel wavelength $\lambda_{z}=10$ Mm (always the
same constant in all setups), the initial amplitudes of magnetic field $b$ and
velocity $u=b/\sqrt{4\pi\rho}$ are either 3.33% or 10% of $B_{0}$ and $v_{A}$,
respectively, and $v_{z}\sim
A_{p}=0.25\left(v_{A}^{2}-C_{S}^{2}\right)^{-1}u^{2}v_{A}$ represents the
ponderomotive component of the speed (its order is $10^{-3}$). Boundary
conditions at other boundaries are periodic. Introduction of $A_{p}\neq$ $0$
insures a smooth solution of the MHD equations at the boundaries; its
influence is studied in Sect. 3.1. The physical configuration is shown in Fig.
1.

The described above forcing is applied during 1 period for the forward wave
and 3 periods for the backward wave, which we call the main setups hereafter
(see Table 1 for setup parameters). Beside the main setups, we run several
complementary simulations without backward wave, or with different amplitudes
of counter-propagating waves, or setups with a single period in both waves.

As suggested by the nonlinear term
$\left(\mathbf{z}^{\pm}\cdot\mathbf{\nabla}\right)\mathbf{z}^{\mp}$ in
Elsässer form of MHD equations, in order to allow for effective interactions,
the counter-propagating AWs should have different polarizations. In our
setups, the forward wave is polarized along $x$ and its wavevector
$\mathbf{k}_{\perp}^{-}\parallel\mathbf{\hat{y}}$; the backward wave is
polarized along $y$-axis and $\mathbf{k}_{\perp}^{+}\parallel\mathbf{\hat{x}}$
(see Fig. 1).

Figure 1: Top: physical setup, early phase. The green rectangle denotes the
numerical box with equilibrium plasma, the red and blue areas represent
velocity perturbations (positive and negative) of the $z^{-}$ wave (far
boundary) and $z^{+}$ wave (near boundary). Polarization planes and non-
uniformity directions are annotated. Bottom: longitudinal sketches, late phase
after AWs collision: main setups with 1 period in $z^{-}$ and 3 periods in
$z^{+}$ wave (left), and complimentary setups with 1 periods in both waves
(right). The grey areas denote cross-sections taken for further analysis.
Table 1: Parameters of the numerical setups.  Main setups:
---
High-resolution
Numerical box |  $256\times 256\times 512$ pixels
$L_{z}$ | 50 Mm
$L_{x}$, $L_{y}$ |  equal to $\lambda_{\perp}$
Number of periods |  $z^{-}$ – 1 period
|  $z^{+}$ – 3 periods
$u$ |  0.1 (same for $z^{-}$ and $z^{+}$)
$\lambda_{z}$ (or $\lambda_{\parallel}$) | 10 Mm
$\lambda_{\perp}$ | from 0.4 to 25.0 (10 configurations)
$k_{\perp}$ | from 15.7 to 0.25
$k_{\perp}/k_{\parallel}$ | from 25.0 to 0.4
Low-resolution
Numerical box |  $128\times 128\times 256$ pixels
$L_{z}$ | 50 Mm
$L_{x}$, $L_{y}$ |  equal to $\lambda_{\perp}$
$u$ |  0.033 (same for $z^{-}$ and $z^{+}$)
$\lambda_{z}$ (or $\lambda_{\parallel}$) | 10 Mm
$\lambda_{\perp}$ | from 0.16 to 25.0 (12 configurations)
$k_{\perp}$ | from 39.3 to 0.25
$k_{\perp}/k_{\parallel}$ | from 62.5 to 0.4
Non-zero cross-helicity:
Numerical box |  $256\times 256\times 364$ pixels
$L_{z}$ | 30 Mm
$L_{x}$, $L_{y}$ |  equal to $\lambda_{\perp}$
Number of periods |  $z^{-}$, $z^{+}$ – 1 period
$u^{-}$ | 0.1
$u^{+}$ | 0.03
$\lambda_{z}$ (or $\lambda_{\parallel}$) | 10 Mm
$\lambda_{\perp}$, $k_{\perp}$, $k_{\perp}/k_{\parallel}$ | same as in high-resolution main setups
Perpendicular and longitudinal structure:
Numerical grid |  $256\times 256\times 364$ pixels
$L_{z}$ | 30 Mm
$L_{x}$, $L_{y}$ |  equal to $\lambda_{\perp}$
Number of periods |  $z^{-}$, $z^{+}$ – 1 period
$u$ |  0.1 (same for $z^{-}$ and $z^{+}$)
$\lambda_{z}$ (or $\lambda_{\parallel}$) | 10 Mm
$\lambda_{\perp}$, $k_{\perp}$, $k_{\perp}/k_{\parallel}$ | same as in high-resolution main setups

The numerical box has physical $z$-length $L_{z}$ either 50 Mm (main setups)
or 30 Mm (complementary setups). The sizes along $x$ and $y$ are set equal to
the perpendicular wavelength $\lambda_{\perp}$ (hence change from setup to
setup). The numerical box for the main setups has either $256\times 256\times
512$ pixels (high-resolution) or $128\times 128\times 256$ pixels (low-
resolution). We have verified that the decrease of numerical resolution does
affect the results: the waves start to decay during their propagation and the
wave profiles get distorted. However, this effect is small even for the case
of low-resolution setups. In complimentary setups the numerical box always has
$256\times 256\times 384$ pixels, thus its spatial resolution coincides with
that of the high-resolution setups. We compared various numerical schemes and
parameters of MPI-AMRVAC and chose the best settings (`powel` scheme for the
$\mathbf{\nabla\cdot B}$ corrector, high-resolution numerical box etc.). We
also pay special attention to distinguish the physical phenomena from
numerical artifacts.

##  3 Results

###  3.1 Nonlinear effects in a single AW

First we verify the effect of nonlinear self-interaction within a single
Alfvén wave. In Fig. 2 we show the longitudinal (along $z$) profiles of
$v_{x}$, $v_{y}$, and $v_{z}$ of the forward Alfvén wave, initiated via
boundary conditions described in Sect. 2, during its developed phase, but
before the collision with the backward wave. For visualization the quantities
are normalized by the following constants: the mother wave $v_{x}$ by the
initial amplitude $u=0.10$, the horizontal component $v_{y}$ and the
ponderomotive component $v_{z}$ by $A_{p}=2.53\cdot 10^{-3}$.

The amplitude and spatial structure of the ponderomotive component $v_{z}$
perfectly reproduces theoretical predictions: its wavenumbers are two times
larger than in the mother wave and its amplitude varies from 0 to 2
(McLaughlin et al., 2011; Zheng et al., 2016). We also observed a self-
consistent generation of $v_{y}$ that appears only in oblique waves with
$\lambda_{\perp}\neq 0$ ($v_{y}=0$ at $\lambda_{\perp}=0$). Our preliminary
simulations (two-dimensional setups were sufficient there) have shown the
following trend in the variation of $v_{y}$ with varying cross-field
wavelength: the amplitude of $v_{y}$ grows proportionally to
$1/\lambda_{\perp}$ at the larger scales $\lambda_{\perp}>\lambda_{z}$, this
growth slows down at $\lambda_{\perp}\sim\lambda_{z}$, and eventually $v_{y}$
becomes a constant independ on $\lambda_{\perp}$ at smaller scales
$\lambda_{\perp}\ll\lambda_{z}$. The spatial extension of $v_{y}$ in both
parallel and perpendicular directions is two times shorter than of the mother
wave. The amplitude of $v_{y}$ is always smaller than that of $v_{z}$. Similar
perturbations of perpendicular velocity were observed also in torsional waves
(Shestov et al., 2017).

The observed perturbations of $v_{y}$ and $v_{z}$ propagate along the magnetic
field with the Alfvén speed $v_{A}$ and are natural companions of AWs not
caused by the numerical effects or boundary conditions for $v_{z}$. The
perturbations always develop in AWs regardless of the ways how the waves are
initiated – by boundary or initial conditions, with or without boundary
perturbations given by Eqs. 7 and 10. In other words, the observed propagating
wave is the eigenmode of the compressible nonlinear MHD.

Figure 2: Velocities $v_{x}$ (left), $v_{y}$ (middle), and $v_{z}$ (right) in
a single Alfvén wave. The color tables (inlines in top) have different
amplitudes in different panels to reflect the range of velocities

. Initially in the Alfvén wave only $v_{x}$ and $B_{x}$ (not shown) and
$v_{z}$ are driven; the $v_{y}$ component of the velocity is generated self-
consistently due to nonlinear self-interaction within the Alfvén wave.

We thus observe typical characteristics of AWs before they collide.

###  3.2 AWs collision

To study effects of the AW collisions, we let the two counter-propagating
waves to fully propagate through each other, and analyze perpendicular
profiles of $v_{x}$ of the forward-propagating $z^{-}$ wave in its leading
maximum – $x-y$ plane with $z=42.5$ at instant $t=45$, see Fig. 1, bottom left
panel (main setups with $u=0.1$ are used). In Fig. 3 the panels show $v_{x}$
of three different setups with $\lambda_{\perp}=0.5$, $0.8$, and $3.0$ Mm. The
perturbations of the wave profiles depend on the perpendicular scale: they are
significant for smallest $\lambda_{\perp}=0.5$, moderate for
$\lambda_{\perp}=0.8$, and weak for the largest $\lambda_{\perp}=3.0$. Similar
perturbations are also observed in the $z^{+}$ wave. However, in setups with
only one wave present, such perturbations do not appear, and hence their
development can be attributed to AWs collision.

Figure 3: Comparison of perpendicular profiles of $v_{x}$ in the leading
maximum of the $z^{-}$ wave in setups with $\lambda_{\perp}=0.5$ (left), $0.8$
(middle), and $3.0$ Mm (right). The axes are in Mm and are different in every
panel.

Appearance of such small-scale perturbations propagating with Alfvén velocity
can be treated as generation of new AWs at smaller perpendicular scales
$\lambda_{\perp}^{\prime}<\lambda_{\perp}$.

###  3.3 Dependence on perpendicular scales

In order to distinguish the nonlinearly generated waves from the mother wave,
we further analyze the wave profiles in the perpendicular cross-section of
$z^{-}$ wave. We extract the perturbed velocity $\Delta v_{x}$ by subtracting
the initial harmonic profile of $v_{x}$, $\Delta
v_{x}=v_{x}-A\sin(k_{\perp}^{-}y)$, where the amplitude $A$ is adjusted to
cancel the perturbation in the wave maximum. The results are shown in Fig. 4
for the setups with $\lambda_{\perp}=0.8$ (top panels) and
$\lambda_{\perp}=3.0$ (bottom panels). The left panels show $v_{x}$, middle
$\Delta v_{x}$, and right $v_{y}$. The induced velocities $\Delta v_{x}$ and
$v_{y}$ have amplitudes $\sim(0.02\div 0.03)u$ and are non-uniform in both $x$
and $y$ directions.

To evaluate numerical effects, we made the similar analysis for $z^{-}$ wave
in the absence of $z^{+}$ waves. Here the perturbations $\Delta v_{x}$ are
observed as well; but they have at least factor 10 smaller amplitude and are
uniform along $x$. It means that numerical effects produce significantly
weaker perturbations with different spatial profiles. On the contrary, after
collisions with counter-propagating $z^{+}$ waves, the perturbations co-
propagating with $z^{-}$ waves have both $v_{x}$ and $v_{y}$ components,
larger amplitudes, and profiles non-uniform both along $y$ and $x$, which
cannot be ascribed to numerical effects. Furthermore, the perturbations of
$v_{y}$ generated by the AW collisions can not be attributed solely to the
single AW self-interaction where perturbations of the $v_{y}$ are zero at the
original wave maximum.

The spatial patterns of the induced velocities fall in two distinct groups:
all spatial patterns at $\lambda_{\perp}<\lambda_{\perp}^{tr}$ are similar to
that shown on the top panels in Fig. 4, and all patterns at
$\lambda_{\perp}\geq\lambda_{\perp}^{tr}$ are similar to that shown on the
bottom panels (the transition scale $\lambda_{\perp}^{tr}=3.0$ for $u=0.1$
used in this figure). The perturbations in the former group have a current-
sheet structuring, similar to that reported by Verniero et al. (2018) for the
strong turbulence regime. The perturbations in the second group have symmetric
structure. The same two groups of spatial structures are also observed in the
setups with different amplitudes $u$, but with different transition scales,
such that $\lambda_{\perp}^{tr}$ is larger for smaller $u$ (for example,
$\lambda_{\perp}^{tr}=4.0$ for $u=0.033$).

Figure 4: Perpendicular profiles of velocities in setups with
$\lambda_{\perp}=0.8$ (top) and $\lambda_{\perp}=3$ (bottom). Left: measured
$v_{x}$; middle: difference $\Delta v_{x}$ between the measured $v_{x}$ and a
harmonic function; right: measured $v_{y}$. Perturbations of $\Delta v_{x}$
and $v_{y}$ are produced in result of AWs collision. In each panel the color
table matches the maximum amplitude of the measured quantity.

The dependence of the amplitudes of induced waves on the perpendicular scales
is shown in Fig. 5. The diamonds correspond to $\Delta v_{x}$ and asterisks
correspond to $v_{y}$. For $u=0.1$ the symbols are blue and green, for
$u=0.033$ they are orange and red. Gray and pink regions indicate the
wavenumber ranges where the wave collisions should generate the strong
(critically balanced) turbulence with $k_{\perp}/k_{\parallel}\sim v_{A}/u$
for $u=0.1$ and $u=0.033$, respectively. In both these cases the amplitude
behavior is similar. At largest $\lambda_{\perp}$ the amplitudes of the
induced waves are much smaller than the amplitudes of the original waves and
the resulting AW turbulence should be weak. As $\lambda_{\perp}$ decreases,
the induced amplitudes first increase slowly and reach a maximum. This maximum
is still much smaller than the initial AW amplitude and is reached at
$\lambda_{\perp}=\lambda_{\perp\mathrm{max}}$ that is still much larger than
the perpendicular scale given by the critically balance condition,
$\lambda_{\perp\mathrm{max}}\gg$
$\lambda_{\perp\ast}=\lambda_{\parallel}u/v_{A}$
($\lambda_{\perp\mathrm{max}}$, $\lambda_{\perp\ast}$ and other characteristic
perpendicular scales are shown in Fig. 5). When $\lambda_{\perp}$ decreases
further beyond $\lambda_{\perp\mathrm{max}}$, the induced amplitudes decrease
and reach a minimum at $\lambda_{\perp}=$ $\lambda_{\perp\mathrm{min}}$ that
is still larger than $\lambda_{\perp\ast}$. After this minimum, a strong
increase of induced perturbations occurs in the region where $\lambda_{\perp}$
becomes several times shorter than $\lambda_{\perp\ast}$. Amplitudes of
generated perturbations become there comparable to the amplitudes of initial
waves and such collisions can generate strong turbulence.

While the observed strengthening of the nonlinear interaction with decreasing
$\lambda_{\perp}$ is expected taking into account that the responsible
nonlinear term is $\sim\left(\mathbf{z}\cdot\nabla\right)\mathbf{z\sim}$
$\lambda_{\perp}^{-1}$, the depression observed at $\lambda_{\perp}\gtrsim$
$\lambda_{\perp\ast}$ and its influence on the transition from weak to strong
turbulence need further investigations. At present we can only state that this
depression should result in a shift of the weak-strong turbulence transition
to the perpendicular scales significantly shorter than that prescribed by the
standard critical balance condition.

Figure 5: Dependence of the amplitudes of the induced waves on
$k_{\perp}/k_{\parallel}$ of the original wave for the $u=0.1$ (top) and
$u=0.033$ (bottom) setups. Diamonds and asterisks correspond to the $\Delta
v_{x}$ and $v_{y}$ respectively. Filled areas correspond to the regions of
classic critical balance in incompressible MHD calculated for particular $u$.
In the top panel the values of $\lambda_{\perp\mathrm{max}}$,
$\lambda_{\perp\mathrm{min}}$, $\lambda_{\perp}^{tr}$, and $\lambda_{\perp*}$
are shown.

###  3.4 Influence of several collisions

Since the initiated $z^{-}$ and $z^{+}$ waves contain 1 and 3 periods,
respectively, the $z^{-}$ wave can interact with 3 periods of the counter-
propagating wave, whereas each period of the $z^{+}$ wave can interact with
only one period of $z^{-}$. We thus expect different amplitudes of the induced
perturbation propagating in $z^{-}$ and $z^{+}$ directions. To verify this, we
measure the perturbations accompanying the $z^{+}$ wave using the same
technique as for $z^{-}$ wave (remember that in $z^{+}$ wave the roles of
$v_{x}$ and $v_{y}$ are exchanged). Comparison of the corresponding
perturbations in the $z^{-}$ and $z^{+}$ waves is given in Fig. 6. The top
panel shows the amplitudes of perturbations accompanying $z^{-}$ (blue and
green symbols) and $z^{+}$ (orange and red symbols), the bottom panel shows
the ratio $z^{-}/z^{+}$ of the perturbations with the corresponding
(orthogonal) polarizations. In both panels the diamonds denote perturbations
with the same polarization as in the original waves ($\Delta v_{x}$ in
$z^{-}$, $\Delta v_{y}$ in $z^{+}$), and the asterisks denote the
complimentary polarization.

The behavior of $z^{+}$ perturbations as function of $\lambda_{\perp}$ is
qualitatively similar to that of $z^{-}$ perturbations. At smallest
$\lambda_{\perp}$ the ratio of the $-/+$ perturbations is about 2, then
approaches 3 with the scale increase, then increases significantly at
$\lambda_{\perp}/\lambda_{\parallel}\sim 1$, and finally drops again to 2 at
large perpendicular scales $\lambda_{\perp}/\lambda_{\parallel}>1$. In the
region of (super-)strong turbulence the observed ratio $z^{-}/z^{+}<3$ means
inapplicability of the perturbation theory: already after the first
interaction the wave profiles are distorted significantly and the following
collisions do not add much.

Figure 6: Dependence of the amplitudes of induced waves on the wavenumber
ratio for different number of wave collisions ($z^{-}$ wave collides 3 times
and $z^{+}$ wave collides once). Top: amplitudes of the perturbations. Bottom:
ratio of the amplitudes $z^{-}/z^{+}$, diamonds correspond to original
polarization ($v_{x}$ in $z^{-}$, and $v_{y}$ in $z^{+}$) and asterisks
correspond to perpendicular polarization.

###  3.5 Non-zero cross-helicity case

In this section we analyze the effects of non-zero cross-helicity (imbalance)
when the counter-propagating initial waves have different amplitudes. This
situation is common in the fast solar wind (Tu et al., 1990; Lucek & Balogh,
1998) and also occurs in numerical simulations in local subdomains of the
simulation box (Perez & Boldyrev, 2009).

We run dedicated setups with initial amplitudes $u^{-}=0.1$ in $z^{-}$ wave
and $u^{+}=0.033$ in $z^{+}$ wave, both waves have one period. We compare
measured perturbations with our main setups in Fig. 7. The black symbols
denote $z^{-}$ perturbation and orange and red symbols denote $z^{+}$
perturbations in imbalanced setups, and blue and green symbols denote main
setups ($u=0.1$, 1 period in $z^{-}$ wave and 3 periods in $z^{+}$ wave).

The perturbations observed in imbalanced cases are smaller then in the main
setups. At the same time the perturbations (expressed in initial amplitudes
$u$) in the $z^{-}$ wave are $\sim 3$ times smaller then in the $z^{+}$ wave.

Figure 7: Amplitude of the induced waves for the case of non-zero cross-
helicity: $u^{-}=0.1$, $u^{+}=0.033$, each wave has single period. Black
symbols denote perturbations in the $z^{-}$ wave, orange and red symbols
denote perturbations in the $z^{+}$ wave. Blue and green symbols represent the
main setups ($u=0.1$, 1 and 3 periods respectively).

###  3.6 Perpendicular Fourier spectra

In order to understand the spectral transport generated by the AW collisions,
we analyze the spatial Fourier spectra of the induced waves. The spectra of
the $v_{x}$ and $v_{y}$ velocities at the leading maximum of $z^{-}$ are given
in Fig. 8 for $\lambda_{\perp}=0.8$ in the top row, $\lambda_{\perp}=2.0$ in
the middle row, and $\lambda_{\perp}=3.0$ in the bottom row. On the left
panels, the spectra of $v_{x}$ in a single $z^{-}$ wave are shown, on the
middle and right panels the spectra of $v_{x}$ and $v_{y}$ after the AW
collision are shown. In each panel the $(x,y)$-coordinates represent
corresponding Fourier wavenumbers and the color shows intensity of a given
spectral component. The quasi-logarithmic color scale is normalized to the
intensity of an ideal harmonic function $u\sin(k_{\perp}^{-}y)$. This function
would have only two peaks with spectral coordinates $(0,\pm 1)$ that
correspond to the brightest components in the left and middle panels. In what
follows, we will drop the $\pm$ sign keeping in mind the inherent symmetry.

The higher-wavenumber spectral components $(0,\left|y\right|>1)$ accompanying
the single $z^{-}$ wave without collisions (left panels) are due to numerical
effects. Note the low level of these components and their uniform
distribution. On the contrary, the real spectral components with higher
wavenumbers are generated by the AW collisions (middle and right panels).

The strongest induced components at $\lambda_{\perp}=0.8$ have spectral
coordinates $(1,2)$ corresponding to the perpendicular wavevector
$\mathbf{k_{\perp}}=k_{\perp}^{+}\mathbf{\hat{x}}+2k_{\perp}^{-}\mathbf{\hat{y}}$.
Generation of waves with such wavevectors supports the mechanism proposed by
Howes & Nielson (2013) (see their Fig. 2 explaining appearance of such
“tertiary” waves). This mechanism is summarized in the Introduction.

The same Fourier components $(1,2)$ of $v_{x}$ are also seen in the middle row
Fig. 8 in the case of intermediate scale $\lambda_{\perp}=2$; in addition, the
spectral components of $v_{x}$ with coordinates $(1,1)$ corresponding to
$\mathbf{k_{\perp}}=k_{\perp}^{+}\mathbf{\hat{x}}+k_{\perp}^{-}\mathbf{\hat{y}}$
are significant as well.

The spatial spectra of $v_{x}$ and $v_{y}$ at the largest scale
$\lambda_{\perp}=3$ are qualitatively different: the strongest induced
components have coordinates $(1,1)$ while the others are negligible. These
spectral components might be formed by a different mechanism than in the
$\lambda_{\perp}=0.8$ case.

In general, the spectral dynamics observed in our simulations, i.e. generation
of higher-wavenumber spectral components, supports scenarios with direct
turbulent cascades generated by AW collisions.

Figure 8: Perpendicular Fourier spectra of $v_{x}$ and $v_{y}$ velocities
measured in the leading maximum of the $z^{-}$ wave. Top row:
$\lambda_{\perp}=0.8$; middle row: $\lambda_{\perp}=2$; bottom row:
$\lambda_{\perp}=3$. Left: $v_{x}$ of the setup with the $z^{-}$ wave only;
middle and right denote $v_{x}$ and $v_{y}$ of the setups after AWs collision.
In each panel the $(x,y)$-coordinates represent corresponding Fourier
wavenumbers and the color shows intensity of a given spectral component. The
quasi-logarithmic color table is normalized to the amplitude of a harmonic
wave.

###  3.7 Field-aligned structure of the induced Alfvén waves

Longitudinal behavior of the Fourier components of the induced Alfvén waves is
studied using the following approach: we Fourier-analyze perpendicular cross-
sections at multiple $z$-coordinates, covering the distance of slightly more
then one full wavelength $\lambda_{\parallel}$ along $z$ (see bottom right
sketch in Fig. 1). In Fig. 9 we show longitudinal behaviour of the spectral
components of $v_{x}$ with coordinates $(0,1)$, $(1,1)$, $(1,2)$, $(2,1)$, and
$(2,2)$ with different colors. The mother wave with spectral coordinates
$(0,1)$ is shown with black. The top panel shows the setups with
$\lambda_{\perp}=0.5$, and the bottom panel shows the setup with
$\lambda_{\perp}=2.0$. The intensity of the spectral components is multiplied
by factor 10 in the top panel, and by factor 100 in the bottom panel.

We observe drastically different behavior of the spectral components in
different setups. While we do not see any regularity in the larger-scale
setup, in the setup with $\lambda_{\perp}=0.5$ the growth of $(1,1)$ and
$(1,2)$ components is highly correlated and their parallel scales are somehow
shorter than in initial AWs. In addition, the energy of the induced waves tend
to concentrate near the center of the mother wave.

Figure 9: Longitudinal dependence of amplitudes of spatial Fourier components
of $v_{x}$. Top: setup with $\lambda_{\perp}=0.5$; bottom: setup with
$\lambda_{\perp}=2$ Mm. Different colors correspond to particular spectral
coordinates.

##  4 Discussion and Application

Results of our simulations revealed several new properties of AW collisions in
compressional plasmas, which can affect Alfvénic turbulence and anisotropic
energy deposition in plasma species. The most striking new property is the
modified relation between the parallel and perpendicular scales in the strong
turbulence regime where the energy is efficiently transferred to the smaller
scale during one collision.

The turbulence strength is usually characterized by the nonlinearity parameter
$\chi_{k}\equiv\tau_{k}^{\mathrm{L}}/\tau_{k}^{\mathrm{NL}}=\left(k_{\perp}v_{k}\right)/\left(k_{z}v_{A}\right)$,
where
$\tau_{k}^{\mathrm{NL}}=\lambda_{\perp}/v_{k}=2\pi/\left(k_{\perp}v_{k}\right)$
is the nonlinear mixing time,
$\tau_{k}^{\mathrm{L}}=\lambda_{z}/v_{A}=2\pi/\left(k_{z}v_{A}\right)$ is the
linear (correlation) crossing time of colliding AWs, and $v_{k}$ is the
velocity amplitude of the colliding AWs. Denote by $\delta v_{k}$ the velocity
amplitude of generated waves. When the classic critical balance condition of
incompressible MHD is satisfied,

$\chi_{k}=1,$ (11)

the nonlinear mixing becomes as fast as the linear crossing and the turbulence
is believed to be strong, $\delta v_{k}/v_{k}\sim 1$ (Goldreich & Sridhar,
1995).

However, as follows from our simulations (see e.g. Fig. 5 showing $\delta
v_{k}/v_{k}$ as function of $k_{\perp}/k_{z}$ for two fixed amplitudes,
$u\equiv v_{k}/v_{A}=0.1$ and $0.033$, and $\delta v_{k}=\sqrt{\Delta
v_{x}^{2}+v_{y}^{2}}$), the spectral transport in compressible MHD is
strongly, about one order of magnitude, suppressed at $k_{\perp}/k_{z}$
satisfying Eq. 11. Namely, $\delta v_{k}/v_{k}\ll 1$ at $k_{\perp}/k_{z}=10$
for $u=0.1$ and at $k_{\perp}/k_{z}=30$ for $u=0.033$. At $k_{\perp}/k_{z}$
increasing further, the spectral transport eventually becomes fast and the
turbulence strong, $\delta v_{k}/v_{k}\sim 1$, which happens at
$k_{\perp}/k_{z}$ obeying the modified critical balance condition

$\tilde{\chi}_{k}=\alpha\chi_{k}=1,$ (12)

where $\alpha<1$ is the factor reducing efficiency of the nonlinear mixing (in
other words, the effective nonlinear time increases by the factor $1/\alpha$).
Consequently, the turbulence becomes strong at perpendicular wavenumbers that
are larger than in the classic critically balanced case.

The origin and nature of $\alpha$ need further clarification. Since
$\alpha\neq 1$ arises when the plasma compressibility is taken into account,
it should depend on the relative content of thermal energy, e.g. on the plasma
$\beta$. For parameters adopted in our simulations, $\alpha\approx 0.3$ in the
critically-balanced state where the scale ratio $k_{\perp}/k_{z}$ obeys
$\tilde{\chi}_{k}\sim 1$. Such departure from the classic critical balance
affects dynamics of the strong AW turbulence (see below). In the general case
of arbitrary scales the functional dependence
$\alpha=\alpha\left(\beta,v_{k},k_{\perp}/k_{z}\right)$ is complex; in
particular, $\delta v_{k}/v_{k}$ becomes a decreasing function of
$k_{\perp}/k_{z}$ in some interval (in Fig. 5 it happens at
$k_{\perp}/k_{z}\lesssim 10$), which should greatly affect the weak AW
turbulence. We do not exclude that $\alpha$ may also depend on other
plasma/wave parameters.

Let us consider the Alfvénic turbulence driven by the fluctuating velocity
$v_{k0}$ at the wavenumber ratio $k_{\perp 0}/k_{z0}$ obeying the critical
balance condition $\tilde{\chi}_{k0}=1$, in which case the turbulence is
already strong at the driving scales. The spectral energy flux in the inertial
range is

$\epsilon_{s}=\frac{\rho v_{k}^{2}}{\tau_{k}^{\mathrm{L}}}=\alpha\frac{\rho
v_{k}^{2}}{\tau_{k}^{\mathrm{NL}}}\approx\alpha\frac{\rho
v_{k}^{3}k_{\perp}}{2\pi}=\mathrm{const}\equiv\alpha_{0}\frac{\rho
v_{k0}^{3}k_{\perp 0}}{2\pi},$ (13)

where $\tau_{k}^{\mathrm{L}}=\lambda_{z}/v_{A}$ is the AW collision time and
$\rho\approx n_{0}m_{i}$ is the mass density. Note that the spectral flux
$\epsilon_{s}$ from Eq. 13 is reduced as compared to the incompressible strong
turbulence driven at the same perpendicular scale, but remains the same for
the turbulence driven at the same parallel scale.

Assume that there is a weak dependence
$\alpha=\alpha_{0}\left(v_{k}/v_{k0}\right)^{\delta}$, where $0<\delta<3/4$.
Such dependence is suggested by the following semi-empirical considerations.
As the observed spectra are power-law, the scaling of $\alpha$ with $v_{k}$
should be power law as well. Furthermore, the index $\delta$ of the power-law
dependence should be small positive to reproduce the observed in simulations
mismatch between the classic and real critical balances (which is larger for
larger wave amplitude). Moreover, such positive values of $\delta$ appear to
be compatible with the observed spectral indexes of turbulence in the quasi-
stationary solar wind, which are slightly larger than $-5/3$ (up to $-3/2$).

The kinetic energy spectrum is then flatter than the Kolmogorov one,

$W_{s\perp}\sim v_{k}^{2}/k_{\perp}\sim k_{\perp}^{-5/3+2\delta/9},$ (14)

and its spectral index varies between $-5/3$ and $-3/2$, as is typically
observed in the solar wind turbulence. In the case of $\alpha$ constant along
the critical balance path, $\delta=0$, the spectrum reduces to the classic
Kolmogorov $W_{s\perp}\sim v_{k}^{2}/k_{\perp}\sim k_{\perp}^{-5/3}.$ The
parallel wavenumber spectrum is, as usual, $W_{sz}\sim v_{k}^{2}/k_{z}\sim
k_{z}^{-2}$.

If the turbulence is weak at injection, $\tilde{\chi}_{k0}=\delta
v_{k0}/v_{k0}<1$, the cascade time increases from the strong turbulence value
$\tau_{k}^{\mathrm{TC}}\sim\tau_{k}^{\mathrm{L}}$ to the weak turbulence value
$\tau_{k}^{\mathrm{TC}}\sim\left(\delta
v_{k}/v_{k}\right)^{2}\tau_{k}^{\mathrm{L}}$. The resulting weakly turbulent
energy flux $\epsilon_{w}$ decreases as compared to the strongly turbulent
energy flux (13), $\epsilon_{w}=\tilde{\chi}_{k}\epsilon_{s}$:

$\epsilon_{w}=\frac{\rho v_{k}^{2}}{\tau_{k}^{\mathrm{TC}}}=\frac{\rho
v_{k}^{2}}{\tau_{k}^{\mathrm{L}}}\tilde{\chi}_{k}^{-2}=\mathrm{const}\equiv\frac{\rho
v_{k0}^{2}}{\tau_{k0}^{\mathrm{L}}}\tilde{\chi}_{k0}^{-2}.$ (15)

The weakly turbulent spectrum is problematic to calculate because of a complex
dependence of $\delta v_{k}/v_{k}$ upon $k_{\perp}$ and $v_{k}$ (see Fig. 5),
which is unknown and difficult to guess. At present, we can only note that the
strength $\tilde{\chi}_{k0}$ of the compressional weak turbulence is much
(about one order, as is demonstrated by Fig. 5) smaller than the
incompressional one, $\tilde{\chi}_{k0}\ll\chi_{k0}$, which drastically
decreases the weakly turbulent energy flux.

Although the large-scale MHD AWs do not dissipate directly, the turbulent
cascade transfers their energy to small scales where dissipative effects come
into play heating plasma. MHD Alfvénic turbulence has been employed as the
mechanism for plasma heating in the solar corona and solar wind, both from the
theoretical/modeling perspective (Van Ballegooijen et al., 2011; Verdini, A.
et al., 2012) and based on experimental observations of quiescent (Morton et
al., 2016; De Moortel et al., 2014; Xie et al., 2017) and flaring loops
(Doschek et al., 2014; Kontar et al., 2017). Here we discuss how the new
properties of AW collisions observed in our simulations can affect models of
quasi-steady turbulent plasma heating in coronal loops.

Recently, Xie et al. (2017) analyzed as many as 50 loops in active regions
using observations of Extreme-ultraviolet Imaging Spectrometer (EIS) (Culhane
et al., 2007) on board the Hinode satellite. They observed non-thermal widths
of spectral lines and found corresponding non-thermal velocities in the range
$v_{\mathrm{nt}}=30\div 40$ km s-1, magnetic field in the loop apexes up to 30
G, loop widths $L_{\perp}\sim 2\div 4$ Mm and loop lengths $L_{z}\sim 100$ Mm.
Brooks & Warren (2016) also used spectroscopic data from EIS and evaluated
non-thermal velocities in loops in 15 active regions. The typical values were
somewhat smaller, with typical values $v_{\mathrm{nt}}\sim 20$ km s-1; the
authors however did not provide any other parameters. Furthermore, Gupta et
al. (2019) analyzed non-thermal widths of spectral lines in high coronal loops
(with heights up to $1.4R_{\odot}$) measured by EIS and found the non-thermal
velocities in the range $20\div 30$ km s-1. The above values can be used to
evaluate the turbulent heating of coronal loops.

We assume that there are AW sources at the loop footpoints. These source can
be due to magnetic reconnection and/or photospheric motion (we will not
specify their origin in more details here). The perpendicular AW wavelengths
$\lambda_{\perp 0}$ are limited by the cross-$\mathbf{B}_{0}$ scale
$l_{\perp}$ of density filaments comprising the loops, $\lambda_{\perp
0}\lesssim l_{\perp}$ (wavenumber $k_{\perp 0}\gtrsim 2\pi/l_{\perp}$). Note
that $l_{\perp}$ can be significantly smaller than the visible loop width $r$.
On the contrary, the coronal plasma is quite homogeneous along
$\mathbf{B}_{0}$ and the possible parallel wavelength $\lambda_{z0}$ are
restricted by the loop length $L_{z}$, $\lambda_{z0}=2\pi/k_{z0}\leq L_{z}$
(wavenumber $k_{z0}\gtrsim 2\pi/L$).

For the wavelengths within the mentioned above limits, a large spectral flux,
and hence a strong plasma heating, can be established by the strong turbulence
driven at the critically-balanced anisotropy $\alpha k_{\perp 0}/k_{z0}=$
$\alpha\lambda_{z0}/\lambda_{\perp 0}=$ $v_{A}/v_{k0}$. The corresponding
energy flux injected in the unit volume is
$\epsilon_{cor}\approx\alpha_{0}\rho v_{k0}^{3}/\lambda_{\perp 0}$. Assuming
that the turbulent velocity at injection $v_{k0}$ is observed as the
nonthermal velocity, $v_{k0}\approx v_{\mathrm{nt}}$, and taking from Xie et
al. (2017) $v_{\mathrm{nt}}=30$ km s-1, magnetic field $B=30$ G, and density
$n_{e}=2\times 10^{9}$ cm-3, we obtain the energy flux
$\epsilon_{cor}\sim\alpha_{0}\rho v_{\mathrm{nt}}^{3}/l_{\perp}\sim$ $3\times
10^{-4}\left(l_{\perp}/L_{\perp}\right)^{-1}$ erg cm-3 s-1. For sufficiently
small $l_{\perp}\lesssim 0.1L_{\perp}$, the energy flux
$\epsilon_{cor}\gtrsim$ $3\times 10^{-3}$ erg cm-3 s-1 is enough to heat
typical coronal loops. The corresponding parallel wavelengths at injection are
$\lambda_{z0}\sim\left(\alpha v_{\mathrm{nt}}/v_{A}\right)^{-1}\lambda_{\perp
0}\lesssim 0.5L_{z}$. Therefore, the turbulent cascade and related plasma
heating can be effective if the perpendicular length scales of the loop
substructures are about 10 times smaller than the loop width, which implies
that the loops should be structured more than was required by previous
turbulent heating models.

##  5 Conclusions

In the framework of compressional MHD, we studied numerically the spectral
transport produced by the collisions between counter-propagating Alfvén waves.
The initial two waves are linearly polarized in two orthogonal planes and
their cross-field profiles vary normally to their polarization planes.
Polarization and spectral characteristics of the perturbations generated after
single and multiple collisions between such AWs are analyzed in detail. The
main properties of the resulting spectral transfer are as follows:

  * •

the perturbations generated by AW collisions have smaller scales than the
original waves, which supports turbulence scenarios based on the direct
turbulent cascade generated by AW collisions;

  * •

we observed two regimes of the AW interaction: the first one is typical for
the case of strong turbulence, and the second one is governed by a different
mechanism;

  * •

the spectral transfer generated by the AW collisions is strongly suppressed at
the scales satisfying the classic critical balance condition (11) of
incompressional MHD, which makes the turbulence weak at these scales;

  * •

the strong turbulence is re-established at significantly smaller perpendicular
scales satisfying the modified critical balance condition (12);

We used these properties to re-evaluate the turbulent heating of the solar
coronal loops. The main conclusion is that the turbulent cascade can heat the
loop plasma provided the loop is structured and the characteristic widths of
the loop sub-structures are more than 10 times smaller than the loop width.

## References

  * Belcher & Davis (1971) Belcher, J. W. & Davis, L. 1971, Journal of Geophysical Research, 76, 3534
  * Beresnyak (2014) Beresnyak, A. 2014, The Astrophysical Journal, 784, L20
  * Beresnyak (2015) Beresnyak, A. 2015, The Astrophysical Journal, 801, L9
  * Beresnyak & Lazarian (2008) Beresnyak, A. & Lazarian, A. 2008, ApJ, 682, 1070
  * Brooks & Warren (2016) Brooks, D. H. & Warren, H. P. 2016, The Astrophysical Journal, 820, 63
  * Bruno & Carbone (2013) Bruno, R. & Carbone, V. 2013, Living Reviews in Solar Physics, 10, 2
  * Chandran & Perez (2019) Chandran, B. D. G. & Perez, J. C. 2019, Journal of Plasma Physics, 85, 905850409
  * Culhane et al. (2007) Culhane, J. L., Harra, L. K., James, A. M., et al. 2007, Sol. Phys., 243, 19
  * De Moortel et al. (2014) De Moortel, I., McIntosh, S. W., Threlfall, J., Bethge, C., & Liu, J. 2014, Astrophysical Journal Letters, 782
  * Doschek et al. (2014) Doschek, G. A., McKenzie, D. E., & Warren, H. P. 2014, The Astrophysical Journal, 788, 26
  * Drake et al. (2016) Drake, D. J., Howes, G. G., Rhudy, J. D., et al. 2016, Physics of Plasmas, 23, 022305
  * Drake et al. (2013) Drake, D. J., Schroeder, J. W. R., Howes, G. G., et al. 2013, Physics of Plasmas, 20, 072901
  * Drake et al. (2014) Drake, D. J., Schroeder, J. W. R., Shanken, B. C., et al. 2014, IEEE Transactions on Plasma Science, 42, 2534
  * Galtier et al. (2000) Galtier, S., Nazarenko, S. V., Newell, A. C., & Pouquet, A. 2000, Journal of Plasma Physics, 63, 447
  * Gogoberidze & Voitenko (2016) Gogoberidze, G. & Voitenko, Y. M. 2016, Ap&SS, 361, 364
  * Goldreich & Sridhar (1995) Goldreich, P. & Sridhar, S. 1995, The Astrophysical Journal, 438, 763
  * Gupta et al. (2019) Gupta, G. R., Del Zanna, G., & Mason, H. E. 2019, Astronomy & Astrophysics, 627, A62
  * Howes & Nielson (2013) Howes, G. G. & Nielson, K. D. 2013, Physics of Plasmas, 20, 072302
  * Iroshnikov (1963) Iroshnikov, P. S. 1963, Astron. Zh., 40, 742
  * Kontar et al. (2017) Kontar, E. P., Perez, J. E., Harra, L. K., et al. 2017, Physical Review Letters, 118, 155101
  * Kraichnan (1965) Kraichnan, R. H. 1965, Physics of Fluids, 8, 1385
  * Lucek & Balogh (1998) Lucek, E. A. & Balogh, A. 1998, ApJ, 507, 984
  * Magyar et al. (2019) Magyar, N., Van Doorsselaere, T., & Goossens, M. 2019, The Astrophysical Journal, 873, 56
  * Mallet et al. (2015) Mallet, A., Schekochihin, A. A., & Chandran, B. D. 2015, Monthly Notices of the Royal Astronomical Society: Letters, 449, L77
  * McLaughlin et al. (2011) McLaughlin, J. A., De Moortel, I., & Hood, A. W. 2011, Astronomy & Astrophysics, 527, A149
  * Montgomery & Matthaeus (1995) Montgomery, D. & Matthaeus, W. H. 1995, The Astrophysical Journal, 447, 706
  * Morton et al. (2016) Morton, R. J., Tomczyk, S., & Pinto, R. F. 2016, The Astrophysical Journal, 828, 89
  * Ng & Bhattacharjee (1996) Ng, C. S. & Bhattacharjee, A. 1996, The Astrophysical Journal, 465, 845
  * Nielson et al. (2013) Nielson, K. D., Howes, G. G., & Dorland, W. 2013, Physics of Plasmas, 20, 072303
  * Perez et al. (2020) Perez, J. C., Azelis, A. A., & Bourouaine, S. 2020, Phys. Rev. Research, 2, 023189
  * Perez & Boldyrev (2009) Perez, J. C. & Boldyrev, S. 2009, Phys. Rev. Lett., 102, 025003
  * Pezzi et al. (2017a) Pezzi, O., Malara, F., Servidio, S., et al. 2017a, Phys. Rev. E, 96, 023201
  * Pezzi et al. (2017b) Pezzi, O., Parashar, T. N., Servidio, S., et al. 2017b, Journal of Plasma Physics, 83, 705830108
  * Porth et al. (2014) Porth, O., Xia, C., Hendrix, T., Moschou, S. P., & Keppens, R. 2014, The Astrophysical Journal Supplement Series, 214, 4
  * Shestov et al. (2017) Shestov, S. V., Nakariakov, V. M., Ulyanov, A. S., Reva, A. A., & Kuzin, S. V. 2017, The Astrophysical Journal, 840, 64
  * Sridhar & Goldreich (1994) Sridhar, S. & Goldreich, P. 1994, The Astrophysical Journal, 432, 612
  * Tu et al. (1990) Tu, C. Y., Marsch, E., & Rosenbauer, H. 1990, Geochim. Res. Lett., 17, 283
  * Van Ballegooijen et al. (2011) Van Ballegooijen, A. A., Asgari-Targhi, M., Cranmer, S. R., & DeLuca, E. E. 2011, Astrophysical Journal, 736, 28
  * Verdini, A. et al. (2012) Verdini, A., Grappin, R., & Velli, M. 2012, A&A, 538, A70
  * Verniero & Howes (2018) Verniero, J. L. & Howes, G. G. 2018, Journal of Plasma Physics, 84, 905840109
  * Verniero et al. (2018) Verniero, J. L., Howes, G. G., & Klein, K. G. 2018, Journal of Plasma Physics, 84, 905840103
  * Xie et al. (2017) Xie, H., Madjarska, M. S., Li, B., et al. 2017, The Astrophysical Journal, 842, 38
  * Zheng et al. (2016) Zheng, J., Chen, Y., & Yu, M. 2016, Phys. Scr., 91, 015601


# Towards High-Payload Admittance Control for Manual Guidance with
Environmental Contact

Kevin Haninger∗, Marcel Radke∗, Axel Vick, Jörg Krüger * Equal contribution,
alphabetical order. Authors affiliated with the Department of Automation at
the Fraunhofer Institute for Production Systems and Design Technology. Email:
<EMAIL_ADDRESS>This project has received funding from
the European Union’s Horizon 2020 research and innovation programme under
grant agreement No 820689 — SHERLOCK.

###### Abstract

Force control enables hands-on teaching and physical collaboration, with the
potential to improve ergonomics and flexibility of automation. Established
methods for the design of compliance, impedance control, and  collision
response can achieve free-space stability and acceptable peak contact force on
lightweight, lower payload robots. Scaling collaboration to higher payloads
can allow new applications, but introduces challenges due to the more
significant payload dynamics and the use of higher-payload industrial robots.

To achieve high-payload manual guidance with contact, this paper proposes and
validates new mechatronic design methods: standard admittance control is
extended with damping feedback, compliant structures are integrated to the
environment, and a contact response method which allows continuous admittance
control is proposed. These methods are compared with respect to free-space
stability, contact stability, and peak contact force. The resulting methods
are then applied to realize two contact-rich tasks on a 16 kg payload (peg in
hole and slot assembly) and free-space co-manipulation of a 50 kg payload.

##  I Introduction

Figure 1: Manual guidance with a heavy payload (50 kg) on a COMAU AURA robot,
with 170 kg payload capacity and safety touch sensitive skin.

Physical interaction with a robot allows intuitive position teaching and
collaboration, which can improve task programming and online flexibility.
Force control, impedance control, and other variants [1, 2, 3] can enable safe
physical human-robot interaction (HRI), allowing robot position to be modified
by human or environmental forces, which can reduce cycle time [4] and human
exertion [5].

However, contact couples the robot with unknown environment or human dynamics,
where stability and maximum force bounds should be maintained. Robots with
joint torque sensing can implement impedance control with inner-loop torque
control [6], reducing rendered impedance to improve safety and transient
contact force. However, the payload currently achievable with joint torque
sensing is limited – commercially, to 14 kg. Impedance control via current
monitoring, where joint torque is estimated from motor current, also does not
scale to the higher gear reductions currently required for high payloads, as
gearbox friction dominates external forces.

While robots for physical HRI today have limited payload and reach, there are
benefits to higher-payload collaboration. The ergonomic benefits of HRI
increase at higher payloads, where the robot can reduce the weight carried by
the human [7]. Robot support of bulky objects can allow tasks which previously
required two workers to be done with one [8, 9]. However, increased payload
inertia introduces challenges for responsiveness, stability, and safety [10,
4, 11, 12].

Human safety is critical, especially at higher payloads and with high
impedance industrial robots. Hands-on-payload interaction makes it difficult
to use an enabling switch (which the human must press to enable motion), so
safety via power and force limiting must play a larger role. Here, methods for
the reduction of transient and quasi-static contact forces are proposed and
compared, towards verifiable power and force limiting. However, additional
aspects for certification (e.g. Performance Level of the sensor and
communication protocols) are outside the scope of this paper.

High payload HRI has been demonstrated by separating human input from the
payload with either a separate force/torque (F/T) sensor [4, 12] or a low-
impedance input channel [11]. However, direct co-manipulation of the payload
as seen in Figure 1 simplifies the system, allowing a single sensor for
environmental contact, human input, and collision detection, as well as
increasing the flexibility in how the operator grips the payload [9]. The
ability to set the rotation center or virtual constraints can also improve
ergonomics and task performance [13].

Using an industrial robot for admittance control can simplify mechanical
integration, but the high-impedance hardware and inner-loop position control
makes contact stability more difficult [14] and increases collision forces
[15, 16]. Collision is especially critical with larger payloads as it provides
impulsive input with a broad frequency range which can excite payload
resonance modes.  Many control architectures proposed for high-payload co-
manipulation do not experimentally test contact [17, 3].

These issues can be partially addressed through control. Peak contact force
and contact stability can be improved through feedforward control, which can
be the inverse of the robot dynamics [18], a PID compensator [19] or a lead
controller [3]. In hands-on-payload applications, the compensation of inertial
payload forces has been proposed to increase sensitivity [10, 3, 20].
Alternative force control architectures can be used, such as disturbance
observer-based control [21] or an arbitrary finite impulse response force
controller [22]. However, standard admittance control has well-studied
stability properties and allows methods to adapt the admittance parameters
online [23, 24].

However, control has limitations. The ability to reduce payload oscillation
through control is limited as well as the ability to reduce the transient
contact force [15]. As such, this paper also considers the use of collision
detection and environmental compliance. While environmental compliance is
often used informally in admittance control [21, 25], we show it can be
integrated into compliant tables to allow practical contact-rich tasks.

This paper identifies and partially addresses new challenges of high-payload
co-manipulation, mostly resulting from the significant payload inertia and
resulting lower frequency payload resonance. Compared to published hands-on-
payload co-manipulation [9, 10] (14 and 12 kg payloads), we demonstrate co-
manipulation up to a 50 kg payload. Compared to state of the art admittance
control [3], we add damping feedback, collision response, and validate with
physical experiments. First, the dynamic model for robot, payload and control
architecture are proposed and analyzed. A collision response method which
allows continuous co-manipulation in contact is proposed. Contact experiments
compare  the proposed methods with established variable admittance control and
payload acceleration compensation. Two contact-rich co-manipulation tasks are
demonstrated (double peg-in-hole and a slot assembly task) with a 16 kg
payload, and free-space manipulation of a 50 kg payload.  Finally, design
conclusions are drawn for high-payload admittance control with contact.

##  II Admittance Architecture, Analysis, and Implementation

This section introduces the dynamic model, admittance controller, and analyzes
passivity and contact resonance.

###  II-A Control Architecture

The controller, robot dynamics and environment can be seen in Figure 2, shown
for a single DOF in operational space. The rotational DOFs realize the desired
admittance about the X/Y/Z axis of the TCP, and $F$ is transformed into this
coordinate system before payload compensation, both of which occur in the
block $S$.

In green is the closed position-controlled robot, which receives TCP position
commands $x^{d}$ at $1250$ Hz and uses a classical cascade position controller
with hierarchical position/velocity/current control loops. We simplify the
robot controller $R_{c}$ and dynamics $R_{d}$ by considering the closed-loop
robot dynamics  in Laplace domain as $R(s)=V(s)/X^{d}(s)$, where the physical
impact of forces $F$ on the robot are ignored, assumed to be rejected by the
high-gain position control.

The environment can be seen in Figure 2 in blue, where a stiffness $K_{s}$
couples a payload $P=(M_{p}s+B_{p})^{-1}$ to the robot flange. The environment
$E$ and human $H$ apply forces to the payload in parallel  and are both shown
in Figure 2 for completeness. In the sequel we consider $E$ as a pure
stiffness ($E=K_{e}/s$), and do not consider $H$ in closed-loop analysis as
typical human arm stiffnesses are $<3$ N/mm [26], so are typically dominated
by environment stiffness (here, $>20$ N/mm).

Figure 2: Admittance control architecture, where the implemented controller
(red) sends commands to a position controlled robot (green). The robot has a
F/T sensor and gripper with stiffness $K_{s}$ which couples the payload $P$ to
the robot.

We extend the control architecture of [3] with a feedback controller from
robot position $C_{fb}$ as an additional damping term  based on the real robot
velocity. We note that damping in the admittance controller is based on the
virtual velocity of the admittance integrator ($\dot{x}^{d}$), which, due to
delay and robot dynamics, does not correspond to the real robot velocity.
However, damping  based on the real robot velocity is critical for passivity
and stable interactive control [27]. The following control blocks  from Figure
2 are defined  in Laplace domain as:

$\displaystyle A(s)$ $\displaystyle=\frac{1}{M_{a}s^{2}+B_{a}s}$ (1)
$\displaystyle C_{ff}(s)$ $\displaystyle=K_{l}s$ (2) $\displaystyle C_{fb}(s)$
$\displaystyle=B_{fb}\frac{\omega_{v}s}{s+\omega_{v}}$ (3) $\displaystyle
O_{a}(s)$
$\displaystyle=\frac{\omega_{v}\omega_{a}s^{2}}{(s+\omega_{a})(s+\omega_{v})}$
(4)

where  parameters $M_{a}$ and $B_{a}$ are the admittance mass and damping,
$K_{l}$ a lead controller gain, $B_{fb}$ a feedback damping term, $\omega_{v}$
and $\omega_{a}$ filter constants for the velocity and acceleration
estimation.  The Figure 2 block $\hat{P}=M_{p}$ is the estimated payload
inertia.

The following transfer functions can be written, omitting the argument $(s)$
for compactness:

$\displaystyle\frac{V}{F}=G$
$\displaystyle=-\frac{A(1+C_{ff})R}{1-A(1+C_{ff})R\hat{P}^{-1}O_{a}s^{-1}+RC_{fb}s^{-1}},$
(5) $\displaystyle\frac{F}{V}=\overline{E}$
$\displaystyle=\frac{K_{s}s^{-1}(1+PE+PH)}{1+PE+PH+PK_{s}s^{-1}},$ (6)
$\displaystyle\frac{F}{F^{d}}$
$\displaystyle=\frac{-G\overline{E}}{1+G\overline{E}}.$ (7)

###  II-B Controller Design

The dynamics of (5) and (7) will be used for controller design. Stability with
environmental contact is a foundational problem for admittance control [28],
and various stability conditions have been used: passivity [2, 3], mixed
small-gain/pasisivity [29], or assuming environmental dynamics and showing
closed-loop stability [30, 31].

In the following analysis, control parameters are $M_{a}=10$, $B_{a}=1000$,
$K_{l}=0.02$, $B_{fb}=2E^{-5}$, $\omega_{a}=200$, $\omega_{v}=120$, and
dynamic parameters $K_{s}=3E^{5}$, $E=5E^{5}/s$, $H=0$, $P=(16s+630)$ and
$R=(100s^{2}+400s)(5s^{2}+120s+400)^{-1}$, where the robot dynamics are the
derivative of a well-damped position-controlled inertia with $3$ Hz bandwidth.
An input time delay of $T_{d}=2.4e-3$ is considered on $R$. These values of
$\omega_{a}$ and $\omega_{v}$ reflect those which allowed acceptable noise
characteristics on the hardware for the acceleration compensation, seen in
Figure 3. Idealized payload acceleration compensation is considered where
$O_{a}=s^{2}$ instead of with a low-pass filter as in (4).

Figure 3: Force measurements during jogging, with and without compensation of
payload acceleration. The payload compensation substantially increases noise
in the force, from RMS value of $0.14$ to $0.40$ N at rest.

####  II-B1 Rendered admittance in free space

The rendered admittance of the various controllers are compared with $E=H=0$,
comparing $F\rightarrow V$ as seen in Figure 4a. From the magnitude plot of
Figure 4a, it can be seen that the lead compensator increases the admittance
and ideal payload acceleration compensation increases it further, improving
responsiveness.  The damping feedback controller decreases admittance in the
range 1-20 Hz, reducing responsiveness.

(a) Measured force $F$ to robot velocity $v$ in free space ($E=H=0$) with
delay $T_{d}=2.4E^{-3}$.

(b) Desired force $F^{d}$ to force $F$ in contact ($E=5E^{+5}/s,H=0$)

Figure 4: Bode magnitude plot for transfer function $F\rightarrow v$ in free
space (a) and $F^{d}\rightarrow F$ in environmental contact (b) under various
controllers.

####  II-B2 Force tracking in contact

We then compare the magnitude of $F^{d}\rightarrow F$ when $E=5E^{5}/s$ as
seen in Figure 4b. Standard admittance control has a high resonance peak,
which is eliminated by the feedforward lead controller and slightly reduced by
the feedback controller. Ideal payload acceleration compensation reduces
resonance, but the low-pass filtered version  is not as effective. From this,
we conclude that the benefits of the acceleration compensation depend on the
bandwidth of the acceleration estimation.

Contact stability is checked at a range of parameter values, and the impact of
the controller components on minimum stable admittance damping $B_{a}$ can be
seen in Fig. 6(a) over a range of environment stiffness. At higher $K_{e}$,
payload compensation requires higher damping, but the feedback and feedforward
control allow lower stable $B_{a}$, improving contact performance.

(a) $16$ kg payload, instability with admittance control and $K_{l}=0.025$

(b) $50$ kg payload, force oscillation during jogging

Figure 5: Payload oscillation (a) admittance control with $K_{l}=0.025$ and
$16$ kg payload, unstable mode at $25$ Hz (b) during jogging with $50$ kg
payload, dominant mode at around $9$ Hz.

(a)

(b)

Figure 6: Impact of  (a) Environment stiffness and controller components on
minimum stable damping $B_{a}$ and (b) payload natural frequency
$\omega=\sqrt{M_{p}/K_{s}}$ and corresponding damping ratio $\zeta$ on stable
parameter bounds.

####  II-B3 Instability from high lead gain

While lead gain $K_{l}$ improves both free-space admittance and contact
resonance, it can drive instability due to the resonance in the payload, even
in free space as seen in Figure 5(a). This instability does not appear
analytically from passivity analysis of $F\rightarrow V$, but it can be found
through classical linear stability when delay and payload resonance are
modelled. We note that this problem is similar to that of non-collocation [28,
30], but here the resonance is after the sensor, not between sensor and
actuator.

The impact of payload resonance  on stability bounds is shown in Figure 6,
where payload mass and damping $M_{p},B_{p}$, $P=M_{p}s+B_{p}$ are varied to
achieve the resonant natural frequency $\omega=\sqrt{M_{p}/K_{s}}$ and damping
ratio $\zeta=B_{p}/(2\sqrt{M_{p}K_{s}})$. The limit of stable lead gain
$K_{l}$ and damping $B_{a}$ are found numerically, checking linear stability
of (7) in contact with $E=5E^{5}/s$. As resonant frequency decreases, the
minimum stable damping increases and the maximum stable lead gain decreases,
both limiting performance. This trend is exaggerated if the payload is poorly
damped (red). Thus, the low-frequency and lightly-damped payload dynamics
typical with a higher payload set limits on the stable admittance parameters.

###  II-C Controller Implementation

The transfer functions (1)-(4) are discretized with a Tustin transformation at
a sample time of $h=8E^{-4}$, giving the following equation for admittance and
feedforward controller

$\begin{split}(1+\frac{B_{a}h}{2M_{a}})x^{d}_{t+1}=2x^{d}_{t}-(1-&\frac{B_{a}h}{2M_{a}})x^{d}_{t-1}\\\
+\frac{h^{2}}{4M_{a}}\left(\left(1+2h^{-1}K_{l}\right)F_{t+1}\right.&\left.+2F_{t}+\left(1\texttt{-}2h^{-1}K_{l}\right)F_{t\texttt{-}1}\right),\end{split}$
(8)

with discrete time samples $x^{d}_{t}$ and $F_{t}$. Note that Euler
discretization was found to have worse contact stability, as shown in Figure
9.

An ME-Meßsysteme KD6110 F/T sensor is used, with rated limits of $5$ kN and
$250$ Nm. The GSV-8DS amplifier is used, with an internal analog low-pass
filter with a cutoff frequency of 2500 Hz. A deadband of $5$ N and $1$ Nm is
applied after the payload gravity compensation - all forces below this value
are set to zero. This deadband keeps small gravitational compensation errors
from causing robot motion. The RMS sensor noise at rest is $0.14$ N without
acceleration compensation.

##  III Collision Detection

Contact can be made with pure force control, where the contact forces are
regulated according to the admittance, reducing the quasi-static force in
collision. Alternatively, contact can be detected and used to change robot
behavior, which has additional complexity but can reduce peak contact force
[15, 32].

There are two challenges in co-manipulation contact detection. First, human
and environmental contact forces must be distinguished [15], where a common
strategy is to assume human input is mostly low frequency, and use a threshold
on a high-pass filtered force signal to distinguish contact [25]. A second
challenge is that in co-manipulation, often the task is not finished at the
first contact - fine positioning must still be carried out. Thus, the
collision response of  a permanent robot stop may be unsuitable. In co-
manipulation, the collision response of increasing of damping has been
proposed [33, 5, 25], which allows continuous co-manipulation. However, this
has a minimal impact on the initial peak contact force, which can be critical
in high payloads.

Algorithm 1 Proposed collision response

Initialize:

$\tt was\textunderscore contact=False$

while true do

${\tt contact}=F_{t+1}^{j}\mathtt{sign}(v^{j})>\overline{F}\,\,{\tt
for\,\,any}\,\,j=1\dots 3$

if $\tt contact\,and\,not\,was\textunderscore contact$ then

$x^{d,k}_{t+1}\leftarrow x^{d,k}_{t}$ for $k=1\dots 6$

$x^{d,k}_{t-1}\leftarrow x^{d,k}_{t}$ for $k=1\dots 6$

$\tt was\textunderscore contact=True$

else

$x^{d,k}_{t+1}\leftarrow$ (8) for $k=1\dots 6$

end if

end while

To examine the continuity of this response, consider
$x^{d}_{t}=x^{d}_{t-1}+\delta$, then (8) simplifies to

$x^{d}_{t+1}=x^{d}_{t}+\frac{2M_{a}-B_{a}h}{2M_{a}+B_{a}h}\delta+g(f_{t+1},f_{t},f_{t-1}),$
(9)

where $g(\cdot)$ is the contribution from forces and is unchanged from (8). It
can be seen that when $M_{a}>\frac{1}{2}B_{a}h$, a positive $\delta$ is
carried forward to the next step due to the admittance integrator. When
$\delta=0$, this term drops and
$x^{d}_{t+1}=x^{d}_{t}+g(f_{t+1},f_{t},f_{t-1})$, so setting $\delta=0$ by
modifying $x^{d}_{t-1}$ does not cause a discontinuity.

##  IV Contact Experiments

To validate the controller analysis, contact experiments are done without
human contact to improve repeatability. Unless otherwise noted, experiments
use linear mass and damping of $M_{a}=13$, $B_{a}=1000$, $K_{l}=0.013$, and
$B_{fb}=7.5E^{-6}$, and contact is made by setting $F^{d}=30$ N in the
direction of desired contact.

###  IV-A Controller Component Study

The impact of the feedforward/feedback controllers are compared in contact
with a compliant table, seen with yellow foam in Figure 7.

In Figure 8, the improvement of the lead and feedback controllers can be seen
- they reduce peak contact force and have reduced resonance, as predicted in
Figure 4b. The acceleration compensation is shown in Figure 9, where the
damping is decreased from acceleration compensation, and the peak contact
force slightly increased, again as predicted by Figure 4b.

###  IV-B Collision Detection

Contact responses with and without the collision response method proposed in
Section III are carried out on a smaller payload robot (a COMAU Racer 7, with
a 6.3 kg payload attached), reducing damping to $B_{a}=800$. The collision
detection threshold of $\overline{F}=7$ N is used, and the reference force of
$F^{d}=30$ N is applied for the no collision detection and proposed collision
detection. Additionally, the collision response of making a full stop, where
$x^{d}_{t}$ is not updated after detection (thus ignoring $F$ and $F^{d}$
after collision).  The full stop keeps the robot motion control active, but
keeps the position command $x^{d}$ at the same value. This halts the external
admittance control, which is typically unsuitable for co-manipulation, and is
thus included for only comparison.

The results can be seen in Figure 10, where the proposed collision detection
reduces peak contact force from $44$ to $36$ N. The admittance desired
velocity $\dot{x}^{d}$ is also displayed, calculated by the numerical
difference of $x^{d}_{t}$ and smoothed with a moving average filter of length
$5$. The immediate zeroing of $\dot{x}^{d}$ in contact can be seen, and the
continuity in $x^{d}$ is maintained. As can be seen in Figure 10(a), the full
stop eliminates force overshoot and oscillation, but cannot be continuously
moved after contact without resetting the control. The improvement of the
collision detection depends on the control parameters and environment
stiffness. When $K_{l}$ is increased from $0.013$ to $0.02$ in (b), there is
almost no force overshoot to reduce. When the environment is made more
compliant, $K_{e}$, from $32$ to $15$ N/mm in (c), the collision detection
improves performance, as the lower stiffness takes more time to build forces
to stop the pure force control system.

###  IV-C Environment Compliance

Compliance is critical for transient forces in contact, and is often used
informally for admittance control tests on industrial robots [21, 25].
However, integrating compliance at the contact surface is not always feasible
for certain tasks. Here, we integrate compliant joints [34] into the table
structure for the task, as can be seen in Figure 7. The compliant joints form
a four-bar mechanism which allows the contact surface to move laterally,
reducing the stiffness seen at the contact surface. The location of compliance
has been shown to impact the ability to detect collision [32],  and a hard
contact surface here can excite the payload resonance. However, integrating
compliance deeper into the environment structure allows reduction in peak
contact force as seen in Table I. As seen in Fig. 10, the benefits of the
proposed collision detection increase at lower environment stiffness.

$K_{l}$ | $B_{r}$ |  $\hat{M}$ (kg) |  $K_{e}$ (N/mm) |  $\max F$ (N) |  $\omega$ (Hz)
---|---|---|---|---|---
0 | 0 | 0 | $32.3$ | 77.5 | 2.1
0.02 | 0 | 0 | $32.3$ | 72.5 | 2.5
0.02 | 7.5$E^{-6}$ | 0 | $32.3$ | 58.2 | 2.6
0.02 | 7.5$E^{-6}$ | 16 | $32.3$ | 66.4 | 3.2
0.02 | 7.5$E^{-6}$ | 0 | $11.2$ | 55.6 | 3.1
0.02 | 7.5$E^{-6}$ | 0 | $5.1$ | 39.6 | 2.3
0.02 | 7.5$E^{-6}$ | 0 | $2.7$ | 34.4 | 14.7
Table I: Control parameters, environment stiffness $K_{e}$, resulting maximum
contact force, and frequency of dominant frequency mode $\omega$. Figure 7:
Test environment for slot assembly (left), peg-in-hole (middle) and a
compliant table (right) for contact tests. Figure 8: Impact of feedforward and
feedback on contact performance with the yellow compliant table with
$M_{a}=13$, $B_{a}=1000$. Figure 9: The impact of acceleration compensation
and discretization method on contact performance with the yellow compliant
table with $M_{a}=10$, $B_{a}=1000$.

(a) $K_{l}=0.013$, $K_{e}=32$ N/mm

(b) $K_{l}=0.020$, $K_{e}=32$ N/mm (c) $K_{l}=0.013$, $K_{e}=15$ N/mm

Figure 10:  Collision response method in contact, compared with standard
admittance control and a full stop, at various environment stiffness $K_{e}$.
Position $x^{d}$ is scaled 5x for visibility. Peak force is reduced by
collision response, but makes less difference with higher $K_{l}$ and higher
$K_{e}$.

###  IV-D Variable Admittance

Adaptation of the admittance parameters has been proposed to balance free
space responsiveness with good contact properties. We implement the adaptation
rule proposed in [24], where $B_{a}={\tt max}(800e^{-8|v|},550)$, where $v$ is
robot speed in m/s.  This adaptation rule is similar to many others proposed
for co-manipulation [23] and is easily implemented. $B_{a}$ is updated at 100
Hz with a limit in change of $20$ between time steps, and the admittance mass
is fixed $M_{a}=13$. The robot is started at a position 30 cm away from
contact, and an operator was asked to bring the robot into contact without
knowing which control mode was active. A typical result can be seen below in
Fig. 11. This was repeated 4 times per configuration, and one experiment
dropped to bring the average contact velocity closest together. The various
configurations and results can be seen in Table II, evaluated according to:
(1) peak contact force, (2) contact velocity, (3) total energy $\sum|F*v|$ in
the direction of motion for contact, and (4) Time from first force to contact.

$K_{l}$ | $B_{fb}$ | C Det | Var $B_{a}$ | $\max F$ | $v_{c}$ | En. | Time
---|---|---|---|---|---|---|---
|  |  |  | 95.3 | 8.4 | 2.97 | 5.1
x |  |  |  | 93.9 | 7.8 | 3.05 | 5.0
x | x |  |  | 81.5 | 8.5 | 3.17 | 4.3
x | x | x |  | 76.8 | 7.8 | 3.44 | 4.2
|  |  | x | 121 | 9.7 | 2.72 | 4.2
x | x |  | x | 81.8 | 7.8 | 3.33 | 3.5
|  | x | x | 116 | 9.3 | 2.83 | 3.4
x | x | x | x | 76.1 | 7.2 | 3.51 | 3.5
Table II: Manual guidance to contact, where C Det is collision
detection/response and Var $B_{a}$ is the variable admittance. $v_{c}$ is
contact velocity (cm/s), En$=\sum_{i}|F_{i}v_{i}|$ (kJ), and Time (sec) is
time from start to contact.

(a)

(b)

Figure 11: Trace of experiment with variable admittance (a), and zoom of force
in contact (b), comparing variable admittance, normal admittance, and the
proposed controller.

From these results, we conclude that variable impedance can improve the free-
space responsiveness of the robot, reducing average completion time from $4.6$
seconds to $3.7$. However, variable admittance does not significantly improve
contact performance, at least in this implementation. As can be seen in 11 (a,
right), the damping is increased after the initial collision peak. While this
result is implementation-specific, it may be difficult to address as
collisions occur in 10-80 ms. Fig. 11(b) shows the proposed approaches improve
the post-contact performance compared with standard and adaptive admittance.

##  V Co-manipulation Demonstrations

The design conclusions above are applied to realize high-payload co-
manipulation tasks.

###  V-A 16 kg Payload

Two contact-rich assembly tasks are carried out with the setup in Figure 7, a
double peg-in-hole assembly and a slot assembly. In both cases, the tolerances
are running (holes are oversized 0.5mm). Compliance is integrated into the
table structure, reducing total environment stiffness while still enabling the
hard surface contact of the task. The linear parameters used are
$K_{l}=0.013$, $B_{r}=7.5E^{-6}$, $\hat{P}=0$, $M_{a}=13$, $B_{a}=1250$ and
rotational $M_{a}=3$, $B_{a}=150$ (rotation units are radian). For the slot
assembly, rotational admittance is decreased $M_{a}=12$, $B_{a}=600$ to allow
manipulation from one side. Stable contact transitions are achieved, as seen
in the attached video 111Video available at: https://youtu.be/ty-oVZs-n1Q.

###  V-B 50 kg Payload

The same control architecture is used with a COMAU AURA robot for a manual
guidance task in free-space, seen in Figure 1. The total payload is 50 kg,
where compliance in the suction gripper system results in a payload resonance
of about $9$ Hz, giving an effective $K_{s}\approx 1.6E^{+5}$. Due to this
resonance, the admittance gains had to be increased and $K_{l}$ decreased.
Additionally, the payload identification routine was adapted to stop for 3
seconds before collecting measurements. Even after this modification, the
payload calibration had larger magnitude of residual errors, the RMS
calibration error was $23.3$ N ($1.1$ was achieved on the 16 kg payload).
While the F/T sensor is rated for $5$ kN, it was factory calibrated with 10 kg
test weights, resulting in more significant nonlinearities with the $50$ kg
payload.

The control parameters for this setup were tested as linear $K_{l}=0.01$,
$B_{r}=0$, $\hat{P}=0$, $M_{a}=20$, $B_{a}=2000$ and rotational $M_{a}=16$,
$B_{a}=1850$. The force clipping threshold was also increased to $20$ N and
$3$ Nm. The resulting performance can be seen in Figure 12 and the attached
video.

Figure 12: Forces and Cartesian position during manual guidance with 50 kg
payload.

##  VI Conclusion

This paper took a mechatronic approach to high-payload manual guidance with
environment contact. A major challenge of high-payload applications is payload
resonance.  The following design conclusions are drawn:

Feedforward lead $K_{l}$ decreases peak contact force and increasing free-
space admittance. However, maximum $K_{l}$ is limited by payload resonance
(Fig. 4).

Feedback damping $B_{fb}$ decreases peak contact force (albeit less
effectively than feedforward), but also reduces the admittance - this resulted
in $16$% increase in human effort on experiments here.

Collision response by zeroing $\dot{x}^{d}$ in the admittance integrator
allows continuous admittance control, and reduces peak collision force by
$42\%$ in Fig. 10(a). However, it is more effective with more environment
compliance and requires additional integration work.

Environment compliance reduces peak contact force (Table I) and changes the
efficacy of control (Fig. 6a) and collision detection (Fig. 10b). It can be
integrated deeper in the environment structure to improve contact task
feasibility (Fig. 7).

Variable damping can reduce the time taken for a comanipulation task, but did
not improve the collision characteristics here.

Payload acceleration theoretically improves free space admittance and contact,
but when the required low-pass filters are used, it increases peak contact
force and reduces damping.

## References

  * [1] D. E. Whitney, “Historical Perspective and State of the Art in Robot Force Control,” _The International Journal of Robotics Research_ , vol. 6, no. 1, pp. 3–14, Mar. 1987.
  * [2] N. Hogan, “Impedance control: An approach to manipulation: Part i—theory,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 107, no. 1, pp. 1–7, 1985.
  * [3] A. Q. Keemink, H. van der Kooij, and A. H. Stienen, “Admittance control for physical human–robot interaction,” _The International Journal of Robotics Research_ , vol. 37, no. 11, pp. 1421–1444, Sept. 2018.
  * [4] A. Lecours, B. Mayer-St-Onge, and C. Gosselin, “Variable admittance control of a four-degree-of-freedom intelligent assist device,” in _2012 IEEE International Conference on Robotics and Automation_. IEEE, 2012, pp. 3903–3908.
  * [5] F. Dimeas and N. Aspragathos, “Online Stability in Human-Robot Cooperation with Admittance Control,” _IEEE Transactions on Haptics_ , vol. 9, no. 2, pp. 267–278, Apr. 2016.
  * [6] A. Albu-Schäffer, C. Ott, and G. Hirzinger, “A unified passivity-based control framework for position, torque and impedance control of flexible joint robots,” _The International Journal of Robotics Research_ , vol. 26, no. 1, pp. 23–39, 2007.
  * [7] C. Gosselin, _et al._ , “A friendly beast of burden: A human-assistive robot for handling large payloads,” _IEEE Robotics & Automation Magazine_, vol. 20, no. 4, pp. 139–147, 2013.
  * [8] P. Franceschi, _et al._ , “Precise Robotic Manipulation of Bulky Components,” _IEEE Access_ , vol. 8, pp. 222 476–222 485, 2020.
  * [9] L. Roveda, _et al._ , “Human-Robot Cooperative Interaction Control for the Installation of Heavy and Bulky Components,” in _2018 IEEE International Conference on Systems, Man, and Cybernetics (SMC)_ , Oct. 2018, pp. 339–344.
  * [10] F. Aghili, “Robust impedance control of manipulators carrying a heavy payload,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 132, no. 5, 2010.
  * [11] J.-M. Audet and C. Gosselin, “Intuitive Physical Human-Robot Interaction using an Underactuated Redundant Manipulator with Complete Spatial Rotational Capabilities,” _Journal of Mechanisms and Robotics_ , pp. 1–15, May 2021.
  * [12] D. Surdilovic, _et al._ , “Kooperation und Kollaboration mit Schwerlastrobotern – Sicherheit, Perspektive und Anwendungen,” in _Mensch-Roboter-Kollaboration_. Springer Fachmedien Wiesbaden GmbH, Edited Hans Buxbaum, 2020.
  * [13] S. A. Bowyer, B. L. Davies, and F. R. y Baena, “Active Constraints/Virtual Fixtures: A Survey,” _IEEE Transactions on Robotics_ , vol. 30, no. 1, pp. 138–157, Feb. 2014.
  * [14] E. Colgate and N. Hogan, “An analysis of contact instability in terms of passive physical equivalents,” in _, 1989 IEEE International Conference on Robotics and Automation, 1989. Proceedings_ , May 1989, pp. 404–409 vol.1.
  * [15] S. Haddadin, A. De Luca, and A. Albu-Schäffer, “Robot collisions: A survey on detection, isolation, and identification,” _IEEE Transactions on Robotics_ , vol. 33, no. 6, pp. 1292–1312, 2017.
  * [16] K. Haninger and D. Surdilovic, “Bounded Collision Force by the Sobolev Norm: Compliance and Control for Interactive Robots,” in _2019 IEEE International Conference on Robotics and Automation (ICRA)_ , 2019, pp. 8259–8535.
  * [17] A. Sidiropoulos, _et al._ , “A variable admittance controller for human-robot manipulation of large inertia objects,” in _2021 30th IEEE International Conference on Robot & Human Interactive Communication (RO-MAN)_. IEEE, 2021, pp. 509–514.
  * [18] D. Surdilovic, “Robust control design of impedance control for industrial robots,” in _Intelligent Robots and Systems, 2007. IROS 2007\. IEEE/RSJ International Conference On_. IEEE, 2007, pp. 3572–3579.
  * [19] H. Cao, _et al._ , “Smooth adaptive hybrid impedance control for robotic contact force tracking in dynamic environments,” _Industrial Robot: the international journal of robotics research and application_ , vol. 47, no. 2, pp. 231–242, Jan. 2020.
  * [20] S. Farsoni, _et al._ , “Compensation of load dynamics for admittance controlled interactive industrial robots using a quaternion-based kalman filter,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 2, pp. 672–679, 2017.
  * [21] S. Katsura, Y. Matsumoto, and K. Ohnishi, “Modeling of force sensing and validation of disturbance observer for force control,” _IEEE Transactions on industrial electronics_ , vol. 54, no. 1, pp. 530–538, 2007.
  * [22] H. Pham and Q.-C. Pham, “Convex Controller Synthesis for Robot Contact,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 3330–3337, Apr. 2020.
  * [23] F. J. Abu-Dakka and M. Saveriano, “Variable impedance control and learning – A review,” _arXiv:2010.06246 [cs]_ , Oct. 2020.
  * [24] F. Ficuciello, L. Villani, and B. Siciliano, “Variable impedance control of redundant manipulators for intuitive human–robot physical interaction,” _IEEE Transactions on Robotics_ , vol. 31, no. 4, pp. 850–863, 2015.
  * [25] F. Ferraguti, _et al._ , “A variable admittance control strategy for stable physical human–robot interaction,” _The International Journal of Robotics Research_ , vol. 38, no. 6, pp. 747–765, May 2019.
  * [26] T. Tsumugiwa, R. Yokogawa, and K. Hara, “Variable impedance control based on estimation of human arm stiffness for human-robot cooperative calligraphic task,” in _Robotics and Automation, 2002. Proceedings. ICRA’02. IEEE International Conference On_ , vol. 1. IEEE, 2002, pp. 644–650.
  * [27] J. E. Colgate and G. G. Schenkel, “Passivity of a class of sampled-data systems: Application to haptic interfaces,” _Journal of robotic systems_ , vol. 14, no. 1, pp. 37–47, 1997.
  * [28] N. Hogan, “On the stability of manipulators performing contact tasks,” _Robotics and Automation, IEEE Journal of_ , vol. 4, no. 6, pp. 677–686, 1988\.
  * [29] K. Haninger and D. Surdilovic, “Identification of Human Dynamics in User-Led Physical Human Robot Environment Interaction,” in _2018 27th International Symposium on Robot and Human Interactive Communication (RO-MAN)_ , 2018, pp. 509–514.
  * [30] S. D. Eppinger and W. P. Seering, “Three dynamic problems in robot force control,” _IEEE Transactions on Robotics and Automation_ , vol. 8, no. 6, pp. 751–758, 1992.
  * [31] S. P. Buerger and N. Hogan, “Relaxing Passivity for Human-Robot Interaction,” in _2006 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , Oct. 2006, pp. 4570–4575.
  * [32] K. Haninger, _et al._ , “Contact Information Flow and Design of Compliance,” in _On Review, ICRA 2022_.
  * [33] V. Okunev, T. Nierhoff, and S. Hirche, “Human-preference-based control design: Adaptive robot admittance control for physical human-robot interaction,” in _RO-MAN, 2012 IEEE_. IEEE, 2012, pp. 443–448.
  * [34] R. Hartisch and K. Haninger, “Co-design of Environmental Compliance for High-speed Contact Tasks,” in _On Review, 2021_.