\section{Introduction}
\label{sec:intro}

Image super-resolution (SR) refers the process to recover high-resolution (HR)
images from low-resolution (LR) inputs. It is an important image
processing technique to enhance image quality which
subsequently helps to improve higher-level computer vision
tasks~\cite{dai_wacv_2016} \cite{haris_arxiv_2018}.
Over the years, many classical SR methods have been
proposed to successfully use various levels
of features like statistics~\cite{pickup_cj_2009},
edges~\cite{sun_cvpr_2008} \cite{tai_cvpr_2010}
and patches~\cite{wang_iccv_2005} \cite{yang_cvpr_2008} to restore HR images from LR inputs.
While there are also methods developed for SR
using multiple frames~\cite{khattab_iccic_2018},
the scope of introduction here is limited to
single image super-resolution (SISR).

More recently, the powerful deep learning
techniques have led to developments of
many deep learning based SR models
\cite{dong_eccv_2014} \cite{kim_cvpr_2016_2} \cite{ledig_cvpr_2017} \cite{lim_cvprw_2017} \cite{zhang_cvpr_2018} \cite{zhang_eccv_2018}.
These deep learning models commonly rely on a large set of synthetic
training image pairs, where the LR input is downsampled from
the HR reference image using bicubic
interpolation with antialiasing filters.  Common image quality metrics
used to assess performance of these SR models include peak
signal-to-noise ratio (PSNR) and the structural
similarity index (SSIM) \cite{wang_tip_2004},
both emphasizing image restoration fidelity by comparing to the HR reference.
This may lead to SR results of high PSNR values but lack
of HR details perceptually.
Lately, a new metric LPIPS \cite{zhang_r_cvpr_2018}
is proposed to apply image features extracted from pretrained
Alexnet \cite{krizhevsky_nips_2012} to compare two images.
The smaller LPIPS is, the closer the generated SR image is to
the HR reference perceptually.
With advancements in Generative Adversarial Nets (GAN)
\cite{goodfellow_nips_2014}, SR models trained using GAN
\cite{wang_eccv_2018} \cite{ma_cvpr_2020} \cite{ji_cvpr_2020}
have achieved the best performance of image perceptual
quality as compared to LPIPS.

In the past few years, neural architecture search (NAS) that aims to find the optimal network
structure has received a lot of attention\cite{liu_arxiv_2017}
\cite{zoph_arxiv_2016} \cite{zoph_cvpr_2018} \cite{liu_eccv_2018}. It effectively boosts the
SOTA in many typical computer vision problems such as image classification\cite{liu_arxiv_2018},
object detection\cite{liang_iclr_2020}, segmentation\cite{chen_pami_2017} and so on.
Most recently, some researches also begin to apply NAS for image SR
problems \cite{chu_xx_arxiv_2019} \cite{song_dh_aaai_2020} \cite{guo_y_arxiv_2020} with impressive results using
efficient SR models.

In general, the SR problem is ill-posed as there are
multiple HR images corresponding to a single LR image even
when the LR image is constructed from the HR reference using bicubic interpolation
without added noise.
This ambiguity increases when the blur kernel and noise statistics of the LR
are not known, and is even more prominent
in real image SR problems where the LR image is not constructed
from the HR reference.  With the increased uncertainty,
it is common to see different deep learning SR models lead
to different versions of the
restored HR images for a single LR image, especially
when the network architectures are quite different.
To achieve the best performance of the real image SR problem
set forth by the AIM 2020 challenge,
a new fusion scheme is proposed in this study to
generate the final SR output using multi-level ensemble from a suite of
heterogeneous deep learning models that are obtained by applying NAS approach.
The main contributions of the proposed method include:

\begin{itemize}

\item[$\bullet$] A Gaussian Process based NAS (GP-NAS) is first utilized for super-resolution
  with specially designed search space, which can efficiently search and obtain the key
  architecture related parameters and can yield multiple candidate models.

\item[$\bullet$] A multi-level ensemble scheme is proposed in testing, including self-ensemble
  for patches, as well as patch-ensemble and model-ensemble for full-size images.

\item[$\bullet$] The proposed method was applied for the AIM 2020 Real Image Super-Resolution Challenge
  and won the first place in all three tracks (upscaling factors of $\times$2, $\times$3
  and $\times$4) with a comfortable margin in both PSNR and SSIM.

\end{itemize}

\section{Related Works}
\label{sec:rwork}

\noindent\textbf{Deep Learning for Single Image Super-Resolution.} As
the first successful application, Dong \textit{et al}.~\cite{dong_eccv_2014}
proposed a deep CNN model for end-to-end LR to HR mapping and showed
that the training of the neural network is equivalent to global
optimization of traditional sparse-coding-based SR methods.
Kim \textit{et al}. \cite{kim_cvpr_2016_2} designed a deeply recursive
neural network to raise SR performance without increasing
parameters for additional convolutions.  Ledig \textit{et al}. \cite{ledig_cvpr_2017}
were the first to use GAN for SR, introducing a perceptual
loss function to generate photo-realistic SR images from LR inputs.
Inspired by other SR models using deep residual networks
\cite{kim_cvpr_2016_1} \cite{kim_cvpr_2016_2},
Lim \textit{et al}. \cite{lim_cvprw_2017} simplified the network structure
by removing BN layers and optimized the training process to
achieve the best restoration fidelity at that time.
Zhang \textit{et al}. first applied dense skip-connections \cite{zhang_cvpr_2018}
and later channel attention module \cite{zhang_eccv_2018}
in deep residual network for further advancing of SOTA.
Most recently, Guo \textit{et al}. \cite{guo_cvpr_2020} proposed
a dual-regression method by adding a second downsampling
model and corresponding loss to make sure the restored SR
image can best match the LR input after downsampled
by the co-trained secondary model.

\noindent\textbf{NAS for Single Image Super-Resolution.} As the first attempt to apply NAS for SR,
Chu \textit{et al}. \cite{chu_xx_arxiv_2019} made use of an elastic search method on both micro and macro
level with a hybrid controller that profits from evolutionary computation and reinforcement
learning (RL), achieving comparable performance of PSNR with light model.  Based on different types
of residual blocks and evolutionary algorithm, Song \textit{et al}. \cite{song_dh_aaai_2020} proposed a search
method for better and more efficient network for image SR. Guo \textit{et al}. \cite{guo_y_arxiv_2020} put
forward a novel hierarchical NAS approach that considers both the cell-level and network-level
design based on a RL controller.  While all the three works are promising at searching for
efficient SR models where resources like model size or FLOPS are limited,
they are not able to achieve the high PSNR or SSIM values comparing to
other SOTA methods using manually designed residual networks with
dense skip connections.
Aiming at AIM 2020 challenge that does not take model efficiency in consideration, we
mainly concentrate on the macro level structure design for a better network structure
that can achieve the best SR performance in terms of PSNR and SSIM. Moreover, instead of
using RL or evolutionary based search method that tend to be very time consuming, we apply
GP-NAS approach to search the key network structure parameters such as the number of dense
residual block, the block size and the number of features.

\section{Problem Formulation}
\label{sec:pf}

Learning based image SR methods often rely on a large number of image
pairs, including low-res image $\mathbf{I}_{LR}$ and reference high-res
$\mathbf{I}_{HR}$.  For real image SR, as shown in Eq.~\ref{eq:ds},
$\mathbf{I}_{LR}$ could be modeled from $\mathbf{I}_{HR}$ using three steps:
convolution with a kernel $\mathbf{k}$, downsampling $D_s(\cdot)$
and addition of noise $\mathbf{n}$.

\begin{equation}
     \mathbf{I}_{LR}=D_s(\mathbf{I}_{HR}*\mathbf{k})+\mathbf{n}
\label{eq:ds}
\end{equation}

The goal of image SR is to reverse this process, finding the matching $\mathbf{I}_{HR}$ from a known $\mathbf{I}_{LR}$.
This problem is challenging as there are many versions of
$\mathbf{I}_{HR}$ that could generate the same $\mathbf{I}_{LR}$
following the process in Eq.~\ref{eq:ds}, even when the
kernel $\mathbf{k}$ is known and there is no noise
$\mathbf{n}$.  Learning based SR model $f(\cdot)$ use a total of $n$
image pairs to minimize the average error as in Eq.~\ref{eq:sr}.
It is common that only $\mathbf{I}^i_{HR}$ is a real image
and the corresponding $\mathbf{I}^i_{LR}$ is constructed from
$\mathbf{I}^i_{HR}$ following the process described in Eq.~\ref{eq:ds}.

\begin{equation}
     \arg\min_{f} \sum \|f(\mathbf{I}^i_{LR})-\mathbf{I}^i_{HR} \|,
     i \in \{1, 2, \cdots, n\}
\label{eq:sr}
\end{equation}

In the case that both $\mathbf{I}_{HR}$ and $\mathbf{I}_{LR}$
are real images collected separately, the relationship between the image pair
is much more complicated.  To get a digital image $\mathbf{I}$ from
an object $\mathbf{O}$, there are three transformations included in general
as shown in Eq.~\ref{eq:imaging}.  The optical transformation
$\mathbb{O}(\cdot)$ refers to the process of photons reflected
from the object passing through the lens of the camera.  Optical
characteristics of the lens like modulation transfer function (MTF)
and lens settings like aperture are key variables here.  The second
transformation $\mathbb{D}(\cdot)$ refers to the analog to digital
converter (ADC) that turns photons to digital numbers, where
the noise is introduced.  The last one
$\mathbb{I}(\cdot)$ refers to the image signal processor (ISP)
which transforms noisy raw images to end result of sRGB images.
This step is the most complicated of three,
including multiple processes like denoising and color balancing
at both global and local levels.

\begin{equation}
     \mathbf{I} = \mathbb{I}(\mathbb{D}(\mathbb{O}(\mathbf{O})))
\label{eq:imaging}
\end{equation}

As each of the three transformations could be different for HR
and LR images, the relationship between $\mathbf{I}_{LR}$
and $\mathbf{I}_{HR}$ can be illustrated as in
Equations~\ref{eq:lr}-\ref{eq:lrhr}
where they are linked indirectly by the downsampling
of $\mathbf{O}_{HR}$ to $\mathbf{O}_{LR}$ using $D_s(\cdot)$.
With all these added variations, the real image SR becomes
more challenging.  For example, $\mathbf{I}_{HR}$ is
more clear in general compared to $\mathbf{I}_{LR}$.
But for background objects at a further distance,
they could be more blurry in $\mathbf{I}_{HR}$ if its lens
has a smaller f-number which lead to smaller depth-of-field.
The motivation to use heterogeneous model ensemble is based on
the observation that different model could lead to optimization results
biased towards different variation factors
even when using the same set of training image pairs.

\begin{align}
    & \mathbf{I}_{LR} = \mathbb{I}_{LR}(\mathbb{D}_{LR}(\mathbb{O}_{LR}(\mathbf{O}_{LR}))) \label{eq:lr} \\
    & \mathbf{I}_{HR} = \mathbb{I}_{HR}(\mathbb{D}_{HR}(\mathbb{O}_{HR}(\mathbf{O}_{HR}))) \label{eq:hr} \\
    & \mathbf{O}_{LR} = D_s(\mathbf{O}_{HR}) \label{eq:lrhr}
\end{align}

\section{Proposed Real Image SR Method}
\label{sec:method}

