# The Effect of Nonlinear Landau Damping on Ultrarelativistic Beam Plasma
Instabilities

Philip Chang11affiliation: Department of Physics, University of Wisconsin-
Milwaukee, 1900 E. Kenwood Boulevard, Milwaukee, WI 53211, USA Avery E.
Broderick22affiliation: Perimeter Institute for Theoretical Physics, 31
Caroline Street North, Waterloo, ON, N2L 2Y5, Canada 33affiliation: Department
of Physics and Astronomy, University of Waterloo, 200 University Avenue West,
Waterloo, ON, N2L 3G1, Canada Christoph Pfrommer44affiliation: Heidelberg
Institute for Theoretical Studies, Schloss-Wolfsbrunnenweg 35, D-69118
Heidelberg, Germany Ewald Puchwein55affiliation: Institute of Astronomy and
Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge CB3 0HA, UK Astrid Lamberts11affiliation: Department of Physics,
University of Wisconsin-Milwaukee, 1900 E. Kenwood Boulevard, Milwaukee, WI
53211, USA Mohamad Shalaby22affiliation: Perimeter Institute for Theoretical
Physics, 31 Caroline Street North, Waterloo, ON, N2L 2Y5, Canada
33affiliation: Department of Physics and Astronomy, University of Waterloo,
200 University Avenue West, Waterloo, ON, N2L 3G1, Canada 66affiliation:
Department of Physics, Faculty of Science, Cairo University, Giza 12613,
Egypt.

###### Abstract

Very-high energy gamma-rays from extragalactic sources pair-produce off of the
extragalactic background light, yielding an electron-positron pair beam. This
pair beam is unstable to various plasma instabilities, especially the
“oblique” instability, which can be the dominant cooling mechanism for the
beam. However, recently, it has been claimed that nonlinear Landau damping
renders it physically irrelevant by reducing the effective damping rate to a
low level. Here, we show with numerical calculations that the effective
damping rate is $8\times 10^{-4}$ of the growth rate of the linear
instability, which is sufficient for the “oblique” instability to be the
dominant cooling mechanism of these pair beams. In particular, we show that
previous estimates of this rate ignored the exponential cutoff in the
scattering amplitude at large wavenumber and assumed that the damping of
scattered waves entirely depends on collisions, ignoring collisionless
processes. We find that the total wave energy eventually grows to approximate
equipartition with the beam by increasingly depositing energy into long
wavelength modes. As we have not included the effect of nonlinear wave-wave
interactions on these long wavelength modes, this scenario represents the
“worst-case” scenario for the oblique instability. As it continues to drain
energy from the beam at a faster rate than other processes, we conclude that
the “oblique” instability is sufficiently strong to make it the physically
dominant cooling mechanism for high-energy pair beams in the intergalactic
medium.

###### keywords:

BL Lacertae objects: general – gamma rays: general – plasmas – instabilities –
magnetic fields

##  1 Introduction

The very-high energy gamma-ray (VHEGR, $E\geq 100\,{\rm G}{\rm eV}$)
extragalactic sky is dominated by a subset of blazars, which we refer to as
TeV blazars. These extragalactic VHEGR emitters produce TeV photons that are
greatly attenuated via annihilation upon soft photons in the extragalactic
background light (EBL) and produce pairs (see, e.g., Gould & Schréder, 1967;
Salamon & Stecker, 1998; Neronov & Semikoz, 2009).

It has been assumed that these ultrarelativistic pairs produced by VHEGR
annihilation lose energy exclusively through inverse-Compton scattering off of
the cosmic microwave background (CMB), transferring the energy of the original
VHEGR to gamma-rays with energies $\lesssim 100\,{\rm G}{\rm eV}$. As the
gamma-rays are in the LAT bands of Fermi, this is an important target for
Fermi observations.

The absence of observed secondary IC emission leads a number of authors to
argue that this lack of observed emission places lower bounds upon the
intergalactic magnetic field (IGMF; see, e.g., Neronov & Vovk, 2010; Tavecchio
et al., 2010, 2011; Dermer et al., 2011; Taylor et al., 2011; Takahashi et
al., 2012; Dolag et al., 2011) with typical numbers ranging from
$10^{-19}\,{\rm G}$ to $10^{-15}\,{\rm G}$. In addition, the persistent belief
in this IC emission has led other workers to argue that based on the IC
contribution to Fermi extragalactic gamma-ray background (EGRB), the comoving
number density of gamma-ray bright blazars must grow slowly with increasing
redshift, if at all (Kneiske & Mannheim, 2008; Venters, 2010; Abazajian et
al., 2011; Inoue & Ioka, 2012), implying that these TeV blazars cosmologically
evolve qualitatively differently compared to other active galactic nuclei
(AGNs).

These two conclusions depend on IC cooling dominating the evolution of the
ultra-relativistic pairs. However, it was recently found that plasma
instabilities driven by the ultrarelativistic pair beams likely is the
dominant cooling mechanism (Broderick et al., 2012; Schlickeiser et al., 2012,
2013), depositing this energy as heat in the intergalactic medium (Chang et
al., 2012; Pfrommer et al., 2012). Therefore, the lack of an observed IC halo
emission from TeV blazars does not imply the existence of the IGMF as previous
groups have argued (Broderick et al., 2012; Schlickeiser et al., 2012, 2013).

This excess heating of the IGM may resolve a variety of cosmological puzzles,
including naturally explaining anomalies in the statistics of the high-
redshift Ly$\alpha$ forest (Puchwein et al., 2012) and potentially explaining
a number of the X-ray properties of groups and clusters and anomalies in
galaxy formation on the scale of dwarfs (Pfrommer et al., 2012; Lu et al.,
2013). We have recently shown that if the IC halos are ignored, it is possible
to quantitatively reproduce the redshift and flux distributions of nearby hard
gamma-ray blazars and the extragalactic gamma-ray background spectrum above 3
GeV simultaneously with a unified model of AGN evolution (Broderick et al.,
2013a, b). All of these empirical successes provide circumstantial evidence
for the presence of virulent plasma beam instabilities.

Recently, Miniati & Elyiv (2013, hereafter ME13) argued that these
instabilities are physically irrelevant for the cooling of these pair beams.
First, the “oblique” instability would saturate at a very low level due to
nonlinear Landau damping (NLD). ME13 argue that this process occurs when
$3\times 10^{-6}$ of the electron-positron beam energy is contained within the
waves and significantly limiting the instability cooling rate. Second,
inhomogeneities in the IGM prevent this linear instability from even growing
by shifting the unstable waves out of the region of resonance. The second
point relates to the properties of the linear instability in the kinetic
regime, which has been shown by Schlickeiser et al. (2013) to be physically
relevant, contrary to ME13’s claims. We will address this second objection of
ME13 in future work.

In this paper, we will consider the effects of NLD. We begin by discussing the
physics in §2. Using a numerical calculation, we derive the saturation level
of the “oblique” instability in the non-linear regime in §3. We discuss why
our results differ from those presented by ME13 in §4 and present the
implications in §5. Finally, we close with a summary of results and pathway to
future work in §6.

##  2 The Physics of Nonlinear Landau Damping

We consider an unstable wave that is driven by a beam of electrons and
positrons in a background plasma of electrons and protons (or other ions). As
the unstable wave grows in amplitude, it becomes subject to nonlinear wave-
particle and wave-wave interactions. In nonlinear particle-wave interactions,
the most important interaction is induced scattering by thermal ions (Kaplan &
Tsytovich, 1968; Smith & Fung, 1971; Breǐzman et al., 1972; Melrose, 1986),
which is also referred to as NLD (Melrose, 1986).

