\section{Introduction}

Run I of the Large Hadron Collider (LHC)  was marked by the discovery and initial characterization of the Higgs boson.
The comparison of Standard Model (SM) predictions with data from Run I of the LHC was limited by the statistical precision of the experimental data.  This will no longer be the case during Run II, and systematic errors will dominate.  The largest systematic error currently hindering our understanding of Higgs properties is the theoretical precision of the SM predictions.
This is the case for the well-measured di-boson modes~\cite{Higgserrors} which dominate the overall signal-strength determination.
The theoretical uncertainties must be reduced in order to sharpen our understanding of the mechanism of electroweak symmetry breaking in Nature. Calculations through next-to-next-to-leading order (NNLO) in perturbative QCD have become increasingly necessary to match the
precision of LHC measurements. In particular improvements in both the overall production rate of the Higgs boson and in the modeling of its kinematic distributions are needed to match the expected experimental precision of Run II.

NNLO calculations for scattering processes with final-state jets at hadron colliders possess a complex singularity structure.
Partial results for $2 \to 2$  processes such as inclusive jet production~\cite{Ridder:2013mf}, Higgs+jet production~\cite{Boughezal:2013uia,Boughezal:2015dra} and Z+jet production~\cite{Ridder:2016rzm} are available (full results for Z+j using antenna subtraction will be available soon~\cite{Ridder:upcoming}). More recently complete results for Higgs+jet production~\cite{Boughezal:2015aha}, $W^{+}$+jet
production~\cite{Boughezal:2015dva} and $Z$+jet production~\cite{Boughezal:2015ded} were achieved using N-jettiness subtraction method~\cite{Boughezal:2015dva,Boughezal:2015eha,Stewart:2010tn,Gaunt:2015pea}. The purpose of this contribution is to summarize the method and the results achieved for these processes.

\section{Description of N-jettiness subtraction}
\label{sec:jettinessDef}

We review here the salient features of the N-jettiness subtraction scheme for NNLO calculations, which was recently
introduced in the context of the NNLO computation of $W^{+}$ boson and Higgs boson production in association with
a jet~\cite{Boughezal:2015aha,Boughezal:2015dva,Boughezal:2015ded}.  We begin with the definition of $N$-jettiness$,{\cal T}_N$, a global event shape variable designed to veto final-state jets~\cite{Stewart:2010tn}:
\begin{equation}
\label{eq:taudef}
{\cal T}_N = \sum_k \rm{min}_i \left\{ \frac{2 p_i \cdot q_k}{Q_i} \right \}.
\end{equation}
The subscript $N$ denotes the number of jets desired in the final state, and is an input to the measurement.  For the Higgs + jet,
$W^{+} $+ jet and $Z $+ jet processes considered here, we have $N=1$.  Values of ${\cal T}_1$ near zero indicate a final state containing a single narrow energy deposition, while larger values denote a final state containing two or more well-separated energy depositions.
The $p_i$ are light-like vectors for each of the initial-state beams and final-state jets in the problem, while the $q_k$ denote the four-momenta of any final-state radiation.  The $Q_i$ are dimensionful variables that characterize the hardness of the beam-jets and final-state jets.  We set $Q_i = 2 E_i$, twice the energy of each jet.
The cross section for ${\cal T}_N$ less than some value ${\cal T}_N^{cut}$ can be expressed in the form~\cite{Stewart:2010pd,Stewart:2009yx}
\begin{equation} \label{eq:fact}
 \sigma({\cal T}_N < {\cal T}_N^{cut})=  \int H \otimes B \otimes B \otimes S \otimes   \left[ \prod_n^{N} J_n \right] +\cdots .
\end{equation}
The function $H$ contains the virtual corrections to the process.  The beam function $B$ encodes the effect of radiation collinear to one of the two initial beams.  It can be written as a perturbative matching coefficient convoluted with a parton distribution function.  $S$ describes the soft radiation, while $J_n$ contains the radiation collinear to a final-state jet.  The ellipsis denotes power-suppressed terms which become negligible for ${\cal T}_N \ll Q_i$. Each of these functions obeys a renormalization-group equation that allows logarithms of ${\cal T}_N$ to be resummed.  If this expression is instead expanded to fixed-order in the strong coupling constant, it reproduces the cross section for low ${\cal T}_N$. The derivation of this factorization theorem in the small-${\cal T}_N$ limit relies upon the machinery of Soft-Collinear Effective Theory~\cite{Bauer:2000ew}.

The basic idea behind N-jettiness subtraction is that ${\cal T}_N$ fully captures the singularity structure of QCD amplitudes with final-state partons.  This allows us to calculate the NNLO corrections to processes such as Higgs + jet, $W^{+} $+ jet and Z + jet in the following way.  We divide the phase space according to whether ${\cal T}_N$ is greater than or less than ${\cal T}_N^{cut}$.  For ${\cal T}_N > {\cal T}_N^{cut}$ there are at least two hard partons in the final state, since all singularities are controlled by N-jettiness.  This region of phase space can therefore be obtained from, for example, a NLO calculation of Higgs production in association with two jets in the case where the born process is Higgs + jet.  Below ${\cal T}_N^{cut}$, the cross section is given by the factorization theorem of Eq.~(\ref{eq:fact}) expanded to second order in the strong coupling constant.  As long as ${\cal T}_N^{cut}$ is smaller than any other kinematic invariant in the problem, power corrections below the cutoff are unimportant.

All ingredients of Eq.~(\ref{eq:fact}) are known to the appropriate order to describe the low ${\cal T}_N$ region through second order in the strong coupling constant.  The two-loop virtual corrections for the processes discussed here are known~\cite{Gehrmann:2011aa,Gehrmann:2011ab}.
The beam functions are known through NNLO~\cite{Gaunt:2014xga,Gaunt:2014cfa}, as are the jet functions~\cite{Becher:2006qw,Becher:2010pd} and soft function~\cite{Boughezal:2015eha}.  It is therefore possible to combine this information to provide the full NNLO calculation of Higgs + jet
and $W^{+}/Z $+ jet.

A full NNLO calculation requires as well the high ${\cal T}_N$ region above ${\cal T}_N^{cut}$.  However, a finite value of ${\cal T}_N$ implies that there are actually
$N+1$ resolved partons in the final state.  This is the crucial observation; ${\cal T}_N$ completely describes the singularity structure of QCD amplitudes that contain $N$ final-state partons at leading order.  The high ${\cal T}_N$ region of phase space is therefore described by a NLO calculation with $N+1$ jets.  We must choose ${\cal T}_N^{cut}$ much smaller than any other kinematical invariant in the problem in order to a
void power corrections to Eq.~(\ref{eq:fact}) below the cutoff.

\section{Numerical Results} \label{sec:numerics}

We now present some numerical results for Higgs, $W^{+}$ and Z production in association with a jet. For validation checks of the results presented here we refer the reader to the detailed description in~\cite{Boughezal:2015aha,Boughezal:2015dva,Boughezal:2015ded}.  We focus on 8 TeV proton-proton collisions.  Jets are reconstructed using the anti-$k_T$ algorithm~\cite{Cacciari:2008gp} with $R=0.5$.  For the Higgs+jet process
we show results using the NNPDF~\cite{Ball:2012cx} parton distribution functions (PDFs), for the $W^{+}$+jet we use
CT10 PDFs~\cite{Gao:2013xoa} while for Z + jet we use CT14 PDFs~\cite{Dulat:2015mca}.  We use the perturbative order of the PDFs that is consistent with the partonic cross section under consideration: LO PDFs with LO partonic cross sections, NLO PDFs with NLO partonic cross sections, and NNLO PDFs with NNLO partonic cross sections.  We set the renormalization and factorization scales equal to the mass of the Higgs boson, $\mu_R = \mu_F = m_H$ for Higgs+jet, $\mu=M_W$ for $W^{+}$+jet and $\mu=\mu_0=\sqrt{m_{ll}^2+\sum (p_{T}^{jet})^2}$  for Z + jet. For the latter scale choice the sum is over the transverse momenta of all final-state jets, and $m_{ll}$ is the invariant mass of the di-lepton pair arising from the $Z$-boson decay.  To estimate the residual theoretical error, we vary these scales simultaneously around the central value by a factor of two. We set the mass of the Higgs boson as $m_H = 125$ GeV. We impose the following cuts on the final-state jet: $p_T^{jet} > 30$ GeV, $|\eta_{jet}|<2.4$ for Higgs+jet and $|\eta_{jet}|<2.5$ for $W^{+}$+jet. For the Z + jet process we show a plot for the dependence of the ratio $\sigma_{NNLO}/\sigma_{NLO}$
on the power corrections as a function of $\tau_1^{cut}$ and refer the reader to the corresponding paper for more phenomenological studies~\cite{Boughezal:2015ded}.

We begin by showing few distributions in the Higgs plus jet production case.  In Fig.~\ref{fig:pTjet} we show the transverse momentum distribution of the leading jet. There is a shape dependence to the corrections, with the $K$-factor decreasing as $p_T^{jet}$ is increased.
This trend is visible when going from LO to NLO in perturbation theory, and also when going from NLO to NNLO.  We note that the NNLO
result is entirely contained within the NLO scale-variation band.  The shape dependence and magnitude of the NNLO corrections for the $p_T^{jet}$ distribution are in agreement with the results of Ref.~\cite{Boughezal:2015dra}, obtained using sector-improved residue subtraction
scheme~\cite{Czakon:2010td,Boughezal:2011jf}. In Fig.~\ref{fig:pThiggs} we show the transverse momentum of the Higgs boson. The NLO corrections range from 40\% to 120\% near $p_{T}^H=60$ GeV, depending on the scale choice.  The magnitude of this correction decreases as the transverse momentum of the Higgs increases.  The NNLO corrections are more mild, reaching only 20\% at most for the central scale choice $\mu=m_H$.  They also decrease slightly as the transverse momentum of the Higgs increases.  The shape dependence and magnitude of the NNLO corrections for the $p_T^{H}$ distribution are in agreement with the results of Ref.~\cite{Boughezal:2015dra}.  We note that we have combined the two bins closest to the boundary $p_T^H=30$ GeV to avoid the well-known Sudakov shoulder effect~\cite{Catani:1997xc}.

\begin{figure}[!tbp]
  \centering
  \begin{minipage}[b]{0.47\textwidth}
    \includegraphics[width=\textwidth]{h1j_pTJ.pdf}
    \caption{The transverse momentum of the leading jet for the Higgs+jet process at LO, NLO, and NNLO in the strong coupling constant. The lower inset shows the ratios of NLO over LO cross sections, and NNLO over NLO cross sections.  Both shaded regions in the upper panel and the lower inset indicate the scale-variation errors.} \label{fig:pTjet}
  \end{minipage}
  \hfill
  \begin{minipage}[b]{0.47\textwidth}
    \includegraphics[width=\textwidth]{h1j_pTH.pdf}
    \caption{The transverse momentum of the Higgs boson at LO, NLO, and NNLO in the strong coupling constant. The lower inset shows the ratios of NLO over LO cross sections, and NNLO over NLO cross sections.  Both shaded regions in the upper panel and the lower inset indicate the scale-variation errors.} \label{fig:pThiggs}
  \end{minipage}
