# Non-Neutralized Electric Current of Active Regions Explained as a Projection
Effect

Xudong Sun  Mark C. M. Cheung  Institute for Astronomy, University of Hawai‘i
at Mānoa, Pukalani, HI 96768, USA  Lockheed Martin Solar and Astrophysics
Laboratory, Palo Alto, CA 94304, USA

###### keywords:

Active Regions, Magnetic Fields; Electric Currents and Current Sheets;
Magnetic fields, Models; Magnetic fields, Photosphere

##  1 Introduction

sec:intro

###  1.1 Definitions of Current Neutralization

subsec:definition

Magnetic field in the solar interior is believed to exist in a “fibril” state:
isolated flux tubes generated from dynamo actions are embedded in a relatively
field-free background (Fan, 2009). Strong field strength and twist keep the
tubes coherent against vigorous convection (Schüssler, 1979). This places an
interesting constraint on the net electric current $I$. Using Ampère’s law and
Stokes theorem, we have

$\ilabel{eq:stokes}\begin{split}{\bm{j}}&=\mathbf{\nabla}\times{\bm{B}},\\\
I&=\int_{C}{\bm{j}}\cdot\hat{\bm{n}}\,\mathrm{d}A=\oint_{\partial
C}{\bm{B}}\cdot\mathrm{d}{\bm{l}},\\\ \end{split}$ (1)

where ${\bm{j}}$ is the electric current density, ${\bm{B}}$ is the magnetic
field, $A$ is an area on the tube cross section, $\hat{\bm{n}}$ is its unit
normal vector, $C$ is any surface with area $A=\int_{C}dA$, $\partial C$ is
the perimeter of $C$, and ${\bm{l}}$ is its tangent vector. For simplicity, we
assume unity magnetic permeability. Suppose a magnetic flux tube is compact in
space. More specifically, there is a closed loop $\partial C$ such that $B=0$
all along the loop. This shows that a compact flux tube has $I=0$, i.e. it is
current-neutralized.

The argument implies neutralized electric currents in the axial direction. For
a tube with coherent twist, ${\bm{j}}$ near the tube center will point in one
direction when projected onto the axis. Conversely, ${\bm{j}}$ in the tube
periphery must point in the opposite direction, forming a sheath or a skin
layer. The former is termed the “direct current” (DC), and the latter the
“return current” (RC; e.g., Melrose, 1995). The flux tube is “current
neutralized” when DC and RC add to zero.

Solar active regions (ARs), formed through emergence of subsurface flux tubes,
are known to harbor electric currents in the corona. In many cases, their
sheared or sigmoidal loops are incompatible with a current-free, potential
field morphology. Moreover, they must contain sufficient free energy to power
flares and coronal mass ejections (CMEs). This picture is quite different from
the convection zone. Here, the magnetic field is volume-filling and dominates
the plasma dynamics.

In the thin layer of the photosphere, the plasma density and pressure decrease
rapidly with height; the magnetic field begins to transition from a fibril to
a volume-filling state. Maps of the photospheric ${\bm{B}}$, magnetograms, are
routinely inferred from spectropolarimetric observations. Assuming a local
Cartesian geometry and the photosphere as a thin boundary at $z=0$, the
vertical current density $j_{z}$ and the net current $I_{z}$ can be estimated
from the horizontal field component ${\bm{B}}_{h}$, similar to Equation
(eq:stokes):

$\ilabel{eq:stokesph}\begin{split}j_{z}&=\mathbf{\nabla}\times{\bm{B}}_{h},\\\
I_{z}&=\int_{C}j_{z}\,\mathrm{d}A=\oint_{\partial
C}{\bm{B}}_{h}\cdot\mathrm{d}{\bm{l}},\\\ \end{split}$ (2)

where $C$ is now a finite area lying entirely on $z=0$.

For an isolated, spatially compact sunspot, there exists a closed loop
$\partial C$ lying on $z=0$ that surrounds the sunspot such that
${\bm{B}}(z=0)=0$. This implies $I_{z}=0$. That is, all currents that enter
the photosphere must return to the convection zone, i.e., they are “balanced”.
In practice, however, this is not how net currents are measured for observed
ARs. First of all, sunspots are rarely completely isolated. Furthermore,
eruptive ARs tend to harbor sunspots of opposite polarities in close
proximity. To measure the net current, standard practice is to take $C$ as the
contiguous area of a magnetic polarity partially bounded by the polarity
inversion line (PIL; i.e. $B_{z}=0$). Figure 1(a) of Török et al. (2014) shows
a good example. DC and RC are expected to reside in the core region and the
periphery, respectively.

###  1.2 Non-Neutralized Currents in ARs

subsec:literature

Ground-based vector magnetograms provided early evidence for non-neutralized
currents in the AR photosphere. Leka et al. (1996) followed the evolution of
an emerging AR over several days. A significant net current appears in each
magnetic polarity, whose increase cannot be explained by the shearing flows in
the photosphere alone and must be carried by the emerging flux. Wheatland
(2000) analyzed a sample of 21 ARs. In most cases, the net current for the
entire AR is consistent with zero. Within each magnetic polarity, however, it
is significantly different from zero.

Using high-resolution vector magnetograms from the Spectropolarimeter (SP) on
board the Hinode satellite, Ravindra et al. (2011) and Georgoulis, Titov, and
Mikić (2012) found that the non-neutralized currents are always accompanied by
well-formed, sheared PILs. Meanwhile, Venkatakrishnan and Tiwari (2009) and
Gosain, Démoulin, and López Fuentes (2014) found isolated sunspots without a
sheared PIL to be well neutralized. This is a natural consequence of Equation
(eq:stokesph): the sheared component of ${\bm{B}}_{h}$ along the PIL should
contribute significantly to the line integral around a single magnetic
polarity.

The current distribution has perceived important implications for solar
eruptions. Many CME models employ DC-dominated flux ropes as the driver (e.g.,
Titov and Démoulin, 1999). The inclusion of RC will reduce the outward
magnetic “hoop force” that drives the eruption (Forbes, 2010).
Observationally, non-neutralized currents indeed correlate with eruptive
activities. Liu et al. (2017), Vemareddy (2019), and Avallone and Sun (2020)
compared ARs with and without major eruptions using data from the Helioseismic
and Magnetic Imager (HMI) aboard the Solar Dynamics Observatory (SDO). The
quiescent ARs are mostly current-neutralized, while the active ones are not.
Falconer (2001) and Kontogiannis et al. (2019) demonstrated the net current as
a useful index for space weather predictions. On a related note, electric
currents can also cause an apparent line-of-sight flux imbalance because of
the directionality of the magnetic field they produce (Gary and Rabin, 1995).

The origin of the non-neutralized current in ARs is not entirely clear. If
flux emergence is responsible, it must explain how it can arise from a
current-neutralized subsurface flux tube. Longcope and Welsch (2000) proposed
an analytical model where only a fraction of the current passes into the
corona. The rest, including most of the weaker RC, is shunted to the surface
layers and hidden from the observer. Török et al. (2014) demonstrated a
similar effect using a magnetohydrodynamic (MHD) simulation from Leake,
Linton, and Török (2013). As a subsurface, neutralized flux tube emerges, non-
neutralized current develops simultaneously with the flux emergence and PIL
shear. The degree of non-neutralization becomes more pronounced in higher
layers as RC remains trapped in lower ones.

Shear and vortex surface flows are often used to induce coronal electric
currents in MHD models. Török and Kliem (2003) and Dalmasse et al. (2015)
showed that different photospheric line-tied motions can generate different
degrees of current neutralization. Only motions that directly shear the PIL
will generate a non-neutralized component. Simulations by Fan (2001) and
Manchester et al. (2004) showed that the shear flows on the PIL are driven by
the magnetic tension force during flux emergence, which comes from the twisted
tube itself as it expands in the corona.

###  1.3 Outline

subsec:outline

Here we explore an alternative origin of the observed, non-neutralized
electric current in emerging ARs. Using an analytical flux tube model, we show
that a neutralized subsurface flux tube can lead to a significant non-
neutralized vertical current on the surface due to a geometric projection
effect.

We described two versions current neutralization in Section subsec:definition:
one in the context of the subsurface flux tube itself, and the other in the
context of the photospheric observation. As we shall see, they refer to
different components of the current system when the tube axis obliquely
intersects the surface. Such a distinction is seldom discussed in the
literature, but is central to this work.

Below, we describe the flux tube model in Section sec:model, present our
results in Section sec:result, and discuss the implications in Section
sec:discussion. Additional details of the flux tube model are available in the
appendix.

Figure 1.: Magnetic field of a partially emerged FG03 flux tube ($R=0.5$,
$a=0.12$, $q=0.8$, $h=-0.5$). The tube axis (dark gray curve) just starts to
touch the photosphere at $z=0$. Maps on two horizontal slices at $z=-0.5$ and
$0$ show $B_{z}$ distribution. Selective field lines are traced from a
distance of $w=0.08$ (yellow) and $0.18$ (green) from the tube axis in the
$x=0$ plane. They are located within the DC core and RC sheath, respectively.
The black and red arrows show the radial vector for an arbitrary point in the
spherical and torus coordinate system, respectively (see appendix). The inset
illustrates the spherical coordinate $(r,\theta,\phi)$, and the torus
coordinate $(w,\psi,\phi)$. The poloidal angle $\psi$ here is negative.

f:mfr

##  2 Flux Tube Model

sec:model

###  2.1 Model Setup

subsec:setup_fg03

We consider a twisted, semi-torus flux tube model from Fan and Gibson (2003,
hereafter FG03). In a global Cartesian coordinate with unit vectors
$(\hat{\bm{x}},\hat{\bm{y}},\hat{\bm{z}})$, the subsurface tube is centered at
${\bm{r}}_{0}=(0,0,h)$, where $h<0$ (Figure f:mfr).

The semi-torus tube axis has a radius of $R$; it is axisymmetric with a polar
axis parallel to $\hat{\bm{y}}$. In a local spherical coordinate centered at
${\bm{r}}_{0}$ with unit vector
$(\hat{\bm{r}},\hat{\bm{\theta}},\hat{\bm{\phi}})$, the magnetic field
${\bm{B}}$ inside the tube is defined as

$\ilabel{eq:bsph}{\bm{B}}=\mathbf{\nabla}\times\left[\dfrac{A(r,\theta)}{r\sin{\theta}}\,\hat{\bm{\phi}}\right]+B_{\phi}(r,\theta)\,\hat{\bm{\phi}},$
(3)

where

$\ilabel{eq:bsph0}\begin{split}A(r,\theta)&=\dfrac{qa^{2}B_{t}}{2}\exp\left(-w^{2}/a^{2}\right),\\\
B_{\phi}(r,\theta)&=\dfrac{aB_{t}}{r\sin{\theta}}\exp\left(-w^{2}/a^{2}\right).\\\
\end{split}$ (4)

Here, $r$ is the radial distance from the origin; $\theta$ is the polar angle
from the polar axis; $\phi$ is the azimuth angle where 0 is along $+z$ and
increases clockwise toward $+x$ ($-\pi/2\leq\theta,\phi\leq\pi/2$). Moreover,
$w$ is the distance from the tube axis. $B_{t}$, $q$, and $a$ are constants:
$B_{t}$ is the characteristic field strength; $q$ is a measure of the twist;
$a$ is the characteristic length scale for the tube radius. Field lines near
the tube axis wind about the axis at a rate of $q/a$ per unit length. We
truncate ${\bm{B}}$ to 0 outside the surface $w=3a$. The tube is thus isolated
with a circular cross section and a radius of $3a$. Detailed expressions of
${\bm{B}}$ and ${\bm{j}}$ in various coordinate systems are available in the
appendix.

In this work, we consider the following default parameters: $B_{t}=1.0$,
$R=0.5$, and $a=0.12$. Additionally, we limit ourselves to the cases where
$q>0$, such that the flux tube has a right-handed twist.

To emulate flux emergence, we displace the flux tube upward by increasing $h$
toward zero, and evaluate ${\bm{B}}$ and ${\bm{j}}$ on the horizontal slices
$z=0$. This is equivalent to sampling the variables at different heights if
the flux tube is fixed in space. The process is demarcated by several key
values of $h$:

  * •

The emergence starts at $h=-(R+3a)=-0.86$.

  * •

The axis of the tube reaches the surface at $h=-R=-0.5$.

  * •

The crest of the tube detaches from the surface at $h=-(R-3a)=-0.14$.

  * •

The emergence ends at $h=0$.