The proposed real image SR method using dense residual network, GP-NAS and heterogeneous model ensemble
is explained in this section.  First, the primary dense residual network
and the search space of different hyperparameters are introduced.
Then, method to find heterogeneous models using GP-NAS is explained, followed by
the multi-level ensemble that is used to generate full-size SR images for the AIM 2020 challenge.

\subsection{Dense Residual Network (DRN)}
\label{sec:drn}

 \begin{figure}
 %
 \begin{center}
     \includegraphics[width=\linewidth]{drn-cropped}
 \end{center}
 %
 \caption{The deep dense residual network architecture for image super resolution.}
 \vspace{-10pt}
 \label{fig:drn}
 %
 \end{figure}

The backbone model of the proposed method is
a deep dense residual network originally developed for raw image demosaicking and denoising.
As depicted in Fig.~\ref{fig:drn},
in addition to the shallow feature convolution at the front and the upsampler at the end, the proposed
network consists of a total depth of $D$ dense residual blocks (DRB).  The input convolution layer
converts the $3$-channel LR input to a total of $F$-channel shallow features.
For the middle DRB blocks, each one includes $L$ stages
of double layers of convolution and the outputs of all $L$ stages are concatenated together
before convoluted from $F\times L$ to $F$ channels.  An additional channel-attention layeris
included at the end of each block, similar to RCAN~\cite{zhang_eccv_2018}.  There are two types of skip connections
included in each block, the block skip connection (BSC) and inter-block skip connection (IBSC).
The BSC is the shortcut between input and output of block $B_i$, while the IBSC includes
two shortcuts from the input of block $B_{i-1}$ to the two stages inside block $B_i$ respectively.  The
various skip connections, especially IBSC, are included to combine features with a large range of receptive fields.
The last block is an enhanced upsampler that transforms all $F$-channel LR features
to the estimated $3$-channel SR image.
This dense residual network has three main hyperparameters:
$F$ is the number of feature channels, $D$ is the number of DRB layers
and $L$ is the number of stages for each DRB. All these three hyperparameters will greatly affect
the performance of SR.  Previous efforts mainly use professional expertise or experience to choose
them based on, which is laborious. To overcome this issue, we apply NAS to search for
the optimal network structure, which will be elaborated in the subsequent subsection.

\subsection{Gaussian Process based Neural Architecture Search}
\label{sec:nas}

Since most NAS methods are still time consuming, we had proposed Gaussian Process based Neural
Architecture Search (GP-NAS) \cite{li_cvpr_2020} to accelerate the searching process.
Figure \ref{fig:gp-nas} illustrates the framework of the GP-NAS. The GP-NAS formulates NAS
from a Bayesian perspective. Specifically, given the hyper-parameters of GP-NAS, we are capable of
predicting the performance of any architectures in the search space effectively. Then, the NAS
process is converted to hyperparameters estimation. By mutual information maximization, we can
efficiently sample networks. Accordingly, based on the performances of sampled networks,
the posterior distribution of hyperparameters can be gradually and efficiently updated.
Based on the estimated hyperparameters, the architecture with best performance can be obtained.
More details about our proposed GP-NAS can be found in \cite{li_cvpr_2020}.

 \begin{figure}[ht!]
 %
 \begin{center}
     \includegraphics[width=\linewidth]{GP_NAS.png}
 \end{center}
 %
 \caption{The framework of the GP-NAS.}
 \vspace{-10pt}
 \label{fig:gp-nas}
 %
 \end{figure}

\subsection{Multi-Level Ensemble}
\label{sec:mle}

Targeting for the AIM 2020 challenge, where the test images are much larger than the training patches, a
multi-levelensemble scheme is designed to achieve optimal image restoration quality. First, for model-ensemble,
the input LR image is processed by a suite of heterogeneous models separately
and the output HR images are averaged to get the final output.  Additionally, each full size LR input is
cropped to patches, with each has an overlapping buffer with neighbouring patches.  A patch-ensemble method
is then used to blend all restored HR patches together, using different weights for each pixel which are
correlated to the distance between the patch center and corresponding pixels.  The most commonly used
self-ensemble is also applied by flipping and/or transposing the input patch before restoration.

\begin{table}[t!]
 \small
 \begin{center}
  \begin{tabular}{rcccccc}
  \toprule
      {} & \multicolumn{2}{c}{Upscaling $\times 2$} & \multicolumn{2}{c}{Upscaling $\times 3$} & \multicolumn{2}{c}{Upscaling $\times 4$} \\
      {} & {PSNR} & {SSIM} & {PSNR} & {SSIM} & {PSNR} & {SSIM} \\
 \midrule
 DRN\textsuperscript{*} & 32.51\textsuperscript{1} & 0.9209\textsuperscript{1} & 31.07\textsuperscript{1} & 0.8796\textsuperscript{1} & 30.26\textsuperscript{3} & 0.                     8401\textsuperscript{3} \\
 RCAN\textsuperscript{*} & 32.31\textsuperscript{3} & 0.9188\textsuperscript{3} & 30.97\textsuperscript{3} & 0.8787\textsuperscript{3} & 30.31\textsuperscript{2} & 0.                    8403\textsuperscript{1} \\
 RCAN & 32.40\textsuperscript{2} & 0.9200\textsuperscript{2} & 31.01\textsuperscript{2} & 0.8792\textsuperscript{2} & 30.32\textsuperscript{1} & 0.8402\textsuperscript{1} \\
 \midrule
 DRN\textsuperscript{*}+RCAN\textsuperscript{*} & 32.56\textsuperscript{2} & 0.9210\textsuperscript{2} & 31.21\textsuperscript{2} & 0.8811\textsuperscript{3} & 30.47\textsuperscript{3}  & 0.8423\textsuperscript{1} \\
 RCAN\textsuperscript{*}+RCAN & 32.49\textsuperscript{3} & 0.9204\textsuperscript{3} & 31.17\textsuperscript{3} & 0.8807\textsuperscript{3} & 30.48\textsuperscript{1} & 0.8421\textsuperscript{3} \\
 DRN\textsuperscript{*}+RCAN & 32.62\textsuperscript{1} & 0.9215\textsuperscript{1} & 31.24\textsuperscript{1} & 0.8814\textsuperscript{1} & 30.48\textsuperscript{1} & 0.8423\textsuperscript{1} \\
 \midrule
 3-Model Ensemble & \textbf{32.63} & \textbf{0.9218} & \textbf{31.28} & \textbf{0.8822} & \textbf{30.55} & \textbf{0.8435} \\
  \bottomrule
  \end{tabular}
 \end{center}
      \caption{Quantitative results of single-model and model-ensemble methods. Best results are in bold and the ranks in each category are superscripted.}
      \vspace{-10pt}
  \label{tab:as}
 \end{table}

\section{Experimental Results}

The AIM 2020 challenge aims to find a generic model to super-resolve LR images
captured in practical scenarios.  To achieve this goal, paired LR and HR images were
taken by various DSLR cameras.  However, images used for training, validation and testing are
captured in the same way with the same set of cameras, so the transformation processes described in
Equations~\ref{eq:lr}-\ref{eq:lrhr} are not changed among images
of the same upscaling factor, meaning the learned transformation from training data
is expected to achieve similar results on the validation and test images.
For this new real image SR dataset,
a total of $19,000$ LR-HR pairs are available for model training for each
of the $\times 2$, $\times 3$ and $\times 4$ upscaling factors.
The LR image resolutions are $380\times 380$ for $\times 3$, $272\times 272$ for
$\times 2$ and $194\times 194$ for $\times 4$ respectively.

For our experiments of each upscaling factor, $600$ of the $19,000$ pairs
are reserved for validation while the remaining ones are used for training.
Note that any LR-HR pairs that are not perfectly aligned,
those with normalized cross-correlation (NCC) less than 0.99,
were excluded from both training and validation.
For each epoch, a $120\times 120$ patch is randomly cropped and augmented
with flipping and transposing from each training image.  A mixed loss of $L1$
and multi-scale structural similarity (MS-SSIM) is taken for training.
For the experiment, the new model candidate search scheme using
GP-NAS was implemented in PaddlePaddle~\cite{ma_fdd_2019}
and the final-training of searched models were conducted using
PyTorch~\cite{paszke_nips_2019}

\begin{figure}[hb!]
\captionsetup[subfigure]{labelformat=empty}
\begin{center}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_395_LR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_395_SR1}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_395_SR2}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_395_SR3}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_395_ESR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_395_HR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_455_LR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_455_SR1}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_455_SR2}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_455_SR3}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_455_ESR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_455_HR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_109_LR}
      \caption{LR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_109_SR1}
      \caption{DRN\textsuperscript{*}}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_109_SR2}
      \caption{RCAN\textsuperscript{*}}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_109_SR3}
      \caption{RCAN}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_109_ESR}
      \caption{Ensemble}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x2_109_HR}
      \caption{HR}
  \end{subfigure}
\end{center}
    \vspace{-10pt}
    \caption{Visual and quantitative comparison of $\times 2$ SR results.}
    \vspace{-10pt}
\label{fig:x2}
\end{figure}

\begin{figure}[ht!]
\captionsetup[subfigure]{labelformat=empty}
\begin{center}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_46_LR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_46_SR1}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_46_SR2}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_46_SR3}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_46_ESR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_46_HR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_88_LR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_88_SR1}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_88_SR2}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_88_SR3}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_88_ESR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_88_HR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_400_LR}
      \caption{LR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_400_SR1}
      \caption{DRN\textsuperscript{*}}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_400_SR2}
      \caption{RCAN\textsuperscript{*}}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_400_SR3}
      \caption{RCAN}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_400_ESR}
      \caption{Ensemble}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x3_400_HR}
      \caption{HR}
  \end{subfigure}
\end{center}
    \vspace{-10pt}
    \caption{Visual and quantitative  comparison of $\times 3$ SR results.}
    \vspace{-10pt}
\label{fig:x3}
\end{figure}

\begin{figure}[ht!]
\captionsetup[subfigure]{labelformat=empty}
\begin{center}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_107_LR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_107_SR1}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_107_SR2}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_107_SR3}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_107_ESR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_107_HR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_173_LR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_173_SR1}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_173_SR2}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_173_SR3}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_173_ESR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_173_HR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_35_LR}
      \caption{LR}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_35_SR1}
      \caption{DRN\textsuperscript{*}}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_35_SR2}
      \caption{RCAN\textsuperscript{*}}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_35_SR3}
      \caption{RCAN}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_35_ESR}
      \caption{Ensemble}
  \end{subfigure}
  \begin{subfigure}[b]{0.16\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=false, clip=true]{x4_35_HR}
      \caption{HR}
  \end{subfigure}
\end{center}
    \vspace{-10pt}
    \caption{Visual and quantitative comparison of $\times 4$ SR results.}
    \vspace{-10pt}
\label{fig:x4}
\end{figure}

\subsection{Ablation Study}
\label{sec:as}

Ablation studies were conducted to examine the effectiveness of
the proposed heterogeneous model ensemble. Three models
were trained for all three upscaling factors.
DRN\textsuperscript{*} is the selected DRN model using GP-NAS
where $F=128, D=18, L=3$.  RCAN\textsuperscript{*} is the GP-NAS
selected variation of RCAN~\cite{zhang_eccv_2018} with 128 features,
5 residual groups and 10 residual blocks in each group.
And RCAN uses the original settings in~\cite{zhang_eccv_2018}
with 64 features, 10 residual groups and 20 residual blocks per group.

