# Turbulent convection as a significant hidden provider of magnetic helicity
in solar eruptions

Shin Toriumi  Institute of Space and Astronautical Science, Japan Aerospace
Exploration Agency, 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa 252-5210,
Japan<EMAIL_ADDRESS>Hideyuki Hotta  Department of Physics, Graduate
School of Science, Chiba University, 1-33 Yayoi-cho, Inage-ku, Chiba 263-8522,
Japan  Institute for Space-Earth Environmental Research, Nagoya University,
Nagoya 464-8601, Japan  Kanya Kusano  Institute for Space-Earth Environmental
Research, Nagoya University, Nagoya 464-8601, Japan

###### Abstract

Solar flares and coronal mass ejections, the primary space weather
disturbances affecting the entire heliosphere and near-Earth environment,
mainly emanate from sunspot regions harbouring high degrees of magnetic twist.
However, it is not clear how magnetic helicity, the quantity for measuring the
magnetic twist, is supplied to the upper solar atmosphere via the emergence of
magnetic flux from the turbulent convection zone. Here, we report state-of-
the-art numerical simulations of magnetic flux emergence from the deep
convection zone. By controlling the twist of emerging flux, we find that with
the support of convective upflow, the untwisted emerging flux can reach the
solar surface without collapsing, in contrast to previous theoretical
predictions, and eventually create sunspots. Because of the turbulent twisting
of magnetic flux, the produced sunspots exhibit rotation and inject magnetic
helicity into the upper atmosphere, amounting to a substantial fraction of
injected helicity in the twisted cases that is sufficient to produce flare
eruptions. This result indicates that the turbulent convection is responsible
for supplying a non-negligible amount of magnetic helicity and potentially
contributes to solar flares.

## Introduction

Solar flares and coronal mass ejections are the primary sources of space
weather disturbances driving various plasma processes in the entire
interplanetary space, including the near-Earth environment[1, 2, 3]. The
strongest events among these eruptions emanate from the solar active regions,
in which the rotation and shear motions of strongly magnetised sunspots are
often observed[4]. Once a solar flare occurs and a coronal mass ejection is
launched, a helical magnetic flux rope is observed in the interplanetary
space[5, 6], which is the clearest evidence that solar flares are the sudden
release of the excessive magnetic energy accumulated in the helical magnetic
structure in the solar corona through magnetic reconnection and plasma
instability[7, 8, 9].

Magnetic helicity is a measure to quantify the topology of a magnetic field,
such as twists, kinks, and internal linkages:

$\displaystyle H=\int\mbox{\boldmath$A$}\cdot\mbox{\boldmath$B$}\,dV,$ (1)

where $A$ is the vector potential of a magnetic field $B$, i.e.,
$\mbox{\boldmath$B$}=\nabla\times\mbox{\boldmath$A$}$. It is well conserved
even in resistive magnetohydrodynamic (MHD) processes but gauge-invariant only
if the magnetic field is fully contained in a closed volume. Otherwise, the
relative magnetic helicity,[10, 11]

$\displaystyle H_{\rm R}=\int(\mbox{\boldmath$A$}+\mbox{\boldmath$A$}_{\rm
p})\cdot(\mbox{\boldmath$B$}-\mbox{\boldmath$B$}_{\rm p})\,dV.$ (2)

is widely used because it is gauge invariant in any case, and the potential
magnetic field $\mbox{\boldmath$B$}_{\rm p}$ is often adopted as a reference.
The total amount of injected helicity, $H_{\rm R}$, and the helicity injection
rate (helicity flux),

$\displaystyle\frac{dH_{\rm R}}{dt}=2\int_{S}\left[(\mbox{\boldmath$A$}_{\rm
p}\cdot\mbox{\boldmath$B$}_{\rm h})V_{z}-(\mbox{\boldmath$A$}_{\rm
p}\cdot\mbox{\boldmath$V$}_{\rm h})B_{z}\right]\,dS,$ (3)

where the subscripts h and $z$ denote the horizontal and vertical directions,
respectively, are widely used in analyses of flare-productive active regions.
For instance, observations show that active regions with a larger amount of
magnetic helicity tend to produce stronger flares[12, 13, 14] and that the
activity level is enhanced in association with the total injected helicity or
the temporal variation of helicity flux[15, 16, 17, 18, 19]. This is why
magnetic helicity attracts broad attention of the solar and heliophysics
community in the context of flare prediction and forecasting.

Based on a copious amount of observational evidence, it is widely believed
that the injection of magnetic helicity into the corona is mainly due to the
emergence of twisted magnetic flux from the convection zone and the
consequential motions of sunspots, such as rotation and shearing[20, 21, 22,
23, 24, 16, 25, 26, 27]. Therefore, investigating how magnetic helicity is
accumulated in the corona as the magnetic flux emerges and builds up active
regions is an important factor in understanding the occurrence of solar
flares. However, it has almost never been considered how and to what extent
the background convection affects the helicity injection, which is probably
because we cannot probe the solar interior using direct optical observations.
The detection of flux emergence in the convection zone using helioseismology
is a promising technique but still in the development stage.[28, 29, 30]

To solve this problem, we perform numerical simulations in which a twist-free
magnetic flux tube emerges from the turbulent convection zone and examine
whether the helicity injection into the upper atmosphere is negligibly small
or comparable to the observations. Previous theoretical studies suggested that
an untwisted emerging flux cannot reach the photosphere because it experiences
counteracting aerodynamic drag (see[31] and the references therein). Under the
ideal condition with no background convection, even if the non-twisted flux
tube reaches the photosphere, the helicity injection will be zero (as far as
the emerging bipole keeps its geometrical symmetry). In this study, we explore
the occurrence of helicity injection by calculating the case with convection
and, if so, to what extent.

We perform our computations using radiative MHD code R2D2[32] to reproduce the
realistic turbulent thermal convection in the Sun. We use the convection model
that is in a statistically steady state as the initial condition ($t=0$). We
do not consider solar rotation in our model; therefore, the net kinetic
helicity of the background flow field is negligibly small. The simulations
with no solar rotation and thus infinitesimal net kinetic helicity enable
investigation of the magnetic helicity injection that is caused purely by
turbulence.

At $t=0$, a magnetic flux tube with an axial field strength of 12 kG and a
typical radius of 8 Mm is placed at a depth of 22 Mm in the rectangular
computational box that covers the entire solar convection zone, with the box
spanning over $0\leq x\leq 98.3\ {\rm Mm}$, $0\leq y\leq 98.3\ {\rm Mm}$, and
$-201\ {\rm Mm}\leq z\leq 676\ {\rm km}$ (bottom panel of Figure 1). The
bottom boundary ($-201\ {\rm Mm}$) is deeper than most previous convective
flux emergence simulations, which were at most $-30\ {\rm Mm}$[33, 34, 35].
Thus, the present model allows us to investigate the effects of large-scale
convection on the emergence of magnetic flux and the resultant sunspot
formation[36, 37, 38, 39]. In the initial state, a mechanical balance between
the flux tube and surroundings is achieved by lowering the entropy inside the
tube. Therefore, the magnetic flux starts rising in response to the background
velocity field without artificial buoyancy. For the purpose of comparison, we
calculate two additional cases where the flux tubes are weakly and strongly
twisted (right-handed twist), but the background convection field remains the
same. The twist strengths of the weakly- and strongly-twisted tubes are
$q_{\rm cr}/4$ and $q_{\rm cr}/2$, respectively, where $q_{\rm cr}\,(=0.125\
{\rm Mm}^{-1})$ is the critical twist for the kink instability[40], namely,
these tubes are stable against the instability. In all three cases, the total
magnetic flux in the axial direction is on the order of $2\times 10^{22}\ {\rm
Mx}$.

## Results

### General evolution

Figure 1: The vertical magnetic field strength ($B_{z}$) at the photospheric
surface, emergent intensity normalised by the average quiet-Sun intensity
($I/I_{0}$), three-dimensional (3D) volume rendering of the total magnetic
field ($|\mbox{\boldmath$B$}|$) and initial magnetic field lines for the non-
twisted, weakly-twisted and strongly-twisted cases. The white square
emphasises the region where the two sunspots collide and rotate. The square
also corresponds to the field of view of Figure 3. A movie for the vertical
field and intensity is available in the online Supplementary Information.

Figure 1 shows the vertical magnetic field $B_{z}$ in the photospheric
surface, the emergent intensity and the magnetic field strength and field
lines in the 3D space for the three runs, i.e., the non-, weakly- and
strongly-twisted cases. In this study, we define the layer where the optical
depth is unity ($\tau=1$) as the photospheric surface and measure physical
quantities to make direct comparisons with observations. Figure 1 shows that
the flux tube is levitated by a convective upflow located in the centre of the
computational box in all three cases. Even the non-twisted flux tube
successfully reaches the photosphere, in contrast to the previous prediction
that a flux tube needs some twisting to maintain its cohesion against
aerodynamic drag[31].

From $t=20\ {\rm hr}$, the positive and negative magnetic elements distributed
in the photosphere gradually assemble and build up sunspots with positive and
negative polarities. The positive sunspot moves beyond the side periodic
boundary and at around $t=30\ {\rm hr}$, the two sunspots collide with each
other, eventually creating a strongly-packed bipolar sunspot
($\delta$-sunspot), which is known to be highly flare-active [41]. This
situation agrees with the scenario where a single flux system emerges at
multiple locations to build up a colliding bipolar sunspot [42, 43, 36]. The
magnetic energy of the flux tube at the initial state is set to be the same
for the three cases. However, once the sunspots are established in the
photosphere, their decay is more rapid for the weaker twist cases because the
local convection cells can more easily intrude into the sunspots and break
them into pieces. In the twist-free case, the sunspots disappear by around
$t=60\ {\rm hr}$, i.e., the lifetime in the photosphere is approximately 40
hr.

One remarkable feature of the developed sunspots is their continued rotation.
In the strong twist case, the two sunspots rotate in the same clockwise
direction because the flux tube, endowed initially with right-handed torsion,
releases its twist as it appears in the photosphere[44, 45, 36]. Observations
show that many more flaring active regions show rotations of bipolar sunspots
in the same directions than in the opposite directions[46]. Of particular
note, however, is that the sunspots in the no-twist case also exhibit
rotations, with the negative sunspot in a clockwise direction and the positive
sunspot in a counterclockwise direction (the square-framed zone in Figure 1).
Because this flux tube is not given any twist at the initial state, the
observed sunspot rotations are presumed to be driven by the background
turbulence beneath the solar surface.

### Magnetic helicity injection and sunspot rotation

Figure 2: Temporal evolutions of total unsigned magnetic flux ($\Phi$),
relative magnetic helicity ($H_{\rm R}$) and helicity flux ($dH_{\rm R}/dt$)
measured in the photosphere for the three cases.

Figure 2 shows the temporal evolutions of injected magnetic helicity, $H_{\rm
R}$, and helicity flux, $dH_{\rm R}/dt$, as compared with the total unsigned
magnetic flux,

$\displaystyle\Phi=\int|B_{z}|\,dS,$ (4)

all measured in the photosphere for the three cases. The initial total
magnetic flux in the flux tubes in the axial direction is of the order of
$2\times 10^{22}\ {\rm Mx}$. Thus, if a flux tube emerges bodily as a whole to
form a pair of sunspots, the total photospheric flux would be $4\times
10^{22}\ {\rm Mx}$. However, in all cases it is on the order of $3\times
10^{22}\ {\rm Mx}$, indicating that some fractions of the flux tubes remain
below the photosphere. The total fluxes peaked around $t=40\ {\rm hr}$ and
then gradually decreased as the sunspots decayed.

The middle panel of Figure 2 illustrates that positive magnetic helicity is
injected in the cases of strong (black line) and weak (blue line) twists. This
is reasonable considering that these flux tubes were initially given right-
handed twists and that the sunspots displayed clockwise rotations. What is
noteworthy about this plot is that the positive helicity injection also occurs
in the non-twisted case (red line). The accumulated helicity amounts as much
as to $2.9\times 10^{43}\ {\rm Mx}^{2}$, which is about 20% to 50% of those in
the twisted cases. This result reveals that a non-negligible amount of
magnetic helicity can be injected even when a twist-free flux tube emerges in
the convective zone with null net kinetic helicity. It should be noted that
the magnetic helicity normalised by the square of the total photospheric flux,
$H_{\rm R}/\max{(\Phi)}^{2}$, for the three cases is 0.029, 0.058 and 0.094
(peak values), which are comparable to the recent observations of flaring
active regions[47, 48]. Therefore, the present simulations provide the
reasonable reproduction of solar active regions. The bottom panel of Figure 2
shows the evolution of helicity flux (see Equation (3)).

Figure 3: 3D view of the magnetic structures below the two rotating sunspots,
whose field of view is indicated in Figure 1. The reddish volume rendering
shows the total field strength ($|\mbox{\boldmath$B$}|$), while the yellow
lines are the selected magnetic field lines. The 2D slice at the top presents
the photospheric vertical field ($B_{z}$) at the $\tau=1$ surface, with the
grayscale saturating at $\pm 2\ {\rm kG}$. The three panels on the right show
the velocity fields ($V_{z}$ and $\mbox{\boldmath$V$}_{\rm h}$ vector) at
$-5$, $-10$ and $-20\ {\rm Mm}$. Contours (black or green depending on the
background) show the 50% and 80% levels of the maximum total field strength
($|\mbox{\boldmath$B$}|$) of the entire horizontal (i.e. $98.3\ {\rm Mm}\times
98.3\ {\rm Mm}$) plane at each depth, representing the two magnetic pillars.
Cyan arrows highlight the directions of the horizontal flows. All physical
values are obtained by averaging the data sets over a period between $t=35.0\
{\rm hr}$ and 38.0 hr.

