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

The next generation of large ground-based telescopes hold the potential to open a new realm of discovery at the Galactic center, both in fundamental physics and astrophysics. The precision of the measurement of stellar orbits has  been anticipated to be powerful enough to detect very small deviations in the orbital parameters that can be used to test General Relativity and alternative theories of gravity.

While the potential for discovery has been long recognized at the Galactic center, realistic estimates the sensitivity of such experiments has not been well studied for the next generation of extremely-large ground-based telescopes (TMT, GMT, ELT). The difficulty stems from the fact that the measurement of orbits at the Galactic center requires understanding 5 major components: (1) the astrophysical content near the black hole (e.g. the predicted number of stars as well as their physical properties); (2) estimates of the performance of the adaptive optics system; (3) estimates of the performance of the science instrument and the specific configuration that will be used for observations; (4) a model for the cadence of observations and uncertainties in positional and radial velocity measurements; (5) a framework for estimating the uncertainty from fitting many orbits simultaneously to constrain physical parameters such as post-Newtonian effects on the orbits. The lack of a comprehensive study has hindered our ability to quantitatively estimate and optimize the sensitivity to tests of various physical models of gravity.

Here, we present some of the first realistic simulations of observations of stars at the Galactic center with TMT, taking into account all 5 major components listed above. We use TMT as an example of this work because we have extensive simulations from the adaptive optics system (NFIRAOS) and integral-field spectrograph and imager (IRIS). A comparable study could be done with GMT and ELT as well. In Section \ref{sec:ao}, we present the estimates of the performance and point spread function delivered by NFIRAOS. We discuss the simulations of the instrumental performance of IRIS on Galactic center sources. Using our existing knowledge of the stellar population and distribution of stars, we predict statistically the number and types of stars that will be detected by TMT in Section \ref{sec:astrophysics}. We present our observational cadence tool and the model for the effect of source confusion used to determine the sensitivity to physical models in Section \ref{sec:models}. We discuss our resulting constraints on models such as the precession of the periaspe of the orbit of S0-2 in Section \ref{sec:results}. Section \ref{sec:conclusions} states our conclusions and prospects for future improvements.

\section{ADAPTIVE OPTICS AND IRIS INSTRUMENT}
\label{sec:ao}

We use simulations of the AO point spread function (PSF) from the TMT NFIRAOS team in order to simulate the observations at the Galactic center. These PSFs are described in \citenum{2014AJ....147...93D} and \citenum{2013aoel.confE..83Y}. Bases on observing conditions at the Galactic center from Mauna Kea (airmass, atmospheric conditions, guide star configuration, etc.), the PSF is predicted to have Strehl ratio of about 0.7 at K-band. Since the simulations are limited to the central $1^{\prime\prime}\times1^{\prime\prime}$, spatial PSF variations should be very small so we use the same PSF for all sources.

We use the IRIS imaging and spectroscopic simulator to produce images and spectra at the Galactic center. The properties of the simulator is described in \citenum{2014AJ....147...93D}. The simulator takes into account different noise sources as well as the throughput and sensitivity predictions for the AO system and the IRIS instrument in order to simulate realistic observations.

\subsection{Astrometric and Radial Velocity Precision}
\label{sec:errors}
With its larger aperture and better AO correction, the astrometric precision of TMT is expected to be much improved compared to current 8-10 m telescopes. For the simulations here, we use the results of astrometric experiments performed by \citenum{2013aoel.confE..83Y}. We expect that the astrometric performance of TMT will likely be in the range of 25 to 50 $\mu$as for well measured stars, compared to about 150 $\mu$as for observations today \cite{2016ApJ...830...17B}.

The radial velocity precision is a function of both the signal to noise ratio (SNR) as well as the spectral resolution. We determine the expected radial velocity precision by simulating spectra at different SNR and at three different spectral resolution: R = 4000, 8000, and 10,000, corresponding to the potential modes for IRIS. We find that early-type stars (A-type and earlier) K-band have higher radial velocity uncertainties than late-type stars (K and M-type) because early-type stars have much fewer lines than late-type stars. At a given spectral resolution and SNR, late-type stars generally have a factor of $\sim5$ times smaller uncertainties (Fig. \ref{fig:rv_err}). We use the predictions for the radial velocity uncertainties in our orbit simulations below. One of the biggest advantage of TMT will be the increase in spectroscopic sensitivity. Today, we are spectroscopically about 50\% complete at Kp = 16.0 in this region, with the faintest sources with a measured radial velocity is about K = 17 \cite{2013ApJ...764..154D,2016ApJ...830...17B,2017ApJ...837...30G}, TMT will likely measure radial velocities down to K = 20-22. Likely, all the known short period stars without radial velocity measurements today will have measurements with TMT.

\begin{figure}[hbt]
\center
\includegraphics[width=4.0in]{fig1.pdf}
\caption{The result of simulations showing expected radial velocity as a function of SNR at R = 4000, 8000, and 10,000 at K-band. Observations of early-type stars (blue) will have larger uncertainties compared to late-type stars (red) because they have fewer spectral lines in K-band.}
\label{fig:rv_err}
\end{figure}

\section{WHAT LIES AROUND A SUPERMASSIVE BLACK HOLE?}
\label{sec:astrophysics}

While current 8-10 m telescopes have observed over 100 stars within just 1$^{\prime\prime}$ from supermassive black hole, future extremely large telescopes should detect many more stars -- both fainter and closer to the black hole. Currently, the detection of sources are limited by stellar confusion with both other observed stars as well as unresolved sources. This confusion limit is set by both the angular resolution of the telescope and the performance of the AO system. Typical imaging observations of stars at the Galactic center from the Keck Telescopes have Strehl ratios of 20 to 35\% with a core of about 60 mas at Kp (2.2 $\mu$m). These images represent the state of our knowledge in terms of the number and types of stars in this region.

In order to simulate observations with TMT, we need to extrapolate our current knowledge of the stellar population and distribution of stars immediately around the black hole to simulate the properties of the underlying stars. We know today that there are two major populations of stars around the supermassive black hole: an old population from about 1 Gyr ago, and a young population formed just 4-6 Myr ago\cite{2009ApJ...703.1323D,2011ApJ...741..108P,2013ApJ...764..155L,2015ApJ...808..106S}. These populations are currently best separated spectroscopically using integral-field spectroscopy with AO, since they have about the same luminosity and color in the near-infrared (Fig. \ref{fig:cmd}). Based on the luminosity function of the currently observed young stars, \citenum{2013ApJ...764..155L} was able to fit the mass function and age of the young stellar cluster. The observations are complete only to about early B-type stars so they represent the high-mass end ($> 15 M_\odot$) of the luminosity function. Similarly, for the old stars, we can only spectroscopically observe the evolved stars which have high luminosities. Based on these stars, \citenum{2011ApJ...741..108P} derived a star formation history and mass function for the old population. Using a stellar population model, we can predict the luminosity function at fainter magnitudes for both of these populations, normalized such that the currently observed number of stars match the simulated stars at magnitudes where current observations are complete. All currently known stars, regardless of whether they have been spectroscopically identified are included in the simulations.  We simulate to a depth of K = 22 mag, about 8 magnitudes fainter than the current spectroscopic limit (Kp = 15.5) and about 4 magnitudes fainter than the imaging limit (Kp = 18.0).

\begin{figure}[hbt]
\center
\includegraphics[width=6.0in]{fig2.pdf}
\caption{\textbf{Left:} Theoretical color magnitude diagram of stars at the distance of the Galactic center. The H-K color is the intrinsic color of the stars. The small different in near-infrared colors and the large variations in dust extinction makes it difficult to differentiate between stars with photometry alone. Current spectroscopy from 8-10 m telescopes are limited only to the most massive young stars or red giants. TMT will provide many magnitudes of improvement in spectroscopic sensitivity. \textbf{Right:} An example of a spectrum of a short period star, S0-102 observed with Keck today. We are unable to obtain a radial velocity measurement or a spectral type of this star. With TMT, we predict that high SNR spectra will be possible for these stars. Obtaining radial velocity measurements of these stars are critical for accurate orbit determination.}
\label{fig:cmd}
\end{figure}

In order determine the distribution of stars in this region, we also use the currently observed density profiles of these two population to extrapolate the density of young and old stars as a function of distance from the supermassive black hole. Inherent in such an assumption is that the density profile of the unobserved sources follow the same distribution as the brighter sources. This may not be the case because of effects such as mass segregation or tidal stripping or collisions of red giants could modify the brighter population relative to the fainter one. However, since these factors are currently observationally unconstrained, we will make the simplifying assumption that the fainter population follows the brighter one. The density profile within this region of interest has the young stars steeply rising with an observed surface density profile of $\Sigma_{yng}(R) \propto R^{-1.7}$ compared to the old population with $\Sigma_{old}(R) \propto R^{-0.5}$, where R is the projected distance of the stars from the black hole. This steep density profile and the luminosity function results in more young stars than old stars within 0.04 pc of the black hole.

The final astrophysical component is the orbital parameters for the simulated stars. Current observations of stars within the region of interest (about 0.04 pc from the black hole) show that they have orbits with orientations ($i$ - inclination, $\Omega$ - longitude of the ascending node, $\omega$ - argument of periapsis) that are isotropically distributed on the plane of the sky. This result is expected because the dynamical relaxation times for stars around the supermassive black hole should be fairly fast. The relevant time scales in this region are the timescale from General Relativistic precession and vector resonant relaxation, both of which should relax the stars within several million years \cite{2007ApJ...666..919Y}. We therefore uniformly draw the stars with randomly oriented orbital planes and from an eccentricity ($e$) distribution of the form $P(e) \propto e^{2}$ \cite{2017ApJ...837...30G}. For the semi-major axis ($a$), we draw from the radial distribution of stars discussed in the previous paragraph. We also randomly choose $T_o$, the time of closest approach of the star to the black hole.

When combined, the astrophysical simulations of stars around the supermassive black hole contains stellar parameters such as the effective temperature, surface gravity, and age, and 6 orbital parameters ($e$, $i$, $\Omega$, $\omega$, $a$, and $T_o$). Deviations from the Keplerian orbits are included as additional free parameters, which depend on the theory of gravity being tested.

\begin{figure}[hbt]
\center
\includegraphics[width=6.0in]{fig3.pdf}
\caption{Simulation of the central 0.5 arcsecond around the supermassive black hole Sgr A* (white cross) with current adaptive optics with the Keck Telescopes (left) compared to that with TMT (right). The greater resolving power of TMT along with better AO correction will likely allow us to detect a factor of $> 10$ times more stars in this region. We simulate how these orbit can be used simultaneously in a joint analysis to constraint tests of gravity and other science cases.}
\label{fig:keck_tmt}
\end{figure}

\begin{figure}[hbt]
\center
\includegraphics[width=5.0in]{fig4.pdf}
\caption{Overlay of a sample orbits from the simulations described in Section \ref{sec:astrophysics}. TMT will likely discover stars with orbital periods below 5 years given its increased sensitivity and angular resolution compared to telescopes today.}
\label{fig:tmt_orbits}
\end{figure}

\section{ORBIT SCIENCE WITH TMT}
\label{sec:models}
We test several representative science cases that will greatly benefit from the sensitivity and angular resolution of the next generation of large ground-based telescopes: (1) measurement of the precession of the periapse predicted by General Relativity; (2) measurement of the extended dark cusp of stars and stellar remnants around the black hole; (3) measuring the black hole mass and distance to the Galactic center.