In FG03, the tube field at $z=0$ was used to prescribe a time-dependent
boundary condition that drives the coronal dynamics. Here, we focus solely on
the layer $z=0$ as a proxy of the photospheric observations, ignoring any
coronal evolution or its feedback. This approach has been used to help
interpret the emerging non-potential field structure observed in the
photosphere with success (Lites et al., 1995; Gibson et al., 2004; Luoni et
al., 2011).

###  2.2 Poloidal-Toroidal Decomposition

subsec:ptd_fg03

The axial direction of a torus is known as the toroidal direction, which is
just $\hat{\bm{\phi}}$ in our spherical coordinate. The non-axial direction,
on the other hand, includes both $\hat{\bm{r}}$ and $\hat{\bm{\theta}}$. We
may define a torus coordinate where a non-axial vector can be decomposed into
radial and poloidal components (see Figure f:mfr and appendix). In a constant
$\phi$ plane, the former directs away from the tube axis, while the latter
follows a circle around the tube axis.

We show in the appendix that the non-axial ${\bm{B}}$ and ${\bm{j}}$ in the
FG03 flux tube are purely poloidal. Using the subscript P (T) to denote the
poloidal (toroidal) components, we have in spherical coordinate

$\ilabel{eq:btp}\begin{split}{\bm{B}}_{\mathrm{P}}&=B_{r}\hat{\bm{r}}+B_{\theta}\hat{\bm{\theta}},\\\
{\bm{B}}_{\mathrm{T}}&=B_{\phi}\hat{\bm{\phi}},\\\ \end{split}$ (5)
$\ilabel{eq:jtp}\begin{split}{\bm{j}}_{\mathrm{P}}&=j_{r}\hat{\bm{r}}+j_{\theta}\hat{\bm{\theta}},\\\
{\bm{j}}_{\mathrm{T}}&=j_{\phi}\hat{\bm{\phi}}.\\\ \end{split}$ (6)

Their respective contributions to the surface observation vary at different
stages of the flux emergence, when the tube axis intersects the surface at
different angles. In particular, the projected values in the vertical
direction can be evaluated by taking the dot product with $\hat{\bm{z}}$,
e.g., ${\bm{B}}_{\mathrm{P}}\cdot\hat{\bm{z}}$.

###  2.3 Two Versions of Current Neutralization

subsec:dcrc_fg03

Considering the flux tube by itself, we define DC and RC as components of the
net toroidal current. For $q>0$, DC (RC) refers to the net current for
$j_{\phi}>0$ ($j_{\phi}<0$). We denote the total DC as
$I_{\mathrm{T}}^{\mathrm{D}}$, the total RC as $I_{\mathrm{T}}^{\mathrm{R}}$,
and the net toroidal current $I_{\mathrm{T}}$. We can evaluate them by
integrating $j_{\phi}$ on a cross section $\phi=0$:

$\ilabel{eq:itdr}\begin{split}I^{\mathrm{D}}_{\mathrm{T}}&=\int_{\phi=0}\left.j_{\phi}\right|_{j_{\phi}>0}\,\mathrm{d}A,\\\
I^{\mathrm{R}}_{\mathrm{T}}&=\int_{\phi=0}\left.j_{\phi}\right|_{j_{\phi}<0}\,\mathrm{d}A,\\\
I_{\mathrm{T}}&=I_{\mathrm{T}}^{\mathrm{D}}+I_{\mathrm{T}}^{\mathrm{R}}.\\\
\end{split}$ (7)

For a photospheric vector magnetogram, we denote the total DC in a single
magnetic polarity as $I^{\mathrm{D}}_{z}$, the total RC as
$I^{\mathrm{R}}_{z}$, and the net vertical current $I_{z}$. For $q>0$, we can
evaluate them by integrating $j_{z}$ on $z=0$ where $B_{z}>0$:

$\ilabel{eq:izdr}\begin{split}I^{\mathrm{D}}_{z}&=\int_{z=0;B_{z}>0}\left.j_{z}\right|_{j_{z}>0}\,\mathrm{d}A,\\\
I^{\mathrm{R}}_{z}&=\int_{z=0;B_{z}>0}\left.j_{z}\right|_{j_{z}<0}\,\mathrm{d}A,\\\
I_{z}&=I^{\mathrm{D}}_{z}+I^{\mathrm{R}}_{z}.\\\ \end{split}$ (8)

We can separately evaluate the projected poloidal and toroidal contribution by
integrating ${\bm{j}}_{\mathrm{P}}\cdot\hat{\bm{z}}$ and
${\bm{j}}_{\mathrm{T}}\cdot\hat{\bm{z}}$ instead of $j_{z}$.

The degree of current neutralization are often quantified by the relative
magnitude of DC versus RC. In this study, we use the index

$\ilabel{eq:rt}R_{\mathrm{T}}=\left|\dfrac{I^{\mathrm{D}}_{\mathrm{T}}}{I^{\mathrm{R}}_{\mathrm{T}}}\right|$
(9)

for the toroidal current and

$\ilabel{eq:rz}R_{z}=\left|\dfrac{I^{\mathrm{D}}_{z}}{I^{\mathrm{R}}_{z}}\right|$
(10)

for the observed vertical current (e.g., Liu et al., 2017; Avallone and Sun,
2020). We will show that a geometric projection effect is capable of producing
large $I_{z}$ and $R_{z}$ from a flux tube with small $I_{\mathrm{T}}$ and
unity $R_{\mathrm{T}}$.

##  3 Result

sec:result

###  3.1 Flux Tube Cross Section

subsec:res_ft

For our fiducial case with $q=0.8$, the toroidal current density on the tube
cross section features a compact DC core and a diffuse RC sheath (Figure
f:itor(a)). Though not strictly axisymmetric because of the bending of the
flux tube axis, the sign change occurs near the radius $w\approx a=0.12$
(Equation (eq:jsph3)). The net toroidal current integrated within a circle of
radius $w$ peaks at $w=0.12$ due to the contribution from the DC core. It
becomes quite neutralized once the RC sheath is included (Figure f:itor(b)).
For the entire tube $w=3a=0.36$, there is
$I_{\mathrm{T}}=0.004I_{\mathrm{T}}^{\mathrm{D}}$ and $R_{\mathrm{T}}=1.004$.
Strictly speaking, an additional RC skin layer with a net current of
$-0.004I_{\mathrm{T}}^{\mathrm{D}}$ is required. We ignore this layer as it
does not affect any of our conclusions. Because $I_{\mathrm{T}}^{\mathrm{D}}$
and $I_{\mathrm{T}}^{\mathrm{R}}$ both linearly scale with $q$,
$R_{\mathrm{T}}$ is independent of the twist.

In contrast, the poloidal current density has a single sign, rotating counter-
clockwise around the axis in the same direction as the poloidal field
(Equation (eq:jtoi3); Figure f:itor(a)). Its magnitude is independent of $q$
by construct, and is comparable to the toroidal current density for $q=0.8$.

Figure 2.: Electric current on a tube cross section ($q=0.8,\phi=0$). (a)
Distribution of ${\bm{j}}$ near the tube axis. The background color shows the
toroidal current density $j_{\mathrm{T}}$; the arrows show poloidal current
density vectors ${\bm{j}}_{\mathrm{P}}$. The tube axis is marked with a cross.
The coordinate $(\Delta y,\Delta z)$ indicates the offset from the tube axis
in $y$ and $z$ direction. (b) Net toroidal current $I_{\mathrm{T}}$ within a
tube radius $w$, normalized by $I^{\mathrm{D}}_{\mathrm{T}}$ of the entire
tube.f:itor

Figure 3.: Surface magnetic maps ($q=0.8$, $z=0$). (a)–(c) Vertical magnetic
field $B_{z}$ when the tube is centered at $h=-0.5$, $-0.3$, and $0$,
respectively. Dotted curves indicate regions where $B_{z}>0$. (d)–(f)
Contribution from the poloidal field ${\bm{B}}_{\mathrm{P}}\cdot\hat{\bm{z}}$.
(g)-(i) Contribution from the toroidal field
${\bm{B}}_{\mathrm{T}}\cdot\hat{\bm{z}}$. The top row is the sum of the middle
and bottom rows. For corresponding current density maps, see Figure
f:mapsj.f:mapsb

Figure 4.: Surface current density maps ($q=0.8$, $z=0$). (a)–(c) Vertical
electric current density $j_{z}$ when the tube is centered at $h=-0.5$,
$-0.3$, and $0$, respectively. (d)–(f) Contribution from the poloidal current
${\bm{j}}_{\mathrm{P}}\cdot\hat{\bm{z}}$. (g)-(i) Contribution from the
toroidal current ${\bm{j}}_{\mathrm{T}}\cdot\hat{\bm{z}}$. Dotted curves
indicate regions where $B_{z}>0$. The top row is the sum of the middle and
bottom rows.f:mapsj

###  3.2 Synthetic Surface Observations

subsec:res_surf

The vertical magnetic field $B_{z}$ shows two elongated polarities when the
crest of the tube intersects the surface (Figure f:mapsb(a)). These are known
observationally as “magnetic tongues”, which are a manifestation of the
magnetic twist (López Fuentes et al., 2000). As $h$ increases, the tongues
retract and transition into two circles, which correspond to the two legs of
the flux tube (Figure f:mapsb(b)-(c)). The poloidal and toroidal contributions
to the vertical field are more important in the early and late stages,
respectively (Figure f:mapsb(d)-(i)).

The vertical current density $j_{z}$ shows two elongated, ribbon-shaped
polarities early on (Figure f:mapsj(a)). They mainly arise from the projected
poloidal component (Figure f:mapsj(d)). For $q=0.8$, $j_{z}$ is mostly
positive in the positive magnetic polarity, so DC dominates. As the tube
continues to rise, the RC sheath becomes visible on the surface, and the
$j_{z}$ distribution starts to resemble $j_{\phi}$ on the tube cross section
(Figure f:mapsj(g)-(i)).

Figure 5.: Integrated surface observations for different $q$. (a) Net current
$I_{z}$. (b) Ratio between DC and RC $R_{z}$. (c) Projected poloidal (thin
solid) and toroidal (dash-dotted) components in $I_{z}$. (d) Surface magnetic
flux $\Phi$ ($B_{z}>0$ only). Gray bands indicate $h\geq-(R-3a)$. Dotted lines
show the three heights in Figure f:mapsj.f:profile

Figure 6.: Surface distribution of $j_{z}$ for $h=-0.5$, when the apex of the
flux rope axis is at the photospheric level. (a)-(c) For $q=0.2$, $0.5$, and
$0.8$, respectively. Dotted curves indicate regions where $B_{z}>0$. Arrows
show ${\bm{B}}_{h}$ on the PIL.f:jz

For $q=0.8$, the surface vertical current is highly non-neutralized when the
axis just starts to emerge ($h=-R$). $I_{z}$ and $R_{z}$ peak at $h=-0.44$ and
$-0.50$, with a maximum of $I_{z}=2.6I_{\mathrm{T}}^{\mathrm{D}}$ and
$R_{z}=8.8$, respectively (Figure f:profile(a)-(b)). The dominance of DC is
clearly owing to the poloidal contribution (Figure f:profile(c)). At $h=-0.5$,
the contribution of the poloidal component in $I_{z}$ is about $18$ times that
of the toroidal counterpart.

To evaluate the general properties of the model, we consider two additional
cases for $q=0.2$ and $0.5$. Several notable features are as follows.

Patterns in surface map.— For $h=-0.5$, the $j_{z}$ and $B_{z}$ polarities
become better aligned as $q$ increases (Figure f:jz). This is an interesting
feature of the FG03 flux tube, where ${\bm{j}}_{\mathrm{P}}$ and
${\bm{B}}_{\mathrm{P}}$ depend differently on $q$. The predominant poloidal
current contribution ${\bm{j}}_{\mathrm{P}}\cdot\hat{\bm{z}}$ is independent
of $q$ (Equation (eq:jtpz)), whereas the toroidal contribution
${\bm{j}}_{\mathrm{T}}\cdot\hat{\bm{z}}$ merely appends a curved tail to the
$j_{z}$ ribbon, creating a “Yin-Yang” pattern. Therefore, the $j_{z}$ patterns
remain almost unchanged with $q$. Contrarily, the poloidal field contribution
${\bm{B}}_{\mathrm{P}}\cdot\hat{\bm{z}}$ linearly scales with $q$ (Equation
(eq:btpz)). As $q$ increases, the $B_{z}$ magnetic tongues rotate clockwise to
better match the current ribbons. The PIL becomes largely parallel to the
$j_{z}$ ribbons for $q=0.8$; the horizontal field ${\bm{B}}_{h}$ is also more
sheared. As $h$ increases further, U-shaped field lines in the lower portion
of the tube will graze the surface to create “bald patches”, where
${\bm{B}}_{h}$ on the PIL directs from the negative to the positive side
(Titov, Priest, and Demoulin, 1993).