Wave-particle interactions induce the transformation of one plasma wave,
characterized by a frequency and wavevector $(\omega,{\bf k})$, into another
$(\omega^{\prime},{\bf k}^{\prime})$ via nonlinear scattering on the particles
that constitute the plasma. The kinetic equation for these waves in the
presence of wave-particle interactions is (Kaplan & Tsytovich, 1968; Breǐzman
et al., 1972)

$\frac{dW_{{\bf k}}}{dt}=2\Gamma_{{\bf k}}W_{{\bf k}}-\frac{W_{{\bf
k}}\omega_{p}}{8(2\pi)^{5/2}n_{e}m_{e}v_{e}^{2}}\int\frac{({\bf k}\cdot{\bf
k}^{\prime})^{2}}{k^{2}k^{\prime 2}}\phi({\bf k},{\bf k}^{\prime})W_{{\bf
k}^{\prime}}d{\bf k}^{\prime}$ (1)

where $W_{{\bf k}}$ is the spectral energy density of the waves, normalized
such that the total energy density is given by

$W=\frac{1}{(2\pi)^{3}}\int W_{{\bf k}}d{\bf k}\,,$ (2)

$\Gamma_{{\bf k}}=\Gamma_{\rm gr}({\bf k})+\Gamma_{\rm LD}({\bf k})$ is the
sum of the unstable wave growth rate, $\Gamma_{\rm gr}({\bf k})$, and the
linear Landau damping rate, $\Gamma_{\rm LD}({\bf k})$,
$\omega_{p}\equiv\sqrt{4\pi n_{e}e^{2}/m_{e}}$ is the electron plasma
frequency of the background plasma, given in terms of the proper electron
density ($n_{e}$) and rest mass ($m_{e}$), and $v_{e}$ and $v_{i}$ are the
electron and ion thermal velocities, respectively. We also note that wave-
particle interactions also convert these electrostatic waves into
electromagnetic waves as noted in Kaplan & Tsytovich (1968). However, we
ignore electromagnetic modes in this work to focus on the the electrostatic
waves. Including these modes as well other nonlinear processes is the subject
of future work.

The growth rate is given by the oblique growth rate, which is

$\Gamma_{\rm gr}({\bf k})\equiv\frac{1}{\tau_{\rm gr}}\approx
0.4\frac{n_{b}}{n_{e}}\gamma_{b}\omega_{p}\Theta(1-v_{\rm ph}(k)/c),$ (3)

where $n_{b}$ is the beam density and $\gamma_{b}$ is the Lorentz factor of
the ultrarelativistic beam, $\Theta$ is the Heaviside function, and $v_{\rm
ph}(k)=\omega/k$ is the phase speed of the Langmuir wave. Equation (3) was
first found by Bret et al. (2010) by fitting the maximum growth rate in the
kinetic regime. We have confirmed this growth rate in the electrostatic
approximation for the kinetic regime (Broderick et al., in prep., see also
Schlickeiser et al. 2013). More importantly, this result holds true for a
large range in ${\bf k}$. The reason for this is that electrostatic waves
oscillate at $\omega_{p}$ almost independently of ${\bf k}$ and their phase
speed along the z-axis (arbitrarily defined) is given by $v_{\rm
ph}=\omega/k\cos\theta$, where $\theta$ is the angle between the direction in
question and the wave vector. Hence for any $k\geq\omega_{p}/c$, there exist
some $\cos\theta$ such that the $v_{\rm ph}\approx c$ and hence these waves
are in resonance with a relativistic beam. In other words, the oblique
instability grows for any $k\geq\omega_{p}/c$, but the angle between the most
unstable wavevector and beam varies.

Damping in the linear regime is given by linear Landau damping, whose rate is

$\Gamma_{\rm LD}({\bf
k})=-\left(\frac{\pi}{8}\right)^{1/2}\omega_{p}\left(\frac{v_{\rm
ph}}{v_{e}}\right)^{3}\exp\left(-\frac{v_{\rm ph}^{2}}{2v_{e}^{2}}\right).$
(4)

The overlap integral $\phi$ is given by (Kaplan & Tsytovich, 1968)

$\phi({\bf k},{\bf k}^{\prime})=\frac{3v_{e}^{2}\left(k^{2}-k^{\prime
2}\right)}{4\omega_{p}|{\bf k}-{\bf
k}^{\prime}|v_{i}}\exp\left[-2\left(\frac{3v_{e}^{2}\left(k^{2}-k^{\prime
2}\right)}{4\omega_{p}|{\bf k}-{\bf k}^{\prime}|v_{i}}\right)^{2}\right].$ (5)

An important feature of equation (5) is the dependence on
${k^{\prime}}^{2}-k^{2}$, which sets the sign of $\phi(k,k^{\prime})$.
Scattering of a wave with wavevector ${\bf k}$ into another wave with
wavevector ${\bf k}^{\prime}$ can only proceed if $\phi(k,k^{\prime})>0$,
i.e., the wave energy in the ${\bf k}$ wave is damped, while the ${\bf
k}^{\prime}$ wave energy grows. This requires that $k^{\prime}<k$, i.e., the
scattered wave has a longer wavelength than the incident wave. The demand that
induced scattering drives waves to longer wavelength arises from the transfer
of some momentum from the incident wave into the polarization clouds
surrounding the ions (Smith & Fung, 1971).

##  3 Numerical Studies

We now solve equation (1) numerically assuming that linear growth and
nonlinear damping via NLD are the two mechanisms that control the initial
evolution of Langmuir waves. However, because ${\bf k}$ is three-dimensional,
we adopt the simplifying assumption that the Langmuir waves are isotropic,
i.e., $W_{{\bf k}}=W_{k}$. This simplifying assumption is reasonable as long
as the induced scattering processes are sufficiently rapid that it isotropizes
the waves111Further support for this approximation emerges if the number of
TeV blazars that contribute to the heating of a given patch of the
intergalactic medium exceeds 100. (Kaplan & Tsytovich, 1968). Equation (1)
reduces then to222Equation (5) can also be simplified if we set ${({\bf
k}_{1}\cdot{\bf k}_{2})^{2}}/{k_{1}^{2}k_{2}^{2}}$ to the angle-averaged value
of $1/3$ as done by Kaplan & Tsytovich (1968). We computed this integral both
ways and found little difference in the saturation amplitude.

$\displaystyle\frac{dW_{k}}{dt}$ $\displaystyle=$ $\displaystyle
2\Gamma_{k}W_{k}-\frac{W_{k}\omega_{p}}{8(2\pi)^{3/2}n_{e}m_{e}v_{e}^{2}}$ (6)
$\displaystyle\times\int k^{\prime 2}\cos^{2}\theta\,\phi({\bf k},{\bf
k}^{\prime})W_{k^{\prime}}dk^{\prime}d\cos\theta$

where $\theta$ is the angle between ${\bf k}$ and ${\bf k}^{\prime}$. Without
loss of generality, we have fixed ${\bf k}$ along the z-axis. Here $\phi({\bf
k},{\bf k}^{\prime})$ can be simplified to

$\displaystyle\phi({\bf k},{\bf k}^{\prime})$ $\displaystyle=$
$\displaystyle\frac{3v_{e}^{2}\left(k^{2}-k^{\prime
2}\right)}{4\omega_{p}v_{i}\sqrt{k^{2}+k^{\prime 2}-2kk^{\prime}\cos\theta}}$
(7)
$\displaystyle\times\exp\left[-2\left(\frac{3v_{e}^{2}\left(k^{2}-k^{\prime
2}\right)}{4\omega_{p}v_{i}\sqrt{k^{2}+k^{\prime
2}-2kk^{\prime}\cos\theta}}\right)^{2}\right].$