Then, what effect causes a finite positive helicity injection in the non-
twisted case? Figure 3 shows the 3D rendering below the sunspot collision area
(corresponding to the square-framed zone in Figure 1) and the plots of
velocity fields at three different depths. We average the data over the period
from $t=35\ {\rm hr}$ to 38 hr, when the helicity flux is peaked (bottom panel
of Figure 2). Below the positive and negative sunspots, a pair of vertical
magnetic fluxes (represented by yellow magnetic field lines) extend downward
into the deep convection zone. The velocity plots reveal that the two magnetic
pillars (indicated by contours) reside in the strong downflowing plumes, to
which surrounding plasmas stream in, leaving local vortices (highlighted by
cyan arrows). These structures are created because some portions of the
initially horizontal flux tube are dragged into the strong downflow plumes and
become vertical magnetic concentrations[35, 37]. At the same time, the plumes
also drive the surrounding plasmas to flow into them, accompanied by the
vortices. The directions of the vortices agree with those of the sunspot
rotations in the photosphere, i.e., clockwise (counterclockwise) for the
negative (positive) sunspot. This indicates that the local vortices streaming
into the downflow plumes spin the magnetic pillars and drive the sunspot
rotations in the photosphere. If the sunspot rotation is because of the
release of magnetic twists, the two magnetic pillars (sunspots in the
photosphere) should rotate in the same direction[44, 45]. However, this is not
the case. Thus, even if the flux tube is initially endowed with no magnetic
twist, it is possible that the surrounding turbulent flow, in which the local
vortices reside, exerts a spinning effect on the flux tube and, accordingly,
the sunspot rotation occurs in the photosphere, leading to the significant
injection of magnetic helicity into the upper atmosphere.

### Flare productivity

Figure 4: (Top) Non-potential magnetic field ($B_{\rm
np}=|\mbox{\boldmath$B$}-\mbox{\boldmath$B$}_{\rm p}|$) for the non-twisted
and strongly-twisted cases. The cyan and orange lines indicate the smoothed
contours of $B_{z}=-2\ {\rm kG}$ and 2 kG, respectively. (Middle) Magnetic
twist flux density ($\tau_{\rm twist}=T_{\rm W}|B_{z}|$). Contours are
identical to those in the top panels. (Bottom: flare phase diagram) Critical
radius of the circular reconnection region required to satisfy the double-arc
instability (DAI) condition, $r_{\rm c}$, and the magnetic energy that can be
released, $E_{\rm r}$, for all points on the polarity inversion lines in the
identified high free-energy regions (HiFERs). The colours denote the different
HiFER patches, with the reddish (blueish) colours for the larger (smaller)
patches. Observations show that X2-class flares or larger are likely if there
are points in the area of $E_{\rm r}>4\times 10^{31}\ {\rm erg}$ and $r_{\rm
c}<1\ {\rm Mm}$[49]. The black curve represents the self-similar scaling of
$r_{\rm c}\propto E_{\rm r}^{1/3}$.

One of the most promising methodologies that have been suggested to predict
the flare eruptions is to analyse the spatial distribution of high free-energy
regions (HiFERs), where the non-potential magnetic field, $B_{\rm
np}=|\mbox{\boldmath$B$}-\mbox{\boldmath$B$}_{\rm p}|$, exceeds the critical
value of 1 kG, and examine the occurrence of the double-arc instability
(DAI[50]) for each HiFER patch near the polarity inversion line[49]. Top and
middle panels of Figure 4 show the distributions of $B_{\rm np}$ and the
magnetic twist flux density $\tau_{\rm twist}$ at the representative timings
for the non-twisted and strongly-twisted cases. The bottom panels are the
flare phase diagram, which is the scatterplot of the minimum critical radius
of the circular reconnection region required to satisfy the DAI condition,
$r_{\rm c}$, versus the magnetic energy that can be released, $E_{\rm r}$, for
the identified HiFERs. Observations show that large X-class flares likely
occur when there is a HiFER with $r_{\rm c}<1\ {\rm Mm}$ and $E_{\rm
r}>4\times 10^{31}\ {\rm erg}$ on the flare phase diagram[49].

In the strongly-twisted case, the HiFER patches ($B_{\rm np}>1\ {\rm kG}$) are
coherent and the positive twist is distributed at the polarity inversion line
between the two major sunspots, which is in agreement with the analysis
results in the previous subsections. In the flare phase diagram, several
locations are very close to the X-flare criticality. Therefore, we can predict
that M-class flares are possible with a slight chance of X flares for the
strong twist case. In the non-twisted case, however, although there are
regions where $B_{\rm np}$ exceeds 2 kG, the HiFERs are more fragmented. The
twists around the polarity inversion lines are a mixture of positive and
negative signs, suggesting that the twists of both signs can be generated,
probably randomly, because of the nature of turbulent convection, when a
twist-less magnetic flux emerges. No HiFERs satisfy the X-flare threshold in
the flare phase diagram, but lower-class (e.g. C-class) flares may still
occur.

## Discussion

We investigate the physical mechanisms that provide magnetic helicity into the
solar corona, which is critically important in understanding the occurrence of
flare eruptions, by comparing the flux emergence simulations with and without
magnetic twists under the existence of turbulent thermal convection. It is
widely believed that the emergence of twisted flux tubes and associated
sunspot rotations supply helicity into the active region corona. However, this
is not trivial because of the lack of optical observations below the
photosphere, which is the primary motivation of this study. We summarise the
key results as follows.

  1. 1.

A non-twisted flux tube can reach the photosphere and develop an active
region. This is contrary to the previous simulation results without thermal
convection during the emergence, a rising untwisted flux tube splits into two
vortex rolls and quickly loses its identity because of the counteracting
aerodynamic drag[31]. With the support of convection, a non-twisted flux tube
may complete its journey to the solar surface.

  2. 2.

The injected magnetic helicity is non-zero even in the case of no magnetic
twist. In fact, the amount of injected helicity was quite large, reaching
about 50% as in the twisted cases. Our analysis reveals that in the convection
zone below the developed sunspots, inflows into the strong downwelling plumes
produce local vortex motions, which spin the vertically standing magnetic flux
tubes that extend along the plumes. As a result, the sunspots in the
photosphere show rotations in the directions determined by the subsurface
vortices, and a considerable amount of helicity injection into the upper
atmosphere occurs.

  3. 3.

The instability analysis suggests that the generated non-potential magnetic
field in the non-twisted case can produce eruptions, albeit weaker than the
twisted cases. This indicates that the “twist” observed in a variety of forms
in the flare-productive active regions of the Sun is not only supplied by the
magnetic “twist” of the emerging flux and the associated sunspot rotations, as
previously believed, but also includes a non-negligible fraction of the
“twist” brought by the background turbulent convection.

In the present simulations, although the net kinetic helicity of the initial
background convection was infinitesimally small, the resultant helicity
injection in the photosphere is comparable to that provided by the emergence
of twisted flux tubes. Considering the nature of the turbulence where we do
not consider the solar rotation, the directions of the rotations of the local
vortices that contributed to the spinning of the flux tubes may probably be
determined by chance. Because the downflow plumes accompanied by the vortices
are a coherent structure owing to the large time scale in the deeper
convection zone, although the directions of the vortices are randomly
determined, the sign of helicity injection may remain the same (positive in
the present runs) over time. This indicates that the sign of helicity
injection in simulations with different background turbulence can be negative,
which will be examined in the future studies.

Asymmetry in sunspot rotation is reported in bipolar active regions[51, 46],
and this could be due to the difference in vortex motions in the convection
zone. Also, a considerable fraction of solar active regions do not obey the
hemispheric helicity sign rule, with the fraction of violators being 20–40%
[52, 53]. This randomness may be the result of the stochastic imposing of
magnetic helicity by the background turbulence to the twist of the emerging
flux that is determined by the solar dynamo[54].

## Methods

In this study, we use the radiative MHD code R2D2[32], which realistically
simulates the thermal convection over the entire solar convection zone, to
model the emergence of magnetic flux tubes and the spontaneous sunspot
formation[36, 37, 38, 39]. This code computes the 3D MHD equations with
realistic radiation transfer and equation of state, and by implementing
RSST[55, 56, 57, 58], it effectively suppresses the high sound speed in the
deep solar interior and relaxes the Courant-Friedrichs-Lewy condition.

We use a rectangular Cartesian box as the computational domain, spanning over
$0\leq x\leq 98.3\ {\rm Mm}$, $0\leq y\leq 98.3\ {\rm Mm}$, and $-201\ {\rm
Mm}\leq z\leq 676\ {\rm km}$, resolved by a $1024\times 1024\times 384$ grid.
The grid spacing for the horizontal directions was $\Delta x=\Delta y=96\ {\rm
km}$ (uniform), while that for the vertical direction was uniform at $\Delta
z=48\ {\rm km}$ from the top boundary to $z=-5.6\ {\rm Mm}$, which linearly
increased to the bottom boundary up to $\Delta z=1486\ {\rm km}$. The
horizontal boundaries assumed the periodic boundary condition, while the
magnetic field at the top boundary was connected to the potential field. The
initial background convection was the same as that in[38], in which a magnetic
flux tube emerges to successfully produce a bipolar sunspot group. The net
kinetic helicity is negligibly small because the calculation box did not
consider the solar rotation. The normalised net kinetic helicity at $t=0$,
when the flux tube was embedded, was;

$\displaystyle\frac{\displaystyle\int_{z<0}\mbox{\boldmath$V$}\cdot(\nabla\times\mbox{\boldmath$V$})\,dV}{\displaystyle\int_{z<0}|\mbox{\boldmath$V$}\cdot(\nabla\times\mbox{\boldmath$V$})|\,dV}=0.00523\%.$
(5)

The magnetic flux tube was embedded at a depth of 22 Mm. Unlike the previous
R2D2 simulations, which used the force-free magnetic flux tubes, we adopt
commonly used Gaussian-type flux tubes, where the longitudinal and azimuthal
components of the magnetic field are given by

$\displaystyle B_{x}(r)=B_{\rm tb}\exp{\left(-\frac{r^{2}}{R_{\rm
tb}^{2}}\right)}$ (6)

and

$\displaystyle B_{\phi}(r)=qrB_{x}(r),$ (7)

where $B_{\rm tb}$, $r$, $R_{\rm tb}$ and $q$ are the axial field strength,
radial distance from the tube axis, typical radius and twist intensity,
respectively. We calculate three cases with varying $q$ (Table 1). Here,
$B_{\rm tb}$ was also varied so that the magnetic energies of the flux tubes,
$E_{\rm mag}(=\int B^{2}/(8\pi)\,dV)$, are the same in all three cases. The
total axial magnetic flux $\Phi_{x}\,(=\int B_{x}\,dS)$ is of the order of
$2\times 10^{22}\ {\rm Mx}$. To achieve the mechanical balance with no
magnetic buoyancy, we adjust the internal pressure and density by reducing the
entropy by $\frac{B^{2}}{8\pi}/\left(\frac{\partial p}{\partial
s}\right)_{\rho}$. Therefore, the flux tubes started emergence only because of
advection exerted by the background turbulence. In all three cases, we set $q$
to be 0 or below the critical value for the kink instability $q_{\rm
cr}\,(=1/R_{\rm tb})$[40], where $q>0$ indicates that the tube has a right-
handed twist.

Case | $R_{\rm tb}$ | $q$ | $B_{\rm tb}$ | $E_{\rm mag}$ | $\Phi_{x}$
---|---|---|---|---|---
Non-twisted | 8 Mm | 0 | 12.2 kG | $5.9\times 10^{34}\ {\rm erg}$ | $2.4\times 10^{22}\ {\rm Mx}$
Weakly-twisted | 8 Mm | $q_{\rm cr}/4$ | 12.1 kG | $5.9\times 10^{34}\ {\rm erg}$ | $2.4\times 10^{22}\ {\rm Mx}$
Strongly-twisted | 8 Mm | $q_{\rm cr}/2$ | 11.5 kG | $5.9\times 10^{34}\ {\rm erg}$ | $2.3\times 10^{22}\ {\rm Mx}$
Table 1: Summary of the simulation cases. The critical twist for the kink
instability corresponds to $q_{\rm cr}=1/R_{\rm tb}=0.125\ {\rm Mm}^{-1}$.

We measure several physical quantities including the relative magnetic
helicity ($H_{\rm R}$), helicity flux ($dH_{\rm R}/dt$) and total unsigned
magnetic flux ($\Phi$) at the photospheric surface, here defined as where the
optical depth is unity ($\tau=1$). For measuring the helicity, there is a
degree of freedom in choosing the vector potential $\mbox{\boldmath$A$}_{\rm
p}$ and, in this study, we select vector potentials satisfying
$\mbox{\boldmath$A$}_{\rm p}\cdot\hat{\mbox{\boldmath$z$}}=0$ and calculate it
using the method by[59],

$\displaystyle\mbox{\boldmath$A$}_{\rm
p}=\frac{1}{2\pi}\hat{\mbox{\boldmath$z$}}\times\int_{S^{\prime}}B_{z}(\mbox{\boldmath$x$}^{\prime})\frac{\mbox{\boldmath$r$}}{r^{2}}\,dS^{\prime},$
(8)

where $\mbox{\boldmath$r$}=\mbox{\boldmath$x$}-\mbox{\boldmath$x$}^{\prime}$.
In this process, the grid number was reduced from $1024\times 1024$ to
$256\times 256$ to accelerate the computation speed. The helicity flux was
then calculated using Equation (3), and the injected magnetic helicity over
the course of flux emergence (Equation (2)) was obtained by integrating the
helicity flux over time.

To investigate the flare productivity of the generated active regions, the
instability analysis is performed based on the DAI theory, which states that a
flux rope becomes unstable if it is sufficiently twisted and has a sufficient
amount of magnetic flux against the overlying confinement field[50]. This
instability is characterised by the parameter