Non-neutralized current vs. twist.— The maximum $I_{z}$ and $R_{z}$ both
positively correlate with $q$ (Figure f:profile(a)-(b)). This is because the
$B_{z}$ and $j_{z}$ polarities become better aligned at higher $q$: the area
occupied by DC increases. The $R_{z}$ values appear unexpectedly low for
$q=0.8$ early on ($h<-0.8$). We find that the toroidal contribution in $I_{z}$
is initially negative (Figure f:profile(c)). This RC component is larger for
$q=0.8$, leading to an overall smaller $R_{z}$. It is soon offset by the
faster increasing DC component in the poloidal current.

Evolution of current and magnetic flux.— As $h$ increases, $I_{z}$ and $R_{z}$
initially increase along with the surface magnetic flux $\Phi$ (Figure
f:profile(d)). The former peak at $h\approx-R$, when a significant portion of
the toroidal magnetic flux is still below the surface. As emergence proceeds,
$I_{z}$ and $R_{z}$ will decrease due to the toroidal contribution, whereas
$\Phi$ will continue to increase. For an untwisted flux tube, the maximum
$\Phi$ is expected to occur at $h=-(R-3a)=-0.14$, when the crest of the tube
fully emerges. Such is also the case for both $q=0.2$ and $0.5$. For $q=0.8$,
however, $\Phi$ peaks significant earlier at $h=-0.37$ with a maximum 5%
higher than the total toroidal flux. This is owing to the strong poloidal
field contribution (Equation (eq:btpz)). A portion of highly twisted magnetic
field lines will intersect the $z=0$ plane more than two times (appendix).

We conclude that the observed, non-neutralized vertical current (large $I_{z}$
and $R_{z}$) can arise from a neutralized flux tube (zero $I_{\mathrm{T}}$ and
unity $R_{\mathrm{T}}$) due to the projected poloidal component. The
conclusion holds for the early stages of emergence only, when the crest of the
tube is partially emerged, and the axis obliquely intersects the surface. For
the later stages, the tube axis becomes more perpendicular to the surface. The
observed current becomes more neutralized as the toroidal component takes
over.

##  4 Discussion

sec:discussion

For an idealized, torus-shaped flux tubes, as long as $B_{\mathrm{T}}$
decreases away from the tube axis, the surface
${\bm{j}}_{\mathrm{P}}\cdot\hat{\bm{z}}$ distribution is expected to be
bipolar like FG03. Whether the surface current is neutralized will depend on
the patterns of ${\bm{B}}_{\mathrm{P}}\cdot\hat{\bm{z}}$, which is sensitive
to the global magnetic twist (Figure f:jz). Other conditions being equal, a
flux tube with higher twist should possess higher magnetic helicity and free
energy, which are known to drive solar eruptions. This, combined with the
positive correlation between $R_{z}$ and $q$, provides a natural explanation
to the increased eruptivity in ARs with significant non-neutralized currents.

We temper our conclusion with some caveats. MHD simulations have shown that it
is difficult for a twisted flux tube to rise bodily into the corona (e.g.,
Fan, 2001; Archontis and Hood, 2010). The concave portions of the field lines
are loaded with dense plasma, which prevents the tube axis from rising for
more than a few scale heights above the surface. Indeed, the tongue-shaped
$B_{z}$ and $j_{z}$ polarities (Figures f:mapsb(a)) frequently appear in
observations (e.g., Poisson et al., 2016) without developing further into two
circular polarities. If so, the $h\leq-R$ phase in our model should be more
relevant than $h>-R$ in realistic settings. In these early stages of flux
emergence, the poloidal component is expected to play a prominent role. We
note that counterexamples also exist where the axial flux fully emerges in
simulation (e.g., Hood et al., 2009).

In observations, CME-productive ARs frequently exhibit long-lived, narrow DC
ribbon pairs that are well aligned with the PIL similar to the $q=0.8$ case
(see, e.g., Figure 1 of Liu et al., 2017 and Figure 3 of Avallone and Sun,
2020). The overall patterns of $j_{z}$, however, are filamentary except for
the DC ribbons (if present); evidence for a single shell of RC surrounding the
DC is also lacking. The observed $R_{z}$ evolution sometimes show a fast
increase with $\Phi$ and peaks earlier than $\Phi$ (see, e.g., Figure 2 of
Vemareddy, 2019 and Figure 3 of Avallone and Sun, 2020). Their absolute
values, however, are considerably lower than in our model (Figure
f:profile(b)). The maximum $R_{z}$ is typically $\sim$2 for CME-productive
ARs, and at best reaches $\sim$4 when the central flux rope structure is
largely isolated (Liu et al., 2017).

It is worth mentioning that in ideal MHD, the evolution of ${\bm{j}}$ is
fundamentally different from that of ${\bm{B}}$. There is no “frozen-in”
theorem for ${\bm{j}}$ so $I$ is not conserved in a flux tube;
$\partial{\bm{j}}/\partial t$ will depend not only on the velocity field
${\bm{v}}$ and ${\bm{j}}$ itself, but ${\bm{B}}$ as well (e.g., Inverarity and
Titov, 1997). Predicting the behavior of the electric current can be rather
difficult without solving the MHD equations for the full system (Parker,
1996).

Our simplistic approach neglects a rich variety of physics during flux
emergence. We discuss several possible effects below.

  * •

As the flux tube approaches the surface, it will be flattened into a sheet-
like structure due to the diminishing pressure scale height at the top of the
convection zone (Spruit, Title, and van Ballegooijen, 1987; Cheung et al.,
2008). The internal ${\bm{B}}$ and ${\bm{j}}$ distributions will change.

  * •

The drastic expansion of the emerged flux due to magnetic buoyancy
instabilities (section 3.3 of Cheung and Isobe, 2014) can shunt the RC to the
interface layer, as demonstrated by Longcope and Welsch (2000) and Török et
al. (2014).

  * •

The expansion will also induce a Lorentz force, which drives twist or shear
flows on the surface (Longcope and Welsch, 2000; Manchester et al., 2004; Fang
et al., 2012a). These flows are capable of producing non-neutralized currents
if they act close to the PIL (Dalmasse et al., 2015).

  * •

The emerged flux may undergo internal or external reconnection to produce a
new flux rope in the corona (Archontis and Török, 2008; Leake, Linton, and
Török, 2013; Takasao et al., 2015; Syntelis et al., 2019; Toriumi and Hotta,
2019; Toriumi et al., 2020). This evolution in the magnetic field is reflected
in changes of the current distribution.

  * •

Interactions between different emerging flux regions (e.g. Chintzoglou et al.,
2019), and between existing sunspots and newly emerged regions (Fan, 2016;
Cheung et al., 2019), are known to cause flares. The foot prints of the
reconnecting quasi-separatrix layers (QSLs; e.g., Démoulin, Priest, and Lonie,
1996), demarcated by flare ribbons, often overlap with the pre-flare current
ribbons. Observations haven shown that the local $j_{z}$ can greatly enhance
during the eruption (Janvier et al., 2014).

  * •

The buffeting effect from the convective flows may fragment the magnetic flux
(Cheung et al., 2008; Martínez-Sykora, Hansteen, and Carlsson, 2008; Stein and
Nordlund, 2012; Chen, Rempel, and Fan, 2017) and contribute to the filamentary
structure of $j_{z}$ (Yelles Chaouche et al., 2009; Fang et al., 2012b).

Due to these effects, the observed surface ${\bm{B}}$ and ${\bm{j}}$ are
unlikely to be identical to those on a subsurface horizontal slice. A
distinction between poloidal and toroidal components with respect to their
subsurface origin can be difficult in practice. Nevertheless, the simple model
here is capable of reproducing several key observational features. The
projection effect may be important to the photospheric current distribution,
in particular during the early stages of flux emergence.

In FG03, the flux tube is kinematically lifted from subsurface ($z<0$) at a
constant velocity ${\bm{v}}_{0}=v_{0}\hat{\bm{z}}$. As the tube intersects the
$z=0$ plane, its magnetic field ${\bm{B}}_{0}={\bm{B}}(z=0)$ and
${\bm{v}}_{0}$ induce an ideal electric field
${\bm{E}}\propto-{\bm{v}_{0}}\times{\bm{B}}_{0}$, which is used as an evolving
boundary condition for the domain of interest ($z>0$). The initial potential
field in the $z>0$ volume interacts with the newly emerged flux, resulting in
a sigmoidal flux rope that eventually erupts.

For simplicity, we place the origin at ${\bm{r}}_{0}=(0,0,0)$ in the Cartesian
coordinate. A horizontal slice at $z_{0}>0$ represents the photospheric
condition when the tube is centered at a depth of $h=-z_{0}$ (Figure f:mfr).

We define two auxiliary variables: $\rho$ for the projected distance of a
point ${\bm{r}}$ to the origin on the $y=0$ plane, and $w$ for the distance to
the tube axis:

$\ilabel{eq:aux}\begin{split}\rho&=\left(x^{2}+z^{2}\right)^{1/2},\\\
w&=\left[(\rho-R)^{2}+y^{2}\right]^{1/2}.\\\ \end{split}$ (11)

In our convention, both the Cartesian and spherical coordinate systems are
right-handed (Figure f:mfr). The coordinates $(x,y,z)$ and $(r,\theta,\phi)$
are related as

$\ilabel{eq:car2sp}\begin{split}x&=r\sin{\theta}\sin{\phi},\\\
y&=r\cos{\theta},\\\ z&=r\sin{\theta}\cos{\phi}.\\\ \end{split}$ (12)

Their unit vectors $(\hat{\bm{x}},\hat{\bm{y}},\hat{\bm{z}})$ and
$(\hat{\bm{r}},\hat{\bm{\theta}},\hat{\bm{\phi}})$ are related as

$\ilabel{eq:ecar2sph}\begin{split}\hat{\bm{x}}&=\sin{\theta}\sin{\phi}\,\hat{\bm{r}}+\cos{\theta}\sin{\phi}\,\hat{\bm{\theta}}+\cos{\phi}\,\hat{\bm{\phi}},\\\
\hat{\bm{y}}&=\cos{\theta}\,\hat{\bm{r}}-\sin{\theta}\,\hat{\bm{\theta}},\\\
\hat{\bm{z}}&=\sin{\theta}\cos{\phi}\,\hat{\bm{r}}+\cos{\theta}\cos{\phi}\,\hat{\bm{\theta}}-\sin{\phi}\,\hat{\bm{\phi}}.\\\
\end{split}$ (13)

In the spherical coordinate, the three components of the magnetic field
${\bm{B}}$ are

$\ilabel{eq:bsph3}\begin{split}B_{r}&=B_{t}\dfrac{qR\cos{\theta}}{r\sin{\theta}}\exp\left(-w^{2}/a^{2}\right),\\\
B_{\theta}&=B_{t}\dfrac{q(r-R\sin{\theta})}{r\sin{\theta}}\exp\left(-w^{2}/a^{2}\right),\\\
B_{\phi}&=B_{t}\dfrac{a}{r\sin{\theta}}\exp\left(-w^{2}/a^{2}\right).\\\
\end{split}$ (14)

The three Cartesian components are

$\ilabel{eq:bcar3}\begin{split}B_{x}&=B_{t}\dfrac{qxy+az}{\rho^{2}}\exp\left(-w^{2}/a^{2}\right),\\\
B_{y}&=B_{t}q\dfrac{R-\rho}{\rho}\exp\left(-w^{2}/a^{2}\right),\\\
B_{z}&=B_{t}\dfrac{qyz-ax}{\rho^{2}}\exp\left(-w^{2}/a^{2}\right).\end{split}$
(15)

In the spherical coordinate, the three components of the current density
${\bm{j}}$ are

$\ilabel{eq:jsph3}\begin{split}j_{r}&=\dfrac{2B_{t}R\cos\theta}{a\rho}\exp\left(-w^{2}/a^{2}\right),\\\
j_{\theta}&=\dfrac{2B_{t}(r-R\sin{\theta})}{a\rho}\exp\left(-w^{2}/a^{2}\right),\\\
j_{\phi}&=-\dfrac{B_{t}q}{\rho}\left(\dfrac{2w^{2}}{a^{2}}-\dfrac{R}{\rho}-1\right)\exp\left(-w^{2}/a^{2}\right).\end{split}$
(16)

The three Cartesian components are