We calculate equation (6) numerically for $N_{\rm modes}=300$ logarithmically
spaced modes from $k=10^{-6}\omega_{p}/c$ to $10^{3}\omega_{p}/c$ and have
verified this calculation using $N_{\rm modes}=1000$. The lower limit of
$k_{\rm min}$ was chosen to fulfill the requirement $k_{\rm
min}\ll\omega_{p}/c$. The upper limit of $k_{\rm max}=10^{3}\omega_{p}/c$ was
set because it is significantly larger than the estimated $k$ where linear
Landau damping would suppress the instability. We have found that our
calculations are not affected by extending the upper and lower limits on k. We
also set the initial $k^{3}W_{k}$ to a small value of the initial beam energy,
i.e., $10^{-15}$ and confirmed that our results are independent of this
initial value.

In Figure 1, we plot the wave energy $W$ in units of the initial beam energy
density $n_{b}\gamma_{b}m_{e}c^{2}$ as a function of growth e-folding times
$\Gamma_{\rm gr}t$ for a $1$ TeV beam and a $10$ TeV beam, where the beam
density is what is expected at $z=0$ for a TeV blazar with equivalent
isotropic luminosity $EL_{E}=10^{45}~{}\mathrm{erg~{}s}^{-1}$. The wave energy
grows exponentially up to a time $\Gamma_{\rm gr}t\approx 15$, where
exponential growth ends and transitions to a slow linear growth in $W$. The
wave energy, $W$, is equal to the energy of the initial beam when $\Gamma_{\rm
gr}t=2000$. Therefore, in the absence of a significant back-reaction, the beam
experiences one e-folding reduction in energy at $\Gamma_{\rm gr}t\approx
1300$, giving a damping rate of

$\Gamma_{\rm NLD}\approx 8\times 10^{-4}\Gamma_{\rm gr}.$ (8)

fig1.pdf

Figure 1: Wave energy, $W$, in units of the beam energy,
$n_{b}\gamma_{b}m_{e}c^{2}$, as a function of growth e-folding times
$\Gamma_{\rm gr}t$. The dot-dashed curve is a line defined by $W\propto
10^{-3}\Gamma_{\rm gr}t$.

In Figure 2 we plot the wave energy as a function of wavevector, $k$, for
three different times, $\Gamma_{\rm gr}t=$ 12 (solid line), 15 (dotted line),
and 30 (dashed line). These times are also marked in Figure 1 by vertical
lines of the same type as in Figure 2. For $\Gamma_{\rm gr}t=$ 12 (solid
line), NLD is not important and the instability grows for all
$k\geq\omega_{p}/c$. As the unstable waves grow, the effect of NLD begins to
become important and long wavelength modes $k<\omega_{p}/c$ begin to grow at
the expense of short wavelength modes. This is clearly seen for $\Gamma_{\rm
gr}t=$ 15 (dotted line) and 30 (dashed line). However, the largest wavevector
modes are not suppressed by the effect of NLD and remain at a level of
$\approx 10^{-3}$ of the beam energy density. These large wavevector modes
survive for the duration of the calculation and continually pump energy into
long wavelength modes.

fig2.pdf

Figure 2: Wave energy, $k^{3}W_{k}$, in units of the beam energy,
$n_{b}\gamma_{b}m_{e}c^{2}$, as a function of wavevector $k$ for three
different times $\Gamma_{\rm gr}t=$12 (solid line), 15 (dotted line), and 30
(dashed line).

To understand the origin of the survival of these large wavevector modes and
the development of the empty region between $k\approx 10^{-2}\omega_{p}/c$ and
$\sim 20\omega_{p}/c$, it is helpful to return to the coupling term (5).
Because the argument in the exponent is dominated by $k^{\prime}\approx k$,
modes that are closely spaced to each other are strongly scattered. When the
difference is large, i.e., $k^{2}-k^{\prime 2}$ is large, the scattering is
exponentially suppressed. Therefore, NLD is strongest on modes near each
other. As a result, as the long wavelength modes grow, they quickly sap energy
from the nearest modes, leaving the short wavelength modes untouched due to
the exponential suppression. This suppression becomes important where the
exponent is order unity or

$\frac{v_{e}^{2}\left(k^{2}-k^{\prime 2}\right)}{\omega_{p}|{\bf k}-{\bf
k}^{\prime}|v_{i}}\approx\frac{v_{e}^{2}k}{\omega_{p}v_{i}}\gtrsim 1$ (9)

This gives the condition

$k\gtrsim\sqrt{\frac{m_{e}}{m_{i}}}\frac{\omega_{p}}{c}\frac{c}{v_{e}}\approx
20\left(\frac{T}{10^{4}\,{\rm K}}\right)^{-1/2}\frac{\omega_{p}}{c},$ (10)

which is of order where the suppression occurs (see Figure 2).

##  4 Relationship to Previous Work

###  4.1 Comparison with Miniati & Elyiv (2013)

Our result differs from ME13, who found that the effect of NLD on the
“oblique” instability is to drive the saturation of the excited Langmuir waves
to a physically irrelevant amplitude. ME13 estimates the excited wave energy
to be $W/(n_{b}\gamma m_{e}c^{2})\approx 3\times 10^{-6}$, whereas our
numerical calculation finds a value that is over two orders of magnitude
larger. There are two crucial differences. First, ME13 based their estimate on
the order-of-magnitude calculation which ignores the effect of the suppression
of NLD at large wavenumber. Second, ME13 assumes that the damping of these
scattered waves is completely due to collisions, which are extraordinarily
slow. In this case, the wave energy of the excited waves is reduced to
${W_{r}}/{\gamma_{b}n_{b}m_{e}c^{2}}={\Gamma_{\rm NLD}}/{\Gamma_{\rm c}},$
where $W_{r}$ is the wave energy of resonant waves, $\Gamma_{\rm c}$ is the
electron collision rate and is typically given by electron-ion collisions, and
$\Gamma_{\rm NLD}$ is the maximum growth rate of the resonant waves that are
unstable to growth due to NLD of the linearly unstable mode.

The choice of collisional damping likely underestimates the true damping rate.
In particular, collisionless processes like wave-wave scattering and wave-
particle interaction are likely to produce a damping rate for these scattered
waves that far exceeds collisional damping. Indeed, collisionless simulations
(Davidson, 1972) have shown that nonlinear wave-wave interactions can lead to
particle heating and wave damping in the absence of collisions.

Here we presume that the damping rate of the scattered waves are sufficiently
rapid to give

$\frac{W_{r}}{\gamma_{b}n_{b}m_{e}c^{2}}=\frac{\Gamma_{\rm NLD}}{\Gamma_{\rm
gr}}\approx 8\times 10^{-4}.$ (11)

We are motivated to use this estimate by comparison with the work of Ziebell
et al. (2008b, a) who found while NLD effects dominate the scattering of
Langmuir waves, the rate of three-wave interactions is competitive with NLD
and results in the quasi-isotropic heating of the electrons. The results of
Ziebell et al. (2008b, a) suggest that collisional processes are irrelevant to
the damping of scattered waves. Moreover, nonlinear wave-wave coupling will
likely lead to a more equitable distribution of mode energies among the
different wavevectors whereas the effect of NLD moves energy toward smaller
wavevectors. Thus, inclusion of nonlinear wave-wave coupling will lead to a
stronger damping rate for the oblique instability by countering the effect of
NLD. As a result, the calculation that we present here likely represents the
worst case scenario for the “oblique” instability.

###  4.2 Comparison with Sironi & Giannios (2014)

Recent numerical simulations have been brought to bear upon these dilute beam-
plasma processes to study the physics of nonlinear saturation by Sironi &
Giannios (2014, hereafter SG14). However, due both to their spectral
resolution and constraints on the parameters simulated these appear to be
unable to capture the physics of NLD.

The range of $k$ over which NLD redistributes energy is limited by the
particular form of the overlap integral, eq. 5. To estimate this range, we
note that the overlap integral has a maximum width for ${\bf k}^{\prime}$
antiparallel to ${\bf k}$, i.e., the backward scattering case; in this case,
we take ${\bf k}^{\prime}=-(k-\Delta k)\hat{{\bf k}}$ and apply this to
equation (5) to give

