# Wave Propagation at Oblique Shocks: How Did Tycho Get Its Stripes?

J. Martin Laming11affiliation: Space Science Division, Naval Research
Laboratory, Code 7684, Washington DC 20375<EMAIL_ADDRESS>
###### Abstract

We describe a new model for the “stripes” of synchrotron radiation seen in the
remnant of Tycho’s supernova. In our picture, cosmic rays streaming ahead of
the forward shock generate parallel propagating (with respect to the local
magnetic field direction) circularly polarized Alfvén waves that are almost
free of dissipation, and due to being circularly polarized exhibit no spatial
variation of magnetic field strength. Following interaction with the SNR shock
with nonzero obliquity, these parallel propagating waves become obliquely
propagating, due the the wave refraction (different in principle for the
different plane wave components), and dissipation sets in. The magnetosonic
polarization decays faster, due to transit time damping, leaving only the
Alfvén mode. This surviving mode now exhibits a spatial variation of the
magnetic field, leading to local maxima and minima in the synchrotron
emission, i.e. the stripes. We attribute the initial wave generation to the
Bell instability, which in contrast to the resonant generation of upstream
Alfvén waves, gives rise to a preferred wavelength, and hence the single wave
period at which the stripes are seen. Based on estimates for damping rates due
to turbulent cascade and transit time damping, we estimate the dependence of
the visibility of the stripes on the shock obliquity, and determine a maximum
cosmic ray energy in Tycho’s SNR in the range $6\times 10^{14}-1\times
10^{15}$ eV.

###### Subject headings:

acceleration of particles — cosmic rays — magnetic fields — shock waves — ISM:
supernova remnants

††slugcomment: Draft

##  1\. Introduction

The advent of the Chandra and XMM-Newton X-ray astronomy missions has
revitalized the field of cosmic ray acceleration. Chandra’s combination of
arcsecond resolution imaging and the moderate energy spectral resolution
afforded by CCD X-ray detectors has proven well suited to the study of
supernova remnants, and specifically the X-ray synchrotron radiation emitted
by cosmic ray electrons. This often reveals itself as thin rims of continuum
emission in the 4 - 6 keV waveband (a region essentially free of quasi-thermal
line emission from shocked plasma) spatially coinciding with the supernova
remnant forward shock. This can be understood (e.g. Vink & Laming, 2003) as
cosmic ray electrons radiating in the strong (i.e. amplified) magnetic field
at the shock. The radial extent is limited by either the radiative loss time
of the electrons, or the decay of magnetic field postshock. The remnant of
Tycho’s supernova shows an even more intriguing pattern of synchrotron
radiation; a series of “stripes” or ripples, (Eriksen et al., 2011), observed
most clearly in the 4 - 6 keV waveband where line emission from the quasi-
thermal plasma in the SNR is absent. The wavelength corresponds to the
gyroradius of cosmic ray protons with energies in the range $10^{14}-10^{15}$
eV, but the precise mechanism by which these structures form is not known.

The Bell instability (Bell, 2004, 2005), which amplifies magnetic field in the
shock precursor region, generates linearly polarized structures in a near
perpendicular shock geometry. Considering the effect of these short wavelength
fluctuations on the cosmic ray current, Vladimirov et al. (2009) and Bykov et
al. (2011a) have shown that long wavelength upstream structures can result,
with spatial variations of the magnetic field strength and hence synchrotron
emissivity. Bykov et al. (2011b) argue that these long wavelength structures
are responsible for the “stripes”. A number of conditions must be met. Most
importantly, the shock region where the stripes appear must be “nearly
perpendicular” (Bykov et al., 2011b, don’t specify how close to 90∘ they
require), and that in this nearly perpendicular region, shock acceleration
must still be efficient. But as discussed elsewhere (Zank et al., 2006; Laming
et al., 2013), the efficiency of shock acceleration at quasi-perpendicular
shocks is open to question.

Malkov et al. (2012) offer an alternative idea, that is appears most promising
at parallel shocks. They find fully nonlinear exact ideal MHD solutions
supported by the cosmic ray return current, which comprise pulses of Alfvén
waves that can propagate ahead of the main shock. Malkov et al. (2012) argue
that these Alfvén pulses when visible in X-rays should appear as quasi-
periodic stripes, with spacing similar to that of the observed stripes. In
their equations, Malkov et al. (2012) neglect thermal and cosmic ray pressure
gradients, assuming that the ponderomotive force of the turbulence is much
stronger. In the opposite limit high frequency sound waves are generated by
the Drury instability (Drury & Falle, 1986). Caprioli & Spitkovsky (2013)
suggest that the forward shock may push forwards in upstream cavities created
at the saturation of the Bell instability, also at parallel shocks. Rakowski
et al. (2011) discuss a similar idea in connection with shock structures
(although not stripes) seen in SN 1006.

In this paper we pursue a different model for these synchrotron stripes in
Tycho’s SNR. Cosmic rays drifting ahead of a quasi-parallel shock generate
upstream circularly polarized Alfvén waves. We assume these to be parallel
propagating (see e.g. Bamert et al., 2004; Gargaté & Spitkovsky, 2012), and
with no spatial variation in magnetic field pressure. Upon passage through the
shock, both the magnetic field and the wave propagation change direction, and
by different amounts, so that the formerly parallel propagating waves are now
obliquely propagating. The circular polarization decomposes into its linearly
polarized constituents. The magnetosonic polarization becomes compressive, and
quickly decays by transit time damping. The Alfvén polarization survives
longer, and its spatially varying magnetic field gives rise the regions of
enhanced synchrotron emission, observable as the stripes or ripples seen by
Eriksen et al. (2011). The following sections treat the wave refraction at the
shock, the wave transmission and reflection coefficients, and the postshock
damping of the magnetosonic and Alfveń modes. Section 5 puts these topics
together to explain the origin of the stripes, and to derive physical
implications from this identification.

Figure 1.— Schematic diagram of wave propagation at the oblique shock, in the
shock rest frame. Upstream (left), plasma moves towards the shock front with
velocity ${\bf u}_{1}$, carrying with it magnetic field ${\bf B}_{1}$ at angle
$\theta_{1}$ to ${\bf u}_{1}$ (or the shock normal). Postshock, the magnetic
field is ${\bf B}_{2}$ at angle $\theta_{2}$ to the shock normal and the flow
velocity is ${\bf u}_{2}$ at angle $\alpha$ to the shock normal. In the
preshock medium, waves are assumed to be parallel propagating along ${\bf
B}_{1}$. Postshock, the waves refract and travel at an angle $\theta_{2W}$ to
the shock normal. They are no longer parallel propagating. In principle,
Alfvén and magnetosonic polarizations will refract at slightly different
angles. The difference between these two angles of refraction tends to zero as
the Alfvén Mach number $M_{A}\rightarrow\infty$, the approximation taken in
this paper.

##  2\. Wave Refraction at Shocks

McKenzie & Westphal (1970) and Webb et al. (1999) give formal accounts of wave
properties at oblique shocks. We follow and extend slightly for our particular
application the more pedagogical treatment of Achterberg & Blandford (1986). A
schematic diagram of the wave propagation at an oblique shock is shown in Fig
1. In the preshock medium, we take both Alfvén and magnetosonic polarizations
to be parallel propagating (with respect to the upstream magnetic field) with
phase velocity $V_{A1}$, the upstream Alfvén speed. Downstream, the Alfvén
polarization has phase speed $V_{A2}\cos\left(\theta_{2}-\theta_{2W}\right)$,
where $\theta_{2}$ is the angle from the shock normal to the postshock
magnetic field, and $\theta_{2W}$ is the angle from the shock normal to the
wave propagation direction. The phase speed of the magnetosonic wave depends
on the plasma beta (the ratio of gas pressure to magnetic pressure), and is
given by (e.g Melrose, 1986)

$\displaystyle v^{2}$ $\displaystyle=$ $\displaystyle{1\over
2}\bigg{(}V_{A2}^{2}+c_{s2}^{2}\pm$ (2)
$\displaystyle\left[\left(V_{A2}^{2}+c_{s2}^{2}\right)^{2}-4V_{A2}^{2}c_{s2}^{2}\cos^{2}\left(\theta_{2}-\theta_{2W}\right)\right]^{1/2}\bigg{)}$
$\displaystyle\simeq$ $\displaystyle
V_{A2}^{2}\cos^{2}\left(\theta_{2}-\theta_{2W}\right)-{V_{A2}^{4}\over
c_{s2}^{2}}\big{[}\cos^{2}\left(\theta_{2}-\theta_{2W}\right)$
$\displaystyle+\cos^{4}\left(\theta_{2}-\theta_{2W}\right)\big{]}+\ldots\quad\beta>>1$
$\displaystyle\simeq
V_{A2}^{2}+c_{s2}^{2}\sin^{2}\left(\theta_{2}-\theta_{2W}\right)\quad\beta<1.$

