# Efficiency of MHD Wave Generation in Weakly Ionized Atmospheres

Paul S. Cally School of Mathematics, Monash University
Victoria 3800, Australia

###### Abstract

Generation of Alfvén and slow magneto-acoustic waves in weakly ionized
atmospheres by excitation of the charges-only component of the two fluid
(charges and neutrals) plasma is shown to be more or less efficient depending
on the energy fraction initially allocated to the three stationary “flow
differential” modes which characterize the inter-species drift. This is
explained via detailed analysis of the full ten-dimensional spectral
description of two-fluid linear magnetohydrodynamics. Excitation via the
velocity of the charges only is found to be very inefficient, in accord with
previous results, whilst excitation via the magnetic field perturbation alone
is highly efficient. All ten eigenvalues and eigenvectors are presented
analytically in the high collision frequency regime.

Solar atmosphere(1477) — Plasma astrophysics(1261) —
Magnetohydrodynamics(1964)

##  1 Introduction

The low atmospheres of cool stars are generally very weakly ionized. For
example, the quiet solar photosphere has an ionization fraction as low as
$10^{-4}$ (Khomenko et al., 2014). This has led to some discussion about
whether Alfvén waves can be excited there. Alfvén waves are often invoked as
important drivers of the solar wind and coronal heating (Cranmer & van
Ballegooijen, 2005; McIntosh et al., 2011), so this is an important issue.

To date, the argument has revolved around timescales, and less obviously the
wave initialization mechanism. The electron-ion elastic scattering collision
frequency at the quiet Sun photospheric base is of order $10^{10}$ $\rm
s^{-1}$, which allows the charges to be modelled as a single fluid at
frequencies much lower than this. The coupling between the neutrals fluid and
the charges fluid is then described by the charges-neutrals collision
frequency $\nu_{\mathrm{cn}}$ of around $10^{9}$ $\rm s^{-1}$ and the neutral-
charges collision frequency $\nu_{\mathrm{nc}}$ of about $10^{6}$ $\rm
s^{-1}$, dropping to a few tens of thousands per second at the top of the
photosphere. The two fluids (neutrals and charges) are coupled only via these
collisions.

The standard one-fluid (1F) ideal magnetohydrodynamic (MHD) formula for Alfvén
wave energy flux density (energy per unit area per unit time) is

$F=Ea$ (1)

directed along the field lines, where $E=E_{\mathrm{kin}}+E_{\mathrm{mag}}$ is
the wave energy density made up equally of kinetic and magnetic contributions,
and $a$ is the Alfvén speed. In terms of the plasma velocity amplitude $v$ and
equilibrium density $\rho$,
$E_{\mathrm{mag}}=E_{\mathrm{kin}}={\textstyle\frac{1}{4}}\rho v^{2}$. The
extra factor of ${\textstyle\frac{1}{2}}$ comes from RMS averaging. Thus
$F={\textstyle\frac{1}{2}}\rho v^{2}\,a$ overall. Vranjes et al. (2008) argued
that this should be reduced by the factor $\chi^{-2}$ in a two-fluid (2F,
charges and neutrals) plasma, where the ratio of neutrals to charges density
$\chi=\rho_{\mathrm{n}}/\rho_{\mathrm{c}}$ is typically of order
$10^{3}-10^{4}$. This is because any velocity originally on the charges alone
must quickly be shared with the neutrals, thereby greatly reducing $v$. Based
on this insight, they concluded that a significant flux of Alfvén waves could
not be generated in the solar photosphere.

On the other hand, Tsap et al. (2011) found no such diminution of flux in
their 2F model. The cause of the discrepancy was discussed by Soler et al.
(2013b), who concluded that the issue hinges on whether or not the neutrals
fluid is given initial velocity ${\boldsymbol{v}_{\mathrm{n}}}$ matching that
of the charges fluid ${\boldsymbol{v}_{\mathrm{c}}}$. If it is, then there is
no collisional quenching of the velocity and the full 1F energy flux formula
applies. If on the other hand the neutrals are left unperturbed
(${\boldsymbol{v}_{\mathrm{n}}}=\mathbf{0}$) by the driver, then the
$\chi^{-2}$ reduction must be applied. Soler et al. judged that both
viewpoints are valid; it just depends upon the nature of the perturbation that
initiates the wave.

The nature of the dependence on initial conditions is explored further here
using a spectral description.

###  1.1 Spectral Overview

For specified wavevector $\boldsymbol{k}$, it is shown that the linearized
system is fully described by a tenth order eigenvalue system, corresponding to
six slowly decaying MHD wave modes (slow, Alfvén and fast, each propagating in
the $\pm\boldsymbol{k}$ directions); one stationary very slowly decaying
isobaric mode with zero total pressure; and three very rapidly decaying
(nanoseconds) stationary flow differential modes whose main characteristic is
that they exhibit very different velocities of the neutrals and charges
fluids. Only the MHD modes transport energy, and only the flow differential
modes have significantly discrepant velocities between the two species. The
eigenmodes are independent of each other, so their energies are simply
whatever they were initially given by the excitation process, subject to their
respective fast or slow collisional decays.

Any energy given to the flow differential modes is immediately lost to
collisions, and represents an inefficiency of the excitation. For many simple
excitations, this is around 50%. Energy deposited in the isobaric mode also
does not contribute to energy flux, but typically very little of this mode is
generated. Energy density $E_{i}$ given to each MHD mode (where $i$ represents
one of the six modes) plays a full role in energy transport, carrying flux
${\boldsymbol{F}}_{i}=E_{i}\mathbf{V}_{i}$, where $\mathbf{V}_{i}$ is the
corresponding group velocity for that mode. Typically, realistic excitation
results in nearly equal energies in the positively and negatively directed
versions of each mode type, and so _net_ flux is small or zero, but flux in
each direction can independently be large.

However, pure excitation via just a magnetic field perturbation does not place
significant energy in the flow differential modes, and so essentially all
energy is efficiently allocated to travelling MHD modes, 50% in each
direction. Therefore, a purely magnetic excitation can in fact be a most
efficient generator of MHD waves. The same holds true for full-plasma
excitations where charges and neutrals fluids are given the same initial
velocities, as in Tsap et al. (2011), because again the flow differential
modes are not significantly excited.

The eigensystem decomposition is described in detail in Section 2. Section 3
discusses the quadratic expressions for wave energy and flux. An analytic
analysis of three different charges-only Alfvén wave excitation mechanisms is
presented in Section 4. Numerical examples of both Alfvén and magneto-acoustic
excitation corresponding to different levels in the low solar atmosphere are
set out in Section 5. The results are discussed and summarized in Section 6.

##  2 Mathematical Formulation

###  2.1 Governing Equations

Consider a partially ionized hydrogen plasma. Extension to a more realistic
chemical mixture does not change the arguments to be presented. The linearized
coupled two-fluid (2F) wave equations take the form (Soler et al., 2013a)

$\displaystyle\rho_{\mathrm{n}}\frac{\partial{\boldsymbol{v}_{\mathrm{n}}}}{\partial
t}=-\mbox{\boldmath$\nabla$}p_{\mathrm{n}}+\alpha_{\mathrm{cn}}({\boldsymbol{v}_{\mathrm{c}}}-{\boldsymbol{v}_{\mathrm{n}}}),$
(2a) $\displaystyle\frac{\partial p_{\mathrm{n}}}{\partial t}=-\gamma
P_{\mathrm{n}}\mbox{\boldmath$\nabla$}{\boldsymbol{\cdot}}\hskip
1.5pt{\boldsymbol{v}_{\mathrm{n}}},$ (2b) for the neutrals, and
$\displaystyle\rho_{\mathrm{c}}\frac{\partial{\boldsymbol{v}_{\mathrm{c}}}}{\partial
t}=-\mbox{\boldmath$\nabla$}p_{\mathrm{c}}+\frac{1}{\mu}(\mbox{\boldmath$\nabla$}{\boldsymbol{\times}}{\boldsymbol{b}}){\boldsymbol{\times}}{\boldsymbol{B}}-\alpha_{\mathrm{cn}}({\boldsymbol{v}_{\mathrm{c}}}-{\boldsymbol{v}_{\mathrm{n}}}),$
(2c) $\displaystyle\frac{\partial p_{\mathrm{c}}}{\partial t}=-\gamma
P_{\mathrm{c}}\mbox{\boldmath$\nabla$}{\boldsymbol{\cdot}}\hskip
1.5pt{\boldsymbol{v}_{\mathrm{c}}},$ (2d)
$\displaystyle\frac{\partial{\boldsymbol{b}}}{\partial
t}=\mbox{\boldmath$\nabla$}{\boldsymbol{\times}}({\boldsymbol{v}_{\mathrm{c}}}{\boldsymbol{\times}}{\boldsymbol{B}}),$
(2e)

for the charges, where the equilibrium magnetic field ${\boldsymbol{B}}$ and
charges and neutrals gas pressures $P_{\mathrm{c}}$ and $P_{\mathrm{n}}$ are
perturbed respectively by ${\boldsymbol{b}}$, $p_{\mathrm{c}}$ and
$p_{\mathrm{n}}$. These equations follow from linearizing Equations (1) of
Popescu Braileanu et al. (2019) under the assumption of temperature equality
between charges and neutrals. The equilibrium charges and neutrals mass
densities are $\rho_{\mathrm{c}}$ and $\rho_{\mathrm{n}}$,
$\alpha_{\mathrm{cn}}$ is the friction coefficient, and $\mu$ is the vacuum
permeability. The fluid velocities of the charges and neutrals are
respectively
${\boldsymbol{v}_{\mathrm{c}}}=(u_{\mathrm{c}},v_{\mathrm{c}},w_{\mathrm{c}})$
and
${\boldsymbol{v}_{\mathrm{n}}}=(u_{\mathrm{n}},v_{\mathrm{n}},w_{\mathrm{n}})$.
Without loss of generality, ${\boldsymbol{B}}$ is arbitrarily oriented in the
$z$-direction and the wavevector
$\boldsymbol{k}\equiv-i\,\mbox{\boldmath$\nabla$}=(\sin\theta,0,\cos\theta)k$
in the $x$-$z$ plane. Since ${\boldsymbol{b}}$ is perpendicular to
$\boldsymbol{k}$ (which also follows from
$\mbox{\boldmath$\nabla$}{\boldsymbol{\cdot}}\hskip 1.5pt{\boldsymbol{b}}=0$),
we need only retain two components, $b_{y}$ and $b_{\perp}=b_{x}\cos\theta-
b_{z}\sin\theta$.

These 2F energy equations above are adiabatic, and do not explicitly feed
energy lost via collisions back into the thermal state of the plasma. However,
it is to be understood that this is where it goes (Cally & Gómez-Míguez,
2023).

Although Equations (2) are commonly used in two-fluid studies, Vranjes et al.
(2008) (and more recently Alharbi et al., 2022) questioned the significance of
the Lorentz force term on the right hand side of Equation (2c) on the basis
that the collisional frequency greatly exceeds the ion gyrofrequency
$\Omega_{\mathrm{i}}$ and hence that ions traverse only a tiny portion of
their gyration path before scattering. However, Tsap et al. (2011) demonstrate
that the ratio of the Lorentz force to the net ion-neutral collision force is
of order

$\frac{\Omega_{\mathrm{i}}}{\nu_{\mathrm{in}}}:\frac{|{\boldsymbol{v}_{\mathrm{c}}}-{\boldsymbol{v}_{\mathrm{n}}}|}{|{\boldsymbol{v}_{\mathrm{c}}}|}.$

The left hand side is typically $10^{-3}-10^{-2}$ in the low solar atmosphere
(see Khomenko et al., 2014, Fig. 1). However, rapid collisions easily ensure
that the right hand side is smaller than this, so the Lorentz force dominates
or is at least comparable to the collisional force. Indeed, as will be shown,
the collisions impose two timescales on the system: (i) a very rapid scale of
order nanoseconds over which
$|{\boldsymbol{v}_{\mathrm{c}}}-{\boldsymbol{v}_{\mathrm{n}}}|$ is reduced
almost to zero, and (ii) a much longer ambipolar diffusion timescale which
typically exceeds Alfvén wave periods of interest. The Lorentz force is
unimportant over the former scale, but clearly significant over the latter, so
it is retained here.

###  2.2 Modal Eigenfrequencies and Eigenvectors

Let
$\boldsymbol{X}=(u_{\mathrm{n}},v_{\mathrm{n}},w_{\mathrm{n}},\psi_{\mathrm{n}}=p_{\mathrm{n}}/\rho_{\mathrm{n}},u_{\mathrm{c}},v_{\mathrm{c}},w_{\mathrm{c}},\psi_{\mathrm{c}}=p_{\mathrm{c}}/\rho_{\mathrm{c}},\alpha_{\perp}=b_{\perp}/\sqrt{\mu\,\rho},\alpha_{y}=b_{y}/\sqrt{\mu\,\rho})^{T}$,
where $\rho=\rho_{\mathrm{c}}+\rho_{\mathrm{n}}$ is the total equilibrium
density. Then Equations (2) can be written in matrix form
$d\boldsymbol{X}/dt={\mathrm{M}}\,\boldsymbol{X}$ with exact solution
$\boldsymbol{X}(t)=e^{{\mathrm{M}}\,t}\boldsymbol{X}(0)$. The matrix
${\mathrm{M}}$ is not defective so the matrix exponential may be constructed
by direct exponentiation of the eigenvalues,
$e^{{\mathrm{M}}\,t}={\mathrm{P}}^{T}\mathop{\rm
diag}[e^{\lambda_{1}t},\ldots,e^{\lambda_{10}t}]{\mathrm{P}}^{-1}$, where the
$\lambda_{m}$ are the eigenvalues of ${\mathrm{M}}$ and the columns of the
$10\times 10$ matrix ${\mathrm{P}}$ are the corresponding eigenvectors.