\end{figure}

In Fig.~\ref{fig:ptspectrum} we show the transverse momentum spectrum of the leading jet for $W^{+}$ + jet at LO, NLO and
NNLO in perturbation theory.  The ratios of the NLO cross section over the LO result, as well as the NNLO cross section over the NLO one,
are shown in the lower inset.  The shaded bands in the upper inset indicate the theoretical errors at each order estimated by varying
the renormalization and factorization scales by a factor of two around their central value, as do the shaded regions in the lower inset.  In the
 lower inset we have shown the results for both ${\cal T}_N^{cut}=0.05$ GeV and ${\cal T}_N^{cut}=0.07$ GeV, for the scale choice
 $\mu=2 M_W$, to demonstrate the ${\cal T}_N^{cut}$ independence in every bin studied.  The NLO corrections are large and positive for this scale choice, increasing the cross section by 40\% at $p_T^{jet}=40 $ GeV and by nearly a factor of two at $p_T^{jet}=180$ GeV.  The scale variation at NLO reaches approximately $\pm 20\%$ for $p_T^{jet}=180$ GeV.  The shift when going from NLO to NNLO is much more mild, giving only a percent-level decrease of the cross section that varies only slightly as $p_T^{jet}$ is increased.  The scale variation at NNLO is at the percent level and is nearly invisible on this plot.

\begin{figure}[!tbp]
  \centering
  \begin{minipage}[b]{0.47\textwidth}
    \includegraphics[width=\textwidth]{w1j_pTJ_tauband.pdf}
    \caption{The transverse momentum spectrum of the leading jet for $W^{+}$ + jet at LO, NLO and NNLO in perturbation theory.  The bands indicate the estimated theoretical error. The lower inset shows the ratios of the NLO over the LO cross section, and the NNLO over the NLO cross section.  Both shaded regions in the upper panel and the lower inset indicate the scale-variation errors.  The dashed and solid black lines in the lower inset respectively show the distribution for ${\cal T}_1^{cut}=0.05$ GeV and ${\cal T}_1^{cut}=0.07$ GeV, for the scale choice $\mu=2 M_W$.} \label{fig:ptspectrum}
  \end{minipage}
  \hfill
  \begin{minipage}[b]{0.475\textwidth}
    \includegraphics[width=\textwidth]{w1j_pTW_tauband.pdf}
    \caption{The transverse momentum spectrum of the $W^{+}$-boson at LO, NLO and NNLO in perturbation theory.  The bands indicate the estimated theoretical error. The lower inset shows the ratios of the NLO over the LO cross section, and the NNLO over the NLO cross section. Both shaded regions in the upper panel and the lower inset indicate the scale-variation errors.  The dashed and solid black lines in the lower inset respectively show the distribution for ${\cal T}_1^{cut}=0.05$ GeV and ${\cal T}_1^{cut}=0.07$ GeV, for the scale choice $\mu=2 M_W$.}  \label{fig:ptWspectrum}
  \end{minipage}
\end{figure}

The transverse momentum spectrum of the $W$-boson is shown in Fig.~\ref{fig:ptWspectrum}.  The NLO
corrections are again 40\% for $p_T^W \geq 50$ GeV with a sizable scale dependence, while the NNLO corrections
are flat in this region and decrease the cross section by a small amount.  The phase-space region
$p_T^{W}<30$ GeV only opens up at NLO, leading to a different pattern of corrections for these transverse momentum values.
The instability of the perturbative series in the bins closest to the boundary $p_T^W=30$ GeV is caused by the well-known Sudakov-shoulder
effect~\cite{Catani:1997xc}.

Finally, in Fig.~\ref{fig:taucheck} we show the dependence of the sum of the cross sections above and below the N-jettiness cutoff
${\cal T}_N^{cut}$ and the effect of power corrections. The validation is done for the ratio $\sigma_{\text{NNLO}}/\sigma_{\text{NLO}}$ in
 13 TeV proton-proton collisions.  We have checked that the NLO cross section obtained with $N$-jettiness subtraction agrees with
the result obtained with standard techniques.  These cross sections are obtained using CT14 PDFs at the same order in perturbation theory as the partonic cross section, and contain the following fiducial cuts on the leading final-state jet and the two leptons from CMS~\cite{Khachatryan:2014zya}: $p_T^{jet} > 30$ GeV, $|\eta_{jet}|<2.4$, $p_T^{l} > 20$ GeV, $|\eta_{l}|<2.4$ and $71\, \text{GeV} ~< m_{ll} ~< 111 \, \text{GeV}$. The ATLAS analysis is similar but with slightly different cuts.  We reconstruct jets using the anti-$k_T$ algorithm~\cite{Cacciari:2008gp} with $R=0.5$.
A dynamical scale $\mu_0=\sqrt{m_{ll}^2+\sum p_{T}^{jet,2}}$ is chosen to describe this process, where the sum is over the transverse
 momenta of all final-state jets, and $m_{ll}$ the invariant mass of the di-lepton pair arising from the $Z$-boson decay. In this
validation plot we have set the renormalization and factorization scales to $\mu_R = \mu_F = 2 \times \mu_0$; since the corrections are larger for this scale choice, it is easier to illustrate the important aspects of the ${\cal T}_1^{cut}$ variation.

\begin{figure}[h]
\centering
\includegraphics[width=3.3in]{validation.pdf}%
\caption{Plot of the NNLO cross section over the NLO result, $\sigma_{\text{NNLO}}/\sigma_{\text{NLO}}$, as a function of ${\cal T}_1^{cut}$, for the scale choice $\mu= 2 \times \mu_0$.  The vertical bars accompanying each point indicate the integration errors.} \label{fig:taucheck}
\end{figure}

A few features can be seen in Fig.~\ref{fig:taucheck}.  First, in the region ${\cal T}_1^{cut} < 0.2$ GeV the result becomes independent of the particular value of the cut chosen within the numerical errors.  The NNLO correction for $\mu = 2 \times \mu_0$ corresponds to a $+3\%$ shift in the cross section.  The plot makes clear that we have numerical control over the NNLO cross section to the per-mille level, completely sufficient for phenomenological predictions.  We observe an approximately linear dependence of $\sigma_{\text{NNLO}}$ on ${\rm ln}\left({\cal T}_1^{cut}\right)$ in the region $ 0.2 \, {\rm GeV} < {\cal T}_1^{cut} < 0.5 \, {\rm GeV}$, indicating the onset of the neglected power corrections.  These power corrections have the form $({\cal T}_N / Q) \, {\rm ln}^n({\cal T}_N / Q)$, where $n \leq 3$ at NNLO~\cite{Gaunt:2015pea} and $Q$ is a hard scale such as $p_T^{jet}$.

\section{Conclusions}

We have presented in this proceedings the complete NNLO calculation of $W^{+}/Z$ + jet and Higgs boson production in association with a jet in hadronic collisions.  To perform this computation we have used a new subtraction scheme based on the $N$-jettiness event-shape variable  ${\cal T}_N$.
We will further study the phenomenological impact of our NNLO result in future work, including the prediction for the exclusive one-jet bin, where an intricate interplay between various sources of higher-order corrections was recently pointed out~\cite{Boughezal:2015oga}.
We will in addition use the full NNLO result to improve upon the resummation of jet-veto logarithms that occur when Higgs production is measured in exclusive jet bins.  Previous work has indicated that this resummation has an important effect in reducing the theoretical uncertainties that plague the predictions for exclusive jet multiplicities~\cite{Boughezal:2013oha}.


\section{Introduction}
	The experimentally measured \cite{Nat460_1101, Sci325_178}  Fermi surfaces for the surface states of flat semi-infinite slabs of the three-dimensional topological insulator \ce{Bi2Te3} differs at moderate Fermi energies from the circular cross sections of Dirac cones predicted by the Dirac fermion Hamiltonian $H = v(\vec{p}\times\hat{z})\cdot\vec{\sigma}$ where $\hat{z}$ is the surface normal to the slabs. In order to reproduce the experimentally observed Fermi surfaces, Fu introduced a hexagonal warping term to the Hamiltonian based on the symmetries of the underlying crystal structure \cite{PRL103_266801}. The modified effective Hamiltonian for the surface states of the flat slab now reads as
\begin{equation}
	H = v(\vec{p}\times\hat{z})\cdot\vec{\sigma} + \lambda \sigma_z (k_x^3-3k_xk_y^2).  \label{ham}
\end{equation}

	Liu \textit{et al} subsequently derived a four-band effective Hamiltonian \cite{PRB82_045122} which describes both the bulk as well as surface states near the Dirac point of the \ce{Bi2Se3} family of TIs, of which \ce{Bi2Te3} is a member. They noted that Fu's hexagonal warping term can also obtained from the model Hamiltonian.

	We showed in an earlier work \cite{SciRep4_5062} that the additional hexagonal warping term in a flat \ce{Bi2Te3} slab leads to the distortion of the Fermi surface and the opening of a bandgap by an \textit{in-plane} magnetization, both of which would not have occured in the absence of the hexagonal warping term. These and the displacement of the Fermi surfaces by the in-plane magnetization lead to a complex profile of the transmission from an unmagnetized source segment to a magnetized drain segment, as a function of the magnetization direction and the electron energy.

	Besides flat slabs and flakes, another common morphology for experimentally grown TI nanostructures is the nanowire \cite{SciRep3_1212, SciRep3_1564, APL103_193107}, which can approximately be described as a cylinder. The effective Hamiltonian for the surface states of a TI cylinder \textit{without} the hexagonal warping term has been derived to be \cite{PRB79_245331,PRL105_206601,PRL105_136403}
\begin{equation}
		H = v(p_\phi\sigma_z - p_z\sigma_\phi).  \label{Hcy0}
\end{equation}

	The surface states of a TI cylinder differ from those of the flat TI slab in two important aspects. First, the circumference of the cylinder imposes periodic boundary conditions on the surface states along the azimuthal direction. Second, the curvature along the circumference of the cylinder leads to a position dependent normal vector to the surface. This in turn leads to a position dependent spin orbit interaction field, which is reflected by the presence of the $\sigma_\phi\equiv (\cos(\phi)\sigma_x+\sin(\phi)\sigma_y)$ term in the Hamiltonian.

	The above observations motivate us to analyze the combined effects of the  of the hexagonal warping term together with a cylindrical geometry in \ce{Bi2Te3} nanowires. Moreover, whilst the effective Hamiltonian for the surface states of a cylindrical TI nanowire has already been derived,\cite{PRB84_195406} the hexagonal warping terms have been neglected in the derivation. Thus in this work, we derive the effective Hamiltonian for the surface states of a \ce{Bi2Te3} nanowire which now includes the hexagonal warping term starting from the four-band Hamiltonian of Liu \textit{et al}. To test the accuracy of the effective surface state Hamiltonian, we compare the bandstructure from the Hamiltonian with that from the full four-band Hamiltonian, and show reasonable matching between the two bandstructures.

	As an example applicatin of the surface state Hamiltonian, we study the transmission between two \ce{Bi2Te3} magnetized in different directions perpendicular to the cylinder axis.  The hexagonal warping term leads to a dependence of the energy-$k_z$ relation on the magnetization direction in the plane of the cylinder circumference. This in turn results in a sawtooth transmission profile in the transmission between two \ce{Bi2Te3} cylinders magnetized in different directions, a feature which may be of utility in achieving multiple-state magnetic memory bits.