$\ilabel{eq:jcar3}\begin{split}j_{x}&=-\dfrac{B_{t}}{\rho^{2}}\left\\{q\left[\dfrac{2w^{2}}{a^{2}}-\dfrac{R}{\rho}-1\right]z-\dfrac{2xy}{a}\right\\}\exp\left(-w^{2}/a^{2}\right),\\\
j_{y}&=\dfrac{2B_{t}(R-\rho)}{a\rho}\exp\left(-w^{2}/a^{2}\right),\\\
j_{z}&=\dfrac{B_{t}}{\rho^{2}}\left\\{q\left[\dfrac{2w^{2}}{a^{2}}-\dfrac{R}{\rho}-1\right]x+\dfrac{2yz}{a}\right\\}\exp\left(-w^{2}/a^{2}\right).\\\
\end{split}$ (17)

A torus coordinate $(w,\psi,\phi)$ is related to the Cartesian coordinate as

$\ilabel{eq:car2toi}\begin{split}x&=(R+w\cos{\psi})\sin{\phi},\\\
y&=-w\sin{\psi},\\\ z&=(R+w\cos{\psi})\cos{\phi}.\\\ \end{split}$ (18)

Here $w$ is the radius with respect to the tube axis, and $\phi$ is the
toroidal angle with respect to the center of the torus. They are identically
defined as in the spherical coordinate. The new coordinate $\psi$ defines the
poloidal angle, which is the polar angle with respect to the tube axis. It is
0 in the $x$-$z$ plane, and increase in the same direction as $\theta$ (Figure
f:mfr).

The unit vectors $(\hat{\bm{w}},\hat{\bm{\psi}})$ and
$(\hat{\bm{r}},\hat{\bm{\theta}})$ are related as

$\ilabel{eq:etoi2sph}\begin{split}\hat{\bm{w}}&=\dfrac{r^{2}-R\rho}{wr}\,\hat{\bm{r}}-\dfrac{Ry}{wr}\,\hat{\bm{\theta}},\\\
\hat{\bm{\psi}}&=\dfrac{Ry}{wr}\,\hat{\bm{r}}+\dfrac{r^{2}-R\rho}{wr}\,\hat{\bm{\theta}}.\\\
\end{split}$ (19)

In this definition, the axial component is toroidal ($\hat{\bm{\phi}}$). The
non-axial components include the radial ($\hat{\bm{w}}$) and poloidal
($\hat{\bm{\psi}}$) contributions. Using Equations (eq:bsph3), (eq:jsph3), and
(eq:etoi2sph), we find

$\ilabel{eq:btoi3}\begin{split}B_{w}&=0,\\\
B_{\psi}&=\dfrac{qwB_{t}}{\rho}\exp\left(-w^{2}/a^{2}\right),\\\ \end{split}$
(20) $\ilabel{eq:jtoi3}\begin{split}j_{w}&=0,\\\
j_{\psi}&=\dfrac{2wB_{t}}{a\rho}\exp\left(-w^{2}/a^{2}\right).\\\ \end{split}$
(21)

This indicates that the non-axial ${\bm{B}}$ and ${\bm{j}}$ are purely
poloidal.

Using Equations (eq:btp), (eq:jtp), (eq:ecar2sph), (eq:bsph3), and (eq:jsph3),
we can evaluate the contributions of poloidal and toroidal components in
$B_{z}$ and $j_{z}$:

$\ilabel{eq:btpz}\begin{split}{\bm{B}}_{\mathrm{P}}\cdot\hat{\bm{z}}&=B_{t}\dfrac{qyz}{\rho^{2}}\exp\left(-w^{2}/a^{2}\right),\\\
{\bm{B}}_{\mathrm{T}}\cdot\hat{\bm{z}}&=-B_{t}\dfrac{ax}{\rho^{2}}\exp\left(-w^{2}/a^{2}\right),\\\
\end{split}$ (22)
$\ilabel{eq:jtpz}\begin{split}{\bm{j}}_{\mathrm{P}}\cdot\hat{\bm{z}}&=\dfrac{2B_{t}yz}{a\rho^{2}}\exp\left(-w^{2}/a^{2}\right),\\\
{\bm{j}}_{\mathrm{T}}\cdot\hat{\bm{z}}&=\dfrac{B_{t}qx}{\rho^{2}}\left(\dfrac{2w^{2}}{a^{2}}-\dfrac{R}{\rho}-1\right)\exp\left(-w^{2}/a^{2}\right).\\\
\end{split}$ (23)

Figure 7.: Properties of a FG03 flux tube on a cross section ($q=0.8$,
$\phi=x=0$). (a) Map of toroidal current density $j_{\phi}$ (left) and Lorentz
force ${\bm{f}}$ (right). The contours are for $j_{\phi}=-0.45$, $0$, and $2$.
The background colors show the force amplitude; the arrows show the force
vectors at $w=0.12$. The tube axis is marked with a cross. (b) Map of the
length $L$ (left) and the number of turns $T$ around the axis (right) for
field lines passing through the $x=0$ plane and end at $\Delta z=-0.5$.
Horizontal dotted line shows the photosphere for $h=-0.37$ when the observed
magnetic flux reaches maximum. A portion of field lines will cross the
photosphere four times; one such field line is shown (for $x>0$). The wedge-
shaped shaded region denotes their intersections with $x=0$.

f:cutx0

The semi-torus tube is axisymmetric about the $y$-axis. However, on a cross
section with constant $\phi$, variables are not axisymmetric about the tube
axis. For example, for $j_{\phi}$, the inner, lower portion of the RC sheath
appears to be stronger than other parts (Figure f:cutx0(a)).

We evaluate the Lorentz force density ${\bm{f}}={\bm{j}}\times{\bm{B}}$ on the
cross section $\phi=0$. As both the radial field ($B_{w}$) and radial current
density ($j_{w}$) are zero, ${\bm{f}}$ must be purely radial. The force points
away from the axis so the tube tends to expand. It is the strongest at the
lower portion of the DC core (Figure f:cutx0(a)).

We trace field lines in both directions from the $\phi=0$ plane to investigate
their lengths $L$ and turns of twist $T$ around the axis (Figure f:cutx0(b)).
The latter is evaluated from the change of $\psi$ in the torus coordinate. For
$q=0.8$, the axis of the semi-torus has $L=\pi R=1.57$; field lines passing
through the lower edge are more than twice in length due to twist. Field lines
close to the axis has $T=qR/2a=1.67$; those passing slightly below the axis
have a larger $T$ due to increased $L$.

For the $q=0.8$ case, we find that a portion of the field lines may intersect
the photosphere four times (Figure 8(b)). These field lines cross the $\phi=0$
plane in a wedge-shaped region below the tube axis, where $L$ and $T$ are both
large. This is a result of competition between the curvatures given by the
flux tube axis and by the twist in the flux tube. No such field lines are
present for the $q=0.2$ and $0.5$ cases.

#### Acknowledgments

XS is partially supported by NSF award 1848250 and NASA award 80NSSC190263.
MCMC acknowledges support from NASA award 80NSSC19K0088. This work was
initiated during a scientific visit to the Institute for Astronomy by MCMC
with travel support from NSF award 1848250. We thank T. Török, Y. Liu, and V.
S. Titov for their valuable comments. Visualization softwares include the
Plotly Python package.