In the following we will only consider the leading terms in equations 1 and 2,
which will lead to the same angle of refraction $\theta_{2W}$ for both
polarizations. In principle, however, the different polarizations can refract
at slightly different angles. At the shock transition, we demand that the
phases of upstream and downstream waves are equal, i.e. $\omega_{1}t-{\bf
k}_{1}\cdot{\bf r}=\omega_{2}t-{\bf k}_{2}\cdot{\bf r}$, and that the
components of wavevectors perpendicular to the shock velocity are continuous.
For magnetosonic waves at low $\beta$ these constraints give the equations

$\displaystyle
k_{1}V_{A1}-k_{1}\cos\theta_{1}u_{1}=k_{2}V_{A2}-k_{2}u_{2}\cos\left(\theta_{2W}-\alpha\right)$
(3) $\displaystyle k_{1}\sin\theta_{1}=k_{2}\sin\theta_{2W},$ (4)

where $\tan\alpha=u_{2\perp}/u_{2\|}=ru_{2\perp}/u_{1\|}$, the tangent of the
angle the postshock flow makes to the shock normal with $r$ being the shock
compression ratio. Dividing (3) by (4) to eliminate $k_{1}$ and $k_{2}$ leads
to

${V_{A1}\over\sin\theta_{1}}-u_{1}\cot\theta_{1}={V_{A2}\over\sin\theta_{2W}}-u_{2}\cot\theta_{2W}\cos\alpha-
u_{2}\sin\alpha,$ (5)

which can be rearranged to give an equation for $\theta_{2W}$;

${M_{A2}\cos\theta_{2W}-1\over
M_{A2}\sin\theta_{2W}}=r{M_{A1}\cos\theta_{1}-1\over
M_{A1}\sin\theta_{1}}-\tan\alpha.$ (6)

For $M_{A1}\rightarrow\infty$, $M_{A2}\rightarrow\infty$, (and hence
$\tan\alpha\rightarrow 0$) and well away from the precisely perpendicular
shock where $\cos\theta_{2W}=\cos\theta_{1}=0$, this reduces to

$\cot\theta_{2W}=r\cot\theta_{1}={r^{2}\over\tan\theta_{2}}$ (7)

where we have also used the result $\tan\theta_{2}=r\tan\theta_{1}$ from
equation A4. Writing $\tan\theta_{2W}=\tan\left(\theta_{2}-\Delta\right)$
where $\Delta$ is the angle between the magnetosonic wave propagation
direction and the postshock magnetic field direction, we find

$\tan\Delta={\left(r^{2}-1\right)\over r}\cos\theta_{1}\sin\theta_{1}.$ (8)

The maximum deviation of the wave propagation from the magnetic field, and
hence the maximum wave compression is achieved for $\theta_{1}=45^{\circ}$,
where $\tan\Delta=-15/8$ (assuming $r=4$) and $\Delta=62^{\circ}$. Solving for
$\theta_{2}$ we find $\theta_{2}=76^{\circ}$ and $\theta_{2W}=14^{\circ}$.

For magnetosonic waves at high $\beta$ and for Alfvén waves equation 3 above
becomes

$\displaystyle
k_{1}V_{A1}-k_{1}\cos\theta_{1}u_{1}=k_{2}V_{A2}\cos\left(\theta_{2W}-\theta_{2}\right)$
$\displaystyle-k_{2}u_{2}\cos\left(\theta_{2W}-\alpha\right)$ (9)

accounting for the different dispersion properties of the Alfvén wave.
Following the same approach as above, we find

$\displaystyle{M_{A2}-\cos\theta_{2}\over
M_{A2}}\cot\theta_{2W}=r\cot\theta_{1}-{r\over M_{A1}\sin\theta_{1}}$
$\displaystyle-\tan\alpha+{\sin\theta 2\over M_{A2}},$ (10)

which can be recast as

$\displaystyle\left(1-{r^{1/2}\cos\theta_{1}\over
M_{A1}}\right)\cot\theta_{2W}=r\cot\theta_{1}-{r\over M_{A1}\sin\theta_{1}}$
$\displaystyle-\tan\alpha+{r^{3/2}\sin\theta_{1}\over M_{A1}}.$ (11)

In the limit $M_{A1}\rightarrow\infty$, ensuring $\beta>>1$ in the postshock
medium, this reduces to $\cot\theta_{2W}=r\cot\theta_{1}$ as above. Away from
this limit, $\cot\theta_{2W}$ will be different for magnetosonic and Alfvén
modes. In supernova remnants, this is unlikely to be a significant effect, but
in the lower $M_{A1}$ shocks driven by solar coronal mass ejections, this
could be an important consideration. The limit $M_{A1}\rightarrow\infty$ also
captures the case when $\omega_{1}\rightarrow 0$, appropriate for the Bell
nonresonant instability where modes grow at zero frequency.

##  3\. Wave Transmission and Reflection at Shocks

Vainio & Schlickeiser (1999) give transmission and reflection coefficients for
forward and backward propagating Alfvén waves at a parallel shock. Using the
jump conditions for the tangential electric field,

$\left[u_{n}B_{t}-B_{n}u_{t}\right]=0,$ (12)

the continuity of the transverse momentum,

$\left[\rho u_{n}u_{t}-B_{n}B_{t}/4\pi\right]=0,$ (13)

and the continuity of the mass flux

$\left[\rho u_{n}\right]=0,$ (14)

their result is

${T\atop R}={\left(M_{A1}+H\right)\left(\sqrt{r}\pm 1\right)\sqrt{r}\over
2\left(M_{A1}\pm H\sqrt{r}\right)}.$ (15)

Here $M_{A}$ is the Alfvén Mach number of the shock, and $H$ is the cross
helicity of the waves. $H=1$ for forward and $H=-1$ for backward propagating
waves. At a parallel shock, the two polarizations, Alfvén and magnetosonic,
behave the same. At an oblique shock, differences emerge. For the Alfvén mode,
which perturbs magnetic field and velocity in the direction perpendicular to
the plane in which the wave refraction occurs, results for $T$ and $R$ are
obtained from equation 15 with the replacement $M_{A1}\rightarrow
M_{A1}/\cos\theta_{1}$.

The magnetosonic mode is much more involved since perturbed magnetic field and
velocity vectors lie in the plane of refraction. It appears to be tractable
only in the limit $M_{A1}\rightarrow\infty$, when the motion of the shock
front itself in response to the passing wave can be neglected. We evaluate the
jump conditions for the tangential electric field and the continuity of
transverse momentum using

$\displaystyle B_{n}=-B_{2}\cos\theta_{2}+\delta B\sin\theta_{2W}$ (16)
$\displaystyle B_{t}=B_{2}\sin\theta_{2}+\delta B\cos\theta_{2W}$ (17)
$\displaystyle u_{n}=-u_{2\|}+\delta u\sin\theta_{2W}$ (18) $\displaystyle
u_{t}=u_{2\perp}+\delta u\cos\theta_{2W}.$ (19)

Substituting into the first jump condition gives

$\displaystyle\left(-u_{1}+\delta
u_{1}^{f}\sin\theta_{1}\right)\left(B_{1}\sin\theta_{1}+\delta
B_{1}^{f}\cos\theta_{1}\right)$
$\displaystyle-\left(-B_{1}\cos\theta_{1}+\delta
B_{1}^{f}\sin\theta_{1}\right)\delta u_{1}^{f}\cos\theta_{1}$ (20)
$\displaystyle=\left(-u_{2\|}+\delta
u_{2}\sin\theta_{2W}\right)\left(B_{2}\sin\theta_{2}+\delta
B_{s}\cos\theta_{2W}\right)$ $\displaystyle-\left(-B_{2}\cos\theta_{2}+\delta
B_{s}\sin\theta_{2W}\right)\left(u_{2\perp}+\delta
u_{s}\cos\theta_{2W}\right).$

Multiplying out and assuming
$-u_{1}B_{1}\sin\theta_{1}=-u_{2\|}B_{2}\sin\theta_{2}-u_{2\perp}B_{2}\cos\theta_{2}$
yields

$\displaystyle\delta u_{1}B_{1}-u_{1}\delta B_{1}\cos\theta_{1}=\delta
u_{2}B_{2}\cos\left(\theta_{2}-\theta_{2W}\right)$ $\displaystyle-
u_{2\|}\delta B_{2}\cos\theta_{2W}-u_{2\perp}\delta B_{2}\sin\theta_{2W}.$
(21)

We put $\delta u_{1}=\delta u_{1}^{f}$ or $\delta u_{1}^{b}$ to represent an
initially forward or backward propagating wave, and then put $\delta
u_{2}=\delta u_{2}^{f}+\delta u_{2}^{b}=\left(-\delta B_{2}^{f}+\delta
B_{2}^{b}\right)/\sqrt{4\pi\rho_{2}}$ and rearrange to find