\section{Derivation of Effective Hamiltonian}
We take as our starting point the model four-band Hamiltonian in Eqs. 16 and 17 of Ref. [\onlinecite{PRB82_045122}],
\begin{eqnarray}
	H_{(4B)} &=& (M_0 + M_1 k_z^2 + M_2 k_\parallel^2)\Gamma_5 + B\Gamma_4k_z   \nonumber \\
	&&  + A(\Gamma_1k_y-\Gamma_2k_x) + R_1\Gamma_3(k_x^3-3k_xk_y^2) \nonumber \\
	&& + R_2\Gamma_4(3k_x^2k_y-k_y^3) \label{HCa}
\end{eqnarray}
where
\begin{eqnarray*}
	&& \Gamma_1 = \sigma_x\otimes\tau_1,	\Gamma_2 = \sigma_y\otimes\tau_1, \\
	&& \Gamma_3 = \sigma_z\otimes\tau_1, \Gamma_4 = \mathbb{I}_\sigma\otimes\tau_2, 	\Gamma_5 = \mathbb{I}_\sigma\otimes\tau_3,
\end{eqnarray*}
and $k_\parallel^2 \equiv k_x^2+k_y^2$. The $\tau$'s can be interpreted as describing the orbital degree of freedom and the $\sigma$ terms the real spins.  This Hamiltonian can be derived based on the symmetries exhibited by the underlying crystal structure of the \ce{Bi2Se3} family of TIs, and describes both the bulk as well as surface states.  (The subscript (4B) stands for ``4-band''. )

We first derive the two-band effective Hamiltonian for the surface states of a TI (\ce{Bi2Te3}) cylinder with its axis parallel to the $z$ direction without the $R_1$ and $R_2$ terms (i.e. without the hexagonal warping), and incorporate the effect of curvature by an approach modified from Imura \textit{et al} \cite{PRB84_195406}. This would yield a Hamiltonian compatible to the usual Dirac Hamiltonian of Eq. \eqref{Hcy0}.  The general idea is to first split the four-band Hamiltonian $H_{(4B)}$, which describes both the bulk and surface states, into two parts which we call the perpendicular Hamiltonian $H_{(4B),\perp}$ and parallel Hamiltonian $H_{(4B),\parallel}$ for short. The first part  contains momentum operators perpendicular to the surface, and the second part momentum operators parallel to the surface.

We know that the surface states exist near the zero-energy Dirac point where the energy contribution due to the perpendicular Hamiltonian is minimal. In other words, the energy of these states should come only from the momentum of the particles moving on the surface. We hence find the zero energy ground states of the perpendicular Hamiltonian which are localized near the surface, and then treat the parallel Hamiltonian as a perturbation to the perpendicular Hamiltonian ground states. The ground states are degenerate due to the presence of the orbital and spin degrees of freedom. According to degenerate perturbation theory, the first order effects of a perturbation are given by projecting the perturbing Hamiltonian in the basis of the degenerate ground states. The parallel Hamiltonian projected into the basis of the degenerate ground states thus yields the effective Hamiltonian for the states localized near the surface.

 In detail, we first identify $k_y\rightarrow(-i\partial_y)$, $(-i\partial_y) = - i(\frac{\partial r}{\partial_y}\partial_r + \frac{\partial \phi}{\partial_y}\partial_\phi)$, etc., to rewrite Eq. \eqref{HCa} in cylindrical coordinates $(r, \phi)$. We treat the terms containing $r$ derivatives resulting from the $(M_0 - M_2\partial_r^2)\Gamma_5 + A(\Gamma_1 k_y-\Gamma_2 k_x)$ as our base Hamiltonian and the remaining terms as perturbations.

\subsection{Zero energy eigenstates of $H_{(4B),\perp}$}
The perpendicular Hamiltonian in cylindrical coordinates reads as
\begin{eqnarray}
	\langle\vec{r}|H_{(4B),\perp} &=& (M_0 - M_2\partial_r^2)(\mathbb{I}_\sigma\otimes\tau_3) + A ( - (\sigma_\phi)(-i\partial_r) )\otimes\tau_1 \nonumber  \\
	&\equiv& M_\perp(\mathbb{I}_\sigma\otimes\tau_3) + A ( - (\sigma_\phi)(-i\partial_r) )\otimes\tau_1
	\label{H0Cy}
\end{eqnarray}
where, following Ref. [\onlinecite{PRB84_195406}], we dropped the term resulting from $M_2(\frac{1}{r})(-\partial_r)$. This is justifiable on the basis that the edge states we are seeking are confined to around $r=R$ where in the $R\rightarrow\infty$ limit the contributions of terms proportional to $\frac{1}{R}$ go to 0. As mentioned above, we look for eigenstates of  Eq. \eqref{H0Cy} with eigenenergy 0 and with spatial $r$ dependence of the form $\exp(\lambda(r-R))$.

It turns out to be physically more instructive to diagonalize the $\sigma$ part of Eq. \eqref{H0Cy} first. The $\sigma$ eigenstates $|\sigma_\pm\rangle$ can be  solved to be $|\sigma_\pm\rangle \equiv \frac{1}{\sqrt{2}} (1, \pm i \exp(i \phi))^\mathrm{T}$. Substituting the eigenstates into Eq. \eqref{H0Cy}, we then have an equation for the $\tau$ degree of freedom --
\begin{equation}
	\langle\vec{r}, \sigma_\pm|H_{(4B), \perp} = (M_\perp \tau_3  \pm i A\lambda \tau_1)\langle\vec{r},\sigma_\pm|.
	\label{H1Cy}
\end{equation}

We now solve for the zero energy eigenstates of Eq. \eqref{H1Cy}.  Writing the (unnormalized) zero energy eigenstate of Eq. \eqref{H1Cy} as $(1, \alpha)^\mathrm{T}$ and solving for $\alpha$ yields the eigenstates $(1, \pm \frac{i M_\perp}{A\lambda})^\mathrm{T}$. Setting the secular equation for the eigenvalues of Eq. \eqref{H1Cy} to be 0 then yields $(M_\perp)^2 = (A\lambda)^2$.

By definition $M_\perp = M_0 - \lambda^2 M_2$. A condition for the existence of topological edge states is that $M_0$ and $M_2$ need to be of differing signs. It is conventional to choose $M_0$ to be negative, in which case $M_\perp$ is negative. As $M_\perp$ is negative and $A$ is positive, the only consistent solution for  $(M_\perp)^2 = (A\lambda)^2$ is to choose $M_\perp=-A\lambda$. The (normalized) zero-energy eigenstates of Eq. \eqref{H1Cy} are then
\[
		|\tau_\pm\rangle = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 \\ \mp i \end{pmatrix}.
\]

We note that the equation $M_\perp=-A\lambda$ actually expands out to $(M_0-\lambda^2M_2)=-A\lambda$ which is quadratic in $\lambda$ and hence yields two solutions which we denote as $\lambda_{(\pm)}$. Note that the bracketed $(\pm)$ in this subscript differs from the $\pm$ index which pertains to the spin ($\sigma$) eigenstate of Eq. \eqref{H0Cy}.
 Explicitly, we have
\begin{equation}
	\lambda_{(\pm)} = \frac{A\pm\sqrt{A^2+4M_0M_2}}{2M_2}. \label{lambdaPM}
\end{equation}

Putting everything together, a zero-energy eigenstate of Eq. \eqref{H1Cy} which vanishes at the cylinder boundaries $r=R$ is given by
\begin{eqnarray*}
		&& \langle \vec{r}|\pm\rangle \\
		&\propto&  \frac{\exp(\lambda_{(+)} (r-R)) - \exp(\lambda_{(-)}(r-R))}{2} \begin{pmatrix} 1 \\ \mp i \\ \pm i\exp(i\phi) \\ \exp(i\phi) \end{pmatrix}.
\end{eqnarray*}

To obtain the normalization constant, we evaluate the following integral in the limit of  $\lambda_{(\pm)}R \gg 1$,
\begin{eqnarray*}
	&& \int^R_0 \ \mathrm{d}r\  2\pi r (\exp(\lambda_{(+)} (r-R))-\exp(\lambda_{(-)} (r-R)))^2\\
	&\approx& \frac{\pi R (\lambda_{(-)}-\lambda_{(+)})^2 }{\lambda_{(-)}\lambda_{(+)}(\lambda_{(+)}+\lambda_{(-)})}.
\end{eqnarray*}
(Ref. [\onlinecite{PRB84_195406}] writes of the $\lambda_{(\pm)} \ll R$ limit but we find that $\lambda_{(\pm)}R \gg 1$ is a better description of what we end up operationally doing.) The normalized zero-energy eigenstates thus read as
\begin{equation}
	\langle \vec{r}|\pm\rangle = \rho(r)\begin{pmatrix} 1 \\ \mp i \\ \pm i\exp(i\phi) \\ \exp(i\phi) \end{pmatrix}
\end{equation}
where $\rho(r)  = \sqrt{\frac{\lambda_{(-)}\lambda_{(+)}(\lambda_{(+)}+\lambda_{(-)})}{\pi R (\lambda_{(-)}-\lambda_{(+)})^2 }}[\exp(\lambda_{(+)} (r-R))-\exp(\lambda_{(-)} (r-R))]$.