$\displaystyle\kappa=\left|\frac{\displaystyle\int_{\rm rec}T_{\rm
W}\,d\Phi}{\displaystyle\Phi_{\rm over}}\right|,$ (9)

or equivalently,

$\displaystyle\kappa=\left|\frac{\displaystyle\int_{S_{\rm rec}}\tau_{\rm
twist}\,dS}{\displaystyle\Phi_{\rm over}}\right|,$ (10)

where $T_{\rm W}$ is the amount of twist integrated over each magnetic field
line in a flux rope, $S_{\rm rec}$ is the footpoint area of the magnetic field
lines that reconnect to form the flux rope, $\Phi_{\rm over}$ is the magnetic
flux of the overlying field, and $\tau_{\rm twist}=T_{\rm W}|B_{z}|$. In DAI,
$\kappa$ is usually larger for larger $S_{\rm rec}$; therefore, there is a
critical value $S_{\rm c}$ at which the instability occurs:
$\kappa>\kappa_{0}\sim 0.1$. It is shown that the ratio of the magnetic
helicity of current-carrying magnetic field ($|H_{\rm j}|$) to the total
relative helicity ($|H_{\rm R}|$) well discriminates whether a flare event
becomes eruptive or not[60]. While the $\kappa$ parameter is the critical
parameter which determines the onset of DAI (i.e., whether a flare occur), the
helicity ratio $|H_{\rm j}|/|H_{\rm R}|$ may be capable of distinguishing the
eruptivity when the flare occurs.

In the $\kappa$-scheme[49], the coronal magnetic field is first calculated by
the non-linear force-free field extrapolation based on the vector magnetic
field of the photosphere. For the computational resources, the grid number was
reduced from $1024\times 1024$ to $512\times 256$ in this process. Then,
HiFERs are identified as the regions where the non-potential field $B_{\rm
np}(=|\mbox{\boldmath$B$}-\mbox{\boldmath$B$}_{\rm p}|)$ exceeds the threshold
value of 1 kG, which is based on the observational result that $B_{\rm np}=1$
kG sufficiently encompasses the the distribution of non-potential magnetic
fields that drive large flares[49]. For each HiFER, the critical area $S_{\rm
c}$ is measured as the minimum circular area that satisfies the DAI condition
($\kappa>\kappa_{0}$), and $r_{\rm c}$ is obtained as the radius of $S_{\rm
c}$, i.e., $S_{\rm c}=\pi r_{\rm c}^{2}$. The releasable energy for each HiFER
is estimated as

$\displaystyle E_{\rm r}=\frac{S_{\rm r}^{1/2}}{8\pi}\int_{S_{\rm r}}B_{\rm
np}^{2}\,dS,$ (11)

where $S_{\rm r}$ is the area of the footpoint of the magnetic flux that pass
over the circular area $S_{\rm c}$.

## Data Availability

The data are available from the corresponding author upon reasonable request.

## References

  * [1] Shibata, K. & Magara, T.  Solar Flares: Magnetohydrodynamic Processes.  _Living Reviews in Solar Physics_ 8, 6, DOI: 10.12942/lrsp-2011-6 (2011).
  * [2] Chen, P. F.  Coronal Mass Ejections: Models and Their Observational Basis.  _Living Reviews in Solar Physics_ 8, 1, DOI: 10.12942/lrsp-2011-1 (2011).
  * [3] Fletcher, L. _et al._ An Observational Overview of Solar Flares.  _Space Science Reviews_ 159, 19–106, DOI: 10.1007/s11214-010-9701-8 (2011).
  * [4] Toriumi, S. & Wang, H.  Flare-productive active regions.  _Living Reviews in Solar Physics_ 16, 3, DOI: 10.1007/s41116-019-0019-7 (2019).
  * [5] Zurbuchen, T. H. & Richardson, I. G.  In-Situ Solar Wind and Magnetic Field Signatures of Interplanetary Coronal Mass Ejections.  _Space Science Reviews_ 123, 31–43, DOI: 10.1007/s11214-006-9010-4 (2006).
  * [6] Marubashi, K. & Lepping, R. P.  Long-duration magnetic clouds: a comparison of analyses using torus- and cylinder-shaped flux rope models.  _Annales Geophysicae_ 25, 2453–2477, DOI: 10.5194/angeo-25-2453-2007 (2007).
  * [7] Priest, E. R. & Forbes, T. G.  The magnetic nature of solar flares.  _The Astronomy and Astrophysics Review_ 10, 313–377, DOI: 10.1007/s001590100013 (2002).
  * [8] Gibson, S. E. & Fan, Y.  Coronal prominence structure and dynamics: A magnetic flux rope interpretation.  _Journal of Geophysical Research (Space Physics)_ 111, A12103, DOI: 10.1029/2006JA011871 (2006).
  * [9] Kusano, K. _et al._ Magnetic Field Structures Triggering Solar Flares and Coronal Mass Ejections.  _The Astrophysical Journal_ 760, 31, DOI: 10.1088/0004-637X/760/1/31 (2012).
  * [10] Berger, M. A. & Field, G. B.  The topological properties of magnetic helicity.  _Journal of Fluid Mechanics_ 147, 133–148, DOI: 10.1017/S0022112084002019 (1984).
  * [11] Finn, J. M. & Antonsen, J., Thomas M.  Magnetic helicity: What is it and what is it good for?  _Comments on Plasma Physics and Controlled Fusion_ 9, 111–126 (1985).
  * [12] Nindos, A. & Andrews, M. D.  The Association of Big Flares and Coronal Mass Ejections: What Is the Role of Magnetic Helicity?  _The Astrophysical Journal_ 616, L175–L178, DOI: 10.1086/426861 (2004).
  * [13] LaBonte, B. J., Georgoulis, M. K. & Rust, D. M.  Survey of Magnetic Helicity Injection in Regions Producing X-Class Flares.  _The Astrophysical Journal_ 671, 955–963, DOI: 10.1086/522682 (2007).
  * [14] Tziotziou, K., Georgoulis, M. K. & Raouafi, N.-E.  The Magnetic Energy-Helicity Diagram of Solar Active Regions.  _The Astrophysical Journal Letters_ 759, L4, DOI: 10.1088/2041-8205/759/1/L4 (2012).
  * [15] Moon, Y. J. _et al._ Flare Activity and Magnetic Helicity Injection by Photospheric Horizontal Motions.  _The Astrophysical Journal_ 574, 1066–1073, DOI: 10.1086/340975 (2002).
  * [16] Kusano, K., Maeshiro, T., Yokoyama, T. & Sakurai, T.  Measurement of Magnetic Helicity Injection and Free Energy Loading into the Solar Corona.  _The Astrophysical Journal_ 577, 501–512, DOI: 10.1086/342171 (2002).
  * [17] Park, S.-H. _et al._ The Variation of Relative Magnetic Helicity around Major Flares.  _The Astrophysical Journal_ 686, 1397–1403, DOI: 10.1086/591117 (2008).
  * [18] Park, S.-H. _et al._ The Occurrence and Speed of CMEs Related to Two Characteristic Evolution Patterns of Helicity Injection in Their Solar Source Regions.  _The Astrophysical Journal_ 750, 48, DOI: 10.1088/0004-637X/750/1/48 (2012).
  * [19] Jing, J. _et al._ Evolution of Relative Magnetic Helicity and Current Helicity in NOAA Active Region 11158.  _The Astrophysical Journal Letters_ 752, L9, DOI: 10.1088/2041-8205/752/1/L9 (2012).
  * [20] Zirin, H. & Liggett, M. A.  Delta spots and great flares.  _Solar Physics_ 113, 267–283, DOI: 10.1007/BF00147707 (1987).
  * [21] Tanaka, K.  Studies on a very flare-active $\delta$ group: Peculiar $\delta$ spot evolution and inferred subsurface magnetic rope structure.  _Solar Physics_ 136, 133–149, DOI: 10.1007/BF00151700 (1991).
  * [22] Leka, K. D., Canfield, R. C., McClymont, A. N. & van Driel-Gesztelyi, L.  Evidence for Current-carrying Emerging Flux.  _Astrophysical Journal_ 462, 547, DOI: 10.1086/177171 (1996).
  * [23] López Fuentes, M. C., Demoulin, P., Mandrini, C. H. & van Driel-Gesztelyi, L.  The Counterkink Rotation of a Non-Hale Active Region.  _The Astrophysical Journal_ 544, 540–549, DOI: 10.1086/317180 (2000).
  * [24] Chae, J.  Observational Determination of the Rate of Magnetic Helicity Transport through the Solar Surface via the Horizontal Motion of Field Line Footpoints.  _The Astrophysical Journal Letters_ 560, L95–L98, DOI: 10.1086/324173 (2001).
  * [25] Brown, D. S. _et al._ Observations of Rotating Sunspots from TRACE.  _Solar Physics_ 216, 79–108, DOI: 10.1023/A:1026138413791 (2003).
  * [26] Toriumi, S., Schrijver, C. J., Harra, L. K., Hudson, H. & Nagashima, K.  Magnetic Properties of Solar Active Regions That Govern Large Solar Flares and Eruptions.  _The Astrophysical Journal_ 834, 56, DOI: 10.3847/1538-4357/834/1/56 (2017).
  * [27] Toriumi, S. & Park, S.-H.  Solar Flares and Magnetic Helicity.  _arXiv e-prints_ arXiv:2204.06010, DOI: 10.48550/arXiv.2204.06010 (2022).
  * [28] Ilonidis, S., Zhao, J. & Kosovichev, A.  Detection of Emerging Sunspot Regions in the Solar Interior.  _Science_ 333, 993, DOI: 10.1126/science.1206253 (2011).
  * [29] Birch, A. C., Braun, D. C., Leka, K. D., Barnes, G. & Javornik, B.  Helioseismology of Pre-emerging Active Regions. II. Average Emergence Properties.  _The Astrophysical Journal_ 762, 131, DOI: 10.1088/0004-637X/762/2/131 (2013).
  * [30] Toriumi, S., Ilonidis, S., Sekii, T. & Yokoyama, T.  Probing the Shallow Convection Zone: Rising Motion of Subsurface Magnetic Fields in the Solar Active Region.  _The Astrophysical Journal Letters_ 770, L11, DOI: 10.1088/2041-8205/770/1/L11 (2013).
  * [31] Fan, Y.  Magnetic fields in the solar convection zone.  _Living Reviews in Solar Physics_ 18, 5, DOI: 10.1007/s41116-021-00031-2 (2021).
  * [32] Hotta, H., Iijima, H. & Kusano, K.  Weak influence of near-surface layer on solar deep convection zone revealed by comprehensive simulation from base to surface.  _Science Advances_ 5, 2307, DOI: 10.1126/sciadv.aau2307 (2019).
  * [33] Cheung, M. C. M., Rempel, M., Title, A. M. & Schüssler, M.  Simulation of the Formation of a Solar Active Region.  _The Astrophysical Journal_ 720, 233–244, DOI: 10.1088/0004-637X/720/1/233 (2010).
  * [34] Rempel, M. & Cheung, M. C. M.  Numerical Simulations of Active Region Scale Flux Emergence: From Spot Formation to Decay.  _The Astrophysical Journal_ 785, 90, DOI: 10.1088/0004-637X/785/2/90 (2014).
  * [35] Chen, F., Rempel, M. & Fan, Y.  Emergence of Magnetic Flux Generated in a Solar Convective Dynamo. I. The Formation of Sunspots and Active Regions, and The Origin of Their Asymmetries.  _The Astrophysical Journal_ 846, 149, DOI: 10.3847/1538-4357/aa85a0 (2017).
  * [36] Toriumi, S. & Hotta, H.  Spontaneous Generation of $\delta$-sunspots in Convective Magnetohydrodynamic Simulation of Magnetic Flux Emergence.  _The Astrophysical Journal Letters_ 886, L21, DOI: 10.3847/2041-8213/ab55e7 (2019).
  * [37] Hotta, H. & Iijima, H.  On rising magnetic flux tube and formation of sunspots in a deep domain.  _Monthly Notices of the Royal Astronomical Society_ 494, 2523–2537, DOI: 10.1093/mnras/staa844 (2020).
  * [38] Hotta, H. & Toriumi, S.  Formation of superstrong horizontal magnetic field in delta-type sunspot in radiation magnetohydrodynamic simulations.  _Monthly Notices of the Royal Astronomical Society_ 498, 2925–2935, DOI: 10.1093/mnras/staa2529 (2020).
  * [39] Kaneko, T., Hotta, H., Toriumi, S. & Kusano, K.  Impact of subsurface convective flows on the formation of sunspot magnetic field and energy build-up.  _Monthly Notices of the Royal Astronomical Society_ 517, 2775–2786, DOI: 10.1093/mnras/stac2635 (2022).
  * [40] Linton, M. G., Longcope, D. W. & Fisher, G. H.  The Helical Kink Instability of Isolated, Twisted Magnetic Flux Tubes.  _Astrophysical Journal_ 469, 954, DOI: 10.1086/177842 (1996).
  * [41] Sammis, I., Tang, F. & Zirin, H.  The Dependence of Large Flare Occurrence on the Magnetic Structure of Sunspots.  _The Astrophysical Journal_ 540, 583–587, DOI: 10.1086/309303 (2000).
  * [42] Toriumi, S., Iida, Y., Kusano, K., Bamba, Y. & Imada, S.  Formation of a Flare-Productive Active Region: Observation and Numerical Simulation of NOAA AR 11158.  _Solar Physics_ 289, 3351–3369, DOI: 10.1007/s11207-014-0502-1 (2014).
  * [43] Toriumi, S. & Takasao, S.  Numerical Simulations of Flare-productive Active Regions: $\delta$-sunspots, Sheared Polarity Inversion Lines, Energy Storage, and Predictions.  _The Astrophysical Journal_ 850, 39, DOI: 10.3847/1538-4357/aa95c2 (2017).
  * [44] Fan, Y.  The Emergence of a Twisted Flux Tube into the Solar Atmosphere: Sunspot Rotations and the Formation of a Coronal Flux Rope.  _The Astrophysical Journal_ 697, 1529–1542, DOI: 10.1088/0004-637X/697/2/1529 (2009).
  * [45] Sturrock, Z., Hood, A. W., Archontis, V. & McNeill, C. M.  Sunspot rotation. I. A consequence of flux emergence.  _Astronomy & Astrophysics_ 582, A76, DOI: 10.1051/0004-6361/201526521 (2015).
  * [46] Yan, X. L., Qu, Z. Q. & Kong, D. F.  Relationship between rotating sunspots and flare productivity.  _Monthly Notices of the Royal Astronomical Society_ 391, 1887–1892, DOI: 10.1111/j.1365-2966.2008.14002.x (2008).
  * [47] Moraitis, K., Sun, X., Pariat, É. & Linan, L.  Magnetic helicity and eruptivity in active region 12673.  _Astronomy & Astrophysics_ 628, A50, DOI: 10.1051/0004-6361/201935870 (2019).
  * [48] Thalmann, J. K. _et al._ Magnetic Helicity Budget of Solar Active Regions Prolific of Eruptive and Confined Flares.  _The Astrophysical Journal_ 887, 64, DOI: 10.3847/1538-4357/ab4e15 (2019).
  * [49] Kusano, K., Iju, T., Bamba, Y. & Inoue, S.  A physics-based method that can predict imminent large solar flares.  _Science_ 369, 587–591, DOI: 10.1126/science.aaz2511 (2020).
  * [50] Ishiguro, N. & Kusano, K.  Double Arc Instability in the Solar Corona.  _The Astrophysical Journal_ 843, 101, DOI: 10.3847/1538-4357/aa799b (2017).
  * [51] Yan, X. L. & Qu, Z. Q.  Rapid rotation of a sunspot associated with flares.  _Astronomy and Astrophysics_ 468, 1083–1088, DOI: 10.1051/0004-6361:20077064 (2007).
  * [52] Pevtsov, A. A., Berger, M. A., Nindos, A., Norton, A. A. & van Driel-Gesztelyi, L.  Magnetic Helicity, Tilt, and Twist.  _Space Science Reviews_ 186, 285–324, DOI: 10.1007/s11214-014-0082-2 (2014).
  * [53] Park, S.-H., Leka, K. D. & Kusano, K.  Magnetic Helicity Flux across Solar Active Region Photospheres. I. Hemispheric Sign Preference in Solar Cycle 24.  _The Astrophysical Journal_ 904, 6, DOI: 10.3847/1538-4357/abbb93 (2020).
  * [54] Hotta, H., Rempel, M. & Yokoyama, T.  Large-scale magnetic fields at high Reynolds numbers in magnetohydrodynamic simulations.  _Science_ 351, 1427–1430, DOI: 10.1126/science.aad1893 (2016).
  * [55] Rempel, M.  Solar Differential Rotation and Meridional Flow: The Role of a Subadiabatic Tachocline for the Taylor-Proudman Balance.  _The Astrophysical Journal_ 622, 1320–1332, DOI: 10.1086/428282 (2005).
  * [56] Hotta, H., Rempel, M., Yokoyama, T., Iida, Y. & Fan, Y.  Numerical calculation of convection with reduced speed of sound technique.  _Astronomy & Astrophysics_ 539, A30, DOI: 10.1051/0004-6361/201118268 (2012).
  * [57] Hotta, H., Rempel, M. & Yokoyama, T.  High-resolution Calculation of the Solar Global Convection with the Reduced Speed of Sound Technique. II. Near Surface Shear Layer with the Rotation.  _The Astrophysical Journal_ 798, 51, DOI: 10.1088/0004-637X/798/1/51 (2015).
  * [58] Iijima, H., Hotta, H. & Imada, S.  Semiconservative reduced speed of sound technique for low Mach number flows with large density variations.  _Astronomy & Astrophysics_ 622, A157, DOI: 10.1051/0004-6361/201834031 (2019).
  * [59] Pariat, E., Démoulin, P. & Berger, M. A.  Photospheric flux density of magnetic helicity.  _Astronomy and Astrophysics_ 439, 1191–1203, DOI: 10.1051/0004-6361:20052663 (2005).
  * [60] Pariat, E. _et al._ Relative magnetic helicity as a diagnostic of solar eruptivity.  _Astronomy & Astrophysics_ 601, A125, DOI: 10.1051/0004-6361/201630043 (2017).