The following quantities are introduced to define ${\mathrm{M}}$: neutral to
charges mass ratio $\chi=\rho_{\mathrm{n}}/\rho_{\mathrm{c}}$, (total) Alfvén
velocity ${\boldsymbol{a}}={\boldsymbol{B}}/\sqrt{\mu\rho}$, and (total) sound
speed $c$ defined in terms of the charges and neutrals squared sound speeds
$c_{\mathrm{c}}^{2}=\gamma P_{\mathrm{c}}/\rho_{\mathrm{c}}$ and
$c_{\mathrm{n}}^{2}=\gamma P_{\mathrm{n}}/\rho_{\mathrm{n}}$ by
$c^{2}=(c_{\mathrm{c}}^{2}+\chi c_{\mathrm{n}}^{2})/(1+\chi)$. Thermal
equilibrium is assumed between the species,
$c_{\mathrm{c}}^{2}=2c_{\mathrm{n}}^{2}$, so
$c^{2}=c_{\mathrm{n}}^{2}(2+\chi)/(1+\chi)$. The total Alfvén speed
$a=B/\sqrt{\mu\rho}$ is related to the charges-only Alfvén speed
$a_{\mathrm{c}}=B/\sqrt{\mu\rho_{\mathrm{c}}}$ by
$a_{\mathrm{c}}=a\sqrt{1+\chi}$. It is natural to use $a$ rather than
$a_{\text{c}}$ as the measure of magnetic influence in the strongly coupled
regime. We also introduce the charges-neutral collision frequency
$\nu_{\mathrm{cn}}=\alpha_{\mathrm{cn}}/\rho_{\mathrm{c}}$ and the neutral-
charges collision frequency
$\nu_{\mathrm{nc}}=\alpha_{\mathrm{cn}}/\rho_{\mathrm{n}}$, related by
$\nu_{\mathrm{cn}}=\chi\nu_{\mathrm{nc}}$.

The matrix ${\mathrm{M}}$ is written down explicitly in Appendix A. In the
absence of collisions, the top left $4\times 4$ and bottom right $6\times 6$
submatrices of ${\mathrm{M}}$ decouple, yielding two oppositely directed
acoustic waves and two steady (zero frequency) incompressive velocity shears
on the neutrals, and the six oppositely directed slow, Alfvén and fast waves
on the charges.

On the other hand, the ten eigenfrequencies $\omega=i\,\lambda$ and the
accompanying eigenvectors of the full collisionally coupled system are quite
different in nature. All eigenfrequencies are complex for finite non-zero
$\nu_{\mathrm{nc}}$ and may be calculated numerically for any combination of
parameters. It is instructive to calculate their asymptotic values in the
large collision rate regime $\nu_{\mathrm{nc}}\gg\omega$ that applies to waves
of interest in the low solar atmosphere.

A perturbation method for deriving these formulae is sketched in Appendix B.

#### Alfvén waves

$\omega_{\text{A}}=a\,k\cos\theta\left(\pm
1-\frac{i\,a\,k\,\chi\cos\theta}{2(1+\chi)\nu_{\mathrm{nc}}}+\mathcal{O}\left((\omega/\nu_{\mathrm{nc}})^{2}\right)\right),$
(3a)

which accords with Equation (21) of de Pontieu & Haerendel (1998) for the two-
dimensional case $\theta=0$. The corresponding eigenvectors are

$A_{\pm}=\left(0,\mp
1-\frac{i\,a\,k_{z}(\chi+2)}{2(\chi+1)\nu_{\mathrm{nc}}}+\mathcal{O}\left(\nu_{\mathrm{nc}}^{-2}\right),0,0,0,\mp
1+\frac{i\,a\,k_{z}\,\chi}{2(\chi+1)\nu_{\mathrm{nc}}}+\mathcal{O}\left(\nu_{\mathrm{nc}}^{-2}\right),0,0,0,1\right),$
(3b)

which verifies that the velocities on the neutrals and charges (second and
sixth components) differ only at $\mathcal{O}(\nu_{\mathrm{nc}}^{-1})$.

#### Slow waves

$\omega_{\text{s}}=\pm\frac{k\sqrt{a^{2}+c^{2}-\Delta}}{\sqrt{2}}-\frac{ik^{2}\chi
H_{\text{s}}}{2\nu_{\mathrm{nc}}(\chi+1)(\chi+2)^{2}\left(a^{2}+c^{2}-\Delta\right)\Delta}+O\left((\omega/\nu_{\mathrm{nc}})^{2}\right),$
(3c)

where $\Delta=\sqrt{a^{4}+c^{4}-2a^{2}c^{2}\cos 2\theta}$ satisfies
$|c^{2}-a^{2}|\leqslant\Delta\leqslant a^{2}+c^{2}$, and

$H_{\text{s}}=a^{2}c^{2}\cos
2\theta\left(a^{2}(\chi+2)(\chi+3)+c^{2}(\chi+3)-\Delta(\chi+2)\right)+\Delta\left(a^{4}(\chi+2)^{2}+a^{2}c^{2}(\chi+2)+c^{4}\right)\\\
-\left(a^{2}(\chi+2)+c^{2}\right)\left(a^{4}(\chi+2)+c^{4}\right).$ (3d)

The corresponding eigenvector up to and including terms of
$\mathcal{O}(\nu_{\mathrm{nc}}^{-1})$ is far too long to present here, so we
give it to leading order only. For the positively (upper sign) and negatively
(lower sign) directed modes, the eigenvector is

$S_{\pm}=\left(\mp\frac{\sqrt{a^{2}+c^{2}-\Delta}}{\sqrt{2}a},0,\pm\frac{\cot\theta\left(a^{2}-c^{2}+\Delta\right)}{\sqrt{2}a\sqrt{a^{2}+c^{2}-\Delta}},\frac{c^{2}(\chi+1)\csc\theta\left(a^{2}\cos
2\theta-c^{2}+\Delta\right)}{a(\chi+2)\left(a^{2}+c^{2}-\Delta\right)},\mp\frac{\sqrt{a^{2}+c^{2}-\Delta}}{\sqrt{2}a},0,\right.\\\
\left.\pm\frac{\cot\theta\left(a^{2}-c^{2}+\Delta\right)}{\sqrt{2}a\sqrt{a^{2}+c^{2}-\Delta}},\frac{2c^{2}(\chi+1)\csc\theta\left(a^{2}\cos
2\theta-c^{2}+\Delta\right)}{a(\chi+2)\left(a^{2}+c^{2}-\Delta\right)},1,0\right)+\mathcal{O}\left(\nu_{\mathrm{nc}}^{-1}\right).$
(3e)

Again the velocities on the neutrals and charges differ only at
$\mathcal{O}(\nu_{\mathrm{nc}}^{-1})$.

#### Fast waves

$\omega_{\text{f}}=\pm\frac{k\sqrt{a^{2}+c^{2}+\Delta}}{\sqrt{2}}-\frac{ik^{2}\chi
H_{\text{f}}}{2\nu_{\mathrm{nc}}(\chi+1)(\chi+2)^{2}\left(a^{2}+c^{2}+\Delta\right)\Delta}+O\left((\omega/\nu_{\mathrm{nc}})^{2}\right),$
(3f)

where

$H_{\text{f}}=-a^{2}c^{2}\cos
2\theta\left(a^{2}(\chi+2)(\chi+3)+c^{2}(\chi+3)+\Delta(\chi+2)\right)+\Delta\left(a^{4}(\chi+2)^{2}+a^{2}c^{2}(\chi+2)+c^{4}\right)\\\
+\left(a^{2}(\chi+2)+c^{2}\right)\left(a^{4}(\chi+2)+c^{4}\right).$ (3g)

The eigenvector $F_{\pm}$ for the positively directed fast mode is exactly as
for $S_{\pm}$ given in Equation (3e) but with $\Delta$ replaced by $-\Delta$.

#### Flow differential mode rapid decay (triple)

$\omega_{\text{diff}}=-i\,(1+\chi)\nu_{\mathrm{nc}}+\mathcal{O}\left(a^{2}k^{2}/\nu_{\mathrm{nc}}\right).$
(3h)

The eigenvectors are

$\begin{split}f_{x}=&\left(-1,0,0,-\frac{ic^{2}k\sin\theta}{\nu_{\mathrm{nc}}(\chi+2)},\chi,0,0,\frac{2ic^{2}k\chi\sin\theta}{\nu_{\mathrm{nc}}(\chi+2)},-\frac{iak\chi}{\nu_{\mathrm{nc}}(\chi+1)},0\right)+\mathcal{O}(\nu_{\mathrm{nc}}^{-2}),\\\\[8.0pt]
f_{y}=&\left(0,-1,0,0,0,\chi,0,0,0,-\frac{iak\chi\cos\theta}{\nu_{\mathrm{nc}}(\chi+1)}\right)+\mathcal{O}(\nu_{\mathrm{nc}}^{-2}),\\\\[8.0pt]
f_{z}=&\left(0,0,-1,-\frac{ic^{2}k\cos\theta}{\nu_{\mathrm{nc}}(\chi+2)},0,0,\chi,\frac{2ic^{2}k\chi\cos\theta}{\nu_{\mathrm{nc}}(\chi+2)},0,0\right)+\mathcal{O}(\nu_{\mathrm{nc}}^{-2}),\end{split}$
(3i)

corresponding to flows in the $x$, $y$ and $z$ directions respectively. To
leading order, all three of these modes have oppositely directed velocities
with ratio
$|{\boldsymbol{v}_{\mathrm{c}}}|/|{\boldsymbol{v}_{\mathrm{n}}}|=\chi$ and no
thermal or magnetic components, so their energy is purely kinetic to this
order. All other eigenmodes have near-equal velocities in the highly
collisional regime. The rapid decay of these flow eigenmodes represents the
decay of the flow differentials.

#### Isobaric mode slow decay

$\omega_{\text{pr}}=-\frac{2\,i\,c^{2}k^{2}(1+\chi)}{\nu_{\mathrm{nc}}(2+\chi)^{2}}+\mathcal{O}\left((c\,k/\nu_{\mathrm{nc}})^{3}\right)c\,k.$
(3j)

The eigenvector is

$\left(\frac{ik\sin\theta}{\nu_{\mathrm{nc}}},0,-\frac{ik\cos\theta\left(\chi\tan^{2}\theta-2\right)}{\nu_{\mathrm{nc}}(\chi+2)},-1,0,0,-\frac{ik\chi\sec\theta}{\nu_{\mathrm{nc}}(\chi+2)},\chi,0,0\right)+\mathcal{O}(\nu_{\mathrm{nc}}^{-2}).$
(3k)

This mode balances the pressures on the neutrals and charges,
$p_{\mathrm{c}}+p_{\mathrm{n}}=0$, with small
$\mathcal{O}(\nu_{\mathrm{nc}}^{-1})$ velocity ($u$ and $w$) and magnetic
($b_{\perp}$) perturbations. Hence it is stationary apart from the slow
diffusive decay.

To the order shown, the three classic MHD waves inherit a slow temporal decay
rate of order $1/\nu_{\mathrm{nc}}$, even for $\chi\gg 1$. The decay rate of
the erstwhile stationary isobaric mode is even slower, of order
$1/(\chi\nu_{\mathrm{nc}})$ in the weakly ionized regime. In stark contrast
though, the three modes representing an imbalance between
${\boldsymbol{v}_{\mathrm{c}}}$ and ${\boldsymbol{v}_{\mathrm{n}}}$ decay
extremely rapidly, with rate exactly $\nu_{\mathrm{nc}}+\nu_{\mathrm{cn}}$
independent of ionization fraction. For Alfvén waves, this is in accord with
the analysis of Section 4.1.3 of Soler et al. (2013b).

The four modes with zero real component of eigenfrequency are called “entropy
modes” by Soler et al. (2013a), and “vortex modes” by Zaqarashvili et al.
(2011) based on the natures of zero-frequency modes in 1F plasmas. However, in
the 2F context, we prefer the usage introduced above. Only the isobaric mode
has any thermal signature, and so could reasonably be called an “entropy”
mode, but in view of its pressure balance character, we adopt “isobaric” as a
more precise descriptor. The flow differential modes are most notable for
their discrepant velocities, and so are best named accordingly, rather than
with reference to the unrelated vorticity
$i\,\boldsymbol{k}{\boldsymbol{\times}}{\boldsymbol{v}}$.

###  2.3 Cutoffs

From the exact expression for the Alfvén eigenfrequency, it can be shown that
there are Alfvén cutoff wavenumbers at