\subsection{Projection of $H_{(4B),\parallel}$ into zero-energy $(2\times2)$ subspace of $H_{(4B),\perp}$}
We now proceed to derive the effective Hamiltonian on the surface state of the \ce{Bi2Te3} cylinder. We first consider $H_{(4B),0}=(M_0 + M_1 k_z^2 + M_2 k_\parallel^2)\Gamma_5+B\Gamma_4 k_z + A(\Gamma_1 k_y - \Gamma_2 k_x)$ [where the ``0"' in the subscript denotes the absence of the $R_1$ and $R_2$ hexagonal warping terms in Eq. \eqref{HCa}]. From this, we extract the Hamiltonian $H_{(4B),\parallel; 0}$ that corresponds to motion parallel to the surface of the cylinder [i.e., the terms which are not already included in $H_{(4B), \perp}$ in Eq. \eqref{H0Cy}]:
\[
	\langle\vec{r}|H_{(4B),\parallel; 0} = (M_1 k_z^2)\Gamma_5+B\Gamma_4 k_z - i\frac{A}{r}\sigma_r\otimes\tau_1\partial_\phi.
\]

We then project $H_{(4B),\parallel; 0}$ into the subspace spanned by $|\pm\rangle$, the zero-eigenstates of $H_{(4B),\perp}$. In other words, we treat $H_{(4B),\parallel;0}$ as a perturbation to $H_{(4B),\perp}$. An explicit calculation of $\int \mathrm{d}r\ r \langle\vec{r},\pm'|H_{(4B),\parallel; 0}|\vec{r},\pm\rangle$ for the four possible combinations of $\pm'$ and $\pm$ yields
\begin{eqnarray*}
	\tilde{H}_{(2B);0}  &\equiv&   \begin{pmatrix} \langle  +|H_{(4B),\parallel;0}|+ \rangle &  \langle  +|H_{(4B),\parallel;0}|- \rangle \\  \langle -|H_{(4B),\parallel;0}|+\rangle &  \langle  -|H_{(4B),\parallel;0}|-\rangle \end{pmatrix} \\
	&=& \begin{pmatrix} B k_z & A p_\phi  \\
	 A p_\phi & - B k_z \end{pmatrix}.
\end{eqnarray*}

The spin components of the matrix above are in the $| \sigma_\pm\rangle$ basis. They can be converted to the original spin basis of Eq. \eqref{HCa} through the unitary transformation $H_{(2B);0} = U\tilde{H}_{(2B);0}U^\dagger$ with
\[
U = \sum_{\sigma=(\uparrow/\downarrow),\tilde{\sigma}=\sigma_{\pm}} |\sigma\rangle\langle\sigma|\tilde{\sigma}\rangle\langle\tilde{\sigma}|.
\]
This yields
\begin{equation}
	H_{(2B);0} = A \left(p_\phi\sigma_z+\frac{1}{2R}\mathbb{I}_\sigma\right) - B k_z\sigma_\phi \label{H1ti},
\end{equation}
where $p_\phi\equiv(-i\partial_\phi)/R$. Eq. \eqref{H1ti} coincides with the usual Rashba form of the Dirac fermion Hamiltonian $v (\vec{\sigma}\times\vec{p})\cdot\hat{n}$ on the surface of the cylinder expanded out in cylindrical coordinates with the identification of $A=B/\hbar=v$ except for the presence of the additional $\frac{A}{2R}$ term. The latter has been interpreted as a manifestation of the spin Berry phase \cite{PRB84_195406} but can be ignored for most transport calculations as a constant shift in energy.

Finally, we consider the remaining terms in Eq. \eqref{HCa} that have not yet been dealt with, and which correspond to the warping terms, i.e.
\[
	 H_{(4B),\parallel;R}=R_1\Gamma_3 (k_x^3-3k_xk_y^2) + R_2\Gamma_4(3k_x^2k_y-k_y^3).
\]

These terms require a bit more subtlety to handle as they contain $r$ and its derivatives. We often see the momentum in the $x$ direction being written as $\hat{p}_x\simeq-i\partial_x$, so that the expectation value of the momentum at position $x$ gets written as $\psi(x)^*(-i\partial_x \psi(x))$. We should, however, recognize that $\psi(x)^*(-i\partial_x \psi(x))$ is actually $\langle \psi| \Big( |x\rangle\langle x|p_x\Big)|\psi\rangle$ where we are taking the expectation value of the \textit{non-Hermitian} operator $|x\rangle\langle x|p_x$ with respect to the state $|\psi\rangle$. This does not lead to any problems when $\psi(x)=\langle x|\psi\rangle = \exp(i k x)$ where $k$ is real, but yields a complex value where $k$ is complex. This contravenes the fact that an expectation value which corresponds to a physical measurement should yield a real value. The most natural way around this is to modify the expression for the momentum to $\mathrm{Re} (\psi^*(-i\partial_x \psi)) = \frac{1}{2} ( (i\partial_x\psi^*)\psi + \psi^*(-i\partial_x\psi)) = \frac{1}{2}\langle \psi| \Big( \{ |x\rangle\langle x|, p_x\} \Big) |\psi\rangle$. The preceding discussion is relevant for the case of a TI cylinder because the $r$ dependence of the $|\pm\rangle$ states take the form of $\exp(- r\lambda)$ which can be interpreted as a wave with a complex wavevector. Following our discussion,  the matrix elements of the projected $(2\times2)$ Hamiltonian $H_{(2B);R}$ which represents the warping terms acting on the surface states, are given by
\begin{eqnarray}
	 \left[H_{(2B);R}\right]_{\pm',\pm} &=& \frac{1}{2} \langle \pm'|\{ R_1 \Gamma_3 (p_x^3-3p_xp_y^2)\nonumber \\
	&& + R_2 \Gamma_4 (3p_x^2p_y-p_y^3), |\vec{r}\rangle\langle\vec{r}|\}|\pm\rangle.\nonumber
\end{eqnarray}
Evaluating the integrals, dropping terms containing $\exp(-\lambda R)$, and performing the unitary transformation to the real spin basis as before, gives the final expression
\begin{eqnarray}
	H_{(2B)} &\equiv& H_{(2B;0)} +H_{(2B;R)} \nonumber\\
	&=&-B k_z \sigma_\phi  + \frac{A}{R}\left(\sigma_z(-i\partial_\phi)+\mathbb{I}_\sigma\frac{1}{2}\right) \nonumber \\
	&-&\frac{3i\lambda_{(-)}\lambda_{(+)}}{2R} \{ \partial_\phi,(R_1 \sin(3\phi)\sigma_r - R_2 \cos(3\phi)\sigma_\phi)\}.\nonumber\\  \label{Heff}
\end{eqnarray}
This effective Hamiltonian is the key result of this work. The steps outlined in the previous section can generally be applied to reduce a full 4-band Hamiltonian to a $(2\times2)$ surface state Hamiltonian. For instance, the surface state Hamiltonian for a flat TI slab with warping term in Eq. \eqref{ham}, as obtained by Fu, can be obtained from the full Hamiltonian in Eq. \eqref{HCa} by following the same procedure, but with the TI cylinder replaced by a TI crystal truncated perpendicular to the $z$ direction, as we briefly outline in the Appendix.

We then find that the parameters $v$ and $\lambda$ in Eq. \eqref{ham} correspond to $A$ and $R_1$ in Eq. \eqref{HCa}, respectively. The values of $A$ in both Refs. [\onlinecite{PRL103_266801}] and [\onlinecite{PRB82_045122}] match within 10\% ($2.55\ \mathrm{eV}$\AA  versus $2.87\ \mathrm{eV}$\AA, respectively)	 and we shall use the latter value for our numerical calculations. The values of $R_1$, however, differ by a factor of 5.5 ($250~\mathrm{eV}$\AA$^3$ versus $45.02~\mathrm{eV}$\AA$^3$) between the two works. We assume the value of  $R_1$ to be 5 times the latter value in our numerical calculations. This value is chosen as it is intermediate between the values quoted in the two works, but closer to towards that in Ref. [\onlinecite{PRL103_266801}] which was obtained by comparison to actual experimental data, and hence may be deemed to be more reliable.

\section{Comparison between two-band and four-band Hamiltonians}
In this section, we compare how well the effective two-band Hamiltonian of Eq. \eqref{Heff} for the surface states replicates the surface states obtained via the full four-band Hamiltonian of Eq. \eqref{HCa}. To accomplish this, we write the spatial parts of the latter in the basis of the normalized eigensolutions to the two-dimensional cylindrical coordinate Laplacian equation which vanish at $r=R$, thereby satisfying the boundary condition that the eigenstate wavefunctions go to 0 at the surface of the cylinder.  These eigensolutions $|m,l\rangle$ are indexed by the quantum numbers $m$ and $l$ which carry the usual quantum mechanical interpretations of the $z$ and total angular momentum respectively, and are given by
\[
	\langle \vec{r}|m, l\rangle = \frac{\exp(i m\phi) J_m\left(\frac{r}{R} j_{m,l}\right)j_{m,l}^2}{\sqrt{\pi}R^3J_{m+1}(j_{m,l})}
\]
where $J_{m}$ is the Bessel function of the first kind of order $m$, and $j_{m,l}$ its $l$th zero. (A theorem in elementary quantum mechanics states that the eigenstates for a Hermitian Hamiltonian form a complete basis set, so we do not need to consider states proportional to $\exp(i m\phi)J_n(\frac r{R} j_{n,l})$ for $m\neq n$.)  The resulting matrix is then diagonalized numerically to yield the energy eigenvalues.

There is a technical caveat that has to be considered in evaluating the matrix elements involving the third order derivatives in $r$, which occur in the $R_1$ and $R_2$ hexagonal warping terms. For example, the terms in the Hamiltonian containing $R_1$ terms acting on arbitrary states $\langle \Psi|$ and $|\Phi \rangle$ can be most succinctly be expressed by
\begin{eqnarray*}
	&& R_1 \Psi^\dagger \Gamma_3(k_x^3-3k_xk_y^2) \Phi \\
	=&&R_1 \mathrm{Im} \Psi^\dagger\Gamma_3\Big[\cos3\phi\Big(-\frac{3 \partial_\phi^2}{r^3}-\frac{3 \partial_r}{2 r^2}+\frac{3 \partial_r\partial_\phi^2}{2 r^2}+\frac{3 \partial_r^2}{2 r}\\
	&-&\frac{1}{2} \partial_r^3 \Big)+\sin3\phi\Big(\frac{4 \partial_\phi}{r^3}-\frac{\partial_\phi^3}{2 r^3}-\frac{9 \partial_r\partial_\phi}{2 r^2}+\frac{3 \partial_r^2\partial_\phi}{2 r}\Big)\Big]\Phi.
\end{eqnarray*}
One might then naively evaluate the $\partial^3_r$ matrix elements $\langle m',l'\partial^3_r |m,l\rangle$ as
\begin{equation}
	\int^R_0 \mathrm{d}r\  r J_{m'}(r/R j_{m',l'})\partial^3_r J_m(r/R j_{m,l}). \label{r30}
\end{equation}
It turns out, however, that such a na\"{i}ve calculation would not yield a Hermitian Hamiltonian. A more proper interpretation of the matrix element should be
\begin{equation}
\int^\infty_0 \mathrm{d}r\ r J_{m'}(r/R j_{m',l'})\partial^3_r (\Theta(R-r)J_m(r/R j_{m,l})),\label{r31}
\end{equation}
where $\Theta$ is the Heaviside step function. Compared to Eq. \eqref{r30}, Eq. \eqref{r31} has extra terms after integration due to the Dirac delta functions that arise from differentiating the step function. \footnote{These extra terms are also present in the evaluation of integrals involving the second $r$ derivatives but evaluate to 0 as they involve a product of a $J_m(r/R j_{m,l})$ and the derivative of another one at $r=R$.} For matrix elements involving third $r$ derivatives, the extra terms include products of derivatives of Bessel functions which do not evaluate to 0 at $r=R$, and have to be accounted for.

	The calculation of the dispersion relations with the full four-band Hamiltonian enables us to compare how well our effective two-band Hamiltonian approximates the surface states of the four-band Hamiltonian. The calculation also establishes the energy and wavevector ranges above which bulk states start to emerge, and would invalidate any transport calculations based on the two-band Hamiltonian, which considers only the surface states.

	Fig. \ref{gcombEK1} shows the dispersion relations calculated for \ce{Bi2Se3} and \ce{Bi2Te3}  using the four-band Hamiltonian parameters given in Ref. [\onlinecite{PRB82_045122}] for TI nanocylinders except for the $R_1$ parameter for $\ce{Bi2Te3}$, for which we assume a value 5 times that quoted in the paper as discussed in the previous section.

\begin{figure}[ht!]
\centering
\includegraphics[scale=0.5]{T2combEK1.eps}
\caption{  Dispersion relations for the edge states and some bulk states for  (a), (b) \ce{Bi2Se3} and (c), (d) \ce{Bi2Te3} for TI nanocylinders of radii  (a),(c) $15~\mathrm{nm}$  and (b),(d) $30~\mathrm{nm}$. }
\label{gcombEK1}
\end{figure}

	As the radius of the cylinders increases, the energy spacings between the bulk and edge bands decrease. In \ce{Bi2Se3} the edge states are more cleanly separated in energy from the bulk states compared to \ce{Bi2Te3} where some of the higher energy edge states are found at energy ranges where bulk states are present as well. Somewhat surprisingly, the energy values (and for \ce{Bi2Te3}, the wavevectors) at which the bulk states emerge and render the two-band effective Hamiltonian for surface states inappropriate, is only very weakly dependent on the cylinder radii. The results indicate that the bulk states may be ignored up to around $0.25\ \mathrm{eV}$ in \ce{Bi2Se3} and $0.1\ \mathrm{eV}$ in \ce{Bi2Te3} cylinders, for the range of cylinder radii considered.

\begin{figure}[ht!]
\centering
\includegraphics[scale=0.55]{T2efEkCompB.eps}
\caption{  The dispersion relations for a TI cylinder of radius (a) $15\ \mathrm{nm}$ and (b) $30\ \mathrm{nm}$ calculated using the full four-band Hamiltonian (blue lines), the two-band effective Hamiltonian (red lines) and the two-band Hamiltonian without the hexagonal warping terms (green dots) .}
\label{gefEkComp}
\end{figure}

Next, we compare the dispersion relations for \ce{Bi2Te3} obtained from the two-band and four-band Hamiltonians. Fig. \ref{gefEkComp} shows the dispersion relations for two TI cylinders of different radii calculated using the full and the effective Hamiltonian, as well as the effective Hamiltonian without the hexagonal warping term $H_{(2B);0} = A \left(p_\phi\sigma_z+\frac{1}{2R}\mathbb{I}_\sigma\right) - B k_z\sigma_\phi$. The depicted edge states are two-fold degenerate, i.e., each distinct line in the figure represents two states.  In both panels of the figure, there is evidently a significant mismatch between the dispersion relations of the two-band and four-band Hamiltonians if the hexagonal warping terms of the two-band Hamiltonians were neglected (green dots). On the other hand, the dispersion relations obtained by the effective surface state Hamiltonian of Eq. \eqref{Heff} show a reasonably close fit with those obtained by the full Hamiltonian of Eq. \eqref{HCa} at small $k_z$ and at energy values lower than the threshold at which the bulk bands emerge. At larger values of $k_z$ below $0.1\ \mathrm{eV}$, the effective Hamiltonian tends to underestimate the eigenenergy. However the qualitative trend that the bottom-most pairs of degenerate particle states tend to cluster together at large $k_z$ is reproduced by the effective Hamiltonian. The agreement between the effective and full Hamiltonian $E-k$ curves improves with increasing cylinder radius. This is to be expected as the effective Hamiltonian was derived under the $R\rightarrow\infty$ limit.

\section{Dispersion relations and magnetoresistance}
Having now established that our effective two-band Hamiltonian is an adequate approximation for the full four-band Hamiltonian for a \ce{Bi2Te3} nanowire, we now use the former to study a \ce{Bi2Te3} cylinder with an magnetization. We use the numerical value of $50\ \mathrm{nm}$ for the cylinder radius in our calculations. This choice is motivated by two reasons. The results of the previous section suggests a reasonably good match between  the eigenstates calculated by the four-band Hamiltonian and the two-band effective Hamiltonian can be obtained at $50\ \mathrm{nm}$ radius.  This radius value also falls within the range reported for experimentally fabricated \ce{Bi2Te3} nanowires \cite{SciRep3_1212,  SciRep3_1564}.

\begin{figure}[ht!]
\centering
\includegraphics[scale=0.67]{T2muXcombEKb.eps}
\caption{   The dispersion relations for an magnetization applied in the $x$ direction for a \ce{Bi2Te3} nanowire of radius $500$ \AA, and (a) $M_x = 0.01\ \mathrm{eV}\hbar^{-1}$, (b) $M_x = 0.02\ \mathrm{eV}\hbar^{-1}$ and (c) $M_x = 0.04\ \mathrm{eV}\hbar^{-1}$ respectively plotted in solid lines, and the dispersion relations without the magnetizations plotted in dotted lines for comparison. }
\label{gT2muXcombEK}
\end{figure}

The panels of Fig. \ref{gT2muXcombEK} shows the dispersion relations for increasing magnetizations in the $x$ direction. Similar to the TI nanowires studied in Ref. [\onlinecite{JAP117_17C749}],  the energies of the particle states at $k_z=0$ increase with increasing magnetizations, and the low energy particle bands become convex at small $|k_z|$ at large magnetizations.

\begin{figure}[ht!]
\centering
\includegraphics[scale=0.48]{T2combMuX2b.eps}
\caption{   (a) The dispersion relations for a \ce{Bi2Te3} nanowire of radius $500$ \AA with magnetization  $\vec{M} = M (\cos(\phi_M), \sin(\phi_M), 0)$ for $\phi_M = 0$ (red lines), $\frac{\pi}{12}$ (green lines) and $\frac{\pi}{6}$ (blue lines). (b) The energy of the first band whose band bottom is visible in panel (a) plotted as a function of $k_z$ and $\phi_M$.  }
\label{gT2combMuX2}
\end{figure}

Differing from the TI nanocylinders studied in Ref. [\onlinecite{JAP117_17C749}],  the hexagonal warping terms in \ce{Bi2Te3} result in a dependence of the dispersion relations on the direction of the applied magnetization perpendicular to the cylinder axis. Panel (a) of Fig. \ref{gT2combMuX2} shows that for a given band, the value of $\phi_M$ which gives the lowest energy depends on both the band as well as the value of $k_z$. Panel (b) in turn shows that the energy of each band varies periodically with a period of $\pi / 3$ as $\phi_M$ is varied. This is reminiscent of the behavior of Fig. 2a of Ref. [\onlinecite{SciRep4_5062}]  in which the bandgap in the flat $+z$ terminated \ce{Bi2Te3} slab subjected to an in-plane magnetization exhibits a similar 6 fold periodic variation with the in-plane magnetization angle.

The combined effects of the band bending at small $|k_z|$ and the magnetization angle dependence give rise to a rich transmission profile as the energy and magnetization angles are varied.  Fig. \ref{gcombMrX} shows the transmission from a  source cylinder magnetized in the $+x$ direction to a coaxial drain cylinder  magnetized to the same magnitude shown in Fig. \ref{gT2mr}, as a function of the drain magnetization angle and the energy for two values of source magnetization.

\begin{figure}[ht!]
\centering
\includegraphics[scale=0.48]{T2mr}
\caption{   The system considered in this section consisting of a source \ce{Bi2Te3} cylinder magnetized in the $+x$ direction and a coaxial drain \ce{Bi2Te3} cylinder, both of radius $50 \mathrm{nm}$, coaxial with each other.    }
\label{gT2mr}
\end{figure}

\begin{figure}[ht!]
\centering
\includegraphics[scale=0.33]{T2combMrX}
\caption{   The left panels show transmission from a $50\ \mathrm{nm}$ TI cylinder magnetized in the $+x$ direction to a co-axial drain cylinder magnetized to the same magnitude and in the $(\cos(\phi_M),\sin(\phi_M),0)$ direction at (a) $M=0.02\ \mathrm{eV}\hbar^{-1}$ and (b) $M=0.03\ \mathrm{eV}\hbar^{-1}$ as a function of energy and drain magnetization angle $\phi_M$. The panels on the right show the dispersion relations at various magnetization angles at the same magnetization magnitude as the accompanying plots on the left, as a function of the wave vector $k_z$ and the energy. The energy values of the transmission profiles (left) and dispersion relations (right) are aligned to facilitate comparison. We have also drawn lines linking the energy at which key features in the transmission profile emerge with their corresponding features in the band diagrams.    }
\label{gcombMrX}
\end{figure}

We first consider panel (a) of Fig. \ref{gcombMrX} which corresponds to the smaller magnetization magnitude.  At energies below slightly more than $0.015\ \mathrm{eV}$, theew is only one propagating mode in both the source and drain. The transport in this regime exhibits the conventional magnetoresistance where the transmission is highest for parallel magnetizations between the source and drain, and lowest for the antiparallel configuration.  At energies lying between slightly more than $0.015\ \mathrm{eV}$ to about $0.025\ \mathrm{eV}$ the transmission exhibits the anomalous magnetoresistance in Ref. [\onlinecite{JAP117_17C749}]. The band-bending near  $|k_z=0|$ results in the transmission at antiparallel configuration being higher than at intermediate configurations between parallel and anti-parallel source and drain magnetization configurations. Brieftly, the anomalously high transmission at antiparallel transmission occurs because for a given direction of current flow, the states whose bands bend for a given magnetization direction, are localized at the same angular positions around the cylinder circumference as those states whose bands do not bend for the opposite magnetization direction. The overlap between the angular postions of where the wavefunctions have the largest ampltidue of the bent bands on one side, and the unbent bands on the other side, of the source-drain interface lead to increased transmission from the source to the drain and hence the anomalously high transmission \cite{JAP117_17C749}.  In contrast, at higher energies above approximately $0.023\ \mathrm{eV}$ in panel (a) where band bending does not occur we have again the conventional magnetoresistance.

A `sawtooth' profile of high transmission at $\phi_M = n \pi/3$  appears in a narrow energy range between about $0.019\ \mathrm{eV}$ to slightly less than $0.021\ \mathrm{eV}$. The emergence of this sawtooth pattern is due to the variation of the energy of the band bottom of the band emerging at this energy range with the drain magnetization direction -- the narrow tip of each  high transmission `tooth' at each integer $n$ $\phi_M = n\pi/3$ corresponds to the drain magnetization direction at which the band bottom has the lowest value of energy. At higher values of energy the band bottoms of the bands which subsequently emerge do not vary as much with the drain magnetization direction, and the bands no longer bend upwards at small $|k_z|$. This leads to a conventional magnetoresistance behavior at high energies. Comparing now between panels (b) and (a), we see that the larger magnetization  in (b) results in more extensive bending of the lower energy bands at small $|k_z|$. This results in a larger energy range over which the anomalous magnetoresistance appears. The larger variation of the energy of each band with the magnetization angle also leads to a larger energy range for the sawtooth transmission profile.

The sawtooth transmission profile may possibly be of utility as a four-state magnetic memory element. The magnetic memory element may be biased and gated so that the total current is given by the energy integral over an energy range falling within the sawtooth transmission profile regime (about $0.024\ \mathrm{eV}$ to $0.027\ \mathrm{eV}$ ), and the 4 states represented by the drain magnetization lying in $\phi_M = 0$, $\pi/6$, $2\pi/6$ and $3\pi/6$. The sawtooth transmission profile then leads to a relatively large difference in current between memory states represented by  adjacent values of the drain magnetization directions  (e.g. between $\phi_M=\pi/6$ and $\phi_M=2\pi/6$) while retaining a current difference large enough between the $\phi_M=0$ and $\phi_M=2\pi/6$, and $\phi_M =\pi/6$ and $\phi_M=3\pi/6$ pairs, to be able to robustly distinguish between the four states.

\section{Conclusion}
In this work we derived the two-band effective Hamiltonian for the surface states of a \ce{Bi2Te3} nanocylinder incorporating the hexagonal warping terms. We calculated the dispersion relations from the underlying four-band Hamiltonian and saw that the energies at which the bulk states emerge is only weakly dependent on the cylinder diameter. The dispersion relations calculated from the two-band Hamiltonian match those of the four-band Hamiltonian reasonably well. The dependence of the \ce{Bi2Te3} dispersion relation on the angle at which the magnetization is applied perpendicular to the cylinder axis gives a rich transmission profile between two cylinders magnetized in different directions. In particular, the sawtooth profile which exists at some values of energy and large magnetizations may be of utility in a four-state magnetic memory bit.

\section*{Acknowledgment}
The authors acknowledge the Singapore National Research Foundation for support under NRF Award Nos. NRF-CRP9-2011-01 and NRF-CRP12-2013-01, and MOE under Grant No. R263000B10112.

\section*{Appendix}
Here, we briefly outline the derivation of the effective Hamiltonian for the surface states of a flat \ce{Bi2Te3} slab of infinite dimensions in the $x$ and $y$ directions, semi-infinite thickness in the $z$ direction, and terminated by a flat surface with the outward normal in the $+z$ direction. For convenience, we set the surface to be at $z=0$, so that the surface states which decay exponentially into the surface as $z\rightarrow-\infty$ have the form $\exp(\lambda^{(z)} z)$.

In this case, the perpendicular Hamiltonian containing derivatives in the $z$ direction reads, in Cartesian coordinates, as
\begin{equation}
	\langle \vec{r}|H^{(z)}_{(4B),\perp} = M^{(z)}_\perp(\mathbb{I}_\sigma\otimes\tau_3) + B (\mathbb{I}_\sigma\otimes\tau_2) k_z \label{HzPt},
\end{equation}
where $M^{(z)}_\perp \equiv (M_0 - M_2\partial_z^2)$, and the $(z)$ superscripts indicate that these quantities pertain to the $+z$ terminated slab. We solve for the zero-energy eigenstates of $H^{(z)}_{(4B),\perp}$. These zero-energy eigenstates which are proportional to $\exp(\lambda^{(z)} z)$ and hence, consistent with our boundary condition, exist only when $M^{(z)}_\perp = -\lambda^{(z)} /B$ (similar to the case of the cylindrical geometry considered in the main paper). Solving $M^{(z)}_\perp = -\lambda^{(z)} /B$ for $\lambda^{(z)}$ in turn gives two solutions, i.e., $\lambda^{(z)}_{\pm} = (B \pm \sqrt{B^2+4M_0M_1})/(2M_1)$. The zero-energy eigenstates of Eq. \eqref{HzPt} can then be chosen to be proportional to
\begin{eqnarray}
	\langle \vec{r} |+^{(z)}\rangle &\propto&(\exp(\lambda^{(z)}_+z) - \exp(\lambda^{(z)}_-z) )   \begin{pmatrix} 1 \\ -1 \\ 0 \\ 0 \end{pmatrix}, \\
	\langle \vec{r}  |-^{(z)}\rangle &\propto& (\exp(\lambda^{(z)}_+z) - \exp(\lambda^{(z)}_-z) ) \begin{pmatrix} 0 \\ 0 \\ 1 \\ -1 \end{pmatrix}.
\end{eqnarray}

Taking the projection of the remaining terms of Eq. \eqref{ham} not already included in  Eq. \eqref{HzPt}, onto the $|\pm^{(z)}\rangle$ states then yields Fu's surface Hamiltonian of Ref. [\onlinecite{PRL103_266801}]
\begin{equation}
H = A (k_x\sigma_y - k_y\sigma_x) + R_1 (k_x^3 - 3 k_x k_y^2) \sigma_z,
\end{equation}
upon identifying $A \rightarrow v_f$ and $R_1 \rightarrow \lambda $.


\section{Conclusion}
We present a simple non-parametric model for clustering short documents (such as tweets) into storylines, which are conceptually coherent and temporally focused.
Future work may consider learning more flexible temporal distance functions, which could potentially represent temporal periodicity or parametric models of content popularity.
\section{TREC Evaluation}
To test the efficacy of this approach, we evaluate on the Twitter Timeline Generation (TTG) task in the Microblog track of TREC 2014. It involves taking tweets based on a query $Q$ at time $T$ and returning a summary that captures relevant information. We perform the task on 55 queries with different timestamps and compare our results with 13 groups that submitted 50 runs for this task in 2014.

We consider the following systems:
\begin{description}
\setlength\parskip{0pt}
\setlength\parsep{0pt}
\item[Baseline] We replace the distance-dependent prior with a standard Dirichlet prior. The number of clusters is heuristically set to 20. Annealed Gibbs sampling is employed for inference.
\item[Offline inference] The dd-CRP model with offline inference procedure (described in \autoref{sec:inference}).
\item[Online inference] The dd-CRP model with online inference procedure (described in \autoref{sec:online}).
\end{description}

For the online inference implementation, we set the size of window and number of iterations to five days and 500 respectively. For the baseline, the parameter of the Dirichlet prior was set to a vector  of $0.5$ for each cluster. These values were chosen through 10-fold cross validation.

To measure the quality of the clusterings obtained by these models, we compare the average weighted and unweighted F-measures for 55 TREC topics, using the evaluation scripts from the TREC TTG task. Overall results are shown in \autoref{tab:sim}. The \textsc{online model} has the best weighted F1 score, outperforming the offline version of the same model, even though its inference procedure is an approximation to the \textsc{offline model}. It may be that its approximate inference procedure discourages long-range linkages, thus placing a greater emphasis on the temporal dimension. Both models were trained over 500 iterations, and the \textsc{online model} was 30\% faster to train than the offline model.

Compared to the other 2014 TREC TTG systems, our dd-CRP models are competitive. Both models outperform all but one of the fourteen submissions on the unweighted $F_1$ metric, and would have placed fourth on the weighted $F^{w}_1$ metric. Note that the TREC evaluation scores both clustering quality and retrieval. We use only the baseline retrieval model, which achieved a mean average precision of 0.31. The competing systems shown in \autoref{tab:sim} all use retrieval models that are far superior: the retrieval model for top-ranked PKUICST team (line 4) achieved a mean average precision (MAP) of 0.59~\cite{lv2014pkuicst}, and the QCRI~\cite{magdy2014qcri} and and hltcoe~\cite{xu2014hltcoe} teams (lines 5 and 6) used retrieval models with MAP scores of at least 0.5. Bayesian dd-CRP storyline clustering was competitive with these timeline generation systems despite employing a far worse retrieval model, so improving the retrieval model to achieve parity with these alternative systems seems the most straightforward path towards better overall performance.

\begin{table*}
\centering
\begin{tabular}{l l l l l l}
    \toprule
    Model & Rec. & Rec.$^{w}$ & Prec. & $F_1$ & $F_1^{w}$ \\ \midrule
\textit{dd-CRP clustering models} \\
1. \textsc{baseline} & 0.14 & 0.27 & 0.33 & 0.20 & 0.30 \\
2. \textsc{offline} & 0.32 & 0.47 & 0.27 & 0.29 & 0.34 \\
3. \textsc{online} & 0.34 & 0.55 & 0.26 & 0.29 & 0.35 \\[8pt]
\textit{Top systems from Trec-2014 TTG}\\
4. \texttt{TTGPKUICST2}~\cite{lv2014pkuicst} & 0.37 & 0.58 & 0.46 & 0.35 & 0.46 \\
5. \texttt{EM50}~\cite{magdy2014qcri} & 0.29 & 0.48 & 0.42 & 0.25 & 0.38 \\
6. \texttt{hltcoeTTG1}~\cite{xu2014hltcoe} & 0.40 & 0.59 & 0.34 & 0.28 & 0.37 \\
    \bottomrule
\end{tabular}
\caption{Performance of Models in the TREC 2014 TTG Task. Weighted recall and $F_1$ are indicated as Rec.$^{w}$ and $F_1^w$.}
\label{tab:sim}
\end{table*}

\section{Inference}
\label{sec:inference}
The key sampling equation for the dd-CRP is the posterior likelihood,
\begin{align*}
\Pr(c_i = j \mid \vec{c}_{-i}, \vw) \propto & \Pr(c_i = j) P(\vw \mid \vec{c}).
\end{align*}
The prior is defined in Equation~\ref{eq:prior}. Let $\ell$ represent the likelihood under the partitioning induced when the link $c_i$ is cut. Now, the likelihood term has two cases: in the first case, $j$ is already in the same connected component as $i$ (even after cutting the link $c_i$), so no components are merged by setting $c_i = j$. In this case, the likelihood $P(\vw \mid c_i = j)$ is exactly equal to $\ell$. In the second case, setting $c_i = j$ causes two clusters to be merged. This gives the likelihood,
\begin{align*}
P&(\vw \mid c_i = j, \vec{c}_{-i}) \\
& \propto \frac{P(\{\vw_k : \zc_k = \zc_j \vee \zc_k = \zc_i\})}
{P(\{\vw_k : \zc_k = \zc_i\})
P(\{\vw_k : \zc_k = \zc_j\}) },
\label{eq:likelihood}
\end{align*}
where the constant of proportionality is exactly equal to $\ell$. Each of the terms in the likelihood ratio is a Dirichlet Compound Multinomial likelihood. This likelihood function is itself a ratio of gamma functions; by eliminating constant terms and exploiting the identity $\Gamma(x+1) = x \Gamma(x)$, we can reduce the number of Gamma function evaluations required to compute this ratio to the number of words which appear in \emph{both} clusters $\zc_i$ and $\zc_j$. Words that occur in neither cluster can safely be ignored, and the gamma functions for words which occur in exactly one of the two clusters cancel in the numerator and denominator of the ratio. Note also that we only need compute the likelihood for $c_i$ with respect to each cluster, not for every possible follower link.

\subsection{Online inference}
\label{sec:online}
While we make every effort to accelerate the computation of individual Gibbs samples, the complexity of the basic algorithm is superlinear in the number of instances. This is due to the fact that each sample requires computing the probability of instance $i$ joining every possible cluster, while the number of clusters itself grows with the number of instances (this growth is logarithmic in the Chinese Restaurant Process). Scalability to the streaming setting therefore requires more aggressive optimizations.

To get back to linear time complexity, we employ a fixed-lag sampling procedure~\cite{doucet2000sequential}. After receiving instance $i$, we perform Gibbs sampling only within the fixed window $[t_i - \tau, t_i]$, leaving $c_j$ fixed if $t_j < t_i - \tau$. This approximate sampling procedure implicitly changes the underlying model, because there is no possibility of linking $i$ to a later message $j$ if the time gap $t_j - t_i > \tau$.

Since we are only interested in obtaining a single storyline clustering --- rather than a full Bayesian distribution over clusterings --- we perform annealing for samples towards the end of the sampling window. Specifically, we set the temperature to $\gamma = 2.0$  and exponentiate the sampling likelihood by the inverse temperature~\cite{geman1984stochastic}. This has the effect of interpolating between probabilistically-correct Gibbs sampling and a hard coordinate-ascent procedure.

\subsection{Hyperparameter estimation}
\label{sec:hyper}
The model has three parameters to estimate:
\begin{itemize}
\setlength\parskip{0pt}
\setlength\parsep{0pt}
\setlength\itemsep{0pt}
\item $\alpha$, the concentration parameter of the dd-CRP
\item $a$, the offset of the distance function
\item $\eta$, the scale of the symmetric Dirichlet prior.
\end{itemize}
We interleave maximization-based updates to these parameters with sampling, in a procedure inspired by Monte Carlo Expectation Maximization~\cite{wei1990monte}. Specifically, we compute gradients on the likelihood $P(\vc)$ with respect to $\alpha$ and $a$, and take gradient steps after every fixed number of samples. For the symmetric Dirichlet parameter $\eta$, we employ the heuristic from~\newcite{minka2000estimating} by setting the parameter to $\eta = \frac{(K-1)/2}{\sum_k\log p_k}$, where $K$ is the number of words that appear exactly once, and $p_k$ is the probability of choosing the $k^{th}$ word from the vocabulary under the unigram distribution for the entire corpus.

\section{Introduction}
A long-standing goal for information retrieval and extraction is to identify and group textual references to ongoing events in the world~\cite{allan2002topic}. Success on this task would have applications in personalized news portals~\cite{gabrilovich2004newsjunkie}, intelligence analysis, disaster relief~\cite{vieweg2010microblogging}, and in understanding the properties of the news cycle~\cite{leskovec2009meme}. This task attains a new importance in the era of social media, where citizen journalists can document events as they unfold~\cite{lotan2011arab}, but where repetition and untrustworthy information can make the reader's task especially challenging~\cite{becker2011beyond,marcus2011twitinfo,petrovic2010streaming}.

A major technical challenge is in fusing information from two heterogeneous data sources: textual content and time. Two different documents about a single event might use very different vocabulary, particularly in sparse social media data such as microblogs; conversely, two different sporting events might be described in nearly identical language, with differences only in the numerical outcome. Temporal information is therefore critical: in the first case, to find the commonalities across disparate writing styles, and in the second case, to identify the differences. A further challenge is that unlike in standard document clustering tasks, the number of events in a data stream is typically unknown in advance. Finally, there is a high premium on scalability, since online text is produced at a high rate.

Due to these challenges, existing approaches for combining these modalities have been somewhat heuristic, relying on tunable parameters to control the tradeoff between textual and temporal similarity. In contrast, the Bayesian setting provides elegant formalisms for reasoning about latent structures (e.g., events) and their stochastically-generated realizations across text and time. In this paper, we describe one such model, based on the distance-dependent Chinese Restaurant Process (dd-CRP; Blei and Frazier, 2011)\nocite{blei2011distance}. This model is distinguished by the neat separation that it draws between textual content, which is treated as a stochastic emission from an unknown Multinomial distribution, and time, which is modeled as a prior on graphs over documents, through an arbitrary \emph{distance function}. However, straightforward implementations of the dd-CRP are insufficiently scalable, and so the model has been relatively underutilized in the NLP literature~\cite{titov2011bayesian,kim2011accounting,sirts2014pos}. We describe improvements to Bayesian inference that make the application of this model feasible, and present encouraging empirical results on the Tweet Timeline Generation task from TREC 2014~\cite{lin2014overview}.

\section{Model}
The basic task that we address is to group short text documents into an unknown number of storylines, based on their textual content and their temporal signature. The textual content may be extremely sparse --- the typical Tweet is on the order of ten words long --- so leveraging temporal information is crucial. Moreover, the temporal signal is multiscale: in the 24-hour news cycle, some storylines last for less than an hour, while others, like the disappearance of the Malaysian Airlines 370 plane in 2014, continue for weeks or months. In some cases, the temporal distribution of references to a storyline will be unimodal and well-described by a parametric model~\cite{marcus2011twitinfo}; in other cases, it may be irregular, with bursts of activity followed by periods of silence~\cite{he2007analyzing}. Finally, it will be crucial to produce an implementation that scales to large corpora.

The distance-dependent Chinese Restaurant Process (dd-CRP) meets many of these criteria~\cite{blei2011distance}. In this model, the key idea is that each instance (document) $i$  ``follows'' another instance $c_i$ (where it is possible that $c_i = i$), inducing a graph. We can compute a partitioning over instances by considering the connected components in the undirected version of the follower graph; these partitions correspond to ``tables'' in the conventional ``Chinese Restaurant'' analogy~\cite{aldous1985exchangeability}, or to clusters. The advantage of this approach is that it is fundamentally non-parametric, and it introduces a clean separation between the textual data and the covariates: the text is generated by a distribution associated with the partition, while the covariates are associated with the following links, which are conditioned on a distance function.

The distribution over follower links for document $i$ has the following form,
\begin{equation}
\Pr(c_i = j) \propto
\begin{cases}
  f(d_{i,j}), & i \neq j\\
  \alpha, & i = j,\\
\end{cases}
\label{eq:prior}
\end{equation}
where $d_{i,j}$ is the distance between units $i$ and $j$, and $\alpha > 0$ is a parameter of the model. Large values of $\alpha$ induce more self-links and therefore more fine-grained partitionings. Since we are concerned with temporal covariates, we define the distance function as follows:
\begin{equation}
f(d_{i,j}) = e^{\frac{-|t_i - t_j|}{a}}.
\end{equation}
Thus, the likelihood of document $i$ following document $j$ decreases exponentially as the time gap $|t_i - t_j|$ increases.

The text of each document $i$ is represented by a vector of word counts $\vw_i$. The likelihood distribution is multinomial, conditioned on a parameter $\theta$ associated with the partition to which document $i$ belongs. By placing a Dirichlet prior on $\theta$, we can analytically integrate it out. Writing $\vzc$ for the cluster membership induced by the follower graph $\vc$, we have:
\begin{small}
\begin{align}
P(\vw \mid \vc; \eta) = & \prod_k P(\{ \vw_i : \zci = k \}; \eta)\\
 = &\prod_k \int_\theta P(\{ \vw_i : \zci = k \} \mid \theta) P(\theta ; \eta) d\theta
\label{eq:likelihood}
\end{align}
\end{small}

Given a multinomial likelihood $P(\vw \mid \theta)$ and a (symmetric) Dirichlet prior $P(\theta \mid \eta)$, this integral has a closed-form solution as the Dirichlet-Multinomial distribution (also known as the multivariate Polya distribution). The joint probability is therefore equal to the product of \autoref{eq:prior} and \autoref{eq:likelihood},
\begin{align}
P(\vw, \vc) = \prod_i P(c_i ; \alpha, a) \prod_k P(\{\vw_i : \zci = k\}; \eta).
\label{eq:joint}
\end{align}
The model has three hyperparameters: $\alpha$, which controls the likelihood of self-linking, and therefore affects the number of clusters; $a$, which controls the time scale of the distance function, and therefore affects the importance of the temporal dimension to the resulting clusters; and $\eta$, which controls the precision of the Dirichlet prior, and therefore the importance of rare words in the textual likelihood function. Estimation of these hyperparameters is described in \autoref{sec:hyper}.
\section{Related work}
Topic tracking and first-story detection are very well-studied tasks; space does not permit a complete analysis of the related work, but see~\cite{allan2002topic} for a summary of ``first generation'' research. More recent non-Bayesian approaches have focused on string overlap~\cite{suen2013nifty}, submodular optimization~\cite{shahaf2012trains}, and locality-sensitive hashing~\cite{petrovic2010streaming}.
In Bayesian storyline analysis, the seminal models are Topics-Over-Time~\cite{wang2006topics}, which associates a parametric distribution over time with each topic~\cite{ihler2006adaptive}, and the Dynamic Topic Model~\cite{blei2006dynamic}, which models topic evolution as a linear dynamical system~\cite{nallapati2007multiscale}. Later work by \newcite{diao2012finding} offers a model for identifying ``bursty'' topics, with inference requiring dynamic programming. All these approaches require the number of topics to be identified in advance. \newcite{kim2011accounting} apply a distance-dependent Chinese Restaurant \emph{Franchise} for temporal topic modeling; they evaluate using predictive likelihood rather than comparing against ground truth, and do not consider online inference.

The Infinite Topic-Cluster model~\cite{ahmed2011unified} is non-parametric over the number of storylines, through the use of the recurrent Chinese Restaurant Process (rCRP). The model is substantially more complex than our approach. Unlike the dd-CRP, the rCRP is Markovian in nature, so that the topic distribution at each point in time is conditioned on the previous epoch (or, at best, the previous $K$ epochs, with complexity of inference increasing with $K$). This Markovian assumption creates probabilistic dependencies between the topic assignment for a given document and the documents that follow in subsequent epochs, necessitating an inference procedure that combines sequential Monte Carlo and Metropolis Hastings, and a custom data structure; this inference procedure was complex enough to warrant a companion paper~\cite{ahmed2011online}. The rCRP is also employed by Diao and Jiang (2013, 2014)\nocite{diao2013unified,diao2014recurrent}.
In contrast, the dd-CRP makes no Markovian assumptions, and efficient inference is possible through relatively straightforward Gibbs sampling in a fixed window.


\section{Introduction}
The origin of high-energy emission has long been a key question in Active Galactic Nuclei (AGN)
physics. While they represent only a minority of AGN, blazars constitute a unique laboratory to
probe the acceleration processes in relativistic jets emanating from the central engine.
Since its launch in 2008, the {\it Fermi}-LAT (Large Area Telescope) has revolutionized
our knowledge of the $\gamma$-ray sky with a combination of high sensitivity, wide
field-of-view, and a large energy range (about 20 MeV to more than 300 GeV). Its nominal
sky-survey operating mode has enabled a continuous monitoring of the complete $\gamma$-ray
sky at an unprecedented level, making remarkable discoveries. The third {\it Fermi}-LAT
catalog \citep[3FGL,][]{acero2015}, based on the first 48 months of data, consists of more than 3000
$\gamma$-ray  sources including $\sim$60$\%$ AGN - most of which are blazars \citep[$\sim$98$\%$][]{ackermann2015}.
Blazars therefore form a major population of the $\gamma$-ray sky and are potential targets
for understanding the physical processes responsible for production of high-energy emission.
Moreover, blazars also account for the majority of the high energy $\gamma$-ray background
above 100~GeV \citep{ajello2015}.

BL Lacertae objects (BL~Lacs) and flat spectrum radio quasars (FSRQs) are clubbed together and called blazars.
In spite of the dissimilarity of their optical spectra -- FSRQs show strong broad emission lines, while
BL~Lacs have only weak or no emission lines in their optical spectra -- they share the same peculiar continuum
properties (strong variability and polarization properties). They are characterized by powerful non-thermal
emission ranging from radio to the $\gamma$-ray bands and exhibit strong variability, over a variety
of time scales from minutes to months and often the radio components detected in VLBI observations
exhibit superluminal motion. These properties are interpreted as resulting from the emission of
high-energy particles accelerated within a relativistic jet aligned with the direction of sight.
In many cases, the bulk of the observed flux emerges in $\gamma$-ray photons in GeV (TeV for some sources)
energy range populating a major fraction of the high-energy sky \citep{ackermann2015, acero2015}. One of the most
intriguing and challenging quests of current astrophysics is to understand the physical conditions and processes
responsible for the  production of high-energy particles, and the emission of $\gamma$-rays.

\section{Particles accelerated in relativistic jets}

Detailed investigation of multi-wavelength flux and spectral variability of individual sources including
cross-band (radio, optical, X-ray, $\gamma$-ray, and polarization), relative timing analyses of
outbursts, and/or VLBI component ejection/kinematics suggest that $\gamma$-rays are associated with
the compact regions of relativistic jets energized by the central SMBHs \citep[e.g.,][and references therein]{agudo2011, jorstad2001, jorstad2010, marscher2010, schinzel2012, rani2013a, rani2015}. For
individual objects, the close association of the $\gamma$-ray flare
with a continuous change of the optical polarization angle provides evidence for the presence of highly ordered
magnetic fields in the regions of $\gamma$-ray production \citep{abdo2010c, marscher2008}.
Recent studies (at least for some blazars) suggest a
significant correlation between the variations of $\gamma$-ray flux and the direction of the jet outflow
in the sub-parsec scale region \citep{rani2014}, which implies that the observed inner jet morphology has
a strong connection with the observed $\gamma$-ray flares.  However, the fine details of acceleration
mechanisms and radiation processes are still missing.

Blazars exhibit a characteristic bimodal spectral energy distribution (SED). The SED low-energy
peak (from radio to optical/UV and X-rays for some sources) is commonly interpreted as synchrotron emission
from high-energy electrons, while no consensus on the origin of the high-energy peak (from X-rays to
$\gamma$-rays) has been reached. In the framework of leptonic models, this second peak corresponds
to photons produced via inverse-Compton scattering of soft seed photons by the same electrons, with
these seed photons originating either from within the jet (Synchrotron Self-Compton, SSC) or from regions
external to the jet, such as the accretion disk, the dust torus or the broad emission line clouds
(External Compton, EC). In hadronic models, the high energy peak is due to photons produced in
photo-hadron interactions or by proton or charged $\mu$/$\pi$ synchrotron emission and subsequent
cascades \citep[e.g.,][]{bottcher2007, bottcher2012, bottcher2013}.
Concerning the high-energy emission, leptonic models are more favored in the literature;
however, hadronic processes remain a viable option.

No clear consensus has been achieved on the location of emission regions. The radio--$\gamma$-ray
correlations \citep{fuhrmann2014, marscher2008, marscher2010, rani2013a}  suggests the latter
leading the former, which puts the location of high-energy emission regions closer to the central
black hole.   The observed $\gamma$-ray spectral breaks at a few GeVs require the emitting region to
be located within the BLR if the photon-photon pair production scenarios are responsible for their
origin \citep{abdo2009, finke2010, rani2013b, rani2013_3c273, poutanen2010}. On the other hand, coincidence
of $\gamma$-ray flares with the appearance of new components from the base of the jet
indicates distances longer than few parsecs \citep{jorstad2013, agudo2011, schinzel2012}.
Interaction of moving features with the stationary features can amplify the magnetic field and
accelerate particles \citep{marscher2014}.  Coincidence of $\gamma$-ray flares with the passage of moving components through
stationary features in jets has been observed in several sources \citep{schinzel2012, rani2015}, which implies
that $\gamma$-ray flares are produced downstream of the core at a distance of more than a few parsecs.

Particle acceleration in relativistic jets is also an open question. Relativistic shocks seem to be
effective in low magnetization regimes; however in the case of high magnetization, shocks are weak and poor
accelerators. Particle acceleration through second order stochastic processes is rather slow in
highly magnetized flows; but it is still a viable option \citep{blanford1973}. In highly magnetized environments,
particle acceleration via magnetic reconnection is very likely to occur \citep{kagan2014}. Most likely
all of these processes are in action on different scales in jets. Closer to the central engines, magnetic
reconnection could inject energetic particles in compact regions to produce rapid $\gamma$-ray flares.
Further downstream the jet  in low magnetized environments, the observed $\gamma$-ray variations
could be a natural consequence of relativistic shocks.  Magnetoluminescence and/or electromagnetic detonation
is lately proposed to explain the fastest observed TeV flares.  The basic idea is that the tightly
wound magnetic field in the inner jet regions is subject to dynamical instability which could be responsible
for the rapid $\gamma$-ray variability. In other words, impulsive particle acceleration can take place while
the magnetic flux ropes untangle \citep{blanford2015}, which could be responsible for high-energy flares.

 \begin{figure}
\begin{minipage}{0.5\linewidth}
 \includegraphics[scale=0.45,angle=0, trim=130 0 120 2, clip]{s50716_map.eps}
\end{minipage}\hfill
\begin{minipage}{0.5\linewidth}
\includegraphics[scale=0.323,angle=0,trim=0 0 0 2, clip]{PA_var.eps}
\includegraphics[scale=0.323,angle=0,trim=0 0 1.8 2, clip]{core_flx_var.eps}
\includegraphics[scale=0.323,angle=0,trim=5 0 0 2, clip]{gamma_var.eps}
\end{minipage}\hfill
\caption{Left: The 86~GHz VLBI image of S5 0716+714 observed on February 18, 2013. Right: (a)
Orientation variations in inner region ($\leq$0.2~mas) of the jet, (b) the 7~mm VLBI core flux density
variations,  and (c) the monthly averaged $\gamma$-ray photon flux light curve.     }
\label{fig1}
\end{figure}

\section{Causal connection between gamma-ray emission and jet morphology}
Our recent studies on the BL Lac S5 0716+714 suggest a significant correlation between the
variations of $\gamma$-ray flux and the direction of the jet outflow in sub-parsec scale region
\citep[see][for details]{rani2014}. The core flux density and inner-jet orientation (position angle) variations
were investigated using the 7~millimeter Very Long Baseline Interferometry (VLBI) observations of
S5 0716+714 for a time period between August 2008 and September 2013. In Fig.\ \ref{fig1} (left),
we show an example map of the source. The variations
in the core flux density and orientation of the sub-parsec scale jet i.e.\ position angle (PA) are
shown in Fig.\ \ref{fig1} (a and b).
The $\gamma$-ray flux variations (Fig.\ \ref{fig1} c) over the same time period were examined using the
observations obtained  by the {\it Fermi}-LAT (Large Area Telescope).

The discrete cross-correlation function \citep[DCF,][]{edelson1988} method was used to investigate the correlation between
the observed $\gamma$-ray and jet morphology (core flux density and PA) variations. The DCF analysis
results are presented in Fig.\ \ref{fig2}. The analysis suggests that the high-energy radiation has
a tight correlation between the VLBI core flux density and the direction of the jet flow flow (Fig.\ \ref{fig2}).
We observed a time lag of 82$\pm$39 days between$\gamma$-ray and core flux density variations, which
puts the location of $\gamma$-ray emission region upstream of the 7~mm VLBI core by 3.8$\pm$1.9~pc
(de-projected). Given the location of the 7~mm core at a distance of $\sim$6.5~pc \citep{rani2015},
we conclude that the $\gamma$-ray emission region is at a distance of 2.7$\pm$1.9~pc from the central engine. The
correlated variations observed in the innermost PA and the high-energy flux of
S5 0716+714 add a unique new point to constrain the location and origin of high-energy emission in
blazars.

  \begin{figure}
\begin{minipage}{0.7\linewidth}
\includegraphics[scale=0.22,angle=0,trim=0 0 0 0, clip]{dcf_gamma_core.eps}
\includegraphics[scale=0.22,angle=0,trim=0 0 0 0, clip]{dcf_gamma_pa.eps}
\end{minipage}
\begin{minipage}{0.3\linewidth}
  \caption{Cross-correlation analysis curves: $\gamma$-ray vs.\ core flux density (a) and $\gamma$-ray vs.\ position
angle (b). The dashed lines show the 99$\%$ confidence levels.    }
\end{minipage}
\label{fig2}
    \end{figure}

\section{Future Perspectives}
High-frequency VLBI is a promising technique to pinpoint the high-energy emission sites.
With the current GMVA (Global mm-VLBI Array) resolution at 86~GHz frequency ($\sim$50~$\mu$as), one could
scale down to less than a thousand Schwarzschild radii. The Event Horizon Telescope (EHT) will offer an
angular resolution of $\sim$20~$\mu$as, which means that the future of high-resolution VLBI is very bright.
Participation of ALMA will bring a revolutionary leap in capabilities for mm-VLBI. The {\it Fermi} mission with its
extraordinary capabilities will keep observing the $\gamma$-ray sky at least for the next few years.
In the upcoming years, TeV astronomy will also grow making several more new and major discoveries.
Not to mention, high-energy polarization missions (Astro-H, GEMS, X-calibur, PoGOLite, Polar, Harpo,
and many more) are also on their way to teach us more about high-energy emission models.

\ack
The $Fermi$-LAT Collaboration acknowledges support from a number of agencies and institutes for both
development and the operation of the LAT as well as scientific data analysis. These include NASA and
DOE in the United States, CEA/Irfu and IN2P3/CNRS in France, ASI and INFN in Italy, MEXT, KEK, and JAXA
in Japan, and the K.~A.~Wallenberg Foundation, the Swedish Research Council and the National Space Board
in Sweden. Additional support from INAF in Italy and CNES in France for science analysis during the
operations phase is also gratefully acknowledged.~~This study makes use of 43 GHz VLBA data from the
VLBA-BU Blazar Monitoring Program (VLBA-BU-BLAZAR; http://www.bu.edu/blazars/VLBAproject.html).
This research has made use of data obtained with the Global Millimeter VLBI Array (GMVA), which consists of
telescopes operated by the MPIfR, IRAM, Onsala, Mets\"ahovi, Yebes and the VLBA. The data were correlated at
the correlator of the MPIfR in Bonn, Germany. The VLBA is an instrument of the National Radio Astronomy
Observatory, a facility of the National Science Foundation operated under cooperative agreement by Associated
Universities, Inc.