## Acknowledgements

The initial convection data were obtained using the Supercomputer Fugaku
provided by the RIKEN Center for Computational Science. Numerical computations
were in part carried out on Cray XC50 at Center for Computational
Astrophysics, National Astronomical Observatory of Japan. This work was
supported by JSPS KAKENHI Grant Nos. JP20KK0072 (PI: S. Toriumi), JP20K14510
(PI: H. Hotta), JP21H01124 (PI: T. Yokoyama), JP21H04492 (PI: K. Kusano),
JP21H04497 (PI: H. Miyahara), and by MEXT as a Program for Promoting
Researches on the Supercomputer Fugaku.

## Author contributions statement

S.T. conducted the simulations, analysed the results, and wrote the
manuscript. H.H. developed the simulation code and prepared the initial
convection data. K.K. conducted the $\kappa$-scheme analysis. All authors
reviewed the manuscript.

## Additional information

Correspondence and requests for materials should be addressed to S.T.

## Competing interests

The authors declare no competing interests.


# A discrete analogue of Terrell’s characterization of rectangular
distributions

Nickos Papadatos111 e-mail<EMAIL_ADDRESS>url:
users.uoa.gr/$\sim$npapadat/

(  | National and Kapodistrian University of Athens, Department of Mathematics,
---|---
| Section of Statistics and Operations Research, Panepistemiopolis, 157 84
Athens, Greece.
)

###### Abstract

George R. Terrell (1983, Ann. Probab., vol. 11(3), pp. 823–826) showed that
the Pearson coefficient of correlation of an ordered pair from a random sample
of size two is at most one-half, and the equality is attained only for
rectangular (uniform over some interval) distributions.

In the present note it is proved that the same is true for the discrete case,
in the sense that the correlation coefficient attains its maximal value only
for discrete rectangular (uniform over some finite lattice) distributions.

MSC: Primary 60E15; 62E10; Secondary 62G30.
Key words and phrases: discrete rectangular distribution; order statistics;
Hahn polynomials; Pearson coefficient of correlation.

##  1 Introduction and main result

For independent, identically distributed, random variables $X_{1}$, $X_{2}$,
from a probability distribution function $F$, the corresponding order
statistics will be denoted by $X_{1:2}\leq X_{2:2}$, that is,
$X_{1:2}=\min\\{X_{1},X_{2}\\}$, $X_{2:2}=\max\\{X_{1},X_{2}\\}$, and the
Pearson coefficient of correlation by

$\rho_{12}:=\frac{\mbox{\rm\hskip 0.86108ptCov\hskip
0.86108pt}(X_{1:2},X_{2:2})}{\sqrt{\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}X_{1:2}}\sqrt{\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}X_{2:2}}}.$

For a random pair $(X,Y)$, the correlation coefficient, $\rho(X,Y)$, is well
defined (and belongs to the interval $[-1,1]$) if and only if both $X$, $Y$
are non-degenerate with finite second moment. Thus, for $\rho_{12}$ to be well
defined it is necessary and sufficient that $F$ is non-degenerate and posses
finite second moment (so that $0<\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}X_{j:2}<\infty$, $j=1,2$).

It is known for a long time that $\rho_{12}>0$; this follows immediately if we
take expectations to $X_{1:2}X_{2:2}=X_{1}X_{2}$ and
$X_{1:2}+X_{2:2}=X_{1}+X_{2}$, yielding $\mbox{\rm\hskip 0.86108ptCov\hskip
0.86108pt}(X_{1:2},X_{2:2})=\left(\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}X_{2:2}-\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}X_{1}\right)^{2}>0$, since $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}X_{1}<\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}X_{2:2}$ for any
non-degenerate $F$. At the time of 80’s, David Scott and Robert Bartoszyński,
in connection with a problem in cell division, proposed to George R. Terrell
the conjecture that $\rho_{12}$ is never greater then one-half. This was
proved true:

######  Theorem 1.1.

(Terrell, 1983). If $F$ is a continuous distribution with finite variance then
$\rho_{12}\leq 1/2$, with equality if and only if $F$ is a rectangular
distribution.

Terrell’s result is based on the development of a function in a series of
Legendre polynomials, it is quite complicated, and imposes the unnecessary
restriction that $F$ is continuous. Székely and Móri (1985) interpreted
Terrell’s result as a maximal correlation problem – see Gebelein (1941), Rényi
(1959) – that is,

$R(X_{1:2},X_{2:2}):=\sup_{g_{1},g_{2}}\rho(g_{1}(X_{1:2}),g_{2}(X_{2:2}))$

where the supremum is taken over non-constant functions $g_{1}\in
L^{2}(X_{1:2})$, $g_{2}\in L^{2}(X_{2:2})$. In this way, Székely and Móri
improved Terrell’s result in four directions. First, they removed the
restriction that $F$ is continuous; second, they simplified Terrell’s proof;
third, they showed that the correlation coefficient of $g_{1}(X_{1:2})$ and
$g_{2}(X_{2:2})$ is less than $1/2$, for any distribution and any functions
$g_{1}(X_{1:2})$ and $g_{2}(X_{2:2})$; fourth, and most important, they
extended these results to any sample size $n\geq 2$, obtaining the inequality
(denote by $X_{1:n}\leq\cdots\leq X_{n:n}$ the order statistics from a random
sample of size $n$ from $F$)

$\rho(X_{i:n},X_{j:n})\leq\sqrt{\frac{i(n+1-j)}{j(n+1-i)}},\ \ \ 1\leq i<j\leq
n,\ \ n\geq 2,$ (1.1)

which is valid for all non-degenerate distribution functions $F$ with
$\mbox{\rm\hskip 0.86108ptVar\hskip 0.86108pt}X_{i:n}+\mbox{\rm\hskip
0.86108ptVar\hskip 0.86108pt}X_{j:n}<\infty$. The equality in (1.1), for a
single value of $(i,j,n)$, characterizes the rectangular distributions.

Subsequently, Nevzorov (1992) obtained a similar inequality for upper records,
in which the equality characterizes the location-scale family of the standard
exponential random variable. Later on, López-Blázquez and Castaño-Martínez
(2006) proved that the TSM inequality (1.1) is also valid when the order
statistics are based on a without-replacement sample of size $n$, taken from a
finite population ${\cal P}=\\{x_{1},\ldots,x_{N}\\}$ with $N>n$ distinct
elements.

It is shown in Papadatos and Xifara (2013) that all presenting results are
based, essentially, in the following polynomial regression property (PRP):

$\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}(X^{k}|Y)=A_{k}Y^{k}+P_{k-1}(Y),\
\ \ \mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}(Y^{k}|X)=B_{k}X^{k}+Q_{k-1}(X),\ \ k=1,2,\ldots,$ (1.2)

where $P_{k-1}$ and $Q_{k-1}$ are polynomials of degree at most $k-1$. If a
random pair $(X,Y)$ satisfies (1.2) then, under mild conditions, its maximal
correlation, $R(X,Y)$, equals to $\sup_{k\geq 1}\sqrt{A_{k}B_{k}}$, so one
simply has to calculate the principal coefficients $A_{k}$ and $B_{k}$ in
(1.2), and choose the value of $k=k_{0}$ that maximizes the products
$A_{k}B_{k}$. If this value of $k_{0}$ is unique, then the equality in the
inequality

$\rho(g_{1}(X),g_{2}(Y))\leq R(X,Y)$

is attained if and only if $\mbox{\rm\hskip 0.86108ptIP\hskip
0.86108pt}[g_{1}(X)=c_{1}+\lambda_{1}\phi_{k_{0}}(X)]=1$, $\mbox{\rm\hskip
0.86108ptIP\hskip 0.86108pt}[g_{2}(Y)=c_{2}+\lambda_{2}\psi_{k_{0}}(Y)]=1$,
and $\lambda_{1}\lambda_{2}\mbox{sign}(A_{k_{0}})>0$, where
$\\{\phi_{k}\\}_{k=0}^{\infty}$ is the complete, orthonormal polynomial system
in $L^{2}(X)$ (with the convention that each $\phi_{k}$ has positive principal
coefficient) and $\\{\psi_{k}\\}_{k=0}^{\infty}$ the corresponding system in
$L^{2}(Y)$. Let us denote by $U_{1:n}<\cdots<U_{n:n}$ the order statistics
either from uniform in the interval $(0,1)$, or from the discrete uniform in
$\\{1,\ldots,N\\}$ (the latter in the without-replacement case). It is a
simple exercise to verify that the pair $(U_{i:n},U_{j:n})$ posses the PRP
(1.2), and it is easy to calculate $A_{k}$ and $B_{k}$ (in the without-
replacement case this is slightly more complicated); also, it is plain to
check the PRP for upper records $(W_{n},W_{m})$ ($n<m$) from the standard
exponential, and to obtain the constants $A_{k}$ and $B_{k}$ in (1.2). Then,
the PRP method yields all the presenting results at once and, in the discrete
case, it provides the exact characterization of those populations ${\cal P}$
that attain the upper bound in the TSM inequality (1.1), even if ties are
allowed; see Section 3 in Papadatos and Xifara (2013). In all of the forgoing
results, the sequence $A_{k}B_{k}$ is uniquely maximized by its first term,
hence, $R=\sqrt{A_{1}B_{1}}$, and the equality is attained by linear functions
$g_{1}$, $g_{2}$; equivalently, $R(X,Y)=|\rho(X,Y)|$. While the equality
$|\rho|=R$ often appears to problems regarding maximal correlation under PRP,
it is not always true; see Papadatos (2014).