$\displaystyle M_{A1}\cos\theta_{1}\pm 1$ $\displaystyle=$
$\displaystyle{T^{f}\atop
R^{b}}\bigg{\\{}{B_{2}\over\sqrt{r}B_{1}}\cos\left(\theta_{2}-\theta_{2W}\right)$
(22) $\displaystyle+$ $\displaystyle{M_{A1}\over
r}\cos\theta_{2W}+M_{A1}{u_{2\perp}\over u_{1}}\sin\theta_{2W}\bigg{\\}}$
$\displaystyle-$ $\displaystyle{R^{f}\atop
T^{b}}\bigg{\\{}{B_{2}\over\sqrt{r}B_{1}}\cos\left(\theta_{2}-\theta_{2W}\right)-{M_{A1}\over
r}\cos\theta_{2W}$ $\displaystyle-$ $\displaystyle M_{A1}{u_{2\perp}\over
u_{1}}\sin\theta_{2W}\bigg{\\}}.$

The +ve sign in the left hand side corresponds to $T_{f}=\delta
B_{2}^{f}/\delta B_{1}^{f}$ and $R^{f}=\delta B_{2}^{b}/\delta B_{1}^{f}$,
while the -ve sign corresponds to $R^{b}=\delta B_{2}^{f}/\delta B_{1}^{b}$
and $T^{b}=\delta B_{2}^{b}/\delta B_{1}^{b}$. A similar procedure for the
second jump condition gives

$\displaystyle\pm M_{A1}\cos\theta_{1}$ $\displaystyle+$ $\displaystyle\cos
2\theta_{1}={T^{f}\atop R^{b}}\bigg{\\{}{M_{A1}\over\sqrt{r}}\cos\theta_{2W}$
(23) $\displaystyle-$ $\displaystyle M_{A1}\sqrt{r}{u_{2\perp}\over
u_{1}}\sin\theta_{2W}+{B_{2}\over
B_{1}}\cos\left(\theta_{2}+\theta_{2W}\right)\bigg{\\}}$ $\displaystyle+$
$\displaystyle{R^{f}\atop
T^{b}}\bigg{\\{}-{M_{A1}\over\sqrt{r}}\cos\theta_{2W}+M_{A1}\sqrt{r}{u_{2\perp}\over
u_{1}}\sin\theta_{2W}$ $\displaystyle+$ $\displaystyle{B_{2}\over
B_{1}}\cos\left(\theta_{2}+\theta_{2W}\right)\bigg{\\}}.$

Taking $M_{A1}\rightarrow\infty$, so that also $u_{2\perp}\rightarrow 0$ in
equations 22 and 23 we find

$\displaystyle{T^{bf}\atop R^{bf}}$ $\displaystyle=$
$\displaystyle{\cos\theta_{1}\over 2\cos\theta_{2W}}\sqrt{r}\left(\sqrt{r}\pm
1\right)$ (24) $\displaystyle=$ $\displaystyle\left(\sqrt{r}\pm 1\over
2\sqrt{r}\right)\sqrt{\sin^{2}\theta_{1}+r^{2}\cos^{2}\theta_{1}},$

where we have substituted from $\cot\theta_{2W}=r\cot\theta_{1}$ in the final
step. This agrees with Vainio & Schlickeiser (1999) and equation 15 in the
appropriate limits ($\theta_{1}\rightarrow 0$ and $M_{A1}\rightarrow\infty$).

The forgoing has treated Alfvén and magnetosonic waves of nonzero frequency in
the upstream and downstream shock regions. As mentioned above, when being
driven by the cosmic ray current, these waves grow in the upstream region at
zero frequency, and hence $\delta u_{1}=\rightarrow 0$. In equations 22 and 23
this leads to the terms $\pm 1$ and $\pm M_{A1}\cos\theta_{1}$ on the left
hand sides of equation 22 and 23 respectively being dropped. We give the
corresponding expression for the transmission and reflection coefficients for
Alfvén (A) and magnetosonic (M) modes

$\displaystyle{T\atop R}$ $\displaystyle=$
$\displaystyle{M_{A1}r/\cos\theta_{1}\pm H\sqrt{r}\over
2\left(M_{A1}/\cos\theta_{1}\pm H\sqrt{r}\right)}\rightarrow{r\over
2}\quad{\rm(A)}$ (25) $\displaystyle=$ $\displaystyle{1\over
2}\sqrt{\sin^{2}\theta_{1}+r^{2}\cos^{2}\theta_{1}}\quad{\rm(M)},$ (26)

which are independent of the cross-helicity $H$ in the limit that
$M_{A1}\rightarrow\infty$.

##  4\. Wave Damping

Behind the shock, the newly obliquely propagating waves are subject to various
damping mechanisms. Both Alfvén and magnetosonic polarizations decay by
turbulent cascade, and the magnetosonic wave, now compressive due to its
oblique propagation also decays by transit time damping by the shocked quasi-
thermal ions.

Many formulations exist for the damping by turbulent cascade. We follow
Boldyrev (2005) who gives an expression to cover the cases of both the weak
turbulence (Goldreich & Sridhar, 1995) and strong turbulence (Iroshnikov,
1963; Kraichnan, 1965). In these limits the approximate damping rates for the
large scale turbulence are

$\displaystyle\gamma_{GS}$ $\displaystyle=$ $\displaystyle k_{2\perp}\delta
u_{2}=\omega\tan\left(\theta_{2}-\theta_{2W}\right)\delta u_{2}/V_{A2}$ (27)
$\displaystyle\gamma_{IK}$ $\displaystyle=$ $\displaystyle k_{2\perp}\delta
u_{2}^{2}/V_{A2}=\omega\tan\left(\theta_{2}-\theta_{2W}\right)\delta
u_{2}^{2}/V_{A2}^{2},$ (28)

where $\omega=k_{2\|}V_{A2}$ from equation 1. The strong turbulence expression
for $\gamma_{IK}$ is probably the most applicable. These simplest expressions
refer to balanced turbulence, with equal wave amplitudes propagating in each
direction. In our case, and decay rate of the transmitted waves will depend on
the intensity of counter propagating reflected waves, and vice versa. We will
assume cosmic rays streaming ahead of the shock generate wave travelling in
one direction only (away from the shock), and then we will take $\delta
u_{2}=\delta u_{1}R$, where $R$ is the reflection coefficient calculated above
in equation 24, in the damping rates.

Figure 2.— Schematic diagram of a plane parallel wave front encountering the
shock. Upstream the wave is a parallel propagating circularly polarized wave,
with no spatial variation in the magnetic field strength. Behind the shock the
magnetosonic polarization is damped quickly, leaving a spatial variation in
magnetic field strength due to the surviving Alfvén polarization, with maxima
every half wavelength along the wave propagation direction, $\lambda_{1}/2$,
as indicated. When projected onto the shock front the distance between
magnetic field maxima becomes $\lambda_{1}/2\sin\theta_{1}$, where
$\theta_{1}$ is the shock obliquity, the angle of the upstream unperturbed
magnetic field to the shock normal.

The magnetosonic polarization is also subject to transit time damping by
shocked but quasi-thermal ions. This damping rate in the small gyroradius
limit is given by (Melrose, 1986)

$\displaystyle\gamma_{TTD}=-\int\int{4\pi^{2}q^{2}V_{A2}^{2}\over\hbar\omega
c^{2}}\left(v_{\perp}^{2}k_{\perp}\over\Omega\right)^{2}\delta\left(\omega-
k_{\|}v_{\|}\right)$ $\displaystyle\times\hbar k_{\|}{\partial f\over\partial
p_{\|}}2\pi p_{\perp}dp_{\perp}dp_{\|},$ (29)

where $\omega$ and $\Omega$ are the wave angular frequency and thermal ion
gyrofrequency respectively, $V_{A2}$ is the postshock Alfvén speed, and $v$
and $p$ are the ion velocity and momentum respectively, with subscripts
$\perp$ and $\|$ indicating components perpendicular or parallel to the
ambient magnetic field. We represent the ion distribution function $f$ as a
“kappa” distribution;

$f={n\over\left(2\pi\right)^{3/2}m^{3}v_{Th2}^{3}}{\Gamma\left(\kappa\right)\over\Gamma\left(\kappa-3/2\right)\kappa^{3/2}}\left[1+{p^{2}\over
2\kappa^{2}v_{Th2}^{2}}\right]^{-\kappa},$ (30)

where $m$ is the ion mass, $v_{Th2}$ is its thermal velocity,
$\Gamma\left(\kappa\right)$ is the Gamma function of with argument $\kappa$,
and $\kappa$ is the index of the distribution. For $\kappa\rightarrow\infty$,
$f$ tends to a Maxwellian. For finite $\kappa$, the kappa distribution has
extended wings compared to a Maxwellian, and may be taken to represent
suprathermal ions in addition to the Maxwellian core. Normalizing the kappa
distribution over an infinite momentum range, we require $\kappa>3/2$ to keep
the number density $n$ finite, and $\kappa>5/2$ to keep the energy finite.
Smaller values of kappa require a high-momentum cutoff to keep particle number
and energy densities finite. Substituting equation 29 into equation 28 we get