$\phi(\Delta k)=\frac{3v_{e}^{2}\Delta
k}{4\omega_{p}v_{i}}\exp\left[-2\left(\frac{3v_{e}^{2}\Delta
k}{4\omega_{p}v_{i}}\right)^{2}\right],$ (12)

for $\Delta k/k\ll 1$. This peaks at ${3v_{e}^{2}\Delta
k}/{4\omega_{p}v_{i}}=1/2$, falling shortly thereafter. Thus, to marginally
resolve NLD requires a spectral resolutions in excess of

$\frac{\Delta
k}{\omega_{p}/c}=\frac{2v_{i}c}{{3}v_{e}^{2}}=\frac{2}{{3}}\sqrt{\frac{m_{e}c^{2}}{kT_{e}}\frac{m_{e}}{m_{p}}\frac{T_{i}}{T_{e}}},$
(13)

SG14 found that the temperature of the background electrons approached
relativistic temperatures (see Figure 3 of SG14) in their simulations and
hence $\Delta k\approx 0.02\omega_{p}/c$, assuming $T_{i}=T_{e}$.

The spectral resolution of a numerical computation is set by the simulated
domain’s physical size. In SG14 the reported size is $128\,c/\omega_{P}$,
implying a spectral resolution of $\Delta k_{\rm
min}=(2\pi/128)\omega_{p}/c\approx 0.05\omega_{p}/c$, larger than the minimum
required to reolve the NLD. Therefore, the calculations described in SG14 are
unable to capture the impact of NLD even for the most optimistic case of
backward scattering. Moreover, it is likely that $T_{i}\ll T_{e}$ as the
collisions needed to maintain this equilibrium are absent and electromagnetic
interactions are inefficient because of the large mass ratio; a significantly
lowered $T_{i}$ would make this disparity even more substantial.

In addition, the high temperatures reached in SG14’s simulations suppresses
the effectiveness of NLD. To see this, we estimate the NLD term in equation
(6) for the backward scattering case discussed above. Again taking ${\bf
k}^{\prime}=-(k-\Delta k)\hat{{\bf k}}$ and integrating over $\Delta k$, we
find equation (6) becomes

$\frac{dW_{k}}{dt}\approx
2\Gamma_{k}W_{k}-\frac{W_{k}\omega_{p}}{8(2\pi)^{5/2}n_{e}m_{e}v_{e}^{2}}\frac{4\pi}{3}k^{2}\Delta
kW_{k},$ (14)

where we have approximated the $\phi({\bf k},{\bf
k}^{\prime})\approx\phi(\Delta k)\approx 1/3$ over the interval $\Delta k$
given by equation (13). Substituting $\Delta k$ by equation (13) and
$\Gamma_{k}$ by equation (3), the ratio between the second and first terms,
indicating the importance of NLD, is

$\Gamma_{k}^{-1}\frac{4\pi\omega_{p}k^{2}\Delta
kW_{k}}{48(2\pi)^{5/2}n_{e}m_{e}v_{e}^{2}}\approx
10^{-4}\frac{\omega_{p}}{kc}\frac{k^{3}W_{k}}{\gamma_{b}n_{b}m_{e}c^{2}}\left(\frac{m_{e}c^{2}}{k_{B}T}\right)^{3/2}.$
(15)

For conditions relevant to SG14’s simulations,
$k^{3}W_{k}/\gamma_{b}n_{b}m_{e}c^{2}\approx 0.1$, $k_{B}T/m_{e}c^{2}\approx
1$ and $\omega_{p}/kc\approx 1$, implying NLD is suppressed by five orders of
magnitude compared to linear growth.

It remains unclear if quasilinear relaxation plays an important role. In the
simulations performed by SG14 no more than 10% of the energy of the original
beam is drained by the oblique instability, which they attribute to
quasilinear relaxation processes. Assuming this to be the case, extrapolating
to the parameters of intergalactic TeV-driven beams results in even more
stringent limits on the efficiency of plasma instabilities. However, the
accuracy and applicability of these extrapolations, from $\gamma\approx
10^{2}$ and $n_{b}/n_{e}\approx 10^{-2}$ to $\gamma\approx 10^{6}$–$10^{7}$
and $n_{b}/n_{e}\approx 3\times 10^{-18}$, is far from clear, and depends
critically on the identification of the physical and potentially numerical
causes of the instability saturation. An exhaustive study of these effects,
quasilinear relaxation and nonlinear wave-wave coupling in conjunction with
NLD, is left for future work.

##  5 Implications

We now compare the damping rate given by equation (8) to the current (i.e., at
$z=0$) damping rate due to inverse Compton scattering off cosmic microwave
background (CMB) photons, which is given by Broderick et al. (2012):

$\Gamma_{\rm IC}=\frac{4\sigma_{T}u_{\rm CMB}}{3m_{e}c}\gamma_{b}\approx
2.2\times 10^{-13}\left(\frac{E}{{\rm
TeV}}\right)\left(\frac{1+z}{2}\right)^{4}\,{\rm s}^{-1}.$ (16)

This sets a maximum beam density (paper I), which we can use to get the
effective maximum damping rate of the beam. For $\Gamma_{\rm NLD}=8\times
10^{-4}\Gamma_{\rm gr}$, we find

$\Gamma_{\rm NLD,max}\approx 6.9\times
10^{-12}\left(\frac{1+z}{2}\right)^{3\zeta-5.5}\left(\frac{EL_{E}}{10^{45}\,{\rm
ergs\,s}^{-1}}\right)\left(\frac{E}{\rm TeV}\right)^{2}\,{\rm s}^{-1},$ (17)

where $\zeta$ is our parameterization of the extragalactic background light
and is given by $\zeta=4.5$ for $z<1$ and 0 otherwise.

However, the above calculation is not self-consistent as we have presumed that
the beam density is limited only by inverse Compton scattering. If instead we
assume that the damping of the beam is driven by these plasma instabilities
(Broderick et al., 2012) with the damping rate given by equation (8), we find
a self-consistent effective damping rate of

$\Gamma_{\rm NLD}\approx
10^{-12}\left(\frac{1+z}{2}\right)^{(6\zeta-3)/4}\left(\frac{EL_{E}}{10^{45}\,{\rm
ergs\,s}^{-1}}\right)^{1/2}\left(\frac{E}{\rm TeV}\right)^{3/2}\,{\rm s}^{-1}$
(18)

Hence, we find that NLD while important does not appear to be sufficiently
strong to prevent the oblique instability from dominating the cooling of the
pair beam at energies $\gtrsim 0.8$ TeV. This domination is not complete at
$z=0$: a substantial fraction of the beam energy can now be lost to inverse
Compton scattering, of order 46% at 1 TeV and 20% at 10 TeV. We caution that
this value may be significantly reduced if nonlinear wave-wave interactions
reduce the effectiveness of NLD.

Figure 3 shows contours of the ratio between the damping rates due to plasma
effects in equation (18) and due to inverse Compton scattering in equation
(16). Lines denoting ratios of 0.1, 0.5, 1.0, 2.0 and 10.0 are shown and the
shaded green region denotes where inverse Compton scattering dominates plasma
effects. For larger photon energies, the dominance of plasma effects becomes
more pronounced. In addition, the dominance of plasma effects are also more
pronounced for increasing redshifts up to $z=1$. Beyond $z>1$, the constant
physical density of the extragalactic background light implies that the
stronger $(1+z)^{4}$ scaling of inverse Compton scattering will become more
and more important for the energy budget of these TeV beams.

fig3.pdf