As shown in Table~\ref{tab:as}, the average performances of three
individual models are very close.  In general, DRN\textsuperscript{*}
is the best with RCAN as the close second.  In comparison, results from
model ensemble is always better than the individual models of the ensemble
no matter it is a 2-model ensemble or 3-model one.
Overall, the best results are always from the 3-model ensemble.
One thing interesting to note is, out of three 2-model ensemble results,
the one combining RCAN and RCAN\textsuperscript{*} is consistently
worse than others.  It could be explained that, for
heterogeneous model ensemble, differences in network architecture
is more beneficial than differences in hyperparameters of
the same network architecture.

Some image examples are shown in Figures~\ref{fig:x2}-\ref{fig:x4},
where PSNR values of SR results are also annotated for quantitative
comparison too.  For small upscaling factor $\times 2$ as in
Fig.\ref{fig:x2}, there is no big difference among SR
results visually.  In the top example, individual SR models
are all able to super-resolve the fine line features which are
blurred in the LR input.  By combining the models together,
the ensemble result has higher PSNR than individual ones.
But the other two examples show that the ensemble PSNR is lower
than the best individual one.  For the middle one, RCAN\textsuperscript{*}
seems to be the outlier comparing with the other two.
Removing outlier before averaging would have increased
the performance of ensemble.  However, the outlier DRN\textsuperscript{*}
in the bottom example is the most accurate and removing it would lead
to worse ensemble result.

Visual difference starts to stand out more in $\times 3$
examples as represented in Fig.\ref{fig:x3}.  For the
first example on the top, all three
individual models are able to super-resolve the small
dot pattern which are blurred in LR input, but each is
slightly different in details.  By combining the models together,
the ensemble result has higher PSNR than individual ones.
The middle example also show ensemble result is higher
than all individual ones and the visual difference
are around the feather-like structures, showing different color
and contrast.  The bottom example shows DRN\textsuperscript{*}
has sharper image than the other two and its PSNR is higher
than the ensemble as a result.

More interesting visual differences are observed in $\times 4$
examples as shown in Fig.\ref{fig:x4}.  For the first two examples,
individual models resolve the blurred details differently,
like the lines near right bottom right corner in the first
and the small text in the second.  While the ensemble
result doesn't increase the image clarity, it has higher
PSNR than all individual ones.  The last one is worth noting
as it demonstrates the difference in the optical transformation
$\mathbb{O}(\cdot)$ for different cameras as explained in
Section~\ref{sec:pf}.  The front object in HR image is clearer
comparing to LR due to higher spatial resolution, but its background
plant is more blurry than the LR counterpart, probably because
the camera used for HR images has a smaller depth-of-field.
The individual SR models handle this quite differently, with
DRN\textsuperscript{*} mimics the depth-of-field of the LR camera
while the other two closer to the HR version. The ensemble method
results in a significant increase in PSNR comparing to individual ones.

\begin{table}[ht!]
\small
\begin{center}
 \begin{tabular}{rcccccc}
 \toprule
     {} & \multicolumn{2}{c}{Upscaling $\times 4$} & \multicolumn{2}{c}{Upscaling $\times 3$} & \multicolumn{2}{c}{Upscaling $\times 2$} \\
     {} & {PSNR} & {SSIM} & {PSNR} & {PSNR} & {PSNR} & {SSIM} \\
 \midrule
 \textbf{Baidu (ours)} & \textbf{31.3960\textsuperscript{1}} & \textbf{0.8751\textsuperscript{1}} & \textbf{30.9496\textsuperscript{1}} & \textbf{0.8762\textsuperscript{1}} & \textbf{33.4460\textsuperscript{1}} & \textbf{0.9270\textsuperscript{1}} \\
 ALONG & 31.2369\textsuperscript{2} & 0.8742\textsuperscript{3} & 30.3745\textsuperscript{5} & 0.8661\textsuperscript{6} & 33.0982\textsuperscript{8} & 0.9238\textsuperscript{6} \\
 CETC-CSKT & 31.1226\textsuperscript{4} & 0.8744\textsuperscript{2} & 30.7651\textsuperscript{2} & 0.8714\textsuperscript{2} & 33.3140\textsuperscript{2} & 0.9245\textsuperscript{2} \\
 SR-IM & 31.1735\textsuperscript{3} & 0.8728\textsuperscript{7} & - & - & - & - \\
 DeepBlueAI & 30.9638\textsuperscript{6} & 0.8737\textsuperscript{4} & 30.3017\textsuperscript{7} & 0.8665\textsuperscript{4} & 33.1771\textsuperscript{7} & 0.9236\textsuperscript{7} \\
 JNSR & 30.9988\textsuperscript{5} & 0.8722\textsuperscript{8} & - & - & - & - \\
 OPPO\_CAMERA & 30.8603\textsuperscript{8} & 0.8736\textsuperscript{5} & 30.5373\textsuperscript{4} & 0.8695\textsuperscript{3} & 33.3091\textsuperscript{3} & 0.9242\textsuperscript{4}  \\
 Kailos & 30.8659\textsuperscript{7} & 0.8734\textsuperscript{6} & 30.1303\textsuperscript{8} & 0.8664\textsuperscript{5} & 32.7084\textsuperscript{12} & 0.9196\textsuperscript{11} \\
 SR\_DL & 30.6045\textsuperscript{9} & 0.8660\textsuperscript{12} & - & - & - & - \\
 Noah\_TerminalVision & 30.5870\textsuperscript{10} & 0.8662\textsuperscript{11} & 30.5641\textsuperscript{3} & 0.8661\textsuperscript{7} & 33.2888\textsuperscript{4} & 0.9228\textsuperscript{8} \\
 Webbzhou & 30.4174\textsuperscript{12} & 0.8673\textsuperscript{10} & - & - & - & - \\
 TeamInception & 30.3465\textsuperscript{13} & 0.8681\textsuperscript{9} & - & - & 33.2322\textsuperscript{6} & 0.9240\textsuperscript{5} \\
 lyl & 30.3191\textsuperscript{14} & 0.8655\textsuperscript{13} & 30.3654\textsuperscript{6} & 0.8642\textsuperscript{8} & 32.9368\textsuperscript{10} & 0.9210\textsuperscript{9} \\
 MCML-Yonsei & 30.4201\textsuperscript{11} & 0.8637\textsuperscript{15} & - & - & 32.9032\textsuperscript{11} & 0.9186\textsuperscript{12} \\
 MoonCloud & 30.2827\textsuperscript{15} & 0.8644\textsuperscript{14} & - & - & - & - \\
 qwq & 29.5878\textsuperscript{17} & 0.8547\textsuperscript{16} & 29.2656\textsuperscript{9} & 0.8521\textsuperscript{9} & 31.64\textsuperscript{13} & 0.9126\textsuperscript{13} \\
 SrDance & 29.5952\textsuperscript{16} & 0.8523\textsuperscript{17} & - & - & - & - \\
 MLP\_SR & 28.6185\textsuperscript{18} & 0.8314\textsuperscript{18} & - & - & - & - \\
 RRDN\_IITKGP & 27.9708\textsuperscript{19} & 0.8085\textsuperscript{20} & - & - & 29.8506\textsuperscript{14} & 0.8453\textsuperscript{14} \\
 congxiaofeng & 26.3915\textsuperscript{20} & 0.8258\textsuperscript{19} & - & - & - & - \\
 AiAiR & - & - & 18.1903\textsuperscript{10} & 0.8245\textsuperscript{10} & 33.2633\textsuperscript{5} & 0.9243\textsuperscript{3} \\
 GDUT-SL & - & - & - & - & 32.9725\textsuperscript{9} & 0.9204\textsuperscript{10} \\
 \bottomrule
 \end{tabular}
\end{center}
     \caption{Quantitative results of all three tracks for the AIM 2020 Real Image Super-Resolution Challenge. The superscript number indicates ranking of each metric.}
     \vspace{-10pt}
 \label{tab:aimsr}
\end{table}

\subsection{AIM 2020 Challenge Results}
\label{sec:aimc}

To generate the full-size SR images for the AIM 2020 challenge, all three
models from the ablation study were used for all three upscaling tracks.
For $\times 4$ track, a double regression model~\cite{guo_cvpr_2020} was also trained
to include in the model ensemble.  Each full-size LR test images were
cropped to $120\times 120$ patches and self-ensemble ($\times 8$) was applied.
The cropping window was slided at 60-pixel spacing and the overlapping
patches were average using weights correlated with the distance between
the patch center and each pixel.  With ensemble applied to all three levels,
the generated full-size images were submitted to the challenge and won
the first place at all three tracks.  As shown in Table~\ref{tab:aimsr},
both our PSNR and SSIM values lead the second place with a comfortable margin.

One set of representative $\times 4$ images are shown in the
Fig.~\ref{fig:fs_x4}.  The full-size SR output are located at
the top, with selected areas zoomed in to compare with
bicubic interpolation results.

\section{Conclusions}
In this paper, based on the models searched via GP-NAS, we have introduced a new heterogeneous model
ensemble method for real image super resolution. Since network architecture greatly affects the results
of SR, we first apply GP-NAS approach to search the key factors such as the number of residual network
block, block size and the number of features in our network structure.  Then, different models selected
using GP-NAS are fused together to boost the performance of SR.
Combined with patch-ensemble and self-ensemble, the proposed new scheme is validated to be
highly effective, generating impressive testing results on all three
tracks ($\times 2$, $\times 3$ and $\times 4$) of the AIM 2020 challenge in terms of both PSNR
and SSIM.

\begin{figure}[ht!]
\captionsetup[subfigure]{labelformat=empty}
\begin{center}
  \begin{subfigure}[b]{0.48\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=true]{SR_032_A_20}
  \end{subfigure}
  \begin{subfigure}[b]{0.48\textwidth}
    \centering
      \includegraphics[width=\textwidth, interpolate=true]{SR_040_A_20}
  \end{subfigure}
  \begin{subfigure}[b]{0.24\textwidth}
    \centering
      \vspace{-0.25em}
      \includegraphics[width=\textwidth, interpolate=false]{LR_032_1}
  \end{subfigure} \hspace*{-0.5em}
  \begin{subfigure}[b]{0.24\textwidth}
    \centering
      \vspace{-0.25em}
      \includegraphics[width=\textwidth, interpolate=false]{SR_032_1}
  \end{subfigure}
  \begin{subfigure}[b]{0.24\textwidth}
    \centering
      \vspace{-0.25em}
      \includegraphics[width=\textwidth, interpolate=false]{LR_040_1}
  \end{subfigure} \hspace*{-0.5em}
  \begin{subfigure}[b]{0.24\textwidth}
    \centering
      \vspace{-0.25em}
      \includegraphics[width=\textwidth, interpolate=false]{SR_040_1}
  \end{subfigure}
  \begin{subfigure}[b]{0.24\textwidth}
    \centering
      \vspace{-0.25em}
      \includegraphics[width=\textwidth, interpolate=false]{LR_032_2}
      \caption{Bicubic}
  \end{subfigure} \hspace*{-0.5em}
  \begin{subfigure}[b]{0.24\textwidth}
    \centering
      \vspace{-0.25em}
      \includegraphics[width=\textwidth, interpolate=false]{SR_032_2}
      \caption{SR\_Ensaemble}
  \end{subfigure}
  \begin{subfigure}[b]{0.24\textwidth}
    \centering
      \vspace{-0.25em}
      \includegraphics[width=\textwidth, interpolate=false]{LR_040_2}
      \caption{Bicubic}
  \end{subfigure} \hspace*{-0.5em}
  \begin{subfigure}[b]{0.24\textwidth}
    \centering
      \vspace{-0.25em}
      \includegraphics[width=\textwidth, interpolate=false]{SR_040_2}
      \caption{SR\_Ensaemble}
  \end{subfigure}
 \end{center}
    \vspace{-10pt}
    \caption{Examples of full-size test images from the AIM 2020 challenge ($\times 4$).}
    \vspace{-10pt}