$\displaystyle\gamma_{TTD}$ $\displaystyle=$
$\displaystyle{2\pi^{2}q^{2}V_{A2}^{2}\over
m^{7}v_{Th2}^{5}\Omega^{2}c^{2}}{k_{\perp}^{2}\over
k_{\|}}\sqrt{2\over\pi}{n\Gamma\left(\kappa\right)\over\Gamma\left(\kappa-3/2\right)\kappa^{3/2}}$
$\displaystyle\times\left(2\kappa m^{2}v_{Th2}^{2}\right)^{\kappa+1}$
$\displaystyle\times\int_{0}^{\infty}\left[2\kappa
m^{2}v_{Th2}^{2}+{m^{2}\omega^{2}\over
k_{\|}^{2}}+p_{\perp}^{2}\right]^{-\kappa-1}p_{\perp}^{5}dp_{\perp}.$

The integral is performed by substituting
$p_{\perp}=\sqrt{m^{2}\omega^{2}/k_{\|}^{2}+2\kappa m^{2}v_{Th2}^{2}}\tan\psi$
and evaluating the resulting integral,
$\int_{0}^{\pi/2}\sin^{5}\psi\cos^{2\kappa-5}\psi
d\psi=\int_{0}^{1}\left(1-\cos^{2}\psi\right)^{2}\cos^{2\kappa-5}\psi
d\left(\cos\psi\right)$ to find

$\displaystyle\gamma_{TTD}={2\pi^{2}q^{2}nV_{A2}^{2}\over
m\Omega^{2}c^{2}}{k_{\perp}^{2}v_{Th2}\over
k_{\|}}\sqrt{2\over\pi}\left[1+{V_{A2}^{2}\over 2\kappa
v_{Th2}^{2}}\right]^{2-\kappa}$
$\displaystyle\times{\Gamma\left(\kappa\right)\over\Gamma\left(\kappa-3/2\right)\kappa^{3/2}}{8\kappa^{2}\over\left(\kappa-1\right)\left(\kappa-2\right)}.$
(32)

As $\kappa\rightarrow\infty$, $\gamma_{TTD}\rightarrow
4\sqrt{2\pi}k_{\perp}^{2}v_{Th2}\exp\left({-1/\beta}\right)/k_{\|}$. For
application below, where $\beta>>1$ and $\omega\simeq k_{\|}V_{A2}$, we
rewrite as

$\gamma_{TTD}\simeq
4\sqrt{2\pi}\omega\tan^{2}\left(\theta_{2}-\theta_{2W}\right)v_{Th2}/V_{A2}.$
(33)

For a fast shock neglecting energy losses to cosmic rays,
$v_{Th2}\simeq\sqrt{3}v_{s}/4\sim 2000$ km s-1, where the shock velocity in
Tycho $v_{s}\simeq 5000$ km s-1.

##  5\. How Did Tycho Get Its Stripes?

Eriksen et al. (2011) observe “stripes” in synchrotron emission from Tycho’s
SNR. The stripes are separated by about 10”, which corresponds to $6\times
10^{17}$ cm at the presumed 4.0 kpc distance (Hayato et al., 2010), or
$5\times 10^{17}$ cm at 3.2 kpc (Slane et al., 2014). We argue that these
structures arise at oblique shocks, where upstream parallel circularly
polarized waves which are undamped become obliquely propagating damped waves
postshock. The key point is that the magnetosonic polarization damps faster
than the Alfvén polarization, and so a spatial variation in the magnetic field
strength will emerge associated with the surviving Alfvén mode. Along the
direction of propagation the absolute magnitude of the magnetic field strength
has a maximum every half wavelength ($\lambda_{1}/2$). When projected onto the
shock front, the maxima are separated by a distance of
$\lambda_{1}/2\sin\theta_{1}$, as shown schematically in Fig 2. Although less
damped than the magnetosonic polarization, the Alfvén modes still damp within
a wavelength or so postshock. These magnetic field variations are illuminated
by synchrotron radiation from cosmic ray electrons, the radiative cooling time
for which typically lies between the magnetosonic and Alfvén mode damping
times. This hierarchy ensures that the damped magnetosonic polarization does
not contribute significantly to the emission, while the Alfvén mode only
contributes close to the shock, while the electrons are still radiating in
X-rays.

Figure 3.— Variation with shock obliquity $\theta_{BN}$ of the intensity
contrast $\Delta I$ between stripe maxima and minima, for the cases of
Iroshnikov/Kraichnan (IK) or Goldreich-Sridhar (GS) turbulence.

In the case of Iroshnikov-Kraichnan turbulence, the difference in synchrotron
intensity between maxima and minima in the stripes may be written

$\displaystyle\Delta I$ $\displaystyle\propto$
$\displaystyle\left[\exp\left(-\gamma_{IK}t\right)^{s+1}-\exp\left(-\gamma_{IK}t-\gamma_{TTD}t\right)^{s+1}\right]$
(34) $\displaystyle\times\cos^{s+1}\theta_{1}$

where $s$ is the index of the synchrotron photon spectrum and $s+1$ is the
power law dependence of the magnetic field strength on the synchrotron
emission. The final factor of $\cos\theta_{1}$ gives the approximate shock
obliquity dependence of the preshock magnetic field amplification. Goldreich-
Sridhar turbulence simply requires $\gamma_{GS}$ in place of $\gamma_{IK}$.
This intensity contrast is maximized at
$t=\ln\left(1+\gamma_{TTD}/\gamma_{IK}\right)/\gamma_{TTD}/\left(s+1\right)$
with value

$\displaystyle\Delta I$ $\displaystyle\propto$
$\displaystyle\bigg{[}\exp\left(-\gamma_{IK}/\gamma_{TTD}\right)\ln\left(1+\gamma_{TTD}/\gamma_{IK}\right)$
(35)
$\displaystyle-\exp\left(-\gamma_{IK}/\gamma_{TTD}-1\right)\ln\left(1+\gamma_{TTD}/\gamma_{IK}\right)\bigg{]}$
$\displaystyle\times\cos^{s+1}\theta_{1}.$

Figure 3 shows the dependence of $\Delta I$ on the shock obliquity
$\theta_{1}$, for both Iroshnikov-Kraichnan (IK; strong) turbulence and
Goldreich-Sridhar (GS; weak) turbulence, calculated from equation 34, using
equations 26, 27, and 32. We take $\delta u_{2}=\delta u_{1}R$, and $\delta
u_{1}\simeq 180$ km s-1, corresponding to $R\simeq 2$ from equations 25, and
postshock magnetic field and density of 180 $\mu$G and 1.2 amu cm-3
respectively (Slane et al., 2014). In both cases the emission is relatively
strongly peaked over a restricted range of angles, possibly suggesting a
reason why the stripes are not ubiquitous over Tycho, but only seen in certain
“special” regions.

The separation at which the stripes preferably appear can be used to estimate
the maximum cosmic ray ion energy in Tycho. The Bell instability grows
structures at a parallel wavevector in the ambient magnetic field given by
(Bell, 2004, 2005; Laming et al., 2014)

$\displaystyle k_{\|}$ $\displaystyle=$ $\displaystyle 1.5\times 10^{-8}\eta
n_{i}\left(u_{1}\over 5000~{}{\rm km~{}s}^{-1}\right)^{3}\left(3\mu~{}{\rm
G}\over B\right)$ (36)
$\displaystyle\times{\gamma_{max}-3\gamma_{1}/4\over\gamma_{max}\gamma_{1}\left(\ln\gamma_{max}-1\right)}={2\pi\over\lambda_{1}}.$

Here $\eta$ is the fraction of shock energy going into cosmic rays, $n_{i}$ is
the ion number density in the preshock medium. The highest energy cosmic rays
have Lorentz factor $\gamma_{max}$, while the highest energy magnetized cosmic
rays have Lorentz factor $\gamma_{1}$. The cosmic ray current with particles
with $\gamma_{1}<\gamma<\gamma_{max}$ drives the Bell instability. We have
also set $\delta B/B=1$ in equation 3 of Laming et al. (2014) for
$k_{\|}r_{g}$, appropriate for the far upstream region of the cosmic ray
precursor, and divided through by the cosmic ray gyroradius.