Figure 3:  Contour plot of the ratio between damping due to plasma effects
($\Gamma_{\rm NLD}$) and damping due to inverse Compton scattering
($\Gamma_{\rm IC}$) as a function of redshift and photon energy, $E$, for a
blazar with equivalent isotropic luminosity, $EL_{E}=10^{45}\,{\rm
ergs\,s}^{-1}$ (black solid lines), and $10^{46}\,{\rm ergs\,s}^{-1}$ (red
dashed lines). Lines denoting ratios of $0.1$, 0.5, 1.0, 2.0 and 10.0 are
shown. The shade regions denotes where inverse Compton scattering dominates
plasma effects for $EL_{E}=10^{45}\,{\rm ergs\,s}^{-1}$ (light green) and
$10^{46}\,{\rm ergs\,s}^{-1}$ (dark green). As discussed in the text, this is
likely the “worse case” scenario for heating due to plasma effects as the
inclusion of nonlinear wave-wave damping likely improves the efficiency of the
“oblique” instability.

The limits derived above assume that NLD alone limits the growth rate of the
linear “oblique” instability. We have already noted that this neglects
nonlinear wave-wave coupling and quasilinear effects. However, the modulation
instability may also play an important role in limiting the growth of the
long-wavelength modes that are fed by NLD. The modulation instability is a
result of the ponderomotive force that results from a spatially inhomogeneous
distribution of Langmuir waves. In particular, it can be shown that the change
on the time-averaged electron density is $\delta
n_{e}\propto-|E|^{2}\propto-W$, i.e., regions of high wave density correspond
to regions with lower electron density (Boyd & Sanderson, 2003). Such a shift
in the electron density leads to a shift in the plasma frequency that is of
order

$\frac{\delta\omega_{p}}{\omega_{p}}=\frac{\delta n_{e}}{2n_{e}}.$ (19)

As a result sufficiently long-wavelength Langmuir waves that propagate freely
in regions of high wave density (low electron density) can lie below the
plasma resonance outside, and therefore are trapped. Which modes become
trapped depend upon the shift in the plasma frequency and therefore the wave
density. Assuming the Langmuir wave dispersion relation,
$\omega(k)=\omega_{p}(1+3k^{2}\lambda_{D}^{2}/2$), where
$\lambda_{D}=v_{e}/\omega_{p}$ is the Debye length, the condition for trapped
modes is given by

$\frac{W}{n_{e}k_{B}T}\gg 3(k\lambda_{D})^{2}.$ (20)

The late-time development of the modulation instability in multiple dimensions
is currently believed to result in strong turbulence, which then rapidly damps
the participating Langmuir waves and results in direct heating of the
background plasma (ME13, Schlickeiser et al. 2012).

For the long wavelength modes generated by NLD, $k\approx
10^{-2}\omega_{p}/c$, equation (20) implies

$\frac{f_{w}n_{b}\gamma_{b}m_{e}c^{2}}{n_{e}k_{B}T}\gg\left(\frac{3\times
10^{-2}v_{e}}{c}\right)^{2}\approx 6\times 10^{-10}\left(\frac{T}{10^{4}\,{\rm
K}}\right)$ (21)

where $f_{w}$ is the saturation amplitude of these waves relative to the beam
energy $n_{b}\gamma_{b}m_{e}c^{2}$. For $f_{w}=10^{-3}$ given by the
saturation amplitude due to NLD, this gives

$\frac{f_{w}n_{b}\gamma_{b}m_{e}c^{2}}{n_{e}k_{B}T}\approx 2\times
10^{-9}\left(\frac{f_{w}}{10^{-3}}\right)\left(\frac{n_{b}/n_{e}}{3\times
10^{-18}}\right)\left(\frac{\gamma_{b}}{10^{6}}\right)\left(\frac{T}{10^{4}\,{\rm
K}}\right)^{-1},$ (22)

which is above the criterion given in equation (21) and can allow these waves
to rapidly heat the background electrons (Schlickeiser et al., 2012). In doing
so, the modulation instability may limit the effectiveness with which NLD can
drive long-wavelength modes, and hence the linear growth of the instability,
as well as provide a natural mechanism for converting the wave energy into
heat. However, more work is required on this question.

##  6 Conclusions

We have calculated the effects of NLD on the saturation amplitude of the
“oblique” instability to which high energy pair beams in the IGM are linearly
unstable to. Using a numerical calculation, we find that the “oblique”
instability remains the most powerful cooling mechanism for these pair beams
contrary to the earlier claims of ME13. In particular, we find that the beam
saturates at a rate that is $\approx 0.1\%$ of kinetic growth rate of the
“oblique” instability. The damping of the beam leads to the transfer of beam
energy into long-wavelength non-resonant waves. When comparing to the estimate
of ME13, we find that our damping rate exceeds their estimate by two orders of
magnitude. Using this damping rate, we conclude that the oblique instability
is effective in quenching the beam at $z\approx 1$, but it less effective at
different redshifts.

We caution that the results that we present here are limited to the effects of
NLD. As we argue above, the inclusion of nonlinear wave-wave coupling will
likely lead to a more equitable distribution of energy among wavevectors. The
nonlinear damping rate of the “oblique” instability is likely to increase
under these conditions. Thus, the calculation that we present here represents
the “worst-case” scenario for the “oblique” instability where the magnitude of
plasma effects is comparable to the effects of inverse Compton scattering.

###### Acknowledgements.

P.C. and A.L. gratefully acknowledges support from the NASA ATP program
through NASA grant NNX13AH43G, and the NSF through grant AST-1255469. A.E.B.
and M.S. receive financial support from the Perimeter Institute for
Theoretical Physics and the Natural Sciences and Engineering Research Council
of Canada through a Discovery Grant. Research at Perimeter Institute is
supported by the Government of Canada through Industry Canada and by the
Province of Ontario through the Ministry of Research and Innovation. C.P.
gratefully acknowledges financial support of the Klaus Tschira Foundation.
E.P. acknowledges support by the ERC grant “The Emergence of Structure during
the epoch of Reionization”.