Our goal is to perform a sensitivity analysis using realistic simulations of short-period stars in the TMT era using the following methodology:
\begin{itemize}
	\item Identify and simulate the orbital parameters of stars that will be used in our simulations.
	\item Simulate an observation plan.
	\item Model the observational astrometric and radial velocity uncertainties.
	\item Model the stellar confusion of sources and their impact on positional and radial velocity measurements.
	\item Using all the observations epochs and related uncertainties, compute an estimate of the sensitivity to physical effects on the orbits.
	\end{itemize}
The time baseline considered in this work starts in 1995, when high-angular resolution measurements began at Keck, up to 2047 -- after 20 years of TMT operations. We consider observations using a combination of two types of telescopes: (1) observations made using 8-10 m telescopes today  and (2) observations made using TMT--like telescope that starting in 2027. The steps of our simulations are described with more details below.

\subsection{Number of stars simulated, orbital parameters}
The sample of stars used in this simulation includes 17 well measured stars whose orbital parameters are reported by \citenum{2012Sci...338...84M,2016ApJ...830...17B,2017ApJ...837...30G} as well as 105 simulated stars that TMT is expected to detect in the central region (Section~\ref{sec:astrophysics}). All these stars are included in the source confusion analysis which is described in Section~\ref{sec:sourceconfusion}. For the purpose of the sensitivity analysis, in the Keck era, we consider the impact of combining the orbit of 17 stars in a joint orbital-analysis. In the TMT-era, we include the orbits of all 122 stars (both currently detected and simulated).

\subsection{Cadence of observations}

In order to simulate science cases, we need to model when the observations occur, the number of observations, and the observational uncertainties. As mentioned above, we consider two classes of telescopes. The first class corresponds to existing telescopes and the epochs of observation considered consist of a combination of existing observations (between 1995 and 2017) and of future simulated observations (after 2017). From 1995-2017, we use the actual observed epochs of astrometry from Keck, and observed epochs of radial velocities from Keck and VLT \cite{2016ApJ...830...17B,2017ApJ...837...30G}. From 2017-2027, we include three astrometric and three spectroscopic observations per year.

We assume a cadence of ten astrometric and ten spectroscopic observations per year for TMT after 2027. Such a high observations frequency is required to capture the rapid motion of the fast-moving stars with short period that are expected to be discovered with the next generation of telescopes. For each year, the exact epochs of observation are randomly sampled from a uniform distribution in the observational window which spans from May 1st to September 1st.

\subsection{Observational uncertainties}\label{sec:obs_sigma}

\subsubsection{Astrometric uncertainty}
We model the astrometric uncertainty as a function of magnitude using a constant at bright magnitudes and power-law at fainter magnitudes. This is intended to simulate the observed behavior of current astrometric errors, which flatten at bright magnitudes due to systematics such as tip-tilt jitter. At faint magnitudes, the background and photon noise dominate resulting in a power-law dependence. We approximate the astrometric errors with:
\begin{equation}\label{eq:sigma_astro}
	\sigma_\textrm{astro}(m)= A \times 10^{\alpha\left(m-m_0\right)}  \quad \textrm{for} \quad m>m_0   \qquad \textrm{and} \qquad \sigma_\textrm{astro}(m)= A \quad \textrm{for} \quad m\leq m_0 \, .
\end{equation}
This relation depends on three constants: the constant $A$ which corresponds to the typical uncertainty reached for bright stars; the slope of the power-law $\alpha$ and the breaking point of the power-law $m_0$. These three parameters depends on the type of telescopes. The values used in this work are presented in Table~\ref{tab:sigma_astro}.
\begin{table}[htb]
\begin{center}
	\begin{tabular}{c c c c }
		\multicolumn{4}{c}{Astrometric uncertainty} \\
		\hline
		Telescope	&	A   [$\mu$as] &  $\alpha$	& $m_0$	 \\ \hline\hline
		Keck	    &	100	        &   0.2     &  15    \\
		TMT	        &	25          &	0.2     &  17    \\
		\hline
	\end{tabular}
\caption{Parameters characterizing the astrometric uncertainties for current and future telescopes parametrized by Eq.~(\ref{eq:sigma_astro}).}
\label{tab:sigma_astro}
\end{center}
\end{table}

\subsubsection{Radial velocity uncertainty}\label{sec:RV_sigma}
The radial velocity uncertainty depends on the spectral-type of a star (early-type or late-type) as well as its magnitude (see Section \ref{sec:errors}). We model this dependence using:
\begin{equation}\label{eq:sigma_RV}
	\sigma_\textrm{RV}(m)= B \times 10^{\beta\left(m-m_0\right)}  \quad \textrm{for} \quad m>m_0   \qquad \textrm{and} \qquad \sigma_\textrm{RV}(m)= B \quad \textrm{for} \quad m\leq m_0 \, .
\end{equation}
The constant $B$ corresponds to the typical RV uncertainty for bright stars, the parameter $\beta$ corresponds to the slope of the power-law and $m_0$ correspond to the breaking point in the power-law. The value of these three constants depend on the spectral-type and on the characteristic of telescope. They are given in Table~\ref{tab:sigma_RV}. Typically, the next generation of extremely large telescopes is expected to improve the current RV uncertainty by one order of magnitude.

\begin{table}[htb]
\begin{center}
	\begin{tabular}{cc c c }
		\multicolumn{4}{c}{Early-type star RV uncertainty} \\
		\hline
		Telescope	&	$B$   [km/s]  &  $\beta$	& $m_0$	 \\ \hline\hline
		Keck	    &	30	        &   0.17    &  14    \\
		TMT	        &	 3          &	0.17    &  14    \\
		\hline
	\end{tabular}\hspace*{2cm}
	\begin{tabular}{c c c c }
		\multicolumn{4}{c}{Late-type stars RV uncertainty} \\
		\hline
		Telescope	&	$B$   [km/s]  &  $\beta$  & $m_0$	 \\ \hline\hline
		Keck	    &	7	        &   0.37    &  14    \\
		TMT	        &	0.7         &	0.37    &  14    \\
		\hline
	\end{tabular}