## References

  * Archontis and Hood (2010) Archontis, V., Hood, A.W.: 2010, Flux emergence and coronal eruption. A&A 514, A56. DOI. ADS.
  * Archontis and Török (2008) Archontis, V., Török, T.: 2008, Eruption of magnetic flux ropes during flux emergence. A&A 492, L35. DOI. ADS.
  * Avallone and Sun (2020) Avallone, E.A., Sun, X.: 2020, Electric Current Neutralization in Solar Active Regions and Its Relation to Eruptive Activity. ApJ 893, 123\. DOI. ADS.
  * Chen, Rempel, and Fan (2017) Chen, F., Rempel, M., Fan, Y.: 2017, Emergence of Magnetic Flux Generated in a Solar Convective Dynamo. I. The Formation of Sunspots and Active Regions, and The Origin of Their Asymmetries. ApJ 846, 149\. DOI. ADS.
  * Cheung and Isobe (2014) Cheung, M.C.M., Isobe, H.: 2014, Flux Emergence (Theory). Living Rev. Sol. Phys. 11, 3\. DOI. ADS.
  * Cheung et al. (2008) Cheung, M.C.M., Schüssler, M., Tarbell, T.D., Title, A.M.: 2008, Solar Surface Emerging Flux Regions: A Comparative Study of Radiative MHD Modeling and Hinode SOT Observations. ApJ 687, 1373\. DOI. ADS.
  * Cheung et al. (2019) Cheung, M.C.M., Rempel, M., Chintzoglou, G., Chen, F., Testa, P., Martínez-Sykora, J., Sainz Dalda, A., DeRosa, M.L., Malanushenko, A., Hansteen, V., De Pontieu, B., Carlsson, M., Gudiksen, B., McIntosh, S.W.: 2019, A comprehensive three-dimensional radiative magnetohydrodynamic simulation of a solar flare. Nature Astronomy 3, 160\. DOI. ADS.
  * Chintzoglou et al. (2019) Chintzoglou, G., Zhang, J., Cheung, M.C.M., Kazachenko, M.: 2019, The Origin of Major Solar Activity: Collisional Shearing between Nonconjugated Polarities of Multiple Bipoles Emerging within Active Regions. ApJ 871, 67\. DOI. ADS.
  * Dalmasse et al. (2015) Dalmasse, K., Aulanier, G., Démoulin, P., Kliem, B., Török, T., Pariat, E.: 2015, The Origin of Net Electric Currents in Solar Active Regions. ApJ 810, 17\. DOI. ADS.
  * Démoulin, Priest, and Lonie (1996) Démoulin, P., Priest, E.R., Lonie, D.P.: 1996, Three-dimensional magnetic reconnection without null points 2. Application to twisted flux tubes. JGR 101, 7631\. DOI. ADS.
  * Falconer (2001) Falconer, D.A.: 2001, A prospective method for predicting coronal mass ejections from vector magnetograms. J. Geophys. Res. 106, 25185\. DOI. ADS.
  * Fan (2001) Fan, Y.: 2001, The Emergence of a Twisted $\Omega$-Tube into the Solar Atmosphere. ApJ 554, L111. DOI. ADS.
  * Fan (2009) Fan, Y.: 2009, Magnetic Fields in the Solar Convection Zone. Living Reviews in Solar Physics 6, 4\. DOI. ADS.
  * Fan (2016) Fan, Y.: 2016, Modeling the Initiation of the 2006 December 13 Coronal Mass Ejection in AR 10930: The Structure and Dynamics of the Erupting Flux Rope. ApJ 824, 93\. DOI. ADS.
  * Fan and Gibson (2003) Fan, Y., Gibson, S.E.: 2003, The Emergence of a Twisted Magnetic Flux Tube into a Preexisting Coronal Arcade. ApJ 589, L105. DOI. ADS.
  * Fang et al. (2012a) Fang, F., Manchester, I. Ward, Abbett, W.P., van der Holst, B.: 2012a, Buildup of Magnetic Shear and Free Energy during Flux Emergence and Cancellation. ApJ 754, 15\. DOI. ADS.
  * Fang et al. (2012b) Fang, F., Manchester, I. Ward, Abbett, W.P., van der Holst, B.: 2012b, Dynamic Coupling of Convective Flows and Magnetic Field during Flux Emergence. ApJ 745, 37\. DOI. ADS.
  * Forbes (2010) Forbes, T.: 2010, In: Schrijver, C.J., Siscoe, G.L. (eds.) Models of coronal mass ejections and flares, 159\. ADS.
  * Gary and Rabin (1995) Gary, G.A., Rabin, D.: 1995, Line-of-Sight Magnetic Flux Imbalances Caused by Electric Currents. Sol. Phys. 157, 185\. DOI. ADS.
  * Georgoulis, Titov, and Mikić (2012) Georgoulis, M.K., Titov, V.S., Mikić, Z.: 2012, Non-neutralized Electric Current Patterns in Solar Active Regions: Origin of the Shear-generating Lorentz Force. ApJ 761, 61\. DOI. ADS.
  * Gibson et al. (2004) Gibson, S.E., Fan, Y., Mandrini, C., Fisher, G., Demoulin, P.: 2004, Observational Consequences of a Magnetic Flux Rope Emerging into the Corona. ApJ 617, 600\. DOI. ADS.
  * Gosain, Démoulin, and López Fuentes (2014) Gosain, S., Démoulin, P., López Fuentes, M.: 2014, Distribution of Electric Currents in Sunspots from Photosphere to Corona. ApJ 793, 15\. DOI. ADS.
  * Hood et al. (2009) Hood, A.W., Archontis, V., Galsgaard, K., Moreno-Insertis, F.: 2009, The emergence of toroidal flux tubes from beneath the solar photosphere. A&A 503, 999\. DOI. ADS.
  * Inverarity and Titov (1997) Inverarity, G.W., Titov, V.S.: 1997, Formation of current layers in three-dimensional, inhomogeneous coronal magnetic fields by photospheric motions. J. Geophys. Res. 102, 22285\. DOI. ADS.
  * Janvier et al. (2014) Janvier, M., Aulanier, G., Bommier, V., Schmieder, B., Démoulin, P., Pariat, E.: 2014, Electric Currents in Flare Ribbons: Observations and Three-dimensional Standard Model. ApJ 788, 60\. DOI. ADS.
  * Kontogiannis et al. (2019) Kontogiannis, I., Georgoulis, M.K., Guerra, J.A., Park, S.-H., Bloomfield, D.S.: 2019, Which Photospheric Characteristics Are Most Relevant to Active-Region Coronal Mass Ejections? Sol. Phys. 294, 130\. DOI. ADS.
  * Leake, Linton, and Török (2013) Leake, J.E., Linton, M.G., Török, T.: 2013, Simulations of Emerging Magnetic Flux. I. The Formation of Stable Coronal Flux Ropes. ApJ 778, 99\. DOI. ADS.
  * Leka et al. (1996) Leka, K.D., Canfield, R.C., McClymont, A.N., van Driel-Gesztelyi, L.: 1996, Evidence for Current-carrying Emerging Flux. ApJ 462, 547\. DOI. ADS.
  * Lites et al. (1995) Lites, B.W., Low, B.C., Martinez Pillet, V., Seagraves, P., Skumanich, A., Frank, Z.A., Shine, R.A., Tsuneta, S.: 1995, The Possible Ascent of a Closed Magnetic System through the Photosphere. ApJ 446, 877\. DOI. ADS.
  * Liu et al. (2017) Liu, Y., Sun, X., Török, T., Titov, V.S., Leake, J.E.: 2017, Electric-current Neutralization, Magnetic Shear, and Eruptive Activity in Solar Active Regions. ApJ 846, L6. DOI. ADS.
  * Longcope and Welsch (2000) Longcope, D.W., Welsch, B.T.: 2000, A Model for the Emergence of a Twisted Magnetic Flux Tube. ApJ 545, 1089\. DOI. ADS.
  * López Fuentes et al. (2000) López Fuentes, M.C., Demoulin, P., Mandrini, C.H., van Driel-Gesztelyi, L.: 2000, The Counterkink Rotation of a Non-Hale Active Region. ApJ 544, 540\. DOI. ADS.
  * Luoni et al. (2011) Luoni, M.L., Démoulin, P., Mandrini, C.H., van Driel-Gesztelyi, L.: 2011, Twisted Flux Tube Emergence Evidenced in Longitudinal Magnetograms: Magnetic Tongues. Sol. Phys. 270, 45\. DOI. ADS.
  * Manchester et al. (2004) Manchester, I. W., Gombosi, T., DeZeeuw, D., Fan, Y.: 2004, Eruption of a Buoyantly Emerging Magnetic Flux Rope. ApJ 610, 588\. DOI. ADS.
  * Martínez-Sykora, Hansteen, and Carlsson (2008) Martínez-Sykora, J., Hansteen, V., Carlsson, M.: 2008, Twisted Flux Tube Emergence From the Convection Zone to the Corona. ApJ 679, 871\. DOI. ADS.
  * Melrose (1995) Melrose, D.B.: 1995, Current Paths in the Corona and Energy Release in Solar Flares. ApJ 451, 391\. DOI. ADS.
  * Parker (1996) Parker, E.N.: 1996, Comment on “Current Paths in the Corona and Energy Release in Solar Flares”. ApJ 471, 489\. DOI. ADS.
  * Poisson et al. (2016) Poisson, M., Démoulin, P., López Fuentes, M., Mandrini, C.H.: 2016, Properties of Magnetic Tongues over a Solar Cycle. Sol. Phys. 291, 1625\. DOI. ADS.
  * Ravindra et al. (2011) Ravindra, B., Venkatakrishnan, P., Tiwari, S.K., Bhattacharyya, R.: 2011, Evolution of Currents of Opposite Signs in the Flare-productive Solar Active Region NOAA 10930. ApJ 740, 19\. DOI. ADS.
  * Schüssler (1979) Schüssler, M.: 1979, Magnetic buoyancy revisited: analytical and numerical results for rising flux tubes. A&A 71, 79\. ADS.
  * Spruit, Title, and van Ballegooijen (1987) Spruit, H.C., Title, A.M., van Ballegooijen, A.A.: 1987, Is there a weak mixed polarity background field? Theoretical arguments. Sol. Phys. 110, 115\. DOI. ADS.
  * Stein and Nordlund (2012) Stein, R.F., Nordlund, Å.: 2012, On the Formation of Active Regions. ApJ 753, L13. DOI. ADS.
  * Syntelis et al. (2019) Syntelis, P., Lee, E.J., Fairbairn, C.W., Archontis, V., Hood, A.W.: 2019, Eruptions and flaring activity in emerging quadrupolar regions. A&A 630, A134. DOI. ADS.
  * Takasao et al. (2015) Takasao, S., Fan, Y., Cheung, M.C.M., Shibata, K.: 2015, Numerical Study on the Emergence of Kinked Flux Tube for Understanding of Possible Origin of $\delta$-spot Regions. ApJ 813, 112\. DOI. ADS.
  * Titov and Démoulin (1999) Titov, V.S., Démoulin, P.: 1999, Basic topology of twisted magnetic configurations in solar flares. A&A 351, 707\. ADS.
  * Titov, Priest, and Demoulin (1993) Titov, V.S., Priest, E.R., Demoulin, P.: 1993, Conditions for the appearance of “bald patches” at the solar surface. A&A 276, 564\. ADS.
  * Toriumi and Hotta (2019) Toriumi, S., Hotta, H.: 2019, Spontaneous Generation of $\delta$-sunspots in Convective Magnetohydrodynamic Simulation of Magnetic Flux Emergence. ApJ 886, L21. DOI. ADS.
  * Toriumi et al. (2020) Toriumi, S., Takasao, S., Cheung, M.C.M., Jiang, C., Guo, Y., Hayashi, K., Inoue, S.: 2020, Comparative Study of Data-driven Solar Coronal Field Models Using a Flux Emergence Simulation as a Ground-truth Data Set. ApJ 890, 103\. DOI. ADS.
  * Török and Kliem (2003) Török, T., Kliem, B.: 2003, The evolution of twisting coronal magnetic flux tubes. A&A 406, 1043\. DOI. ADS.
  * Török et al. (2014) Török, T., Leake, J.E., Titov, V.S., Archontis, V., Mikić, Z., Linton, M.G., Dalmasse, K., Aulanier, G., Kliem, B.: 2014, Distribution of Electric Currents in Solar Active Regions. ApJ 782, L10. DOI. ADS.
  * Vemareddy (2019) Vemareddy, P.: 2019, Degree of electric current neutralization and the activity in solar active regions. MNRAS 486, 4936\. DOI. ADS.
  * Venkatakrishnan and Tiwari (2009) Venkatakrishnan, P., Tiwari, S.K.: 2009, On the Absence of Photospheric Net Currents in Vector Magnetograms of Sunspots Obtained from Hinode (Solar Optical Telescope/Spectro-Polarimeter). ApJ 706, L114. DOI. ADS.
  * Wheatland (2000) Wheatland, M.S.: 2000, Are Electric Currents in Solar Active Regions Neutralized? ApJ 532, 616\. DOI. ADS.
  * Yelles Chaouche et al. (2009) Yelles Chaouche, L., Cheung, M.C.M., Solanki, S.K., Schüssler, M., Lagg, A.: 2009, Simulation of a flux emergence event and comparison with observations by Hinode. A&A 507, L53. DOI. ADS.


# Closed-form Tight Bounds and Approximations for the Median of a Gamma
Distribution

Richard F. Lyon
<EMAIL_ADDRESS>
###### Abstract

We show how to find upper and lower bounds to the median of a gamma
distribution, over the entire range of shape parameter $k>0$, that are the
tightest possible bounds of the form $2^{-1/k}(A+Bk)$, with closed-form
parameters $A$ and $B$. The lower bound of this form that is best at high $k$
stays between 48 and 50 percentile, while the uniquely best upper bound stays
between 50 and 55 percentile. We show how to form even tighter bounds by
interpolating between these bounds, yielding closed-form expressions that more
tightly bound the median. Good closed-form approximations between the bounds
are also found, including one that is exact at $k=1$ and stays between 49.97
and 50.03 percentile.

Keywords: median upper bound, median lower bound, tight bounds, closed-form,
arctan, rational function, interpolation

MSC 2020: 62E17, 62E10

##  1 Introduction

The known Laurent series for the median of a gamma distribution produces some
upper and lower bounds, and very accurate approximations at high enough values
of the shape parameter $k$, but poor results near $k=1$ and below. The known
approximate and bounding results for $k$ close to zero do not extend to high
$k$. We find new bounds that correspond to percentiles close to 50 across the
entire range of shape parameter, and apply them to making even tighter bounds
as well as better approximation formulae. An asymptotic approximation and
bound at low $k$, as opposed to the usual focus on high $k$, provides the key
functional form needed for this approach.

##  2 Problem formulation

The gamma distribution PDF is
$\frac{1}{\Gamma(k)\theta^{k}}x^{k-1}e^{-\frac{x}{\theta}}$, but we’ll use
$\theta=1$ because both the mean and median simply scale with this parameter.
Thus we use this PDF with just the shape parameter $k$, with $k>0$ and $x\geq
0$:

$p(x)=\frac{1}{\Gamma(k)}x^{k-1}e^{-x}$

The mean of this distribution, $\mu$, is well known to be $\mu=k$. The median
$\nu$ is the value of $x$ at which the CDF equals one-half:

$\frac{1}{2}=\int_{0}^{\nu}p(x)dx=\int_{0}^{\nu}\frac{x^{k-1}}{\Gamma(k)}e^{-x}dx$

This equation has no easy solution, but the median is well known to be a bit
below the mean, bounded by [6]

$k-\frac{1}{3}<\nu<k\quad\mathrm{and}\quad 0<\nu$

Bounds that are tighter in some part of the shape parameter range can be
obtained from the known Laurent series partial sums, or from the low-$k$
asymptote and bounds of Berg and Pedersen [3].

We seek upper and lower bounds that are tighter, especially in the middle part
of the $k$ range, than are previously known. Further, we seek simple
approximation formulae for the median, leveraging these bounds.

##  3 Prior work

A Laurent series for $\nu$ with rational coefficients has been discovered,
with deep connections to some math by Ramanujan. Choi [7] applied Ramanujan’s
work to this particular question, providing 4 coefficients (through the
$k^{-3}$ term). Berg and Pedersen [3], based on work by Marsaglia [8],
extended this to 10 coefficients. Neither commented on the radius of
convergence, which appears to be in the neighborhood of $k=1$. So for large
enough $k$ and $N$ the series yields excellent approximations, but for $k<1$
it is useless.

$\nu\approx k+\sum_{j=0}^{N}a_{j}k^{-j}$

with $a_{j}=\\{\frac{-1}{3},\frac{2^{3}}{3^{4}\cdot 5},\frac{2^{3}\cdot
23}{3^{6}\cdot 5\cdot 7},\frac{2^{3}\cdot 281}{3^{9}\cdot 5^{2}\cdot
7},\frac{-2^{3}\cdot 17\cdot 139753}{3^{13}\cdot 5^{3}\cdot 7\cdot
11},\frac{-2^{3}\cdot 708494947}{3^{15}\cdot 5^{3}\cdot 7^{2}\cdot 11\cdot
13}\dots\\}$. Thus, where Choi [7] had 144 instead of the correct $2^{3}\cdot
23=184$:

$\nu=k-\frac{1}{3}+\frac{8}{405k}+\frac{184}{25515k^{2}}+O\left(\frac{1}{k^{3}}\right)$

Partial sums of the Laurent series are not generally bounds, but the first two
($k$ and $k-1/3$) are upper and lower bounds, respectively, and the sums
ending with $-3$ and $-5$ powers of $k$ are also upper and lower bounds,
respectively.

Berg and Pedersen [3] also derived an asymptote for small $k$, which we call
$\nu_{0}$:

$\nu\approx\nu_{0}=e^{-\gamma}\,2^{-1/k}$