## References

  * Abazajian et al. (2011) Abazajian, K. N., Blanchet, S., & Harding, J. P. 2011, Phys. Rev. D, 84, 103007
  * Ackermann et al. (2012) Ackermann, M., et al. 2012, Science, 338, 1190
  * Boyd & Sanderson (2003) Boyd, T. J. M. and Sanderson, J. J., 2003, The Physics of Plasmas.
  * Bret et al. (2010) Bret, A., Gremillet, L., & Dieckmann, M. E. 2010, Physics of Plasmas, 17, 120501
  * Breǐzman et al. (1972) Breǐzman, B. N., Ryutov, D. D., & Chebotaev, P. Z. 1972, Soviet Journal of Experimental and Theoretical Physics, 35, 741
  * Broderick et al. (2012) Broderick, A. E., Chang, P., & Pfrommer, C. 2012, ApJ, 752, 22
  * Broderick et al. (2013a) Broderick, A. E., Pfrommer, C., Puchwein, E., & Chang, P. 2014a, ApJ, in press, ArXiv:1308.0340
  * Broderick et al. (2013b) —. 2014b, ApJ, 790, 137
  * Chang et al. (2012) Chang, P., Broderick, A. E., & Pfrommer, C. 2012, ApJ, 752, 23
  * Davidson (1972) Davidson, R. C. 1972, Methods in nonlinear plasma theory.
  * Dermer et al. (2011) Dermer, C. D., Cavadini, M., Razzaque, S., Finke, J. D., Chiang, J., & Lott, B. 2011, ApJ, 733, L21
  * Dolag et al. (2011) Dolag, K., Kachelriess, M., Ostapchenko, S., & Tomàs, R. 2011, ApJ, 727, L4
  * Domínguez et al. (2013) Domínguez, A., Finke, J. D., Prada, F., Primack, J. R., Kitaura, F. S., Siana, B., & Paneque, D. 2013, ApJ, 770, 77
  * Gould & Schréder (1967) Gould, R. J., & Schréder, G. P. 1967, Physical Review, 155, 1408
  * Inoue & Ioka (2012) Inoue, Y., & Ioka, K. 2012, Phys. Rev. D, 86, 023003
  * Kaplan & Tsytovich (1968) Kaplan, S. A., & Tsytovich, V. N. 1968, Soviet Ast., 11, 956
  * Kneiske & Mannheim (2008) Kneiske, T. M., & Mannheim, K. 2008, A&A, 479, 41
  * Lu et al. (2013) Lu, Z., Mo, H., Lu, Y., Katz, N., Weinberg, M. D., van den Bosch, F. C., & Yang, X. 2014, MNRAS, 439, 1294
  * Melrose (1986) Melrose, D. B. 1986, Instabilities in Space and Laboratory Plasmas
  * Miniati & Elyiv (2013) Miniati, F., & Elyiv, A. 2013, ApJ, 770, 54
  * Neronov & Semikoz (2009) Neronov, A., & Semikoz, D. V. 2009, Phys. Rev. D, 80, 123012
  * Neronov & Vovk (2010) Neronov, A., & Vovk, I. 2010, Science, 328, 73
  * Pfrommer et al. (2012) Pfrommer, C., Chang, P., & Broderick, A. E. 2012, ApJ, 752, 24
  * Puchwein et al. (2012) Puchwein, E., Pfrommer, C., Springel, V., Broderick, A. E., & Chang, P. 2012, MNRAS, 423, 149
  * Salamon & Stecker (1998) Salamon, M. H., & Stecker, F. W. 1998, ApJ, 493, 547
  * Schlickeiser et al. (2012) Schlickeiser, R., Ibscher, D., & Supsar, M. 2012, ApJ, 758, 102
  * Schlickeiser et al. (2013) Schlickeiser, R., Krakau, S., & Supsar, M. 2013, ApJ, 777, 49
  * Sironi & Giannios (2014) Sironi, L., & Giannios, D. 2014, ApJ, 787, 49
  * Smith & Fung (1971) Smith, D. F., & Fung, P. C. W. 1971, Journal of Plasma Physics, 5, 1
  * Takahashi et al. (2012) Takahashi, K., Mori, M., Ichiki, K., & Inoue, S. 2012, ApJ, 744, L7
  * Tavecchio et al. (2011) Tavecchio, F., Ghisellini, G., Bonnoli, G., & Foschini, L. 2011, MNRAS, 414, 3566
  * Tavecchio et al. (2010) Tavecchio, F., Ghisellini, G., Foschini, L., Bonnoli, G., Ghirlanda, G., & Coppi, P. 2010, MNRAS, 406, L70
  * Taylor et al. (2011) Taylor, A. M., Vovk, I., & Neronov, A. 2011, A&A, 529, A144
  * Venters (2010) Venters, T. M. 2010, ApJ, 710, 1530
  * Ziebell et al. (2008a) Ziebell, L. F., Gaelzer, R., Pavan, J., & Yoon, P. H. 2008a, Plasma Physics and Controlled Fusion, 50, 085011
  * Ziebell et al. (2008b) Ziebell, L. F., Gaelzer, R., & Yoon, P. H. 2008b, Physics of Plasmas, 15, 032303


mm

11institutetext:  School of Engineering, University of Edinburgh, Edinburgh
EH9 3FB, United Kingdom
11email<EMAIL_ADDRESS>22institutetext: Canon Medical Research Europe Ltd., Edinburgh EH6 5NP, United
Kingdom
22email<EMAIL_ADDRESS>
# Group Distributionally Robust Knowledge Distillation

Konstantinos Vilouras 11 Xiao Liu 1122 Pedro Sanchez 11 Alison Q. O’Neil 1122
Sotirios A. Tsaftaris 11

###### Abstract

Knowledge distillation enables fast and effective transfer of features learned
from a bigger model to a smaller one. However, distillation objectives are
susceptible to _sub-population shifts_ , a common scenario in medical imaging
analysis which refers to groups/domains of data that are underrepresented in
the training set. For instance, training models on health data acquired from
multiple scanners or hospitals can yield subpar performance for minority
groups. In this paper, inspired by distributionally robust optimization (DRO)
techniques, we address this shortcoming by proposing a group-aware
distillation loss. During optimization, a set of weights is updated based on
the per-group losses at a given iteration. This way, our method can
dynamically focus on groups that have low performance during training. We
empirically validate our method, GroupDistil on two benchmark datasets
(natural images and cardiac MRIs) and show consistent improvement in terms of
worst-group accuracy.

###### Keywords:

Invariance Knowledge Distillation Sub-population Shift Classification

##  1 Introduction

The rapid success of deep learning can be largely attributed to the
availability of both large-scale training datasets and high-capacity networks
able to learn arbitrarily complex features to solve the task at hand. Recent
practices, however, pose additional challenges in terms of real-world
deployment due to increased complexity and computational demands. Therefore,
developing lightweight versions of deep models without compromising
performance remains an active area of research.

Knowledge distillation To this end, knowledge distillation is a promising
technique aiming to guide the learning process of a small model (student)
using a larger pre-trained model (teacher). This guidance can be exhibited in
various forms; for example, Hinton et al. [3] propose to match the soft class
probability distributions predicted by the teacher and the student. Tian et
al. [10] augment this framework with a contrastive objective to retain the
dependencies between different output dimensions (since the loss introduced in
[3] operates on each dimension independently). On the contrary, Romero et al.
[7] consider the case of matching intermediate teacher and student features
(after projection to a shared space since the original dimensions might be
mismatched). Another interesting application involves distilling an ensemble
of teachers into a single student model, which has been previously explored in
the case of histopathology [9] and retinal [4] images.

Distillation may lead to poor performance Although knowledge distillation has
been widely used as a standard approach to model compression, our
understanding of its underlying mechanisms remains underexplored. Recently,
Ojha et al. [6] showed, through a series of experiments, that a distilled
student exhibits both useful properties such as invariance to data
transformations, as well as biases, e.g., poor performance on rare subgroups
inherited from the teacher model. Here, we attempt to tackle the latter issue
which can be especially problematic, for instance, for medical tasks where we
have access to only a few data for certain groups in the training set. Note
that there also exist concurrent works [11] similar to ours that consider the
case of distillation on long-tailed111i.e., on sets where the majority of data
comes from only a few classes datasets.

In this paper, we consider the task of knowledge distillation under a common
type of distribution shift called _sub-population shift_ where both training
and test sets overlap, yet per-group proportions differ between sets. In sum,
our contributions are the following:

  * •

We propose a simple, yet effective, method that incorporates both the original
distillation objective and also group-specific weights that allow the student
to achieve high accuracy even on minority groups.

  * •

We evaluate our method on two publicly available datasets, i.e., on natural
images (Waterbirds) and cardiac MRI data (M&Ms) and show improvements in terms
of worst-group accuracy over the commonly used knowledge distillation loss of
[3].

##  2 Methodology

Preliminaries Let $D$ denote the data distribution from which triplets of data
instances $x\in\mathcal{X}$, labels $y\in\mathcal{Y}$ and domains
$d\in\mathcal{D}$ are sampled, respectively. Let also the teacher
$f_{T}:\mathcal{X}\rightarrow\mathcal{Z}_{T}$ and student
$f_{S}:\mathcal{X}\rightarrow\mathcal{Z}_{S}$ map input samples to class
logits. The vanilla knowledge distillation loss introduced in Hinton et al.
[3] is defined as

$\mathcal{L}_{KD}=(1-\alpha)\cdot H(y,\sigma(z_{S}))+\alpha\tau^{2}\cdot
D_{KL}(\sigma(z_{T}/\tau),\sigma(z_{S}/\tau)),$ (1)