\caption{Parameters characterizing the RV uncertainty for early-type stars (left table) and late-type stars (right table) measured with a Keck and TMT. The RV uncertainty depends on the star's magnitude $m$ through  Eq.~(\ref{eq:sigma_RV}).}
\label{tab:sigma_RV}
\end{center}
\end{table}

\subsection{Source confusion}\label{sec:sourceconfusion}
Confusion arises when a star is observed close (typically at an angular distance smaller than the angular resolution of the telescope) to another source. In such cases, the observations can be biased and in the worst case, faintest stars are completely hidden by bright stars and cannot be observed at all. The near-infrared emission from the black hole is also an important source of confusion for some stars at their closest approach to the black hole.

The model source confusion we identified all epochs when two (or more) orbits cross. This is done by considering the Keplerian orbits for all the 122 stars used in this work. For all these confusion events, we determine if a star can be observed or if it will be hidden by the other sources nearby. We also model an increase in uncertainty from the effects of confusion. This confusion uncertainty is then added in quadrature to the observational uncertainty in Section~\ref{sec:obs_sigma}.

We describe the model for astrometric and radial velocity confusion in more details below.

\subsubsection{Astrometric source confusion}\label{subsec_astroconfusion}
We consider that two sources are astrometrically confused at a given epoch when their angular distance is smaller than the astrometric confusion distance limit $\delta_\text{a}$. This confusion limit depends on the typical point spread function of the detector and depends on the telescope. The values used in this work are given in Table \ref{tab_telescopes}.
\begin{table}[htb]
\begin{center}
\begin{tabular}{cc}
\hline
Telescopes	&	$\delta_\text{a}$ [mas]	\\
\hline\hline
Keck	&	63 	    \\
TMT	    &   18     \\
\hline
\end{tabular}
\caption{Astrometric confusion distance limits $\delta_\text{a}$ for Keck and TMT.}
\label{tab_telescopes}
\end{center}
\end{table}

When several stars are predicted to be observed within an astrometric radius of $\delta_\textrm{a}$, we consider only the brightest star is observed and all the fainter stars will be unobservable. In addition, we model the observed astrometric position of the brightest star to be the center of light of all the stars included within a distance $\delta_\textrm{a}$. The center of light is defined as
\begin{equation}
	x^j_{CL} = \frac{F_b x_b^j + \sum_{k\neq b} w_{kb} F_k x_k^j  }{F_b + \sum_{k\neq b} w_{kb} F_k} \, ,
\end{equation}
where the subscript $j$ denotes the two different astrometric directions, the sum of all the stars $k$ that are closer than $\delta_\textrm{a}$ with the exception of the brightest star (denoted by a subscript $b$), $x^j_k$ are the astrometric positions of the star $k$, $F_k$ is its flux and $w_{kb}$ is a weighting function to model the point spread function of the detector. In this work, we use a simple model of the core of the point spread function using a Gaussian:
\begin{align}
w_{kb} = \frac{1}{\sqrt{2 \pi} \delta_\text{a}} \exp{\left( -\frac{d_{kb}^2}{2\delta_\text{a}^2} \right)} \, ,
\end{align}
where $d_{kb}$ is the astrometric separation between the star $k$ and the brightest star $b$.

We model the source confusion for the bright star by determining a confusion uncertainty which is given by the absolute value of the confusion bias defined as $x^j_b-x^j_{CL}$. The expression of the confusion uncertainty for the brightest star $b$ is therefore given by
\begin{equation}
	\sigma^j_\textrm{astro;conf} =  \frac{ \sum_{k\neq b} w_{kb} \frac{F_k}{F_b} \left|x_k^j-x_b^j\right|  }{1 + \sum_{k\neq b} w_{kb} \frac{F_k}{F_b}} \, ,
\end{equation}
where the subscript $j$ denotes the two different astrometric directions.

\subsubsection{Spectroscopic source confusion}\label{subsec_rvconfusion}
If two or more sources are close to each other on the plane of the sky with approximately the same radial velocities, measuring the RV of each star with high accuracy is difficult. We have developed a methodology that identifies when this occurs; for these epochs, there are two choices: increase the uncertainty or remove the observation from epoch for that star entirely. In order to model this we consider the spectral features of the stars, classified into three categories: featureless, early-type, or late-type. Based on a star's classification, location, and brightness, we implement a series of decisions to model the effect of confusion its radial velocity measurement:
\begin{enumerate}
	\item We do not consider radial velocity measurements for the stars that are featureless in the K-band or is fainter than the spectroscopic limit of the telescope $\text{m}_\text{f}$ (see Table~\ref{tab_rv_confusion_keck}).

	\item For each observation of a given star $s$, we examine nearby sources whose angular distance is closer than an astrometric distance threshold $\xi_\text{a}$ (the values used for this threshold are given in Table~\ref{tab_rv_confusion_keck}). If no neighbor is found, this observation does not suffer from source confusion.

	\item For all the epochs of observations when a given star has astrometric neighbor(s), we compare the magnitude of the different stars. If there is a neighboring source brighter by $\Delta m>2$, we consider the observation to be lost. On the other hand, if the star is brighter by $\Delta m>2$  than all its neighbors, we consider this measurement as unaffected by confusion. Finally, if nearby stars are within 2 magnitudes ($|\Delta m|\leq 2$) of the star $s$ under consideration, we go on with the next step.

	\item We identify nearby stars that are featureless. While these stars do not have spectral features, they add noise to the spectrum of the given star. We model the effect of this noise will affect the radial velocity uncertainty as:
	\begin{align}
	\sigma_\textrm{RV;conf;featureless} =& \frac{F_\textrm{featureless}}{F_s} \ \sigma_\textrm{s,RV} \, ,
	\end{align}
	where $F_\textrm{featureless}/F_s$ denotes the flux ratio of the featureless source to the star $s$ and $\sigma_\textrm{s,RV}$ is the observational RV uncertainty of the star described in Section~\ref{sec:RV_sigma}. This uncertainty is added  to the observational uncertainty in quadrature.

	\item We now consider the radial velocities of nearby stars that may be early or late-type stars. If the difference in radial velocities between the nearby sources and the source under consideration is larger than the threshold, $\xi_\textrm{RV}$, given in Table~\ref{tab_rv_confusion_keck}, we do not consider the radial velocity measurement to be affected by confusion.

	\item In the rare case where nearby sources have radial velocities that are not well separated from the source under consideration, we need to take into account the spectral-types of the sources involved:
	\begin{enumerate}
		\item if the star $s$ under consideration is an early-type star:
		\begin{itemize}
			\item If it is close to least one late-type star, the observation is considered missed. This is due to that the many spectral features of a late-type star will overwhelm our ability to detect the few weak lines in the early-type star.
			\item If it is close to at least one early-type star that is brighter, the observation is considered as missed, as the spectral feature of the brighter star will dominate.
			\item  If it is surrounded by early-type stars that are all fainter, we consider an observation possible, but with additional uncertainties. The fainter stars will bias the RV measurement of the observed star according to their brightness as well as their radial velocities. We model this uncertainty as:
			\begin{equation}\label{eq:RV_conf}
				\sigma_\textrm{RV;conf}=\frac{\sum_k F_k/ F_s \left|  RV_{k}-RV_s\right| }{1 + \sum_k F_k/F_\text{s}}  \, ,
			\end{equation}
			where the sum is of the nearby early-type stars $k$, $F_k$ is the flux of the stars and $RV_k$ is the radial velocity of the stars.
 		\end{itemize}
		\item If the star $s$ under consideration is a late-type star:
		\begin{itemize}
			\item If it is surrounded only by early-type stars, this observation is considered as successful and no additional confusion uncertainty is included.
			\item If it is surrounded by at least one late-type star that is brighter, the observation is considered as missed.
			\item If it is surrounded by at least one late-type star, with all the late-type stars in the neighborhood being fainter, the observation is considered as successful. In this case, the late-type stars that are astrometrically and spectroscopically confused with the star $s$ under consideration will impact its uncertainty. The confusion uncertainty is computed using Eq.~(\ref{eq:RV_conf}) where the sum $k$ is of the late-type stars near the star $s$.
		\end{itemize}
	\end{enumerate}
\end{enumerate}

\begin{table}[htb]
\begin{center}
\begin{tabular}{c  c  c  c  c  c  c  }
\hline
	& \multicolumn{3}{c} {Keck}  & \multicolumn{3}{c} {TMT} \\
stars			& $\text{m}_\text{f}$ 	&	$\xi_\text{a}$	[mas] &	$\xi_\text{rv}$	[km/s] & $\text{m}_\text{f}$ &	$\xi_\text{a}$ [mas]	& $\xi_\text{rv}$ [km/s]\\
\hline\hline
early-type 		& 16.0	&	80 	&	200 	& 22.0	&	26.7 		&	200 \\
late-type	    &	16.5	&	80	&	75.0 		&	22.5	&	26.7 		&	75 	\\
\hline
\end{tabular}
\caption{ The threshold magnitude $m_\text{f}$, the confusion distance limit $\xi_\text{a}$ and the radial velocity limit $\xi_\text{rv}$ for observation with both Keck and TMT telescopes.}
\label{tab_rv_confusion_keck}
\end{center}
\end{table}

\subsection{Sensitivity Analysis}
The different steps described previously allows us to derive a list of epochs for astrometric and spectroscopic observations and their related uncertainties for the 122 stars considered in this analysis. In this section, we will describe how we estimate the uncertainties on different scientific measurements using these observations in an orbital fit.

We re-iterate that in the Keck era, only the 17 stars that are currently well measured are considered in a joint analysis. From 2027, we start to include observations from the other 105 simulated stars.

In order to estimate the uncertainties, we compute the covariance matrix related to the linearized least-square problem and to use this covariance matrix as an estimator of the uncertainty that can be reached with a given set of observations. This covariance matrix $\Sigma$ is computed as the inverse of the normal matrix $\Sigma = N^{-1}$ defined by $N= P^T \cdot P$. The matrix $P$ contains the partial derivatives of the model with respect to the parameters. More precisely
\begin{equation}\label{eq:partial}
	P_{ij}=\frac{1}{\sigma_i}\frac{\partial M(t_i,\bm p)}{\partial p_j}\, .
\end{equation}
where $M(t,\bm p)$ represents the model to be adjusted to the observations, $\bm p$ are the parameters of the model, $t_i$ are the epochs of observations and $\sigma_i$ is the uncertainty related to the observation from the epoch $t_i$. We note that the covariance matrix obtained with this method corresponds exactly to the one that would be obtained in a Bayesian framework under the assumptions of a Gaussian likelihood, independent measurements and if the problem is sufficiently linear in the neighborhood of its minimum (in practice, if the posterior is Gaussian, the two methods will agree). Under these assumptions, the normal matrix corresponds to the Fisher information matrix. The expected uncertainty on the scientific parameters are given by the diagonal of the covariance matrix.

In this analysis, the orbital model that we use includes the Newtonian and the first post-Newtonian correction from the central SMBH and the Newtonian contribution from the extended mass. The stellar-mass density of the extended mass follows a power-law $r^{-\gamma}$, which leads to the following expression for the mass enclosed within a radius $r$
\begin{equation}
	M_\textrm{\tiny EM}(r)=M_\textrm{\tiny EM}(r_0)\left(\frac{r}{r_0}\right)^{3-\gamma} \, .
\end{equation}
The equations of motion integrated are therefore given by
\begin{equation}
	\ddot{\bm a} = - \frac{G M \bm{x}}{r^3}-G M_\textrm{EM}(r_0) \frac{r^{-\gamma}}{r_0^{3-\gamma}}\bm x +\eta \frac{G M }{c^2 r^3} \left(4\frac{G M}{r}-v^2 \right)\bm{x}  + 4 \eta\frac{G M (\bm x \cdot \bm v)}{c^2 r^3} \bm{v} \nonumber \, .
\end{equation}
The first part of this equation of motion is the standard Newtonian acceleration from the SMBH, the second part is the Newtonian acceleration due to the extended mass, the rest of this equation is the first post-Newtonian contribution from the SMBH. We included a parameter $\eta$ to allow us to parametrize a deviation from GR (i.e. in GR, $\eta=1$ whereas Newtonian gravity has $\eta = 0$). The uncertainty related to this parameter allows us to infer the signal to noise ratio (SNR) related to the detection of the GR advance of periastron ($SNR=1/\sigma_\eta$).

The result of the integration is then transformed into astrometric and RV observations by a standard Euler rotation using the Thiennes-Innes constant. The distance to the Galactic center $R_0$ is used to transform the distance in the plane of the orbit into observed angular position in the plane of the sky. We emphasize that in this work, we include 4 parameters that are related to the astrometric position and velocity of the SMBH and one parameter that is related to the velocity of the black hole along the line of sight. We allow these parameters to be free; fixing these parameters will improve the uncertainties in the GR measurement.

In total, the model used includes the following parameters: 6 orbital parameters per star (period, eccentricity, time of closest approach and the three regular angles), the SMBH $GM$, the distance to the Galactic Center $R_0$, the astrometric position and velocity of the SMBH ($x_0$, $y_0$, $v_{x_0}$ and $v_{y_0}$), the velocity of the SMBH along the line of sight ($v_{z_0}$), the amount of extended mass $M_\textrm{EM}(r_0)$ and eventually a parameter to characterize the test of GR: $\eta$.

\subsection{Results}\label{sec:results}
Using the procedure presented above, we run several simulations to measure the time evolution of the uncertainty related to several parameters of scientific interest.

The first case that we consider is a case where we assume GR to be valid (i.e. we do not consider the parameter $\eta$ in the fit). From these fits, we can see how the uncertainties on the SMBH mass, on the distance to the Galactic Center $R_0$ and on the amount of extended mass included within S0-2 apoapse evolve with time. These evolutions are presented in Figures~\ref{fig:gm_r0} and \ref{fig:gr_extm}. We consider two sets of simulations: one with no TMT observations and one where TMT observations begin in 2027. The improvement expected from TMT observations are spectacular for all parameters. More precisely, TMT observations will allow us quickly to improve the estimate on all parameters by at least one order of magnitude, reaching the impressive relative accuracy of 1 per 1000 for $GM$ and $R_0$. Regarding the extended mass, TMT allows us to probe values below $10^3$ solar masses, which will allows to distinguish between different astrophysical model\cite{2005PhR...419...65A}.

In the second simulation, we consider the parameter $\eta$ in the simulation as well. This parameter allows us to quantify the expected SNR for the relativistic advance of the periastron. This parameter can also be considered as a parameter encoding deviation from GR. As can be seen from the left panel from Figure~\ref{fig:gr_extm}, TMT will also improve the detection of relativistic effects by one order of magnitude  and allows us to further improve tests of GR in a strong-gravity regime around a SMBH\cite{2017PhRvL.118u1101H}.

In addition to the scientific cases presented here, TMT will allow us to measure the SMBH spin. A combined measurement of the SMBH spin\cite{2017arXiv170904492G} and quadrupole moment would also lead to a test of the no-hair conjecture as proposed by \citenum{2008ApJ...674L..25W}. Simulations including these effects are currently on-going.

\begin{figure}[hbt]
\center
\includegraphics[width=6.0in]{fig5.pdf}
\caption{The left panel shows the uncertainty on the black hole mass parameter in the unit of mass of the sun. The left panel displays the uncertainty on the distance of SgrA* in parsec. The red curve shows the evolution with Keck and TMT data. The blue curve on the other hand is with only Keck data.}
\label{fig:gm_r0}
\end{figure}

\begin{figure}[hbt]
\center
\includegraphics[width=6.0in]{fig6.pdf}
\caption{The left panel shows the signal to noise ratio for the GR precession signal and the left panel displays the uncertainty on the extended mass parameter at the Galactic Center. The red curve shows the evolution with Keck and TMT data. The blue curve on the other hand is with only Keck data. The grey dashed horizontal line in the right panel is the extended mass value at the distance of semi--major axis of S0--2.}
\label{fig:gr_extm}
\end{figure}

\section{CONCLUSIONS}
\label{sec:conclusions}

The measurement of stellar orbits at the Galactic center holds the potential to push the frontiers of our understanding of the physics and astrophysics of supermassive black holes with the advent of future extremely large ground-based optical and near-infrared telescopes. While the potential tests of General Relativity and other theories of gravity at the Galactic center have been anticipated for years, this is the first study to quantitatively study the impact on different science cases using realistic models of the telescope, instrument, observing cadence, and simulations of the underly stellar cluster. We use a model of the TMT AO system, the IRIS instrument, and our current knowledge of the orbits and stellar population in our simulations. We find that in all cases, within a few years of TMT operation, there will be a large leap in the sensitivity of the orbits to physical effects such as the GR precession of the orbit and the extended mass distribution. The measurement of the precession of the periapse from General Relativity will be measured at $> 10\sigma$ significance, while fundamental parameters such as the distance to the Galactic center can be measured to better than 5 pc, which is 20 times better precision than today. The science cases outlined here are only a small sample of the potential science that is anticipated. There are also many areas that will only be discovered once this region is revealed at with unprecedented sensitivity and angular resolution.


\section{Introduction}
The realization of topological aspects of the band theory of solids has opened new avenues in the study of insulators. A particularly dramatic manifestation of the topology is the presence of \emph{boundary modes}, i.e, single particle states exponentially localized at surface/edge of a topologically nontrivial insulator with energies lying in the bulk gap. These modes are protected ``by the topology'' and cannot be removed by local perturbations\cite{bhbook, shen_book, hasan-kane_TI, kane-mele, fu-kane-mele}. The topological invariants associated with the boundary modes have been studied in diverse contexts\cite{tknn, bellissard-baldes_noncomm_IQHE, baldes_disorder_IQHE, prodan_disorder_TI, volovik_book, essin-gurarie_bulk_bdry}, and related to the topological invariants of the bulk under the \emph{bulk-boundary correspondence}.

Many of these ideas have also been generalized to semimetals, which are gapped almost everywhere in the bulk Brillouin zone(BZ). The intuition for insulators can thus often be carried over, albeit with momenta defined only on a punctured BZ. A notable example are the 3+1 dimensional \emph{Weyl semimetals}(WSMs)\cite{turner-ashvin2011, witten_wsm}, gapless at a finite (even) number of isolated, topologically protected \emph{Weyl nodes} in the bulk BZ. The surface exhibits \emph{Fermi arcs}, curves in the surface BZ of topologically protected zero energy modes connecting the projections of the Weyl nodes. Recently, WSMs have attracted much attention recently owing to various experimental realizations\cite{hasan_wsm, hasan_wsm2, hasan_wsm3, marin_wsm_photonic,marin_wsm_photonic2}.

WSMs are particularly interesting because the \emph{nodal configuration}, i.e, the positions and chiralities of Weyl nodes (or equivalently, the linearized theory of the bulk) does not completely specify the surface spectrum for a fixed boundary condition. Instead, one also needs to know the \emph{Fermi arc connectivity}, i.e, the pairs of nodes that are connected by Fermi arcs on a given surface, which may depend on the microscopic details of the bulk\cite{vd-str_wsm}. Mathematically, the connectivity can be described in terms of a relative homology group\cite{mathai_wsm_top} or \emph{Euler structures}\cite{mathai_semimetal_top}, which encode the additional global information required to completely characterize a WSM.

The arguments for existence of boundary modes in topological phases can be readily generalized to \emph{interface modes} localized at junctions between topological phases, where one would expect localized modes whenever the two phases forming the junction are topologically distinct\footnote{
  Indeed, the boundary modes for topological phases can be thought of as a special case corresponding to one of these phases being the vacuum, a topologically trivial insulator with infinite energy gap.
}.
Since the formation of a junction can be intuitively visualized as the coupling (due to proximity) of the surface modes of the two topological phases involved, the interface modes can be expected to inherit aspects of topology from both phases. Various such cases involving topological insulators (TIs) and superconductors (TSCs) have already been studied, notable examples being \mbox{TI-TI} junctions\cite{takahashi_ti_opp_vel, beule_ti_junc}, \mbox{TI-WSM} junctions\cite{adolfo_ti_wsm,juergens_ti_wsm} and \mbox{WSM-TSC} junctions\cite{sumathi_wsm_sc, chen_wsm_sc}.

In this article, we consider the interface modes for a WSM-WSM junction. For two WSMs with different nodal configurations, one generically expects zero modes corresponding to the ``difference'' between the Fermi arcs on the two sides. However, a more intriguing possibility is to consider two WSMs with identical nodal configurations, but different Fermi arc connectivities, as schematically illustrated in \fig{fig:sch}. These Fermi arcs cannot be coupled and hence gapped out by hopping terms between surface modes that maintain the translation invariance in the transverse directions, since they lie on different transverse momenta(in the surface BZ).

We find that the coupling between the Fermi arcs of the two surfaces leads to a \emph{Fermi arc reconstruction}. In the tight binding models considered here, the Fermi arcs generically combine to form closed loops of zero modes localized at the interface. In certain cases, these loops detach from the Weyl nodes, leading to a stable Fermi surface of effectively 2D modes. Given tight-binding models for WSMs, these modes can be obtained by a numerical exact diagonalization(ED) of the real space Hamiltonian.

To study the Fermi arc reconstruction analytically, we generalize the transfer matrix method to include interfaces. This approach, being partly a real space method, is ideally suited to probing spatially inhomogeneous systems.
It has also been used to analytically study the bulk and boundary states of noninteracting toplogical phases\cite{hatsugai_cbs, hatsugai_cbs_PRL, dh_lee_cbs, schulman_cbs, tauber-delplace_tm, vd-vc_tm}. Building up on this literature, we derive a closed form condition for the existence of interface modes in terms of the transfer matrices for systems on the two sides of the interface, which is applicable to interfaces between systems described by tight-binding models.

The rest of this article is organized as follows: In Sec.~\ref{sec:wsm}, a generic two band tight-binding model for WSMs is introduced and the associated transfer matrix is computed. In Sec.~\ref{sec:int}, general conditions for the presence of interface modes are derived, which are then applied to specific WSM models in Sec.~\ref{sec:eg}. We conclude with a discussion in Sec.~\ref{sec:conc}, and some of the tedious calculations are relegated to the appendices.

\begin{figure}
	\centering
	\includegraphics[width=0.85\columnwidth]{schematic.pdf}
	\caption{
	  A schematic depiction of an interface between WSMs with different Fermi arc connectivities. The hopping between the two surfaces, denoted by a gray arrow, leads to hybridization of the localized modes on the two surfaces, i.e, a \emph{Fermi arc reconstruction}.
	}
	\label{fig:sch}
\end{figure}

\section{Weyl semimetals}     \label{sec:wsm}
In this section, we study a general family of lattice models which can realize a WSM with an arbitrary (even) number of Weyl nodes lying in a plane in the bulk BZ. We construct the transfer matrix for translations normal to this plane and use it to study the surface modes.

\subsection{The lattice model}   \label{sec:wsm_model}
Consider the family of 2-band models described by the Bloch Hamiltonian
\begin{equation}
  \mathcal{H}(\mathbf{k}) = \mathcal{H}_x(k_x) + \mathcal{H}_\perp(\mathbf{k}_\perp),   \label{eq:hlt_orig}
\end{equation}
where we take
\begin{align} \label{eq:4nodewsm1}
  \mathcal{H}_x(k_x) = & \; \sin k_x \sigma^x + (1 - \cos k_x) \sigma^z, \nonumber \\
  \mathcal{H}_\perp(\mathbf{k}_\perp) = & \; \eta_y(\mathbf{k}_\perp) \sigma^y + \eta_z(\mathbf{k}_\perp) \sigma^z.
\end{align}
Here, ${\boldsymbol{\eta}} = \left( \eta_y, \eta_z \right)^T : \mathbb{T}^2 \to \mathbb{R}^2$ are functions of the transverse momentum $\mathbf{k}_\perp \equiv (k_y, k_z) \in \mathbb{T}^2$, and the Pauli matrices describe the internal (\emph{pseudospin}) degrees of freedom. The spectrum is gapless in the $k_x = 0$ plane for
\begin{equation}
  \eta_y(\mathbf{k}_\perp^\ast) = \eta_z(\mathbf{k}_\perp^\ast) =  0.   \label{eq:nodes}
\end{equation}
We preclude the possibility of gapless modes in the $k_x = \pi$ plane by demanding that $\eta_z > -2$.

The Hamiltonian can be linearized near a particular gapless point $\mathbf{k}^\ast = (0, \mathbf{k}_\perp^\ast)$ as
\[
 \mathcal{H} \approx \delta k_a V_{ab} \sigma^b, \quad a,b \in \{x, y,z\},
\]
where
\begin{equation}
  V = \begin{pmatrix}
	1 & 0 \\
	0 & \mathbb{J}
      \end{pmatrix};
      \quad \mathbb{J}(\mathbf{k}_\perp^\ast) = \left. \frac{\partial{\boldsymbol{\eta}}}{\partial\mathbf{k}_\perp} \right|_{\mathbf{k}_\perp^\ast}
\end{equation}
being the Jacobian of ${\boldsymbol{\eta}}$. The system exhibits a Weyl node at $\mathbf{k}^\ast$ iff $V$ (and hence $\mathbb{J}$) is nonsingular, in which case the associated chirality is given by
\begin{equation}
  \chi = \sgn{ \det{V(\mathbf{k}^\ast)} } = \sgn{ \det{\mathbb{J}(\mathbf{k}_\perp^\ast)} }.
\end{equation}

Choosing a suitable ${\boldsymbol{\eta}}$, WSMs with arbitrary (even) number of Weyl nodes can be realized. We expect Fermi arcs on a surface normal to $x$, with the surface BZ parametrized by $\mathbf{k}_\perp$ and the surface projections of the Weyl nodes simply the solutions of \eq{eq:nodes}. To study the surface spectrum analytically, we next construct the transfer matrix for translations along the $x$-direction.

\subsection{Constructing the transfer matrix}    \label{sec:wsm_tm}
Consider the model of \eq{eq:hlt_orig} on a cylinder $[0,N] \times \mathbb{T}^2$, finite along $x$. Since $k_x$ is not a good quantum number anymore, inverse Fourier transform the Hamiltonian of \eq{eq:hlt_orig} along $x$ to get
\begin{align}
  \mathcal{H}_{\text{real}}(\mathbf{k}_\perp) = & \; \sum_{n}  \bigg[  \mathbf{c}^\dagger_{n+1} \left( \frac{ \sigma^x - i \sigma^z}{2} \right) \mathbf{c}_n + \text{H.c} \nonumber \\
  & + \mathbf{c}^\dagger_n \left( \eta_y \sigma^y + (1 + \eta_z) \sigma^z \right) \mathbf{c}_n \bigg],    \label{eq:hlt_real}
\end{align}
where $\mathbf{c}_n (\mathbf{k}_\perp) = (\c_{n,1}, \c_{n,2})^T$ are the annihilation operators for the two pseudospin states at the layer $n$ along $x$.
The Schr\"odinger equation, $\mathcal{H} \ket{\Psi} = \mathbf{e} \ket{\Psi}$, can be written as a recursion relation\cite{hatsugai_cbs, vd-vc_tm}
\begin{equation}
  J \Psi_{n+1} + M \Psi_n + J^{\dagger} \Psi_{n-1} = \mathbf{e} \Psi_n,    \label{eq:recur}
\end{equation}
where $\Psi_{n} \in \mathbb{C}^2$
and
\begin{align}
  J = \frac{1}{2i} \left( \sigma^x - i \sigma^z \right), \quad   M = \eta_y \sigma^y + (1 + \eta_z) \sigma^z
\end{align}
are the \emph{hopping} and \emph{on-site} matrices, respectively.

The conventional transfer matrix construction works by rewriting \eq{eq:recur} as an equation for $\Psi_{n+1}$ in terms of $\Psi_{n}$ and $\Psi_{n-1}$ by inverting $J$. However, since $J$ is singular ($\text{rank}\, J = 1$), we follow the alternative construction introduced in Ref~\onlinecite{vd-vc_tm}, where the recursion relation is instead rewritten as
\begin{equation}
  \Psi_n = \mathcal{G} \left( J \Psi_{n+1} + J^{\dagger} \Psi_{n-1} \right),   \label{eq:recur_G}
\end{equation}
where
\begin{align}
  \mathcal{G} \equiv & \;  (\mathbf{e} \mathbbm{1} - M)^{-1} \nonumber \\
  = & \; \frac{1}{\mathbf{e}^2 - \Gamma^2} \left[ \mathbf{e} \mathbbm{1} + \eta_y \sigma^y + (1 + \eta_z) \sigma^z  \right],
\end{align}
can be identified as an \emph{on-site Green's function}, with
\begin{equation}
  \Gamma^2 = \eta_y^2 + (1 + \eta_z)^2.  \label{eq:gamma_def}
\end{equation}
Next, we take the reduced singular value decomposition of $J$, \emph{viz}, $J = \mathbf{v} \cdot \mathbf{w}^{\dagger}$, where
\begin{equation}
  \mathbf{v} = \frac{1}{\sqrt{2}} \vecenv{-i}{1}, \quad
  \mathbf{w} =  \frac{1}{\sqrt{2}} \vecenv{i}{1},
\end{equation}
which satisfy
\[
  \mathbf{v}^{\dagger} \mathbf{v} = \mathbf{w}^{\dagger} \mathbf{w} = 1, \quad \mathbf{v}^{\dagger} \mathbf{w} = 0,
\]
since $J^2 = 0$. Thus, $\{\mathbf{v}, \mathbf{w}\}$ forms an orthonormal basis of $\mathbb{C}^2 \ni \Psi_n$, and we may expand $\Psi_n$ as
\begin{align}
  \Psi_n = & \; \left( \mathbf{v}^{\dagger} \Psi_n \right) \mathbf{v} + \left( \mathbf{w}^{\dagger} \Psi_n \right) \mathbf{w}  \nonumber \\
  \equiv & \;  \alpha_n \mathbf{v} + \beta_n \mathbf{w},
\end{align}
The recursion relation of \eq{eq:recur_G} thus becomes
\begin{equation}
  \Psi_n = \mathcal{G} \mathbf{v} \beta_{n+1} + \mathcal{G} \mathbf{w} \alpha_{n-1},
\end{equation}
so that the coefficients in the $\{\mathbf{v}, \mathbf{w}\}$ basis are
\begin{align}
 \alpha_n = & \; \mathbf{v}^{\dagger} \mathcal{G} \mathbf{v} \, \beta_{n+1} + \mathbf{v}^{\dagger} \mathcal{G} \mathbf{w} \, \alpha_{n-1},  \nonumber \\
 \beta_n = & \; \mathbf{w}^{\dagger} \mathcal{G} \mathbf{v} \, \beta_{n+1} + \mathbf{w}^{\dagger} \mathcal{G} \mathbf{w} \, \alpha_{n-1},
\end{align}
where
\begin{align}
 & \mathbf{v}^{\dagger} \mathcal{G} \mathbf{v} = \frac{\mathbf{e}  + \eta_y}{\mathbf{e}^2 - \Gamma^2}, \quad \mathbf{w}^{\dagger} \mathcal{G} \mathbf{w} = \frac{\mathbf{e} - \eta_y}{\mathbf{e}^2 - \Gamma^2}, \nonumber \\
 & \mathbf{v}^{\dagger} \mathcal{G} \mathbf{w} = \mathbf{w}^{\dagger} \mathcal{G} \mathbf{v} = - \frac{1 + \eta_z}{\mathbf{e}^2 - \Gamma^2}.
\end{align}
These equations can finally be rearranged to get a transfer matrix equation:
\begin{equation}
  \Phi_{n+1} = T \Phi_n, \quad \Phi_n \equiv \vecenv{\beta_n}{\alpha_{n-1}},
\end{equation}
where the transfer matrix is
\begin{equation}
  T(\mathbf{e}, \mathbf{k}_\perp) = \frac{1}{1 + \eta_z}\left(
  \begin{array}{cc}
    \mathbf{e}^2 - \Gamma^2 & \quad  -(\mathbf{e} - \eta_y) \\
    \mathbf{e} + \eta_y & \quad -1
  \end{array} \right).    \label{eq:tmat}
\end{equation}
Note that the vectors $\Phi_n$ that are actually translated by the transfer matrix consist of the physical wavefunctions at two adjacent sites.

\begin{figure}[b!]
	\centering
	\includegraphics[width=0.65\columnwidth]{surface_modes.png}
	\caption{
	  The geometric picture of the Weyl nodes and Fermi arcs for the tight-binding model of \eq{eq:hlt_orig} (with ${\boldsymbol{\eta}}$ defined by \eq{eq:4nodewsm} for $\varphi = 4\pi/3$). The Weyl nodes occur in the the $k_x = 0$ plane at the intersection of the curves defined by $\eta_y = 0$ and $\eta_z = 0$, depicted by solid red lines and dashed blue lines, respectively. The surface has localized modes for $\mathbf{k}_\perp = (k_y, k_z)$ in the shaded region, and the subset of these modes at zero energy, viz, the Fermi arcs, are denoted by the bold red lines.
	}
	\label{fig:wsm_engg_sch}
\end{figure}

\subsection{Surface states}    \label{sec:wsm_surf}
The spectrum of the transfer matrix (for a given $\mathbf{e}$ and $\mathbf{k}_\perp$) encodes the bulk/edge states of the system. More explicitly, note that the transfer matrix can be used to write wavefunctions as $\Phi_n = T^n \Phi_0$, so that knowing $\Phi_0$ completely describes a wavefunction. It follows that given an eigenvector,
\begin{equation}
  T \Phi_0 = \rho \Phi_0 \implies \Phi_n = \rho^n \Phi_0,
\end{equation}
so that $\Phi$ describes a bulk (Bloch) state if $\abs{\rho}=1$ and an edge state exponentially localized at the left(right) edge for $\abs{\rho}<1$ ($\abs{\rho}>1$). Furthermore, the boundary conditions are used to obtain the allowed $\Phi_0$'s for the physical eigenstates of the Hamiltonian.

For surface states, one typically considers the Dirichlet boundary condition ($\Psi_{-1} = \Psi_{N+1} = 0$), which translates to demanding that
\begin{equation}
  \Phi_0 \propto \vecenv{1}{0}, \quad
  \Phi_{N+1} \propto \vecenv{0}{1}, \quad
\end{equation}
at the left and the right edge, respectively. For the left edge, we explicitly compute
\begin{equation}
  T \Phi_0 = \rho \Phi_0 \implies
  \begin{cases}
    \mathbf{e} = -\eta_y & \\
    \rho = -(1 + \eta_z). &
  \end{cases}
\end{equation}
Thus, for a given $\mathbf{k}_\perp$, there exists a surface mode with energy $\mathbf{e} = -\eta_y(\mathbf{k}_\perp)$ exponentially localized\footnote{
  We have assumed a semi-infinite system.
} at the left surface if $\abs{\rho}<1$, i.e, if
\begin{equation}
  -1 \leq 1 + \eta_z(\mathbf{k}_\perp) \leq 1 \implies \eta_z(\mathbf{k}_\perp) \leq 0,
\end{equation}
since we have assumed that $\eta_z > -2$. Finally, for the left surface, the Fermi arc, i.e, the locus of the surface modes with $\mathbf{e} = 0$, are given by
\begin{equation}
  \eta_y(\mathbf{k}_\perp) = 0, \quad   \eta_z(\mathbf{k}_\perp) < 0,    \label{eq:fermi_arc}
\end{equation}
with the surface projection of the Weyl nodes given by the end points $\eta_y(\mathbf{k}_\perp) = 0 = \eta_z(\mathbf{k}_\perp)$. An identical calculation for the right surface (i.e, using $\Phi_{N+1}$) yields surface modes dispersing in the opposite direction, i.e, $\mathbf{e} = \eta_y$. The corresponding Fermi arcs are given by
\begin{equation}
  \eta_y(\mathbf{k}_\perp) = 0, \quad   \eta_z(\mathbf{k}_\perp) < 0,    \label{eq:fermi_arcR}
\end{equation}
which is identical to those on the left surface.

Eqns \ref{eq:nodes} and \ref{eq:fermi_arc} have a straightforward geometric interpretation, which provides useful intuition for the construction of WSM lattice models. The level sets $\eta_y = 0$ and $\eta_z = 0$ generically define two sets of curves on $\mathbb{T}^2$, whose intersections are the Weyl nodes (in the $k_x = 0$ plane). The Fermi arcs then lie along the former ($\eta_y = 0$) curves, stretching between the intersections with the latter curve. A particular example is shown in \fig{fig:wsm_engg_sch}.

\section{Interfaces}    \label{sec:int}
In this section, we consider interfaces between noninteracting topological phases described by tight binding models, and derive a general analytical condition for the existence of a localized mode at the interface in terms of the associated transfer matrices. We apply these conditions to the transfer matrix derived in \sect{sec:wsm} for simple interfaces where the analytical solutions are particularly illuminating.

\subsection{A qualitative picture}  \label{sec:int_qual}
Before we get into explicit computations, we seek to qualitatively understand the circumstances under which one might get modes localized at an interface between two topological phases. To arrive at an intuitive picture, consider the two phases being completely decoupled, with their individual surface modes with energies in the corresponding bulk gaps. As one turns on a hopping between the two surfaces (hereafter referred to as \emph{gluing}), these modes may gap out in pairs. Since the interface does not break the translation symmetry in the transverse directions, the corresponding momenta are conserved, so that only hoppings between modes with identical $\mathbf{k}_\perp$ are allowed.

Explicitly, let there be localized modes with energies $\mathbf{e}_{1,2} (\mathbf{k}_\perp)$ on the two surfaces before they are glued. Then, the most general effective Hamiltonian for the interface can be written as
\begin{equation}
  \mathcal{H}_{\text{eff}}(\mathbf{k}_\perp) =
  \begin{pmatrix}
    \mathbf{e}_1(\mathbf{k}_\perp) & t(\mathbf{k}_\perp) \\
    t^\ast(\mathbf{k}_\perp) & \mathbf{e}_2(\mathbf{k}_\perp),
  \end{pmatrix}
\end{equation}
whose spectrum is
\[
  \mathbf{e} = \frac{1}{2} \left( \mathbf{e}_1 + \mathbf{e}_2 \pm \sqrt{(\mathbf{e}_1 - \mathbf{e}_2)^2 + 4 \abs{t}^2} \right).
\]
The locus of zero energy modes is then given by
\begin{equation}
  (\mathbf{e}_1 + \mathbf{e}_2)^2 -(\mathbf{e}_1 - \mathbf{e}_2)^2 =  4 \abs{t}^2.       \label{eq:glued_zero_modes}
\end{equation}
Given $t(\mathbf{k}_\perp)$, this equation describes the zero energy modes localized at the interface.

One can immediately consider one trivial case, \emph{viz}, when $\mathbf{e}_1 = -\mathbf{e}_2$, for which there are no zero modes for any choice of $t$. This is the case, for instance, if the two systems being glued are identical, since the left and right surface states being glued have surface states that disperse in the opposite directions. This is to be expected, since after gluing one is left with a single uniform system, which does not have any localized modes in the bulk.

While good to obtain some intuition, this picture is usually not complete, since the relation between $t(\mathbf{k}_\perp)$ and the actual boundary conditions at the interface are not obvious. Furthermore, there is no direct way to study the associated wavefunction, to compute, for instance, the penetration depth, without resorting to a real space ED calculation. Therefore, in the next subsection, we shall instead use the transfer matrix approach to study the interface modes.

In the case of WSMs, one is usually interested in the zero energy surface modes, \emph{viz}, Fermi arcs. At low energies, the corresponding linearized spectrum is that of a chiral fermion dispersing in the direction normal to the arc\footnote{
  Recall that since a WSM can be thought of as stacked Chern insulators, each point of the Fermi arc is the zero crossing of the edge spectrum of a Chern insulator.
}.
Let the two surfaces have Fermi arcs with normal vectors(lying in the surface BZ) $\mathbf{n}_{1,2}$ near some $\mathbf{k}_\perp$, so that in a neighborhood, the associated linearized spectra are $\mathbf{e}_i = \delta \mathbf{k}_\perp \cdot \mathbf{n}_i$. Using \eq{eq:glued_zero_modes}, the zero modes after gluing are given by
\begin{equation}
  \left[ \delta\mathbf{k}_\perp \cdot (\mathbf{n}_1 + \mathbf{n}_2) \right]^2 - \left[ \delta\mathbf{k}_\perp \cdot (\mathbf{n}_1 - \mathbf{n}_2) \right]^2 = 4 \abs{t}^2.
\end{equation}
Generically, this is an equation for a hyperbola, which can be interpreted as a \emph{level repulsion}, albeit in the $\mathbf{k}_\perp$ space. One extreme case is to have $\mathbf{n}_1 = \mathbf{n}_2$, corresponding to parallel Fermi arcs running in the same directions. The interface modes are then given by  $\left( \delta\mathbf{k}_\perp \cdot \mathbf{n}_1 \right)^2 = \abs{t}^2$, which denotes two lines separated by $2 \abs{t}$. We shall revisit this case in \sect{sec:eg_2}.

Finally, we conclude with a conundrum for WSMs: Consider an interface between two WSMs with identical nodal configuration, for which surface modes exist in disjoint regions of the surface BZ. In general, one would not expect modes at the interface, since the two WSMs are topologically identical. However, the gluing mechanism described above precludes gapping out the original Fermi arcs, since there are no localized modes at the same $\mathbf{k}_\perp$ to couple to.  We study this case in more detail in \sect{sec:eg_4}, and discover that there are, in fact, modes localized at the interface.

\begin{figure*}
	\centering
	\includegraphics[ width=0.95\textwidth]{int.pdf}
	\caption{
	  The basic setup for describing the interface in the transfer matrix formalism. The left and right bulk layers are depicted by green and yellow, while the interface layers are depicted in gray. The matrices $J$'s and $M$'s are generically functions of the transverse momentum.
	}
	\label{fig:int}
\end{figure*}

\subsection{The general condition}     \label{sec:int_gen}
The transfer matrix approach, being inherently a real space method, is ideally suited for the analytical study of interfaces. In particular, the transfer matrix description describes the boundary modes as the solution of two conditions, \emph{viz}, the \emph{decay condition}, associated with the eigenvalues of the transfer matrix and the \emph{boundary condition} corresponding to the associated eigenvectors. The physical states are naturally required to satisfy both, but choosing either one separately, one may often extract additional topological information about the system\cite{vd-vc_tm}.

We now adapt this approach to study the interface modes. As the interface naturally breaks the translation invariance in the $x$ direction, consider the transfer matrices translating along $x$ for the phases to the left ($T_L$) and right ($T_R$) of the interface. Furthermore, there can generally be additional on-site terms at the interface, so we assume that the interface is $W\geq 0$ layers wide. Setting $n=0$ to be its left end, we can use the transfer matrices to write the wavefunctions for $n \notin [1,W-1]$ as
\begin{equation}
  \Phi_n = \begin{cases}
             T_R^{n-W} \Phi_W, & n \geq W,  \\
             T_L^{n} \Phi_0, & n \leq 0.    \label{eq:tmatLR}
           \end{cases}
\end{equation}
Thus, a mode localized at the interface must satisfy
\begin{align}
 & T_R \Phi_W = \rho_R \Phi_W,  & \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \abs{\rho_R} < 1,  &  \quad \nonumber \\
 & T_L \Phi_0 = \rho_L \Phi_0,  & \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \abs{\rho_L} > 1.  &  \quad   \label{eq:rhoLR_cond}
\end{align}
Furthermore, let the effect of the interface on a wavefunction be encoded in an invertible matrix $S$, so that
\begin{equation}
 \Phi_W = S \Phi_0.   \label{eq:S_def}
\end{equation}
Physically, this matrix can be thought of as the transmission matrix for a scattering problem with the interface being the scatterer. This is inherently a real space quantity, and depends on the geometry of the interface and the on-site potential(s). In the next subsection, we shall compute it explicitly for various cases.

The conditions for the existence of a localized interface mode can now be separated into a \emph{matching} and a \emph{decay} condition. This is a discrete analogue of the strategy for computing the bound states of a 1D potential well in quantum mechanics. To wit, one obtains general solutions (with certain unknown constants) to the left and the right of the potential well which exhibit the desired asymptotic behavior at $\pm \infty$, and then fixes the constants by \emph{matching} the solutions across the potential well. We shall follow a similar strategy here, however, we shall find unlike the case of potential wells, the matching conditions are easier to compute in terms of the transfer matrices.

We note that using Eq.~\ref{eq:tmatLR} and \ref{eq:S_def}, a wavefunction of the system is completely determined by the knowledge of $\Phi_0$. Thus, we shall seek to derive a condition only in terms of $\Phi_0$. In particular, \eq{eq:rhoLR_cond} reduces to
\begin{equation}
  S^{-1} T_R S \Phi_0 = \rho_R \Phi_0, \quad T_L \Phi_0 = \rho_L \Phi_0,
\end{equation}
alongwith the condition on $\rho_{L,R}$.

For the matching condition, we ignore the constraints on $\rho_{L,R}$ and simply demand that $\Phi_0 \neq 0$ be a simultaneous eigenvector of $S^{-1} T_R S$ and $T_L$, which implies that
\begin{equation}
  [T_L, S^{-1} T_R S]\Phi_0 = 0.    \label{eq:cond_comm}
\end{equation}
As it turns out\cite{shemesh_common_eigs}, for a $2\times 2$ transfer matrix\footnote{
  For larger transfer matrices, the latter condition is necessary, but not sufficient for a shared eigenvector. A more general condition for arbitrarily-sized transfer matrices is derived in Appendix \ref{app:cond}.
}, this is a necessary as well as \emph{sufficient} condition for the existence of a simultaneous eigenvector. Finally, \eq{eq:cond_comm} has a nontrivial solution for $\Phi_0$ iff the commutator satisfies the Cramers' condition:
\begin{equation}
  \det\left( [T_L, S^{-1} T_R S] \right) = 0.   \label{eq:cond_det}
\end{equation}
The LHS is a function of $\mathbf{e}$ and $\mathbf{k}_\perp$, so that we get an implicit equation for the points that satisfy the matching condition.

Given a particular $(\mathbf{e}, \mathbf{k}_\perp)$ satisfying \eq{eq:cond_det}, one can then compute $\Phi_0 \in \ker \left( [T_L, S^{-1} T_R S] \right)$, so that $\Phi_0$ is an eigenvector of $T_L$ while $S \Phi_0$ is an eigenvector of $T_R$. For the corresponding eigenvalues, we have the following three possibilities:
\begin{enumerate}
 \item $\abs{\rho_R} < 1, \, \abs{\rho_L} > 1$.
 \item $\abs{\rho_R} > 1, \, \abs{\rho_L} < 1$.
 \item $\abs{\rho_A} \lessgtr 1$ for both $A = L,R$.
 \label{blah}
\end{enumerate}
The first case corresponds to a mode decaying into the bulk on both sides of the interface, i.e, the (\emph{physical}) decay condition. The remaining cases are also relevant, since they can be thought of as the physical interface modes for a different choice of the gluing matrix $S$. For instance, the second case, with a mode growing in the bulk on both sides, can be obtained by swapping the two phases ($L \leftrightarrow R$) and replacing the gluing matrix $S$ by $S^{-1}$. The relevance of these ``unphysical'' cases is discussed in more detail in Appendix \ref{app:unphys}.

Finally, we note that the eigenvalues $\rho_{L,R}$ must vary continuously with $\mathbf{k}_\perp$ and cross the unit circle only in pairs\cite{vd-vc_tm}, these labels cannot change unless the curves intersect\footnote{
  Parts of these ``curves'' can actually be single points. This is the case, for instance, for the projections of the Weyl nodes on the interface BZ.
}. We can then label each segment of the curves obtained from the matching conditions by the eigenvalues $\rho_{L,R}$. Thus, in practice, given a pair of phases described by $T_L$ and $T_R$ and a gluing matrix $S$, we shall compute the curves given by the matching condition (\eq{eq:cond_det}) in the interface BZ. For each section of the curve,  we then compute the null eigenvector of $[T_L, S^{-1}T_R S]$ to assign it one of the three possibilities discussed above.

\subsection{Computations for specific interfaces}     \label{sec:int_spc}
The conditions for interface modes derived in the previous subsection do not refer to the specific form of the transfer matrices or interfaces. We now consider the specific case of the transfer matrix derived in \sect{sec:wsm}. We begin by recalling that for a translation invariant system, the Schr\"odinger equation can be written as a recursion relation\cite{vd-vc_tm}, as shown in \eq{eq:recur}. For spatially inhomogeneous case, this generalizes to
\begin{equation}
  J_n \Psi_{n+1} + M_n \Psi_n + J_{n-1}^\dagger \Psi_{n-1} = \mathbf{e} \Psi_n,       \label{eq:schr_pos_dep}
\end{equation}
where both the hopping and on-site matrices are position dependent. Consider the setup depicted in \fig{fig:int}, which explicitly corresponds to setting
\begin{align}
  J_n = & \; \begin{cases}
	  J_L, & n \leq -1, \\
	  \wt{J}_n, & 0 \leq n \leq W-2, \\
	  J_R, & n \geq W-1,
        \end{cases}
  \nonumber \\
  M_n = & \; \begin{cases}
	  M_L, & n \leq -1, \\
	  \wt{M}_n, & 0 \leq n \leq W-1, \\
	  M_R, & n \geq W,
        \end{cases}
  \label{eq:}
\end{align}
where $\wt{J}_{n}$ and $\wt{M}_{n}$ encode the interface of width $W \geq 0$. The case with $W=0$ is not well defined unless $J_R = J_L$. So far, this setup is valid for interfaces between any systems described by tight-binding models.

In the rest of this article, we restrict ourselves to interfaces between two WSMs described by \eq{eq:hlt_orig} with ${\boldsymbol{\eta}} = {\boldsymbol{\eta}}^{L/R}$ to the left/right of the interface, so that $J_L = J_R = J$ and
\begin{equation}
  M_A = \eta_y^A \sigma^y + (1 + \eta_z^A) \sigma^z, \quad A = L, R.
\end{equation}
The transfer matrices $T_{L,R}$ for the two sides were computed in \sect{sec:wsm_tm}. We now seek to apply the results of \sect{sec:int_gen} to study the zero modes, if any, localized at the interfaces. The simplest case of an interface is is an abrupt one, with no additional local potential. We shall keep the coupling between the two sides to be tunable, so that the hopping matrix connecting the two sides is simply $\kappa J$, with $\kappa \in [0,\infty)$, the two limits corresponding to the uncoupled and the strongly coupled case. A quick glance at \fig{fig:int} reveals that for $\wt{J}_0 = \kappa J$, we must take $W = 2$; however, for $\kappa = 1$, we may take $W = 0$, rendering it particularly simple to analyze. This is the case that we start with.

\begin{figure}[t]
	\centering
	\includegraphics[ width=0.7\columnwidth]{int0.pdf} \\ \ \\
	\includegraphics[ width=0.8\columnwidth]{int1.pdf}
	\caption{
	  Schematic of the interfaces discussed in \sect{sec:int_spc}: (Top) A hard interface with a tunable coupling $\kappa$ between the surfaces. (Bottom) A softer interface with an on-site potential described by $\wt{M}_0$.
	}
	\label{fig:int_eg}
\end{figure}

\subsubsection{Hard interface for $\kappa = 1$}
We set $W = 0$, so that
\[
  J_n = J, \quad
  M_n = \begin{cases}
	  M_L, & n \leq -1, \\
	  M_R, & n \geq 0.
        \end{cases}
\]
The matrix $S$ is required to satisfy $S\Phi_0 = \Phi_W$, so that we need $S = \mathbbm{1}$. The matching condition for the interface then becomes
\begin{equation}
  \det\left( [T_L, T_R] \right) = 0.
\end{equation}
Using the transfer matrices $T_{L,R}$ computed in \eq{eq:tmat}, we show in Appendix \ref{app:calc} that
\begin{align}
   \det\left( [T_L, T_R] \right) = & \; \det\left(  a \, \mathbf{e} \, \sigma^x + i (b \, \mathbf{e}^2 + c) \, \sigma^y + d \,\mathbf{e} \, \sigma^z \right)  \nonumber \\
  = & \; a^2 \mathbf{e}^2 - (b \mathbf{e}^2 + c)^2 + d^2 \mathbf{e}^2,     \label{eq:int_cond}
\end{align}
where
\begin{align}
 a = & \; \Gamma_L^2 - \Gamma_R^2, \nonumber \\
 b = & \; \eta_y^R - \eta_y^L,   \nonumber \\
 c = & \; \eta_y^R \left\{ 1 - \Gamma_L^2  \right\} - \eta_y^L \left\{ 1 - \Gamma_R^2  \right\},  \nonumber \\
 d = & \; 2 \left( \eta_y^R - \eta_y^L \right),
\end{align}
with $\Gamma$ as defined in \eq{eq:gamma_def}.

For zero modes we set $\mathbf{e} = 0$, and the matching condition reduces to $c=0$.
A trivial solution to that condition is $\eta_y^R = \eta_y^R = 0$, which corresponds to overlapping Fermi arcs for the two side. For $\eta_y^R, \eta_y^R \neq 0$, the matching condition can be expressed elegantly as
\begin{equation}
 \frac{1 - \Gamma_L^2}{\eta_y^L} = \frac{ 1 - \Gamma_R^2 }{\eta_y^R}.   \label{eq:cond_w0}
\end{equation}
Each of these is a solution of $\det[T_L, T_R] = c^2 = 0$ with multiplicity two. Tuning away from $\kappa=1$ splits these solutions apart, as we discuss next.

\subsubsection{Hard interface for arbitrary $\kappa$}   \label{sec:int_eg_2}
\newcommand{K}{K}
We set $W=2$ and
\[
  \wt{M}_0 = M_L, \quad \wt{M}_1 = M_R, \quad \wt{J}_0 = \kappa J,
\]
so that explicitly,
\[
  J_n = \begin{cases}
	  \kappa J, & n = 0, \\
	  J, & n \neq 0,
        \end{cases}, \quad
  M_n = \begin{cases}
	  M_L, & n \leq 0, \\
	  M_R, & n \geq 1,
        \end{cases}.
\]
The matrix $S$ is required to satisfy $S\Phi_0 = \Phi_2$. To compute it, consider the recursion relations for $n=0,1$ (involving $\kappa$), which can be explicitly written as
\begin{align*}
  \kappa J \Psi_1 + M_L \Psi_0 + J^\dagger \Psi_{-1} = \mathbf{e} \Psi_0,    \nonumber \\
  J \Psi_2 + M_R \Psi_1 + \kappa J^\dagger \Psi_0 = \mathbf{e} \Psi_1.
\end{align*}
These can be solved by a construction similar to that discussed in \ref{sec:wsm}. Explicitly, for $n=0$, we get
\begin{align}
 \alpha_0 = & \; \kappa \mathbf{v}^{\dagger} \mathcal{G}_L \mathbf{v} \, \beta_{1} + \mathbf{v}^{\dagger} \mathcal{G}_L \mathbf{w} \, \alpha_{-1},  \nonumber \\
 \beta_0 = & \; \kappa \mathbf{w}^{\dagger} \mathcal{G}_L \mathbf{v} \, \beta_{1} + \mathbf{w}^{\dagger} \mathcal{G}_L \mathbf{w} \, \alpha_{-1},
\end{align}
where $\mathcal{G}_A = \left( \mathbf{e} \mathbbm{1} - M_A \right)^{-1}, \; A = L,R$. We get a similar expression for $n=1$, and constructing the transfer matrices, we get
\begin{align}
 \Phi_1 & \; = \diag{\frac{1}{\kappa}, 1} \cdot T_L \Phi_0, \nonumber \\
 \Phi_2 & \; = T_R \cdot \diag{1, \kappa} \Phi_1,
\end{align}
so that
\begin{equation}
  S = T_R K T_L, \quad K = \diag{\frac{1}{\kappa}, \kappa}.
\end{equation}
The condition derived in \eq{eq:cond_det} thus simplifies to
\begin{equation}
  \det\left( [T_L, K^{-1} T_R K] \right) = 0.    \label{eq:cond_w2}
\end{equation}
For the transfer matrices of \eq{eq:tmat}, we simplify \eq{eq:cond_w2} using Mathematica\texttrademark $\,$ to get
\begin{align}
 &  \frac{1 - \Gamma_L^2}{\eta_y^L} - \frac{ 1 - \Gamma_R^2 }{\eta_y^R}   \nonumber \\
 &  =  \; \pm 2 \sinh(\lambda) \left[ \frac{\left( 1 - \Gamma_L^2 \right) \left( 1 - \Gamma_R^2 \right)}{\eta_y^L \eta_y^R}  + 4 \cosh^2(\lambda) \right]^{\frac{1}{2}} \!\!\!\!\!,   \label{eq:cond_kappa}
\end{align}
where $\lambda = \ln\kappa$. This reduces to \eq{eq:cond_w0} for $\kappa = 1$, since the RHS vanishes. We can also explicitly note the splitting of each solution to \eq{eq:cond_w0} into two as $\kappa$ is tuned away from 1.

Finally, we remark that \eq{eq:cond_w2} is invariant under the substitution
\begin{equation}
  T_L \to T_R, \quad T_R \to T_L, \quad \kappa \to \frac{1}{\kappa}.
\end{equation}
Since $\kappa$ is the coupling across the interface, this can be thought of as a ``duality'' between the strong and the weak coupling regime, albeit in a noninteracting context. Note that this duality also switches the eigenvalues, so that the physical states are not invariant under it. Instead, it is the curves defined by the matching conditions that are invariant, and the physical states in the two cases correspond to two branches of the same curve.

\begin{figure}[b]
	\centering
	\includegraphics[ width=0.7\columnwidth]{ed_torus.png}
	\caption{
	  A schematic depiction of the system used for the ED. The system is inhomogeneous along $x$ with two interfaces between the two WSMs depicted in orange and green. The transverse direction corresponds to $\mathbf{k}_\perp \in \mathbb{T}^2$.
	}
	\label{fig:ed_sch}
\end{figure}

\subsubsection{Interface with an on-site potential}
Finally, we consider the case of a single layer on-site potential at the interface to illustrate the computation of $S$. Explicitly, we take $W = 1$ and set
\[
  J_n = J, \quad
  M_n = \begin{cases}
	  M_L, & n \leq -1, \\
	  \wt{M}_0, & n = 0, \\
	  M_R, & n \geq 1,
        \end{cases}
\]
as shown in lower panel of \fig{fig:int_eg}. Here, $\wt{M}_0$ represents an on-site potential at layer $0$, which can model irregularities and/or other local effects at the interface. The recursion relation for $n=0$ can be rewritten as
\begin{equation}
  J \Psi_{1} + J^\dagger \Psi_{-1} = (\mathbf{e} \mathbbm{1} - \wt{M}_0) \Psi_0,
\end{equation}
which can be rewritten as $\Phi_1 = \wt{T}_0 \Phi_0$ following the construction of Sec.~\ref{sec:wsm}. Substituting $S = \wt{T}_0$ in \eq{eq:cond_det} results in the matching condition.

\begin{figure*}
	\centering
	\includegraphics[ width=0.16\textwidth]{2node_k000.png}
	\includegraphics[ width=0.16\textwidth]{2node_k014.png}
	\includegraphics[ width=0.16\textwidth]{2node_k033.png}
	\includegraphics[ width=0.16\textwidth]{2node_k100.png}
	\includegraphics[ width=0.16\textwidth]{2node_k300.png}
	\includegraphics[ width=0.16\textwidth]{2node_k700.png}
	\caption{
	  The analytically computed locii of zero modes on the interface BZ, superimposed on the ED result for the 2-node WSM (\sect{sec:eg_2}) for $\kappa = 0, 1/7, 1/3, 1, 3, 7$ (left to right). The background depicts the lowest conduction band on the surface BZ obtained from ED, with darker shades denoting lower energies. The curves in red and green depict the decaying and growing branches of the curve obtained from the matching condition, while the dashed gray line denotes the remaining solutions. Note that the curves for $\kappa$ and $1/\kappa$ are identical, but the growing and decaying labels have exchanged their roles.
	}
	\label{fig:2node_1}
\end{figure*}

\section{Examples}     \label{sec:eg}
In this section, we derive the interface zero modes for specific models for WSMs and the interface described in \sect{sec:int_eg_2}. We compare the analytical results with the interface modes obtained by numerically diagonalizing the Hamiltonian (as a function of the transverse momentum) on a finite lattice in real space. Since a finite system with a single interface necessarily has two open surfaces, it would lead to additional modes in the spectrum obtained from ED. We instead put the system on a 3-torus, with two interfaces between the topological phases, as depicted schematically in \fig{fig:ed_sch}.

We shall restrict ourselves to the case of tunable couplings, in which case the modes localized at the two interfaces have identical zero modes. To see this, note that the RHS of \eq{eq:cond_kappa} is invariant $L \leftrightarrow R$, while the LHS flips sign. Consequently, one switches between the two branches, switching the decay conditions, but that is precisely the correct decay conditions for the other interface.

\subsection{2-node model} \label{sec:eg_2}
We consider the simplest possible WSM, \emph{viz}, one with a single pair of Weyl nodes, described by the Bloch Hamiltonian of \eq{eq:hlt_orig}, with
\begin{equation}
  {\boldsymbol{\eta}} = \left( \nu \sin k_y, 1 - \cos k_y - \cos k_z \right),
\end{equation}
where $\nu \in \mathbb{R}, \, \nu \neq 0$ is a parameter. The two Weyl nodes are located at $\mathbf{k} = (0,0,\pm  \pi/2)$, whose chiralities are $\pm \sgn{\nu}$. The Fermi arcs are given by
\begin{equation}
  \sin k_y = 0, \quad 1 - \cos k_y - \cos k_z < 0,
\end{equation}
which reduces to $k_y = 0,  \abs{k_z} < \pi/2$.

Switching the sign of $\nu$, the only tunable parameter in this model, switches the chiralities of the Weyl nodes and thus the orientation of the Fermi arcs. Consider then an interface between WSMs with $\nu_R = -\nu_L = \nu$, so that the positive and negative chirality nodes swap places across the interface. To evaluate the matching condition explicitly, note that \mbox{$\Gamma_L = \Gamma_R \equiv \Gamma$}. The equation \mbox{$\eta_y^L = \eta_y^R$} is solved for $k_y = 0, \pi$, while the condition of \eq{eq:cond_w0} reduces to
\begin{equation}
  \frac{1 - \Gamma^2}{ -\nu \sin k_y} = \frac{1 - \Gamma^2}{ \nu \sin k_y}  \implies \Gamma^2 = 1,
\end{equation}
or more explicitly,
\begin{equation}
  \nu^2 \sin^2 k_y + \left( 2 + \cos k_0 - \cos k_y - \cos k_z \right)^2 = 1.   \label{eq:arc_2node}
\end{equation}
Tuning away from $\kappa=1$, the matching condition from \eq{eq:cond_kappa} can be simplified to get
\begin{equation}
  \Gamma^2 = 1 \pm 2 \nu \sin k_y \, \sinh \lambda,
\end{equation}
where $\lambda = \ln \kappa$. In \fig{fig:2node_1}, we plot a sequence of these curve as a function of $\kappa$ alongwith the results of ED. The $\kappa \leftrightarrow 1/\kappa$ duality implies that the analytic curves obtained from the matching condition are identical for the two cases. The physical Fermi arcs, however, switch branches under  this duality

For the decoupled WSMs ($\kappa = 0$), the Fermi arcs coincide and run along the same direction. From the intuitive arguments of \sect{sec:int}, we then expect a ``level repulsion''. Indeed, for $\kappa \leq 1$, we get a pair of curves touching at the Weyl nodes. For $\kappa > 1$, these curves expand further away from the Weyl nodes, until they hit the edges and eventually flip around the BZ around $\kappa \approx 6.85$, with additional left over curves running along the noncontractible loops of the interface BZ. In the next example, we see a more dramatic version of this detachment of lines of zero modes from the Weyl nodes.

\begin{figure*}
	\centering
	\includegraphics[ width=0.16\textwidth]{4node_p026.png}
	\includegraphics[ width=0.16\textwidth]{4node_p035.png}
	\includegraphics[ width=0.16\textwidth]{4node_p045.png}
	\includegraphics[ width=0.16\textwidth]{4node_p050.png}
	\includegraphics[ width=0.16\textwidth]{4node_p060.png}
	\includegraphics[ width=0.16\textwidth]{4node_p070.png}
	\ \\ \ \\
	\includegraphics[ width=0.16\textwidth]{4node_k025.png}
	\includegraphics[ width=0.16\textwidth]{4node_k035.png}
	\includegraphics[ width=0.16\textwidth]{4node_k042.png}
	\includegraphics[ width=0.16\textwidth]{4node_k047.png}
	\includegraphics[ width=0.16\textwidth]{4node_k060.png}
	\includegraphics[ width=0.16\textwidth]{4node_k095.png}
	\caption{
	  The analytically computed locii of zero modes on the interface BZ, superimposed on the ED result for the 4-node WSM (\sect{sec:eg_4}) for a fixed $\varphi_0 = \frac{5\pi}{4}$, with (top panel) fixed $\kappa = 1$ and varying $\Delta\varphi = 0.26\pi, 0.35\pi, 0.45\pi, 0.50\pi, 0.60\pi, 0.70\pi$ from left to right, and (bottom panel) fixed $\Delta\varphi = \frac{\pi}{4}$ for $\kappa = 0.25,0.35,0.42,0.47,0.60,0.95$ from left to right. The background depicts the lowest conduction band on the surface BZ obtained from ED, with darker shades denoting lower energies. The curves in red and green depict the decaying and growing branches of the curve obtained from the matching condition, while the dashed gray line denotes the remaining solutions. For the variation of $\Delta\varphi$ as well as $\kappa$, as one tunes these parameters beyond the range explored in the plots, the interface modes gap out completely.
	}
	\label{fig:4node}
\end{figure*}

\subsection{4-node model: tunable Fermi arcs}    \label{sec:eg_4}
We next consider the family of tight-binding models proposed in Ref.~\onlinecite{vd-str_wsm}, which can be written as
\begin{equation}
  {\boldsymbol{\eta}}(\varphi, \mathbf{k}_\perp) = \mathcal{R}(\varphi) \vb(\mathbf{k}_\perp),
\end{equation}
where $\varphi \in [0, 2\pi)$, $\mathcal{R}(\varphi) \in \mathrm{SO}(2)$ is a rotation matrix and $\vb(\mathbf{k}_\perp) = (b_y, b_z)^T$ a reference function. The Weyl nodes are given by the zeros of $\vb(\mathbf{k}_\perp)$ independent of the parameter $\varphi$, while the Fermi arcs, defined by $\eta_y = 0$, manifestly depend on $\varphi$.

Consider an interface across which we change $\varphi$, so that ${\boldsymbol{\eta}}^A = {\boldsymbol{\eta}}(\varphi_A)$. We begin by noting that
\begin{equation}
  \Gamma_A^2 = 1 + \abs{{\boldsymbol{\eta}}}^2 + 2\eta_z^A = 1 + \abs{\vb}^2 + 2\eta_z^A,
\end{equation}
since $\mathcal{R}$ is a rotation matrix. The matching condition of \eq{eq:cond_w0} becomes
\begin{equation}
  \abs{\vb}^2 \left( \eta_y^R - \eta_y^L \right) = 2 \left( \eta_z^R \eta_y^L - \eta_z^L \eta_y^R \right).
\end{equation}
Defining
\[
   \varphi_0 = \frac{\varphi_R + \varphi_L}{2}, \quad
   \Delta \varphi = \frac{\varphi_R - \varphi_L}{2},
\]
we get
\begin{align}
   \eta_y^L - \eta_y^R = & \; 2 \sin(\Delta\varphi) \left( b_y \sin\varphi_0  + b_z \cos\varphi_0 \right),	 \nonumber \\
   \eta_z^R \eta_y^L - \eta_z^L \eta_y^R = & \; 2 \sin(\Delta\varphi) \cos(\Delta\varphi) \abs{\vb}^2.
\end{align}
For an interface $\varphi_R \neq \varphi_L (\text{mod } 2\pi)$, the condition for interface modes becomes
\begin{equation}
  b_y \sin\varphi_0  + b_z \cos\varphi_0 + 2 \cos(\Delta\varphi) = 0.    \label{eq:cond_4node}
\end{equation}
We could cancel off the prefactor $2 \sin(\Delta\varphi) \abs{\vb}^2$, since the zeros of $\abs{\vb}^2$ are the Weyl nodes and the zeros of $\sin(\Delta\varphi)$ satisfy $\varphi_R - \varphi_L = 2n\pi$.

In Ref.~\onlinecite{vd-str_wsm}, a particular choice of $\vb(\mathbf{k}_\perp)$ was considered:
\begin{equation}
  \label{eq:4nodewsm}
  \vecenv{\eta_y}{\eta_z} = \left(
  \begin{array}{cc} \cos\varphi & \;\; -\sin\varphi \\ \sin\varphi & \;\; \;\;\; \cos\varphi  \end{array}
  \right) \vecenv{\cos k_y}{\cos k_z}.
\end{equation}
This model exhibits four Weyl nodes located at momenta $\mathbf{k} = \left( 0, \pm \frac{\pi}{2}, \pm \frac{\pi}{2} \right)$ with chiralities as shown in \fig{fig:wsm_engg_sch}. The nodal configuration is independent of $\varphi$, while the Fermi arcs are given by
\begin{equation}
  0 = \eta_y = \cos\varphi \cos k_y - \sin\varphi \cos k_z,
\end{equation}
which manifestly depends on $\varphi$. One can thus tune the Fermi arc connectivity using a \emph{bulk parameter}, without changing the position and chiralities of the nodes.

We can finally solve the conundrum stated at the end of \sect{sec:int_qual}. For the coupling $\kappa = 1$, the interface modes are given by \eq{eq:cond_4node}. This has a nontrivial solution for $\mathbf{k}_\perp$ iff
\[
  \cos(\Delta\varphi) \leq \frac{1}{\sqrt{2}} \implies \abs{\varphi_R - \varphi_L} \in \left( \frac{\pi}{2}, \frac{3\pi}{2} \right),
\]
in which case one gets a closed curve of zero modes on the surface BZ. Furthermore, note that except for $\Delta\varphi = \frac{\pi}{2}$, these curves do not intersect the projections of the Weyl nodes on the interface BZ, so that they are essentially 2D Fermi surfaces localized deep inside the bulk of a 3D WSM. We plot some of these curves juxtaposed on the corresponding ED results in the top panel of \fig{fig:4node}.

This result is unexpected and particularly intriguing, since the WSMs on the two sides of the interface have identical nodal configurations and are therefore identical from a linearized theory perspective.
We note that although there are modes localized at the interface, they are not topologically protected; however, they are \emph{stable} against ``small'' perturbations. More precisely, they can be gapped out by the addition of an on-site term at the interface, but only by shrinking them to a point. This is analogous to the case of Weyl nodes, which are stable against small perturbations but can be gapped out by a collision with another Weyl node of opposite chirality.

We can also vary the coupling $\kappa$ for a fixed $\Delta\varphi$ to obtain a set of interesting interface nodal configurations. In the bottom panel of \fig{fig:4node}, we plot a few such cases for a fixed $\Delta\varphi = \pi/4$ and a sequence of $0 < \kappa < 1$. The curves for $\kappa > 1$ are identical to those for $1/\kappa$, and the corresponding physical states can be obtained simply by swapping the decaying and the growing branches, as in the 2-node case.

\section{Discussion}     \label{sec:conc}
This article presents a systematic analytic study of localized modes at the interface between two noninteracting topological phases described by Bloch Hamiltonians. A closed form algebraic condition for the existence of such localized modes is derived using the transfer matrix approach. The results are used to study zero modes localized at the interface between WSMs, and compared with the zero modes obtained from a numerical diagonalization of the real space Hamiltonian.

It is worth emphasizing that the results obtained in \sect{sec:int_gen} are rather general, and can be used to study interface modes between any pairs of noninteracting topological phases. This approach, being a real space method, could also be generalized to analytically study other real space structures that might support localized modes, notable examples being defects and dislocations. The essential task is the computation of the scattering matrix $S$ introduced in \sect{sec:int_gen}, which, generically being a single particle scattering problem, should be analytically tractable.

A particular advantage of the transfer matrix method is a complete analytic knowledge of the wavefunctions associated with the bulk as well as the boundary states. Indeed, the study of real space wavefunctions can open interesting avenues for the realization of interesting topological phases in various lattice models\cite{flore_surf_state}. Furthermore, the algebraic nature of various conditions obtained in this picture can lead to alternative perspectives on the topological invariants associated with the well-known topological phases\cite{hatsugai_cbs, vd-vc_tm}.

The geometric information associated with the Fermi arcs has already been a topic of much discussion in the literature. The example discussed in \sect{sec:eg_4} provides yet another manifestation of this aspect, where one gets stable interface modes even when the two sides have identical nodal configuration, and are hence considered ``topologically identical''. The stable 2D Fermi surface of essentially 2D zero modes thus realized on the interface BZ, alongside the inherently 3D zero modes at the Weyl nodes opens up further interesting possibilities, e.g, the superconducting instabilities of this Fermi surface.

The interface modes, unlike the surface modes, lie deep within the 3D bulk and are generally inaccessible to experimental techniques like ARPES. However, the Fermi arcs contribute to the transport, and are needed to \emph{complete the circuit} for the cyclotron orbits in presence of a magnetic field. Signatures of this effect are visible in the quantum oscillations of the density of states\cite{potter_quant_osc}, which might be generalizable to the case of the interface modes.

The study of states obtained at the interface between various topological phases has already led to the realization of many exotic possibilities. It can reasonably be hoped that an analytical approach to this problem will lead to a further understanding of these possibilities, as well as of the topological phases themselves.

\acknowledgments
I was funded by the Deutsche Forschungsgemeinschaft (DFG) with the CRC network TR 183 (Project B03).
I would like to thank Victor Chua, Ciar\'an Hickey and Maria Hermanns for useful discussions.