We identify $\lambda_{1}/2\sin\theta_{1}=5\times 10^{17}$ cm at a distance to
Tycho of 3.2 kpc (Slane et al., 2014, model A), and taking $\theta_{1}$ from
Fig 3, we estimate $\gamma_{max}$. Letting
$\gamma_{1}\rightarrow\gamma_{max}$, and taking $\eta=0.26$, $n_{i}=0.3$cm-3,
we find $\gamma_{max}=6\times 10^{5}-1\times 10^{6}$ corresponding to
$\theta_{1}=10^{\circ}-15^{\circ}$. Although model dependent, this broadly
agrees with the modeling of Slane et al. (2014), who find the cosmic ray
energy spectrum breaking from something close to a $p^{-4}$ power law to a
steeper fall-off just above $\gamma\sim 10^{5}$ (their Fig. 4), and the
original estimate of Eriksen et al. (2011), who simply identified the spacing
of the stripes with a cosmic ray gyroradius. Projected onto the shock velocity
vector, the wave phase varies with distance
$\lambda_{1}/2\times\cos\theta_{1}\simeq 1\times 10^{17}-1.5\times 10^{17}$
cm. Thus with a shock velocity of 5000 km s-1, changes in the projected
position of the stripes should be visible in $2\times 10^{8}-3\times 10^{8}$
seconds, i.e. 6 - 10 years. Eriksen et al. (2011) report an intensity contrast
of a factor of 25 between stripe maxima and minima. Relating this to $B^{s+1}$
with $s\sim 2.11$ we find a variation of $B$ between maxima and minima of
about a factor 3. This is significantly lower than the likely upstream
magnetic field contrast of between a factor of 9 and 15, but probably
reasonable when we consider that both polarizations are amplified upstream by
this amount, and then the magnetosonic polarization is preferentially damped
by transit time damping. In the case that the electron cooling time becomes
longer than the Alfvén mode damping time, the variation of magnetic field
along the shock front will be less clear, but stripes separated by
$\lambda_{2}/2$ oriented normal to ${\bf k}_{2}$ extending further inside the
shock would result, analogous to ripples associated with sound waves.

##  6\. Conclusions

We suggest a model for the synchrotron “stripes” or ripples observed in
Tycho’s SNR (Eriksen et al., 2011) based only on simple ideas about Alfvén
propagation and dissipation at the forward shock, calculated as a function of
its obliquity. We find that the contrast is highest, and therefore that the
stripes would be most visible for a narrow range of shock obliquities close to
15∘ in the case of strong Iroshnikov-Kraichnan turbulence. In our model, the
emission comes from behind the shock front, once the magnetosonic component of
the originally circularly polarized wave has been transit time damped. This is
in contrast to previous ideas which have placed the structures giving rise to
the synchrotron stripes in the shock precursor. In the model advanced here,
some motion of the stripes over a period of one to a few years should be
visible, leading to possibilities for an observational test.

This work was supported by Basic Research Funds of the Chief of Naval
Research. I acknowledge enlightening conversations with Chee Ng, and I am
grateful to Una Hwang and an anonymous referee for reading and commenting on
the paper.

##  Appendix A MHD Shock Relations

Here we reproduce a few relationships relating to MHD shock, following mainly
Melrose (1986). The postshock perpendicular flow speed is

$u_{2\perp}=u_{1}{\left(r-1\right)\sin\theta_{1}\cos\theta_{1}\over
M_{A1}^{2}-r\cos^{2}\theta_{1}}\rightarrow 0\quad{\rm as}\quad
M_{A1}\rightarrow\infty.$ (A1)

The transverse magnetic field is

$B_{2}\sin\theta_{2}={r\left(M_{A1}^{2}-\cos^{2}\theta_{1}\right)\over
M_{A1}^{2}-r\cos^{2}\theta_{1}}B_{1}\sin\theta_{1}\rightarrow
rB_{1}\sin\theta_{1}\quad{\rm as}\quad M_{A1}\rightarrow\infty.$ (A2)

Writing

${B_{2}\sin\theta_{2}\over
B_{1}\sin\theta_{1}}={r\left(M_{A1}^{2}-\cos^{2}\theta_{1}\right)\over
M_{A1}^{2}-r\cos^{2}\theta_{1}}={\cos\theta_{1}\sin\theta_{2}\over\cos\theta_{2}\sin\theta_{1}}$
(A3)

we find

$\tan\theta_{2}=r\tan\theta_{1}+{\left(r-1\right)\sin\theta_{1}\cos\theta_{1}\over
M_{A1}^{2}-r\cos^{2}\theta_{1}}=r\tan\theta_{1}+{1\over r}\tan\alpha$ (A4)

to be compared with

$\cot\theta_{2W}=r\cot\theta_{1}-\tan\alpha$ (A5)

from equation 6, for the wave refraction.

## References

  * Achterberg & Blandford (1986) Achterberg, A., & Blandford, R. D. 1986, MNRAS, 218, 551
  * Bamert et al. (2004) Bamert, K., Kallenbach, R., Ness, N. F., et al. 2004, ApJ, 601, L99
  * Bell (2004) Bell, A. R. 2004, MNRAS, 353, 550
  * Bell (2005) Bell, A. R. 2005, MNRAS, 358, 181
  * Boldyrev (2005) Boldyrev, S. 2005, ApJ, 626, L37
  * Bykov et al. (2011a) Bykov, A. M., Osipov, S. M., & Ellison, D. C. 2011a, MNRAS, 410, 39
  * Bykov et al. (2011b) Bykov, A. M., Ellison, D. C., Osipov, S. M., Pavlov, G. G., & Uvarov, Y. A. 2011b, ApJ, 735, L40
  * Caprioli & Spitkovsky (2013) Caprioli, D., & Spitkovsky, A. 2013, ApJ, 765, L20
  * Drury & Falle (1986) Drury, L. O., & Falle, S. A. E. G. 1986, MNRAS, 223, 353
  * Eriksen et al. (2011) Eriksen, K., Hughes, J. P., Badenes, C., et al. 2011, ApJ, 728, L28
  * Gargaté & Spitkovsky (2012) Gargaté, L., & Spitkovsky, A. 2012, ApJ, 744, 67
  * Goldreich & Sridhar (1995) Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763
  * Hayato et al. (2010) Hayato, A. et al. 2010, ApJ, 725, 894
  * Iroshnikov (1963) Iroshnikov, P. 1963, Astron. Zh., 40, 742
  * Kraichnan (1965) Kraichnan, R. H. 1965,Phys. Fluids, 8, 1385
  * Laming et al. (2014) Laming, J. M., Hwang, U., Ghavamian, P., & Rakowski, C. E. 2014, ApJ, 790, 11
  * Laming et al. (2013) Laming, J. M., Moses J. D., Ko, Y.-K., Murphy, R. J., Ng, C. K., Rakowski, C. E., & Tylka, A. J. 2013, ApJ, 770, 73
  * Malkov et al. (2012) Malkov, M. A., Sagdeev, R. Z., & Diamond, P. H. 2012, ApJ, 748, L32
  * McKenzie & Westphal (1970) McKenzie, J. F., & Westphal, K. O. 1970, Phys. Fluids, 13, 630
  * Melrose (1986) Melrose, D. B. 1986, Instabilities in Space and Laboratory Plasmas (Cambridge: Cambridge University Press)
  * Rakowski et al. (2011) Rakowski, C. E., Laming, J. M., Hwang, U., Eriksen, K. A., Ghavamian, P., & Hughes, J. P. 2011, ApJ, 735, L21
  * Slane et al. (2014) Slane, P. O., Lee, S.-H., Ellison, D. C., Patnaude, D. J., Hughes, J. P., Eriksen, K. A., Castro, D., & Nagataki, S. 2014, ApJ, 783, 33
  * Vainio & Schlickeiser (1999) Vainio, R. & Schlickeiser, R. 1998, A&A, 331, 793
  * Vink & Laming (2003) Vink, J., & Laming, J. M. 2003, ApJ, 584, 758
  * Vladimirov et al. (2009) Vladimirov, A. E., Bykov, A. M., & Ellison, D. C. 2009, ApJ, 703, L29
  * Webb et al. (1999) Webb, G. M., Zakharian, A., Brio, M., & Zank, G. P. 1999, J. Plasma Physics, 61, 295
  * Zank et al. (2006) Zank, G. P., Li, G., Florinski, V., Hu, Q., Lario, D., & Smith, C. W. 2006, J. Geophys. Res. Space Phys., 111, 6108


# Detecting Braess Routes: an Algorithm Accounting for Queuing Delays With an
Extended Graph
††thanks: This project was sponsored in part by the joint NSF-DOT CPS grant
CNS-1545116 - Traffic Operating System for Smart Cities.

Mikhail Burov  UC Berkeley
<EMAIL_ADDRESS>Can Kizilkale  UC Berkeley and LBNL
<EMAIL_ADDRESS>Alexander Kurzhanskiy  UC Berkeley
<EMAIL_ADDRESS>Murat Arcak  UC Berkeley
<EMAIL_ADDRESS>
###### Abstract

The Braess paradox is a counter-intuitive phenomenon whereby adding roads to a
network results in higher travel time at equilibrium. In this paper we present
an algorithm to detect the occurrence of this paradox in real-world networks
with the help of an improved graph representation accounting for queues. The
addition of queues to the network representation enables a closer match with
real data. Moreover, we search for routes causing this phenomenon (‘Braess
routes’) rather than links, and advocate removing such routes virtually from
navigation systems so that the associated links can continue to serve other
routes. Our algorithm relies on a convex optimization problem utilizing
Beckmann potentials for road links as well as queues, and results in a route
reconfiguration with reduced delay. We assume the availability of historical
data to build the optimization model. We also assume the existence of a
centralized navigation system to manage the routing options and remove the
Braess routes. The theoretical solution demonstrates up to 12% delay reduction
in a network from Montgomery County, Maryland. We validate the improvement
with simulations.