In the present note we shall prove the following discrete analogue of Theorem
1.1. First, we provide a useful definition, cf. Balakrishnan et al. (2003).

######  Definition 1.1.

(Uniform distribution on ${\cal P}=\\{x_{1},\ldots,x_{N}\\}$ in the possible
presence of ties). Let $k_{1}$ of the $x_{i}$’s be equal to $y_{1}$, $k_{2}$
of the $x_{i}$’s be equal to $y_{2}$, and so on, where $k_{1}+\cdots+k_{m}=N$
and, with no loss of generality, $y_{1}<\cdots<y_{m}$. We say that $X$ is
uniform on ${\cal P}=\\{x_{1},\ldots,x_{N}\\}$ if $\mbox{\rm\hskip
0.86108ptIP\hskip 0.86108pt}(X=y_{j})=k_{j}/N$, $j=1,\ldots,m$.

The main result is the following.

######  Theorem 1.2.

Let $X_{1}$, $X_{2}$ be independent random variables with uniform distribution
on ${\cal P}=\\{x_{1},\ldots,x_{N}\\}$, where $N\geq 2$ and
$x_{1}\leq\cdots\leq x_{N}$ with $x_{1}<x_{N}$. Then, under the notation of
Theorem 1.1,

$\rho_{12}\leq\frac{1-N^{-2}}{2+N^{-2}},$ (1.3)

and the equality is attained if and only if ${\cal P}$ is a discrete lattice,
that is, $x_{i+1}-x_{i}=\lambda>0$ (constant), $i=1,\ldots,N-1$.

A straightforward computation shows that the upper bound in (1.3) equals to
$\rho_{12}$ when ${\cal P}={\cal P}_{0}:=\\{1,\ldots,N\\}$. If $U_{1:2}\leq
U_{2:2}$ are the order statistics from ${\cal P}_{0}$, it can be checked that
the random pair $(U_{1:2},U_{2:2})$ does not posses the PRP (1.2); in
particular, $U_{1:2}$ does not have linear regression on $U_{2:2}$. Indeed,
for $x,y\in{\cal P}_{0}$ with $x\leq y$, $\mbox{\rm\hskip 0.86108ptIP\hskip
0.86108pt}(U_{1:2}=x|U_{2:2}=y)=2/(2y-1)$ if $x<y$, and $\mbox{\rm\hskip
0.86108ptIP\hskip 0.86108pt}(U_{1:2}=y|U_{2:2}=y)=1/(2y-1)$. Hence, one finds
$\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}(U_{1:2}|U_{2:2}=y)=y^{2}/(2y-1)$, and this regression (with $y$
restricted to ${\cal P}_{0}$) is, clearly, nonlinear, unless $N=2$. Therefore,
(1.2) fails, and this verifies the essentially different nature of the present
problem, compared to the preceding ones. Notice that the case $N=2$ is
trivial, since $\rho_{12}\equiv 1/3$ for every choice of $x_{1}$, $x_{2}$,
with $x_{1}<x_{2}$ (the coefficient of correlation is location-scale
invariant).

##  2 Proof of Theorem 1.2

We shall apply the ordinary method of Terrell (1983) to the case where $U$
follows a discrete uniform in ${\cal P}_{0}=\\{1,\ldots,N\\}$. In this case,
the orthonormal polynomial system is provided by Hahn polynomials,
$\\{\psi_{k}\\}_{k=0}^{N-1}$, where

$\psi_{k}(x)=N^{1/2}{N+k\choose 2k+1}^{-1/2}{2k\choose
k}^{-1/2}\sum_{j=0}^{k}(-1)^{k-j}{k+j\choose j}{N-1-j\choose k-j}{x-1\choose
j}.$

The following two properties will be used in the sequel; cf. López-Blázquez
and Castaño-Martínez (2006). First,

$\psi_{k}(x)=A_{k}x^{k}+B_{k}x^{k-1}+\cdots,$

where

$A_{k}=\frac{\sqrt{N}}{k!}{2k\choose k}^{1/2}{N+k\choose 2k+1}^{-1/2},\ \ \ \
\ B_{k}=-\frac{k(N+1)\sqrt{N}}{k!}{2k-1\choose k}{N+k\choose
2k+1}^{-1/2}{2k\choose k}^{-1/2}.$ (2.1)

Second, the polynomials $\psi_{k}$ are orthonormal in
$L^{2}(U)\equiv\\{g:{\cal P}_{0}\to\mathbb{R}\\}$, that is,

$\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}\psi_{k}(U)\psi_{m}(U)=\frac{1}{N}\sum_{x=1}^{N}\psi_{k}(x)\psi_{m}(x)=\delta_{k,m}=\left\\{\begin{array}[]{ccc}1,&\mbox{if}&k=m,\\\
0,&\mbox{if}&k\neq m,\end{array}\right.$

and this implies that any function $g:{\cal P}_{0}\to\mathbb{R}$ has a unique
development in Hahn polynomial series, namely,

$g(x)=\sum_{k=0}^{N-1}\delta_{k}\psi_{k}(x),\ \ \ \mbox{where}\ \
\delta_{k}=\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}\psi_{k}(U)g(U)\ \
(k=0,\ldots,N-1).$ (2.2)

The constants $\\{\delta_{k}\\}_{k=0}^{N-1}$ are the Fourier coefficients of
$g$.

An arbitrary ${\cal P}=\\{x_{1},\ldots,x_{N}\\}$ (with $x_{1}\leq\cdots\leq
x_{N}$) is the image of some nondecreasing $g:{\cal P}_{0}\to\mathbb{R}$,
namely, $g(i)=x_{i}$, $i=1,\ldots,N$. In other words, if $X$ is uniformly
distributed over ${\cal P}$ (in the sense of Definition 1.1) and $U$ is
uniformly distributed over ${\cal P}_{0}$ then $X$ and $g(U)$ have the same
distribution. Therefore, $(g(U_{1}),g(U_{2}))$ and $(X_{1},X_{2})$ have the
same distribution, where $X_{1},X_{2}$ are independent copies of $X$, and
$U_{1},U_{2}$ are independent copies of $U$. Since $g$ is nondecreasing, it
holds
$(\min\\{g(U_{1}),g(U_{2})\\},\max\\{g(U_{1}),g(U_{2})\\})=(g(\min\\{U_{1},U_{2}\\}),g(\max\\{U_{1},U_{2}\\}))=(g(U_{1:2}),g(U_{2:2}))$.
This shows that the random pairs $(X_{1:2},X_{2:2})$ and
$(g(U_{1:2}),g(U_{2:2}))$ have the same distribution; hence,
$\rho_{12}=\rho(g(U_{1:2}),g(U_{2:2}))$. Expanding $g$ according to (2.2), we
will be able to express $\rho_{12}$ as a function of
$\delta_{1},\ldots,\delta_{N-1}$, and, next, we shall maximize the resulting
function. The calculations, below, do not a-priori impose any monotonicity
restrictions on $g$; the only assumption is that $g$ is non-constant
(otherwise, the variance of $X$ is zero and $\rho_{12}$ is undefined).

######  Lemma 2.1.

(Hahn representation). Expand an arbitrary $g:{\cal P}_{0}\to\R$ in a Hahn
polynomial series as in (2.2), and let $\sigma_{1}^{2}=\mbox{\rm\hskip
0.86108ptVar\hskip 0.86108pt}g(U_{1:2})$, $\sigma_{2}^{2}=\mbox{\rm\hskip
0.86108ptVar\hskip 0.86108pt}g(U_{2:2})$, $\sigma_{12}=\mbox{\rm\hskip
0.86108ptCov\hskip 0.86108pt}(g(U_{1:2}),g(U_{2:2}))$. Then,

$\displaystyle\sigma_{1}^{2}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{N-1}\delta_{k}^{2}-\frac{\delta_{1}^{2}}{3}\Big{(}1-\frac{1}{N^{2}}\Big{)}-\frac{2}{N}\sum_{k=1}^{N-2}\lambda_{k}\delta_{k}\delta_{k+1}$
(2.3) $\displaystyle\sigma_{2}^{2}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{N-1}\delta_{k}^{2}-\frac{\delta_{1}^{2}}{3}\Big{(}1-\frac{1}{N^{2}}\Big{)}+\frac{2}{N}\sum_{k=1}^{N-2}\lambda_{k}\delta_{k}\delta_{k+1}$
(2.4) $\displaystyle\sigma_{12}$ $\displaystyle=$
$\displaystyle\frac{\delta_{1}^{2}}{3}\Big{(}1-\frac{1}{N^{2}}\Big{)},$ (2.5)

where, in (2.3) and (2.4),

$\lambda_{k}=(k+1)\sqrt{\frac{N^{2}-(k+1)^{2}}{(2k+1)(2k+3)}},\ \ \
k=1,\ldots,N-1.$ (2.6)

Proof: If $\delta_{1}=\cdots=\delta_{N-1}=0$ then $g\equiv\delta_{0}$
(constant) and the result is trivial, since
$\sigma_{1}^{2}=\sigma_{2}^{2}=\sigma_{12}=0$. Assume that
$\beta:=\delta_{1}^{2}+\cdots+\delta_{N-1}^{2}>0$, and set
$h(x):=(g(x)-\delta_{0})/\sqrt{\beta}$, so that

$h(x)=\sum_{k=1}^{N-1}d_{k}\psi_{k}(x),\ \ \mbox{ where }\
d_{k}=\delta_{k}/\sqrt{\beta}\ \ (k=1,\ldots,N-1),$

are the Fourier coefficients of $h$. By construction,
$\sum_{k=1}^{N-1}d_{k}^{2}=1$, $\sigma_{1}^{2}=\beta\mbox{\rm\hskip
0.86108ptVar\hskip 0.86108pt}h(U_{1:2})$, $\sigma_{2}^{2}=\beta\mbox{\rm\hskip
0.86108ptVar\hskip 0.86108pt}h(U_{2:2})$, $\sigma_{12}=\beta\mbox{\rm\hskip
0.86108ptCov\hskip 0.86108pt}(h(U_{1:2}),h(U_{2:2}))$. The orthonormality of
$\psi_{k}$ shows that the random variable $Y:=h(U)$ is standardized; indeed,

$\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Y=\sum_{k=1}^{N-1}d_{k}\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}\psi_{k}(U)=0,\ \ \ \ \mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Y^{2}=\sum_{k=1}^{N-1}\sum_{m=1}^{N-1}d_{k}d_{m}\mbox{\rm\hskip
0.86108ptIE\hskip
0.86108pt}\psi_{k}(U)\psi_{m}(U)=\sum_{k=1}^{N-1}d_{k}^{2}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}\psi_{k}(U)^{2}=1,$

where, the first equality is satisfied because $\psi_{k}$ ($k\geq 1$) is
orthogonal to $\psi_{0}\equiv 1$, and the second one follows from the
orthonormality of $\psi_{k}$ and the definition of $d_{k}$. We now write
$Y_{1}=h(U_{1})$, $Y_{2}=h(U_{2})$, $Z_{1}=h(U_{1:2})$, $Z_{2}=h(U_{2:2})$,
with $U_{1}$, $U_{2}$ being independent uniform from ${\cal P}_{0}$; notice
that the inequality $Z_{1}\leq Z_{2}$ may fail, since $h$ has not been assumed
monotonic. According to the argument above, $Y_{1},Y_{2}$ are independent and
standardized. Taking expectations to the obvious identities
$Z_{1}+Z_{2}=Y_{1}+Y_{2}$ and $Z_{1}Z_{2}=Y_{1}Y_{2}$ we see that
$\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}Z_{1}=-\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}Z_{2}$ and $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{1}Z_{2}=0$. Hence, $\mbox{\rm\hskip 0.86108ptCov\hskip
0.86108pt}(Z_{1},Z_{2})=(\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2})^{2}$; equivalently, $\sigma_{12}=\beta(\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}Z_{2})^{2}$. The relation $2=\mbox{\rm\hskip
0.86108ptVar\hskip 0.86108pt}(Y_{1}+Y_{2})=\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}(Z_{1}+Z_{2})=\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}Z_{1}+\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}Z_{2}+2(\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}Z_{2})^{2}$
shows that $\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}Z_{1}=2-\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}Z_{2}-2(\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2})^{2}=2-(\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2})^{2}-\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}Z_{2}^{2}$
and $\sigma_{1}^{2}=\beta\big{(}2-(\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2})^{2}-\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2}^{2}\big{)}$. Since $\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}Z_{2}=\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2}^{2}-(\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}Z_{2})^{2}$
implies $\sigma_{2}^{2}=\beta\big{(}\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2}^{2}-(\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2})^{2}\big{)}$, all quantities $\sigma_{1}^{2}$,
$\sigma_{2}^{2}$, $\sigma_{12}$, (and $\rho_{12}$), are expressed in terms of
$\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}Z_{2}=\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}h(U_{2:2})$ and $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2}^{2}=\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}h(U_{2:2})^{2}$.

The probability mass function of $U_{2:2}$ is given by $\mbox{\rm\hskip
0.86108ptIP\hskip 0.86108pt}(U_{2:2}=j)=(2j-1)/N^{2}$, $j=1,\ldots,N$, and
therefore,

$\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2}=\frac{1}{N^{2}}\sum_{j=1}^{N}(2j-1)h(j)=\frac{1}{N}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}(2U-1)h(U)=\frac{2}{N}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}Uh(U),$

where the last equality follows from $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}h(U)=\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}Y=0$. Since each
$\psi_{k}(U)$ ($k\geq 2$) is orthogonal to any polynomial of degree at most
one, hence to $U$, we have