where $\gamma\approx 0.577216$ is the Euler–Mascheroni constant. This
asymptote is a lower bound, as we will show. Their factor $2^{-1/k}$ is the
key to good approximations and bounds, and we divide by it to reduce the
dynamic range in some of our later plots. For the range $0.01<k<100$, this
reduces the dynamic range we need to work with by nearly 30 orders of
magnitude—or 300 orders of magnitude for $k$ down to 0.001—but we still need
log–log plots to show the bounds and approximations, and their errors, across
these ranges.

Berg and Pedersen [3] also provide an upper bound $\nu<e^{-1/3k}k$ that is
just above $k-\frac{1}{3}$, and a lower bound $\nu>2^{-1/k}k$.

These previous known bounds are illustrated in Figure 1, where upper and lower
bounds and their errors (or “margins”) are distinguished by different line
styles.

Others have shown good bounds and approximations where $k$ is an integer, that
is, for the Erlang distribution [1, 10, 5]. These results do not extend to the
low-$k$ region.

The approach of approximating functions by interpolating between upper and
lower bounds has been discussed by Barry [2], who used minimax optimization to
find numeric parameters in an interpolation-between-bounds approximation to
the exponential integral. We are not aware of an interpolation approach being
used to find improved closed-form bounds.

Figure 1: Previously published bounds (red; lower bounds solid, upper bounds
dashed) for the median of a gamma distribution (black dotted), are good at
high $k$ or low $k$, but not both. Their margins (errors) are shown in blue
(corresponding line styles). At the left, at $k=0.01$, the median is near
$10^{-30}$, and at $k=0.001$ near $10^{-300}$.

##  4 Tight upper and lower bounds

The Berg lower bounds $2^{-1/k}e^{-\gamma}$ and $2^{-1/k}k$ are tight within
their families $2^{-1/k}A$ and $2^{-1/k}Bk$, but not tight within the wider
two-parameter family $2^{-1/k}(A+Bk)$. We show that the sum of these two is
the uniquely tight upper bound in that family, and that there is a range of
tight lower bounds in the family.

To improve on the lower bounds, and to motivate the family that we consider
further, first solve for $\nu$ in a simple approximation to the distribution’s
integral, using $e^{-x}<1$ for $x>0$:

$\frac{1}{2}<\int_{0}^{\nu}\frac{x^{k-1}}{\Gamma(k)}dx$ $\nu>2^{-1/k}\
\Gamma(k+1)^{1/k}$

which is a tight lower bound, and is a good approximation for $k<0.1$, but not
so great at high $k$—and not what we consider a closed form, due to the gamma
function. This expression resembles the “quantile mechanics” boundary
condition for the gamma distribution from Steinbrecher and Shaw [9], and
converges with the Berg asymptote at low $k$, but we can improve Berg’s result
by utilizing another term. A symbolic calculus system finds for us the next
Taylor series terms about $k=0$ for the power of the gamma function:

$\Gamma(k+1)^{1/k}\approx
e^{-\gamma}+\frac{e^{-\gamma}\pi^{2}}{12}k-0.035k^{2}$

Figure 2: The $A$–$B$ parameter space is shaded with dash-dot lines where
$\nu(k)2^{1/k}=A+Bk$, for a set of very small to very large $k$ values in
geometric progression (using numerically computed $\nu(k)$ values). Key values
of $A$ and $B$ are indicated. Points outside (or on the edge) of the shaded
region represent bounds, while points inside the shaded region represent
functions that cross the median function. There is an obvious uniquely tight
upper bound $\nu_{U}$, and a curved locus of tight lower bounds from
$\nu_{L0}$, which is tightest near $k=0$, to $\nu_{L\infty}$, which is
tightest for high $k$. One point (pentagram) on the curved locus represents a
lower bound $\nu_{L1}$ that is tight at $k=1$, for which $A+B=2\log 2$ (which
is the equation of the dashed line). The point $\nu_{1}$ represents a good
asymptotic approximation close to $\nu_{L0}$, but not a bound; see the next
figure. The dotted line from $\nu_{U}$ to $\nu_{L\infty}$ at $B=1$ intersects
the lines for all $k$ in monotonic order, with $A$ decreasing while $k$
increases.

So we have this improved approximation, which has much less relative error
than Berg’s at low $k$, and has a high-$k$ behavior nearly proportional to $k$
(but is no longer a lower bound because we made it larger by ignoring a next
negative term):

$\nu_{1}=2^{-1/k}\left(e^{-\gamma}+\frac{e^{-\gamma}\pi^{2}}{12}k\right)$

Inspired by this asymptotic approximation, we consider members of this family
of functions, with coefficients $A$ and $B$, and analyze which ones are
bounds:

$\tilde{\nu}=2^{-1/k}\,(A+B\,k)$

A graphical characterization of this family is most informative. Given some
values of $k$ and corresponding numerical $\nu(k)$, we can find the lines
$A+Bk=\nu(k)2^{1/k}$ in $A$–$B$ space, and plot them—see Figure 2. Regions
full of lines are not bounds, and regions without lines are where bounds are
found (including some of Berg’s bounds); we’re interested in the boundaries
between these regions, where tight bounds are to be found.

Figure 3: Zooming in to $\nu_{1}$ and $\nu_{L0}$, note that the point with
$B=e^{-\gamma}\pi^{2}/12$, which we got from the Taylor series of the power of
the gamma function, is actually inside the shaded area, so does not represent
a bound; but a point at slightly lower $B=0.45965$ is on the edge, so
represents a lower bound. These points give zero error at approximately
$k=0.1003$ and $k=0.0708$, respectively (see the next figure). We do not have
analytic formulations for these numeric and graphical observations. Figure 4:
The lower bounds $\nu_{L0}$, $\nu_{L1}$, and $\nu_{L\infty}$ (solid), and
upper bound $\nu_{U}$ (dashed) are shown in red over the ideal median (black
heavy dots), with their absolute errors in blue, all premultiplied by
$2^{1/k}$ to reduce the required plot range. The approximation $\nu_{1}$,
which is not a bound, is also shown; note that its error curve changes from
solid to dashed at the cusp, while the errors for $\nu_{L0}$ and $\nu_{L1}$
have $\log(0)$ cusps where the error grazes zero but does not change sign. The
$k$ parameters at these cusps correspond to the sloped lines indicated in the
previous figures.

The improved asymptote $\nu_{1}$ is great at low $k$, but is neither an upper
nor a lower bound, as shown in 3. We can modify it to approach $k-\frac{1}{3}$
at high $k$ by a few adjustments, via this asymptotic approximation that we
get from a symbolic calculus system:

$k\,2^{-1/k}=k-\log 2+O(k^{-1})$

Thus we find this approximation for high $k$, which is a lower bound as
illustrated in Figure 2:

$\nu_{L\infty}=2^{-1/k}\,\left(\log 2-\frac{1}{3}+k\right)$

A compromise approximation for low $k$ mixes these two, differing from the
high-$k$ approximation in only the $A$ coefficient, leaving a result
consistent to the same order as Berg’s at low $k$, and forming an upper bound
as illustrated in 2:

$\nu_{U}=2^{-1/k}\left(e^{-\gamma}+k\right)$

This mixed approximation has absolute and relative errors approaching zero at
low $k$, and relative error approaching zero at high $k$; but the absolute
error remains high, near $\log 2-\frac{1}{3}-e^{\gamma}\approx 0.20$, at high
$k$. These approximations and their errors are illustrated in Figure 4.

To support the graphical/numerical observation that $\nu_{L\infty}$ and
$\nu_{U}$ are lower and upper bounds, respectively, of the true median
($\nu_{L\infty}<\nu<\nu_{U}$), we examine their asymptotic behaviors in more
detail. At low $k$, it is easy to see, using $\log 2-\frac{1}{3}\approx
0.359814<e^{-\gamma}\approx 0.561459$, and
$\frac{e^{-\gamma}\pi^{2}}{12}\approx 0.461781<1$, that these differences are
positive, for $k\rightarrow 0$:

$2^{1/k}\left(\nu-\nu_{L\infty}\right)=e^{-\gamma}-\left(\log
2-\frac{1}{3}\right)+O(k)>0$
$2^{1/k}\left(\nu_{U}-\nu\right)=k\left(1-\frac{e^{-\gamma}\pi^{2}}{12}\right)+O(k^{2})>0$

At high $k$, a symbolic calculus system gives us for $\nu_{L\infty}$:

$2^{-1/k}k=k-\log 2+\frac{\log^{2}2}{2}\,k^{-1}+O(k^{-2})$

which we can use to construct comparisons to the Laurent series terms for
$\nu$ [7, 3]. Again we find positive differences, with $\frac{\log
2}{3}-\frac{\log^{2}2}{2}\approx-0.009177<\frac{8}{405}$, and $\log
2-e^{-\gamma}\approx 0.131688<\frac{1}{3}$, for $k\rightarrow+\infty$:

$\nu-\nu_{L\infty}=\left(\frac{8}{405}-\left(\frac{\log
2}{3}-\frac{\log^{2}2}{2}\right)\right)\,k^{-1}+O(k^{-2})>0$
$\nu_{U}-\nu=e^{-\gamma}-\log 2+\frac{1}{3}+O(k^{-1})>0$

In addition to these high-$k$ and low-$k$ asymptotic results, we can show the
inequalities also hold at $k=1$ where $\frac{1}{2}\left(\log
2-\frac{1}{3}+1\right)<\log 2<\frac{1}{2}(e^{-\gamma}+1)$, but otherwise we’re
relying on the graphical and numerical results. Since there is considerable
margin in the asymptotes, and the median is well behaved (unique, monotonic,
and smooth, with positive second derivative for all $k>0$ [3, 4]), this seems
reliable enough in concluding that these are bounds (but mathematicians are
invited to interpret these as conjectures to be proved). The positivity
constraints would not all hold, since higher-order terms would not cancel, if
$A_{L\infty}$ or $B_{L\infty}$ were any higher, or if $A_{U}$ or $B_{U}$ were
any lower. In that sense, these upper and lower bounds are proved tight.

For the lower bound $\nu_{L1}$ that is tight at $k=1$, both the value and the
slope need to match the true median. The value is the median of the
exponential distribution, $\nu(1)=\log 2$. The slope $v^{\prime}(k)$, which is
somewhat more troublesome to work out, but is tractable at the special point
$k=1$, is:

$v^{\prime}(k)|_{k=1}=\gamma-2\textrm{Ei}(-\log 2)-\log\log 2\approx
0.9680448$

which is a mathematical expression with a definite value, but is not a closed
form due to the exponential integral, so still requires a numerical approach
to evaluate it. Therefore, we describe this bound with approximate numerical
parameters instead of closed-form analytic expressions.

$B=2\left(v^{\prime}(k)|_{k=1}-\log^{2}2\right)\approx 0.9751836$ $A=2\log
2-B\approx 0.4111107$

The lower bound $\nu_{L0}$ is in worse shape, as we have to search for the $k$
value that gives the lowest $B$ value with $A=e^{-\gamma}$. So its $B$
parameter has no concise mathematical expression, but can be computed to high
precision: $B\approx 0.4596507$.

Table 1: Comparison of several median bounds and asymptotes. Version |  | Description
---|---|---
$\nu$ |  |  | True median of gamma distribution
$e^{-1/3k}k$ |  | Berg’s upper bound, high-$k$ asymptote
$2^{-1/k}\Gamma(k+1)^{1/k}$ |  | New lower bound, low-$k$ asymptote
| $A$ | $B$ | parameters for the form $2^{-1/k}(A+Bk)$
$2^{-1/k}k$ | $0$ | 1 | Berg’s lower bound
$\nu_{0}$ | $e^{-\gamma}$ | 0 | Berg’s asymptote, a lower bound
$\nu_{1}$ | $e^{-\gamma}$ | $\frac{e^{-\gamma}\pi^{2}}{12}$ | Improved low-$k$ asymptote; not a bound
$\nu_{L0}$ | $e^{-\gamma}$ | $0.4596507$ | New tight lower bound, best at low $k$
$\nu_{L1}$ | $0.4111107$ | $0.9751836$ | New tight lower bound, tangent at $k=1$
$\nu_{L\infty}$ | $\log 2-\frac{1}{3}$ | $1$ | New tight lower bound, best at high $k$
$\nu_{U}$ | $e^{-\gamma}$ | $1$ | New uniquely tight upper bound

The bounds and asymptotic approximations discussed here are summarized in
Table 1. In subsequent sections we focus primarily on the new upper and lower
bounds with closed-form coefficients, $\nu_{U}$ and $\nu_{L\infty}$, as a
basis for even tighter closed-form bounds.

