# Multivariate Spatial-Temporal Variable Selection with Applications to
Seasonal Tropical Cyclone Modeling

Marcela Alfaro<EMAIL_ADDRESS>[ Montserrat
<EMAIL_ADDRESS>[ Joseph<EMAIL_ADDRESS>[
Lian Xie<EMAIL_ADDRESS>label=u1 [[ url]http://www.cfdl.meas.ncsu.edu
Department of Statistics, North Carolina State Universitym1 and Department of
Marine, Earth and Atmospheric Sciences, North Carolina State University m2

###### Abstract

Tropical cyclone and sea surface temperature data have been used in several
studies to forecast the total number of hurricanes in the Atlantic Basin. Sea
surface temperature (SST) and latent heat flux (LHF) are correlated with
tropical cyclone occurrences, but this correlation is known to vary with
location and strength of the storm. The objective of this article is to
identify features of SST and LHF that can explain the spatial-temporal
variation of tropical cyclone counts, categorized by their strength. We
develop a variable selection procedure for multivariate spatial-temporally
varying coefficients, under a Poisson hurdle model (PHM) framework, which
takes into account the zero inflated nature of the counts. The method differs
from current spatial-temporal variable selection techniques by offering a
dynamic variable selection procedure, that shares information between
responses, locations, time and levels in the PHM context. The model is used to
study the association between SST and LHF and the number of tropical cyclones
of different strengths in 400 locations in the Atlantic Basin over the period
of 1950-2013. Results show that it is possible to estimate the number of
tropical storms by season and region. Furthermore, the model delimits areas
with a significant correlation between SST and LHF features and the occurrence
and strength of TCs in the North Atlantic Basin.

arXiv:0000.0000

, , and

##  1 Introduction

A tropical cyclone is a rotating system characterized by a low-pressure
center, strong winds, and heavy rain that can cause severe damage both out at
sea and inland. Because of the numerous deaths and damages caused in the past
(Xie and Liu (2014)), it is crucial to develop methods to model and predict
the locations and strengths 111Throughout this paper we use the word strength
to refer to the classification of the storm according to its strength. We use
the word intensity to refer to the $\lambda$ parameter in the Poisson hurdle
model. of storms. According to Keith and Xie (2009), a hurricane is an
environmental heat engine driven by sensible and latent heating from the
ocean. Measures of sea surface temperature (SST) and latent heat flux (LHF)
values in certain areas and during some months of the year are correlated with
cyclone activity (Gray et al. (1994), Xie and Liu (2014) and Xie (2005)), and
thus should be important variables to describe the number of tropical cyclones
(TC) in space and time.

Numerous studies have described the relationship between SST and the number of
TCs per season, with independent models per subregions in the North Atlantic
Basin, delimited by $10^{\circ}N$ to $62^{\circ}N$ latitude and $10^{\circ}W$
to $110^{\circ}W$ longitude; see Figure 1. Lehmiller, Kimberlain and Elsner
(1997), Blake and Gray (2004), Xie (2005), Elsner and Jagger (2006), Keith and
Xie (2009), Werner and Holbrook (2011) and Xie et al. (2014) use climate
indices as covariates in models that predict the number of hurricanes in each
of the following subregions: Caribbean Ocean, Gulf of Mexico and the rest of
the North Atlantic Basin. While some of the climate indices are summaries of
SST from specific geographical locations and months of the year, none of the
models in the literature study all possible combinations of months and
locations.

Previous investigations using spatial models have determined that it is
possible to describe the spatial variation of hurricanes using a finer scale
than the regional analyses discussed above. Elsner, Liu and Kocher (2000) show
that the location of hurricanes on a grid can be modeled using spatial models
and climate indices. Jagger, Niu and Elsner (2002) develop a space time
autoregressive model and Hodges, Jagger and Elsner (2014) use a Bayesian
hierarchical spatial model, with total hurricanes by location as a response
and climate indices per location as covariates. Jagger, Niu and Elsner (2002)
investigate forecasting skill for the model when it predicts hurricane
activity for each year in each location, but pointed out that their estimators
failed to converge using a grid size smaller than 5 by 5 degrees.

One important feature of the TC data on a grid is that the number of storms in
a given year is likely to be zero in most of the locations in the Atlantic.
Since the probability of an occurrence is directly related to the expected
number of occurrences in a year, we propose a model that appropriately model
both parameters and its association. Neelon, Ghosh and Loebs (2013) propose
the Poisson Hurdle Model (PHM) for such cases, as they warn that failing to
account for this dependence in zero-inflated models may produce biased
parameter estimates.

In this paper, we model the spatial-temporal variation of the yearly number of
three strengths of TCs in the Atlantic Basin, and estimate their respective
multivariate spatial-temporal association with SST and LHF. Simultaneously
modeling the multivariate spatial-temporal process improves our ability to
find features in the space-time covariates that are significantly associated
with TC activity in specific locations. To this end, we propose a model to
address the following challenges: (1) a multivariate response that is
potentially zero inflated, (2) an association between response and covariates
that varies over space and time, (3) presence of correlation among
coefficients for different responses, locations and times. We use empirical
orthogonal functions (EOFs) to create fields of SST and LHF features, and use
their scores as covariates. We show that the multivariate spatial-temporal
variation of TCs can be explained by spatial-temporal features of SST and LHF,
using high spatial resolution and seasonal patterns in addition to the annual
signal. We use a Bayesian PHM with spatial-temporal dependence, and add a
prior that can facilitate the selection of areas with no association, taking
into account other correlations among coefficients.

Often the effect of covariates on the response is assumed to be constant
across locations and time. However, this assumption is inappropriate in our
case as TC association with SST and LHF features varies with strength,
location and time within a season (Hodges, Jagger and Elsner (2014) and Xie
and Liu (2014)). Thus, we use spatial-temporally varying coefficients (Gelfand
et al. (2003)) that can relax this assumption. Moreover, our proposed method
applies variable selection techniques in order to identify areas with
significant association in the spatial-temporal domain.

Many variable selection methods have been proposed for Bayesian models; for
example, stochastic search variable selection (SSVS) from George and McCulloch
(1993). Smith and Fahrmeir (2007) and Chipman (1996) propose variations from
SSVS, for single predictors and Scheel et al. (2013) propose to assume the
coefficients of multiple predictors of interest as a product of a binary
spatially dependent field and an independent Gaussian field, to perform
spatial variable selection. Reich et al. (2010), Lum (2012) and Cai et al.
(2013) describe SSVS specifically for spatially and spatial-temporally varying
coefficients in Bayesian models. Reich et al. (2010) use the spike and slab
prior from Ishwaran and Rao (2005) and describe a computationally efficient
way to estimate the model. Cai et al. (2013) propose a model that performs
variable selection for spatial-temporally varying coefficients using a
Dirichlet process prior and nonparametric techniques. Depending on the mixed
distributions, the interpretation of these coefficients becomes very
complicated and almost impossible. Boehm Vock et al. (2015) propose to use
copulas to retain interpretability of the coefficient. We propose a new
version of a copula prior that can create a continuous smooth prior surface
and perform multivariate spatial-temporal variable selection.

The novelty in our method is a dynamic variable selection procedure, which
shares information between responses, locations, time and levels in the PHM
context. This skill is essential in cases where we have a limited number of
variables available, with a strong dependence structure, as in TC forecast
models, but also in many other areas. For example, in epidemiology studies,
where one is interested in detecting location and time where covariates have
significant effects on a zero inflated response; or in forestry or ecology
where models for sparse associations of covariates and rare events in space
and time are also common.

This paper is organized as follows. Section 2 describes the data used in the
application. Section 3 introduces the statistical model and its estimation
details. Section 4 describes a simulation study that investigates the
performance of the variable selection technique, and Section 5 presents the
results from real data application of the model. Finally, Section 6 presents
our discussion and conclusions.

##  2 Data

Our main source of data is HURDAT (HURicane DATa), a publicly available
database that contains detailed information for all recorded TCs in the
Atlantic Basin since 1851. We construct a grid cell count using the 6 hour
track data from 1950 to 2013 presented in Figure 1. The polygon grids are
created using rectangles of 2.5 by 2.5 degrees over the Atlantic Basin, and
using the coordinates from other variables used in the model as centroids.

Figure 1: HURDAT 6 hour track data from 1950 to 2013

In meteorology, the 5 category Saffir-Simpson scale is used to describe
hurricanes. In practice, forecasters report tropical cyclones in three
strength groups according to said scale, which are:

  * •

Low: Tropical storms. TCs with sustained winds of less than 33 $m/s$.

  * •

Mid: Hurricanes of category 1 and 2. TCs with sustained winds between 33 and
49 $m/s$.

  * •

Strong: Major hurricanes of category 3, 4 and 5. TCs with sustained winds of
more than 50 $m/s$.

We use this classification to count the number of unique storms per strength
and year in each grid box. These counts are the response variable in our
model.

SST and LHF are used as covariates in our model. For SST, we use data from the
High-resolution Blended Analysis (Kalnay et al. (1996)), which takes
measurements from NOAA’s Advanced Very High Resolution Radiometer (AVHRR). SST
is reported at a resolution of 2.5 by 2.5 degrees. LHF is defined as the flux
of heat in watts per square meter ($W/m^{2}$) from the Earth’s surface to the
atmosphere, which is associated with evaporation of water at the surface and
subsequent condensation of water vapor in the troposphere. We use data from
the National Centers for Environmental Prediction (NCEP) Reanalysis (Kalnay et
al. (1996)), with a 2.5 by 2.5 degree resolution. Both SST and LHF have daily
observations available for the Atlantic Basin. We average the daily data over
three-month periods (trimesters) that coincide with the definition of season
from World Ocean Database (WOD) (Levitus (1983)), i.e. winter starting in
January, spring in April, summer in July and fall in October. The relationship
between trimesters $w$ and the tropical cyclone season is presented in Figure
2. We define $Z_{\ell}(s_{i},t,w)$ as the average of covariate $\ell$, in
trimester $w$, location $s_{i}$ and year $t$, and then anomalies
$X_{\ell}(s_{i},t,w)$ as departures from the temporal mean per trimester and
location ($\overline{Z}_{\ell}(s_{i},\cdot,w)$) in each variable. Anomalies
for some selected years are presented in Figure 3.

Figure 2: Definition of trimesters compared to the Tropical Cyclone Season for
TC data.

Our data consist of observations for $T=64$ years, $M=4$ trimesters per year,
$N=400$ locations, $K=3$ response variables (number of tropical cyclones of
low, mid and strong strength), and $L=2$ covariates: SST and LHF. Figure 3
presents maps of selected years as examples of the spatial locations in the
Atlantic.

Figure 3: North Atlantic Basin: Anomalies in the North Atlantic Basin during
the first trimester for SST (left column) and SST (middle column). Counts per
grid box of Tropical Storms of Strength 1 (right column). Years 1955, 1980,
and 2013 randomly selected to illustrate data.

##  3 Statistical Methods

In this section we describe our model for yearly counts of TCs per location
and strength. The model simultaneously associates the different count
distributions to covariates’ features derived from anomalies.

###  3.1 Pre-processing

Since the association between SST and LHF and tropical cyclones is not
necessarily limited to an association between measures in the same location,
we seek multivariate spatial-temporal features to be used as covariates in an
interpretative way, while still providing a practical fit to the data. We use
Empirical Orthogonal Functions (EOFs) to describe said features. According to
Lorenz (1956), EOFs are spatial components that display space-time modes of
variability of a quantity over a region. Thus, each covariate anomaly, for
each trimester $w$, ${\bf{X}}_{\ell}(w)=X_{\ell}(\cdot,\cdot,w)$ is
approximated by the principal component (PC) scores, calculated using singular
value decomposition (SVD) of ${\bf{X}}_{\ell}(w)$, for each $w$ and $\ell$. In
this way, we can represent the anomalies as:

$X_{\ell}(s_{i},t,w)=\sum\limits_{r=1}^{R}\xi_{\ell,r}(t,w)\phi_{r}(w,s_{i})+\epsilon(s_{i},t,w),$
(1)

where R is the number of scores used to explain more than $70\%$ of its
variation. EOFs $\phi_{r}(w,s_{i})$ are defined for each trimester $w$,
$\xi_{\ell,r}(t,w)$ are the correspondent PC scores for each covariate $\ell$,
and trimester $w$, and $\epsilon(s_{i},t,w)$ is the variation orthogonal to
the first R EOFs.

###  3.2 Model

Let $y_{k}(s_{i},t)$ be the number of unique TCs in year $t$ of strength $k$
in location $s_{i}$. We model $y_{k}(s_{i},t)$ using a hierarchical Poisson
hurdle model (PHM), where the probability mass function is defined as:

$P(y_{k}(s_{i},t)=m|p_{k}(s_{i},t),\lambda_{k}(s_{i},t))=\left\\{\begin{tabular}[]{lr}$1-p_{k}(s_{i},t)$&if
$m=0$\\\
$p_{k}(s_{i},t)\frac{\lambda_{k}(s_{i},t)^{m}}{(\exp^{\lambda_{k}(s_{i},t)}-1)m!}$&if
$m>0$\\\ \end{tabular}\right.$ (2)

where $p_{k}(s_{i},t)$ is the probability of non-zero counts and
$\lambda_{k}(s_{i},t)$ is the intensity parameter. Both $p_{k}(s_{i},t)$ and
$\lambda_{k}(s_{i},t)$ are described as a logit and loglinear function of
multivariate spatial-temporally varying coefficients of the covariates’
scores, respectively. Thus we have:

$\begin{split}logit(p_{k}(s_{i},t))&=\sum\limits_{\ell=1}^{L}\sum\limits_{w=1}^{M}\sum\limits_{r=1}^{R}\beta_{1,k,\ell,r}(s_{i},w)\xi_{\ell,r}(t,w)\\\
log(\lambda_{k}(s_{i},t))&=\sum\limits_{\ell=1}^{L}\sum\limits_{w=1}^{M}\sum\limits_{r=1}^{R}\beta_{2,k,\ell,r}(s_{i},w)\xi_{\ell,r}(t,w).\end{split}$
(3)

Since our model coefficients vary in space by trimesters and by strength they
are able to explain the association between the space-time modes of
variability in the covariates (LHF and SST) and both the probability of non-
zero counts and the intensity, for every strength. This association is
expected to be close to zero in many spatial locations and during some
trimesters, depending on the TC strength. To this end, we propose a method
that can identify where, when, and for which strength this association is
different from zero. Furthermore, we present a prior for $\beta$ that selects
such cases and uses the correlation between responses, locations, time and
levels in the PHM context to improve this estimation.

###  3.3 Multivariate Spatial-Temporal Prior

Let $A$ be a combination of covariate $\ell$, PC score $r$, response $k$ and
level of PHM $j$. Each coefficient $\beta_{A}(s_{i},w)$ is modeled using a
Gaussian copula and a latent variable $\theta_{A}(s_{i},w)$. We define
$\theta_{A}(s_{i},w)\sim\text{MVN}(0,\Sigma)$ and use $\Sigma$ with a
separable covariance structure with the following form:

$\Sigma_{\theta}=\Gamma_{s}\otimes\Gamma_{w}\otimes\Gamma_{k}\otimes\Gamma_{j}.$
(4)

The term separable means that we can factor the covariance structure into
purely spatial ($s$), temporal ($w$), multivariate ($k$), and by PHM level
($j$) covariance structures, which allows for computationally efficient
estimation and inference. Each covariance $\Gamma$ can be written as a
correlation function $C(h)$ where $h$ is the Euclidean distance when the
distance is defined in space ($s$) or in time ($w$), or the numerical absolute
difference when the function is in terms of the ordered response strengths
($k=1,2,3$) or PHM levels ($j=1,2$).

Then, we define $\beta_{A}(s_{i},w)=F^{-1}_{A}\\{\Phi[\theta_{A}(s_{i},w)]\\}$
where $\Phi$ is the cdf of a standard normal distribution and $F_{A}$ is the
marginal cdf of each $\beta_{A}(s_{i},w)$, distributed as a multivariate spike
and slab mixture of normal densities:

$\beta_{A}\sim\pi_{A}N(\alpha_{A},\sigma^{2}_{A})+(1-\pi_{A})N(0,\sigma^{2}_{A}/C),$
(5)

where $\pi_{A}$ represents the expected proportion of coefficients that are
centered at $\alpha_{A}$. In this way, coefficients are smoothed around
$\alpha_{A}$ or 0. $\sigma^{2}_{A}$ controls variability around the means of
the mixed distribution, and C represents the ratio of variance for
coefficients centered in zero to centered in $\alpha_{A}$.

A copula can flexibly model the marginal distributions while introducing
dependence via multivariate normal distributions, creating a continuous smooth
prior surface. The intuition in our case is that if a coefficient is close to
zero in a specific location and time, nearby locations (in space and time) are
unlikely to have large coefficients. Also, coefficients for adjacent
categories from an ordered response and the two levels of a PHM are expected
to be correlated, since the probability of having zero counts in a specific
location and time is likely to be associated with its intensity.

For $\pi_{A}$, we use a flat prior between 0 and 1. Other prior choices can
give more weight to options like $\pi_{A}=0$: all locations and trimesters
have a coefficient centered in zero for $A$, or the opposite $\pi_{A}=1$. For
$\alpha_{A}$, the prior is a normal distribution with standard deviation 1,
which allows the coefficients to be centered in values ranged from $[-4,4]$.
Previous work on mixed models similar to Equation 5 has defined that model fit
can be sensitive to large values of $\sigma_{A}$ because it would be difficult
to differentiate between coefficients centered at 0 from other values. We use
a Gamma(0.1,0.1) prior to encourage small values for $\sigma_{A}$.

We model the temporal and spatial covariance components using AR1 and
exponential structures, respectively, while using a vague Wishart prior to
model the multivariate covariance and the covariance between PHM levels.
Hyperpriors for the spatial range parameter $r$, and the temporal correlation
and variance $\rho_{t}$ and $\sigma_{t}$ are recommended to be informative, in
the sense of restricting the range of possible parameter values to values
which make sense in terms of the application.

To estimate the parameters in the model, we use a Markov chain Monte Carlo
(MCMC) algorithm. Since we are using a mixed distribution as a prior for
$\beta_{A}(s_{i},w)$, we do not have conjugacy, and we need to use rejection
sampling. The scheme requires factoring $24N\times 24N$ matrices that we
simplify by using the Woodbury matrix identity on the Kronecker structure of
$\Sigma_{\theta}$. In the simulation study in Section 4 we generate 15000
iterations, discard the first 5000 as burn-in, and thin the chains by keeping
every fifth iteration.

##  4 Simulation Study

We study the performance of our model compared to two other competing models
for variable selection. Our main focus is to describe the dynamic variable
selection performance; thus we compare performance for varying degrees of
spatial-temporal correlation in the coefficients $\beta_{A}(s_{i},w)$. For
that, we use the following statistics: median absolute deviation (MAD),
ability to correctly identify null coefficients in space and time and MSE for
each coefficient estimated.

We generate $5$ $\beta_{A}(s_{i},w)$ under three different settings. Setting 1
is composed of smooth surfaces created using linear combinations of latitude,
longitude, and $\alpha_{A}=\\{2,1.5,1,0.5,0\\}$, respectively. Figure 1 in
Acknowledgements shows an example realization. In this case
$\beta_{A}(s_{i},w)$ does not follow the same distribution as in any of our
models. Settings 2 and 3 are generated using the copula transformation
described in Section 3.3 and assuming 5 true, with
$\alpha_{A}=\\{2,1.5,1,0.5,0\\}$, $\pi=\\{1,0.8,0.5,0.2,0\\}$,
$\sigma^{2}_{A}=1$ and $C=\infty$. They differ in the degree of spatial-
temporal correlation. For each $A$ and each setting, we generate $M=4$ time
points and $N=100$ locations.

Samples of $\xi_{A}(w,t)\sim\text{N}(0,2^{-1/2})$ for $w=1,2,..,4$ and
$t=1,2,..,30$, with $A=1,2,..,5$ are generated to simulate scores. We apply
Gram-Schmidt orthonormalization to make them independent among $A$. Then, we
use $\beta_{A}(s_{i},w)$ from settings 1-3, to have
$\boldsymbol{\beta}^{T}\boldsymbol{\xi}$, and use it to generate random
samples of $Y(s_{i},t)$ from a Poisson distribution:
$\text{Pois}(\exp{(\boldsymbol{\beta}^{T}\boldsymbol{\xi})})$.

Following the procedures outlined above, we generate B = 50 data sets using
the following settings for $\beta_{A}(s_{i},w)$:

  * •

Setting 1: Weak Temporal Correlation and Strong Spatial Correlation. Smooth
surfaces are generated using linear combinations of latitude and longitude.

  * •

Setting 2: Weak Spatial-Temporal Correlation. The true coefficients for each
location and trimester are exponentially spatially correlated with range = 2
km and temporally correlated with $\rho=0.1$ for all A.

  * •

Setting 3: Strong Spatial-Temporal Correlation. The true coefficients for each
location and trimester are exponentially spatially correlated with range = 100
km and temporally correlated with $\rho=0.9$ for all A.

We test the performance of three models under three different data settings.
Since we concentrate on showing the advantages of sharing information between
locations and times when performing variable selection, we simplify the
response model to be $y(s_{i},t)\sim\text{Pois}(\lambda(s_{i},t)),$ and use
the second Equation in 3 to model the intensity parameter $\lambda(s_{i},t)$
as a function of the scores $\xi_{\ell,r}(w,t)$. Here, $A$ represents scores
for each covariate (component) that is used in our model. The three models we
consider are our model (Model 3) and two other competing models for variable
selection. They differ in the latent variable distribution used in the
transformation for $\theta_{A}(s_{i},w)\sim\text{MVN}(0,\Sigma_{\theta})$:

  * •

Model 1: $\beta_{A}(s_{i},w)$ are independent across sites and times and
follow the distribution in 5. It does not share information, since
$\Sigma_{\theta}=\sigma^{2}_{t}\times I$ is just a diagonal matrix.

  * •

Model 2: Spatial dependence in $\beta_{A}(s_{i},w)$ is induced with a Gaussian
Copula. It shares information across sites, using
$\Sigma_{\theta}=S\otimes\sigma_{t}^{2}\times I$ where $S$ represents an
exponential spatial correlation matrix with range r.

  * •

Model 3: Spatial-temporal dependence in $\beta_{A}(s_{i},w)$ is induced using
a Gaussian Copula. It shares information across sites and times, using
$\Sigma_{\theta}=S\otimes T$, where S is defined as in Model 2 and T has an
autoregressive 1 structure (AR1) with parameters $\rho_{t}$ and
$\sigma^{2}_{t}$ for time.

In all cases, we use the same priors : $\alpha_{A}\sim\text{N}(0,2)$,
$\pi_{A}\sim\text{Beta}(1,1)$, and
$\sigma^{2}_{A}\sim\text{InvGamma}(0.1,0.1)$. We include Model 1 as a baseline
to show the overall improvements in the estimation of $\beta_{A}(s_{i},w)$
when considering the spatial and the spatial-temporal structure that models 2
and 3 provide.

We compute posterior medians and credible intervals for each
$\beta_{A}(s_{i},w)$ and posterior means for $\pi_{A}$ and $\alpha_{A}$.
Median absolute deviation (MAD) is averaged over space, trimesters and
iterations for each $\beta_{A}$. Percentage of correctly estimated zeroes are
calculated only for $\beta_{A}$ and averaged over space, trimesters and
iterations as well. MSE is also calculated and averaged over iterations for
$\alpha_{A}$ and $\pi_{A}$.

MAD is computed in the following way, for each setting and model:

$\text{MAD}_{A}=\text{Median}_{\mathit{BNM}}|\hat{\beta}_{A}(s_{i},w)-\beta_{A}(s_{i},w)|$

where $B=50$ data sets, $N=100$ locations and $M=4$ time points.

Table 1: Mean Acceptance Rate for all parameters. | Weak T Correlation | Weak ST Correlation | Strong ST Correlation
---|---|---|---
| M1 | M2 | M3 | M1 | M2 | M3 | M1 | M2 | M3
$\beta$ | 0.75 | 0.74 | 0.44 | 0.81 | 0.79 | 0.76 | 0.74 | 0.79 | 0.43
$\alpha$ | 0.19 | 0.14 | 0.35 | 0.71 | 0.71 | 0.85 | 0.71 | 0.71 | 0.85
$\pi$ | 0.31 | 0.35 | 0.38 | 0.36 | 0.35 | 0.32 | 0.35 | 0.35 | 0.43

Table 2 presents MAD for each $\beta_{A}(s_{i},w)$. Values were compared using
credible intervals. We can see how the median absolute deviations for Model 3
are statistically smaller in most of the settings generated with strong
correlation, as expected. In some cases, this improvement is reflected in a
$40\%$ decrease of MAD. Nevertheless, this difference is not reflected in
cases where we have weak correlation and $\alpha$ is generated as far from
zero ($\beta_{1}$ and $\beta_{2}$). When comparing results for data generated
with weak temporal correlation, MAD for our model is smaller than models 2 and
3 in all cases. MAD is statistically smaller using Model 3 in cases where the
real $\alpha$ is equal or close to zero ($\beta_{3}$, $\beta_{4}$, and
$\beta_{5}$).

Table 2: Median MAD per $\beta_{A}$ over locations and trimesters. All
standard errors are less than 0.006. Proportion of correctly estimated
$\beta_{A}=0$ All standard errors are less than 0.005. NA in cases where we
did not simulate zeroes.

Weak T Correlation Weak ST Correlation Strong ST Correlation M1 M2 M3 M1 M2 M3
M1 M2 M3 MAD $\beta_{1}$ 0.3073 0.4018 0.3699 0.6235 0.5996 0.7316 0.7594
0.6712 0.6882 $\beta_{2}$ 0.3229 0.4006 0.3686 0.6368 0.5843 0.7327 0.6731
0.6787 0.6112 $\beta_{3}$ 0.3263 0.3440 0.2824 0.6940 0.7460 0.4722 0.9846
1.0040 0.5306 $\beta_{4}$ 0.3349 0.3380 0.1999 0.5110 0.5174 0.3615 0.9221
0.8880 0.3044 $\beta_{5}$ 0.3197 0.3268 0.2042 0.6480 0.6458 0.0154 0.7164
0.7049 0.2473 Proportion $\beta_{1}$ 0.8121 0.7437 0.8862 NA NA NA NA NA NA of
Correctly $\beta_{2}$ 0.8457 0.7594 0.8670 0.5258 0.4419 0.7860 0.7792 0.8430
0.8404 Estimated $\beta_{3}$ 0.7434 0.8369 0.9409 0.8925 0.8754 0.9437 0.6028
0.7602 0.9692 $\beta_{A}=0$ $\beta_{4}$ 0.7184 0.7346 0.8676 0.9586 0.9746
0.9894 0.8568 0.9476 0.9827 $\beta_{5}$ 0.8941 0.8541 0.9285 0.5510 0.5509
0.9980 0.5563 0.5517 0.8882

As a way to evaluate the performance of our dynamic variable selection method,
we calculate the proportion of correctly estimated $\beta_{A}=0$. Our model
detects a statistically higher proportion of real values equal to zero in all
cases when the data is generated with weak and strong correlation. On the
other hand, our model performs similarly to Model 2 in cases where we have
weak temporal correlation or when $\alpha$ is far away from zero. Simulation
results show an average increase of $13\%$, $38\%$, and $30\%$ in the
percentage of correctly detected zeroes, for each of the settings. The best
case scenario of performance for our dynamic variable selection procedure is
when data presents strong spatial-temporal correlation, and it is close to
zero.

Intuitively, MAD result reflects the fact that our model (M3) benefits from
using temporal correlation to get a better fit when there is strong spatial-
temporal correlation present, compared to the fit to weak spatial-temporal or
temporal dependent data. The same benefit is reflected in the percentage of
correctly detected zeroes, where sharing information in space and time gives
the model more tools to correctly detect zeroes. Therefore, when there is
strong spatial-temporal correlation present in the data, the simulation shows
the need of a model that accounts for both space and time dependence.

Overall, our model shows an average increase of $13\%$, $38\%$, and $30\%$ in
the percentage of correctly detected zeroes for weak temporal, weakly and
strongly correlated data, respectively. Also, on average our model has a MAD
per $\beta_{A}$ $38\%$ lower than the competing models for strongly correlated
data, $24\%$ lower for weakly correlated data, and $16\%$ lower with weak
temporal correlation.

##  5 Analysis of Tropical Cyclones in the Atlantic Basin

###  5.1 Model Comparison

We fit the full model (ST) from Section 3, and two other models (Spatial and
Independent) to the tropical cyclone data described in Section 2. The
difference between the full model (ST), spatial (SP) and independent (IN)
model is the specification of $\Sigma_{\theta}$ in 4. In the SP model we use
$\Sigma_{\theta}=\Gamma_{s}\otimes I_{w}\otimes I_{k}\otimes I_{j}$, and in
the IN model we assume independence over all levels
($\Sigma_{\theta}=I_{s}\otimes I_{w}\otimes I_{k}\otimes I_{j}$). We perform
the comparison in order to collect evidence of the need of the full model (or
the lack thereof).

For each trimester $w$, we calculate anomalies $X_{\ell}(s_{i},t,w)$ for both
SST and LHF and then perform a SVD as described in Section 3. We obtain R = 4
EOFs with their respective scores $\xi_{\ell,r}(t,w)$ for each covariate
$\ell$ and trimester $w$, to have 16 time series in total. The maps of the
four resulting EOFs for each covariate and trimester are presented in the
Acknowledgements. The obtained EOFs for each of the seasons and covariates
together describe at least $70\%$ of the variability. The estimated scores
range from $-6$ to $6$ in all cases and are centered close to zero. In the
results subsection we describe these spatial patterns as features and relate
them to the response with the statistically significant factors.

We compare models using DIC and MSE, to evaluate goodness of fit of the model,
but also their ability to estimate correctly the number of storms per location
and per year. The percentage of locations with no TCs in a given year in the
Atlantic Basin ranges from $78\%$ to $98\%$ for low strength, from $83\%$ to
$99\%$ for mid strength and from $89\%$ to $100\%$ for strong strength and the
range of counts per location is from 0 to 4 storms. Table 3 presents all the
goodness of fit measures for each model.

Table 3: Goodness of fit criterion DIC, and MSE overall, per location and per year for the full model (ST), spatial (SP) and independent (IN). DIC bar values (standard error in parenthesis) for models using different grid box sizes, compared to a model using climate indices. |  | MSE
---|---|---
| $\overline{DIC}$ | TS | HU | MH
Independent (IN) | 183.1 (1.2) | 3.16 (0.3) | 2.51 (0.2) | 0.08 (0.001)
Spatial (SP) | 142.9 (1.1) | 0.20 (0.01) | 0.18 (0.02) | 0.04 (0.001)
Spatial Temporal (ST) | 136.2 (4.2) | 0.13 (0.01) | 0.12 (0.001) | 0.03 (0.001)
200 grid boxes | 136.9 (4.5) |  |  |
400 grid boxes | 136.2 (4.2) |  |  |
800 grid boxes | 136.5 (4.3) |  |  |
Climate indices | 173.6 (5.8) |  |  |

The full model (ST) is selected as the best model by DIC and MSE values for
every response. We compare the observed with the expected values of the
posterior distribution and obtain that the full model performance in terms of
MSE per response is superior to the IN and SP models performance for all
responses. Figures 7 and 6 in the Appendix show the MSE per location and year,
respectively for the full model. Our model has a fairly superior performance
per year compared to climatology (MSE per response is $0.14$ for TC, $0.16$
for HU and $0.03$ for MH), with some exceptions like 1966, 1997 and 2005 for
HU. In addition, the full model presents MSE lower than 1 unit in most of the
locations, with a unique exception in front of the African Coast, very close
to the mid TC main development region. However, prediction is not the only
focus of our model. We also use the model to identify potentially
scientifically relevant associations between covariates and responses.

###  5.2 Sensitivity Study

We consider several choices of hyperpriors when fitting the full model. Also,
we estimate the scores using three different rectangle sizes as grid boxes in
space for SST and LHF measures. We evaluate the performance of the different
modeling options using the Deviance Information Criterion (DIC).

We specify different hyperpriors for the spatial range parameter $r$, and the
temporal correlation and variance $\rho_{t}$ and $\sigma_{t}$. Results were
robust in terms of DIC when changing priors for $\rho_{t}$ and $\sigma_{t}$
from informative to non informative. For this analysis we use a uniform prior
for $\rho_{t}$ and a Gamma(0.1,0.1) for $\sigma_{t}$. We observe that DIC is
sensitive to a change in $r$; goodness of fit improves significantly when
using a uniform distribution from $0$ to $300$ km, instead of $30$ km or
$3000$ km as upper limits. This improvement could be explained by the typical
ratio of a tropical cyclone eye, which approximately ranges from $3$ to $300$
km (NHC (2016)).

In addition to these hyperprior choices, we conduct a specific sensitivity
analysis to compare the goodness of fit when using different rectangle grid
box sizes to construct the SST and LHF scores, and at the same time compare it
with a model that uses climate indices instead. The use of climate indices as
base comparison is particularly interesting in this case, since those are the
covariates that are typically used in TC forecast models (Xie et al. (2014)).
DIC was nearly identical among the cases of using 200, 400 and 800 grid boxes
to construct the SST and LHF scores, but it was higher in the model that uses
climate indices instead of scores (Table 3). We set the rectangle size to be
the same as the response (2.5 by 2.5 degrees) and use 400 grid boxes to
calculate the SST and LHF scores. The first two scores per trimester are
presented in Acknowledgements.

In order to improve MCMC convergence, we tune $C$ to give acceptance rates of
around $30\%$ or higher for most parameters. Convergence is monitored using
trace plots of deviance and several representative parameters. For all the
sensitivity analyses and the final analysis, we generate 20,000 samples and
discard the first 5,000 as burn-in and thin the chains by keeping every fifth
iteration.

###  5.3 Results

Table 4 presents the posterior estimates for each response for the overall
mean $\alpha_{A}$, the proportion $\pi_{A}$ of coefficients centered in
$\alpha_{A}\neq 0$ and the variability $\sigma_{A}$ of coefficients centered
in $\alpha_{A}$. The proportion of locations where estimated coefficients are
significantly different from zero for SST and LHF features are presented in
Table 5; significant difference is determined by whether or not credible
intervals include zero.

Table 4: Posterior estimates for each response for the overall mean
$\alpha_{A}$, the proportion $\pi_{A}$ of coefficients centered in
$\alpha_{A}\neq 0$ and the variability $\sigma_{A}$ of coefficients centered
in $\alpha_{A}$, for responses $Y_{1}=$ Tropical Cyclones, $Y_{2}=$ Hurricanes
and $Y_{3}=$ Major Hurricanes. Values with (*) have credible intervals that do
not include zero.

Tropical Cyclones Hurricanes Major Hurricanes SST, $p$ LHF, $p$ SST, $\lambda$
LHF, $\lambda$ SST, $p$ LHF, $p$ SST, $\lambda$ LHF, $\lambda$ SST, $p$ LHF,
$p$ SST, $\lambda$ LHF, $\lambda$ $\alpha_{A},r=1$ -0.44 (*) 0.56 0.13 0.18
(*) 0.01 -0.02 0.11 0.29 -0.02 0.50 (*) 0.09 0.46 (*) $\alpha_{A},r=2$ -0.49
(*) -0.48 (*) -0.71 (*) -0.21 -0.40 (*) -0.36 0.09 -0.31 0.17 -0.33 (*) 0.04
-0.39 (*) $\alpha_{A},r=3$ -0.76 (*) 2.99 (*) -0.59 (*) 0.24 (*) -0.31 0.08
-0.16 -0.50 -0.38 -0.48 (*) -0.18 -0.76 $\alpha_{A},r=4$ -0.64 (*) -0.34 (*)
-0.03 -0.41 (*) 0.26 0.01 0.13 0.15 -0.07 -0.00 -0.07 0.09 $\pi_{A},r=1$ 0.19
(*) 0.16 (*) 0.01 0.08 (*) 0.03 0.52 (*) 0.12 (*) 0.15 (*) 0.04 0.18 (*) 0.09
(*) 0.02 $\pi_{A},r=2$ 0.22 (*) 0.28 (*) 0.33 (*) 0.80 (*) 0.22 (*) 0.06 0.19
(*) 0.17 (*) 0.05 0.28 (*) 0.22 (*) 0.05 $\pi_{A},r=3$ 0.21 (*) 0.14 (*) 0.33
(*) 0.67 (*) 0.04 0.04 0.13 (*) 0.19 (*) 0.25 (*) 0.23 (*) 0.05 (*) 0.02
$\pi_{A},r=4$ 0.24 (*) 0.21 (*) 0.01 0.18 (*) 0.17 0.55 (*) 0.12 (*) 0.17 (*)
0.07 0.05 0.21 (*) 0.16 (*) $\sigma_{A},r=1$ 0.10 (*) 0.09 (*) 0.05 (*) 0.05
(*) 0.11 (*) 0.51 (*) 0.48 (*) 0.64 (*) 0.17 (*) 0.21 (*) 0.49 (*) 0.09 (*)
$\sigma_{A},r=2$ 0.13 (*) 0.20 (*) 0.08 (*) 0.83 (*) 0.17 (*) 0.11 (*) 0.61
(*) 0.49 (*) 0.16 (*) 0.30 (*) 0.70 (*) 0.38 (*) $\sigma_{A},r=3$ 0.16 (*)
0.20 (*) 0.08 (*) 1.34 (*) 0.11 (*) 0.11 (*) 0.51 (*) 0.87 (*) 0.25 (*) 0.28
(*) 0.43 (*) 0.21 (*) $\sigma_{A},r=4$ 0.20 (*) 0.08 (*) 0.05 (*) 0.05 (*)
0.12 (*) 0.59 (*) 0.50 (*) 0.52 (*) 0.17 (*) 0.17 (*) 0.55 (*) 0.57 (*)

We define a significant factor as an EOF for which more than $10\%$ of the
corresponding spatially varying coefficients for a given trimester $w$ are
statistically different from zero (posterior credible interval do not include
zero). We describe here those factors that have a sum of magnitudes greater
than 20 for intensity and greater than 10 for probability of occurrence. SST
features are the only significant factors of coefficients during winter
($w=1$). These coefficients are the most important in terms of determining
predictors to use in the seasonal forecast, because they can be measured
before the TC season starts (Figure 2). For example, the model indicates that
warmer anomalies in SST over the mid TC main development region (in front of
the African Coast) and in the Arctic during the winter are positively
correlated with more major hurricanes in the Caribbean Sea and with more
hurricanes around Florida, Cuba and the Mid Atlantic. Furthermore, the model
shows that negative anomalies in SST over the Western Atlantic present a
negative correlation with the occurrences of tropical storms over the coast of
Florida and Cuba (see Figure 4).

Figure 4: Winter: SST Features (first left column) and their respective
statistically significant associated factors. Warm colors (red, orange) are
positive, and cold colors (purple, green) are negative. EOFs range from -0.2
to 0.2 and coefficients from -1 to 1.

TC seasonal variation is well described using SST and LHF features during
spring and summer ($w=2,3$). While these covariates cannot be used for
seasonal predictions, identifying their relationship with TC occurrence may
provide guidance for future scientific studies. During spring, we have a
significant factor that describes a positive correlation between positive
anomalies in SST over the Arctic and the number of major hurricanes in the
Gulf of Mexico and the Caribbean Sea. Also, Figure 5 shows three significant
factors that relate positive anomalies of LHF over the tropical Atlantic,
Caribbean and Gulf of Mexico as positively associated with both the occurrence
of tropical storms and their number of occurrences over the Atlantic.

During summer (Figure 8 in the Appendix), we find an opposite pattern for SST
anomalies in the Arctic, since they are negatively correlated with the number
of hurricanes and major hurricanes over the Gulf of Mexico, around Cuba and
over the Mid Atlantic. Furthermore, positive anomalies of LHF in front of the
US East Coast are positively correlated with the occurrence of hurricanes in
the same area. Finally, positive LHF anomalies over the Gulf of Mexico and the
tropical Atlantic are negatively correlated with the occurrence of major
hurricanes around the lower Antilles and with the number of tropical storms
over the Gulf of Mexico, the Caribbean and the Mid Atlantic. Lastly, during
the fall (Figure 9 in the Appendix) all the statistically significant factors
are related to LHF.

Figure 5: Spring: SST or LHF Features (first left column) and their respective
statistically significant associated fractors. Warm colors (red, orange) are
positive, and cold colors (purple, green) are negative. EOFs range from -0.2
to 0.2 and coefficients from -1 to 1.

##  6 Discussion

In this article, we develop a statistical model that allows for dynamic
variable selection. The multivariate spatial temporal smooth surface we use in
the prior reflects the fact that categories, trimesters, locations, and levels
of PHM can affect the association of SST and LHF with the response.

A simulation study shows that our model outperforms the independent and
spatial models providing an average increase of $27\%$ in the percentage of
correctly detected zeroes, in all treatments. Also, on average our model has a
MAD per $\beta_{A}$ $26\%$ lower than the competing models. In some cases,
this improvement is reflected in a $40\%$ decrease of MAD. In the data
application, we find that using 200, 400 or 800 grid boxes to construct the
SST and LHF scores has no effect on model goodness of fit (DIC), but each is
better than a traditional climate index model (Table 3). Also, our model has a
significantly lower MSE values than the competing models.

This study provides a better understanding of how observed TC variations in
space and time are related to variations in features of SST and LHF. No
previous study has attempted to relate SST and LHF features to the variation
of storms, using spatial-temporal models and variable selection techniques.
Our results show that warmer SST anomalies during winter and spring over the
mid TC main development region are associated with more major hurricanes in
the Caribbean Sea and the Gulf of Mexico. The results of this study can be
used to improve seasonal forecast models by using the proposed measures as
seasonal predictors.

##  Appendix A Descriptive Statistics for Coefficients

Table 5: Proportion of locations with estimated coefficients statistically
different from zero (credible intervals do not include zero). From those
coefficients different from zero (+) indicates that more than $60\%$ of the
locations have a positive coefficient for that specific trimester, response
and score. In a similar way, (-) indicates that more than $60\%$ of the
locations have a negative coefficient for that specific trimester, response
and score. Bold values are significant factors. Red values are those displayed
in Figures 4, 5, 8 and 9.

Resp, Level Score Sea Surface Temperature Latent Heat Flux $w=1$ $w=2$ $w=3$
$w=4$ $w=1$ $w=2$ $w=3$ $w=4$ $Y_{1}=$ TS, $p$ $r=1$ 0.02 - 0.00 0.25 - 0.00
0.00 0.00 0.00 0.00 $r=2$ 0.10 - 0.07 - 0.08 - 0.00 0.01 - 0.67 - 0.00 0.44 -
$r=3$ 0.45 - 0.22 - 0.00 - 0.00 0.21 + 0.00 0.00 0.00 $r=4$ 0.48 - 0.01 - 0.14
- 0.00 0.00 0.00 0.28 - 0.00 $Y_{1}=$ TS, $\lambda$ $r=1$ 0.00 0.01 + 0.00
0.00 0.00 0.00 0.12 + 0.00 $r=2$ 0.77 - 0.11 - 0.64 - 0.43 - 0.05 + 0.37 -
0.10 + 0.20 - $r=3$ 0.81 - 0.71 - 0.17 - 0.27 - 0.01 + 0.40 + 0.35 - 0.04 +
$r=4$ 0.01 - 0.00 0.00 0.00 0.12 - 0.00 0.36 - 0.09 - $Y_{2}=$ HU, $p$ $r=1$
0.00 0.00 0.00 0.00 0.02 - 0.00 0.16 + 0.17 - $r=2$ 0.32 - 0.00 0.00 0.00 0.00
0.00 0.00 0.00 $r=3$ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 $r=4$ 0.00 0.00
0.00 0.00 0.04 + 0.18+ 0.18 - 0.06 - $Y_{2}=$ HU, $\lambda$ $r=1$ 0.25 + 0.23
- 0.16 - 0.03 - 0.09 - 0.04 + 0.20 + 0.21 - $r=2$ 0.35 - 0.08 + 0.22 - 0.05 +
0.01 + 0.14 + 0.11 - 0.29 - $r=3$ 0.07 - 0.13 - 0.04 - 0.04 - 0.07 - 0.04 +
0.11 + 0.18 + $r=4$ 0.33 - 0.19 + 0.22 + 0.09 + 0.15 + 0.28 + 0.19 - 0.12 -
$Y_{3}=$ MH, $p$ $r=1$ 0.00 0.00 0.00 0.00 0.00 0.00 0.12 + 0.00 $r=2$ 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.32- $r=3$ 0.06 - 0.03 - 0.00 0.00 0.00 0.05-
0.25 - 0.00 $r=4$ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 $Y_{3}=$ MH,
$\lambda$ $r=1$ 0.17 + 0.04 - 0.00 0.04 + 0.08 + 0.04 + 0.35 + 0.03 + $r=2$
0.03 - 0.10 + 0.15 - 0.03 + 0.05 - 0.00 0.00 0.30 - $r=3$ 0.17 - 0.10 - 0.06 +
0.04 + 0.06 - 0.04 - 0.35 - 0.06 - $r=4$ 0.07 - 0.07 + 0.05 - 0.11 + 0.06 +
0.13 + 0.07 - 0.10 -

##  Appendix B MSE per location and per year

Figure 6: MSE per location for each response: (a) Tropical Cyclones (TC), (b)
Hurricanes (HU) and (c) Major Hurricanes (MH). Figure 7: MSE per year for each
response: Tropical Cyclones (TC), Hurricanes (HU) and Major Hurricanes (MH).

##  Appendix C Significant Factors for Summer and Fall

Figure 8: Summer: SST or LHF Features (first left column) and their respective
statistically significant associated factors. Warm colors (red, orange) are
positive, and cold colors (purple, green) are negative. EOFs range from -0.2
to 0.2 and coefficients from -1 to 1. Figure 9: Fall: LHF Features (first left
column) and their respective statistically significant associated factors.
Warm colors (red, orange) are positive, and cold colors (purple, green) are
negative. EOFs range from -0.2 to 0.2 and coefficients from -1 to 1.

## Acknowledgements

We are grateful for the support provided by the National Science Foundation’s
Research Network for Statistical Methods for Atmospheric and Oceanic Sciences,
award DMS-1107046, and National Science Foundation Awards 1406016 and 1613219.
Marcela Alfaro Córdoba is also grateful for the support provided by CONICIT
Costa Rica and the University of Costa Rica.

Supplement A Code for the MCMC sampler [url].zip file Code for the MCMC
sampler used in the Simulation.

Supplement B Supplementary Plots [url]pdf file We present supplementary plots
to the article. They include an example of the simulation settings
realizations for covariate 1, a plot of time series for each
$\xi_{\ell,r}(t,w)$, a complete plot of statistically significant factors and
a summary of significant factors for summer and fall.

## References

  * Blake and Gray (2004) [author] Blake, Eric S.E. S. and Gray, William M.W. M. (2004). Prediction of August Atlantic Basin Hurricane Activity. Weather and Forecasting 19 1044–1060. 10.1175/814.1
  * Boehm Vock et al. (2015) [author] Boehm Vock, Laura FL. F., Reich, Brian JB. J., Fuentes, MontserratM. and Dominici, FrancescaF. (2015). Spatial variable selection methods for investigating acute health effects of fine particulate matter components. Biometrics 71 167–77. 10.1111/biom.12254
  * Cai et al. (2013) [author] Cai, BoB., Lawson, Andrew BA. B., Hossain, MonirM., Choi, JungsoonJ., Kirby, Russell SR. S. and Liu, JihongJ. (2013). Bayesian semiparametric model with spatially-temporally varying coefficients selection. Statistics in medicine 32 3670–85. 10.1002/sim.5789
  * Chipman (1996) [author] Chipman, HughH. (1996). Bayesian variable selection with related predictors. Canadian Journal of Statistics 24 17–36. 10.2307/3315687
  * Elsner and Jagger (2006) [author] Elsner, James B.J. B. and Jagger, Thomas H.T. H. (2006). Prediction Models for Annual U.S. Hurricane Counts. Journal of Climate 19 2935–2952. 10.1175/JCLI3729.1
  * Elsner, Liu and Kocher (2000) [author] Elsner, James B.J. B., Liu, Kam-biuK.-b. and Kocher, BethanyB. (2000). Spatial Variations in Major U.S. Hurricane Activity: Statistics and a Physical Mechanism. Journal of Climate 13 2293–2305. 10.1175/1520-0442(2000)013¡2293:SVIMUS¿2.0.CO;2
  * Gelfand et al. (2003) [author] Gelfand, Alan EA. E., Kim, Hyon-JungH.-J., Sirmans, C. FC. F. and Banerjee, SudiptoS. (2003). Spatial Modeling With Spatially Varying Coefficient Processes. Journal of the American Statistical Association 98 387–396. 10.1198/016214503000170
  * George and McCulloch (1993) [author] George, Edward I.E. I. and McCulloch, Robert E.R. E. (1993). Variable Selection via Gibbs Sampling. Journal of the American Statistical Association 88 881–889. 10.1080/01621459.1993.10476353
  * Gray et al. (1994) [author] Gray, William M.W. M., Landsea, Christopher W.C. W., Mielke, Paul W.P. W. and Berry, Kenneth J.K. J. (1994). Predicting Atlantic Basin Seasonal Tropical Cyclone Activity by 1 June. Weather and Forecasting 9 103–115. 10.1175/1520-0434(1994)009¡0103:PABSTC¿2.0.CO;2
  * Hodges, Jagger and Elsner (2014) [author] Hodges, Robert E.R. E., Jagger, Thomas H.T. H. and Elsner, James B.J. B. (2014). The sun-hurricane connection: Diagnosing the solar impacts on hurricane frequency over the North Atlantic basin using a space–time model. Natural Hazards. 10.1007/s11069-014-1120-9
  * Ishwaran and Rao (2005) [author] Ishwaran, HemantH. and Rao, J. SunilJ. S. (2005). Spike and slab variable selection: Frequentist and Bayesian strategies. The Annals of Statistics 33 730–773. 10.1214/009053604000001147
  * Jagger, Niu and Elsner (2002) [author] Jagger, Thomas H.T. H., Niu, XufengX. and Elsner, James B.J. B. (2002). A space-time model for seasonal hurricane prediction. International Journal of Climatology 22 451–465. 10.1002/joc.755
  * Kalnay et al. (1996) [author] Kalnay, E.E., Kanamitsu, M.M., Kistler, R.R., Collins, W.W., Deaven, D.D., Gandin, L.L., Iredell, M.M., Saha, S.S., White, G.G., Woollen, J.J., Zhu, Y.Y., Leetmaa, A.A., Reynolds, R.R., Chelliah, M.M., Ebisuzaki, W.W., Higgins, W.W., Sepowiak, J.J., Mo, K. C.K. C., Ropelewski, C.C., Wang, J.J., Jenne, RoyR. and Joseph, DennisD. (1996). The NCEP/NCAR 40-Year Reanalysis Project. Bulletin of the American Meteorological Society 77 437–471. 10.1175/1520-0477(1996)077¡0437:TNYRP¿2.0.CO;2
  * Keith and Xie (2009) [author] Keith, ElinorE. and Xie, LianL. (2009). Predicting Atlantic Tropical Cyclone Seasonal Activity in April. Weather and Forecasting 24 436–455. 10.1175/2008WAF2222139.1
  * Lehmiller, Kimberlain and Elsner (1997) [author] Lehmiller, G. S.G. S., Kimberlain, T. B.T. B. and Elsner, J. B.J. B. (1997). Seasonal Prediction Models for North Atlantic Basin Hurricane Location. Monthly Weather Review 125 1780–1791. 10.1175/1520-0493(1997)125¡1780:SPMFNA¿2.0.CO;2
  * Lorenz (1956) [author] Lorenz, E NE. N. (1956). Empirical Orthogonal Functions and Statistical Weather Prediction. Technical report Statistical Forecast Project Report 1 Department of Meteorology MIT 49 1 52\.
  * Lum (2012) [author] Lum, KristianK. (2012). Bayesian variable selection for spatially dependent generalized linear models.
  * Neelon, Ghosh and Loebs (2013) [author] Neelon, BrianB., Ghosh, PulakP. and Loebs, Patrick FP. F. (2013). A Spatial Poisson Hurdle Model for Exploring Geographic Variation in Emergency Department Visits. Journal of the Royal Statistical Society. Series A, (Statistics in Society) 176 389–413. 10.1111/j.1467-985X.2012.01039.x
  * NHC (2016) [author] NHC (2016). Eye, Eyewall and Tropical Cyclone. Glossary of National Hurricane Center Terms.
  * Reich et al. (2010) [author] Reich, Brian JB. J., Fuentes, MontserratM., Herring, Amy HA. H. and Evenson, Kelly RK. R. (2010). Bayesian variable selection for multivariate spatially varying coefficient regression. Biometrics 66 772–82. 10.1111/j.1541-0420.2009.01333.x
  * Scheel et al. (2013) [author] Scheel, IdaI., Ferkingstad, EgilE., Frigessi, ArnoldoA., Haug, OlaO., Hinnerichsen, MikkelM. and Meze-Hausken, ElisabethE. (2013). A Bayesian hierarchical model with spatial variable selection: the effect of weather on insurance claims. Journal of the Royal Statistical Society: Series C (Applied Statistics) 62 85–100. 10.1111/j.1467-9876.2012.01039.x
  * Smith and Fahrmeir (2007) [author] Smith, MichaelM. and Fahrmeir, LudwigL. (2007). Spatial Bayesian Variable Selection With Application to Functional Magnetic Resonance Imaging. Journal of the American Statistical Association 102 417–431. 10.1198/016214506000001031
  * Werner and Holbrook (2011) [author] Werner, AngelikaA. and Holbrook, Neil J.N. J. (2011). A Bayesian Forecast Model of Australian Region Tropical Cyclone Formation. Journal of Climate 24 6114–6131. 10.1175/2011JCLI4231.1
  * Xie (2005) [author] Xie, LianL. (2005). The effect of Atlantic sea surface temperature dipole mode on hurricanes: Implications for the 2004 Atlantic hurricane season. Geophysical Research Letters 32 L03701. 10.1029/2004GL021702
  * Xie and Liu (2014) [author] Xie, KennyK. and Liu, BinB. (2014). An ENSO-Forecast Independent Statistical Model for the Prediction of Annual Atlantic Tropical Cyclone Frequency in April. Advances in Meteorology 2014 1–11. 10.1155/2014/248148
  * Xie et al. (2014) [author] Xie, LianL., Alfaro-Córdoba, MarcelaM., Liu, BinB. and Fuentes, MontserratM. (2014). 2014 Atlantic Tropical Cyclone Outlook Technical Report, North Carolina State University, Raleigh.


# The parity of specular Andreev reflection under mirror operation in zigzag
graphene ribbon

Yanxia Xing1,3, Jian Wang1,∗, and Qing-feng Sun2,† 1Department of Physics and
the Center of Theoretical and Computational Physics, The University of Hong
Kong, Pokfulam Road, Hong Kong, China.
2Beijing National Lab for Condensed Matter Physics and Institute of Physics,
Chinese Academy of Sciences, Beijing 100080, China.
3Department of Physics, Beijing Institute of Technology, Beijing 100081,
China.

###### Abstract

It is known that the parity of reflection amplitude can either be even or odd
under the mirror operation. Up to now, all the parities of reflection
amplitude in the one-mode energy region are even under the mirror operation.
In this paper, we give an example of odd parity for Andreev reflection (AR) in
a three-terminal graphene-supercondutor hybrid systems. We found that the
parity is even for the Andreev retroreflection (ARR) and odd for specular
Andreev reflection (SAR). We attribute this remarkable phenomenon to the
distinct topology of the band structure of graphene and the specular Andreev
reflection involving two energy bands with different parity symmetry. As a
result of odd parity of SAR, the SAR probability of a four-terminal system
with two superconducting leads (two reflection interfaces) can be zero even
when the system is asymmetric due to the quantum interference of two ARs.

###### pacs:

72.80.Vp, 74.45.+c, 73.40.-c, 74.25.F-

Since the experimental realization of grapheneref1 , it has become an exciting
arena for theoretical and technological investigations.ref2  A number of new
phenomena have been predicted and verified experimentally. For instance, in
the presence of magnetic field, it exhibits a distinctive half-integer quantum
Hall effect.ref1  Its quasi-particles obey the Dirac-like equation and have
relativistic-like behaviors.ref2  Due to the relativistic effect, the Klein
tunneling occurs where an incident electron in graphene can pass through a
potential barrier with probability one.Klein  Then a graphene p-n junction can
be used to focus Dirac electron current with a negative refractive
index.altshuler ; lens

Since good contacts between superconducting leads and graphene have been
realized experimentally,ref16  the transport study through graphene based
normal-metal-superconductor (GNS) heterojunction becomes feasible. In the
presence of a normal metal (graphene)-superconducting interface, an incoming
electron converts into a hole and a cooper pair is formed that enters the
superconductor. Due to the relativistic nature of the electron in graphene,
the electron-hole conversion can either be intraband (within conduction or
valence band) or interband (between conduction and valence bands). When the
electron-hole conversion is intraband, it corresponds to the usual Andreev
reflection (AR)ref15  or Andreev retroreflection (ARR) because the reflected
hole is along the incident direction. This ARR occurs for both relativistic
and non-relativistic electrons. When the electron-hole conversion is
interband, the reflected hole is along specular direction and a specular
Andreev reflection (SAR) takes place,Beenakker  which can lead to novel
phenomena as we will discuss below.

It is known that the parity is a fundamental quantity in physics and
reflection is a general physical phenomenon in nature. In this paper, we
discuss the parity of reflection amplitude for graphene in contact with
superconductor leads. In general, the parity of a reflection amplitude can be
either even or odd when the system is under mirror operation. However, for all
previous known reflection events, the reflection amplitudes in the one-mode
energy region have even parity under the mirror operation. It is yet to find
an odd-parity reflection event. In this paper, we found, for the first time,
that the SAR amplitude has an odd parity under the mirror operation for zigzag
graphene ribbons with even number of chains. This means that the phases of SAR
amplitude for a graphene-superconductor hybrid system and its mirror system
differ by $\pi$. We attribute this phenomenon to the unique band structure of
the graphene. Obviously this phase difference does not affect any observable
quantities for each system. When two systems couple together, however, this
$\pi$ phase manifest through quantum interference between two SARs. So this
$\pi$ phase shift has important consequences for a four terminal device with
two superconducting leads (see Fig.2(a)). When two superconducting leads are
symmetrically attached to the device, the quantum interference of the left and
right SAR leads to a destructive or constructive interference depending on
whether the phase difference of superconducting leads is zero or $\pi$.
Importantly, when two superconducting leads are asymmetrically attached to the
device, the same interference pattern occurs provided that the Dirac point
$E_{0}$ is in line with the condensate of superconducting lead. The quantum
interference between pairs of the AR can be tuned by shifting the Dirac point,
the asymmetry of the two superconducting leads, as well as the phase between
two superconducting leads. Due to the odd parity of SAR, the interference
pattern for SAR is phase contrasted to that of ARR where the parity is even.

Before doing numerical calculation, we first prove that the phases of SAR
amplitude of two systems (i) and (ii) in Fig.1(a) differ by $\pi$, i.e., the
parity of SAR is odd under mirror operation. Note that for graphene systems
electrons in valence and conduction band are usually referred as electrons and
holes, respectively. In the presence of superconducting lead the reference
point of electrons and holes is the Fermi level in the superconducting lead.
In the following, we will refer electrons (holes) as electrons above (below)
Fermi level in superconducting lead. Denote $\psi^{+}_{c}$ ($\psi_{v}^{+}$)
the wavefunction of electrons in conduction (valence) band moving in +y
direction and $\psi^{-}_{c}$ ($\psi_{v}^{-}$) in -y direction in the zigzag
graphene nanoribbon lead. It was known that under reflection
${\hat{P}}:x\rightarrow-x$, $\psi_{c}^{\pm}$ is symmetric while
$\psi_{v}^{\pm}$ is anti-symmetric if the energy of electron is in the first
transmission channelselect [see Fig.1(a)], i.e.,

$\displaystyle{\hat{P}}\psi_{c}^{\pm}(x,y)$ $\displaystyle=$
$\displaystyle\psi_{c}^{\pm}(-x,y)$
$\displaystyle{\hat{P}}\psi_{v}^{\pm}(x,y)$ $\displaystyle=$
$\displaystyle-\psi_{v}^{\pm}(-x,y).$ (1)

which is one of the unique features of zigzag edge nanoribbons with even
number of chains. Assuming the incident electron from the terminal-1, the
wavefunctions for SAR $\psi_{1,3}$ in zigzag nanoribbon lead 1 or 3 of the
system (i) can be written as

$\displaystyle\psi_{1}^{(i)}$ $\displaystyle=$
$\displaystyle\psi_{e}^{+}+r_{11}\psi_{e}^{-}+r_{11A}\psi_{h}^{-}$
$\displaystyle\psi_{3}^{(i)}$ $\displaystyle=$ $\displaystyle
t_{13}\psi_{e}^{+}+r_{13A}\psi_{h}^{+}$ (2)

where $r_{11}$ is the normal reflection amplitude, $t_{13}$ is the
transmission amplitude, $r_{11A}$ and $r_{13A}$ are the Andreev reflection
amplitudes with the reflected hole to the terminal-1 and 3, respectively.
Similarly the wavefunctions for the system (ii) are given by

$\displaystyle\psi_{1}^{(ii)}$ $\displaystyle=$
$\displaystyle\psi_{e}^{+}+{\bar{r}}_{11}\psi_{e}^{-}+{\bar{r}}_{11A}\psi_{h}^{-}$
$\displaystyle\psi_{3}^{(ii)}$ $\displaystyle=$
$\displaystyle{\bar{t}}_{13}\psi_{e}^{+}+{\bar{r}}_{13A}\psi_{h}^{+}$ (3)

Since the system (i) is related to (ii) by the reflection operator
${\hat{P}}$, we have $\psi_{\alpha}^{(i)}={\hat{P}}\psi_{\alpha}^{(ii)}$ with
$\alpha=1,3$. Note that for SAR, the electron is in the conduction band while
the hole is in the valence band, i.e., $\psi_{e}=\psi_{c}$ and
$\psi_{h}=\psi_{v}$. From this relation together with Eqs.(1), (2), and (3),
we obtain

$\displaystyle r_{11A}=-{\bar{r}}_{11A},~{}~{}~{}r_{13A}=-{\bar{r}}_{13A}$
$\displaystyle r_{11}={\bar{r}}_{11},~{}~{}~{}t_{13}={\bar{t}}_{13}$ (4)

Note that the origin of this $\pi$ phase shift (odd parity) is the interband
conversion from the electron to the hole. Therefore the $\pi$ phase shift does
not occur for ARR since it involves only intraband conversion. Now we verify
this statement numerically using a tight-binding model (see below for detailed
description of the model and numerical procedure). The numerical results of AR
probability $R_{11A(13A)}=|r_{11A(13A)}|^{2}$ for two systems are shown in
Fig.1(b). As expected the AR probability are exactly the same for two systems.
However, the phase of AR amplitudes $r_{11A(13A)}$ denoted as
$\Phi^{i,ii}_{11(13)}$ are different. It is shown in Fig.1(c) and Fig.1(d)
that ARR amplitudes ($|E_{0}|>|E_{F}|$, with $|E_{F}|=0.5$) are the same for
two systems in Fig.1(a) while the SAR amplitudes ($|E_{0}|<|E_{F}|$) have a
$\pi$ phase shift. It confirms the odd parity for interband electron-hole
conversion, which comes from the distinct topology of the band structure of
graphene.

To see the consequence of the odd parity of SAR, we examine a symmetric four-
terminal device with two superconducting leads depicted in Fig.2(a) (by
setting asymmetry $\delta N=0$ and phase difference $\delta\phi=0$). For this
system, two beams from terminal-1 has a $\pi$ phase shift due to odd parity of
SAR and interferes destructively at terminal-3 giving rise to a vanishing SAR
coefficient. However, we can arrive the same conclusion using symmetry
argument as follows. Since the system is symmetric with respect to $x=0$, we
must have $r_{13A}={\bar{r}}_{13A}$ when the reflection operation along
x-direction is applied. While from Eq.(4), $r_{13A}=-{\bar{r}}_{13A}$. So the
AR probability $R_{13A}=|r_{13A}|^{2}$ for SAR can also be zero from symmetry
point of view.cheng  Therefore we conclude that the symmetric device can not
be used to test the odd parity of SAR. In the following, we demonstrate that
due to the $\pi$ phase shift the destructive interference still occurs in a
four-probe devices with two superconducting leads attached asymmetrically and
hence can be used to test the odd parity of SAR.

For this purpose, we consider an asymmetric four-terminal device consisting of
a zigzag graphene ribbon with two superconducting leads as shown in Fig.2(a).
The Hamiltonian of the graphene isref19  $H_{0}=\sum_{\bf i}\epsilon_{\bf
i}a^{\dagger}_{\bf i}a_{\bf i}-\sum_{<{\bf ij}>}ta_{\bf i}^{\dagger}a_{\bf
j}$. Here $a_{\bf i}$ and $a_{\bf i}^{\dagger}$ are the annihilation and
creation operators at site ${\bf i}$, $\epsilon_{\bf i}$ is the on-site energy
which can be controlled experimentally by the gate voltageref1 , and the
hopping constant $t=2.75eV$ represents the nearest carbon bond energy. The
pair potential (energy gap) of superconducting terminal-$\beta$ with
$\beta=2,4$ is $\tilde{\Delta}_{\beta}=\Delta_{\beta}e^{i\varphi_{\beta}}$
with $\Delta_{2}=\Delta_{4}=\Delta\simeq 1meV$. In numerical
calculations,cheng  we fix Fermi energy $E_{F}$ and tune the Dirac point
$E_{0}$. We have used $\Delta$ as the energy unit.

Now we study the interference between two ARs from GNS junctions as shown in
Fig.2(a) in which two superconducting leads 2 and 4 are asymmetrically
attached to the zigzag nanoribbon. The horizontal distance $\delta N$ between
two GNS junctions measures the asymmetry of two GNS junctions. The scattering
process can be qualitatively understood as follows. For simplicity, we assume
$\phi_{2}=\phi_{4}$ for the moment. As shown schematically in Fig.2(a), for
SAR the particle-like electrons in terminal-1 split into two beams and are
scattered separately by two GNS junctions (green horizontal lines) as holes
that finally recombine at terminal-3. We examine the total phase accumulated
for each beam that involves the following three processes. Before reaching the
first GNS junction (denoted by the left vertical green line) two beams of
electrons propagate with the same momentum $k_{x}$. After reaching the second
GNS junction (denoted by the right vertical green line) two beams of holes
also propagate with the same momentum $k^{\prime}_{x}$. Obviously phases
accumulated in the above two processes for both beams are the same. Between
them two beams propagate with different momenta $k_{x}$ and $k^{\prime}_{x}$.
Hence the phase difference between two beams is
$\phi=(k_{x}-k^{\prime}_{x})\delta x$ with $\delta x=b\delta N$, where
$b=\sqrt{3}a$ and $a$ the lattice constant. This phase difference can be tuned
by varying the Dirac point $E_{0}$ or the asymmetry $\delta N$ giving rise to
a complicated interference pattern (see Fig.2). In particular, this phase
difference can be zero if $(k_{x}-k^{\prime}_{x})=0$ (i.e.,$E_{0}=0$) or
$\delta N=0$. In general, the total phase difference is
$\phi=(k_{x}-k^{\prime}_{x})\delta x+\phi_{2}-\phi_{4}$.

Interference pattern of AR probability $R_{13A}$ for system depicted in
Fig.2(a) with pair potential phase difference of two superconductors
$\delta\varphi=0$ and $\pi$ ($\delta\varphi\equiv\varphi_{2}-\varphi_{4}$) are
then plotted in Fig.2(b) and (c), respectively. For Fig.2(b) following
observations are in order: (1) For the geometrically symmetric system ($\delta
N=0$), the interference is always destructive with zero $R_{13A}$ as long as
$|E_{0}|<|E_{F}|$.cheng  Clearly this is due to the $\pi$ phase shift depicted
in Fig.1(d) and is consistent with the band selection rule.select  (2) When
Dirac point $E_{0}$ is in line with the condensate energy of the
superconductor, i.e., when $E_{0}=0$, $R_{13A}$ is again zero no matter what
value $\delta N$ assumes. This means that there is a completely destructive
interference between two beams scattered by two GNS junctions attached
asymmetrically to the graphene nano-ribbon. This behavior can be understood as
follows. When $E_{0}=0$ the incoming electron and reflected hole have the same
propagating momentum $k_{x}$ and thus path 1 and 2 in Fig.2(a) experience the
same quantum phase $k_{x}\delta x$ except at the superconducting leads. Hence
the total phase difference is only due to the $\pi$ phase shift between two
SARs. (3) $R_{13A}$ is an even function of Dirac point $E_{0}$ because of the
electron-hole symmetry in graphene. Due to the geometric symmetry, $R_{13A}$
is also an even function of asymmetry $\delta N$. (4) For nonzero $E_{F}$, the
closer the Dirac point $E_{0}$ to $E_{F}$, the more rapidly $R_{13A}$
oscillates as we vary $\delta N$. This is because the difference of
propagating momentum $k_{x}-k^{\prime}_{x}$ increases monotonically as $E_{0}$
approaches to $E_{F}$. (5) When $E_{0}$ is in the vicinity of $E_{F}$,
$R_{13A}$ can reach 0.9 which is much larger than that when $|E_{0}|>|E_{F}|$.
This is because when $E_{F}$ is very close to $E_{0}$, the edge states of
zigzag ribbon begin to contribute, then electron is easier to be scattered by
two GNS junctions located also at edges of zigzag ribbon. Considering the
pseudo-spin conservation, large $R_{13A}$ is always found in the region of
$|E_{0}|<|E_{F}|$, i.e., the SAR region. (6) There is an overall fine
oscillation with a period of $\delta N=3b$. Similar behavior was also found in
zigzag ribbons with a p-n junction where the conductance is determined by the
relative displacement $\delta$ along the p-n junction.Beenakker1  In Fig.2(c)
with the superconducting phase difference $\delta\varphi=\pi$, we see that the
interference pattern is contrary to $\delta\varphi=0$ [Fig.2(b)] where the
constructive interference becomes destructive and vice versa.

To further analyze the interference pattern, we plot in Fig.3(a) the total
$R_{13A}$ vs Dirac point $E_{0}$ for different asymmetry $\delta N$ with the
phase difference between two superconducting leads $\delta\varphi=0$ [main
panel of Fig.3(a)] or $\delta\phi=\pi$ [inset of Fig.3(a)]. Clearly the
interference (oscillatory) pattern occurs only for asymmetric systems ($\delta
N\neq 0$) with oscillation frequency proportional to $\delta N$. When pair
potential phase difference $\delta\varphi=\pi$ is introduced, the interference
pattern reverses, and $R_{13A}$ with $\delta N=0$ becomes the envelop function
of $R_{13A}$ for all nonzero $\delta N$. In Fig.3(b) we plot $R_{13A}$ vs
$\delta N$ for different widths $W$ of nanoribbon. It is shown clearly that
$R_{13A}$ is a periodic function of $\delta N$ with larger periodicity for
larger $W$. In the inset of Fig.3(b) we plot this period versus the width for
different $E_{0}$. The period $P$ is obtained in two ways: (1). from the
expression $P=2\pi/(k_{x}-k^{\prime}_{x})$ where the momenta $k_{x}$ and
$k_{x}^{\prime}$ can be obtained from the band structure for a given $E_{0}$
(black symbols). (2). directly from main panel of Fig.3(b) (red solid circle).
From the inset, it clearly shows that two periods are exactly the same giving
strong evidence that the interference pattern of AR probability are indeed
from two reflected hole beams.

Finally, the interference pattern of AR probability $R_{11A}$ is also studied
(not shown). We found that only ARR probability $R_{11A}$
($|E_{0}|>E_{F}=0.2\Delta$) exhibits interference pattern. We note that since
there is no $\pi$ phase shift involved in ARR, when $\delta N=0$ reflected
electrons through two GNS junctions interfere constructively when
$\delta\varphi=0$ and destructively when $\delta\varphi=\pi$ which is in
contrast to SAR in Fig.2. In fact, interference patterns of SAR and ARR are
always phase contrast not only for $\delta N=0$ but also for all other $\delta
N$.

To test the odd parity of SAR experimentally, it relies on the fabrication of
high quality zigzag graphene nanoribbons. It has been achieved by several
laboratories using different methods last year including the method to unzip
the multi-walled carbon nanotube (CNT),ref13  the anisotropic etching by
thermally activated nickel nanoparticles,ref14  and use reconstruction of the
edge to make zigzag graphene nanoribbons.experiment  In view of the above
experimental breakthrough, we expect that the setup to test our predicted
phenomenon can be realized experimentally.

To reduce the experimental challenge, we have considered an unzipped CNT
device, i.e., (n,n) CNT-zigzag graphene-(n,n) CNT, obtained by unzipping a few
unit cells in the central part of an armchair CNT which has been achieved
experimentally.ref13  For this system, the wavefunction in the armchair CNT
has the same symmetry as that of the zigzag graphene ribbon. Following the
same procedure leading to Eq.(4), we have shown that the unzipped CNT in
contact with a superconducting lead has the odd parity under mirror operation.
Similar conclusions drawn from GNS can be obtained for unzipped CNT with two
superconducting leads.

In conclusion, up to now, the parity of reflection amplitude was found to be
even under the mirror operation. Here we have provided an example of odd
parity for the reflection amplitude, the SAR amplitude in the zigzag graphene-
superconductor hybrid system. This odd parity is due to the combination of
unique band structure of the graphene and the electron-hole conversion
involving two energy bands with different parity symmetry. The signature of
odd parity of SAR can be found from the quantum constructive interference in a
four terminal system with two superconducting leads attached asymmetrically.
Furthermore, the interference pattern due to odd parity of SAR is phase
contrasted to that of ARR where the parity is even.

Acknowledgement We gratefully acknowledge the financial support from a RGC
grant (HKU 705409P) from the Government of HKSAR and from NSF-China under
Grant Nos.10974236 and 10821403.

## References

  * (1) ∗e-mail<EMAIL_ADDRESS>
  * (2) †e-mail<EMAIL_ADDRESS>
  * (3) K. S. Novoselov et al, Science 306, 666 (2004); K. S. Novoselov et al, Nature (London) 438, 197 (2005); Y. Zhang et al, Nature (London) 438, 201 (2005).
  * (4) C. W. J. Beenakker, Rev. Mod. Phys. 80, 1337 (2008); A.H. Castro Neto et al, Rev. Mod. Phys. 81, 109 (2009); A. Rycerz et al, Nature Phys. 3, 172 (2007).
  * (5) M. I. Katsnelson et al, Nature Phys. 2, 620 (2006); A.-F. Young et al, Nature Phys. 5, 222 (2009); A. V. Shytov et al, Phys. Rev. Lett. 101, 156804 (2008).
  * (6) V. V. Cheianov et al, Science 315, 1252 (2007).
  * (7) Y. Xing, J. Wang and Q.-F. Sun, Phys. Rev. B 81, 165425 (2010).
  * (8) H. B. Heersche et al, Nature 446, 56 (2007); F. Miao et al, Science 317, 1530 (2007).
  * (9) A. F. Andreev, Sov. Phys. JETP 19, 1228 (1964).
  * (10) C. W. J. Beenakker, Phys. Rev. Lett. bf 97, 067007 (2006).
  * (11) J. Nakabayashi et al, Phys. Rev. Lett 102, 066803 (2009).
  * (12) S. Cheng et al, Phys. Rev. Lett. 103, 167003 (2009); Q.-F. Sun and X.C. Xie, J. Phys.: Condens. Matter 21, 344204 (2009).
  * (13) D. N. Sheng et al, Phys. Rev. B 73, 233406 (2006); Z. Qiao and J. Wang, Nanotechnology 18, 435402 (2007). W. Long et al, Phys. Rev. Lett. 101, 166806 (2008); J. Li and S.-Q. Shen, Phys. Rev. B 78, 205308 (2008).
  * (14) A. R. Akhmerov et al, Phys. Rev. B 77, 205416 (2008).
  * (15) L. Jiao et al, Nature (London), 458, 877 (2009).
  * (16) L. C. Campos et al, Nano. Lett. 9, 2600 (2009).
  * (17) Ç. Ö. Girit et al, ibid, 323 1705 (2009).

Figure 1:  (Color online) Panel (a): zigzag ribbons with even number of chains
(gray honeycomb) attached by a superconducting lead on the left and right
(orange honey comb), respectively. For SAR the incoming electrons (red arrow)
are scattered by the GNS junction (green solid line) as holes (blue arrow).
The corresponding wave functions at sublattice “A” (solid circle) and “B”
(hollow circle) for the lowest subband in conduction band (bottom) and the
highest subband in valence band (top) are shown schematically. Panel (b): AR
probability from terminal-1 to terminal-1 $R_{11A}$ and to terminal-3
$R_{13A}$ vs. Dirac point $E_{0}$. Panels (c) and (d): AR phase
$\Phi^{i,ii}_{11}$ (c) and $\Phi^{i,ii}_{13}$ (d) of two systems in panel (a)
and their phase different $\Phi^{i}_{11(13)}-\Phi^{ii}_{11(13)}$ vs. $E_{0}$.
Figure 2:  (Color online) Panel (a): sketch of AR interferometer in which the
zigzag ribbon is asymmetrically attached by two superconductor lead-2 and 4.
Electrons in terminal-1 can be Andreev reflected into terminal-3 by either top
or bottom GNS junction (horizontal green lines). Panel(b) and (c): the contour
plot of $R_{13A}$ vs Dirac point $E_{0}$ and asymmetry $\delta N$. The phase
difference of two superconductor leads $\delta\varphi$ is zero in panel (b)
and $\pi$ in panel (c). The other parameters: Fermi energy $E_{f}=0.8$, number
of chains in zigzag ribbon $N=40$ corresponding to width $60a$, the width of
superconductor lead $W_{S}=10b$, where $b=\sqrt{3}a$. Figure 3:  (Color
online) Panel (a): With fixed Fermi level $E_{F}=0.8$, total AR probability
$R_{13A}$ vs Dirac point $E_{0}$ for different asymmetry $\delta N$. In the
main panel $\delta\varphi=0$, and while $\delta\varphi=\pi$ in the inset.
Panel (b): $R_{13A}$ vs asymmetry $\delta N$ with $E_{0}=0.3t$ for different
width $W$ from $10\times 3a$ to $38\times 3a$ with the interval $2\times 3a$
along the black arrow. Inset panel: the hollow signs are the period $P$
obtained from the main panel and the solid red circles are the period $P$ from
the energy band with the expression $P=2\pi/(k_{x}-k_{x}^{\prime})$. The other
parameters: $\delta\varphi=0$, $E_{F}=0.8$.