\label{fig:fs_x4}
\end{figure}

\bibliographystyle{unsrt}


\section{Introduction}

Electrides are a class of unconventional compounds that contain excess valence electrons confined in the void or interlayer space that play the role of anions \cite{dye2009electrides}. Chemically, such materials could serve as strong reducing agents or catalysts for chemical syntheses. From a physics point of view, the confined electrons form (partially occupied) bands close to Fermi level, which could lead to a dramatically reduced work function and high electric conductivity. These attractive physical and chemical properties collectively promise many technological applications such as the splitting of carbon dioxide at room temperature \cite{Toda-NC-2013}, synthesis of ammonia from atmospheric nitrogen under mild conditions \cite{Kitano-NChem-2012} and many others \cite{Hosono-PNAS-2017}. Recently, it was also proposed that the confined electrons in the electrides are favorable for achieving band inversions needed for topological phase transition in the electronic states \cite{Hirayama-PRX-2018,Huang-NanoL-2018, Park-PRL-2018, Zhu-PRM-2019}.

To date, a number of electride materials have been investigated both experimentally and computationally. Depending on the connectivity of crystal cavities and channels, the identified electrides could be classified to zero, one and two dimensions (0D, 1D and 2D). Of them, the 2D layered electride materials \cite{Lee-Nature-2013, Inoshita-PRX-2014}, with anionic electrons distributed in the 2D space, are of great interests since this is naturally connected with the 2D materials physics \cite{zhao2014Ca2N}. Indeed, the recent report on the synthesis of the monolayer Ca$_2$N (exfoliated from the parent layered form) has driven electride research into the nano-regime \cite{2D-Ca2N}. The monolayer Ca$_2$N, by combining high surface areas of 2D materials with the exotic properties of anionic electrons, are favorable for applications as transparent conductors in the modern electronic devices \cite{2D-Ca2N, zhao2014Ca2N}.

Despite the many unique properties that the floating anionic electrons bring, they also result in poor stabilities at room temperature and ambient atmosphere, which significantly limit the utilization of electrides as materials. In electride design, the conventional wisdom is to search for candidates with unpaired electrons present in their crystal voids \cite{Tada-IC-2014,burton2018high,Zhang-PRX-2017,ZHU20191293}.
In the electronic band structure, these confined electrons form the distinct interstitial bands crossing the Fermi level. To date, most known electrides are strong reducing metals due to the presence of these characteristic bands. Consequently, they are unstable in water and air. From an application perspective, it is desirable to have the material remain stable under normal conditions and switch to the electride phase under the operation. Therefore, an ideal material should be semiconducting with the lowest conduction bands occupied by the interstitial electrons. Under the working condition, the interstitial band can be then downshifted to the Fermi level. As such, the interstitial band can act as an effective transport agent in electronic devices. To achieve this goal, the key is to identify the materials that can hold an interstitial band close to the Fermi level.

In many known electrides, the interstitial electrons follow the nearly free electron (NFE) model \cite{ZHU20191293} in which the electrons are considered as plane waves that are weakly perturbed by the periodic lattice potential. The dispersion can be roughly described by a parabolic-like dispersion of $E=\hbar^2k^2/2m$. While the NFE states can be found in the 3D systems, it is more natural to study them in 2D materials ~\cite{NFE-graphene,chen2012band,Cuong_2014,BL-InSe-E}. To our knowledge, most of the previous reports demonstrate that the NFE states are several eVs above the Fermi level. The high energy NFE states require common band structure engineering approaches, such as chemical/electron doping \cite{NFE-C_Nanotube, NFE-C-BN_nanotube}, strain engineering \cite{ni2012tunable, song2016room, Zhu-2013-MoS2}, and electric fields \cite{min2007ab, mak2009observation, liu2015switching, blueP-E, InSe2015JPCL, BL-InSe-E, MoS2-E, WSe2-E}, to lower their energy levels. It was only recently reported that the NFE states can be energetically found near the Fermi levels for a few 2D transition metal carbides \cite{NFE-PRB}. Based on first-principle calculations, the authors found that Sc$_2$C(OH)$_2$ was semiconducting with its conduction band occupied by the NFE states 0.45-0.94 eV above the Fermi Level \cite{NFE-JMCC}. They further proposed that the NFE states can be downshifted by external electric fields \cite{NFE-JMCC}. However, Sc$_2$C(OH)$_2$ is a hypothetical compound. To our knowledge, such switchable 2D electride phases have not been reported based on any known material.

In this work, we focus on potential 2D materials that possess accessible NFE states from strain engineering. We performed a high-throughput screening of the available 2D semiconducting materials database by searching for characteristic NFE states. From the screening, we identified 12 2D materials promising to hold the low energy NFE states at their conduction bands. Furthermore, we found that applying compressive strain can effectively push the NFE bands to the Fermi level. Different from the previous report on a hypothetical material \cite{NFE-JMCC}, 11 of the identified 2D materials can likely be exfoliated from the available layered compounds without further surface functionalization. Therefore, our predictions can be readily validated by experiment if one of the materials can be synthesized in the laboratory. Below, we will discuss the computational methods used and present our discovery in the following sections.

In total, 636 2D candidates from material databases were investigated based on the four steps of screening as shown in Fig.~\ref{Fig1}. The first step was to obtain potential 2D insulators/semiconductors from the online JARVIS-DFT database \cite{choudhary2017high}. For each 2D material, we performed the density functional theory (DFT) calculations using a projector augmented-wave method~\cite{PAW-PRB-1994} as implemented in the Vienna ab initio Simulation Package (VASP)~\cite{Vasp-PRB-1996}.  The exchange and correlation potentials were treated with the generalized gradient approximation (GGA) as proposed by Perdew–Burk–Ernzerhof (PBE)~\cite{PBE-PRL-1996}. Preparation of input parameters and data extraction was done by the Python Materials Genomics (Pymatgen) package~\cite{ong2013python}. We performed the geometry relaxation (MPRelaxSet), static (MPStaticSet) and non-self-consistent (MPNonSCFSet) calculations for each structure. A plane-wave cut-off energy of 700 eV was employed. The relaxation convergence criteria were 10$^{-6}$ eV for energy and 0.01 eV/\AA~ for force. In the relaxation, static and non-self-consistent calculations, the reciprocal density was set to be 300. To avoid the interlayer interaction due to periodic boundary conditions, a vacuum layer of at least 15 \AA~ thickness was added to each simulation at the screening stage. The second step was to extract the 2D semiconducting materials by applying the band gap restriction (between 0.2-2.0 eV). Steps 3 and 4 involve the identification of the NFE states. We filtered the materials with eigenstates at the $\Gamma$ point ($E_\Gamma$) of the first conduction band that are less than 2.0 eV. To ensure the conduction bands held the true NFE states, we also checked that the majority of the density of states (PDOS) could not be projected to the atomic orbitals, similar to a recent high-throughput screening work on electrides \cite{burton2018high}. After a few rounds of parameter tuning, we decided to use 0.3 as the upper bound of the PDOS value at the $\Gamma$ point as determining criteria for whether or not it is an NFE state. These screenings successfully identified 12 materials as summarized in Table \ref{all}. To check the electronic band structure further clearly, we recalculated these 12 materials with higher reciprocal density.

\begin{figure}[ht]
\centering
\includegraphics[width=0.48 \textwidth]{Fig1-sketch.png}
\caption{\label{Fig1} (a) Computational Scheme for the screening of 2D swichable electrides. (b) The definition of physical quantities used in the screening. The PDOS denotes the ratio of densities which can be projected to the atomic orbitals. If the densities cannot be projected to any atomic orbitals, it suggests the electrons are loosely bound with the nucleus.}
\end{figure}

\begin{table}[ht]
\caption{The summary of 2D materials investigated in this work. $E_\textrm{NFE}$ denotes the lowest gap between the NFE states and Fermi energy. $\epsilon_c$ denotes the critical biaxial compressive strain values to turn the material metallic. The exfoliation energies are extracted from the online database \cite{choudhary2017high}. Since PbS did not show zero gap even when put under up to a 20\% strain, we did not show its $\epsilon_c$ value.
}\label{all}
\begin{tabular}{ccccccc}
\hline\hline
Material             &~~Gap~~  &$E_\textrm{NFE}$ & VBM at $\Gamma$ & ~~~~$\epsilon_c$~~~~   & $E_\textrm{exfoliation}$ \\
                     & ({eV})       & (eV)  &   & &(meV/atom)              \\
\hline
RbLiS                & 1.228   & 1.228    & Yes    & -13\%     & 94.14   \\
RbLiSe               & 1.326   & 1.326    & Yes    & -12\%     & 100.08  \\
KAgSe                & 0.511   & 0.511    & Yes    & -9\%      & 101.49  \\
K$_4$HgAs$_2$        & 0.793.  & 0.793    & Yes    & -12\%     & 107.36  \\
KMgSb                & 0.905.  & 0.905    & Yes    & -9\%      & 110.65  \\
NaZnAs               & 0.480.  & 0.480    & Yes    & -6\%      & 118.68  \\
NaZnP                & 0.719.  & 0.719    & Yes    & -7\%      & 124.80  \\
Na$_2$Pd$_3$O$_4$    & 0.376.  & 0.497    & No     & -3\%      & 139.05  \\
PbS                  & 1.654   & 1.924    & No     & N/A       & 187.54  \\
Na$_3$As             & 0.285.  & 0.285    & Yes    & -6\%      & 250.04   \\
Rb$_2$Te             & 0.477   & 0.477    & Yes    & -7\%      & 388.51  \\
SiS                  & 0.494.  & 0.494    & Yes    & -13\%     & N/A   \\
\hline\hline
\end{tabular}
\end{table}

\begin{figure}[ht]
\centering
\includegraphics[width=0.49 \textwidth]{Fig2-420-Na2Pd3O4-strain-new.png}
\caption{\label{Fig2} The electronic structures of monolayer Na$_2$Pd$_3$O$_4$. (a)-(c) show the band structures and corresponding electron densities at $\Gamma$ point with 0, -3\%, 9\% strains. In the left column, each energy band is colored by the atoms which hold the largest portion of PDOS. If it cannot be projected to any atomic site, this corresponds to NFE states as denoted by black color. The right columns display the partial charge density distributions of the three lowest conduction bands at $\Gamma$ point along the $z$-axis. The order of band number are corresponding to the circled numbers in left columns. The side views of the structural model are shown in (a) as the inset. The Na, Pd and O atoms are denoted by yellow, grey and red spheres respectively. The inset in (b) is the 2D band decomposed charge density distribution of the lowest conduction band at $\Gamma$ point with the strain.}
\end{figure}