###### Index Terms:

Braess paradox, Optimization, Beckmann potential, BPR functions

##  I Introduction

Road networks suffer from various types of equilibrium inefficiency due to
selfish routing. Of particular interest is the Braess paradox, a counter-
intuitive phenomenon describing scenarios in which building new road links
results in higher traffic delays at equilibrium. Since its introduction in
1968 in [1], the Braess paradox has been studied extensively to find efficient
ways to predict, detect and prevent its existence.

Early results, such as [2], [3], [4], [5], [6], focus on the classic diamond-
shaped, four-node network with a single OD-pair. A later study [7] extends the
analysis to general traffic networks to predict the occurrence of the paradox;
however, the applicability of the results is limited by a restrictive
assumption that all routes with non-zero flows in the original network are
also utilized after the addition of a road link. Moreover, [7] as well as
related theoretical papers [8] and [9] consider the special case when exactly
one road link or route is built or removed from the network.

Another approach to anticipating the Braess paradox is a network topology
analysis. The study [10] shows that a two-terminal network is immune to Braess
paradox if and only if it is series-parallel. In addition, the paper extends
the result to account for any Pareto inefficiency in a two-terminal network.
References [11] and [12] extend the characterization from undirected graphs to
directed graphs and allow for multiple commodities. Another theoretical study,
[13], explores the concept of matroids to identify networks that are immune to
the Braess paradox. Unfortunately, few transportation networks exhibit a
matroid structure; therefore, the applicability is limited in practice. A
further shortcoming of the theoretical studies mentioned above is that they
present structures that are immune to the Braess paradox, but do not provide
tools to modify existing networks to eliminate this paradox and improve the
efficiency of the equilibrium.

The computational study [14] proposes a mathematical programming method to
detect the Braess paradox in a given network. The problem is formulated as a
bi-level structure and then transformed into a single-level mixed integer
program. By setting tolls to links and analysing the resulting network
latency, the algorithm detects links that cause the Braess paradox and
penalizes them to imitate road closure.

Instead of links, in our paper we search for ‘Braess routes’ (routes that
cause the Braess paradox) with a greedy optimization algorithm and remove them
sequentially from the navigation system. This leaves the drivers with a subset
of routes that are immune to the paradox, which they are free to choose from.
We believe that removing routes is advantageous over removing links, as
removing a link adversely affects all other routes going through this link.
Among other scenarios, our approach enables removing through-traffic from
residential roads, while allowing them to continue to serve the residents. One
shortcoming of the route removal approach, however, is that customers might
simply switch to another navigation system if their freedom of choice is
limited. Further research is needed to address this issue, such as splitting
populations into “selfish” and “altruistic” to model different levels of
cooperation or creating a system of benefits to encourage drivers’
participation.

Unlike other methods for detecting and eliminating the Braess paradox, our
model accounts for intersection structures and queues. A more accurate graph
representation achieved by fitting functions for link delays and queue delays
from data allows us to make theoretical methods applicable to real-world
networks. We validate our approach on an extended graph model of a road
network from North Bethesda, Montgomery County, Maryland. We slightly modified
the geometry of this network with additional links to be able to test a bigger
number of OD-pairs and routes. We were able to demonstrate up to 12% delay
reduction on this extended network. Despite this improvement, we do not claim
the proposed route removal strategy is optimal. Indeed, as shown in [15], the
problem of finding the optimal subnetwork is NP-hard. Instead we trade
optimality with computational tractability and applicability to real-world
networks.

##  II Network Specification

We consider a network with several routes available to each OD-pair. Every
route is represented by a set of consecutive links connecting the origin to
the destination.

###  II-A Link delay function

We estimate the time delay that a vehicle experiences on each link with a
Bureau of Public Roads (BPR) function [16]:

$\Phi(z)=t_{0}\left(1+a\left(\frac{z}{cap}\right)^{b}\right),$ (1)

where

$t_{0}$ is the free-flow time ($t_{0}=\frac{link-length}{free-flow-speed}$),

$a$ and $b$ are parameters that depend on the link’s properties,

$cap$ is the link throughput capacity (in $\frac{veh}{hour}$).

To compute the values of parameters $a$ and $b$ for a particular link, we
followed the method described in [17], which requires speed-density data
points for a curve-fitting algorithm that outputs appropriate parameters. We
used Simulation of Urban Mobility (SUMO) open-source software to build a test
case based on the modified North Bethesda, Montgomery County, Maryland network
around the intersections of Montrose Rd and Montrose Pkwy (Fig. 1). In the
absence of historical data, we generated traffic data from the simulation
environment and estimated the delay functions. The derivation is detailed in
Appendix A.

Figure 1: Montgomery County network, slightly modified with additional links
to test a bigger number of OD-pairs and routes. The red box shows how one edge
is represented with a series of links, as discussed in Appendix A.

###  II-B Queue delay estimation

Link delay functions do not account for the queues, which accumulate when the
link throughput capacity is insufficient for the incoming flow. To address
this issue and model queue delay we introduce “phantom” links, i.e. links that
have no physical analogue in the real-world or simulation, but exist solely to
account for additional delay related to queues. We insert “phantom” links into
routes between consecutive edges incoming to intersections and edges leaving
the same intersections (Fig. 2).

Figure 2: “Phantom” link insertion.

The analysis in Appendix B yields queue delay function:

$\Phi_{q}(z)=\begin{cases}d_{0}&\text{if $z<s$}\\\ \alpha
z+\beta&\text{otherwise,}\end{cases}$ (2)

where

$d_{0}$ is expected constant delay due to intersection structure,

$\alpha$ and $\beta$ are parameters that depend on the link’s properties,

$s$ is the link’s saturation rate.

The proposed function is a continuous non-decreasing non-negative piece-wise
linear function. Therefore, it can be readily used in our route detection
algorithm.

##  III Algorithm

###  III-A Wardrop Equilibrium computation

To compute the Wardrop Equilibrium state we use a convex optimization problem
that utilizes Beckmann potentials [18] for delay functions. The objective
function of the optimization problem to minimize is the sum of Beckmann
potentials, i.e. integrals of delay functions, across all links:

$\displaystyle\sum_{i\in L}\phi_{i}(z_{i})$ (3)

where

$L$ is the set of links in the network,

$\phi_{i}$ is the Beckmann potential of the link $i$
($\phi_{i}^{{}^{\prime}}=\Phi_{i}$),

$\Phi_{i}$ is the delay function of the link $i$,

$z_{i}$ is the flow on the link $i$ ($z_{i}=(R^{T}x)_{i}$),

$x$ is the vector of route flows,

$R$ is the routing matrix defined as

$R_{ij}=\begin{cases}1,\text{ if the route $i$ goes through the link $j$}\\\
0,\text{ otherwise.}\end{cases}$ (4)

The first constraint ensures non-negative flows on routes:

$x^{k}_{j}\geq 0,\ \forall j\in P^{k};\forall k\in O,$ (5)

where

$P^{k}$ is the set of routes corresponding to the OD-pair $k$,

$O$ is the set of OD-pairs.

The second constraint guarantees that flows on routes corresponding to the
same OD-pair sum up to the demand for that pair:

$\displaystyle\sum_{j\in P^{k}}x^{k}_{j}=d_{k},\ \forall k\in O,$ (6)

where $d_{k}$ is the demand on the OD-pair $k$.

Therefore, the problem is:

$\underset{x}{\text{min}}\ (\ref{obj})\ \text{subject to}\
(\ref{non_neg})-(\ref{eq_dem}).$ (7)

The solution of the problem (7) is the vector $x^{*}$ of route flows at
equilibrium, which results in the total network delay:

$Y=x^{*T}R\begin{bmatrix}\Phi_{1}(z_{1}^{*})\\\ ...\\\
\Phi_{m}(z_{m}^{*})\end{bmatrix}_{|z_{i}^{*}=(R^{T}x^{*})_{i}}.$ (8)

###  III-B Elimination procedure

The objective of our algorithm is to find subsets of links/routes such that
the remaining network has reduced latency. To achieve this we remove links and
routes that cause the Braess paradox, which we refer to as ‘Braess’
links/routes, from the original system. Removing a route implies removing it
from the set of options suggested by the navigation system (e.g., Google
Maps); removing a link implies either physically closing down the road segment
or reducing its capacity with tolls or signaling. As explained in the
Introduction, route removal is preferable in practice; however, for the
completeness of the study we discuss several approaches for route and link
elimination, as they are easily obtained from the equilibrium computation
method of the previous section. An important constraint in these approaches is
to keep the connectivity of the network unchanged, i.e. every OD-pair from the
original network must retain at least one route in the reduced network.

In each approach we attribute values to the links/routes with the help of the
optimization problem (7). First we obtain the equilibrium network delay, $Y$,
for the original system. Then, we tentatively remove a link/route from the
network (remove corresponding rows and columns from the routing matrix (4))
and solve (7) again, deriving the new equilibrium network delay, $Y_{new}$.
Having both the original and the new network delays, the value $V$ of the
removed link/route is:

$V=Y_{new}-Y.$ (9)

Links/routes with negative values are associated with the Braess paradox,
because removing them reduces latency.

Greedy Single Link Removal

  1. 1.

Find the link with the smallest value $V_{min}$.

  2. 2.

If $V_{min}<0$, remove the link and revert to step (1).

  3. 3.

Otherwise, terminate.

As stated earlier, removing a link affects the entire set of routes going
through this link, which compromises the effectiveness of the algorithm.
Moreover, this approach is slower than some algorithms discussed further.

Link Combination Removal

  1. 1.

Compute the values of all possible combinations of links and find the one with
the smallest value $V_{min}$.

  2. 2.

If $V_{min}<0$, remove the corresponding combination.

  3. 3.

Otherwise, the network is Braess paradox-free.

Since the number of link combinations grows exponentially with the number of
links, this method is computationally expensive and does not scale well to
large networks.

Link-Route Combination Removal This approach modifies the first method and
tries to address the issue with subset of routes elimination. Instead of
removing links completely,

  1. 1.

For each link find the subset of routes utilizing this link that results in
the minimal network delay.

  2. 2.

Remove the routes that were not present in at least one optimal configuration.

Unlike previous algorithms, this method removes routes, but does it indirectly
by working with link-route combinations.

Greedy Single Route Removal

  1. 1.

Find the route with the smallest value $V_{min}$.

  2. 2.

If $V_{min}<0$, remove the route and revert to step (1).

  3. 3.

Otherwise, terminate.

This approach, unlike the first method, deals with routes directly, which is
faster since in practice the number of commonly used routes is smaller than
the number of links. An OD-pair can potentially have a large number of
possible routes; however the dominant part of the drivers uses the a very
limited subset of those routes. It also allows to identify the occurrence of
the Braess paradox after the first iteration: if the route with the minimal
value has zero initial flow at equilibrium, then the network is Braess
paradox-free.

Route Combination Removal

  1. 1.

Compute the values of all possible combinations of routes and find the one
with the smallest value $V_{min}$.

  2. 2.

If $V_{min}<0$, remove the corresponding combination.

  3. 3.

Otherwise, the network is paradox-free.

This approach is similar to the second approach and has the same limitation
due to computational complexity.

All five approaches were tested to identify the fastest and most accurate
method for further implementation and verification. The improvements achieved
on the sample network were comparable. Therefore, we report below the results
for the fourth method (Greedy Single Route Removal) due to its computational
speed and advantages in implementation.

##  IV Simulation

Series of simulations were conducted to demonstrate the improvement, evaluate
the prediction accuracy and estimate the computational precision of our
algorithm when applied to physical networks. The testing was designed to model
the behavior of a real-world traffic system, therefore a complex structure of
OD-pairs and corresponding routes was built upon (Fig. 1). The parameters of
the simulation model are presented in Table I.

Table I: Simulation parameters. Parameters |  Values
---|---
Number of OD-pairs |  12
Number of routes per OD-pair |  1-5
Total number of routes |  38
Number of simulated vehicles |  500 - 3200

A side benefit of the route removal approach is that it may help small side-
roads in residential areas that usually suffer from congestion due to drivers
attempting to avoid busy freeways. In the era of Google maps, Waze and other
navigation systems, which try to discover short-cuts and guide vehicles
through neighborhoods to free up main roads, residents of these regions face
busy traffic consequences and spend significantly more time than usual on
short trips [19].

The testing procedure consists of the following steps:

  1. 1.

Solve the problem (7) for the original network to obtain the equilibrium flow
distribution, $x^{*}$, and equilibrium network delay, $Y$.

  2. 2.

Derive the route configuration resulting in minimal delay via route
elimination algorithm and reduce the network.

  3. 3.

Solve the problem (7) for the reduced network to obtain the new equilibrium
flow, $x^{*}_{new}$, and network delay, $Y_{new}$. The theoretical delay
reduction is: $I_{th}=\frac{Y-Y_{new}}{Y}$.

  4. 4.

Feed $x^{*}$ into the SUMO to obtain the simulation network delay, $Y^{sim}$,
as a sum of travel times of all vehicles.

  5. 5.

Feed $x^{*}_{new}$ into the SUMO to obtain the new simulated network delay,
$Y^{sim}_{new}$. The simulation delay reduction is
$I_{sim}=\frac{Y^{sim}-Y^{sim}_{new}}{Y^{sim}}$.

  6. 6.

Compare the theoretical and simulated delay reductions $I_{th}$ and $I_{sim}$,
and estimated total delays $Y$ and $Y^{sim}$.

Every simulation set was run with different initial demands corresponding to
OD-pairs (last row of Table I). Congestion scenarios were of particular
interest, because they allowed us to test the developed queue delay estimation
model.

It is important to mention that our model does not account for spillbacks
(full occupation of a link that causes the queue propagation to the upstream
edge). Therefore, upper bounds on demands were applied to avoid spillbacks.
Further research is required for an appropriate spillback representation.

##  V Results

In this section we present the results of our algorithm for several scenarios.
The first two rows in Table II show the size of the simulated traffic. We
addressed conditions associated with different times of day.

Table II: Simulation results. |  Set-up 1 |  Set-up 2 |  Set-up 3 |  Set-up 4
---|---|---|---|---
Demand |  Low |  Medium |  High |  High
Number of vehicles |  500 |  1600 |  2600 |  3300
Improvement ($I_{th}$) |  0% |  3.2% |  11.8% |  8.1%
Improvement Diff. |  0% |  5.3% |  15.3% |  10.7%
Network Delay Diff. |  1.7% |  8.1% |  9.5% |  9.2%

We are interested in three main metrics when analyzing the efficiency of our
approach, which are presented in Table II. The first metric (Row 3) is the
theoretical improvement in the network delay, i.e. total travel time saved by
all vehicles after implementing our algorithm. In the low demand scenario, the
network was Braess paradox-free, and no travel time reduction was achieved.
For moderate and congested traffic the improvement ranges from 3.2% to 11.8%,
which is reasonable considering the insignificant effort required to modify
the network configuration.

The second metric (Row 4) shows the difference between the theoretically
predicted delay improvement and the corresponding simulation result. For the
free-flow set-up there was no Braess paradox and no comparison was needed. For
medium and high demands, the difference varies between 5.3% and 15.3%. Taking
into account the improvement value itself, we can conclude that the algorithm
makes a relatively accurate prediction. Additionally, the theoretical solution
always resulted in actual simulation improvement, i.e., $I_{th}>0\Rightarrow
I_{sim}>0$.

The third metric (Row 5) is the difference between the predicted network delay
and the simulated one. According to the results, we managed to keep the
deviation under 10% for all scenarios. Furthermore, some cases demonstrated
almost identical values for the network delays. The accuracy of the prediction
suggests that the proposed graph representation accounting for queue delays
yields a reasonable model to reflect real traffic conditions.

##  VI Conclusion

We presented a greedy optimization algorithm to detect and sequentially remove
the routes that contribute to the Braess paradox. The algorithm assigns values
to existing routes in terms of how much the network delay would increase if
they were removed, so that routes with negative values indicate a Braess
paradox. The algorithm removes the route with the most negative value at each
iteration. The resulting reduced network configuration demonstrates up to 12%
travel time reduction. The proposed extended graph representation to account
for queue delays provides a major improvement in how well the theoretical
predictions match the simulated ones. We are interested in further studying
the possible methods of modeling spillbacks and other congestion-related
phenomena.

##  Appendix A Link Delay Function Derivation

In this section we present a detailed link delay function derivation from
empirical data points. We treated flow data for each link individually,
simulating traffic on one edge at a time. To capture a wide range of flow
values, we generated a new random flow from the interval [0,3000]
$\frac{veh}{h}$ every 200 seconds. The provided data were sufficient to
construct BPR-functions for 90% of links (Fig. 3). Tuning simulation
parameters covered additional 70% of the remaining links (Fig. 4). For the
rest of the failed edges, we used parameter values corresponding to delay
functions of upstream successful links with minor modifications. We assumed
these parameters are likely to have close values and, therefore, can be almost
interchangeable. However, if this assumption is false the discrepancy between
ground truth delay and estimated delay produced by inaccurate parameter
substitution of one link is insignificant in a network scale.

(a) BPR-function fitting.

(b) Capacity estimation.

Figure 3: Sufficient data examples.

Figure 4: Insufficient data examples.

The negligibility of this difference follows from the fact that road links
constituting a network are relatively short in general. One edge in the graph
network representation is usually displayed by a series of short connected
links in the SUMO simulation network (Fig. 1). The travel time contribution of
one link is in the order of few seconds, which is insignificantly small to
deviate from the ground truth. Moreover, since the number of poorly-fitted
links makes up less than 2% of the number of all links, we can conclude that
the proposed design is sufficiently accurate.