Figure 5: The percentiles achieved by four new median bounds of the form
$2^{-1/k}(A+Bk)$ (solid curves) are plotted, along with the linear bounds $k$
and $k-\frac{1}{3}$ (dotted), upper and lower bounds from Berg and Pedersen
[3] (dash-dot), and a pair of closer bounds formed by interpolation between
$\nu_{U}$ and $\nu_{L\infty}$ using a one-parameter rational function
(dashed). The bounds $2^{1/k}k$, $\nu_{U}$, $\nu_{L\infty}$, and the
interpolated bounds converge on 50th percentile at both low and high $k$,
while the other six do not. Both the upper and lower interpolated bounds are
close to $\nu_{U}$ at low $k$ and close to $\nu_{L\infty}$ at high $k$;
tighter such interpolated bounds, developed in a later section, would crowd
the center of the graph.

Figure 5 shows the percentile values achieved by these bounds and
approximations, compared to the ideal 50% that defines the
median—$\nu_{L\infty}$ always comes in between 48% and 50% and $\nu_{U}$
always between 50% and 55% (percentiles are calculated in Matlab as
100*gammainc(x, k), using the normalized lower incomplete gamma function that
is the CDF for our PDF).

The coefficient of $k^{-1}$ for $\nu_{L\infty}$ is negative, so the lower
bound $\nu_{L\infty}$ slightly violates the $k-\frac{1}{3}$ lower bound in
spite of having asymptotically zero absolute error. That is, for $k>3.021$,
it’s a looser lower bound and a worse approximation than $k-1/3$, even though
it is the tightest lower bound of the form we’re considering. On the other
hand, $\nu_{U}$ is a much tighter upper bound than $k$ is.

##  5 Formulae for tighter bounds

Figure 6: Functions $A(k)$ and $B(k)$, either of which can solve
$\nu=2^{-1/k}(A+Bk)$, with the other constant at the limiting values indicated
by the circles, the parameters of $\nu_{U}$; modeling these curves can lead to
better bounds or approximations.

Letting $A$ and $B$ be functions of $k$, rather than constants, allows tighter
bounding expressions (and potentially exact expressions) for $\nu(k)$, but not
enough structure. Allowing only one of them to vary, and tying the other to
values used in the tight bounds above, allows a more constrained space of
bounds.

With $B_{L\infty}=B_{U}=1$, we can express the median exactly as
$\nu=2^{-1/k}(A(k)+k)$, for some smooth positive real function $A(k)$ that
runs from a limit of $A_{U}=e^{-\gamma}$ as $k\rightarrow 0$ to
$A_{L\infty}=\log 2-\frac{1}{3}$ as $k\rightarrow+\infty$; it is apparently
monotonic.

Alternatively, using $A_{L0}=A_{U}=e^{-\gamma}$, the formula
$\nu=2^{-1/k}(e^{-\gamma}+B(k)k)$ has a smooth positive but non-monotonic
function $B(k)$ that runs between limits $\frac{e^{-\gamma}\pi^{2}}{12}$ and
$1$, but drops a little below its low-$k$ limit before increasing.

This approach converts the problem of finding tighter bounds to the median to
the problem of finding closed-form expressions to bound these well-behaved
functions. Calculating $A(k)$ and $B(k)$ numerically to high precision is easy
when the median is calculated; see Figure 6. For the rest of this paper, we
focus on $A(k)$, since it is monotonic and more nearly symmetric on a $\log k$
axis, and because it corresponds to interpolation between closed-form bounds.

##  6 Interpolators

The function $A(k)$ introduced above can be represented in terms of an
interpolation function $g(k)$ that runs monotonically from 0 at low $k$ to 1
at high $k$:

$A(k)=g(k)A_{L\infty}+(1-g(k))A_{U}=A_{U}-g(k)(A_{U}-A_{L\infty})$
$=e^{-\gamma}-g(k)\left(e^{-\gamma}-\log 2+\frac{1}{3}\right)$

And $g(k)$ is therefore also the function that interpolates between the
bounds, allowing us to write the median in these convenient ways:

$\nu=g(k)\nu_{L\infty}+(1-g(k))\nu_{U}$
$\nu=2^{-1/k}\left(e^{-\gamma}-g(k)\left(e^{-\gamma}-\log
2+\frac{1}{3}\right)+k\right)$

The ideal interpolator can be found numerically:

$g(k)=\frac{A_{U}-A(k)}{A_{U}-A_{L\infty}}=\frac{\nu_{U}-\nu}{\nu_{U}-\nu_{L\infty}}$

It can be interesting to bound or otherwise approximate $g(k)$ with various
$\tilde{g}(k)$. In approximating the ideal with an interpolator
$\tilde{g}(k)$, we achieve absolute and relative error of the median estimate
approaching zero at low $k$ if $\tilde{g}(k)=0+O(k)$, and at high $k$ if
$\tilde{g}(k)=1-O(k^{-1})$. But we might want to do better, matching the
asymptotic slopes of the ideal interpolator to match the median to a higher
order; or we might want to match the exact known value $\nu=\log 2$ at $k=1$.
So we analyze these properties of the ideal interpolator, and give them names.
At low $k$:

$P_{0}=\frac{dg}{dk}=\frac{1-\frac{e^{-\gamma}\pi^{2}}{12}}{e^{-\gamma}-\log
2+\frac{1}{3}}\approx 2.66913$

At high $k$:

$P_{\infty}=-\frac{dg}{d\frac{1}{k}}=\frac{\frac{8}{405}+e^{-\gamma}\log
2-\frac{\log^{2}2}{2}}{e^{-\gamma}-\log 2+\frac{1}{3}}-\log 2\approx 0.143472$

And at $k=1$:

$P_{1}=g(1)=\frac{1+e^{-\gamma}-2\log 2}{e^{-\gamma}-\log
2+\frac{1}{3}}\approx 0.868678$

How such approximation goals relate to bounds is not immediately clear. We
construct some examples. Figure 7 shows the ideal interpolator, computed
numerically, and compares it to bounding interpolators of the forms
$\tilde{g}(k)=\frac{k}{k+b_{0}}$ and
$\tilde{g}(k)=\frac{2}{\pi}\tan^{-1}\frac{k}{b}$, with parameters $b_{0}$ and
$b$ chosen to yield tight upper and lower bounds, as discussed in the next two
sections. The effects of these interpolator bounds on the median bounds is
shown in Figure 8.

Figure 7: The ideal interpolator $g(k)$ (heavy dotted sigmoid) is compared
with upper and lower bounds $\tilde{g}(k)$; their margins $\tilde{g}(k)-g(k)$
are also plotted, magnified and displaced, with the same curve styles. The
curves with largest absolute margins (dashed), which correspond to the
interpolated bounds shown in Figure 5, are for the first-order rational-
function interpolator $\frac{k}{k+b_{0}}$, while the curves with smaller
margins (solid) are for arctan interpolators
$\frac{2}{\pi}\tan^{-1}\frac{k}{b}$. In each case, the one parameter ($b_{0}$
or $b$) is chosen to give a tight bound (analytically in closed form in three
of the four cases). Lower bounds of $g(k)$ make upper bounds of $\nu(k)$, and
vice versa. Figure 8: The absolute margins of the interpolated bounds are
smaller than those of the bounds they started from. Compare Figure 4. The
generally smallest margins are for the arctan interpolators, and the
intermediate for the $N=1$ rational-function interpolators.

##  7 Rational-function interpolators

Consider rational functions as interpolators, of the form

$\tilde{g}_{N}(k)=\frac{\sum_{n=1}^{N-1}a_{i}k^{i}+k^{N}}{\sum_{n=0}^{N-1}b_{i}k^{i}+k^{N}}$

For $N=1$, the only parameter is $b_{0}$, so we have a one-parameter family:

$\tilde{g}_{1}(k)=\frac{k}{b_{0}+k}$

For $N=2$ we have three parameters:

$g_{2}(k)=\frac{a_{1}k+k^{2}}{b_{0}+b_{1}k+k^{2}}$

and so forth.

We can easily constrain the coefficients to match the properties of the ideal
interpolator. At low $k$:

$\frac{a_{1}}{b_{0}}=P_{0}$

At high $k$

$b_{N-1}-a_{N-1}=P_{\infty}$

And at $k=1$:

$\frac{\sum_{n=1}^{N-1}a_{i}+1}{\sum_{n=0}^{N-1}b_{i}+1}=P_{1}$

For $N=1$, the low asymptote is tightly approached with
$b_{0}=\frac{1}{P_{0}}$, yielding an upper bound for the median (lower bound
for the interpolator $g_{\mathrm{ideal}}$). Or the high asymptote is tightly
approached with $b_{0}=P_{\infty}$, yielding a lower bound for the median
(upper bound for $g(k)$). These bounds are illustrated in Figure 7. For
$b_{0}$ between these values, the resulting interpolated function is not a
bound, but is exact at one value of $k$, for example at $k=1$ with
$b_{0}=\frac{1}{P_{1}}-1$.

For $N=2$, there are enough parameters to use any or all of the three
constraints, but it’s not immediately clear which sets of constraints can lead
to bounds. Certainly using all three constraints does not lead to a bound, but
to an interesting approximation. As we did with the $A$ versus $B$ space, we
can investigate the locus of parameter solutions for each $k$, and examine the
edges of these locus-filled areas for tight bounds. Reducing the space to 2D
by constraining one asymptote or the other allows a graphical approach, but
the results are not better than a one-parameter arctan interpolator.

For $N\geq 3$, excellent approximations and bounds are possible, but with so
many parameters are not very interesting.

##  8 Arctan interpolators

Getting a tighter bound or better approximation from the rational-function
family requires at least a handful of parameters. An alternative approach is
to find a one-parameter shape that fits better. We have found that the arctan
shape with parameter $b$ (like $b_{0}$ in the one-parameter rational-function
interpolator, corresponding to the $k$ value at the midpoint,
$\tilde{g}(b)=0.5$) does a good job:

$\tilde{g}_{a}(k)=\frac{2}{\pi}\tan^{-1}\frac{k}{b}$

As Figure 7 shows, the shapes of the arctan interpolators are imperfect, but
are much better than the first-order rational function, fitting better in some
regions than in others. Note that both the $N=1$ rational function and the
arctan interpolators are symmetric about their centers (with $\log k$ as the
independent variable); they are logistic sigmoid and Gudermannian shapes,
respectively. The ideal that they are to bound or approximate, however, is not
quite symmetric in $\log k$. So a family of not-quite-symmetric interpolators
can perhaps do better.

Figure 9: Relative errors of arctan-interpolated approximations (dash-dot
curves) between the upper (dashed) and lower (solid) bounds. These are among
the possibly interesting approximations suggested by the tight-bounds
approach. The non-bounding approximations, optimized for different criteria,
all have maximum relative errors below 1%.

Several special $b$ parameter values for the arctan interpolator can be
derived to match each of the ideal properties mentioned above. We can
constrain the approximation to pass through the known value $\nu=\log 2$ at
$k=1$, with $b=\cot\left(\frac{\pi}{2}P_{1}\right)$, but that does not give a
bound. Or we match the true median at high $k$ to within $O(k^{-2})$ with
$b=\frac{\pi}{2}P_{\infty}$. That also does not give a bound. Or we can match
at low $k$ to within $O(2^{-1/k}\,k^{2})$, like $\tilde{\nu}_{1}$ does, with
$b=\frac{2}{\pi}P_{0}^{-1}$, which yields a lower bound to $g$.

To find an upper bound to $g$ (lower bound to $\nu$), we decrease $b$ until
the margin is nonnegative for all $k$, which is at about $b=0.205282$. Finding
an analytic formulation for that tight bound is a challenge left to others.

Several interpolated bounds and approximations are summarized in Table 2. Most
are closed-form analytic expressions; see Figure 9 for the relative errors of
various arctan interpolator versions.