From the 12 identified materials, it is interesting to note that most of them have the alkaline metals (Li, Na, K, Rb) sitting at the outmost layer of the slab, as shown in Fig. S1 \cite{SI}. This is expected since their valence electrons are much more weakly bound than those in inner shells. To proceed, we first ruled out the materials with high exfoliation energies. Among them, PbS, Na$_3$As, Rb$_2$Te were discarded. We did not consider SiS either since this compound is from a theoretical prediction \cite{yang2016two}. Then, we applied a series of biaxial strains on the remaining compounds. Below, we will discuss how to modulate the NFE states by strain engineering on two representative materials: Na$_2$Pd$_3$O$_4$ and NaZnAs.

Fig. \ref{Fig2} displays the changes of the electronic band structures of Na$_2$Pd$_3$O$_4$ under a series of lattice strains. The pristine monolayer Na$_2$Pd$_3$O$_4$ has an orthorhombic unit cell. It is a semiconductor with a direct gap of 0.377 eV that has both the valence band maximum (VBM) and conduction band minimum (CBM) at the X point (0.5, 0, 0). Not surprisingly, the system becomes metallic when the strains are applied. However, the behaviors are drastically different under different types of strains. Under tension, the band gap goes down continuously and becomes closed when it reaches 9\% strain. On the other hand, compressive strain leads to a more dramatic change and it needs only -3\% strain to close the band gap. To understand the distinct behavior, we tracked the electronic band evolution and analyzed the characteristic wavefunctions under 0, -3\% and 9\% strains. As shown in Fig. \ref{Fig2}c, the dispersion of several conduction bands around the $\Gamma$ point follows a parabolic-like relation. We also plotted the band decomposed charge density of the three lowest conduction bands (labelled from low to high) at the $\Gamma$ point. Both bands 1 and 2 have four spikes. There are two smaller charge peaks at the oxygen's nuclear sites, suggesting that O atoms can take some electrons since they are more electronnegative than Na and Pd. More importantly, the majority of charges are symmetrically distributed outside the slab, confirming the presence of NFE states. With compression biaxial strain, the NFE states (0.497 eV at band 1 and 0.654 eV at band 2) drop more quickly than other conduction bands. At -3\% strain, the CBM is shifted to the $\Gamma$ point due to the lowered NFE states while the original CBM state at the X point remains intact. When the tensile strains are applied, the semicondutor-to-metal transition follows a completely different path in which one nearly flat valence band moves upward and the conduction bands are collectively pushed downward. Due to these joint efforts, the gap is finally bridged. However, it is notable that the NFE states actually take higher energies during the metallization process. The total evolution has been shown in Fig. S3 \cite{SI}.

\begin{figure}
\centering
\includegraphics[width=0.48 \textwidth]{Fig3-420-Na2Pd3O4-potential-strain.png}
\caption{\label{Fig3} (a) The calculated electrostatic potential distribution of monolayer Na$_2$Pd$_3$O$_4$ with/without biaxial strains. The band decomposed charge density distribution at $\Gamma$ point of the first conduction band under -3 \% biaxial strain is also shown for comparison. The positions of each atom along $z$-axis are shown as inset. The Na, Pd and O atoms are denoted by yellow, grey and red spheres respectively. (b) The enlarged detail of the dashed box in (a).}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.48 \textwidth]{Fig4-438-Na2Zn2As2-strain-new.png}
\caption{\label{Fig4} The electronic structures of monolayer NaZnAs. The band structures under (a) 0, (b) -6\% and (c) -8\% biaxial strains. The right column plots the corresponding partial charge densities of three representative bands (two highest valence bands and one lowest conduction band) at $\Gamma$ point along the $z$-axis. For comparison, the isosurfaces of charge density plots for each band were also shown as the insets. The Na, Zn and As atoms are denoted by yellow, grey and green spheres respectively.}
\end{figure}

The direction dependent behavior can be understood by examining the nature of the NFE states. We first found that the bulk Na$_2$Pd$_3$O$_4$ did not exhibit low-energy NFE bands (Fig. S2) \cite{SI}. Different from the bulk electrides where the NFE states are occupying the crystal voids, the NFE states of 2D electrides stem from the image potential states. When the metal has a band gap near the vacuum level $E_\textrm{vacuum}$, an electron below $E_\textrm{vacuum}$ may be trapped in the potential well formed by the Coulomb-like attractive image potential and the repulsive surface barrier \cite{hofer1997time, NFE-C_Nanotube, NFE-C-BN_nanotube, NFE-graphene}. Whether or not the electron can be trapped depends on the spatial distribution of the image potential. If the image potential has a long tail, the NFE states are more likely to be stabilized. Fig. \ref{Fig3} displays the electrostatic potential (ESP) distribution of Na$_2$Pd$_3$O$_4$ under 0, -3\% and 9\% strains. Clearly, the tensile strain sharpens the ESP distribution compared to the pristine form. Thus, the NFE states under tension occupy higher energy levels compared to those in the pristine form. On the other hand, the compressive strain makes the ESP more extended by creating a shoulder around -1 eV. Consequently, the electrons can be trapped around the shoulder with a lower eigenvalue. We also confirmed the results by numerical simulation. Solving the 1D Schordinger equation for the potential wells given by Fig. \ref{Fig3}, we indeed found that the NFE states under the negative strain yielded lower eigenvalues as compared to those under the zero strain (Fig. S4-S5) \cite{SI}. Therefore, we conclude that a negative biaxial strain can promote the activity of the NFE states by modifying the ESP distribution.

Different from Na$_2$Pd$_3$O$_4$, many other materials in Table \ref{all} have their VBM at the $\Gamma$ point. For these materials, the semiconductor-metal transitions under strain become more complicated. As shown in Fig. \ref{Fig4}, the pristine monolayer NaZnAs has its VBM at the $\Gamma$ point. This VBM corresponds to the As's $p_x$/$p_y$ orbitals (bands 1 and 2, see more details in Fig. S6 \cite{SI}). The CBM at the $\Gamma$ point of band 3 is 0.476 eV above the VBM. The decomposed charge density plot suggests that four competing maxima exist in band 3 with two corresponding to the NFE states and the other two corresponding to the electrons trapped by the As nucleus. Compared to band 1 of the Na$_2$Pd$_3$O$_4$ in Fig. \ref{Fig2}a, the dispersion of band 3 in Fig. \ref{Fig4}a is less parabolic since it has more densities inside the slab. Clearly, the hybridization between As's $p_x$/$p_y$ orbitals (bands 1 and 2) and the interstitial electrons dominate the electronic properties near the Fermi level. When put under compressive strain, the energy of NFE states got closer to the Fermi level. Consequently, there is a prominent charge transfer from As's $p_x$/$p_y$ orbitals to the interstitial electrons. When it reaches a critical strain near -6\% (Fig. \ref{Fig4}b), the charge transfer is complete and the hybridization effect becomes negligible. Therefore, further applying the strain leads to a quick downshift of the interstitial electrons (band 1 in Fig. \ref{Fig4}c). For instance, at -8\% strain, the energy of interstitial electrons was reduced to -0.362 eV, while the Fermi level is occupied by As's $p_x$/$p_y$ orbitals (bands 2 and 3).
These phenomena suggest that these 2D materials with low-energy NFE states can exhibit very rich electronic behavior by altering their mechanical attributes. In particular, the coexistence of flat band and NFE band (Fig. \ref{Fig4}b) may be useful for exploring novel quantum phenomena.

We also checked the effects of strain on the previously reported Sc$_2$C(OH)$_2$ (Fig. S7) \cite{SI}. It showed a similar semiconductor-metal transition at the -6\% biaxial strain, suggesting that it's easier to trigger the electride phase transition in Na$_2$Pd$_3$O$_4$ than in Sc$_2$C(OH)$_2$. From Fig. \ref{Fig5}, it is clear that the stress-strain relations follow a linear behavior up to the -8\% strain. The -3\% compressive strain falls into the region of elastic deformation. Therefore, the transition between the electride and non-electride phase of Na$_2$Pd$_3$O$_4$ is likely reversible, which is advantageous for switchable device applications. We also systematically investigated other materials in Table \ref{all}. Among them, 5 materials need a critical strain less than 10\% to complete the semiconducting to metal transition. Introducing strains to 2D nanomaterials has been widely applied to tune the material’s properties \cite{deng2018strain}. These materials require energies between 90-130 meV/atom to be exfoliated from their parent materials, which is higher than graphene (70 meV/atom) and MoS$_2$ (76 meV/atom) but comparable to PtSe$_2$ (111 meV/atom) \cite{wang2015monolayer, choudhary2017high}. To our knowledge, none of the materials in Table \ref{all} have been experimentally synthesized. However, several of them have been theoretically investigated \cite{yang2016two, wang2018kagse, zhang2020high}. If one of the 2D forms can be made, their electric measurements should be straightforward.

\begin{figure}
\centering
\includegraphics[width=0.35 \textwidth]{Fig5-NaZnAs.png}
\caption{\label{Fig5} Calculated stress–strain curves of monolayer Na$_2$Pd$_3$O$_4$. The inset is the schematic diagram of biaxial compressive strain. The solid points denote the critical biaxial compressive strain values $\epsilon_c$ to complete the semiconductor-metal transition.}
\end{figure}

In sum, we report the computational design of 2D switchable atomically thin electrides from first-principles calculations in conjunction with high-throughput screening. We have investigated over 600 potential 2D materials which are potentially exfoliatable from the existing layered materials. Among them, we found that a family of 2D materials may hold the NFE states that are close to the Fermi level. Different from the conventional bulk electrides or the recently synthesized 2D monolayer intrincic electride Ca$_2$N, these monolayer electrides are more chemically inert semiconductors at normal conditions. Under a strain less than 10\%, their NFE states can be lowered and form the partially occupied bands crossing the Fermi level. Consequently, the electride phase transition can be reversible via strain engineering, which provides an ideal platform for switchable nano device/sensors for gas detection and electron transports. Compared to recent studies \cite{NFE-JMCC}, the parent materials discussed in this work have already been synthesized, and several of them require less critical strain to activate the electride phase. In addition to strain manipulation, other band structure engineering methods, such as applying the electric field and chemical doping, are also expected to work given that most of these materials have NFE states close to the Fermi level. Thanks to the improved stability, these predicted materials are more amenable for materials fabrication and operation. We hope this work can stimulate further experimental investigation.

QZ and XHY acknowledge the use of computing resources from XSEDE (DMR180040). This work was also supported by the National Natural Science Foundation of China (Grant Nos. 51872048, 21973012). The study of XHY at University of Nevada is supported by China Scholarship Council (CSC) No.201806650022.


\section{\boldmath INTRODUCTION}
The hyperons, as the $SU(3)$-flavour-octet partners of the nucleons that contain one or more strange quarks,
offer crucial additional dimensions to the study of baryon structures~\cite{hyperon1,hyperon2}.
Treating the heavier strange quarks as spectators, hyperons can provide valuable insight into the behaviour of the lighter up and down quarks in different environments.
Electromagnetic form factors~(EMFFs) are fundamental observables of baryons
that are intimately related to their internal structure and dynamics~\cite{chiral, pqcd, lattice1}.
Despite the fact that much work has been done on the EM structures of protons in both the space-like and
time-like regions~\cite{space_ros, space_jlab, time_cleo, time_babar, time_cmd, time_bes3},
experimental information regarding the EMFFs of hyperons remains limited~\cite{hyperon_babar, hyperon_cleo, hyperon_bes3, hyperon_bes32}.
Moreover, the few existing measurements of time-like neutron FFs~\cite{nnbar_fenice, nnbar_bes3} differ from
each other and lead to conflicting conclusions when compared to those for the proton~\cite{ratio_theory1, ratio_theory2}.
A $\Sigma^+$ hyperon is formed by replacing the proton's down quark with a strange quark; likewise a
$\Sigma^-$ is formed by replacing  the neutron's up quark with a strange quark.
Thus the corresponding ratio of FFs between  the $\Sigma^{+}$ and $\Sigma^{-}$ hyperons could provide
guidance for the nucleons.
The $\Sigma$ hyperon FFs have been predicted by lattice QCD~\cite{lattice2} and in different theoretical models, such as
perturbative QCD~(pQCD)~\cite{ratio_theory2}, the Unitary and Analytic~(U\&A) model~\cite{UandA}, chiral perturbative theory~(ChPT)~\cite{ChPT}, and the di-quark correlation model~\cite{diquark}.
Experimental measurements for $\Sigma$ hyperons,
especially the $\Sigma^{-}$, which has never been previously measured in the time-like region, provide essential
tests of these models and produce important input for the understanding of baryonic structures.