$\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}Uh(U)=d_{1}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}U\psi_{1}(U)+\sum_{k=2}^{N-1}d_{k}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}U\psi_{k}(U)=d_{1}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}U\psi_{1}(U)=\frac{d_{1}}{A_{1}}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}\psi_{1}(U)(A_{1}U+B_{1})=\frac{d_{1}}{A_{1}},$

because $A_{1}U+B_{1}=\psi_{1}(U)$ and $\psi_{1}$ is orthogonal to constants.
Substituting the value of $A_{1}$ from (2.1) we obtain $\mbox{\rm\hskip
0.86108ptIE\hskip
0.86108pt}Z_{2}=2d_{1}/(NA_{1})=(1-N^{-2})^{1/2}d_{1}/\sqrt{3}$, and the
relations $\sigma_{12}=\beta(\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2})^{2}$ and $d_{1}=\delta_{1}/\sqrt{\beta}$ imply (2.5).

The computation of the second moment of $Z_{2}$ is more involved. Write

$\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2}^{2}=\frac{1}{N^{2}}\sum_{j=1}^{N}(2j-1)h(j)^{2}=\frac{1}{N}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}(2U-1)h(U)^{2}=\frac{2}{N}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}Uh(U)^{2}-\frac{1}{N},$ (2.7)

because $\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}Y^{2}=\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}h(U)^{2}=1$. It remains to compute

$\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Uh(U)^{2}=\sum_{k,m=1}^{N-1}d_{k}d_{m}\mbox{\rm\hskip
0.86108ptIE\hskip
0.86108pt}U\psi_{k}(U)\psi_{m}(U)=\sum_{k=1}^{N-1}d_{k}^{2}\mbox{\rm\hskip
0.86108ptIE\hskip
0.86108pt}U\psi_{k}(U)^{2}+2\sum_{k=1}^{N-2}d_{k}d_{k+1}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}U\psi_{k}(U)\psi_{k+1}(U);$

all other terms (with $m>k$) vanish, due to the orthogonality of $\psi_{m}(U)$
to the polynomial $U\psi_{k}(U)$ (of degree $k+1$) when $m\geq k+2$. We
proceed to compute $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}U\psi_{k}(U)^{2}$ and $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}U\psi_{k}(U)\psi_{k+1}(U)$. Write
$\psi_{k}(x)=A_{k}x^{k}+B_{k}x^{k-1}+P_{k-2}(x)$, where $P_{k-2}$ is a
polynomial of degree at most $k-2$ ($P_{-1}\equiv 0$). Then,

$\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}U\psi_{k}(U)^{2}=\mbox{\rm\hskip
0.86108ptIE\hskip
0.86108pt}U\psi_{k}(U)(A_{k}U^{k}+B_{k}U^{k-1}+P_{k-2}(U))=A_{k}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}U^{k+1}\psi_{k}(U)+B_{k}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}U^{k}\psi_{k}(U),$

because $\psi_{k}(U)$ is orthogonal to the polynomial $UP_{k-2}(U)$. Next, we
compute

$\displaystyle\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}U^{k+1}\psi_{k}(U)$
$\displaystyle=$ $\displaystyle\frac{1}{A_{k+1}}\mbox{\rm\hskip
0.86108ptIE\hskip
0.86108pt}\psi_{k}(U)\Big{(}A_{k+1}U^{k+1}+B_{k+1}U^{k}+P_{k-1}(U)-B_{k+1}U^{k}-P_{k-1}(U)\Big{)}$
$\displaystyle=$ $\displaystyle\frac{1}{A_{k+1}}\mbox{\rm\hskip
0.86108ptIE\hskip
0.86108pt}\psi_{k}(U)\Big{(}\psi_{k+1}(U)-B_{k+1}U^{k}-P_{k-1}(U)\Big{)}$
$\displaystyle=$ $\displaystyle-\frac{B_{k+1}}{A_{k+1}}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}U^{k}\psi_{k}(U),$

where the last equality follows from the orthogonality of $\psi_{k}$ to both
$\psi_{k+1}$ and $P_{k-1}$. Therefore,

$\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}U\psi_{k}(U)^{2}=A_{k}\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}U^{k+1}\psi_{k}(U)+B_{k}\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}U^{k}\psi_{k}(U)=\left(B_{k}-\frac{A_{k}B_{k+1}}{A_{k+1}}\right)\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}U^{k}\psi_{k}(U).$

The calculation of $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}U^{k}\psi_{k}(U)$ is easy:

$\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}U^{k}\psi_{k}(U)=\frac{1}{A_{k}}\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}\psi_{k}(U)(A_{k}U^{k}+B_{k-1}U^{k-1}+\cdots)=\frac{1}{A_{k}}\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}\psi_{k}(U)^{2}=\frac{1}{A_{k}}.$

Hence, using the fact that $B_{k}/A_{k}=-k(N+1)/2$, see (2.1), we obtain

$\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}U\psi_{k}(U)^{2}=\frac{B_{k}}{A_{k}}-\frac{B_{k+1}}{A_{k+1}}=\frac{N+1}{2};$

surprisingly, this expectation is independent of $k$. Proceeding similarly,
and in view of (2.1), we compute

$\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}U\psi_{k}(U)\psi_{k+1}(U)=\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}U(A_{k}U^{k}+B_{k}U^{k-1}+\cdots)\psi_{k+1}(U)=A_{k}\mbox{\rm\hskip
0.86108ptIE\hskip
0.86108pt}U^{k+1}\psi_{k+1}(U)=\frac{A_{k}}{A_{k+1}}=\frac{\lambda_{k}}{2},$

with $\lambda_{k}$ as in (2.6). Combining the preceding formulae (recall that
$\sum_{k=1}^{N-1}d_{k}^{2}=1$) we obtain

$\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Uh(U)^{2}=\frac{N+1}{2}+\sum_{k=1}^{N-2}\lambda_{k}d_{k}d_{k+1}.$

Then, from (2.7), $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2}^{2}=1+(2/N)\sum_{k=1}^{N-2}\lambda_{k}d_{k}d_{k+1}$. Moreover,
since $(\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2})^{2}=(1-N^{-2})d_{1}^{2}/3$ we have $\mbox{\rm\hskip
0.86108ptVar\hskip
0.86108pt}Z_{2}=1-(1-N^{-2})d_{1}^{2}/3+(2/N)\sum_{k=1}^{N-2}\lambda_{k}d_{k}d_{k+1}$,
and the relation $\sigma_{2}^{2}=\beta\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}Z_{2}$ yields (2.4), because $\beta d_{1}^{2}=\delta_{1}^{2}$ and
$\beta d_{k}d_{k+1}=\delta_{k}\delta_{k+1}$ (recall that
$\beta=\delta_{1}^{2}+\cdots+\delta_{N-1}^{2}$). Similarly, substituting the
preceding expressions for $(\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2})^{2}$ and $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2}^{2}$ to the relation $\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}Z_{1}=2-(\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}Z_{2})^{2}-\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}Z_{2}^{2}$ we
obtain $\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}Z_{1}=1-(1-N^{-2})d_{1}^{2}/3-(2/N)\sum_{k=1}^{N-2}\lambda_{k}d_{k}d_{k+1}$,
and since $\sigma_{1}^{2}=\beta\mbox{\rm\hskip 0.86108ptVar\hskip
0.86108pt}Z_{1}$, (2.3) follows.

######  Lemma 2.2.

We define the rational functions $R_{0}(x)\equiv 1$ and, recurrently,

$R_{k}(x):=1-\frac{k^{2}(1-k^{2}/x)}{(4k^{2}-1)R_{k-1}(x)},\ \ \ \
k=1,2,\ldots.$

Then, the following inequalities hold true:

$\frac{k+1}{2k+1}<R_{k}(x)<1,\ \ x>k^{2},\ \ k=1,2,\ldots.$ (2.8)

Proof: It is easily seen that $R_{1}(x)=(2+x^{-1})/3$, which is strictly
decreasing on $[1,\infty)$, with $R_{1}(1)=1$ and $R_{1}(\infty)=2/3$; hence,
(2.8) is fulfilled for $k=1$. Suppose that (2.8) is true for some $k-1\geq 1$,
that is, $k/(2k-1)<R_{k-1}(x)<1$, $x>(k-1)^{2}$; equivalently,

$1<\frac{1}{R_{k-1}(x)}<\frac{2k-1}{k},\ \ \ x>(k-1)^{2}.$

If $x>k^{2}$, multiplying the preceding inequality by the positive quantity
$k^{2}(1-k^{2}/x)/(4k^{2}-1)$ we get

$L(x):=\frac{k^{2}(1-k^{2}/x)}{4k^{2}-1}<\frac{k^{2}(1-k^{2}/x)}{(4k^{2}-1)R_{k-1}(x)}=1-R_{k}(x)<U(x):=\frac{k(1-k^{2}/x)}{2k+1},\
\ \ x>k^{2}.$

It is easily seen that $L$ is strictly increasing in $x$, $x>k^{2}$, so
$L(x)>L(k^{2})=0$, and $U$ is also strictly increasing with
$U(x)<U(\infty)=k/(2k+1)$. This shows that $0<1-R_{k}(x)<k/(2k+1)$ for
$x>k^{2}$ and verifies the inductional step, completing the proof.

######  Lemma 2.3.

For $N\geq 2$, there exist nonnegative numbers $\\{\alpha_{k}\\}_{k=1}^{N-1}$
and $\\{\beta_{k}\\}_{k=1}^{N-1}$ satisfying the following four properties:

  * (i)

$\alpha_{1}=\frac{1}{3}\Big{(}2+\frac{1}{N^{2}}\Big{)}$, $\beta_{N-1}=0$;

  * (ii)

$\alpha_{k}>\beta_{k-1}>0$, $k=2,\ldots,N-1$;

  * (iii)

$\alpha_{k}+\beta_{k-1}=1$, $k=2,\ldots,N-1$;

  * (iv)

$\alpha_{k}\beta_{k}=\frac{(k+1)^{2}}{(2k+1)(2k+3)}\Big{(}1-\frac{(k+1)^{2}}{N^{2}}\Big{)}$,
$k=1,\ldots,N-1$.

Proof: Define $\alpha_{k}:=R_{k}(N^{2})$ and $\beta_{k}:=1-R_{k+1}(N^{2})$,
$k=1,\ldots,N-1$, where the functions $R_{k}$ are defined in Lemma 2.2. Then,
(i) and (iii) are obvious, and (ii) follows from (2.8), since
$\alpha_{k}=R_{k}(N^{2})>(k+1)/(2k+1)>1/2$, because $N^{2}>k^{2}$, and, by the
same reasoning, $\beta_{k}\in[0,1/2)$ (more precisely, $\beta_{k}>0$ unless
$k=N-1$). Finally, (iv) is trivial for $k=N-1$ (since $\beta_{N-1}=0$), and
for $k=1,\ldots,N-2$,

$\alpha_{k}\beta_{k}=\alpha_{k}(1-\alpha_{k+1})=\alpha_{k}\frac{(k+1)^{2}(1-(k+1)^{2}/N^{2})}{[4(k+1)^{2}-1]\alpha_{k}}=\frac{(k+1)^{2}}{(2k+1)(2k+3)}\Big{(}1-\frac{(k+1)^{2}}{N^{2}}\Big{)},$

due to the recurrent relation of Lemma 2.2, namely,

$1-R_{k+1}(N^{2})=\frac{(k+1)^{2}(1-(k+1)^{2}/N^{2})}{[4(k+1)^{2}-1]R_{k}(N^{2})}.$

Hence, the lemma is proved.

######  Lemma 2.4.

(Terrell-Hahn representation). Consider the numbers
$\\{\alpha_{k}\\}_{k=1}^{N-1}$ and $\\{\beta_{k}\\}_{k=1}^{N-1}$ as in Lemma
2.3. Then, with the convention $\delta_{N}:=0$,

$\sigma_{1}^{2}=\sum_{k=1}^{N-1}\left(\sqrt{\alpha_{k}}\delta_{k}-\sqrt{\beta_{k}}\delta_{k+1}\right)^{2},\
\ \ \
\sigma_{2}^{2}=\sum_{k=1}^{N-1}\left(\sqrt{\alpha_{k}}\delta_{k}+\sqrt{\beta_{k}}\delta_{k+1}\right)^{2},$

where $\sigma_{1}^{2}=\mbox{\rm\hskip 0.86108ptVar\hskip 0.86108pt}g(U_{1:2})$
and $\sigma_{2}^{2}=\mbox{\rm\hskip 0.86108ptVar\hskip 0.86108pt}g(U_{2:2})$
are as in Lemma 2.1.

Proof: Let $s_{1}^{2}$ be the first sum. Expanding the squares we have

$s_{1}^{2}=\sum_{k=1}^{N-1}\alpha_{k}\delta_{k}^{2}+\sum_{k=1}^{N-1}\beta_{k}\delta_{k+1}^{2}-2\sum_{k=1}^{N-1}\sqrt{\alpha_{k}\beta_{k}}\delta_{k}\delta_{k+1}=\alpha_{1}\delta_{1}^{2}+\sum_{k=2}^{N-1}(\alpha_{k}+\beta_{k-1})\delta_{k}^{2}-2\sum_{k=1}^{N-2}\sqrt{\alpha_{k}\beta_{k}}\delta_{k}\delta_{k+1},$

and Lemma 2.3 implies that

$s_{1}^{2}=\frac{1}{3}\Big{(}2+\frac{1}{N^{2}}\Big{)}\delta_{1}^{2}+\sum_{k=2}^{N-1}\delta_{k}^{2}-\frac{2}{N}\sum_{k=1}^{N-2}\lambda_{k}\delta_{k}\delta_{k+1},$