where $H(p,q)=-\mathbb{E}_{p}[\log q]$ refers to cross-entropy,
$D_{KL}(p,q)=H(p,q)-H(p)$ is the Kullback-Leibler divergence metric, $\sigma$
is the softmax function, $\tau$ is the temperature hyperparameter and $\alpha$
the weight that balances the loss terms.

In our case, inspired by the groupDRO algorithm of Sagawa et al. [8], we
propose an alternative objective that incorporates group-specific weights
(termed GroupDistil). An overview of our method is presented in Figure 1. The
final form of the proposed loss function is shown in Equation 2, where $w_{d}$
and $\mathcal{L}_{KD}^{d}$ refer to the weight and the dedicated distillation
loss for group $d$, respectively. This allows the student to emphasize domains
where performance remains low (e.g., rare subgroups within the training set)
by upweighting their contribution to the overall loss. In that sense, the
final optimization objective can be divided into two steps, i.e., first
performing exponentiated gradient ascent on group weights and then minibatch
stochastic gradient descent on student’s weights. Algorithm 1 provides a full
description of the proposed method. The main difference with the groupDRO
algorithm [8] is the choice of loss function, i.e., we use the distillation
loss instead of cross-entropy, and it is highlighted in Algorithm 1 in blue.

$\mathcal{L}_{GroupDistil}=\sum_{d=1}^{\mathcal{|D|}}w_{d}\cdot\mathcal{L}_{KD}^{d},$
(2)

Figure 1: Overview of the proposed method GroupDistil. Data from multiple
domains $d=1,...,|\mathcal{D}|$ are fed to both the teacher and student to
produce class logits. A soft version of the probability distribution over
classes induced by both models is used to calculate the distillation loss,
whereas the hard version ($\tau=1$) of the student’s probability distribution
is used for the classification loss. The final loss is calculated as the dot
product between the per-group losses and the group weights. This way, more
emphasis is placed on domains where the student has low performance.

Input: Data distribution for a given domain $P_{d}$; Learning rate
$\eta_{\theta}$; Group weight step size $\eta_{w}$

Initialize (uniform) group weights $w^{(0)}$ and student weights
$\theta_{s}^{(0)}$

for  _$t=1,...,T$_ do

$d\sim\mathcal{U}(1,...,|\mathcal{D}|)$

// Draw a random domain

$(x,y)\sim P_{d}$

// Sample (data, labels) from domain $d$

$\ell\leftarrow\mathcal{L}_{KD}^{d}(x,y)$

// Calculate distillation loss (Eqn 1)

$w^{\prime}\leftarrow w^{(t-1)}$; $w^{\prime}_{d}\leftarrow
w^{\prime}_{d}\cdot\exp(\eta_{w}\ell)$

// Update weight for domain $d$

$w^{(t)}\leftarrow w^{\prime}/\sum_{d^{\prime}}w^{\prime}_{d^{\prime}}$

// Normalize new weights

$\theta_{s}^{(t)}\leftarrow\theta_{s}^{(t-1)}-\eta_{\theta}w_{d}^{(t)}\nabla\ell$

// Optimize student weights via SGD

end for

Algorithm 1 _GroupDistil_ method

##  3 Datasets

###  3.1 Waterbirds

Waterbirds [5, 8] is a popular benchmark used to study the effect of sub-
population shifts. More specifically, it contains instances from all possible
combinations of $\mathcal{Y}=$ {landbird, waterbird} labels and $\mathcal{D}=$
{land, water} backgrounds (4 domains in total), respectively. However, it is
collected in a way such that uncommon pairs (i.e., landbirds on water and
waterbirds on land) occur less frequently, thus creating an imbalance. Also,
note that the level of imbalance in the training set is different than that of
the test set. Training, validation and test sets consist of 4795, 1199 and
5794 samples, respectively. In our setup, we first resize images to fixed
spatial dimensions ($256\times 256$), then extract the $224\times 224$ center
crop and normalize using Imagenet’s mean and standard deviation. A few
representative samples of this dataset are depicted in Figure 2.

Figure 2: Examples of images from Waterbirds dataset. Left: waterbird on
water. Middle left: waterbird on land. Middle right: landbird on land. Right:
landbird on water.

###  3.2 M&Ms

The multi-centre, multi-vendor and multi-disease cardiac image segmentation
(M$\&$Ms) dataset [1] contains 320 subjects. Subjects were scanned at 6
clinical centres in 3 different countries using 4 different magnetic resonance
scanner vendors (Siemens, Philips, GE, and Canon) i.e., domains A, B, C and D.
For each subject, only the end-systole and end-diastole phases are annotated.
Voxel resolutions range from $0.85\times 0.85\times 10$ mm to $1.45\times
1.45\times 9.9$ mm. Domain A contains 95 subjects. Domain B contains 125
subjects. Both domains C and D contain 50 subjects. We show example images of
the M$\&$Ms data in Fig. 3. Note that, while the dataset was originally
collected for the task of segmentation, we instead use it here for
classification.

Figure 3: Examples of cardiac MRI images from M$\&$Ms dataset [1].

##  4 Experimental Results

We now discuss our results. First, we briefly discuss hyperparameter choices
made on a per-dataset basis. In each case, we compare the performance of three
types of _students_ : an individual student that was trained with a group
robustness method [8] from scratch (_groupDRO_ hereafter), a distilled student
trained with the objective of Equation 1 (_KD_ hereafter) and also a distilled
student trained with our proposed _GroupDistil_ method as described in
Equation 2 and Algorithm 1 (_GroupDistil_ hereafter).

Waterbirds We consider the case of distilling the knowledge from a ResNet-50
teacher model pre-trained with groupDRO to a ResNet-18 student (shorthanded as
R$50\rightarrow$ R$18$). Note that our proposed method does not impose any
restrictions on either the choice of model architectures or the teacher’s pre-
training strategy (groupDRO was merely chosen to ensure that the teacher
learns robust features on the given training set). For our objective, we use
the following hyperparameters: batch size $=128$, total number of epochs
$=30$, Adam optimizer with learning rate $\eta_{\theta}=10^{-4}$, temperature
$T=4$, $\alpha=0.9$ and group weight step size $\eta_{w}=0.01$. Final results
on the official test set are depicted in Table 1.

Table 1: Distillation results on Waterbirds dataset. First two rows show performance of each model trained from scratch. The last two rows refer to vanilla knowledge distillation (KD) and our proposed method (GroupDistil), respectively. In those cases, we fix the teacher (1st row) and initialize each R18 model with ImageNet-pretrained weights. We report both the average (2nd column) and worst-group (3rd column) accuracy on the official test set. Final results have been averaged across 3 random seeds. Setup | Adjusted avg acc. (%) | Worst-group acc. (%)
---|---|---
groupDRO (R50) | 90.9 | 86.8
groupDRO (R18) | 86.0 | 79.6
KD (R$50\rightarrow$ R$18$) | 95.2 $\pm$ 0.1 | 78.0 $\pm$ 1.5
GroupDistil (R$50\rightarrow$ R$18$) | 86.6 $\pm$ 0.8 |  81.7 $\pm$ 0.8

Our _GroupDistil_ method shows consistent improvement ($>$ 2%) compared to the
knowledge distillation objective (KD) of [3] in terms of worst-group accuracy.
This implies that our distilled student does not rely on spurious attributes
(background information) to classify each type of bird. This is also evident
from the fact that worst-group accuracy in our case does not deviate much
(approx. 5%) from the average accuracy in the test set. Note that our
distilled student even outperforms a ResNet-18 model trained from scratch
using groupDRO, showing that smaller models can largely benefit from larger
ones in this setup.