Table 2: Comparison of several one-parameter interpolated bounds and approximations $\tilde{g}(k)\nu_{L\infty}+(1-\tilde{g}(k))\nu_{U}$, some of which have closed-form parameters. Bounds are indicated by U or L. Version | Parameter
---|---
$\tilde{g}_{1}(k)=\frac{k}{k+b_{0}}$ | Symbolic $b_{0}$ | Numeric $b_{0}$ |
best at low $k$ | $\frac{e^{-\gamma}-\log 2+\frac{1}{3}}{1-\frac{e^{-\gamma}\pi^{2}}{12}}$ | 0.374654 | U
exact at $k=1$ |  $\frac{e^{-\gamma}-\log 2+\frac{1}{3}}{1+e^{-\gamma}-2\log 2}-1$ | 0.151175 | –
best at high $k$ | $\frac{\frac{8}{405}+e^{-\gamma}\log 2-\frac{\log^{2}2}{2}}{e^{-\gamma}-\log 2+\frac{1}{3}}-\log 2$ | 0.143472 | L
$\tilde{g}_{a}(k)=\frac{2}{\pi}\tan^{-1}\frac{k}{b}$ | Symbolic $b$ | Numeric $b$ |
best at low $k$ | $\frac{24}{\pi}\left(\frac{e^{-\gamma}-\log 2+\frac{1}{3}}{12-e^{-\gamma}\pi^{2}}\right)$ | 0.238512 | U
best at high $k$ |  $\frac{\pi}{2}\left(\frac{\frac{8}{405}+e^{-\gamma}\log 2-\frac{\log^{2}2}{2}}{e^{-\gamma}-\log 2+\frac{1}{3}}-\log 2\right)$ | 0.225366 | –
minimax relative error | $\operatorname*{argmin}_{r}\max\left|\frac{\nu-\tilde{\nu}}{\nu}\right|$ | 0.21639 | –
minimax absolute error | $\operatorname*{argmin}_{r}\max\left|\nu-\tilde{\nu}\right|$ | 0.21008 | –
exact at $k=1$ |  $\cot\left(\frac{\pi}{2}\cdot\frac{1+e^{-\gamma}-2\log 2}{e^{-\gamma}-\log 2+\frac{1}{3}}\right)$ | 0.209257 | –
tangent at $k\approx 0.4184$ | ? | 0.205282 | L

##  9 Conclusions

Tight upper and lower bounds to the median of the gamma distribution are
introduced. The simplest new lower bound is never below the 48th percentile,
and the simplest new upper bound, of the same form $2^{-1/k}(A+Bk)$, is never
above the 55th percentile, over the entire range of $k>0$. Using arctan and
rational-function interpolators between these bounds, two one-parameter
families of closed-form bounds and approximations to the median of a gamma
distribution are proposed.

The one-parameter rational-function family has simple closed-form formulae for
tightest upper and lower bounds, staying below 50.85 and above 49.69
percentile, respectively; higher-order rational functions can provide tighter
bounds or better approximations.

The one-parameter arctan family of interpolators is a better fit to the ideal
interpolator, and includes a version that is most accurate in the low-$k$ tail
and provides a closed-from tight upper bound, staying below 50.18 percentile.
With different $b$ parameter, several approximations in the family, including
the closed-form version that is exact at $k=1$, stay between 49.97 and 50.03
percentile. We have not found an analytic formula for the parameter that gives
the tight lower bound, which stays above 49.96 percentile, but have shown
where to find it, graphically or numerically.

The approach of interpolating between tight bounds opens the way to finding
tighter bounds and more accurate approximations, and to finding more such
families of bounds and approximations via other interpolator forms.

## Acknowledgements

The author acknowledges the helpful comments of Google colleagues, especially
Pascal Getreuer, Srinivas Vasudevan, and Dan Piponi.

## References

  * [1] J. Adell and P. Jodrá.  On a Ramanujan equation connected with the median of the gamma distribution.  Transactions of the American Mathematical Society, 360(7):3631–3644, 2008.
  * [2] D. Barry, J.-Y. Parlange, and L. Li.  Approximation for the exponential integral (Theis well function).  Journal of Hydrology, 227(1-4):287–291, 2000.
  * [3] C. Berg and H. L. Pedersen.  The Chen–Rubin conjecture in a continuous setting.  Methods and Applications of Analysis, 13(1):63–88, 2006.
  * [4] C. Berg and H. L. Pedersen.  Convexity of the median in the gamma distribution.  Arkiv för Matematik, 46(1):1–6, 2008.
  * [5] C.-P. Chen.  The median of gamma distribution and a related Ramanujan sequence.  The Ramanujan Journal, 44(1):75–88, 2017.
  * [6] J. Chen and H. Rubin.  Bounds for the difference between median and mean of gamma and Poisson distributions.  Statistics & Probability Letters, 4(6):281–283, 1986.
  * [7] K. P. Choi.  On the medians of gamma distributions and an equation of Ramanujan.  Proceedings of the American Mathematical Society, 121(1):245–251, 1994.
  * [8] J. C. Marsaglia.  The incomplete gamma function and Ramanujan’s rational approximation to $e^{x}$.  Journal of Statistical Computation and Simulation, 24(2):163–168, 1986.
  * [9] G. Steinbrecher and W. T. Shaw.  Quantile mechanics.  European Journal of Applied Mathematics, 19(2):87–112, 2008.
  * [10] X. You.  Approximation of the median of the gamma distribution.  Journal of Number Theory, 174:487–493, 2017.


Photomapping Using
Aerial Vehicle in USARSim
Su Kim

###### Abstract

Creating a photomap plays a critical role in navigation. Therefore, flying
vehicles are usually used to create topdown maps of the environment. In this
report we used two different aerial vehicles to create a map in a simulated
environment

##  1 Introduction to USARSim( Unified System for Automation and Robot
Simulation)

USARSim is a high fidelity simulation of urban search and rescue (USAR) robots
and environments intended as a research tool for the study of human-robot
interaction (HRI) and multirobot coordination. USARSim is designed as a
simulation companion to the National Institute of Standards and Technology
(NIST) reference test facility for autonomous mobile robots for urban search
and rescue[1].
As the time passed, the USARSim has been used for broader range of robot
simulations. The acronym was converted to Unified System for Automation and
Robot Simulation. But it is still USARSim.

Aerial vehicles in this environment play a very important role[2, 6]. They act
as the eyes of a bird, flying in the sky and informing the other robots about
the environment. The better they take pictures, the better other robots will
cooperate. In the other hand two main questions arise. The first one is how to
merge these pictures to have a general map. The second question is what to do
with this map. In the following parts we will talk about these issues.

There are two aerial vehicles implemented in the USARSim, a Blimp named as
Passarola and a Quadcopter named as Airobot as shown in Figure 1. Both of
these vehicles have the ability to carry a camera to take some pictures.

The best and easiest way of using sequence of pictures is to create a photomap
from them. This photomap would be very useful for localization purposes.
Creating such a map is not easy and takes a long time to develop an algorithm
to handle it. There are several image processing algorithms used for
registering sequence of images in order to create a photo map. One of these
algorithms that uses Fourier Melin Transform has been developed by Dr Bulow is
used in our system. Even if there is such an algorithm, we still have some
problems.

One of the challenging problems in having a photomap is having a decent
sequence of images that an image processing algorithm can handle[2, 5]. One of
these problems is movements of the camera. According to the movements of
camera the image processing algorithm has to do some scaling, rotation or
translation processes. The less the changes, the better the final map will be.
Intrinsically, Passarola has a huge amount of inertia that makes its movements
smoother. That makes it move without abrupt changes in altitude or position.
In the better word, its movements along X, Y, and Z axises will be smoother.
Thus we chose passarola for carrying the camera.

In the following parts of this report we will explain what we have done in
detail.

Figure 1: A simulated Passarola on the left and an Airobot on the right side

##  2 What we have done with the source code

There is an implementation of FMI in Matlab and another in C++ [6, 7]. The C++
one was in the experements folder of Jacobotics source code and was basically
developed to work with the Philips webcam in order to test some parts of FMI
algorithm. Since that camera was a real camera and there are some issues
regarding calibration, they had to be considered during photomaping. When we
are working in a simulated area, no calibration is required[8]. The first task
we did was to change the FMI source code in C++ language. So we added a new
parameter to its constructor under the name of calibration. If it is set in
the constructor, calibration will be done. Otherwise it won’t do calibration.

The second task that we did was to move the FMI source code from experimnets
to the library folder, and set the pathes and links. Then everybody can use it
easily. At the moment it is in the lib/mapping/FMI folder.

Then we configured a Blimp to have a downlooking camera and some other things
like propellers etc. We added two important sections in Jacobotics source
code. One of them is an actuator to drive our blimp, and the other is an image
grabber to take pictures and create a photomap of them. The actuator is inside
/lib/virtual/actuators/PassarolaDrive. For using this package it is required
to create an object from this type and use setMotorSpeeds method. It receives
three parameters and returns nothing. This method generally sets the motor
speed of one of the propellers. The parameters are xZAngle, thrustPropeller
and tailPropeller
Where
xZAngle, float, is the rotation angle of the support thrust motors bars, that
make possible change the altitude of the robot (i.e up/down). The value is the
absolute rotation angle, in radians per second.
thrustPropeller, float, is the module of the velocity vector to be applied by
the front thrusters, to move the robot in the X0Z plane (i.e forward/backward
and up/down as the value of XZAngle). The value is the absolute linear
velocity, in meters per second.
tailPropeller, float, is the rotational velocity (i.e left/right). The value
is the absolute rotational velocity, in meters per second.

The other section that I added was PassarolaImageGrabber. It is placed in
lib/virtual/autonomy/PassarolaImageGrabber. This part is responsible for
getting pictures from the image server and store them in files. These files
could be used for testing FMI algorithm in offline mode, but it could be
easily disabled. Since we are in the phase of testing this algorithm, we
prefered to enable this feautre.

Figure 2: A sample result of photomaping in OpenCV Figure 3: The real
environment

The next importnat part done in this file is to apply FMI on the sequence of
images captured from the camera. Our FMI algorithm works only on special
format pictures. Since the image server can take picures only in colored
format and rectangular format, we have to convert them to grayscale and square
shape. For doing this task we used GraphicMagick library. So for using this
part the user has to install this package first. After cropping and converting
the picture sequence, we applied FMI algorithm to it. The output of the FMI
algorithm is X, Y, Rotation and Scale. These variables associated with the new
image can create the whole photomap. In order to create the photomap, we put
the new picture in the appropriate place in our GUI. We developed our own GUI
using OpenCV to present the output of the algorithm in online mode. One of the
outputs is shown in figure 2. The real environment that we applied our
algorithm is shown in figure 3. One of the problems of using OpenCV is that it
needs a fixed size for the photomap. Since we are creating the map in time, it
is a growing map and we can’t have a fixed size for it. There are several
solutions for handling it, but everything should be done from scratch. So it
would be easier to use their GUI.

The idea was to send these variables and images to the original GUI . Then we
had to send it on Wireless Simulation Server to be captured by that GUI. We
talked with them and Ravi said they are to change the Image message and it is
not possible to send it to their GUI. The other option is not to use WSS
server and use some parts of their source code and develop our own QTObject in
order to integrate it with their GUI. But lots of things have to be changed,
because their GUI slows down when the number of images is increasing. So it is
not again usable and we have to fix this problem of the GUI first. By default
QT is repainting everything on its canvas. It doesn’t care which part has
recetly been changed. This will force high workload on CPU and memory of the
computer, and slows down the overal progress. The idea is to render and
repaint only the newer part in photomap. Since it was suggested very late, a
few days ago, we haven’t had enough time to do it yet. But we are trying to do
it as soon as possible.

##  3 Conclusion

Having a general map of a disaster situation is very useful in a cooperative
multirobot system. In this report we explained our approach for integrating
FMI algorithm into Jacobotics source code, in order to have a growing map of
the environmet in a disaster situation. We showed that the job was done pretty
well, but there are still some tasks to be done in near future. The most
important task is to take it to a QT GUI from an OpenCV one.

## References

  * [1] USARsim: Unified system for automation and robotics simulator (usarsim) (2006)
  * [2] B. Zitova, J. Flusser, Image registration mehotds: a survey, Elsevier, (2003)
  * [3] A. Jacoff, E. Messina, J. Evans. ”Experiences in deploying test arenas for autonomous mobile robots”, in Proceedings of Performance Metrics for Intelligent Systens, 2001.
  * [4] Nuchter, Andreas, Seyedshams Feyzabadi, Deyuan Qiu, and Stefan May. ”SLAM à la carte-GPGPU for Globally Consistent Scan Matching.” In ECMR, pp. 271-276. 2011.
  * [5] Lucas, Bruce D., and Takeo Kanade. ”An iterative image registration technique with an application to stereo vision.” IJCAI. Vol. 81. 1981.
  * [6] Baumann, Peter, Shams Feyzabadi, and Constantin Jucovschi. ”Putting pixels in place: A storage layout language for scientific data.” Data Mining Workshops (ICDMW), 2010 IEEE International Conference on. IEEE, 2010.
  * [7] Maes, Frederik, et al. ”Multimodality image registration by maximization of mutual information.” Medical Imaging, IEEE Transactions on 16.2 (1997): 187-198.
  * [8] Feyzabadi, S., and J. Schonwalder. ”Identifying TCP congestion control algorithms using active probing.” Passive and Active Measurement Conference (PAM), Poster. 2010.