Due to their short lifetimes, studies of hyperon EMFFs are mostly performed in the time-like region via electron-positron annihilation experiments.
The differential, one-photon exchange cross section for the $e^{+}e^{-}\rightarrow B\bar{B}$ process,
where $B$ is a spin-1/2 baryon, can be expressed in terms of the electric and magnetic FFs $G_{E}$ and $G_{M}$
as~\cite{ee}:
\begin{equation}
\label{equ:eq1}
\frac{d\sigma^{B}(s)}{d\Omega} = \frac{\alpha^{2}\beta C|G_{M}|^{2}}{4s}\left[(1+\cos^{2}\theta) +
  \frac{1}{\tau}|\frac{G_{E}}{G_{M}}|^{2}\sin^{2}\theta\right],
\end{equation}
where $\alpha$ is the fine-structure constant,
$s$ is the square of center-of-mass~(c.m.)~energy,
$\beta=\sqrt{1-4m_{B}^{2}/s}$ is a phase-space factor,
$\tau=\frac{s}{4m_{B}^{2}}$,
$m_{B}$ is the baryon mass, and $\theta$ is its c.m.~production angle.
The Coulomb correction factor $C$~\cite{schwinger,coulomb}
accounts for the electromagnetic interaction of charged point-like fermion pairs in the final state.
It reads  $C=y/(1-e^{-y})$ with $y=\pi\alpha(1+\beta^{2})/\beta$
for a charged point-like fermion pair and $C=1$ for a neutral point-like fermion pair.
For charged point-like fermion pairs, the cross section at threshold is non-zero,
$\sigma(4m^{2}_B)=\pi^2\alpha^3/2m^2_B = 848(m_p/m_B)^2$~pb, where $m_{p}$ is the proton mass~\cite{pdg},
and then grows with increasing $\beta$.
Experimentally, a rapid rise of the $e^{+}e^{-}\rightarrow p\bar{p}$ cross section near threshold
followed by a plateau is observed~\cite{time_babar, time_cmd}.
The at-threshold cross section is consistent with the 848~pb expectation for a point-like charged particle.
However, in this case, the $p\bar{p}$ is produced by a virtual photon with $Q^2=4m_p^2=3.53$~GeV$^2$,
which corresponds to a Compton wavelength of $\sim$0.1~fm,
a scale at which the proton is definitely not point-like.
A similar feature of the cross section for $e^{+}e^{-}\to \Lambda_{c}^{+}\bar{\Lambda}_{c}^{-}$ is observed by the BESIII experiment~\cite{lambdac_bes3},
where the at-threshold cross section is around 240~pb. This is
1.6 times the predicted value for point-like charged particles.
These unexpected threshold effects have been widely discussed in the literature
where they are interpreted as
final state interactions~\cite{int_fsi}, bound states or near-threshold meson resonances~\cite{int_res}, or an attractive Coulomb interaction~\cite{int_coulomb}.
To understand the nature of these threshold effects, experimental measurements of the near threshold charged pair production of other hyperons are expected to be of critical importance.

In this Letter, we present studies of $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$
and $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$ with a data sample of 329.7~pb$^{-1}$
collected at BESIII with c.m.~energies between $2.3864$ and 3.0200~GeV~\cite{dataset}.
The threshold energies for $\Sigma^{+}\bar{\Sigma}^{-}$ and $\Sigma^{-}\bar{\Sigma}^{+}$ pair production are
$2.3787$~GeV and 2.3949~GeV, respectively.
The BESIII detector is described in detail in Ref.~\cite{Ablikim:2009aa}. The critical
elements for the measurements reported here are: the main drift chamber (MDC), which
measures the momenta of charged particles with 0.5\% resolution for 1~GeV/$c$ tracks and
the $dE/dx$ for charged-particle identification (PID); a barrel array of scintillation counters that
measures charged particles' time of flight for additional PID information; and an electromagnetic
calorimeter (EMC) comprising an array of CsI(Tl) crystals that measures photon energies with a
resolution of 2.5\% at 1 GeV.

Simulated event samples produced with a {\sc geant4}-based~\cite{geant4} Monte Carlo~(MC)
package that includes the geometric description of the BESIII detector and its response,
are used to determine the detection efficiency and to estimate the backgrounds.
The signal processes $e^{+}e^{-}\to\Sigma^{\pm}\bar{\Sigma}^{\mp}$ are generated according to
the differential amplitude presented in Ref.~\cite{Anjrej}.
Initial state radiation~(ISR) is simulated with {\sc conexc}~\cite{conexc} and the corresponding correction factors
are calculated for higher order processes. Background from the QED processes $e^{+}e^{-}\to l^{+}l^{-}~(l=e,\mu)$ and $e^{+}e^{-}\to\gamma\gamma$  are investigated  with {\sc babayaga}~\cite{babayaga},
while for $e^{+}e^{-}\to $hadrons and two-photon processes we use
{\sc lundarlw}~\cite{lundarlw} and {\sc bestwogam}~\cite{bestwogam}, respectively.
\begin{figure}[htbp]
\vspace{-0.3cm}
\begin{overpic}[width=8.cm,height=3.5cm,angle=0]{decayv6}
\put(-2,30){\small{a)}}
\put(17.5,20.2){\small{\mycircled{1}}}
\put(17.5,13.5){\small{\mycircled{2}}}
\put(48,30){\small{b)}}
\put(32, 20.2){\small{\mycircled{3}}}
\put(32, 13.5){\small{\mycircled{4}}}
\end{overpic}
\vspace{-1.2cm}
\caption{ a) Decay chain of $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$, with four final state topologies,
and b) Decay chain of $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$ with one final state topology.
}
\label{decaychain}
\end{figure}
The decay chains of $e^{+}e^{-}\to\Sigma^{\pm}\bar{\Sigma}^{\mp}$ are shown in Fig.~\ref{decaychain},
where two decay channels for  $\Sigma^{+}$ (or $\bar{\Sigma}^{-}$) and one channel for
$\Sigma^{-}$ (or $\bar{\Sigma}^{+}$) decay account for more than 99\% of their total decay widths.

\begin{figure*}[htbp]
\begin{center}
\begin{overpic}[width=8.9cm,height=3.5cm,angle=0]{fit_23864v6_50bin}
\put(14,33){\small{a)}}
\put(61,33){\small{b)}}
\put(14,24){\tiny{$\sqrt{s}=2.3864$~GeV}}
\put(63,24){\tiny{$\sqrt{s}=2.3864$~GeV}}
\end{overpic}
\begin{overpic}[width=8.9cm,height=3.5cm,angle=0]{fit_2396v6_50bin}
\put(14,34){\small{c)}}
\put(61,34){\small{d)}}
\put(14,24){\tiny{$\sqrt{s}=2.3960$~GeV}}
\put(63,24){\tiny{$\sqrt{s}=2.3960$~GeV}}
\end{overpic}
\caption{ (color online)
The mass spectra of $M_{\text{bc}}$ (category A) and $M^{\rm rec}_{\bar{n}\pi^{-}}$ (category B)
for $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$ candidate events
at a,b) $\sqrt{s}=2.3864$~GeV and c,d) $\sqrt{s}=2.396$~GeV.
Dots with error bars are the data; histograms are the background events in
MC samples after normalization. Solid curves are the fit results,
dashed curves are the signal, and  dot-dashed curves are the background.}
\vspace{-0.5cm}
\label{stacksigmap}
\end{center}
\end{figure*}

In the process $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$,  there are four final state topologies:
\noindent
\mycircled{1}\mycircled{3}, \mycircled{1}\mycircled{4}, \mycircled{2}\mycircled{3},
and \mycircled{2}\mycircled{4} as shown in Fig.~\ref{decaychain}.
All four configurations are selected in this analysis, significantly improving the statistics.
At BESIII, charged particles are efficiently detected and identified by the MDC and PID systems and
$\pi^{0}$ mesons are reconstructed in the EMC via their $\pi^0\to\gamma\gamma$ decay mode.
The selection criteria for charged tracks, PID, and photon candidates are the same as those used
in Ref.~\cite{pinglamlam}.
Most of the anti-neutrons ($\bar{n}$) annihilate in the EMC and produce several
secondary particles with a total energy deposition that can be as high as 2~GeV;
the position of the $\bar{n}$ interaction and, from this,
the $\bar{n}$ direction can be inferred from the weighted center-of-energy of the shower.
Neutron~($n$)  detection is not done because of its low interaction efficiency and small energy deposition.

The \mycircled{1}\mycircled{3} and \mycircled{2}\mycircled{3} final-state configurations, classified
as category A, can be analyzed by a partial reconstruction technique in which only the detection of
$\bar{\Sigma}^{-}\to\bar{p}\pi^{0}$ is required.
Candidate events are required to have at least one charged track that is identified as a $\bar{p}$ by the PID
system  and at least two good photons that are consistent with originating from $\pi^0\to\gamma\gamma$.
The mass spectrum of $\gamma\gamma$ is required to be from $0.127<M_{\gamma\gamma}<0.139$~GeV/c$^{2}$
to $0.123<M_{\gamma\gamma}<0.14$~GeV/c$^{2}$, depending on c.m.~energies. The $\bar{\Sigma}^{-}$ is reconstructed using all combinations of the
selected $\bar{p}\gamma\gamma$.
The two-body process exploits two variables that are based on energy and momentum conservation: the energy difference
$\Delta E\equiv E-E_{\rm beam}$ and the beam-constrained mass $M_{\text{bc}}\equiv\sqrt{E^{2}_{\text{beam}}-p^{2}}$.
Here, $E(p)$ is the total measurement energy (momentum) of the $\bar{p}\gamma\gamma$ combinations in the c.m.~system,
and $E_{\rm beam}$ is the beam energy.
Candidates are accepted with optimized $\Delta E$ requirements of $-16<\Delta E<7$~MeV to $-24<\Delta E<13$~MeV, depending on c.m.~energies, and with $M_{\text{bc}}>1.15$~GeV/c$^{2}$.