M&Ms In the case of M&Ms, since there is no official benchmark for studying
sub-population shifts, we divide the available data in training and test sets
as follows: First, we consider a binary classification task using the two most
common patient state classes222This choice eliminates Domain B due to lack of
available data., i.e., hypertrophic cardiomyopathy (HCM) and healthy (NOR).
Then, we conduct our experiments using two types of splits as presented in
Table 2. Splits are defined in such a way that allows us to measure
performance on a specific type of scanner (domain) that is underrepresented in
the training set. Thus, for the first (resp. second) split, we keep only 2
patients per class from Domain C (resp. D) in the train set, and the rest in
the test set. Note that, for fair comparison, we ensure that the test set is
always balanced.

Table 2: M&Ms dataset. For each type of experiment, we show the total number of patients per dataset (train or test), domain (A, C, or D) and class (HCM or NOR), respectively. Note that we extract 20 2D frames from each patient for both training and testing. |  | Train | Test
---|---|---|---
Experiment | Class | Domains (A/C/D) | Domains (A/C/D)
Split 1 | HCM | 25/2/10 | -/3/-
NOR | 21/2/14 | -/3/-
Split 2 | HCM | 25/5/2 | -/-/8
NOR | 21/11/2 | -/-/8

For this dataset, we use a ResNet-18 (pre-trained with groupDRO) as teacher
and a DenseNet-121 as a student model (R$18\rightarrow$ D$121$ setup). As a
pre-processing step, we first extract a random 2D frame from each 4D input (in
total, we extract 20 frames from each patient per epoch) and then apply data
augmentations such as intensity cropping, random rotation and flip. The input
for each model is a $224\times 224$ crop from each frame. We also used the
following hyperparameters: batch size $=128$, total number of epochs $=5$,
Adam optimizer with learning rate $\eta_{\theta}=5\cdot 10^{-4}$, temperature
$T=4$, $\alpha=0.9$ and group weight step size $\eta_{w}=0.01$. Results for
the first split are shown in Table 3, whereas for the second split in Table 4.

As in Waterbirds, similar observations can be made for the M&Ms dataset. It is
clear that in both types of splits, our method reaches the highest accuracy.
Also note that the teacher and the student trained from scratch have low
performance on the test set, indicating the challenging nature of this
dataset; yet, distilled students can significantly outperform them. The
results for KD method exhibit high variance (possibly due to the limited
number of available data), indicating that it could be unstable in this setup.
On the contrary, our method remains fairly robust and shows consistent
performance improvements over the rest of the methods.

Table 3: Distillation results on M&Ms dataset, 1st split (testing on Domain C). First two rows show performance of each model trained from scratch. The last two rows refer to vanilla knowledge distillation (KD) and our proposed method (GroupDistil), respectively. We report the average accuracy on the test set. Final results have been averaged across 5 random seeds. Setup | Accuracy (%)
---|---
groupDRO (R18) | 53.3
groupDRO (D121) | 46.7
KD (R$18\rightarrow$ D$121$) | 55.8 $\pm$ 7.0
GroupDistil (R$18\rightarrow$ D$121$) |  62.8 $\pm$ 4.9
Table 4: Distillation results on M&Ms dataset, 2nd split (testing on Domain D). First two rows show performance of each model trained from scratch. The last two rows refer to vanilla knowledge distillation (KD) and our proposed method (GroupDistil), respectively. We report the average accuracy on the test set. Final results have been averaged across 5 random seeds. Setup | Accuracy (%)
---|---
groupDRO (R18) | 62.2
groupDRO (D121) | 56.2
KD (R$18\rightarrow$ D$121$) | 64.2 $\pm$ 7.4
GroupDistil (R$18\rightarrow$ D$121$) |  66.5 $\pm$ 2.8

##  5 Conclusion and Future Work

In this paper we consider the task of knowledge distillation under a
challenging type of distribution shift, i.e., sub-population shift. We showed
that adding group-specific weights to a popular distillation objective
provides a significant boost in performance, which even outperforms the same
student architecture trained from scratch with a group robustness method in
terms of worst-group accuracy. We also made sure that our proposed method
remains fairly general, allowing arbitrary combinations of teacher-student
models.

A limitation of our work is the fact that we assume access to fully labeled
data, i.e., with both label and domain annotations, which is restrictive in
practice. Therefore, for future work, we plan to investigate methods that
infer domains directly from data as in [2, 12].

##  6 Acknowledgements

This work was supported by the University of Edinburgh, the Royal Academy of
Engineering and Canon Medical Research Europe by a PhD studentship to
Konstantinos Vilouras. S.A. Tsaftaris also acknowledges the support of Canon
Medical and the Royal Academy of Engineering and the Research Chairs and
Senior Research Fellowships scheme (grant RCSRF1819\8\25), and the UK’s
Engineering and Physical Sciences Research Council (EPSRC) support via grant
EP/X017680/1.

## References

  * [1] Campello, V.M., Gkontra, P., Izquierdo, C., Martín-Isla, C., Sojoudi, A., Full, P.M., Maier-Hein, K., Zhang, Y., He, Z., Ma, J., et al.: Multi-centre, multi-vendor and multi-disease cardiac segmentation: The m&ms challenge. IEEE TMI (2021)
  * [2] Creager, E., Jacobsen, J.H., Zemel, R.: Environment inference for invariant learning. In: International Conference on Machine Learning. pp. 2189–2200. PMLR (2021)
  * [3] Hinton, G., Vinyals, O., Dean, J.: Distilling the knowledge in a neural network. arXiv preprint arXiv:1503.02531 (2015)
  * [4] Ju, L., Wang, X., Wang, L., Liu, T., Zhao, X., Drummond, T., Mahapatra, D., Ge, Z.: Relational subsets knowledge distillation for long-tailed retinal diseases recognition. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2021: 24th International Conference, Strasbourg, France, September 27–October 1, 2021, Proceedings, Part VIII 24. pp. 3–12. Springer (2021)
  * [5] Koh, P.W., Sagawa, S., Marklund, H., Xie, S.M., Zhang, M., Balsubramani, A., Hu, W., Yasunaga, M., Phillips, R.L., Gao, I., et al.: Wilds: A benchmark of in-the-wild distribution shifts. In: International Conference on Machine Learning. pp. 5637–5664. PMLR (2021)
  * [6] Ojha, U., Li, Y., Lee, Y.J.: What knowledge gets distilled in knowledge distillation? arXiv preprint arXiv:2205.16004 (2022)
  * [7] Romero, A., Ballas, N., Kahou, S.E., Chassang, A., Gatta, C., Bengio, Y.: Fitnets: Hints for thin deep nets. arXiv preprint arXiv:1412.6550 (2014)
  * [8] Sagawa, S., Koh, P.W., Hashimoto, T.B., Liang, P.: Distributionally robust neural networks for group shifts: On the importance of regularization for worst-case generalization. arXiv preprint arXiv:1911.08731 (2019)
  * [9] Tellez, D., Balkenhol, M., Otte-Höller, I., van de Loo, R., Vogels, R., Bult, P., Wauters, C., Vreuls, W., Mol, S., Karssemeijer, N., et al.: Whole-slide mitosis detection in h&e breast histology using phh3 as a reference to train distilled stain-invariant convolutional networks. IEEE transactions on medical imaging 37(9), 2126–2136 (2018)
  * [10] Tian, Y., Krishnan, D., Isola, P.: Contrastive representation distillation. arXiv preprint arXiv:1910.10699 (2019)
  * [11] Wang, S., Narasimhan, H., Zhou, Y., Hooker, S., Lukasik, M., Menon, A.K.: Robust distillation for worst-class performance. arXiv preprint arXiv:2206.06479 (2022)
  * [12] Zhang, M., Sohoni, N.S., Zhang, H.R., Finn, C., Ré, C.: Correct-n-contrast: A contrastive approach for improving robustness to spurious correlations. arXiv preprint arXiv:2203.01517 (2022)