with $\lambda_{k}$ as in (2.6). Observing that $(2+N^{-2})/3=1-(1-N^{-2})/3$,
$s_{1}^{2}$ reduces to the expression of $\sigma_{1}^{2}$, given in (2.3). The
derivation of the expression for $\sigma_{2}^{2}$ is completely similar.

######  Theorem 2.1.

For $N\geq 2$ and any nonconstant $g:{\cal P}_{0}=\\{1,\ldots,N\\}\to\R$,

$\rho(g(U_{1:2}),g(U_{2:2}))\leq\frac{1-N^{-2}}{2+N^{-2}},$

with equality if and only if $g$ is linear.

Proof: Let $g(U)=\sum_{k=0}^{N-1}\delta_{k}\psi_{k}(U)$ be the Hahn
representation of $g$, with $\delta_{1}^{2}+\cdots+\delta_{N-1}^{2}>0$
(otherwise, $g$ is constant). If $\sigma_{i}^{2}=\mbox{\rm\hskip
0.86108ptVar\hskip 0.86108pt}g(U_{i:2})$ ($i=1,2$) and
$\sigma_{12}=\mbox{\rm\hskip 0.86108ptCov\hskip
0.86108pt}(g(U_{1:2}),g(U_{2:2}))$ then, from Lemma 2.4 (and with the
convention $\delta_{N}=0$),

$\displaystyle\sigma_{1}^{2}\sigma_{2}^{2}$ $\displaystyle=$
$\displaystyle\Bigg{\\{}\sum_{k=1}^{N-1}\left(\sqrt{\alpha_{k}}\delta_{k}-\sqrt{\beta_{k}}\delta_{k+1}\right)^{2}\Bigg{\\}}\Bigg{\\{}\sum_{k=1}^{N-1}\left(\sqrt{\alpha_{k}}\delta_{k}+\sqrt{\beta_{k}}\delta_{k+1}\right)^{2}\Bigg{\\}}$
$\displaystyle\geq$
$\displaystyle\Bigg{\\{}\sum_{k=1}^{N-1}\left(\sqrt{\alpha_{k}}\delta_{k}-\sqrt{\beta_{k}}\delta_{k+1}\right)\left(\sqrt{\alpha_{k}}\delta_{k}+\sqrt{\beta_{k}}\delta_{k+1}\right)\Bigg{\\}}^{2}$
$\displaystyle=$
$\displaystyle\Bigg{\\{}\sum_{k=1}^{N-1}\left(\alpha_{k}\delta_{k}^{2}-\beta_{k}\delta_{k+1}^{2}\right)\Bigg{\\}}^{2}$
$\displaystyle=$
$\displaystyle\Bigg{\\{}\alpha_{1}\delta_{1}^{2}+\sum_{k=2}^{N-1}\left(\alpha_{k}-\beta_{k-1}\right)\delta_{k}^{2}\Bigg{\\}}^{2}$
$\displaystyle\geq$ $\displaystyle\alpha_{1}^{2}\delta_{1}^{4},\Bigg{.}$

where the second inequality holds true because $\alpha_{k}-\beta_{k-1}>0$ and
$\alpha_{1}>0$, see Lemma 2.3(ii), while the first one is a simple application
of the well-known Cauchy inequality. On substituting the value of $\alpha_{1}$
from Lemma 2.3(i), the previous argument verifies the inequality

$\left(\sigma_{1}^{2}\sigma_{2}^{2}\right)^{1/2}\geq\frac{1}{3}\Big{(}2+\frac{1}{N^{2}}\Big{)}\delta_{1}^{2},$

in which the equality holds if and only if $\delta_{k}=0$ for all $k\geq 2$,
i.e., when $g$ is linear. From (2.5) of Lemma 2.1,
$\delta_{1}^{2}=3\sigma_{12}/(1-N^{-2})$; thus,
$(\sigma_{1}^{2}\sigma_{2}^{2})^{1/2}\geq\sigma_{12}(2+N^{-1})/(1-N^{-2})$,
and this is equivalent to the desired inequality.

Applying Theorem 2.1 to the particular case where $g$ is given by
$g(i)=x_{i}$, $i=1,\ldots,N$ (which is nondecreasing and non-constant), we
conclude the result of Theorem 1.2.

##  3 Concluding Remarks

######  Remark 3.1.

(An improvement to the covariance/expectation bound for $n=2$). Since the
variance of $X=g(U)$ equals to
$\sigma^{2}=\delta_{1}^{2}+\cdots+\delta_{N-1}^{2}\geq\delta_{1}^{2}$, an
immediate consequence of (2.5) is the Bessel-type inequality $\mbox{\rm\hskip
0.86108ptCov\hskip 0.86108pt}(X_{1:2},X_{2:2})\leq\sigma^{2}(1-N^{-2})/3$;
this is similar to the bound $\mbox{\rm\hskip 0.86108ptCov\hskip
0.86108pt}(X_{1:2},X_{2:2})\leq\sigma^{2}(1-2N^{-1})/(3(1-N^{-1}))$, obtained
by Balakrishnan et al. (2003) in the without-replacement case. Notice that the
equality in any of these two bounds is attained only for lattice populations.
Both bounds converge to $\sigma^{2}/3$ as $N\to\infty$, and the maximizing
population (provided $\mu=\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}X$ and
$\sigma^{2}$ are held fixed), converges weakly to the uniform ${\cal
U}(\mu-\sigma\sqrt{3},\mu+\sigma\sqrt{3})$ distribution. Papathanasiou (1990)
proved that the inequality $\mbox{\rm\hskip 0.86108ptCov\hskip
0.86108pt}(X_{1:2},X_{2:2})\leq\sigma^{2}/3$ holds true for random samples
$X_{1},X_{2}$ from any distribution with finite variance $\sigma^{2}$, and the
equality characterizes the location-family of rectangular distributions with
fixed variance. Later on, Balakrishnan and Balasubramanian (1993) observed
that, in view of the obvious identity $\mbox{\rm\hskip 0.86108ptCov\hskip
0.86108pt}(X_{1:2},X_{2:2})=(\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}X_{2:2}-\mu)^{2}$, the supremum of this covariance equals to the
smallest upper bound for $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}X_{2:2}-\mu$ in terms of $\sigma$. The form of the latter bound,
called HDG bound, is long known for $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}X_{n:n}$ (any $n\geq 2$), and it is proved in two classical papers
(appeared in consecutive pages of the same journal!) by Hartley and David
(1954) and Gumbel (1954). The particular case ($n=2$) of the HDG bound yields
the inequality $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}X_{2:2}\leq\mu+\sigma/\sqrt{3}$, in which the equality characterizes
${\cal U}(\mu-\sigma\sqrt{3},\mu+\sigma\sqrt{3})$. Clearly, this is equivalent
to Papathanasiou’s covariance bound. The improvement on the HDG bound for
uniform populations (with $N$ elements, in the sense of Definition 1.1) is an
immediate consequence of the discrete covariance bound. More precisely,

$\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}X_{2:2}=\mu+(\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}X_{2:2}-\mu)=\mu+\Big{(}\mbox{\rm\hskip
0.86108ptCov\hskip
0.86108pt}(X_{1:2},X_{2:2})\Big{)}^{1/2}\leq\mu+(1-N^{-2})^{1/2}\
\frac{\sigma}{\sqrt{3}},$

and the equality characterizes the lattice population with mean $\mu$ and
variance $\sigma^{2}$, that is,
$x_{k}=\mu+\sigma\sqrt{3}(2k-N-1)/\sqrt{N^{2}-1}$, $k=1,\ldots,N$.

######  Remark 3.2.

(Connection to continuous case). One may define the functional
$\rho_{12}(\cdot)$ with domain ${\cal F}:=\\{F:F$ is a distribution function
with finite non-zero variance$\\}$ as

$\rho_{12}(F):=\rho(X_{1:2},X_{2:2}),\ \ \ \mbox{$X_{1}$, $X_{2}$ independent
with distribution $F$},$

where $X_{1:2}\leq X_{2:2}$ are the order statistics. For $N\geq 2$, we denote
by ${\cal F}_{N}:=\\{{\cal P}:{\cal P}=\\{x_{1},\ldots,x_{N}\\}$ is any
uniform population with at least two distinct elements$\\}$ (uniform in the
sense of Definition 1.1). Then it can be verified that ${\cal
F}_{0}:=\cup_{N\geq 2}{\cal F}_{N}$ is dense in ${\cal F}$ with respect to the
weak topology of probability measures, and the functional $\rho_{12}$ is
continuous. Theorem 2.1 shows that $\sup_{F\in{\cal F}_{0}}\rho_{12}(F)=1/2$,
hence, the same is true for $\sup_{F\in{\cal F}}$, and this implies the bound
in Terrell’s result, Theorem 1.1.

######  Remark 3.3.

(An equivalent minimization problem). To simplify notation, set
$\gamma_{k}=2\lambda_{k}/N$ ($k=1,\ldots,N-1$) with $\lambda_{k}$ as in (2.6),
and $\lambda=\alpha_{1}$ as in Lemma 2.3(i); that is,

$\lambda=\frac{1}{3}\Big{(}2+\frac{1}{N^{2}}\Big{)},\ \ \
\gamma_{k}=\frac{2(k+1)}{N}\sqrt{\frac{N^{2}-(k+1)^{2}}{(2k+1)(2k+3)}},\ \
k=1,\ldots,N-1.$

According to Lemma 2.1, we have

$\Big{(}\frac{1}{\rho_{12}}\Big{)}^{2}=\frac{9N^{4}}{(N^{2}-1)^{2}}\
\frac{\Big{(}\lambda\delta_{1}^{2}+\delta_{2}^{2}+\cdots+\delta_{N-1}^{2}\Big{)}^{2}-\Big{(}\gamma_{1}\delta_{1}\delta_{2}+\gamma_{2}\delta_{2}\delta_{3}+\cdots+\gamma_{N-2}\delta_{N-2}\delta_{N-1}\Big{)}^{2}}{\delta_{1}^{4}},$

where this expression should be treated as $+\infty$ for $\delta_{1}=0$. The
main difficulty in the proof of Theorem 2.1 was to minimize the above ratio
w.r.t.
$(\delta_{1},\ldots,\delta_{N-1})\in\mathbb{R}^{N-1}\setminus\\{\mbox{\boldmath$0$}\\}$,
and the main result showed that the minimizing points are coincide with the
axis $\\{(\delta_{1},0,\ldots,0),\delta_{1}\neq 0\\}$. The ratio is scale
invariant, so, dividing each $\delta_{k}$ by $\delta_{1}\neq 0$ (since
$\delta_{1}=0$ makes the ratio infinite, and, certainly, cannot be a
minimizer) we can assume that $\delta_{1}=1$. Then, with a slight change of
notation (remove the constant scalar and set $n=N-1$, $x_{k}=\delta_{k}$), we
define

$f(x_{2},\ldots,x_{n}):=\Big{(}\lambda+x_{2}^{2}+\cdots+x_{n}^{2}\Big{)}^{2}-\Big{(}\gamma_{1}x_{2}+\gamma_{2}x_{2}x_{3}+\cdots+\gamma_{n-1}x_{n-1}x_{n}\Big{)}^{2}.$

Consequently, it remains to minimize the (four degree polynomial) function $f$
and to prove that its global minimum is uniquely attained at $0$, so that
$f\geq f(\mbox{\boldmath$0$})=\lambda^{2}$. Despite the fact that this
(unrestricted) minimization problem looks like a simple exercise, this is not
so; the minimization points (if exist) depend on the parameters $\lambda$,
$\gamma_{1},\ldots,\gamma_{n-1}$. To see this, set $n=3$, $\lambda_{0}=11/6$
(since $\lambda=\lambda_{0}$ when $n=3$), and write
$f(x_{2},x_{3};\gamma_{1},\gamma_{2})=(\lambda_{0}+x_{2}^{2}+x_{3}^{2})^{2}-\big{(}\gamma_{1}x_{2}+\gamma_{2}x_{2}x_{3}\big{)}^{2}$.
It can be checked that the function $f(x_{2},x_{3};4/3,2/3)$, though positive,
is globally minimized at the points $(x_{2},x_{3})$ with

$x_{2}=\pm\Big{\\{}66\sqrt{6}-81\Big{\\}}^{1/2}\Big{/}16,\ \ \
x_{3}=(3\sqrt{6}-5)\Big{/}16,$

not at $0$. Also, $f(\cdot\ ;1,2)$ is not bounded below; the minimum value of
$f(\cdot\ ;1,8/5)$ is negative; $f(\cdot\ ;1,1)$ attains its global minimum at
$0$. These examples (in the simplest case $n=3$) indicate that the multi-
parameter case is too complicated, to enable an explicit solution. In the
contrary, the proposed Terrell-type method of proof succeed in obtaining
Theorems 1.2, 2.1. Motivated from these examples, a natural question of
mathematical nature would ask for (necessary and) sufficient conditions on the
arbitrary positive parameters $\lambda$, $\gamma_{1},\ldots,\gamma_{n-1}$,
guaranteing that the minimum of $f(x_{2},\ldots,x_{n})$, above, is uniquely
attained at $0$.

######  Remark 3.4.

(Connection to maximal correlation). It is important to make clear that the
assertion of Theorem 2.1 cannot be reduced to a maximal correlation problem
(unless $N=2$), in contrast to all of the existing results mentioned in
Section 1. To see this, it is useful to recall Rényi’s (1959) characterization
of the maximal correlation of $(X,Y)$, namely,

$R(X,Y)^{2}=\mbox{$\sup_{f}$}\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}\Big{[}\big{(}\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}f(X)|Y\big{)}^{2}\Big{]},$

where the supremum is taken over $f$ with $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}f(X)=0$, $\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}f(X)^{2}=1$.
Applying Rényi’s result to $(X,Y)=(U_{1:2},U_{2:2})$, and writing
$f(k)=x_{k}$, $k=1,\ldots,N$, it is seen that