$\begin{split}k_{\text{A,c}}^{\pm}&=\frac{\nu_{\mathrm{nc}}\,|\sec\theta|}{8(\chi+1)a}\,\sqrt{\chi^{2}+20\chi-8\mp(\chi-8)^{3/2}\chi^{1/2}}\\\\[4.0pt]
&=\frac{2\sqrt{2}\,\nu_{\mathrm{nc}}\,|\sec\theta|}{a}\,\frac{1+\chi}{\sqrt{\chi^{2}+20\chi-8\pm(\chi-8)^{3/2}\chi^{1/2}}}\\\\[8.0pt]
&=\begin{cases}\displaystyle\frac{2\nu_{\mathrm{nc}}|\sec\theta|}{a}\left(1-\frac{1}{\chi}+\mathcal{O}\left(\chi^{-3}\right)\right)&\text{for
``$+$'' sign,}\\\\[12.0pt]
\displaystyle\frac{\nu_{\mathrm{nc}}|\sec\theta|\sqrt{\chi}}{2a}\left(1+\frac{3}{2\chi}+\mathcal{O}\left(\chi^{-2}\right)\right)&\text{for
``$-$'' sign,}\end{cases}\quad\text{for $\chi\gg 1$}.\end{split}$ (4)

The second line here accords with Equation (20) of Soler et al. (2013b) (for
the case $\theta=0$), noting that their $c_{\text{A}}$ is our
$a_{\text{c}}=B/\sqrt{\mu\rho_{\mathrm{c}}}=a\sqrt{1+\chi}$, the charges-only
Alfvén speed. As they note, there is a real cutoff wavenumber interval only if
$\chi>8$. For reference, although $\chi\gg 1$ in the low solar atmosphere, it
is only about 1 at the top of the chromosphere. The Alfvén eigenfrequency is
pure imaginary (evanescent) for $k_{\text{A,c}}^{+}\leqslant k\leqslant
k_{\text{A,c}}^{-}$. The slow wave is subject to similar cutoffs (Soler et
al., 2013a). However, these Alfvén and slow cutoffs are at very high
wavenumbers and associated frequencies in the solar context, so are of limited
practical importance. They are illustrated in Figure 1.

###  2.4 Section Summary

In summary, the full solution of the wave equations takes the form

$\boldsymbol{X}(t)=\sum_{m=1}^{10}C_{m}\boldsymbol{X}_{\\!m}\,e^{-i\,\omega_{m}t},$
(5)

where the $\omega_{j}$ are the complex eigenfrequencies and the
$\boldsymbol{X}_{\\!m}$ are the corresponding eigenvectors. In matrix form,
this is

$\boldsymbol{X}(t)={\mathrm{P}}e^{-i\,\Omega\,t}{\boldsymbol{C}}$ (6)

where $\Omega=\mathop{\rm diag}[\omega_{1},\ldots,\omega_{10}]$,
${\boldsymbol{C}}=(C_{1},\ldots,C_{10})^{T}$, and ${\mathrm{P}}$ is the matrix
of eigenvectors of ${\mathrm{M}}$ arranged as columns consistently with the
eigenvalues.

In all cases of interest, the three flow differential modes (arbitrarily
$m=5$, 6 and 7) have large negative imaginary eigenfrequencies, and hence
disappear almost instantaneously. For all intents and purposes therefore, only
the remaining seven eigenvalues and their eigenvectors are significant. We saw
analytically in Equations (3) and shall see graphically in Figure 2 that these
exhibit only minor drift between neutrals and charges.

##  3 Energy and Flux

It is conventional to introduce quadratic expressions for wave-energy density
and wave-energy-density-flux in linear wave studies. Linear energy
contributions such as
${\boldsymbol{B}}{\boldsymbol{\cdot}}{\boldsymbol{b}}/\mu$ for the magnetic
energy average to zero over a period and are ignored. Extending the standard
Eckart (1960) process to the 2F equations (2), it is easily shown (Cally &
Gómez-Míguez, 2023) that

$\frac{\partial E}{\partial
t}+\mbox{\boldmath$\nabla$}{\boldsymbol{\cdot}}\hskip
1.5pt{\boldsymbol{F}}=-\alpha_{\mathrm{cn}}\,\left|{\boldsymbol{v}_{\mathrm{n}}}-{\boldsymbol{v}_{\mathrm{c}}}\right|^{2},$
(7)

where
$E=\rho_{\mathrm{n}}|{\boldsymbol{v}_{\mathrm{n}}}|^{2}/2+p_{\mathrm{n}}^{2}/(2\rho_{\mathrm{n}}c_{\mathrm{n}}^{2})+\rho_{\mathrm{c}}|{\boldsymbol{v}_{\mathrm{c}}}|^{2}/2+p_{\mathrm{c}}^{2}/(2\rho_{\mathrm{c}}c_{\mathrm{c}}^{2})+|{\boldsymbol{b}}|^{2}/2\mu$
may be interpreted as a quadratic energy density, and
${\boldsymbol{F}}=p_{\mathrm{n}}{\boldsymbol{v}_{\mathrm{n}}}+p_{\mathrm{c}}{\boldsymbol{v}_{\mathrm{c}}}-\mu^{-1}\left({\boldsymbol{v}_{\mathrm{c}}}{\boldsymbol{\times}}{\boldsymbol{B}}\right){\boldsymbol{\times}}{\boldsymbol{b}}$
is the corresponding energy density flux. The right hand side represents
collisional loss. The energy density consists of respectively the kinetic
energy of the neutrals, the thermal energy of the neutrals, the kinetic energy
of the charges, the thermal energy of the charges and the magnetic energy. The
flux consists of the rate of working by the neutrals and charges pressure
perturbations and the Poynting flux.

When a complex representation of the variables is in use, the quadratic terms
in $E$ are replaced by their absolute values, so

$E={\textstyle\frac{1}{2}}\rho_{\mathrm{n}}|{\boldsymbol{v}_{\mathrm{n}}}|^{2}+\frac{|p_{\mathrm{n}}|^{2}}{2\rho_{\mathrm{n}}c_{\mathrm{n}}^{2}}+{\textstyle\frac{1}{2}}\rho_{\mathrm{c}}|{\boldsymbol{v}_{\mathrm{c}}}|^{2}+\frac{|p_{\mathrm{c}}|^{2}}{2\rho_{\mathrm{c}}c_{\mathrm{c}}^{2}}+\frac{|{\boldsymbol{b}}|^{2}}{2\mu},$
(8)

and the flux becomes

${\boldsymbol{F}}=\mathop{\rm
Re}\nolimits\left(p_{\mathrm{n}}^{*}{\boldsymbol{v}_{\mathrm{n}}}+p_{\mathrm{c}}^{*}{\boldsymbol{v}_{\mathrm{c}}}-\frac{1}{\mu}\left({\boldsymbol{v}_{\mathrm{c}}}{\boldsymbol{\times}}{\boldsymbol{B}}\right){\boldsymbol{\times}}{\boldsymbol{b}}^{*}\right).$
(9)

Any particular initialization of oscillations corresponds to a particular
partitioning of energy into the ten eigenmodes.

The quadratic wave energy density may be expressed as an Hermitian form in
terms of the spectral coefficients ${\boldsymbol{C}}$:

$E=\rho_{\mathrm{n}}\,{\boldsymbol{C}}^{\dagger}\Psi(t){\boldsymbol{C}}.$ (10)

Similarly, the vertical wave energy flux is

$F_{z}={\textstyle\frac{1}{2}}\rho_{\mathrm{n}}{\boldsymbol{C}}^{\dagger}\Xi(t)\,{\boldsymbol{C}}.$
(11)

The Hermitian matrices $\Psi$ and $\Xi$ are described in Appendix C. Of
course, the flux in any other direction can be calculated similarly, but
$F_{z}$ will serve to illustrate the point.

##  4 Initializing Alfvén Waves on the Charges Alone

Because of its intrinsic interest as well as its relative simplicity, we now
look at the Alfvén case, polarized in the $y$-direction, but initiate it in
three contrasting ways: (i) with an initial velocity $v_{\text{c}}$ only; (ii)
with an initial magnetic perturbation $b_{y}$ only; and (iii) with a fully
formed Alfvén wave on the charges alone. By way of exposition, a further
whole-plasma kinetic excitation is also briefly discussed: (iv) initiation via
$v_{\text{n}}=v_{\text{c}}$ only.

  1. (i)

Ignoring terms of order $\nu_{\mathrm{nc}}^{-1}$ in the eigenvectors set out
above, the velocity-only initiation proposed by Vranjes et al. (2008) is
constructed uniquely from our eigenvectors as

$\color[rgb]{0,0,0}\boldsymbol{X}(0)=(0,0,0,0,0,1,0,0,0,0)={\textstyle\frac{1}{2}}(A_{-}-A_{+}+2f_{y})/(\chi+1).$
(12a)

This consists initially of two oppositely directed full-plasma Alfvén waves
and the $y$-polarized flow differential mode that decays almost instantly.
Based on Equations (C1) and (C2), the total energy density at time $t=0$
associated with unit velocity is
$E(0)=\rho_{\mathrm{n}}/(2\chi)=\rho_{\mathrm{c}}/2$, as expected. Almost
instantly, the $f_{y}$ component decays due to collisions, leaving only the
two oppositely directed Alfvén waves, each with energy density
$\color[rgb]{0,0,0}E_{\pm}=\frac{1}{8}\rho_{\mathrm{c}}/(\chi+1)$, which is
negligible in comparison for $\chi\gg 1$. The corresponding wave energy fluxes
$a\,E_{+}$ and $-a\,E_{-}$ directed along the equilibrium magnetic field are
similarly small. This is very inefficient.

  2. (ii)

Alternatively, initializing the wave with a magnetic perturbation only,

$\boldsymbol{X}(0)=(0,0,0,0,0,0,0,0,0,1/\sqrt{\chi+1})=(A_{+}+A_{-})/(2\sqrt{\chi+1})$
(12b)

similarly injects energy $E(0)=\rho_{\mathrm{c}}/2$. However now the modal
decomposition contains no decaying flow differential component to leading
order. The two energy densities are $E_{\pm}=\rho_{\mathrm{c}}/4$, equally
splitting the energy flux $\pm\,\rho_{\mathrm{c}}a/4$ in the two directions.
This is very efficient.

  3. (iii)

Next, we suppose that a fully developed Alfvén wave is placed on the charges
only at $t=0$, as if the charges and neutrals were collisionally disconnected
on $t<0$ with collisions turning on discontinuously at $t=0$. Then

$\begin{split}\boldsymbol{X}(0)&=\frac{1}{\sqrt{2}}(0,0,0,0,0,-1,0,0,0,(\chi+1)^{-1/2})\\\
&=\left\\{{\textstyle\frac{1}{2}}\left[1+(\chi+1)^{-1/2}\right]A_{+}+{\textstyle\frac{1}{2}}\left[1-(\chi+1)^{-1/2}\right]A_{-}-(\chi+1)^{-1/2}f_{y}\right\\}/(2(\chi+1))^{1/2},\end{split}$
(12c)

with total energy $E(0)=\rho_{\mathrm{c}}/2$. Once again, the flow
differential mode $f_{y}$ vanishes rapidly, this time leaving slightly
unbalanced Alfvén waves propagating in the two field-aligned directions. The
respective energy densities are
$E_{\pm}=\rho_{\mathrm{c}}\left(\sqrt{\chi+1}\pm 1\right)^{2}/8(\chi+1)$ for
the Alfvén waves and $E_{y}=\rho_{\mathrm{c}}\chi/4(\chi+1)$ for the flow
differential mode. In the $\chi\gg 1$ regime, these are asymptotically
$E_{\pm}\sim\rho_{\mathrm{c}}/8$, with a small excess at higher order in the
prograde Alfvén wave compared to the retrograde one, and
$E_{y}\sim\rho_{\mathrm{c}}/4$. So, in effect, this charges-only Alfvén
excitation rapidly loses half its energy to collisions and initiates
oppositely directed Alfvén waves each with a quarter of the total energy, and
carrying fluxes $\pm\rho_{\mathrm{c}}a/8$ in each direction. Despite the
imposed directionality of the original Alfvén wave on the charges, the
resulting Alfvén waves on the full plasma are nearly balanced in the two
directions. This case is 50% efficient in generating Alfvén waves.

  4. (iv)

Finally, we mention the whole-plasma excitation case where the species
velocities $v_{\text{n}}=v_{\text{c}}$ are initiated together, as envisaged
explicitly in Section 6 of Soler et al. (2013b) and implicitly by Tsap et al.
(2011). Then

$\boldsymbol{X}(0)=(0,-1/\sqrt{\chi+1},0,0,0,-1/\sqrt{\chi+1},0,0,0,0)=\left(A_{+}-A_{-}\right)/(2\sqrt{\chi+1}).$
(12d)

This again places no energy on the flow differential modes, with almost
identical consequences to case (ii): the oppositely directed Alfvén waves are
generated equally, with no collisional energy loss.

Overall, generation of Alfvén waves via the charges alone can range from
highly inefficient (case i) to highly efficient (case ii) and in between (case
iii). The distinguishing feature is the energy initially placed in the flow
differential mode, which is rapidly damped and its energy lost to
thermalization. Cases (ii) and (iv) suffer no such losses to leading order in
the inverse collision frequency, and so are highly efficient.