Additionally, the link capacities can be extracted from fundamental diagrams
built upon the collected data (in $\frac{veh}{h}$) (Fig. 3(b)). The first
method is to set the capacity estimate to the maximal recorded flow value.
Another approach is to use a piece-wise-linear approximation to derive the
capacity as a y-coordinate of the intersection of two approximation lines.

##  Appendix B Queue Delay Function Derivation

In this section we present a detailed queue delay function derivation based on
the queue formation analysis. Queues occur when the flow ($z$) on the link
exceeds the saturation rate ($s$) of this link. The saturation rate is the
upper bound on the number of vehicles able to leave the link within a period
of time. Therefore, the queue formation depends on the difference between the
inflow and the outflow on a particular link. Depending on the link type and
its relative position with a specific intersection, we can distinguish several
possible options for saturation rate estimation:

  * •

If the link incomes to a signalized intersection, its corresponding saturation
rate equals to the maximum number of vehicles ($n$) able to pass the
intersection on green light within one cycle normalized to one hour:

$s=\frac{3600n}{D}$ (10)

where, $D$ is the cycle duration (in seconds).

  * •

If there is a STOP sign at the end of the link, its corresponding saturation
rate is:

$s=\frac{3600}{w}$ (11)

where, $w$ is the delay (in seconds) a vehicle experiences when forced to stop
at the STOP sign on an empty road. This value depends primarily on the speed
limit on the link and has a small fluctuation from edge to edge.

  * •

If the link is free from any traffic regulation causing vehicles to stop or
slow down, the saturation rate equals to the link physical throughput
capacity.

  * •

The link is the secondary link, having no priority on an unsignalized
intersection. Up to this point, the queue delay is given as $\Phi_{q}(z)$
which is a function of the flow on the corresponding link. In this case,
however, the delay incurred by the intersection structure on the secondary
link is also a function of the flow of the primary link. Although this
violates the fundamental assumption of delay being a function of the
corresponding link flow only, if the delay was symmetric for both the primary
and secondary links, we could have still implemented this in our optimization
problem and computed the equilibrium points. However, primary link sees no
delay from the intersection hence the delay function is asymmetric which makes
it hard if not impossible to compute equilibrium points the same way.
Therefore, this scenario is not represented in our model and is subject to
further research.

Similar to the link delay function derivation, we need to use real-world or
simulated empirical data to estimate the parameters for queue delay function.
We present the detailed function computation procedure for only the first
scenario, which features a link incoming to a signalized intersection. Other
scenarios utilize a slightly modified approach, which we will mention at the
end of this section.

Knowing the traffic light cycle length ($D$) and keeping in mind the queue
size dependency on the inflow-outflow difference, we estimate the one cycle
queue growth rate as:

$dq=\frac{(z-N)D}{3600},$ (12)

where $N$ is the throughput of the link.

Based on the flow value, we can distinguish two scenarios:

  * •

The flow is smaller than the saturation rate ($z<s$). In this case, $dq=0$,
all vehicles are able to pass the intersection and no queue is forming. To
account for a potential red phase arrival, instead of setting the queue delay
to zero, we choose a specific constant $d_{0}$, which is the expected value of
a delay due to phase change:

$d_{0}=\operatorname{\mathbb{E}}(y)=\frac{1}{D}\displaystyle\sum_{i=1}^{L_{red}}i=\frac{L_{red}(1+L_{red})}{2D},$
(13)

where $L_{red}$ is the duration of the red phase.

  * •

The flow is greater than the saturation rate ($z>s$). In this case, $dq>0$ and
the queue is growing. To find the growth rate, the flow-throughput dependency
and, thus, the saturation rate for a particular link are required.

Based on the queue formation model presented earlier, queue delay function
depends on the saturation rate of the link. To find the one-cycle throughput
capacity of the link we simulate flows of various values in ascending order
for a short period of time (200 seconds) each and record the number of
vehicles leaving the link (entering the intersection) within one cycle. The
flow point at which the throughput linear growth stops (Fig. 5(a)) corresponds
to the one-cycle saturation rate of the link. To obtain the value for one-hour
period, scale the result by $\frac{3600}{D}$. The throughput of the link (in
$\frac{veh}{h}$) is a piece-wise function of the flow which has the form:

$N=\begin{cases}z&\text{if $z<s$}\\\
\frac{2(z-s)}{z}+s&\text{otherwise.}\end{cases}$ (14)

For the flows smaller than the saturation rate, the throughput is a simple
linear function, because all vehicles are able to leave the link. Otherwise,
the throughput equals to the sum of the saturation rate and an additional
term, which never exceeds 2. This term accounts for rare scenarios with one or
multiple abnormally fast-moving vehicles that manage to exceed the usual
saturation rate.

(a) Saturation rate.

(b) Throughput function.

Figure 5: Link throughput at a signalized intersection.

Feeding the derived throughput function back into equation (12), we obtain a
one-cycle gain in numbers of vehicles to the queue due to excessive flow. It
is important to note that the imposed queue delay varies among vehicles within
the same flow, and drivers in the head of a heavy traffic would spend
significantly less time in the queue, than the ones in the tail. To avoid
complications in the optimization problem formulation, an average delay can be
assigned to every vehicle on the link. Plotting the computed delay results in
a piece-wise linear function (Fig. 6(a)).

(a) Queue delay data.

(b) Fitting delay function.

Figure 6: Signalised intersection flow-delay dependency.

The last step is to fit a linear function $\alpha z+\beta$ to the non-constant
area of the graph (Fig. 6(b)) and learn the parameters $\alpha$ and $\beta$
resulting in equation (2).

The presented approach can also be used to determine queue delay functions for
links ending up with the STOP sign. Both the expected delay $d_{0}$ and the
cycle length $D$ would be artificially set to $w$ from equation (11), since
every vehicle without exception is required to stop at the STOP sign, which
takes it exactly $w$ seconds to do.

## References

  * [1] Dietrich Braess.  Über ein paradoxon aus der verkehrsplanung.  Unternehmensforschung, 12(1):258–268, 1968.
  * [2] John D Murchland.  Braess’s paradox of traffic flow.  Transportation Research, 4(4):391–394, 1970.
  * [3] Larry J Leblanc.  An algorithm for the discrete network design problem.  Transportation Science, 9(3):183–199, 1975.
  * [4] Caroline Fisk.  More paradoxes in the equilibrium assignment problem.  Transportation Research Part B: Methodological, 13(4):305–309, 1979\.
  * [5] NF Stewart.  Equilibrium vs system-optimal flow: some examples.  Transportation Research Part A: General, 14(2):81–84, 1980.
  * [6] Marguerite Frank.  The braess paradox.  Mathematical Programming, 20(1):283–302, 1981.
  * [7] Richard Steinberg and Willard I Zangwill.  The prevalence of Braess’ paradox.  Transportation Science, 17(3):301–318, 1983.
  * [8] Stella Dafermos and Anna Nagurney.  On some traffic equilibrium theory paradoxes.  Transportation Research Part B: Methodological, 18(2):101–110, 1984\.
  * [9] Azuma Taguchi.  Braess’paradox in a two-terminal transportation network.  Journal of the Operations Research Society of Japan, 25(4):376–389, 1982.
  * [10] Igal Milchtaich.  Network topology and the efficiency of equilibrium.  Games and Economic Behavior, 57(2):321–346, 2006.
  * [11] Pietro Cenciarelli, Daniele Gorla, and Ivano Salvo.  Graph theoretic investigations on inefficiencies in network models.  arXiv preprint arXiv:1603.01983, 2016.
  * [12] Xujin Chen, Zhuo Diao, and Xiaodong Hu.  Excluding braess’s paradox in nonatomic selfish routing.  In International Symposium on Algorithmic Game Theory, pages 219–230. Springer, 2015.
  * [13] Satoru Fujishige, Michel X Goemans, Tobias Harks, Britta Peis, and Rico Zenklusen.  Matroids are immune to Braess’ paradox.  Mathematics of Operations Research, 42(3):745–761, 2017.
  * [14] Yangbeibei Ji, Wei Mao, and Xiaoning Zhang.  Detecting braess paradox links with a mixed integer linear programme.  International Journal of Industrial and Systems Engineering, 17(3):275–284, 2014.
  * [15] Tim Roughgarden.  On the severity of Braess’s paradox: designing networks for selfish users is hard.  Journal of Computer and System Sciences, 72(5):922–953, 2006.
  * [16] United States. Bureau of Public Roads.  Traffic assignment manual for application with a large, high speed computer, volume 2.  US Department of Commerce, Bureau of Public Roads, Office of Planning, Urban …, 1964.
  * [17] Rafał Kucharski and Arkadiusz Drabicki.  Estimating macroscopic volume delay functions with the traffic density derived from measured speeds and flows.  Journal of Advanced Transportation, 2017, 2017.
  * [18] Martin Beckmann, Charles B McGuire, and Christopher B Winsten.  Studies in the economics of transportation.  Technical report, 1956.
  * [19] Jane Macfarlane.  When apps rule the road: The proliferation of navigation apps is causing traffic chaos. it’s time to restore order.  IEEE Spectrum, 56(10):22–27, 2019.