$R^{2}:=R(U_{1:2},U_{2:2})^{2}=\frac{1}{N^{2}}\max\sum_{j=1}^{N}\frac{1}{2j-1}\left\\{2(x_{1}+\cdots+x_{j-1})+x_{j}\right\\}^{2},$

where the maximum is taken over $x_{1},\ldots,x_{N}$ satisfying

$\sum_{k=1}^{N}(2N-2k+1)x_{k}=0,\ \ \ \ \
\sum_{k=1}^{N}(2N-2k+1)x_{k}^{2}=N^{2}.$

For $N=3$ one finds a maximization point explicitly, namely,

$x_{1}^{*}=-\sqrt{\frac{2}{5}+\frac{13}{10\sqrt{19}}},\ \ \ \ \
x_{2}^{*}=\sqrt{1-\frac{2}{\sqrt{19}}},\ \ \ \ \
x_{3}^{*}=\sqrt{4-\frac{1}{2\sqrt{19}}},$

yielding $R=(2+\sqrt{19})/15>(1-N^{-2})/(2+N^{-2})=8/19$. Setting
$f(k)=x_{k}^{*}$ and $g(Y)=(1/R)\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}(f(X)|Y)$, that is, $g(k)=y_{k}^{*}$ with

$y_{1}^{*}=-\sqrt{4-\frac{1}{2\sqrt{19}}},\ \ \ \ \
y_{2}^{*}=-\sqrt{1-\frac{2}{\sqrt{19}}},\ \ \ \ \
y_{3}^{*}=\sqrt{\frac{2}{5}+\frac{13}{10\sqrt{19}}},$

if is easy to check that the pair $(f,g)$ is maximally correlated:
$\mbox{\rm\hskip 0.86108ptIE\hskip 0.86108pt}f(U_{1:2})=\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}g(U_{2:2})=0$, $\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}f(U_{1:2})^{2}=\mbox{\rm\hskip 0.86108ptIE\hskip
0.86108pt}g(U_{2:2})^{2}=1$, and $\rho(f(U_{1:2}),g(U_{2:2}))=\mbox{\rm\hskip
0.86108ptIE\hskip 0.86108pt}f(U_{1:2})g(U_{2:2})=(2+\sqrt{19})/15$.

The situation is similar for all $N$, namely, it can be shown that
$\rho_{N}:=(1-N^{-2})/(2+N^{-2})<R_{N}<1/2$, $N\geq 3$, where $R_{N}$ is the
maximal correlation of an ordered pair from the uniform population, ${\cal
P}_{0}$, of size $N$. Clearly, for $N\geq 4$ it would be quite surprising if
$R_{N}$ could be expressed in a closed form. Since the maximal correlation is
attainable in our case (the corresponding operator is compact; see Rényi’s,
1959), the TSF inequality (1.1) implies that $R_{N}<1/2$. Certainly,
$R_{N}\geq\rho_{N}$. In order to prove that $R_{N}>\rho_{N}$, it suffices to
find a pair of functions with coefficient of correlation greater than
$\rho_{N}$. Testing small perturbations of linear functions, we were led to
define $f_{0}(x):=x-3x^{2}/N^{3}$, $g_{0}(y):=y+3y^{2}/N^{3}$ (both strictly
increasing). Set $\rho_{0}=\rho(f_{0}(U_{1:2}),g_{0}(U_{2:2}))$,
$\sigma_{1}^{2}=\mbox{\rm\hskip 0.86108ptVar\hskip 0.86108pt}f_{0}(U_{1:2})$,
$\sigma_{2}^{2}=\mbox{\rm\hskip 0.86108ptVar\hskip 0.86108pt}g_{0}(U_{2:2})$.
For $N\geq 3$ we computed (with the help of Mathematica)

$\displaystyle\frac{60^{2}\sigma_{1}^{2}\sigma_{2}^{2}N^{14}(2N^{2}+1)^{2}}{(N^{2}-1)^{2}(N^{2}-4)}\left(\rho_{0}^{2}-\rho_{N}^{2}\right)$
$\displaystyle=$ $\displaystyle-666-1332N+846N^{2}+2952N^{3}+1077N^{4}$
$\displaystyle-828N^{5}-934N^{6}-792N^{7}-485N^{8}+34N^{10}+20N^{12}.$

The substitution $N=n+3$ ($n\geq 0$) shows that the right-hand side, $w(n)$,
is strictly positive:

$\displaystyle w(n)$ $\displaystyle=$ $\displaystyle
7010100+35183016n+72768816n^{2}+86119956n^{3}+66523137n^{4}+35823456n^{5}$
$\displaystyle+13910474n^{6}+3946848n^{7}+815185n^{8}+119820n^{9}+11914n^{10}+720n^{11}+20n^{12}.$

Hence, $R_{N}\geq\rho_{0}>\rho_{N}$ for every $N\geq 3$.

The main result of the present note, although quite specialized, gives rise to
a number of questions that, at least to author’s view, are of some interest.
Let $U_{1:n}\leq\cdots\leq U_{n:n}$ be the ordered sample from the uniform
population ${\cal P}_{0}=\\{1,\ldots,N\\}$. For $1\leq i<j\leq n$ we define

$R_{N}=R_{N}(\mbox{$i$,$j$:$n$}):=\mbox{$\sup_{f,g}$}\
\rho(f(U_{i:n}),g(U_{j:n})),\ \ \ \ \ \
R^{\prime}_{N}=R^{\prime}_{N}(\mbox{$i$,$j$:$n$}):=\mbox{$\sup_{g}$}\
\rho(g(U_{i:n}),g(U_{j:n})),$

where the suprema are taken over non-constant $f$, $g$ (with domain ${\cal
P}_{0}$). The first one is the well-known maximal correlation; the second
appeared in Theorem 2.1 for $i=1$, $j=2$, $n=2$. Moreover, define

$\rho_{N}=\rho_{N}(\mbox{$i$,$j$:$n$}):=\rho(U_{i:n},U_{j:n}),\ \ \ \ \ \
R^{\prime\prime}_{N}=R^{\prime\prime}_{N}(\mbox{$i$,$j$:$n$}):=\mbox{$\sup_{g\nearrow}$}\
\rho(g(U_{i:n}),g(U_{j:n})),$

where the supremum is taken over nondecreasing non-constant $g$. It is obvious
that for all $i$, $j$, $n$, $N$, $0<\rho_{N}\leq R_{N}^{\prime\prime}\leq
R_{N}^{\prime}\leq R_{N}<[i(n+1-j)]^{1/2}[j(n+1-i)]^{-1/2}$; see (1.1). Our
preliminary calculations indicate the relation
$\rho_{N}=R_{N}^{\prime\prime}=R_{N}^{\prime}<R_{N}$ for all $i$, $j$, $n$,
$N$, e.g., $R_{N}^{\prime}(\mbox{$1$,$3$:$3$})=(3-7N^{-2})/(9-N^{-2})$
(checked for small $N$). If the foregoing relation holds true for all choices
of $i$, $j$, $n$, $N$, it will verify that the population ${\cal P}_{0}$
dominates the coefficient of correlation of order statistics from any uniform
population ${\cal P}=\\{x_{1},\ldots,x_{N}\\}$, as in the continuous case (as
well as the without-replacement one). An extra difficulty here is due to the
lack of a simple formula for $\rho_{N}(\mbox{$i$,$j$:$n$})$.

Another natural (but more complicated) extension would be to relax the
uniformity assumption. One possible way to do this is to endow the population
with a fixed probability vector ${\mbox{\boldmath$p$}}=(p_{1},\ldots,p_{N})$
with $p_{k}>0$, $k=1,\ldots,N$. Then ${\cal
P}_{0}({\mbox{\boldmath$p$}})=\\{1,\ldots,N\\}$ is as in the uniform case,
but, now, the point $k$ has probability $p_{k}$. The problem is to obtain the
corresponding results (in an obvious notation) for some of the quantities
$\rho_{N}(\mbox{$i$,$j$:$n$};{\mbox{\boldmath$p$}})$,
$R^{\prime\prime}_{N}(\mbox{$i$,$j$:$n$};{\mbox{\boldmath$p$}})$,
$R^{\prime}_{N}(\mbox{$i$,$j$:$n$};{\mbox{\boldmath$p$}}_{N})$,
$R_{N}(\mbox{$i$,$j$:$n$};{\mbox{\boldmath$p$}})$. The value of
$R^{\prime\prime}_{N}$ (and $R_{N}^{\prime}$) is of special interest, since it
provides the least upper bound of the correlation of order statistics among
all populations ${\cal P}({\mbox{\boldmath$p$}})$ endowed with $p$, i.e.,
containing $N$ points $x_{1}\leq\cdots\leq x_{N}$, $x_{1}<x_{N}$, and the
probability of $x_{k}$ is $p_{k}$. Here, in analogy to Definition 1.1, if ties
appear in the population, e.g., if
$x_{k-1}<x_{k}=x_{k+1}=\cdots=x_{s}<x_{s+1}$, then the probability of
$y=x_{k}$ is defined to be equal to $p_{k}+\cdots+p_{s}$. In this setup, it
would be of interest to characterize the population that maximizes the
correlation of, e.g., the extreme order statistics,
$\rho\big{(}X_{1:n},X_{n:n};{\mbox{\boldmath$p$}}\big{)}$. Clearly, the actual
value of $R^{\prime\prime}_{N}$ depends on $p$, and even the form of
$R^{\prime}_{3}(\mbox{$1$,$2$:$2$};{\mbox{\boldmath$p$}})$ is quite
complicated. To see this, consider the (binomial-type) vectors
${\mbox{\boldmath$p$}}=(1/4,1/2,1/4)$ and
${\mbox{\boldmath$p$}}=(1/16,3/8,9/16)$. For the first case we calculate
$\rho_{3}=R^{\prime}_{3}=R^{\prime\prime}_{3}=9/23$, while for the second case
one finds
$\rho_{3}=169\sqrt{3}/\sqrt{655027}=0.361674...<0.362354...=R^{\prime}_{3}=R^{\prime\prime}_{3}$,
where $R_{3}^{\prime}$ is the unique positive root of an even, 6th-degree,
polynomial, with (large) integer coefficients. If all $p_{k}$ are rational
numbers, $p_{k}=\alpha_{k}/\beta_{k}$, we can write $p_{k}=c_{k}/M$ (in lowest
terms). Then, Theorem 2.1 (applied with ties) yields the inequality
$R_{N}^{\prime}(\mbox{$1$,$2$:$2$};{\mbox{\boldmath$p$}})\leq
R^{\prime}_{M}(\mbox{$1$,$2$:$2$})=(1-M^{-2})/(2+M^{-2})$, a bound that is
(often) poor, since $M$ is usually much larger than $N$. For instance, for the
previous (binomial-type) examples one finds the inaccurate bounds
$R_{3}^{\prime}\leq 5/11$ ($M=4$) and $R_{3}^{\prime}\leq 85/171=0.497...$
($M=16$), respectively.

Finally, fix $N\geq 2$ and define

${\cal M}_{N}:=\Big{\\{}{\mbox{\boldmath$p$}}=(p_{1},\ldots,p_{N}):0\leq
p_{k}<1\ \ (k=1,\ldots,N),\ \ p_{1}+\cdots+p_{N}=1\Big{\\}}.$

It is fairly expected that, as $p$ varies in ${\cal M}_{N}$, the coefficient
of correlation of any two order statistics from ${\cal
P}_{0}({\mbox{\boldmath$p$}})$ will be maximized in the uniform case,
${\mbox{\boldmath$p$}}=(1/N,\ldots,1/N)$, but no proof is available; we also
conjecture that a similar result is true for the without-replacement case.

## References

  * [1] Balakrishnan, N.; Balasubramanian, K. (1993). Equivalence of Hartley–David–Gumbel and Papathanasiou bounds and some further remarks. Statist. Probab. Lett., 16, 39–41.
  * [2] Balakrishnan, N.; Charalambides, C.; Papadatos, N. (2003). Bounds on expectation of order statistics from a finite population., J. Statist. Plann. Inference, 113, 569–588.
  * [3] Gumbel, E.J. (1954). The maxima of the mean largest value and of the range. Ann. Math. Statist., 25, 76–84.
  * [4] Gebelein, H. (1941). Das Statistische Problem der Korrelation als Variation und Eigenwertproblem und sein Zusammenhang mit der Ausgleichrechnung. Angew. Math. Mech., 21, 364–379.
  * [5] Hartley, H.O.; David, H.A. (1954). Universal bounds for mean range and extreme observations. Ann. Math. Statist., 25, 85–99.
  * [6] López-Blázquez, F.; Castaño-Martínez, A. (2006). Upper and lower bounds for the correlation ratio of order statistics from a sample without replacement. J. Statist. Plann. Inference, 136, 43–52.
  * [7] Nevzorov, V.B. (1992). A characterization of exponential distributions by correlations between records. Math. Meth. Statist., 1, 49–54.
  * [8] Papadatos, N. (2014). Some counterexamples concerning maximal correlation and linear regression. J. Multivariate Anal., 126, 114–117.
  * [9] Papadatos, N.; Xifara, T. (2013). A simple method for obtaining the maximal correlation coefficient and related characterizations. J. Multivariate Anal., 118, 102–114.
  * [10] Papathanasiou, V. (1990). Some characterizations of distributions based on order statistics. Statist. Probab. Lett., 9, 145–147.
  * [11] Rényi, A. (1959). On measures of dependence. Acta Math. Acad. Sci. Hungar., 10, 441–451.
  * [12] Székely, G.J.; Móri, T.F. (1985). An extremal property of rectangular distributions. Statist. Probab. Lett., 3, 107–109.
  * [13] Terrell, G.R. (1983). A characterization of rectangular distributions. Ann. Probab., 11(3), 823–826.