Of the three scenarios discussed, case (i), excitation via $v_{\text{c}}$
only, is the least plausible as well as the least efficient. It is difficult
to envisage a mechanism that would inject energy into the charges plasma
velocity only and not the magnetic field $b_{y}$. Specifically, a purely
mechanical mechanism could not distinguish between the charges and the
neutrals, and an electrical impulse would generate a magnetic perturbation via
Faraday induction.

A similar analytical analysis could be performed for the magneto-acoustic
waves, though with greater algebraic complexity. However, we choose to proceed
numerically. In the following section Alfvén as well as the analogous magneto-
acoustic wave excitations are explored via computationally-derived eigenvalues
and eigenvectors at three representative levels in the low solar atmosphere
for charges-only initiations akin to case (iii) above.

##  5 Numerical Results: Three Heights

Consider three heights in the solar atmosphere; lower photosphere ($h=0$ km),
upper photosphere ($h=250$ km) and temperature minimum ($h=560$ km), with
representative sound speed $c$, Alfvén speed $a$, ionization ratio $\chi$ and
collision frequency $\nu_{\mathrm{nc}}$ set out in Table 1. The complex
eigenfrequencies $\omega_{\text{r}}+i\,\omega_{\text{i}}$ of the slow, Alfvén
and fast waves are plotted against $k$ in Figure 1 for the temperature minimum
$h=560$ km model. The first Alfvén cutoff, $k_{\text{A,c}}^{+}$, and its
similar slow counterpart, are apparent. They accurately match the asymptotic
formulae given in Equations (3) until the cutoffs start to intrude.

Table 1: Atmospheric Models and Wave Eigenfrequencies

| $h$ | $c$ | $a$ | $\chi$ | $\nu_{\mathrm{nc}}$ | flow $\tau$ | isobaric $\tau$ | slow | Alfvén | fast
---|---|---|---|---|---|---|---|---|---|---
| (km) | ($\rm km\,s^{-1}$) | ($\rm km\,s^{-1}$) |  | ($\rm s^{-1}$) | (s) | (s) | Freq (Hz) |  $\tau$ (s) | Freq (Hz) |  $\tau$ (s) | Freq (Hz) |  $\tau$ (s)
photospheric base | 0 | 10 | 0.63 | 1400 | $7\times 10^{5}$ |  | $1.2\times 10^{5}$ | 0.5453 | $9\times 10^{4}$ | 0.5456 | $1.2\times 10^{5}$ | 10 | $6\times 10^{7}$
upper photosphere | 250 | 8.6 | 1.4 | 11000 | 1600 |  | $3.0\times 10^{4}$ | 1.208 | 416 | 1.212 | 551 | 8.6 | $6\times 10^{4}$
temperature minimum | 560 | 8.2 | 5.5 | 10500 | 1000 |  | $2.0\times 10^{3}$ | 4.42 | 2.0 | 4.76 | 2.2 | 8.8 | 9.2

Note. — Representative atmospheric parameters at the photospheric base $h=0$,
upper photosphere $h=250$ km and temperature minimum $h=560$ km drawn roughly
from Figures 1 and 2 of Cally & Gómez-Míguez (2023), adapted for a hydrogen
atmosphere from Model C7 of Avrett & Loeser (2008) with 100 G base magnetic
field and a 600 km magnetic scale height fall-off. In the second block,
frequencies $\mathop{\rm Re}\nolimits\omega/2\pi$ (Hz) and decay times
$\tau=1/|\mathop{\rm Im}\nolimits\omega|$ of the slow, Alfvén and fast waves
in these atmospheres with $k=k_{1}$ (1 km wavelength) and $\theta=30^{\circ}$
are listed, as well as the decay times of the
$|{\boldsymbol{v}_{\mathrm{n}}}-{\boldsymbol{v}_{\mathrm{c}}}|$ neutrals-
charges flow differential and the isobaric mode.

For concreteness, we primarily examine these models with the fiducial
wavenumber $k=k_{1}=2\pi/10^{3}=0.0063$ $\rm rad\,m^{-1}$, corresponding to a
wavelength of 1 km, and $\theta=30^{\circ}$. A non-zero propagation angle
$\theta$ removes any ambiguity between Alfvén and slow waves. With these
values, the corresponding numerically calculated frequencies
$\omega_{\text{r}}/2\pi$ (Hz) and decay times $\tau=1/\omega_{\text{i}}$ are
also set out in Table 1. They scale with $k$ as set out in Equations (3).

The slow and Alfvén modes are not followed beyond their cutoffs. The decay
rates of the acoustically dominated fast wave are well below those of the
magnetically dominated slow and Alfvén waves.

Figure 1: Phase speed $\omega_{\text{r}}/k$ ($\rm km\,s^{-1}$, log-linear,
left) and decay rate $|\omega_{\text{i}}|$ ($\rm s^{-1}$, log-log, right)
plotted against wavenumber $k$ (rad $\rm m^{-1}$) for the temperature minimum
model $h=560$ km with $\theta=30^{\circ}$. For comparison, the specific
wavenumber $k_{1}=0.0063$ referenced in Table 1 and various examples below is
comfortably short of the cutoffs and therefore in the asymptotic regime of
Equations (3). Blue: slow wave; red long-dashed: Alfvén wave; green short-
dashed: fast wave. The vertical line indicates the Alfvén cutoff
$k_{\text{A,c}}^{+}$.

Figure 2: Energies in the ten components
$(u_{\mathrm{n}},v_{\mathrm{n}},w_{\mathrm{n}},\psi_{\mathrm{n}},u_{\mathrm{c}},v_{\mathrm{c}},w_{\mathrm{c}},\psi_{\mathrm{c}},\alpha_{\perp},\alpha_{y})$
of the eigenmodes (log scale)

for the three models of Table 1, $h=0$, 250 and 560 km respectively for the
first, second and third columns. In each case, $k=k_{1}$ and
$\theta=30^{\circ}$. Total energy is normalized to 1, and each mode is
labelled with its eigenfrequency. The rows correspond to the different mode
types. From top to bottom: isobaric mode; the three flow differential modes;
slow; Alfvén; fast. The colours correspond to neutrals-kinetic (light fawn);
neutrals-thermal (dark fawn); charges-kinetic (red); charges-thermal (red-
brown); and magnetic (blue).

In the spirit of Vranjes et al. (2008), the cases where either an Alfvén, slow
or fast wave is initiated on the _decoupled_ charges alone are now explored.
This is an extension of Scenario (iii) of Section 4, and represents a sort of
thought experiment where we suppose the neutrals and charges are collisionally
decoupled on $t<0$, and one of these MHD modes is placed on the charges alone.
It is of interest to determine how these charges-only eigenmodes project onto
the 10D eigenspace of the coupled system. Collisions are then turned on at
$t=0$.

Figure 2 physically characterizes the eigenmodes by presenting the energy
distributions associated with each of the ten components of $\boldsymbol{X}$
for each of seven eigenmodes of the coupled system (rows) for each $h$
(columns) with $k=k_{1}$ and $\theta=30^{\circ}$ (the oppositely directed
slow, Alfvén and fast modes are omitted as they are the same as for the
forward-directed cases). Noting that the ordinate is presented
logarithmically, it is seen that:

  * •

The isobaric mode resides predominantly on the charges. The pressure
perturbations in the neutral and charged fluids cancel out,
$p_{\mathrm{n}}+p_{\mathrm{c}}=0$, allowing a static atmosphere in the strong
coupling regime. Therefore, the respective thermal energy densities are in the
ratio $2:\chi$, explaining the predominance of the latter at large $\chi$.
Kinetic and magnetic energies are negligible in comparison to the therml
energy in this mode.

  * •

In the three flow differential modes, there is a clear difference between
${\boldsymbol{v}_{\mathrm{n}}}$ and ${\boldsymbol{v}_{\mathrm{c}}}$. If they
were identical, their respective kinetic energies would be in the ratio
$\chi:1$, with more energy on the neutrals, which is certainly not the case.
This situation is unsustainable in the presence of inter-species collisions,
and hence produces the extremely rapid decay indicated by the
eigenfrequencies.

  * •

The slow modes are dominated by magnetic and neutrals-kinetic energies in the
$x$-$z$ plane, with the neutral and charges velocities almost matching.

  * •

The Alfvén modes are polarized in the $y$ direction only, with equipartition
between kinetic and magnetic energies. The neutrals and the charges velocities
match closely, and hence the neutrals kinetic energy dominates the charges.

  * •

The fast modes are similar to the slow modes, except that thermal energy plays
a much larger role, as is to be expected in predominantly acoustic waves.

Figure 3: Dimensionless velocity differential (logarithmic scale)
$|{\boldsymbol{v}_{\mathrm{c}}}-{\boldsymbol{v}_{\mathrm{n}}}|\,\sqrt{{\textstyle\frac{1}{2}}\rho_{\mathrm{c}}/E}$
associated with each eigenmode (as labelled) for the models at $h=0$ km
(blue); $h=250$ km (yellow); and $h=560$ km (black), all with $k=k_{1}$ and
$\theta=30^{\circ}$. The flow differential
${\boldsymbol{v}_{\mathrm{c}}}-{\boldsymbol{v}_{\mathrm{n}}}$ is dominated by
the three extremely rapidly decaying flow differential modes, with the MHD and
isobaric modes retaining only small remnant drifts that vanish as
$\nu_{\mathrm{nc}}\to\infty$.

Figure 3 shows the distribution of non-dimensionalized
$|{\boldsymbol{v}_{\mathrm{c}}}-{\boldsymbol{v}_{\mathrm{n}}}|$ over the ten
eigenmodes for each $h$. It is seen that this is dominated by the flow
differential modes, with orders of magnitude smaller drift between neutrals
and charges on all the other modes. This shows that when these modes are
essentially quenched over a few nanoseconds, very little species drift
remains, and the wave behaves largely as 1F MHD waves, though with diffusion
operating over times of order $\nu_{\mathrm{nc}}/\omega_{\mathrm{r}}^{2}$.
From Figure 2, nearly all energy in the flow differential modes resides in the
kinetic energy of the charges flow, explaining why
$\sqrt{{\textstyle\frac{1}{2}}\rho_{\mathrm{c}}|{\boldsymbol{v}_{\mathrm{c}}}-{\boldsymbol{v}_{\mathrm{n}}}|^{2}/E}\approx
1$ for those modes.

###  5.1 Projection onto Eigenmodes

Figure 4 surveys the modal energy distributions of all three initialization
cases at $t=0$ and $t=\delta t$ in the three low-atmosphere models ($\delta
t=2$ ns, 10 ns and 150 ns respectively for $h=0$ km, 250 km and 560 km). We
see that one flow differential mode accounts for half the energy at $t=0$, but
this has essentially disappeared by $t=\delta t$. In essence, the oscillations
can be thought to start from this few-nanoseconds state in which the flow
differential eigenmodes have been suppressed.

The modal fluxes in the MHD waves calculated using Equation (9) accord
perfectly with the energies multiplied by their respective ideal MHD group
velocities,