The \mycircled{1}\mycircled{4} and \mycircled{2}\mycircled{4} final states, classified as category B,
are reconstructed by requiring two good charged tracks with one identified as a $\pi^{-}$ and the other
identified as either a $\pi^{+}$ or $p$,
and the most energetic shower in these events is assigned as the
$\bar{n}$ candidate. To discriminate $\bar{n}$-initiated showers from those produced by photons, three variables
are retained for further selection based on c.m.~energy-dependent requirements: the total energy in the $\bar{n}$-assigned  EMC shower,
the second moment of the shower~\cite{hyperon_bes3}, and the
number of crystals with above-threshold signals within a 40$^{\circ}$ cone around the shower.
After that, kinematic fits that include the $\bar{n}$ direction are performed to identify signal events.
Since the $\bar{n}$ shower does not provide a good measure of its total energy, $E_{\bar{n}}$, this is left
as a free parameter in the kinematic fits.
If a $\pi^{+}$ is identified, the fit imposes the $n\bar{n}\pi^{+}\pi^{-}$ hypothesis with a missing $n$.
If a $p$ is identified, the fit imposes the $p\bar{n}\pi^{-}\pi^{0}$ hypothsis with a missing $\pi^{0}$.
In both fits, total energy-momentum conservation is constrained and $M_{\bar{n}\pi^{-}}$ is also constrained to the mass of the $\bar{\Sigma}^{-}$.
The $p\pi^{-}$ invariant mass is required
to be  $|M(p\pi^{-})-m(\Lambda)|>0.005$~GeV/$c^{2}$ to eliminate background from
$e^{+}e^{-}\to\Lambda\bar{\Lambda}\to p\pi^{-}\bar{n}\pi^{0}$.
Furthermore, the $\chi^{2}$ value from the kinematic fit is required to be less than 20.

The reconstruction of $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$ is similar to that for \mycircled{2}\mycircled{4}
in the $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$ analysis since they have the same final states.
The only difference is that
$M_{\bar{n}\pi^{+}}$ is constrained to the mass of the $\bar{\Sigma}^{+}$ in the kinematic fit.

Figure~\ref{stacksigmap}
shows the distributions of $M_{bc}$ for category A and the recoil mass of $\bar{n}\pi^{-}$, $M^{\rm rec}_{\bar{n}\pi^{-}}$,
for category B using selected $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$ candidates,
where significant signals in both categories are observed in data at $\sqrt{s}=2.3864$ and
2.3960~GeV.
Backgrounds are studied with MC
samples and only hadronic final states survive the selection criteria. In category~A,
the backgrounds are from $e^+e^-$ annihilation events with the same final states as the signal
process, with one or more additional $\pi^{0}$, and with an additional $\gamma$-ray.
In category~B,  the backgrounds are from annihilation events with the same final states as the
signal process, multi-$\pi$ processes such as $\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ and processes with one more
$\pi^{0}$ in the final states.
These background processes are mainly from contributions including intermediate states such as $\Delta$,
$\Lambda$ and $\Sigma$ baryons, but none of them produce peaks in the signal regions as shown by the histograms of  Fig.~\ref{stacksigmap}.
Figure~\ref{stacksigmam} shows distributions of $M_{n\pi^{-}}$ for
$e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$ candidate events at $\sqrt{s}=2.3960$ and
2.6444~GeV, respectively, where significant signals in data are observed.
In the background study, no peaking background is observed in the $n\pi^{-}$ mass spectrum.

\begin{figure}[htbp]
\begin{center}
\begin{overpic}[width=8.7cm,height=3.5cm,angle=0]{fit_minusv8}
\put(14,35){\small{a)}}
\put(63,35){\small{b)}}
\put(25,35){\tiny{$\sqrt{s}=2.3960$~GeV}}
\put(75,35){\tiny{$\sqrt{s}=2.6444$~GeV}}
\end{overpic}
\caption{ (color online) The $M_{n\pi^{-}}$ distributions for selected $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$
events at a) $\sqrt{s}=2.3960$~GeV and b) $\sqrt{s}=2.6444$~GeV.
Dots with error bars are the data; histograms are the background events in
MC after normalization. Solid curves are the fit results,
 dashed curves are the signal, and dot-dashed curves are the background.}
\vspace{-0.5cm}
\label{stacksigmam}
\end{center}
\end{figure}

\begin{table*}[htbp]
\caption{Summary of the calculated cross section and effective FFs at each c.m.~energy and the quantities used in the
calculation, $\epsilon=\varepsilon(1+\delta^{r})(1+\delta^{v})\delta^{\text{data/MC}}$, defined in the text.
The energy points with asterisks are combined data samples with c.m~energies weighted by
the luminosities of the subsamples. The 2.7500~GeV is a combined data set of 2.7000 and 2.8000~GeV, and 2.9884~GeV is a combined data set of
 2.9500, 2.9810, 3.0000 and 3.0200~GeV.
}
\vspace{-0.5cm}
\label{sumsigma}
  \footnotesize

  \begin{center}
  \begin{tabular}{ c|c|ccccc|ccccc}
  \hline
  \hline
  &   &  \multicolumn{5}{c|} {$e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$} &  \multicolumn{4}{c} {$e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$}  \\
\hline
  $\sqrt{s}$~(GeV) & $\mathcal{L}$~(pb$^{-1}$) &  $\epsilon_{A}$(\%)  &  $\epsilon_{B}$(\%) &  $\sigma^{\text{Born}}$~(pb) & $|G_{\text{eff}}|$($\times10^{-2})$ &  $|G_{E}/G_{M}|$  & $\epsilon$(\%)  &  $N$ & $\sigma^{\text{Born}}$~(pb) & $|G_{\text{eff}}|$($\times10^{-2})$  \\
\hline
\hline
  2.3864	   & 22.6    	& 5.8  & 12.6  &  $58.2\pm5.9^{+2.8}_{-2.6}$& $16.5\pm0.9\pm0.9$   &-       & \multicolumn{4}{c}{(below threshold)}\\
  2.3960		   & 66.9    	& 9.5 & 14.1 &  $68.6\pm3.4\pm2.3$     & $15.0\pm0.4\pm0.5$      & $1.83\pm0.26\pm0.24$        & 18.8 & $29.6\pm6.7$ & $2.3\pm0.5\pm0.3$ & $3.9\pm0.5\pm0.6$\\
  2.5000		  & 1.10		    & 18.4 & 21.6 &  $130\pm29\pm11$	& $14.0\pm1.6\pm0.6$       &-        & 20.2  &$4.8^{+2.9}_{-2.2}$ & $21.2^{+12.7}_{-9.5}\pm1.4$ & $5.9^{+1.8}_{-1.3}\pm0.2$\\
  2.6444	  & 33.7 	    & 24.4 & 20.5  &  $59.9\pm3.6\pm3.2$    &$8.6\pm0.3\pm0.2$      	& \multirow{2}{*}{$0.66\pm0.15\pm0.11$}  & 16.7  & $33.1\pm7.7$  & $5.8\pm1.4\pm0.4$ & $2.8\pm0.3\pm0.1$\\
  2.6464	& 34.0			& 24.2  & 20.7 &  $58.9\pm3.5\pm2.4$   &$8.5\pm0.3\pm0.2$       &          & 16.8  &$38.0\pm8.4$  & $6.6\pm1.5\pm0.5$ & $2.9\pm0.3\pm0.1$\\
  *2.7500	 & 2.04			& 25.0  & 19.7 &  $36.9\pm12.8\pm3.2$  & $6.7\pm1.2\pm0.3$       & -        & -   & - & - & -  \\
  2.9000   	& 105.	       & 26.5  & 20.6 &  $16.7\pm1.2\pm1.1$   & $4.5\pm0.2\pm0.2$       	&  $1.06\pm0.36\pm0.09$      & 14.2 & $18.0\pm7.1$ & $1.2\pm0.5\pm0.1$ & $1.2\pm0.2\pm0.1$\\
  *2.9884	& 65.2			& 25.5  & 21.4 &  $12.4\pm1.3\pm1.3$    & $3.9\pm0.2\pm0.2$      & -       &14.9 & $9.4^{+5.4}_{-4.6}$ & $1.0^{+0.6}_{-0.5}\pm0.1$ & $1.1\pm0.3\pm0.1$\\
  \hline
  \hline
  \end{tabular}
  \end{center}
  \end{table*}

The Born cross section for $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$ is determined from the relation:
\begin{equation}
\label{calsigma}
\sigma^{\text{B}} = \frac{N_{i}}{\mathcal{L}(1+\delta^{r})\frac{1}{|1-\Pi|^{2}}\delta^{\text{data/MC}}_{i}\mathcal{B}_{i}\varepsilon_{i}},~(i=A,B),
\end{equation}
where $N$ is the signal yield extracted from the fits;
$\mathcal{L}$ is the integrated luminosity;
$1+\delta^{r}$ is the ISR correction factor incorporating the input cross section from
this analysis iteratively;
$\frac{1}{|1-\Pi|^{2}}$ is the vacuum polarization factor~\cite{VP};
$\varepsilon$ is the detection efficiency determined from signal MC events.
The factor $\delta^{\text{data/MC}}$ is a correction factor for efficiency differences between data and MC simulation,
determined from studies of high statistics, low-background control samples of
$J/\psi\to\Sigma^{+}\bar{\Sigma}^{-}$ and $J/\psi\to\Lambda\bar{\Sigma}^{-}\pi^{+}$, respectively.
The decay branching fraction $\mathcal{B}$ accounts for the intermediate states in the $\bar{\Sigma}^{-}$ decay
(51.57\% for $\bar{\Sigma}^{-}\to \bar{p}\pi^{0}$ and  48.31\% for $\bar{\Sigma}^{-}\to \bar{n}\pi^{-}$).

To determine the signal yields, un-binned maximum likelihood fits are performed to the $M_{\text{bc}}$ and $M^{\rm rec}_{\bar{n}\pi^{-}}$
distributions for categories A and B, respectively.
The probability density function (PDF) for the signal is described with a MC-simulated shape
convolved with a Gaussian function to account for mass resolution differences between
data and MC simulation.
The background PDF for category A is an Argus function~\cite{argus};
for category B it is a second order polynomial.
In the fit, the two categories are constrained by the same Born cross section $\sigma^{\text{Born}}$,
and the expected signal yields are calculated from
$N_{i}=\sigma^{\text{Born}}\cdot \mathcal{L}\cdot \varepsilon_{i}\cdot (1+\delta)\cdot \delta^{\text{data/MC}}_{i}\cdot \mathcal{B}_{i}$.
The fit results at $\sqrt{s}=2.3864$ and 2.3960~GeV are shown in Fig.~\ref{stacksigmap}.
Similarly, the signal yield of $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$ is determined by fitting
the $n\pi^{-}$ mass spectrum, where the signal is described with the MC simulated shape convolved with a Gaussian
 function and the background is described with a 2nd-order polynomial.
Fit results at $\sqrt{s}=2.3960$ and 2.6444~GeV are shown in Fig.~\ref{stacksigmam}.

The quantities used in the cross section calculations are summarized in Table~\ref{sumsigma}.
It should be noted that, due to limited statistics,
data at c.m~energies 2.7000 and 2.8000~GeV are combined;
data at   2.9500, 2.9810, 3.0000
and 3.0200 GeV are combined.
Currently, individual measurements on $|G_{E}|$ and $|G_{M}|$ at each energy point
are not possible. Therefore, the effective FFs of $\Sigma^{\pm}$, defined as
$|G_{\text{eff}}|^2 \equiv (|G_E|^2+2\tau|G_M|^2)/(2\tau +1)$~\cite{effectiveff}, are reported here and shown in Table~\ref{sumsigma}.