$\frac{\partial\omega}{\partial
k_{z}}=\begin{cases}\displaystyle\pm\frac{\cos\theta\left(\Delta\left(a^{2}+c^{2}\right)+2a^{2}c^{2}\cos^{2}\theta-a^{4}-c^{4}\right)}{\sqrt{2}\,\Delta\sqrt{a^{2}+c^{2}-\Delta}},&\text{slow}\\\\[4.0pt]
\displaystyle\pm a,&\text{Alfv\'{e}n}\\\\[4.0pt]
\displaystyle\pm\frac{\cos\theta\left(\Delta\left(a^{2}+c^{2}\right)-2a^{2}c^{2}\cos^{2}\theta+a^{4}+c^{4}\right)}{\sqrt{2}\,\Delta\sqrt{a^{2}+c^{2}+\Delta}},&\text{fast}\\\
\end{cases},$ (13)

i.e., ${\boldsymbol{F}}=E\,\partial\omega/\partial\boldsymbol{k}$, and so need
not be presented explicitly. However, this serves as a check on the numerics.

An alternative and more realistic scenario is that some driver operates solely
on the charges over seconds or minutes. During this period, the near-
instantaneous decay of the flow differential modes keeps the neutrals and
charges fluid velocities strongly coupled, thereby resulting in energy
distributions across the components of $\boldsymbol{X}$ as set out in the
first and the last three rows of Figure 2. Across the height range, it is seen
that:

  * •

Initiating with the charges-restricted Alfvén wave (yellow) loses half its
energy to the flow differential mode, which quickly vanishes, and about a
quarter each to the upgoing and downgoing Alfvén modes.

  * •

Slow wave (acoustic on the charges) initiation (brown) also loses half its
energy to the flow differential, with most of the remainder apportioned to the
isobaric mode. Very little ends up in the travelling waves.

  * •

Fast wave initiation (blue) loses half its energy to the flow differential
with around $20-25$% going to each of the upgoing and downgoing slow waves. A
small amount also goes to the two fast waves at $h=560$ km. Recalling that the
fast wave on the charges is primarily magnetic, due to the large
$a_{\mathrm{c}}$, it was to be expected that it would primarily drive the slow
waves (also magnetic) of the coupled plasma.

Overall, Alfvén and slow waves (the two magnetic waves) are quite efficiently
excited and are near-symmetric in direction, whilst generation of the
acoustically dominated fast wave is very inefficient for all three drivers.

Initiating with $b_{\perp}$ or $b_{y}$ alone (not shown), and no velocity,
just splits the energy into two equal but oppositely directed slow or Alfvén
modes respectively, _à la_ d’Alembert. This was described analytically for
Alfvén waves in Section 4. There is no flow differential at any stage, and
hence no energy loss, so this initiation is particularly efficient. Half of
the initialization energy propagates in each direction.

Figure 4: Energy in each of the eigenmodes (negatively directed fast,
negatively directed Alfvén, …, positively directed Alfvén, positively directed
fast, as labelled) at times $t=0$ (left) and $t=\delta t>0$ (right) for model
$h=0$ (top row, $\delta t=2$ ns), $h=250$ km (middle row, $\delta t=10$ ns)
and $h=560$ km (bottom, $\delta t=150$ ns) with initiation by Alfvén (yellow),
slow (brown) and fast (blue) positively directed eigenmodes of the charged
fluid only. Total initial energies are normalized to 1. The three flow
differential modes have almost disappeared in just nanoseconds, whilst the
others have not changed.

Although there is quite efficient full-plasma wave generation in this charges-
only wave initiation scenario, Figure 4 explains why _net_ flux is still very
small. It is _not_ because negligible Alfvén (or slow or fast) wave is
excited. Instead, despite the initial state being an upward (positive
$z$-direction) MHD eigenmode on the charges-only plasma, it projects onto
almost equal amplitudes of positively and negatively directed MHD waves in
each case. Therefore, when calculating fluxes, these nearly cancel, though
each carries up to a quarter of the initial perturbation energy. The initial
but fast-disappearing energy in the flow differential modes was never going
anywhere anyway, since it had zero $\omega_{\mathrm{r}}$, and in any case took
only half the energy with it.

A real wave excitation region will of course be of finite extent. In most
circumstances, it will generate near-equal fluxes of upward and downward
propagating MHD waves. In the case of initiation on the charges only,
resulting waves on the full plasma carry up to 25% of the perturbation energy
in each direction for the Alfvén wave, a little less for slow waves initiated
by fast waves, and very little for fast waves. Net flux is therefore very
small. However, once they emerge above or below the excitation region, their
uni-directional Alfvén and slow wave fluxes will be manifest and substantial.
It is not correct to say that no significant flux is produced in the Vranjes
et al. scenario.

##  6 Discussion

For specified wavevector $\boldsymbol{k}$, the linearized two-fluid
collisionally-coupled MHD equations for a weakly ionized plasma imply a ten-
dimensional space of solutions, consisting of the familiar six MHD wave modes
– slow, Alfvén and fast propagating either forward or backward – as well as a
very slowly decaying stationary isobaric mode and three extremely rapidly
decaying (nanoseconds) flow differential mode. The role of these flow modes is
to support large velocity differentials between neutrals and charges. They
alone carry this responsibility, as the remaining seven eigenmodes
intrinsically support only very small differentials.

With hindsight, it could be no other way. It is well known that collisions
discharge any significant difference between ${\boldsymbol{v}_{\mathrm{c}}}$
and ${\boldsymbol{v}_{\mathrm{n}}}$ much more rapidly than any other timescale
associated with 2F MHD waves, so they must have their own dedicated eigenmodes
with eigenfrequencies that have large negative imaginary parts. When
$\omega\ll\nu_{\mathrm{nc}}$, the remaining 2F MHD waves differ only slightly
from the standard 1F description, $\nu_{\mathrm{nc}}\to\infty$, in which
${\boldsymbol{v}_{\mathrm{c}}}={\boldsymbol{v}_{\mathrm{n}}}$ exactly. The
asymptotic eigenfrequency formulae presented in Equations (3) describe these
processes in detail. For the Alfvén mode specifically, addressed in detail in
Section 4, the asymptotic eigenvector corresponding to eigenvalue
$\omega_{\mathrm{A}}$ is given by Equation (3b), so
$\left(v_{\mathrm{c}}-v_{\mathrm{n}}\right)/|v_{\mathrm{c}}|\sim
i\,a\,k_{z}/\nu_{\mathrm{nc}}$ as $\nu_{\mathrm{nc}}\to\infty$, showing indeed
that the interspecies drift scales as $\nu_{\mathrm{nc}}^{-1}$ and vanishes in
the high-collision limit.

In the cases examined numerically in Section 5, wave initiation in the form of
charges-specific MHD waves results in almost exactly one half of the initial
energy residing on the flow differential modes (see Figure 4), which
disappears immediately. The remaining MHD eigenmodes are all roughly
symmetrically present in energy, i.e., the forward and backward directed modes
of each species have nearly equal energies. This is not an accident. If it
were not so, the remaining net flux would not be small, contrary to the very
small velocities implied by the vanishing of the flow differential modes. The
result is also consistent with our analytical analysis of Section 4.

As in case (ii) of Section 4, an oscillation excited purely via a magnetic
perturbation produces no significant flow differential modes and hence
essentially no energy loss. The flux still splits into two oppositely directed
parts though, now 50:50. This is a very efficient excitation mechanism.

On the other hand, it is also possible in principle to excite the flow
differential modes solely, in which case all energy is immediately lost.
However, it is difficult to imagine a process that would do this in practice,
as it would involve an initial state consisting of one or more of the flow
differential eigenmodes given asymptotically in Equation (3i).

In summary, the ten-dimensional spectral decomposition of the 2F MHD Equations
(2) gives considerable insight into how excitation of the charged fluid alone
rapidly suppresses interspecies drift via the rapid decay of the flow
differential modes. The remaining energy, dependent on the specific initial
state, typically finds itself in long-lived MHD waves. This indicates that it
is indeed plausible to launch substantial wave flux, especially Alfvén wave
flux, upward from a weakly ionized photosphere even if the excitation
mechanism only directly accesses the charges.

Of course, if excitation is mechanical and equally drives both species in
concert, there was never an issue in the first place.

##  Appendix A Matrix M and its Characteristic Polynomial

The fundamental coefficient matrix

${\mathrm{M}}=\\\
\left(\begin{array}[]{cccc|cccccc}-\nu_{\mathrm{nc}}&0&0&-ik_{x}&\nu_{\mathrm{nc}}&0&0&0&0&0\\\
0&-\nu_{\mathrm{nc}}&0&0&0&\nu_{\mathrm{nc}}&0&0&0&0\\\
0&0&-\nu_{\mathrm{nc}}&-ik_{z}&0&0&\nu_{\mathrm{nc}}&0&0&0\\\
-\frac{ic^{2}(\chi+1)k_{x}}{\chi+2}&0&-\frac{ic^{2}(\chi+1)k_{z}}{\chi+2}&0&0&0&0&0&0&0\\\\[4.0pt]
\hline\cr\chi\nu_{\mathrm{nc}}&0&0&0&-\chi\nu_{\mathrm{nc}}&0&0&-ik_{x}&ia(\chi+1)k&0\\\
0&\chi\nu_{\mathrm{nc}}&0&0&0&-\chi\nu_{\mathrm{nc}}&0&0&0&ia(\chi+1)k_{z}\\\
0&0&\chi\nu_{\mathrm{nc}}&0&0&0&-\chi\nu_{\mathrm{nc}}&-ik_{z}&0&0\\\
0&0&0&0&-\frac{2ic^{2}(\chi+1)k_{x}}{\chi+2}&0&-\frac{2ic^{2}(\chi+1)k_{z}}{\chi+2}&0&0&0\\\
0&0&0&0&iak&0&0&0&0&0\\\ 0&0&0&0&0&iak_{z}&0&0&0&0\\\ \end{array}\right),$
(A1)

where $k_{x}=k\sin\theta$ and $k_{z}=k\cos\theta$. Vertical and horizontal
dividers have been included to accentuate the neutrals (top left block) and
charges (bottom right) submatrices and explicate their collisional couplings
(top right and bottom left).

Since the $y$ polarization is decoupled from $x$-$z$, ${\mathrm{M}}$ can be
broken down into separate third and seventh order matrices (it can be
expressed in block-diagonal form under the appropriate reordering of rows and
columns). Hence, the tenth order characteristic polynomial (dispersion
function) of ${\mathrm{M}}$ may be factored into a cubic

$\omega^{2}\left(\omega+i\,(\chi+1)\nu_{\mathrm{nc}}\right)-a^{2}k^{2}(\chi+1)\cos^{2}\theta\left(\omega+i\,\nu_{\mathrm{nc}}\right)=0,$
(A2)

which captures the two Alfvén modes and one flow differential mode (those
oscillations polarized in the $y$-direction), and an equation of the seventh
order (see also Eq. (57) of Zaqarashvili et al., 2011) that contains the four
magneto-acoustic modes, two flow differential modes and the isobaric mode:

$\omega\left(\omega^{4}-(a^{2}+c^{2})k^{2}\omega^{2}+a^{2}c^{2}k^{4}\cos^{2}\theta\right)=\\\
\omega\,\tau^{2}\left(\frac{\omega^{6}}{(\chi+1)^{2}}-\frac{k^{2}\omega^{4}\left(a^{2}(\chi+2)+3c^{2}\right)}{(\chi+1)(\chi+2)}+\frac{c^{2}k^{4}\omega^{2}\left(2\left(a^{2}(\chi+2)+c^{2}\right)+a^{2}(\chi+2)\cos
2\theta\right)}{(\chi+2)^{2}}-\frac{2a^{2}c^{4}k^{6}(\chi+1)\cos^{2}\theta}{(\chi+2)^{2}}\right)\\\
+i\,\tau\left(\frac{2\,\omega^{6}}{\chi+1}-\frac{k^{2}\omega^{4}\left(a^{2}(\chi+2)^{2}+c^{2}(4\chi+5)\right)}{(\chi+1)(\chi+2)}+\frac{c^{2}k^{4}\omega^{2}\left(2a^{2}(\chi+2)\cos
2\theta+a^{2}(\chi+2)(\chi+3)+2c^{2}(\chi+1)\right)}{(\chi+2)^{2}}\right.\\\
\left.-\frac{2a^{2}c^{4}k^{6}(\chi+1)\cos^{2}\theta}{(\chi+2)^{2}}\right),$
(A3)

where $\tau=1/\nu_{\mathrm{nc}}$ is the neutrals-charges collision timescale.
In the fully coupled limit $\tau\to 0$, this dispersion relation reduces to
that of the classic 1F magneto-acoustic modes and the stationary isobaric mode
$\omega=0$. The full tenth order formulation is retained here for unity of
exposition.

##  Appendix B Eigensystem Perturbation Analysis

Although the eigenvalues and eigenvectors of ${\mathrm{M}}$ may be found
numerically, the analytic large-$\nu_{\mathrm{nc}}$ asymptotic formulae set
out in Equations (3) are valuable aids to understanding. The process of
developing these formulae is not straightforward though. A matrix method is
sketched here, beginning with splitting
${\mathrm{M}}={\mathrm{M}}_{0}+\nu_{\mathrm{nc}}{\mathrm{M}}_{1}$. Thus

$\left({\mathrm{M}}_{0}+\nu_{\mathrm{nc}}{\mathrm{M}}_{1}\right)\boldsymbol{X}=\lambda\boldsymbol{X}.$
(B1)

Recall that the eigenfrequencies are $\omega=i\,\lambda$.

At large $\nu_{\mathrm{nc}}$, ${\mathrm{M}}_{1}$ plays the leading role. It is
only rank-3, but is diagonalizable:
${\mathrm{M}}_{1}={\mathrm{S}}\Lambda{\mathrm{S}}^{-1}$, where the columns of
${\mathrm{S}}$ are the eigenvectors and the diagonal matrix
$\Lambda=\mathop{\rm diag}(-(\chi+1),-(\chi+1),-(\chi+1),0,0,0,0,0,0,0)$ is
made up of the eigenvalues. The non-zero entries correspond to the three flow
differential modes. (Alternatively, LU decomposition works just as well in
that it also collects these modes in only the first three rows of the upper
triangular ${\mathrm{U}}$ matrix.) Then Equation (B1) may be recast as

$\left({\mathrm{Q}}+\nu_{\mathrm{nc}}\Lambda\right)\boldsymbol{Y}=\lambda\boldsymbol{Y},$
(B2)

where ${\mathrm{Q}}={\mathrm{S}}^{-1}{\mathrm{M}}_{0}{\mathrm{S}}$ and
$\boldsymbol{Y}={\mathrm{S}}^{-1}\boldsymbol{X}$.

The required procedure then differs between the flow differential and the
other seven modes.

###  B.1 Flow Differential Modes

Substituting
$\lambda=\nu_{\mathrm{nc}}(\lambda_{0}+\nu_{\mathrm{nc}}^{-1}\lambda_{1}+\nu_{\mathrm{nc}}^{-2}\lambda_{2}+\cdots)$
and
$\boldsymbol{Y}=\boldsymbol{Y}_{\\!\\!0}+\nu_{\mathrm{nc}}^{-1}\boldsymbol{Y}_{\\!\\!1}+\nu_{\mathrm{nc}}^{-2}\boldsymbol{Y}_{\\!\\!2}+\cdots$
into Equation (B2) and equating powers of $\nu_{\mathrm{nc}}$, it is found
that

$\Lambda\boldsymbol{Y}_{\\!\\!0}=\lambda_{0}\boldsymbol{Y}_{\\!\\!0},$ (B3)

an eigenvalue equation that determines $\lambda_{0}=-(\chi+1)$ for the three
flow differential modes and the corresponding leading order eigenvectors
$\boldsymbol{Y}_{\\!\\!0}$.

At the next order,

$\left(\Lambda-\lambda_{0}{\mathrm{I}}\right)\boldsymbol{Y}_{\\!\\!1}-\lambda_{1}\boldsymbol{Y}_{\\!\\!0}=-{\mathrm{Q}}\boldsymbol{Y}_{\\!\\!0},$
(B4)

etc. This may be solved for $\lambda_{1}$ and $\boldsymbol{Y}_{\\!\\!1}$ by
writing it as an $11\times 11$ matrix equation for
$(\boldsymbol{Y}_{\\!\\!1},\lambda_{1})$ when supplemented with a
normalization condition. In practice, one component of
$\boldsymbol{Y}_{\\!\\!1}$ corresponding to a nonzero component of
$\boldsymbol{Y}_{\\!\\!0}$ is set to zero keeping that component of
$\boldsymbol{Y}$ unchanged under the perturbation. Equations (3i) result on
recovering $\boldsymbol{X}={\mathrm{S}}\boldsymbol{Y}$. The procedure may be
repeated to higher order if required.

No use is made of the remaining seven eigen-solutions, which are better
attacked as follows.

###  B.2 MHD and Isobaric Modes

Defining ${\mathrm{J}}=\mathop{\rm diag}(0,0,0,1,1,1,1,1,1,1)$ it is clear
that ${\mathrm{J}}\Lambda=0$, and Equation (B2) may be used to derive

$\left(\mathrm{R}+\nu_{\mathrm{nc}}^{-1}{\mathrm{Q}}\right)\boldsymbol{Y}=\lambda\left({\mathrm{J}}+\nu_{\mathrm{nc}}^{-1}{\mathrm{I}}\right)\boldsymbol{Y},$
(B5)

where $\mathrm{R}=\Lambda+{\mathrm{J}}{\mathrm{Q}}$.

With
$\lambda=\lambda_{0}+\nu_{\mathrm{nc}}^{-1}\lambda_{1}+\nu_{\mathrm{nc}}^{-2}\lambda_{2}+\cdots$,
and $\boldsymbol{Y}$ expanded as before, it follows that

$\mathrm{R}\boldsymbol{Y}_{\\!\\!0}=\lambda_{0}{\mathrm{J}}\boldsymbol{Y}_{\\!\\!0}.$
(B6)

This is a generalized eigenvalue equation from which the leading behaviours of
the remaining seven $\lambda$ and $\boldsymbol{Y}$ may be determined.

At the next order,

$\left(\mathrm{R}-\lambda_{1}{\mathrm{J}}\right)\boldsymbol{Y}_{\\!\\!1}-\lambda_{1}{\mathrm{J}}\boldsymbol{Y}_{\\!\\!0}=\left(\lambda_{0}{\mathrm{I}}-{\mathrm{Q}}\right)\boldsymbol{Y}_{\\!\\!0},$
(B7)

which may again be used to solve for $\lambda_{1}$ and
$\boldsymbol{Y}_{\\!\\!1}$ subject to a normalization. The remainder of
Equations (3) result.

##  Appendix C Energy and Flux in Matrix Form

In complex matrix form, we may write

$E=\rho_{\mathrm{n}}\,\boldsymbol{X}^{\dagger}\Phi\boldsymbol{X}$ (C1)

where the dagger indicates conjugate transpose and

$\Phi={\textstyle\frac{1}{2}}\mathop{\rm
diag}\left[1,1,1,\frac{\chi+2}{c^{2}(\chi+1)},\frac{1}{\chi},\frac{1}{\chi},\frac{1}{\chi},\frac{\chi+2}{2c^{2}\chi(\chi+1)},\frac{\chi+1}{\chi},\frac{\chi+1}{\chi}\right].$
(C2)

In terms of the eigenvector decomposition and coefficient vector
${\boldsymbol{C}}$ the energy density is

$E=\rho_{\mathrm{n}}\,{\boldsymbol{C}}^{\dagger}\Psi(t){\boldsymbol{C}},$ (C3)

where
$\Psi(t)=e^{i\,\Omega^{*}t}{\mathrm{P}}^{\dagger}\Phi{\mathrm{P}}e^{-i\,\Omega\,t}$,
and $\Omega=\mathop{\rm diag}[\omega_{1},\ldots,\omega_{10}]$. By
construction, $\Psi$ is Hermitian. In principle, since $\Psi$ is not diagonal,
there is cross-talk between the modes. However, in practice, the off-diagonal
contributions to $E$ appear in purely imaginary complex conjugate pairs which
do not contribute at all to the overall energy as they cancel, and in any case
they are entirely negligible in magnitude. For that reason, we may attribute
all energy to the ten individual eigenmodes and not interactions between them.

The $z$-component of flux $F_{z}$ for example may also be written as a
quadratic form,

$F_{z}={\textstyle\frac{1}{2}}\rho_{\mathrm{n}}\boldsymbol{X}^{\dagger}\Upsilon\boldsymbol{X},$
(C4)

where the real symmetric (and therefore Hermitian) matrix

$\Upsilon=\begin{pmatrix}0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\
0&0&0&1&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&-\frac{a(\chi+1)\cos\theta}{\chi}&0\\\
0&0&0&0&0&0&0&0&0&-\frac{a(\chi+1)}{\chi}\\\
0&0&0&0&0&0&0&\frac{1}{\chi}&0&0\\\ 0&0&0&0&0&0&\frac{1}{\chi}&0&0&0\\\
0&0&0&0&-\frac{a(\chi+1)\cos\theta}{\chi}&0&0&0&0&0\\\
0&0&0&0&0&-\frac{a(\chi+1)}{\chi}&0&0&0&0\\\ \end{pmatrix},$ (C5)

or in terms of ${\boldsymbol{C}}$,

$F_{z}={\textstyle\frac{1}{2}}\rho_{\mathrm{n}}{\boldsymbol{C}}^{\dagger}\Xi\,{\boldsymbol{C}},$
(C6)

where
$\Xi=e^{i\,\Omega^{*}t}\,{\mathrm{P}}^{\dagger}\Upsilon{\mathrm{P}}\,e^{-i\,\Omega\,t}$.

Unlike the energy density, the cross-talk in flux between eigenmodes can be
substantial in magnitude. However, cross-talk between MHD modes averages to
zero over time, as does the interaction between a wave mode and the isobaric
mode. The interaction with a flow differential mode of course quickly
vanishes. So only the diagonal entries, attributed to slow, Alfvén and fast
waves, contribute meaningfully.

###

The author thanks Elena Khomenko and Martin Gómez Míguez for their very useful
comments and suggestions on an initial draft of this paper.

## References

  * Alharbi et al. (2022) Alharbi, A., Ballai, I., Fedun, V., & Verth, G. 2022, MNRAS, 511, 5274, doi: 10.1093/mnras/stac444
  * Avrett & Loeser (2008) Avrett, E. H., & Loeser, R. 2008, ApJS, 175, 229, doi: 10.1086/523671
  * Cally & Gómez-Míguez (2023) Cally, P. S., & Gómez-Míguez, M. M. 2023, The Astrophysical Journal, 946, 108, doi: 10.3847/1538-4357/acbb63
  * Cranmer & van Ballegooijen (2005) Cranmer, S. R., & van Ballegooijen, A. A. 2005, ApJS, 156, 265, doi: 10.1086/426507
  * de Pontieu & Haerendel (1998) de Pontieu, B., & Haerendel, G. 1998, A&A, 338, 729
  * Eckart (1960) Eckart, C. 1960, Hydrodynamics of Oceans and Atmospheres (Oxford: Pergamon)
  * Khomenko et al. (2014) Khomenko, E., Collados, M., Díaz, A., & Vitas, N. 2014, Physics of Plasmas, 21, 092901, doi: 10.1063/1.4894106
  * McIntosh et al. (2011) McIntosh, S. W., de Pontieu, B., Carlsson, M., et al. 2011, Nature, 475, 477, doi: 10.1038/nature10235
  * Popescu Braileanu et al. (2019) Popescu Braileanu, B., Lukin, V. S., Khomenko, E., & de Vicente, Á. 2019, A&A, 630, A79, doi: 10.1051/0004-6361/201935844
  * Soler et al. (2013a) Soler, R., Carbonell, M., & Ballester, J. L. 2013a, ApJS, 209, 16, doi: 10.1088/0067-0049/209/1/16
  * Soler et al. (2013b) Soler, R., Carbonell, M., Ballester, J. L., & Terradas, J. 2013b, ApJ, 767, 171, doi: 10.1088/0004-637X/767/2/171
  * Tsap et al. (2011) Tsap, Y. T., Stepanov, A. V., & Kopylova, Y. G. 2011, Sol. Phys., 270, 205, doi: 10.1007/s11207-011-9727-4
  * Vranjes et al. (2008) Vranjes, J., Poedts, S., Pandey, B. P., & de Pontieu, B. 2008, A&A, 478, 553, doi: 10.1051/0004-6361:20078274
  * Zaqarashvili et al. (2011) Zaqarashvili, T. V., Khodachenko, M. L., & Rucker, H. O. 2011, A&A, 529, A82+, doi: 10.1051/0004-6361/201016326


11institutetext: Orange Innovation, Lannion, France 22institutetext:
Université de technologie de Compiègne

# Viewing the process of generating counterfactuals as a source of knowledge
\- Application to the Naive Bayes classifier

Vincent Lemaire 11 Nathan Le Boudec 1122
Françoise Fessant 11 Victor Guyomard 11

###### Abstract

There are now many comprehension algorithms for understanding the decisions of
a machine learning algorithm. Among these are those based on the generation of
counterfactual examples. This article proposes to view this generation process
as a source of creating a certain amount of knowledge that can be stored to be
used, later, in different ways. This process is illustrated in the additive
model and, more specifically, in the case of the naive Bayes classifier, whose
interesting properties for this purpose are shown.

##  1 Introduction

Machine learning, one of the branches of artificial intelligence, has enjoyed
many successes in recent years. The decisions made by these models are
increasingly accurate, but also increasingly complex. However, it appears that
some of these models are like black boxes: their decisions are difficult, if
not impossible, to explain [2]. This lack of explicability can lead to a
number of undesirable consequences: lack of user confidence, reduced usability
of the models, presence of biases, etc. These needs have given rise to the
field of XAI (eXplainable AI). XAI [16, 1] is a branch of artificial
intelligence that aims to make the decisions made by machine learning models
intelligible, understandable, to users.

Among XAI methods, counterfactual reasoning is a concept from psychology and
the social sciences [14]. It involves examining possible alternatives to past
events [17]. Humans often use counterfactual reasoning by imagining what would
happen if an event had not occurred, and this is exactly what counterfactual
reasoning is. Applied to artificial intelligence, the question is, for
example, “Why did the model make this decision rather than another?” or “How
would the decision have been different if a certain condition had been
changed?”.

Within the framework of counterfactual reasoning, this article proposes to
consider this generation process as a source of knowledge that can be stored
and then exploited in diverse ways. This process is illustrated in the case of
additive models and in particular in the case of the naive Bayes classifier,
whose interesting properties for this purpose will be shown.

The rest of this paper is organised as follows: section 2 presents the key
concepts used in the rest of the paper, so that this paper can be read
independently. Section 3 presents the first contribution of this paper by
showing that it is possible to find “additive trajectories” of counterfactuals
in the case of the naive Bayes classifier. The section 4 presents the second
contribution of this paper by detailing how these trajectories, this
knowledge, can be store in a database and then used. Finally, before
concluding, the section 5 illustrates, by means of an unsubscription problem,
how clustering applied to this database generates new knowledge.

Note: In the remainder of this article we will focus on supervised
classification problems where a predictive model $f$ is trained using $N$
examples, each described by a set of $d$ explanatory variables (a vector
$X=\\{X_{1},....,X_{d}\\}$, derived from a $\mathcal{X}$ distribution) so as
to predict a categorical target variable denoted $Y=\\{y_{1},...,y_{C}\\}$
derived from a $\mathcal{Y}$ distribution.

##  2 Concepts

###  2.1 Countrefactual and semi-factual example

In machine learning, a counterfactual explanation aims to explain why a
particular result was obtained by suggesting hypothetical changes in the input
features $X$ that might have led to a different prediction [12, 19]. In other
words, it identifies the factors that could have influenced a particular
outcome. Counterfactual reasoning can be defined as follows:

###### Definition 1

Let $f$ be $\mathcal{X}\mapsto\mathcal{Y}$ a machine learning model such that,
for a given individual $X$, $f(X)=\hat{y}_{i}$. In this case, a counterfactual
example is a new example $X^{\prime}$ such that $f(X^{\prime})\neq\hat{y}_{i}$
and $X\neq X^{\prime}$. If $X$ had different values for some of its
explanatory variables ($X_{1}^{\prime},X_{2}^{\prime},...$) while all the
other variables remained identical, the class $\hat{y}_{j}\neq\hat{y}_{i}$
would have been returned.

Figure 1: Illustration of a counterfactual and a semi-factual. The red dots
represent initial examples ($X$). The orange dot represents a semi-factual,
the purple dot represents a counterfactual and the white line represents the
decision boundary between the red and green classes.

Above, $\hat{y}_{i}$ is known and factual, while $\hat{y}_{j}$ is the
unexpected result that did not happen, counterfactual. Note, however, that a
change from $X$ to $X^{\prime}$ does not necessarily lead to a change in the
class prediction: this is known as a semi-factual example [4]. The knowledge
of counterfactual or semi-factual examples makes it possible to explain how to
change, to modify, the decisions of the model: “Your bank loan was not
accepted BUT IF you had more seniority in our company, the decision would have
been the opposite (or closer to acceptance)”. These two concepts are
illustrated in Figure 1. The understanding produced by a “counterfactual”
method of explanation is local, because it applies to a particular individual,
and “instance-based”, because it is produced in the form of a new example.

###  2.2 Informativeness and actionability of variables

In decision making, and particularly in the context of counterfactual
reasoning, the identification of informative and actionable variables is
essential. An informative variable is defined as a variable that has a
significant impact on the value of the output of the predictive model.
However, it is not enough to know which variables are informative. It is also
important to identify the actionable variables, which we define as a variable
on which it is possible to act. The most valuable type of variable is the
actionable informative variable, i.e. one that not only has a significant
impact on the output variable, but can also be acted upon to improve or
influence the result.

###  2.3 The concept of trajectories

The scientific literature offers numerous methods for generating
counterfactual examples, either model dependent (depending on the type of
classifier) or agnostic to the predictive model [17, 3]. In this article we
are interested in methods that induce the notion of trajectory. If we take
Figure 1, we are interested in methods which allow us to approach the decision
frontier step by step (until $X$ cross it) by successive modifications of the
initial example. In this case, the succession of $X^{\prime}$ resulting from
these successive operations is called a ‘trajectory’ (see Figure 2).

Figure 2: Illustration of two counterfactuals: one achieved in 1 step, the
second in 3 steps

Over the course of this article, we will be particularly interested in
predictive methods and/or models that allow for a notion of additivity in this
trajectory. That is, the univariate modifications of $X$ add up and make it
possible to approach the decision frontier by variable-by-variable
modification, step by step, as shown in Figure 2. We will show later that this
is the case for the Naive Bayes classifier.

##  3 Optimised search for counterfactuals in the case of the naive Bayes
classifier

###  3.1 Reminders on the naive Bayes classifier

The naive Bayes classifier (NB) is a widely used tool in supervised
classification problems. It has the advantage of being efficient for many real
data sets [7]. However, the naive assumption of conditional independence of
the variables can, in some cases, degrade the classifier’s performance. This
is why variable selection methods have been developed [11]. They mainly
consist of variable addition and deletion heuristics to select the best subset
of variables maximizing a classifier performance criterion, using a wrapper-
type approach [6]. It has been shown that averaging a large number of
selective naive Bayes classifiers, performed with different subsets of
variables, amounts to considering only one model with a weighting on the
variables. Bayes’ formula under the assumption of independence of the input
variables conditionally to the class variable becomes:

$P(C_{k}|X)=\frac{P(C_{k})\prod_{i}P(X_{i}|C_{k})^{W_{i}}}{\sum_{j=1}^{K}(P(C_{j})\prod_{i}P(X_{i}|C_{j})^{W_{i}})}$
(1)

where $W_{i}\in[0,1]$ is the weight of variable $i$. The predicted class is
the one that maximizes the conditional probability $P(C_{k}|X)$. The
probabilities $P(X_{i}|C_{i})$ can be estimated by interval using
discretization for numerical variables. Gaussian naive Bayes could be also
considered. For categorical variables, this estimation can be done directly if
the variable takes few different modalities, or after grouping (of values) in
the opposite case.

###  3.2 Criteria to be optimised in the search for counterfactuals

Lete $X$ as an example, and two classes $C_{1}$ and $C_{2}$. The search for a
counterfactual consists in optimising and increasing the probability of
belonging to the target class $C_{1}$ when $X$ is initially predicted by the
model to belong to $C_{2}$ (and vice versa). To do this, we can develop a
gluttonous algorithm, which is expensive in terms of computation time and does
not necessarily have the additivity properties described above in Section 2.3.
We propose below to pose the problem differently, rewriting the equation 1 and
looking at how to increase the probability of belonging to a particular class
of interest. To achieve this goal, and to maximise $P(C_{j}|X^{\prime})$ with
respect to the initial value of $P(C_{j}|X)$, we will exploit the following
proposition:

###### Proposition 1

If we take $X$ and $X^{\prime}$ as two elements of the input space
$\mathcal{X}$, we show that for a two-class classification problem, searching
for counterfactuals of $X$ amounts to examining the evolution of the value of
$\Delta$ when we change some of the values of $X$ to $X^{\prime}$, such that:

$\begin{split}\Delta(X,X^{\prime})&=\left(\sum_{i=1}^{d}W_{i}(\log{P(X_{i}|C_{1})}-\log{P(X_{i}|C_{2})})\right)\\\
&-\left(\sum_{i=1}^{d}W_{i}(\log{P(X^{\prime}_{i}|C_{1})}-\log{P(X^{\prime}_{i}|C_{2})})\right)\end{split}$
(2)

Proof: If we start again from the equation 1

$P(C_{j}|{X})=\frac{P(C_{j})\prod_{i=1}^{d}P(X_{i}|C_{j})^{W_{j}}}{\sum_{z}[P(C_{z})\prod_{i=1}^{d}P(X_{i}|C_{z})^{W_{j}}]}$

by posing:

$\small
L_{j}({X})=\log\left(P(C_{j})\prod_{i=1}^{d}P(X_{i}|C_{j})^{W_{i}}\right)=\log{P(C_{j})}+\sum_{i=1}^{d}W_{i}\log{P(X_{i}|C_{j})},$

then we have:

$P(C_{j}|{X})=\frac{e^{L_{j}({X})}}{\sum_{z}e^{L_{z}({X})}}=\frac{1}{\sum_{z}e^{L_{z}({X})-L_{j}({X})}}=\frac{1}{1+\sum_{z\neq
j}e^{L_{z}({X})-L_{j}({X})}}$

and so in the case of two classes:

$P(C_{j}|{X})=\frac{1}{1+e^{L_{j^{\prime}}({X})-L_{j}({X})}}$ (3)

We can see that to get closer to the class $C_{j}$, all we have to do is
reduce the quantity $L_{j^{\prime}}({X})-L_{j}({X})$, and thus reduce:

$\log{P(C_{j^{\prime}})}+\sum_{i=1}^{d}W_{i}\log{P(X_{i}|C_{j^{\prime}})}-\log{P(C_{j})}-\sum_{i=1}^{d}W_{i}\log{P(X_{i}|C_{j})}$

Since $P(C_{j})$ and $P(C_{j^{\prime}})$ are constant, this is equivalent to
decreasing:

$\sum_{i=1}^{d}W_{i}\log{P(X_{i}|C_{j^{\prime}})}-\sum_{i=1}^{d}W_{i}\log{P(X_{i}|C_{j})}$

and therefore to take an interest in the distance:

$\begin{split}\Delta({X},{X^{\prime}})&=\sum_{i=1}^{d}W_{i}(\log{P(X_{i}|C_{j^{\prime}})}-\log{P(X_{i}|C_{j})})\\\
&-\sum_{i=1}^{d}W_{i}(\log{P(X^{\prime}_{i}|C_{j^{\prime}})}-\log{P(X^{\prime}_{i}|C_{j})})\end{split}$

If $\Delta$ is positive then we are getting closer to the decision frontier
(or even crossing it) if $\Delta$ is negative then we are moving away from the
decision frontier and therefore away from the desired objective. The
counterfactual search algorithm becomes straightforward. Simply calculate, for
a given example $X$, the value of $\Delta$ for each explanatory variable and
for each value of this explanatory variable. Then, given these values, iterate
the successive changes in order to obtain a counterfactual example. These
variable-by-variable changes have the property of being additive.

Indeed, if we consider four examples $X^{0}$, $X^{{}^{\prime}1}$,
$X^{{}^{\prime}2}$ and $X^{{}^{\prime}3}\in\mathcal{X}$, which are
respectively (i) an initial example $X^{0}$, then the same example for which
we have modified only one explanatory variable $l$ for $X^{{}^{\prime}1}$, $m$
for $X^{{}^{\prime}2}$, and finally an example that cumulates the two
univariate modifications $l$ and $m$ for $X^{{}^{\prime}3}$, such that :

$\exists!\ l\ \text{such as}\ X^{{}^{\prime}1}_{l}\neq X^{0}_{l}$

$\exists!\ m\ \text{such as}\ X^{{}^{\prime}2}_{m}\neq X^{0}_{m}\ \text{and}\
m\neq l$

and

$X^{{}^{\prime}3}_{k}=\begin{cases}X^{{}^{\prime}1}_{l},&\text{if }k=l\\\
X^{{}^{\prime}2}_{m},&\text{if }k=m\\\ X^{0}_{k},&\text{otherwise}\end{cases}$

then it is obvious, from the additivity over all the variables in the equation
2, that we have :
$\Delta({X_{0}},{X^{\prime}_{3}})=\Delta({X_{0}},{X^{\prime}_{1}})+\Delta({X_{0}},{X^{\prime}_{2}})$.
Modifying one variable and then the other is equivalent to modifying them
simultaneously in the calculation of $\Delta$. It should also be noted that
this additivity is demonstrated from the equation 3, so we can be sure of
increasing the value normalised of the probability of the class of interest,
$P(C|X)$, which is a plus.

Note: The list of $\Delta$ values can potentially be very large if the number
of distinct values of the explanatory variables is large. Nevertheless, it is
common for some naive Bayes classifiers (except the Gaussian version) [20, 21]
to discretise the numerical variables and group the modalities of the
categorical variables, in a supervised manner, in order to obtain an estimate
of the conditional densities ($P(X_{i}|C)$) which are then used in the
calculation of $P(C|X)$. These supervised discretisation and grouping
operations often produce a limited number of intervals or groups of
modalities. This makes it possible to obtain a reasonable number of values to
test.

##  4 Creation and use of a knowledge base

###  4.1 Creation a knowledge base

In the foregoing we have shown how to increase the probability of belonging to
a class of interest and quantify this increase using the equation 2. We have
also shown that this quantity is additive as we change the values of the
explanatory variables one by one. We now propose to store these $\Delta$
values in a table in the form shown in Table 1.

It should be remembered that it is assumed that the numerical variables have
been discretised beforehand and that a grouping of modalities has been carried
out for the categorical variables. Each variable is therefore represented by a
limited number of values (corresponding to the values of $P(X|C)$). Table 1
gives, for illustrative purposes, the values stored in the case where the
predictive model uses two explanatory variables $X_{1}$ and $X_{2}$
respectively discretised (or grouped) into 3 and 2 intervals (groups) of
values ($I$). For each individual, $l$, each row of the table, we store the
values of the equation 2 where $\Delta(X^{l}_{i,*\rightarrow m},X^{l})$ is the
value of $\Delta$ when the value of variable $i$ is changed from its initial
value $`*^{\prime}$ to the value of interval (group) $m$ (simplified as
$\Delta(X^{l}_{i,*\rightarrow m})$ in the table).

| $X_{1}$ | $X_{2}$
---|---|---
| $I_{1}$ | $I_{2}$ | $I_{1}$ | $I_{2}$ | $I_{3}$
${X^{1}}$ | $\Delta({X^{1}_{1,*\rightarrow 1}})$ | $\Delta({X^{1}_{1,*\rightarrow 2}})$ | $\Delta({X^{1}_{2,*\rightarrow 1}})$ | $\Delta({X^{1}_{2,*\rightarrow 2}})$ | $\Delta({X^{1}_{2,*\rightarrow 3}}$
${X^{2}}$ | $\Delta({X^{2}_{2,*\rightarrow 1}})$ | $\Delta({X^{2}_{2,*\rightarrow 2}})$ | $\Delta({X^{2}_{2,*\rightarrow 1}})$ | $\Delta({X^{2}_{2,*\rightarrow 2}})$ | $\Delta({X^{2}_{2,*\rightarrow 3}})$
Table 1: Illustration of the knowledge base created in the form of $\Delta$,
here in the case of two variables and two examples.

We detail in the following sections how to exploit the knowledge stored in
this way. In this article, only the naive Bayes classifier is considered, but
any other classifier and/or counterfactual creation method that produces a
similar data table could be used.

###  4.2 Generation of counterfactuals with certain properties

####  4.2.1 Minimising the number of changes -

We can set ourselves the criterion of finding the counterfactual with the
property of involving the smallest number of modified variables. To do this,
we will exploit the knowledge base. For a given individual, $X$, all we have
to do is read the value of the largest $\Delta$ and then, as we have the
additivity property, read the second value of the largest $\Delta$ for a
second variable and so on. At each stage we check whether
($\hat{f}(\mathbf{X^{\prime}})>0.5$). If this is the case, the counterfactual
has been found footnoteWe could, on the other hand, maximise the number of
changes, but this is often of little interest in practice.

####  4.2.2 Taking account of business constraints or criteria -

In other cases, the aim may be to find the closest counterfactual, but under
‘business constraints’ defined by the user. For example, the search for
counterfactuals could be restricted to making changes only in adjacent
intervals (e.g. intervals of close values). Given the table 1, we would be
allowed to move people from interval 1 to 2 for the second variable, but not
from 1 to 3. The user can also constrain the search by requiring that one of
the variables must always be changed to a certain value, and so on. This type
of constraint can easily be considered and incorporated into a counterfactual
search algorithm using the proposed knowledge table.

The literature on counterfactuals sets out some interesting properties on the
subject, such as (i) the notion of minimality: having a counterfactual that
differs as little as possible from the original example; (ii) realism: the
generated counterfactual must not contain changes that do not make sense from
the point of view of the data (e.g. a decrease in the “age” of an individual),
also known as plausibility [15]; (iii) generating counterfactuals that are
similar to real examples or in dense regions of the class of interest, in
order to have robust counterfactuals [5]; in this case, in the case of the
naive Bayes classifier, we can use a Bayesian distance as proposed in [13]; …

All these properties are easily achieved with the Table 1, since the user can
choose the list and the order of the variables he wish to intervene on, as
well as a distance of his choice between $X$ and $X^{\prime}$.

###  4.3 Additional usable knowledge

####  4.3.1 Preventive and reactive actions -

So far, we have mainly talked about creating counterfactuals to explain the
model’s decisions (as mentioned in the introduction to this article), but
potentially also to be able to take reactive actions. For example, if a bank
customer is predicted to ”leave” (churner), the counterfactual example
indicates one or more actions to be taken in order to try to keep him as a
customer: these are known as “reactive” actions.

Conversely, the study of counterfactual trajectories a posteriori is of great
interest, as it also allows us to identify when a trajectory is approaching
the frontier. In such situations, reactive measures can be taken to reverse
the trend and avoid undesirable outcomes. This approach is particularly
relevant when it comes to predicting churn, for example, as it enables us to
identify customers who are “starting” to churn. By being proactive, it is
possible to put in place targeted strategies to retain these customers and
bring them back to a quality service.

Finally, our knowledge base can also be used to carry out “preventive”
actions. Going back to Figure 1, we can try to create a semi-factual which
moves away from the decision frontier: ”The customer is not predicted as
leaving but is nevertheless close to the decision frontier”. In this case, all
we need to do is look at the negative values of $\Delta$ and take steps away
from it according to the user’s wishes. For example, all the people who are
one step away from crossing the decision frontier, who are easily identifiable
in this case, could be concerned.

####  4.3.2 Profile creation -

The last way of using the knowledge base that we will describe here111The
reader can imagine others: descriptive statistics of the table, number of
individuals at 1, 2, 3 … steps from the decision frontier, visualisation of
the trajectories, … ., is to carry out an exploratory analysis using a
clustering technique. Using the knowledge base, it is possible to group
individuals according to the impact of each possible change, i.e. the impact
resulting from each of these changes. Analysis of the clusters created can be
a source of learning. This is illustrated in the next section.

##  5 Illustration of a clustering on an unsubscribe case

1) Dataset and classifier used : This section uses the “Telco Customer Churn”
dataset (widely used in analysing the results of the XAI method) provided by a
fictitious telecommunications company that provided home telephone and
internet services to 7043 customers in California. The aim is to classify
people who may or may not leave the company. Each customer is described by 20
descriptive variables (3 numerical and 17 categorical) plus the class variable
‘churn’ (yes/no) which has two modalities (75% non-churn). This dataset can be
downloaded from Kaagle [9]. We use 80% of the data for learning and 20% for
testing. The naive Bayes classifier is produced using the Khiops library,
which is available on Github [10].

During the learning process, only 10 of the variables were retained in the
model. Below are all the intervals of values or groups of modalities obtained
during the pre-processing process (the value in brackets gives the weight of
the variable in the model, equation 1, values from 0 to 1):

  * •

1 * - Tenure ($W_{1}$=0.67) : [0-0.5], ]0.5-1.5], ]1.5-5.5], ]5.5-17.5],
]17.5-42.5], ]42.5-58.5], ]58.5-71.5], ]71.5-72]

  * •

2 - InternetService ($W_{2}$=0.78) : [Fiberoptic], [DSL], [No]

  * •

3 - Contract ($W_{3}$=0.37) : [Month-to-month], [Twoyear], [Oneyear]

  * •

4 - PaymentMethod ($W_{4}$=0.29) : [Mailedcheck], [Creditcard(automatic),
Electroniccheck, Banktransfer(automatic)]

  * •

5 - OnlineSecurity ($W_{5}$=0.15) : [No], [Yes], [No internet service]

  * •

6 -TotalCharges ($W_{6}$=0.29) : [18.8;69.225], ]69.225;91.2], ]91.2;347.9],
]347.9;1182.8],

]1182.8 ;2065.7], ]2065.7;3086.8], ]3086.8;7859], ]7859;8684.8]

  * •

7 * - PaperlessBilling ($W_{7}$=0.40) : [Yes], [No]

  * •

8 - TechSupport ($W_{8}$=0.04) : [No], [Yes], [No internet service]

  * •

9 * - SeniorCitizen ($W_{9}=0.28$): [0], [1]

  * •

10 * - Dependents ($W_{9}=0.10$): [Yes], [No]

For all 10 variables, there are a total of 36 intervals/groupings and
therefore 26 values of $\Delta$ to calculate in our knowledge base. Indeed,
for each individual and each variable, there is a $\Delta$ value that has a
null value, the value which corresponds to it factually and which therefore
does not need to be calculated.

2) Classifier analysis stage: Before carrying out the clustering stage, it is
important to take an interest in the variables retained during the
classification model training stage. For example, although it may be
interesting to analyse the ‘Tenure’ variable, it is clearly not an actionable
variable. Indeed, it is not possible to change a customer’s seniority in order
to make him potentially less unfaithful. The same applies to the
‘SeniorCitizen’ and ‘Dependents’ variables. We have also removed the
‘PaperlessBilling’ variable, which has very little impact on the clustering
results described below. As a result, these 4 variables are not retained
during the clustering stage below; only the informative, influential and
actionable variables are retained 222All the variables could have been
retained but the clustering would have been biased by uninteresting variables
from the point of view of creating counterfactual examples (see Section 2.2).

3) Clustering performed : The clustering performed is usual: (i) we use the
table of $\Delta$ values calculated on the test set, (ii) we learn a k-means
with the L2 [8] distance for different $k$ values ($k\in\\{2,12\\}$), (iii)
and finally we retain the k-means whose value of $k$ corresponds to the ‘elbow
point’, here $k=4$, [18] of the curve representing the global reconstruction
distance versus the value of $k$.

Figure 3: Average profile of individuals in clusters represented as
histograms. The names of the values on the abscissa refer to the number of the
variables and number of the intervals (groups) described above. For example,
‘3I2’ refers to the third variable (’Contract’) and its second interval /
group of values (’Twoyear’). The ordinate values are the mean values of the
cluster ($\Delta$).

4) The resulting clusters are shown in Figure 3. An analysis of these 4
clusters, combined with the predictions of the classifier, shows that:

  * •

Cluster 1 (10% of the global population and containing 2% of customers
predicted to be churner): individuals who can be made less churner mainly by
means of variable 3 (‘Contract’) - i.e. by trying to get them to take out an
annual contract (‘Twoyear’ or ‘OneYear’); NB - this marketing action is fairly
difficult to carry out.

  * •

Cluster 2 (24% of the global population and containing no customer predicted
to be churner): people who are very insensitive to the fact that they are
becoming less churner (mostly negative $\Delta$ means). They may not be
targeted by a ‘reactive’ marketing campaign (which is in line with the
classifier’s predictions), but rather by a preventive campaign using the
‘contract’ variable or the ‘payment method’ variable (payment by card or
direct debit).

  * •

Cluster 3 (45% of the global population and containing 47% of customers
predicted to be chruner (almost all of the individuals predicted to be
churner)): some similarities with the individuals in cluster 1 for the
‘Contract’ variable. On the other hand, we can see that the 5th
(’OnlineSecurity’) and 8th (’TechSupport’) variables have a ‘leverage effect’
in reducing churn. Offering a security or support option is very attractive to
these individuals.

  * •

Cluster 4 (21% of the population and containing no customers predicted to be
unfaithful): individuals who are partially opposite to those in the first
cluster, for example for the ‘Contract’ variable, who should not be offered a
‘two-year contract’ in this case.

The analysis of the clusters obtained here is not exhaustive due to space
limitations. It is an exploratory analysis where the data scientist and the
business expert will spend the time needed to refine their joint analyses.
However, the analysis carried out here allows us to identify interesting
‘reactive’ actions to be taken with individuals in cluster 3 or preventive
actions with individuals in cluster 2.

##  6 Conclusion

In the context of methods for explaining the results of a machine learning
model, this article has proposed to consider the process of generating
counterfactual examples as a source of knowledge that can be stored and then
exploited in different ways. This process has been illustrated in the case of
additive models and in particular in the case of the naive Bayes classifier,
whose interesting properties for this purpose have been shown. We have also
suggested the quantities that can be stored and the different ways of
exploiting them. Some of the results have been illustrated on a churn problem,
but the approach is equally exploitable in other application domains.

## References

  * [1] Allen, G.I., Gan, L., Zheng, L.: Interpretable machine learning for discovery: Statistical challenges & opportunities. Arxiv preprint:2308.01475 (2023)
  * [2] Bodria, F., Giannotti, F., Guidotti, R., Naretto, F., Pedreschi, D., Rinzivillo, S.: Benchmarking and survey of explanation methods for black box models. Data Mining and Knowledge Discovery 37(5), 1719–1778 (2023)
  * [3] Brughmans, D., Leyman, P., AU Martens, D.: Nice: an algorithm for nearest instance counterfactual explanations. Data Mining and Knowledge Discovery (2023)
  * [4] Fernández, R.R., De Diego, I.M., Moguerza, J.M., Herrera, F.: Explanation sets: A general framework for machine learning explainability. Information Sciences 617, 464–481 (2022)
  * [5] Guyomard, V., Fessant, F., Guyet, T., Bouadi, T., Termier, A.: Generating robust counterfactual explanations. In: European Conference on Macine Learning (2023)
  * [6] Guyon, I., Elisseeff, A.: An introduction to variable and feature selection. J. Mach. Learn. Res. 3, 1157–1182 (2003)
  * [7] Hand, D.J., Yu, K.: Idiot’s bayes-not so stupid after all? International Statistical Review 69(3), 385–398 (2001)
  * [8] Hartigan, J.A., Wong, M.A.: A k-means clustering algorithm. JSTOR: Applied Statistics 28(1), 100–108 (1979)
  * [9] Kaagle: Telco customer churn dataset (2023), [https://www.kaggle.com/datasets/blastchar/telco-customer-churn], last visited 08/22/2023
  * [10] Khiops: Github khiops (2023), [https://github.com/KhiopsML/khiops], last visited 08/22/2023
  * [11] Langley, P., Sage, S.: Induction of selective bayesian classifiers. In: Proceedings of the Tenth International Conference on Uncertainty in Artificial Intelligence. pp. 399–406. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA (1994)
  * [12] Lemaire, V., Hue, C., Bernier, O.: Data Mining in Public and Private Sectors: Organizational and Government Applications, chap. Correlation Analysis in Classifiers, pp. 204–218. IGI Global (2010)
  * [13] Lemaire, V., Ismaili, O.A., Cornuéjols, A., Gay, D.: Predictive k-means with local models. In: Lu, W., Zhu, K.Q. (eds.) Trends and Applications in Knowledge Discovery and Data Mining. pp. 91–103. Springer International Publishing, Cham (2020)
  * [14] Miller, T.: Explanation in artificial intelligence: Insights from the social sciences. Artificial Intelligence 267, 1–38 (2019)
  * [15] Nemirovsky, D., Thiebaut, N., Xu, Y., Gupta, A.: Countergan: Generating counterfactuals for real-time recourse and interpretability using residual gans. In: Conference on Uncertainty in Artificial Intelligence. pp. 1488–1497. Machine Learning Research, PMLR (2022)
  * [16] Saeed, W., Omlin, C.: Explainable ai (xai): A systematic meta-survey of current challenges and future opportunities. Knowledge-Based Systems 263, 110273 (2023)
  * [17] Stepin, I., Alonso, J.M., Catala, A., Pereira-Fariña, M.: A survey of contrastive and counterfactual explanation generation methods for explainable artificial intelligence. IEEE Access 9, 11781–11803 (2021). https://doi.org/10.1109/ACCESS.2021.3051315
  * [18] Thorndike, R.L.: Who belongs in the family? Psychometrika 18, 267–276 (1953), https://api.semanticscholar.org/CorpusID:120467216, the method can be traced to speculation by Robert L. Thorndike
  * [19] Wachter, S., Mittelstadt, B., Russell, C.: Counterfactual explanations without opening the black box: Automated decisions and the gdpr. Harvard Journal of Law and Technology 31(2), 841–887 (2018)
  * [20] Yang, Y., Webb, G.: A comparative study of discretization methods for naive-bayes classifiers. In: Proceedings of PKAW, vol. 2002 (04 2003)
  * [21] Yang, Y., Webb, G.: Discretization for naive-bayes learning: Managing discretization bias and variance. Machine Learning 74, 39–74 (09 2009). https://doi.org/10.1007/s10994-008-5083-5