Systematic uncertainties associated with the cross section measurements include event selection,
cross section line-shape, angular distribution, fitting method, energy scale, and luminosity.
In the nominal results, the differences  of data and MC efficiencies are corrected with control samples.
We vary the data/MC correction factors within their $\pm 1\sigma$ uncertainty
and the resulting  differences in the cross sections are taken as the uncertainty from the event selection.
The uncertainty associated with the cross section line-shape
is 1.0\%, which includes both the theoretical uncertainty and the parameter uncertainty in the line-shape fit.
The uncertainty from the angular distribution is evaluated by varying $|G_{E}/G_{M}|$ ratios within $\pm 1\sigma$ at the
three energy points with the highest statistics. For the energy points with
unknown $|G_{E}/G_{M}|$ values, two extreme cases
$G_{E}=0$ and $G_{M}=0$ are considered and the difference in the efficiencies divided by
a factor of $\sqrt{12}$ is taken as the uncertainty~\cite{systematic}.
Alternative fits are performed to study the uncertainty from the fit procedure. These include varying the fitting range, varying the signal
shape by fixing the resolution of the convolved Gaussian to be $\pm1\sigma$ different from its nominal value,
and changing the background PDF from a second order to a third order polynomial.
The effects of the c.m.~energy and energy resolution uncertainties are studied for energy points near threshold.
The difference of the cross sections in  $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$ is very small and the corresponding uncertainty on
the cross sections can be neglected.
The uncertainty on the effective FFs  are 4.9\% and 2.8\%  at $\sqrt{s}=2.3864$ and 2.396 GeV due to the change
of Coulomb correction factors.
For the $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$ process, the variation of c.m~energy and energy resolution introduce
uncertainties of 12.0\%
and 14.2\% in the cross section and effective FF, respectively, at $\sqrt{s}=2.396$~GeV.
The integrated luminosity is determined with large angle Bhabha events with an uncertainty of 1.0\%~\cite{dataset}.
All sources of systematic uncertainties are treated as uncorrelated and summed in quadrature;
they are in the range between  3.5\% and 13.0\% of the cross sections, depending on the c.m.~energy.

\begin{figure}[htbp]
\begin{overpic}[width=6.5cm,height=4.5cm,angle=0]{fit_xsecv10.eps}
\end{overpic}
\caption{ (color online) The cross section lineshapes for $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$ (circles) and $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$ (squares).
 The solid line is the pQCD fit for $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$  and the dashed line for $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$.
The vertical dashed and dotted lines denote production thresholds for $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$ and $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$.}
\label{comgeff}
\end{figure}

The measured cross section line-shapes of $e^{+}e^{-}\to\Sigma^{\pm}\bar{\Sigma}^{\mp}$
from $\sqrt{s}=2.3864$ to 3.0200~GeV
are shown in Fig.~\ref{comgeff}.
The near threshold cross sections for $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$
and  $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$ are measured to be $58.2\pm5.9^{+2.8}_{-2.6}$
and $2.3\pm0.5\pm0.3$~pb, respectively,
both are inconsistent with the value of 520~pb expected for point-like charged baryons.
Instead, a new feature is observed in which
the cross sections for $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$ are consistently
smaller than those for $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$.
A pQCD-motivated energy power function~\cite{rinaldo,QCD}, given by
\begin{equation}
\label{pqcd}
\sigma^{\text{B}}(s) = \frac{\beta C}{s}(1+\frac{2m_{B}^{2}}{s})\frac{c_{0}}{(s-c_{1})^{4} (\pi^{2}+\text{ln}^{2}(s/\Lambda_{\text{QCD}}^{2}))^{2}}
\end{equation}
is used to fit the line-shapes, where $c_{0}$ is the normalization, $c_{1}$ is the mean effect of a set of intermediate states that mediates the coupling between the virtual photon~\cite{VMD}
 and is regarded as common for the two processes, and $\Lambda_{\text{QCD}}$ is the QCD scale, fixed to 0.3~GeV.
The fit results are shown in Fig.~\ref{comgeff} with a fit quality of $\chi^{2}/{\rm ndof}=9.7/12$,
where $\rm ndof$ is number of degrees of freedom.
The cross section ratio between $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$ and $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$ is obtained from $c_{0}$
to be $9.7\pm1.3$, and $c_{1}$ is $2.0\pm0.2$~GeV$^{2}$.
Since the effective FF is proportional to the square root of the Born cross section,
the ratio of the effective $\Sigma^{+}$ and $\Sigma^{-}$ FFs is consistent with 3,
which is the ratio of the incoherent sum of the squared charges of the $\Sigma^{+}$ and $\Sigma^{-}$ valence quarks, $\sum_{\text q \in \text B}Q_{q}^{2}$.
The results are in disagreement with the prediction from octet baryon wave functions~\cite{ratio_theory2}, where
the typical $SU(3)$-symmetry breaking effects for hyperon FFs are about $10\sim30\%$.
The U\&A model predicts that the $\Sigma^{+}$ FF is  one order of magnitude smaller than that for the
$\Sigma^{-}$~\cite{UandA}, in sharp contrast to what we measure.
In ChPT, the $\Sigma^{+}$ and $\Sigma^{-}$ FFs are both predicted to be larger than that of the proton~\cite{ChPT}.
In the di-quark model, the $\Sigma^{+}$ FFs should be comparable to that of  $\Lambda$~\cite{diquark}.
Our measurements are not consistent with any of these theoretical predictions.
We notice that a recent prediction for the non-resonant cross section of $e^{+}e^{-}\to \Sigma^{\pm}\bar{\Sigma}^{\mp}$ at the $J/\psi$ mass~\cite{rinaldo2}, based on an effective Lagrangian density,
 is consistent with our result when extrapolated to $\sqrt{s}=3.097$~GeV using Eq.~(\ref{pqcd}).

The value of $|G_{E}/G_{M}|$ can be obtained by fitting the differential
angular distribution according to Eq.~(\ref{equ:eq1}). The large statistics
at $\sqrt{s}=2.3960$, 2.6444, 2.6464 and 2.9000~GeV for
 $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$ allow us to perform a study of the
polar angle of $\Sigma^{+}$ in the c.m.~frame.
The  angular distributions for categories~A and~B at $\sqrt{s}=2.3960$~GeV
are shown in Fig.~\ref{figangular2396}.
These angular distributions have been corrected for the detection efficiency and ISR,
which are obtained from signal MC simulation. Additional bin-by-bin corrections
 due to the data/MC detection differences,
for categories~A and~B, respectively, have also been applied.
Simultaneous fits to the two data sets
to the expression in Eq.~(\ref{equ:eq1}) sharing a common value for $|G_{E}/G_{M}|$ are performed.
The result of $|G_{E}/G_{M}|=1.83\pm 0.26$ is significantly higher than 1.
Using the normalized number of events, $|G_{M}|$ is determined to be
$(9.14\pm 1.42)\times10^{-2}$ and
 $(9.30\pm 1.53)\times10^{-2}$ for category A and B, respectively.
Similar angular distribution fits are performed for the combined $\sqrt{s}=2.6444$
and 2.6464~GeV data sets, denoted as 2.6454~GeV, and  $\sqrt{s}=2.90$~GeV and the results are listed in Table~\ref{sumsigma}.
The systematic uncertainties on $|G_{E}/G_{M}|$ considered here are the difference between data and MC efficiency, the bin size, and the
fit range.
For the $\Sigma^{-}$, on the other hand, the statistics only allow for the determination of $|G_{\rm eff}|$; they are not sufficient to extract $|G_{E}/G_{M}|$.

\begin{figure}[htbp]
\begin{center}
\begin{overpic}[width=8.7cm,height=3.5cm,angle=0]{fit_ang_2396.eps}
\put(15, 32){\small{a)}}
\put(62,32){\small{b)}}
\end{overpic}
\caption{(color online) Simultaneous fit of efficiency corrected angular distribution at $\sqrt{s}=2.396$~GeV
for a) category A b) category B for $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$ events.
Dots with error bars are data, solid curves are the fit results, the contributions from $G_{E}$
and $G_{M}$ are indicated by dashed and dotted curves.}
\vspace{-0.5cm}
\label{figangular2396}
\end{center}
\end{figure}

In summary, the data collected by BESIII at c.m.~energies between $2.3864$ and 3.0200~GeV,
are exploited to perform measurements of $e^{+}e^{-}\to\Sigma^{\pm}\bar{\Sigma}^{\mp}$.
This is the first time that cross sections of $e^{+}e^{-}\to\Sigma^{\pm}\bar{\Sigma}^{\mp}$
in the off-resonance region are presented.
The precision has been significantly improved by reconstructing all dominant decay modes of the $\Sigma$.
Nonzero cross sections near threshold are observed for
$e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$ and $e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$
to be $58.2\pm5.9^{+2.8}_{-2.6}$ and $2.3\pm0.5\pm0.3$~pb, respectively.
The values disagree with the point-like expectations near threshold, $848(m_{p}/m_{B})^{2}$~pb, as has been seen for the proton~\cite{time_babar,time_cmd}.
The cross section line-shapes for $e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-}$ and
$e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+}$ are well-described by pQCD-motivated functions.
The ratio of the $\sigma^{\text{Born}}(e^{+}e^{-}\to\Sigma^{+}\bar{\Sigma}^{-})$ to
$\sigma^{\text{Born}}(e^{+}e^{-}\to\Sigma^{-}\bar{\Sigma}^{+})$ is determined to be $9.7\pm1.3$,
which is inconsistent with predictions from various models~\cite{ratio_theory2, UandA, ChPT, diquark}. The EMFF ratio $|G_{E}/G_{M}|$ of the $\Sigma^{+}$ is determined from its production angle dependence
at three high-statistics energy points.
The $|G_{E}/G_{M}|$ of the $\Sigma^{+}$ shows similar features to those of the proton~\cite{time_babar,time_bes3}, $\Lambda$~\cite{hyperon_bes32},
 and $\Lambda_{c}$~\cite{lambdac_bes3},
that is larger than 1 within uncertainties near threshold
 and consistent with 1 at higher c.m.~energies. \\

The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center and the supercomputing center of USTC for their strong support. This work is supported in part by National Key Basic Research Program of China under Contract No. 2015CB856700; National Natural Science Foundation of China (NSFC) under Contracts Nos.~11625523, 11635010, 11735014, 11822506, 11835012, 11935015, 11935016, 11935018, 11961141012, 11335008, 11375170, 11475164, 11475169, 11625523, 11605196, 11605198, 11705192; the Chinese Academy of Sciences (CAS) Large-Scale Scientific Facility Program; Joint Large-Scale Scientific Facility Funds of the NSFC and CAS under Contracts Nos.~U1732263, U1832207, U1532102, U1732263, U1832103; CAS Key Research Program of Frontier Sciences under Contracts Nos.~QYZDJ-SSW-SLH003, QYZDJ-SSW-SLH040; 100 Talents Program of CAS; INPAC and Shanghai Key Laboratory for Particle Physics and Cosmology; ERC under Contract No.~758462; German Research Foundation DFG under Contracts Nos.~443159800, Collaborative Research Center CRC 1044, FOR 2359; Istituto Nazionale di Fisica Nucleare, Italy; Ministry of Development of Turkey under Contract No.~DPT2006K-120470; National Science and Technology fund; STFC (United Kingdom); The Knut and Alice Wallenberg Foundation (Sweden) under Contract No.~2016.0157; The Royal Society, UK under Contracts Nos.~DH140054, DH160214; The Swedish Research Council; Olle Engkvist Foundation under Contract No 200-0605; U. S. Department of Energy under Contracts Nos. DE-FG02-05ER41374, DE-SC-0012069.