\section{Introduction \label{sec:intro}}
A three-dimensional topological insulator (3DTI) is an extraordinary material with an insulating band gap but topologically protected gapless surface states. These states are possible due to spin-orbit interaction and follow a spin-polarized Dirac-type spectrum \cite{Qi_RMP}. 
Among many unique properties of these materials, one has triggered an intense interest in the condensed matter community: the emergence of Majorana fermions when TIs are in electrical contact with a superconductor \cite{Alicea_Majorana,*Beenakker_Majorana}. 
Indeed, the topological insulator-superconductor system can be used to engineer an effective \textit{spinless} triplet superconductor which is expected to host zero-energy topologically protected Majorana bound states \cite{Fu_2008,*Sau_2010,*Lutchyn_2010,*Alicea_2010,*Oreg_2010,*Potter_2010,*Potter_2011,*Prada_2012,*Houzet_2013,*Tiwari_2013,*Tiwari_2014b,*Rex_2014}. 


\begin{figure}[ht!]
	\includegraphics[width=\columnwidth]{figure1.eps}
	\caption{\label{fig:sketch}
		(a) Sketch of a normal-superconductor junction on the surface state of a topological insulator (regions $N$ and $S$, respectively) with an intermediate region $I$ modeled as a square potential barrier (see potential profile). The angle of incidence $\theta$ is measured with respect to the $x$-direction. 
		(b) An out-of-plane magnetization $m_z$ is induced in the system by proximity from a ferromagnetic insulator. 
		(c) An external magnetic field, stemming from the vector potential $\mbf{A}$ (blue arrow) with arbitrary orientation, also generates an effective magnetization (red arrows). (d) Symmetry classification of the anomalous Green function $F_{\sigma\sigma'}(\omega,\mbf{k})$. With respect to the exchange of the time (frequency), spin, and spatial (momentum) coordinates, the Green function can be even (E) or odd (O). The spin part can be divided into a singlet (S) or three triplet (T) components. The Pauli principle allows for four different combinations; using a ``frequency/spin/momentum'' notation: ESE, OSO, ETO, and OTE. }
\end{figure}

Remarkably, the interplay between the spin-momentum locking at the surface of the 3DTI, the conventional pair potential from a superconductor, and magnetic order from either a ferromagnet or a magnetic field creates a very exotic induced pair potential in the 3DTI. 
First, the proximity induced pair potential from a conventional singlet $s$-wave superconductor on the surface of a 3DTI features a superposition of singlet and triplet states \cite{Read_2000,Gorkov_2001,Schnyder_2008,Yokoyama_2012,Tkachov_2012,Tkachov_2013,Tkachov_2013b,Burset_2014}. Owing to the symmetry of Cooper pairs, the singlet (triplet) components are even (odd) under the exchange of the spatial coordinates of the two electrons (even/odd in momentum). At the same time, they are all even under the exchange of time coordinates (even in frequency). 
Second, at the interface between normal and superconducting regions inversion symmetry is broken. Such inhomogeneity generates new pairing components that are odd (even) in momentum if they are a singlet (triplet) state of spin \cite{Tanaka_2007,*Eschrig_2007,*Tanaka_2007b,*Tanaka_2007c,*Tanaka_JPSJ,Yokoyama_2012,Black-Schaffer_2012,*Black-Schaffer_2013,*Black-Schaffer_2013b}. Therefore, the Pauli principle imposes that these new terms must be odd in frequency. 
Finally, as we demonstrate in this paper, magnetic order directly affects all these components and opens the possibility for the exotic odd-frequency terms to become dominant and produce strong changes in transport observables of these systems. 

Among the unconventional components of the pair potential, a new, very rare type of superconducting condensate can be induced in this system: one with a spin-triplet component that is odd in frequency and even in momentum \cite{Berezinskii_1974,*Berezinskii_1974R,*Balatsky_1992}. 
Such exotic \textit{odd-frequency triplet pairing} has been proposed to explain the long-range proximity effect observed in superconductor-ferromagnet hybrids \cite{Bergeret_2001,*Birge_2010}. On such systems, the odd-frequency triplet component can be filtered from the rest of induced components of the superconducting condensate since it is the only one resistant to disorder \cite{[{For more information about odd-frequency triplet pairing on ferromagnet-superconductor hybrids see }][{}] Bergeret_RMP,*Eschrig_2010}. 
In topological systems, spin-polarization can be used to detect this unconventional pairing. Recently, Cr\'{e}pin \textit{et al.} have proposed a detection scheme based on nonlocal conductance measurements on the one-dimensional edge of a two-dimensional TI \cite{Crepin_2015}. 
Moreover, odd-frequency triplet pairing and the emergence of Majorana bound states are related: any system that supports Majorana zero-energy bound states also develops an odd-frequency spin-triplet component of the pair potential \cite{Asano_2013,Galitski_2014,DasSarma_2015,Ebisu_2015,Snelder_2015}. 

In this paper, we consider hybrid systems between 3DTI and superconductors like the one sketched in \fref{fig:sketch}(a). 
We study the proximity effect when such systems are exposed to a magnetization extended over the whole junction. The magnetization can be induced by proximity to a ferromagnetic material (b) or by an external magnetic field (c). 
Previous works on 3DTI-based normal-superconductor (NS) and Josephson (SNS) junctions focused only on transport signatures of such hybrids and the conditions for the formation of Majorana bound states \cite{Tanaka_2009,*Linder_2010b,*Sengupta_2013,Linder_2010,Olund_2012,Snelder_2013,Tkachov_2013,Tkachov_2013b,Tiwari_2014,Lu_2015,Lu_2015b}. 
Here, we study the interplay between topological order and superconducting correlations performing a symmetry analysis of the induced pair potential based on the anomalous Green function. We go beyond previous works \cite{Yokoyama_2012,Black-Schaffer_2012,*Black-Schaffer_2013,*Black-Schaffer_2013b} by including the effect of the magnetization on the superconducting electrode \cite{Tkachov_2015}. 

In the absence of magnetization, singlet and triplet components of the pair potential with even and odd frequency dependences are created at the TI-superconductor interface. However, the even-frequency singlet component is dominant and no low-energy states are formed \cite{Black-Schaffer_2012,*Black-Schaffer_2013,*Black-Schaffer_2013b}. 
For finite magnetization, the odd-frequency triplet components can become dominant, resulting in a rich subgap structure with emerging low-energy peaks. 
We study their effect on the NIS conductance, with I an intermediate region, and the local density of states (LDOS) for two orientations of the magnetization: in- and out-of-plane. 
For the latter, when the strength of the magnetization approaches the magnitude of the chemical potential in the superconductor, the odd-frequency triplet component becomes stronger than the singlet one. Consequently, a zero-energy peak emerges in the NIS conductance. 
We find the same behavior for in-plane case with strength comparable to the superconducting gap and oriented parallel to the NIS interface. 
When the in-plane magnetization is oriented differently, the rotational invariance of the band structure is broken. As a consequence, the NIS conductance becomes asymmetric with the energy; a feature completely different to previous works based on ferromagnet-superconductor junctions on 3DTIs \cite{Tanaka_2009,*Linder_2010b,Linder_2010,Lu_2015}. 
Independently of the orientation of the in-plane magnetization, when its strength is comparable to the superconducting gap, the odd-frequency triplet component becomes similar to or bigger than the singlet one. Consequently, subgap resonances are present in the LDOS and NIS conductance. 

This paper is organized as follows. 
We present the model for the TI-superconductor hybrid system and describe the effect of the magnetization in \sref{sec:model}.  
In \sref{sec:sym}, we describe the symmetry classification of the induced pair potential and study the conditions for the emergence of new symmetry terms for finite magnetization, which we associate with changes in the local density of states. 
In \sref{sec:cond}, we connect the previous results with an experimentally relevant observable like the NIS conductance. 
We present our conclusions in \sref{sec:conc}. 
Finally, in the Appendix, we describe our method to compute the retarded Green function and the transport observables that we examine. 


\section{Model \label{sec:model}}
In our setup, the surface of a 3DTI extends along the $x$-$y$ plane. The superconducting and intermediate regions are located, respectively, at $x>d$ and $0<x<d$ (see \fref{fig:sketch}). The role of the intermediate region is to include any interface scattering between the normal and the superconducting regions. We thus assume perfectly flat and clean interfaces. 
The low-energy electron and hole excitations at the surface of the 3DTI are described by the Bogoliubov-de Gennes (BdG) equations $\check{H}\Psi=E\Psi$, with $E$ the excitation energy. 
Particle-hole symmetry imposes that if $\Psi_E$ is a solution of the BdG equations with excitation energy $E$, $\Psi_{-E}$ must also be a solution. 
In Nambu (particle-hole) and spin space, with basis $\Psi=[\hat{c}_{\uparrow}(\mbf{k}),\hat{c}_{\downarrow}(\mbf{k}),\hat{c}^{\dagger}_{\uparrow}(-\mbf{k}),\hat{c}^{\dagger}_{\downarrow}(-\mbf{k})]^T$, with $\hat{c}_{\sigma}(\mbf{k})$ the annihilation operator for an electron of spin $\sigma=\uparrow,\downarrow$ and momentum $\mbf{k}$, the Hamiltonian reads as
\begin{equation}
 \check{H}=\left(\!\begin{array}{cc}
                       \hat{h}(\mathbf{k})-\mu(x)\hat{\sigma}_{0} & i\hat{\sigma}_{2}\Delta(x) \\ -i\hat{\sigma}_{2}\Delta(x) & \mu(x)\hat{\sigma}_{0}-\hat{h}^{*}(-\mathbf{k})
                      \end{array}\!\right) \quad ,
\label{eq:hamil}
\end{equation}
where $\mu(x)$ is the chemical potential and the Pauli matrices $\hat{\sigma}_{0,1,2,3}$ act on spin space. 
We model the intermediate region as a square potential with thickness $d$ and height $\mu_I$, as it is shown in \fref{fig:sketch}(a). Therefore, we have
\begin{equation*}
 \mu(x)= \left\{ \begin{array}{cl} 
  \mu_N &, x\leq0 \\ 
  \mu_I &, 0<x<d \\ 
  \mu_S &, x\geq d 
 \end{array}\right. \quad .
\end{equation*}
We note that, due to Klein tunneling, the barrier does not affect the normally incident modes (with $k_y\sim0$) \cite{Tkachov_2012,Tkachov_2013}. 
The pair potential in the superconducting region is $\Delta(x)=\Delta\Theta(x-d)$, with $\Delta\ge0$ and $\Theta(x)$ the Heaviside function. Electron-like quasiparticles are described by a Dirac Hamiltonian as
\begin{equation}\label{eq:dirac}
 \hat{h}(\mathbf{k})= v_F\left(k_x\hat{\sigma}_{1}+k_y\hat{\sigma}_{2} \right)-e\mbf{A}\cdot\mbf{\hat{\sigma}} \quad ,
\end{equation}
with $v_F$ the Fermi velocity, $e$ the electron charge, and $\mbf{A}$ the vector potential (we take $\hbar=1$ and $c=1$ for simplicity). 
For a constant magnetic field, the vector potential depends linearly on the spatial coordinates. 
In what follows, we describe two situations where we can neglect the spatial dependence of $\mbf{A}$ and regard it as constant. 
The effect of a uniform constant $\mbf{A}$ has the same form as a Zeeman interaction where the vector potential plays the role of an exchange field. Consequently, the Hamiltonian in \eref{eq:hamil} is equivalent to the one used in, e.g., Ref.~\onlinecite{Tanaka_2009}, when we define an \textit{effective magnetization}
\begin{align}\label{eq:magnet}
 \check{M}={}&\left(\!\begin{array}{cc}
                       \mbf{m}\cdot\mbf{\hat{\sigma}} & 0 \\ 0 & -\mbf{m}\cdot\mbf{\hat{\sigma}}^{*}
                      \end{array}\!\right) \\ 
  ={}&-e\left[ \left(A_x\hat{\sigma}_{1}+A_z\hat{\sigma}_{3}\right)\otimes\hat{\tau}_{3} +A_y\hat{\sigma}_{2}\otimes\hat{\tau}_{0}\right]\quad , \nonumber
\end{align}
where the Pauli matrices $\tau_{0,1,2,3}$ act in Nambu space. 
The main difference with previous works \cite{Tanaka_2009,Linder_2010b,Snelder_2013} is that magnetic order is induced in the whole junction and is not just limited to an intermediate region. 
We consider this a sensible approach to describe 3DTI hybrid junctions where proximity-induced superconductivity on the two-dimensional surface state can be sensitive to magnetic effects. 

With our choice of the reference frame (see \fref{fig:sketch}), the momentum component parallel to the interface $k_y$ is conserved. Therefore, the solutions of the BdG equations are expressed as $\Psi(x,y)=\mathrm{e}^{ik_yy}\Psi(x,\theta)$, with $\sin\theta=k_y/|\mbf{k}|$. A full list of the scattering states for this junction is presented in Appendix \ref{sec:app1}. Next, we provide the main results for two different orientations of the effective magnetization $\check{M}$. 

\subsection{Out-of-plane effective magnetization}
We start by considering an out-of-plane effective magnetization $\check{M}=m_z\hat{\sigma}_{3}\otimes\hat{\tau}_{3}$ (see \fref{fig:sketch}), 
stemming from a magnetic field such that $\mbf{A}=[0,0,A_z(x)]$. Ignoring the spatial dependence of the vector potential is equivalent to assume that the magnetization is a Zeeman term with exchange field $m_z$. 

In the superconducting region, the positive branch of the low-energy spectrum is then given by
\begin{equation}\label{eq:mz-energy} 
 E_{\pm}(\mbf{k}) = \sqrt{ \mu^2+\Delta^2 +v_F^2|\mathbf{k}|^2 + m_z^2 \pm 2 r_Z} \quad ,
\end{equation}
with $r_Z=\sqrt{\mu^2v_F^2|\mathbf{k}|^2+m_z^2 \left( \mu^2 + \Delta^2\right)}$. 
Although the energy spectrum looks rather complicated, its analysis is very simple if we focus on the $\Gamma$ point $\mbf{k}=0$, where we find
\[
 E_{\pm}(0) = |m_z \pm \sqrt{\mu^2+\Delta^2}| \quad .
\]
It is straightforward to see that both $m_z$ and the superconducting gap $\Delta$ independently open a gap in an otherwise gapless Dirac spectrum. When both $m_z$ and $\Delta$ are finite at the same time, a competition between the superconducting and the magnetic gap occurs. As a result, a zero-energy state can be created when $m_z^2=\mu^2+\Delta^2$. This allows us to define the effective superconducting gap
\begin{equation}
 \Delta_z\equiv\Delta\sqrt{1-(m_z/|\mu|)^2} \quad .
\label{eq:gaps_mz}
\end{equation}
For $0<m_z\le\mu$, the gap is real and superconductivity is the dominant process in the superconducting region. On the other hand, superconductivity is suppressed by a magnetic gap when $m_z>\mu>\Delta$. 
The critical point between the two regimes is thus $m_z=\mu>\Delta$. The Fermi energy for 3DTIs based on Bi compounds or strained HgTe can be estimated as $\mu\gtrsim50$ meV. A comparable exchange field can be obtained by proximity from a ferromagnetic insulator, as sketched in \fref{fig:sketch}(b). In a recent experiment, EuS was deposited on top of the topological insulator Bi$_2$Se$_3$. Ferromagnetic order induced onto the surface of the TI was estimated to be around $(60-400)$ meV/area [nm$^2$] per applied Tesla \cite{Moodera_2013}. For small enough devices, this correction can be much larger than the induced superconducting gap on TIs for fields below the critical field of the superconductor. 

\subsection{In-plane effective magnetization}
We now examine \eref{eq:magnet} in the heavily-doped, weak-field approximation with $\mu\gg |E|,\Delta,|\mbf{m}|$. Since the out-of-plane component becomes relevant only when $m_z\sim\mu$, we can ignore it in this approximation. The resulting in-plane effective magnetization is $\hat{M}=m_x\hat{\sigma}_{1}\hat{\tau}_{3}+m_y\hat{\sigma}_{2}\hat{\tau}_{0}\equiv m_{\parallel}\left(\hat{\sigma}_{1}\hat{\tau}_{3}\cos\alpha+\hat{\sigma}_{2}\hat{\tau}_{0}\sin\alpha\right)$, with the angle $\alpha$ measured from the $k_x$ direction [see more details in \fref{fig:sketch}(c)]. 
The in-plane effective magnetization is defined from the vector potential $\mbf{A}=[A_x(z),A_y(z),0]$ which we approximate by its value at the junction's plane. This approximation is also valid for the superconducting region if we let the penetration length of the superconductor be much larger than the superconducting coherence length \cite{Tkachov_2015,Sauls_1997,*Tanaka_2002,*Tanaka_2009b}. 

In the superconducting region, the positive branch of the low-energy spectrum is given by
\begin{equation}
 E_{\pm}(\mbf{k})\approx\mbf{n}\cdot\mbf{m} \!+\!\sqrt{\left(\mu\pm v_F|\mbf{k}|\right)^2+\Delta^2} \quad ,
\label{eq:hdwf_energy}
\end{equation}
where we have defined $\mbf{n}=\mbf{k}/|\mbf{k}|=(\cos\theta,\sin\theta,0)$. The effect of the weak field is a shift in the energy, analogous to the Doppler shift discussed for a SQUID-like geometry of helical edge states in Ref.~\onlinecite{Tkachov_2015}. 
Using that
\[
 \mbf{n}\cdot\mbf{m} = m_{\parallel}\left(\cos\theta\cos\alpha + \sin\theta\sin\alpha\right)= m_{\parallel} \cos\left(\theta-\alpha\right) \quad ,
\]
we can define the effective gaps
\begin{equation}
 \Delta_{\pm}\equiv |\Delta\pm m_{\parallel} \cos\left(\theta-\alpha\right)| \quad .
\label{eq:gaps_mxy}
\end{equation}
The critical value for which the effective gap closes is now $m_{\parallel}\!\sim\!\Delta$, which can be reached with rather weak fields. 
The magnetic field induces a finite Cooper pair (condensate) momentum which affects quasiparticles differently depending on their direction of motion. 
Due to the spin-momentum locking at the surface state of the TI, quasiparticles moving in opposite directions feel reversed magnetic fields. In other words, they travel upstream or downstream with respect to the superconducting condensate. This asymmetric behavior with the direction of motion is not found in non-topological systems \cite{Sauls_1997,*Tanaka_2002,*Tanaka_2009b,Burset_2009}. 

\subsection{Retarded Green function and experimental observables}
We now construct the retarded Green function associated to the Hamiltonian in \eref{eq:hamil}. The retarded Green function is obtained combining all scattering states as \cite{McMillan_1968,*Furusaki_1991,*Tanaka_1996,*Kashiwaya_2000,*Herrera_2010}
\begin{gather}\label{eq:GF-gen}
 \check{G}^r(x,x',E+i0^+;\theta)= \\
 \left\{\!\! \begin{array}{lr}
   \begin{array}{c} 
               \alpha_1\Psi_3(x)\tilde{\Psi}_1^T(x') \!+\! 				\alpha_2\Psi_3(x)\tilde{\Psi}_2^T(x') \\
               + \alpha_3\Psi_4(x)\tilde{\Psi}_1^T(x') \!+\! \alpha_4\Psi_4(x)\tilde{\Psi}_2^T(x')
   \end{array} \, , & x< x'
\\ 
   \begin{array}{c} 
               \beta_1\Psi_1(x)\tilde{\Psi}_3^T(x') \!+\! 				\beta_2\Psi_2(x)\tilde{\Psi}_3^T(x') \\
               + \beta_3\Psi_1(x)\tilde{\Psi}_4^T(x') \!+\! \beta_4\Psi_2(x)\tilde{\Psi}_4^T(x')
   \end{array} \, , & x> x'
               \end{array} \right. \quad . \nonumber
\end{gather}

An incoming electron (hole) from the normal region onto the superconductor can be 
\begin{inparaenum}[(\itshape a\upshape)]
\item Andreev reflected as a hole (electron); 
\item normal reflected as an electron (hole); 
\item transmitted to the superconductor as an electronlike (holelike) quasiparticle; and
\item as a holelike (electronlike) quasiparticle. 
\end{inparaenum}
We label this scattering process with the number $1$ ($2$). The corresponding wave function is thus $\Psi_{1(2)}(x)$. 
Analogously, the equivalent process where the incoming particle is an electronlike (holelike) excitation incident from the superconductor is labeled with the number $3$ ($4$). 
In Appendix \ref{sec:app2}, we show more details of the calculation of the Green function, the reflection amplitudes $a_j$, $b_j$, and the transmission amplitudes $c_j$ and $d_j$, with $j=1,\dots,4$. 
The wave functions $\tilde{\Psi}_i(x)$ correspond to the conjugate scattering processes to $\Psi_i(x)$. 
The coefficients $\alpha_i$ and $\beta_i$ are determined by the continuity equation obtained after integrating the BdG equations around $x=x'$, namely, 
\begin{gather}\label{eq:GF-cont}
 \frac{-i}{v_F}\hat{\sigma}_{1}\otimes \hat{\tau}_{0} = \\ \check{G}^r(x,x+0^+,E+i0^+;\theta)-\check{G}^r(x,x-0^+,E+i0^+;\theta) \quad . \nonumber
\end{gather}

The spectral density of states is then calculated from the retarded Green function as
\begin{equation*}
 A(x,E;\theta)=-\frac{1}{\pi}\imag\left\{\mathrm{Tr}\check{G}_{ee}^r(x,x,E+i0^+;\theta)\right\} \quad ,
\end{equation*}
where the trace is taken on the electron-electron component of the Green function in Nambu space (i.e., the single-particle Green function). Furthermore, the local density of states (LDOS) is given by
\begin{equation}
 \rho(x,E)=\int\limits_{-\pi/2}^{\pi/2}\!\!\! \mathrm{d}\theta\cos\theta A(x,E;\theta) \quad .
\label{eq:ldos-def}
\end{equation}
In the next sections, we normalize the LDOS in the superconductor by the LDOS calculated deep inside the normal region, $\rho_N(E)=|\mu_N+E|/v_F$ (see Appendix \ref{sec:app3} for more details). 

We define the conductance of the NIS junction as \cite{BTK,Tanaka_1995}
\begin{equation}
 G_{NS}(E)=\int\limits_{-\pi/2}^{\pi/2}\!\!\! \mathrm{d}\theta\cos\theta \left[ 1-|b_1(E,\theta)|^2+|a_1(E,\theta)|^2\right] \quad ,
\label{eq:cond_NS}
\end{equation}
where $a_1(E,\theta)$ and $b_1(E,\theta)$ are the Andreev and normal reflection amplitudes, respectively, for an incident electron. 
We normalize the conductance using the constant value $G_0=G_{NS}(E\gg\Delta)$. 

Notice that both the conductance and the LDOS are averaged over all incident angles, labeled by $\theta\in[-\pi/2,\pi/2]$ \footnotemark[1]. 

\footnotetext[1]{The one-dimensional edge of a 2DTI, which is equivalent to the $\theta=0$ case considered here, presents a well-defined spin polarization. For the 2D surface state, the modes with $\theta\neq0$ usually spoil this neat effect. }


\section{Symmetry of the induced pair potential \label{sec:sym}}
Induced superconducting correlations are manifested in the anomalous Green function (the electron-hole component of the Green function in Nambu space), i.e., the $2\times2$ matrix in spin space
\[
 \tilde{F}_{\sigma\sigma'}(t_1,\mbf{r}_1;t_2,\mbf{r}_2)=\left\langle T_t\Psi_{\sigma}(t_1,\mbf{r}_1)\Psi_{\sigma'}(t_2,\mbf{r}_2) \right\rangle \quad ,
\]
with $T_t$ the time-ordering operator and $\Psi_{\sigma}(t_1,\mbf{r}_1)$ the fermion operator for given spin, time, and space coordinates. In homogeneous systems, we can define the relative coordinates $\mbf{r}=\mbf{r}_1-\mbf{r}_2$ and $t=t_1-t_2$ and the Pauli principle imposes that $\tilde{F}_{\sigma\sigma'}(t,\mbf{r})=-\tilde{F}_{\sigma'\sigma}(-t,-\mbf{r})$. After Fourier transformations, we have
\[
 F_{\sigma\sigma'}(\omega,\mbf{k})=-F_{\sigma'\sigma}(-\omega,-\mbf{k}) \quad .
\]
We can write explicitly the spin-singlet and triplet components of the anomalous Green function in the basis of \eref{eq:hamil} as
\[
 \hat{F}(\omega,\mbf{k})= \left[ F_0(\omega,\mbf{k})\hat{\sigma}_{0} + F_j(\omega,\mbf{k})\hat{\sigma}_j\right]i\hat{\sigma}_{2} \quad .
\]
Evidently, the Pauli principle must be fulfilled independently by the singlet and the three triplet components. 
With respect to the spin degree of freedom, the singlet component is odd and the triplet ones even. Under a sign change in either the frequency or the momentum, the functions $F_\nu(\omega,\mbf{k})$, with $\nu=0,1,2,3$, can also be classified as even or odd. As a consequence, only four types of symmetries are allowed, namely, \begin{inparaenum}[(\itshape i\upshape)]
\item \textbf{E}ven-frequency spin-\textbf{S}inglet \textbf{E}ven-parity (ESE); 
\item \textbf{E}ven-frequency spin-\textbf{T}riplet \textbf{O}dd-parity (ETO); 
\item \textbf{O}dd-frequency spin-\textbf{S}inglet \textbf{O}dd-parity (OSO); and
\item \textbf{O}dd-frequency spin-\textbf{T}riplet \textbf{E}ven-parity (OTE) 
\end{inparaenum} [see \fref{fig:sketch}(b) and, for more details, Ref.~\onlinecite{Tanaka_2007,*Eschrig_2007,*Tanaka_2007b,*Tanaka_2007c,*Tanaka_JPSJ}]. 

In our model, we assume perfectly flat and clean interfaces which conserve the parallel component of the momentum ($k_y$ for the axis choice in \fref{fig:sketch}). In the quasi-one dimensional model, transport along these junctions is described by one-dimensional Green functions $\check{G}(x,x',\omega;k_y)$ separately for each channel $k_y$. 
For positive frequencies, we construct these Green functions using \eref{eq:GF-gen}, which considers scattering solutions of the BdG equations with excitation energies $E\ge0$. For negative frequencies, we need to consider the time-reversal of the scattering processes, i.e., the $E<0$ solutions of the BdG equations, which are provided by the advanced Green function. 
Therefore, to study the frequency dependence of the anomalous Green function, we combine the retarded and the advanced Green functions as 
\begin{gather}\label{eq:anom-GF}
 F_{\nu}[E+i\sgn(E)0^+] = \\ \Theta(-E)F^a_{\nu}(E-i0^+)+\Theta(E)F^r_{\nu}(E+i0^+) \quad , \nonumber
\end{gather}
where $E$ is a real variable, that coincides with the excitation energy for positive values, which we associate with the frequency. $k_y$ is parametrized using the angle of incidence $\theta$ and the advanced Green function is obtained using $\check{G}^a(x,x',E-i0^+;k_y)=[\check{G}^r(x',x,E+i0^+;k_y)]^{\dagger}$. 

\subsection{Zero-field analysis}
We now analyze each component of the anomalous Green function. 
We are mainly interested in local observables, like the LDOS which is determined by the spatial variation of the symmetry of the induced pair potential. Therefore, we only consider the case with $x=x'$ here. 
In Appendix \ref{sec:app2}, we show that the Green function is composed of an \textit{edge} and a \textit{bulk} part, i.e., one term given by the scattering at the interface and another by the solutions far away from the interface, respectively. For each component of the anomalous Green function we thus have
\begin{equation*}
 F_{\nu}(E,\theta)= F_{\nu,\infty}(E,\theta)+ F_{\nu,A}(E,\theta) + F_{\nu,B}(E,\theta) \quad .
\end{equation*}
$F_{\nu,\infty}$ corresponds to the Green function at $x\!=\!x'\!\rightarrow\!\infty$, i.e., the \textit{bulk} term inside the superconducting region. On the other hand, $F_{\nu,A}$ and $F_{\nu,B}$ are the \textit{edge} Green functions proportional to Andreev and normal reflection processes, respectively. In what follows, we ignore $F_{\nu,B}$ since it only includes rapid oscillations. 

In the absence of magnetization, for the bulk terms we immediately see that 
\begin{align*}
 F_{0,\infty}(E,\theta)={}& \frac{i\sgn\left(E\right)}{2v_F|\cos\theta|} \frac{\Delta} {\sqrt{E^2-\Delta^2}} \quad , \\
 F_{2,\infty}(E,\theta)={}& \frac{i\sgn\left(E\right)}{2v_F|\cos\theta|} \frac{\Delta} {\sqrt{E^2-\Delta^2}}\sin\theta \quad , \\
 F_{1,\infty}(E,\theta)={}& F^r_{3,\infty}(E,\theta) = 0 \quad ,
\end{align*}
where 
\begin{gather*}
 \sqrt{E^2-\Delta^2} = \\ \left\{ \begin{array}{cr}
               \sgn\left(E\right)\sqrt{[E\!+\!i\sgn\left(E\right)0^{+}]^2-\Delta^2} & , \, |E|>\Delta \\ i\sgn\left(E\right)\sqrt{\Delta^2-[E\!+\!i\sgn\left(E\right)0^{+}]^2} & , \, |E|\le\Delta  \end{array} \right. \quad ,
\end{gather*}
which is an odd-function of the energy. Therefore, we find that $F_{0,\infty}$ and $F_{2,\infty}$ are even-frequency functions; they are thus classified as ESE and ETO terms, respectively. Even in the absence of magnetic order, the combination of a spin-orbit induced Dirac spectrum with conventional spin-singlet superconductivity results in a superposition of singlet and triplet superconducting correlations \cite{Read_2000,Gorkov_2001,Schnyder_2008,Yokoyama_2012,Tkachov_2012,Tkachov_2013,Tkachov_2013b,Burset_2014}. 
The ETO term is proportional to $\sin\theta=k_y/|k|$ because it is odd in momentum. Therefore, after averaging over incident angles, any odd-momentum terms must vanish. 
Consequently, the effect of ETO terms on a local observable in our system is canceled by the average over incident angles. 

At the IS interface with $x=d$, the breaking of spatial inversion generates odd-frequency terms both for the singlet and triplet components \cite{Black-Schaffer_2012,*Black-Schaffer_2013,*Black-Schaffer_2013b,Tanaka_2007,*Eschrig_2007,*Tanaka_2007b,*Tanaka_2007c,*Tanaka_JPSJ}. Angle-averaging cancels any odd-frequency singlet term, but the odd-frequency triplet components remain. 
We thus find
\[
 F_{1,A}(E,\theta)=\frac{i}{v_F}a_3(E,\theta)\mathrm{e}^{i\frac{2\sqrt{E^2-\Delta^2}}{v_F\cos\theta} x} \quad ,
\]
with $a_3(E,\theta)$ the Andreev reflection amplitude for an electron-like quasiparticle incident from the superconductor, which results in an OTE contribution that survives the integration over incident angles. 

To compare the relative weight of the triplet components against the singlet one, we define the magnitude of the angle-averaged triplet vector 
\begin{equation}\label{eq:triplet}
 F_t(E)=\left[ \sum\limits_{j=1}^{3} \left| \int\limits_{-\pi/2}^{\pi/2}\!\!\!\mathrm{d}\theta \cos\theta F_{j}(E,\theta) \right|^2 \right]^{1/2}\quad .
\end{equation}
Due to the integration over incident angles, $F_t(E)$ represents only odd-frequency triplet terms. Therefore, at zero field, $F_t(E)=|F_{1,A}(E)|$. It is interesting to note that, in our spin basis, this term is given by $F_1\propto \left(\uparrow\up-\downarrow\dw\right)$, while the opposite polarization, $F_2\propto \left(\uparrow\up+\downarrow\dw\right)$, vanishes. This is a consequence of our calculation being limited to one surface of the 3DTI; for the other surface the polarization of the Cooper pairs at the interface is reversed. 


\begin{figure}
	\includegraphics[width=\columnwidth]{figure2.eps}
	\caption{\label{fig:fdos_z}
	 (a) Evolution of $\rho(E=0)$, angle-averaged $|F_{0}|$, and $F_t$ with the distance inside the superconducting region from the interface at $x=d$ at zero magnetization and $E=0$. (b) Same quantities as before as a function of  $m_z$, for $E=0$ and $x=d$. The green solid and dashed lines correspond to the two triplet polarizations $|F_{1,2}|\propto(\uparrow\up\mp\downarrow\dw)$, respectively. (c,d) LDOS as a function of the energy for several values of $m_z$ calculated at the IS interface $x=d$ (c), and inside the superconducting region, where $x\rightarrow\infty$ (d). For all plots, we set $\mu_N=10\mu$, $\mu_I=20\mu$, $\mu_S=\mu=10^3\Delta$, and $d=0.1(v_F/\mu)$. }
\end{figure}


We show in \fref{fig:fdos_z}(a) the evolution of the angle-averaged quantities $|F_{0}(E=0)|$ and $F_t(E=0)$ as a function of the distance inside the superconducting region (dashed blue and dotted black lines, respectively). The bulk ESE contribution is dominant and responsible for the suppression of the density of states (red solid line) for $x>\xi$. At the interface, however, we find that $F_{t}(E=0)\neq0$. The OTE contribution vanishes inside the superconducting region but is responsible for the enhancement of the LDOS close to the interface. 

\subsection{Symmetry analysis for out-of-plane effective magnetization}
We now study the local anomalous Green function for an out-of-plane effective magnetization. We present a list of the components of the anomalous Green function in Appendix \ref{sec:app4}, for the perpendicular case treated here and for the in-plane one analyzed in the next section. The presence of the perpendicular magnetization affects all components and, in particular, induces two new bulk terms. $F_{1,\infty}$ is proportional to $k_y=|k|\sin\theta$ (ETO) and becomes zero after integration over $\theta$. On the other hand, $F_{3,\infty}$ is an OTE term proportional to $E m_z$ and it vanishes in the absence of magnetization or for $E=0$. Therefore, the main effect on the bulk region is determined by $F_0$. 

The situation is different at the interface. We plot in \fref{fig:fdos_z}(b) angle-averaged $|F_{0}(E=0)|$ and $F_t(E=0)$ (blue dashed and black dotted lines, respectively), as a function of $m_z$, for $x=d$. We choose $E=0$ to analyze the symmetry of any zero-energy peak in the LDOS (red solid line). 
The \textit{edge} contributions of the triplet components have odd-frequency parts that remain finite after integrating over the incident angles. $F_t(0)$ becomes comparable to $|F_{0}(0)|$ for a magnetization strength $m_z\approx0.7\mu$. 
Additionally, we plot the triplet components $|F_{1,2}(E=0)|\propto(\uparrow\up\mp\downarrow\dw)$ using solid and dashed green lines, respectively. At $m_z=0$, only $|F_1(0)|$ is present. For a finite $m_z$, $|F_2(0)|$ becomes a finite OTE term and increases with the field strength like $|F_1|$ at the same rate that the ESE term is suppressed. The LDOS at the interface is enhanced accordingly. 
With the increase of $m_z$, the feature of a Majorana resonant state becomes more pronounced. Also, here Majorana resonances are always accompanied by odd-frequency pairing. 

The magnetization reduces the bulk gap $\Delta_z$, given in \eref{eq:gaps_mz}, for $m_z\leq\mu$. After the gap closes, the LDOS and the anomalous terms of the Green function are strongly suppressed. 
We show the energy dependence of the LDOS at the interface with $x=d$ in \fref{fig:fdos_z}(c) and in the bulk of the superconducting region with $x\rightarrow\infty$ in \fref{fig:fdos_z}(d). 
At the interface, the LDOS features a strong resonance at $|E|=\Delta_z$ and is finite for $|E|<\Delta_z$. 
For $m_z>\mu$, the LDOS is still finite and maximum at $E=0$ [see \fref{fig:fdos_z}(b)]. However, the value of this finite LDOS at the interface decreases rapidly when increasing $m_z$ over $\mu$ because it is suppressed by a magnetic gap. 
Deep inside the superconducting region, where only the ESE component of the anomalous Green function is finite, we find a BCS-type density of states around the effective gap $\Delta_z$. For $m_z>\mu$, a magnetic gap opens destroying superconductivity. 


\begin{figure}
	\includegraphics[width=\columnwidth]{figure3.eps}
	\caption{\label{fig:LDOS_mxy}
	(a) Energy-momentum plot of the low-energy bands inside the superconducting region for several values of $m_{\parallel}$. (b) Bulk LDOS of the superconducting region of the NIS junction as a function of the excitation energy $E$ for different $m_{\parallel}$. (c,d) LDOS at the interface $x=d$ with $\mu_I/\mu=10$ and $d=0.1(v_F/\mu)$ as a function of the excitation energy $E$ with different in-plane magnetizations. We show $\alpha/\pi=0,0.5$ for (c,d), respectively. 
	For (b,c,d) we set $\mu=10^{3}\Delta$. }
\end{figure}

\subsection{Symmetry analysis for in-plane effective magnetization}
In the previous section, we showed that the effect of a magnetic field in the heavily-doped, weak-field approximation is a Doppler-type shift of the excitation energy. This shift has an important effect on the energy bands: in \fref{fig:LDOS_mxy}(a), we plot the conduction (red) and valence (blue) bands in momentum space for several values of $m_{\parallel}$. For reference, the zero-energy plane is marked by a gray area. At zero field, the spectrum has rotational symmetry and the superconducting gap is the same for $|\mbf{k}|=k_F$. For finite fields, $m_{\parallel}\neq0$, the rotational symmetry of the spectrum is broken. However, particle-hole symmetry is recovered using $\epsilon(\mbf{k})=-\epsilon(-\mbf{k})$. The orientation of the NIS interfaces with respect to this bulk spectrum becomes very important since the energy spectrum can be asymmetric around $k_y=0$, with $k_y$ the momentum component parallel to the interface. Only the case where the in-plane field is oriented 
parallel to the normal-superconductor interface, i.e., for $\alpha=\pi/2$, the energy bands are $k_y$-symmetric. 
Additionally, there is a \textit{critical} magnetization $m_{\parallel}=\Delta$ where both bands touch the zero-energy plane and the superconducting gap closes at some values of $|\mbf{k}|=k_F$, for any orientation of the in-plane field. The orientation of the in-plane field $\alpha$ changes the position in momentum space where this closing occurs \cite{[{The gap closing can greatly affect the edge states of the two-dimensional TI, as recently described in }] [{}]Reinthaler_2015}. 

To study the effect of the gap closing on the bulk of the superconducting region, we show the bulk density of states in \fref{fig:LDOS_mxy}(b) for several values of $m_{\parallel}$. After averaging over the angle of incidence, the bulk density of states is independent of the orientation of the magnetization. When $0<m_{\parallel}<\Delta$, the effective gaps $\Delta_{\pm}$ defined in \eref{eq:gaps_mxy} emerge and the LDOS exhibits a peak at $E\sim|\Delta_+|$ and is gapped for $E\lesssim|\Delta_-|$. When the gap closes for $m_{\parallel}=\Delta$, the bulk density of states becomes V-shaped with two resonances at $E=\pm2\Delta$. 
For this case, the LDOS at the interface features a sharp zero-energy peak for any orientation of the in-plane field, as it is shown in the dashed-dotted blue lines of \fref{fig:LDOS_mxy}(c,d). 
Beyond this point, the zero-energy peak becomes wider if $\alpha<\pi/4$ or splits into two otherwise. For the bulk LDOS, a region emerges around $E=0$ with a flat density of states of the same magnitude as that of the normal region. Outside this energy range, the V-shaped profile becomes wider with increasing the strength of the field. 


\begin{figure}
	\includegraphics[width=\columnwidth]{figure4.eps}
	\caption{\label{fig:fdos}
	 Plots of $\rho(E=0)$, $|F_{0}(E=0)|$, and $F_t(E=0)$ as a function of $m_{\parallel}$. (a) Bulk results for $x\rightarrow\infty$. (b,c,d) Results at the interface with $x=d$ for different orientations of the effective magnetization: $\alpha/\pi=0,0.25,0.5$ for (b,c,d), respectively. The solid and dashed green lines correspond to the polarizations $\uparrow\up-\downarrow\dw$ and $\uparrow\up+\downarrow\dw$ of the OTE term, respectively. For all plots, $\mu_N=\mu_S=\mu=10^{3}\Delta$, $\mu_I/\mu=10$, and $d=0.1(v_F/\mu)$. }
\end{figure}


We now study the components of the anomalous Green function plotting the bulk contribution of $|F_{0}(E=0)|$ and $F_t(0)$ as a function of $m_{\parallel}$ in \fref{fig:fdos}(a). There are two distinct regimes. For $m_{\parallel}<\Delta$, the ESE term is greatly enhanced, while the OTE is finite but small compared to it. The LDOS is suppressed which corresponds to a conventional gapped profile. 
For the critical magnetization $m_{\parallel}=\Delta$, all terms increase. 
At this point, the LDOS jumps from $0$ to $1$. From then on, the triplet vector magnitude is finite and similar to the ESE term for $m_{\parallel}\gtrsim 2\Delta$ while the LDOS remains constant. Therefore, a finite bulk LDOS is obtained because of an equal superposition of OTE and ESE terms in the bulk of the superconducting region. 

We show in \fref{fig:fdos}(b,c,d) the ESE and OTE components at the interface ($x=d$), as a function of $m_{\parallel}$ for different orientations of the in-plane magnetization [$\alpha=0$, $\pi/4$, and $\pi/2$ for Figs. 4(b), 4(c), and 4(d), respectively]. 
For the symmetric orientation $\alpha=\pi/2$, the OTE term becomes dominant over the singlet one for values equal or greater than the critical magnetization [\fref{fig:fdos}(d)]. 
In the opposite case with $\alpha=0$ [\fref{fig:fdos}(b)], the singlet term is always dominant over the triplet one. Independently of the orientation, the LDOS is greatly enhanced for values greater than the critical one. 

Finally, the behavior of $|F_{1,2}|\propto(\uparrow\up\mp\downarrow\dw)$ strongly depends on the orientation of the magnetization (see green solid and dashed lines for $|F_{1,2}|$, respectively). Indeed, for $\alpha=0$ we find that $|F_{1}|=|F_{2}|$ while they are very different for $\alpha\neq0$. Therefore, the spin polarization of the induced Cooper pairs can be controlled by the orientation of the external field \footnotemark[2].  

\footnotetext[2]{The orientation of the in-plane vector potential $\mbf{A}$ and, hence, that of the effective magnetization $\mbf{m}_{\parallel}$, is shifted from that of the magnetic field by $-\pi/2$. }

\section{Conductance spectroscopy \label{sec:cond}}
Our setup, sketched in \fref{fig:sketch}, is ideal for conductance spectroscopy measurements like the ones recently performed in Refs. \onlinecite{Finck_2014,Brinkman_2015}. The conductance of the NIS junction, obtained from \eref{eq:cond_NS}, is normalized to $G_0=G_{NS}(E\gg\Delta)$, as it is commonly done in experiments. 


\begin{figure}
	\includegraphics[width=\columnwidth]{figure5.eps}
	\caption{\label{fig:NScond}
	 Conductance as a function of the excitation energy for several effective magnetizations. For each case we include a sketch of the NIS junction displaying the orientation of the effective magnetization as a blue arrow. (a) A perpendicular magnetization is finite in the whole NIS junction with $\mu_N=10\mu$, $\mu_I=20\mu$, and $\mu_S=\mu$. The thin gray lines are taken at intervals $0.1\mu$ from $m_z=0$. (b,c) In-plane magnetization finite in the whole junction oriented perpendicular (b) or parallel (c) to the NIS interfaces. For both cases, $\mu_N=\mu_S=\mu$ and $\mu_I/\mu=10$. (d) With the same parameters as (b), superposition of the conductance for both surfaces; sketched on the left. The black solid line corresponds to the $m_{\parallel}/\Delta=1.6$ case of (b) while the dashed one is the conductance of the other surface with the same parameters. In all plots, $d=0.1(v_F/\mu)$ and $\mu=10^3\Delta$. }
\end{figure}


We start considering an out-of-plane effective magnetization finite in the whole junction. 
Consequently, when we approach the critical point where $m_z\sim\mu$, a magnetic gap starts to develop in both normal and superconducting regions. In order to have propagating states in the normal region, we choose $\mu_N=10\mu$ and $\mu_S=\mu=10^3\Delta$. We still include a square potential barrier setting $\mu_I=20\mu$. The results for this setup are shown in \fref{fig:NScond}(a) for several values of $m_z$ in the range $0<m_z<\mu$. 
For $m_z=0$, the tunnel conductance features a gapped profile. Due to the Klein tunneling effect, the conductance is not fully gapped for $E<\Delta$, even in the presence of a strong barrier. 
For $0<m_z<\mu$, the effective gap is reduced and, as $m_z$ approaches $\mu$, the subgap conductance increases clearly featuring a zero-bias peak. 
The formation of this zero-bias peak coincides with the regime where the odd-frequency triplet $F_t$ becomes dominant with respect to the ESE term $|F_0|$ [see \fref{fig:fdos_z}(b)]. For $m_z>\mu$, the gap in the excitation spectrum is given by $E_z=|m_z-\sqrt{\mu^2+\Delta^2}|$ and the conductance is zero for $|E|\leq E_z$. 

We now consider an in-plane magnetization in the heavily-doped, weak-field approximation. Under this approximation, if we take $\mu_N\sim\mu_S\gg |E|,\Delta,m_{\parallel}$, the effect of the field in the normal and intermediate regions is almost negligible. Assuming a setup with a symmetric square potential barrier with $\mu_N=\mu_S=\mu$ and $\mu_I=10\mu$, we plot in \fref{fig:NScond}(b,c) the conductance for several values of $m_{\parallel}$. 
We show an in-plane magnetization perpendicular to the NIS interfaces ($\alpha=0$) in \fref{fig:NScond}(b) and parallel to them ($\alpha=\pi/2$) in \fref{fig:NScond}(c). 
In both cases, the effect of the field is to split the superconducting gap into two effective gaps $\Delta_{\pm}$ defined in \eref{eq:gaps_mxy}. 

For $\alpha=\pi/2$, the resonance at $E=\Delta$ for $m_{\parallel}=0$ remains pinned to the smallest of the gaps, but the one at the biggest gap is smeared by the angle average. Due to Klein tunneling, these resonances can not be clearly resolved in the conductance, although they appear in the LDOS, as shown in \fref{fig:LDOS_mxy}(c,d). The smallest of the gaps closes for $m_{\parallel}\gtrsim\Delta$ and a zero-energy peak appears in the conductance. It is in this regime where the odd-frequency triplet term is greater than the singlet one, as shown in \fref{fig:fdos}(d). 

Strikingly, for $\alpha\neq\pi/2$, the conductance is asymmetric with the energy. 
This asymmetry is a consequence of the distortion of the energy bands by the magnetic field, as shown in \fref{fig:LDOS_mxy}(a). The field breaks the rotation symmetry of the energy bands and the modified bands have a point-like symmetry that is determined by the orientation of the in-plane magnetization. Since we study transport in the direction perpendicular to the NIS interfaces, when the magnetization is aligned with the interface, i.e., $\alpha=\pi/2$, the symmetry is recovered [see \fref{fig:NScond}(c)], while it is lost for any other orientation. 

Owing to the special particle-hole symmetry of our system, we find that, for $\alpha\neq\pi/2$, Andreev reflection is very asymmetric with the energy while normal reflection is symmetric. This effect holds independently of the barrier strength provided that the effective magnetization is misaligned with the NIS interface. Thus, as an example, we can consider a junction with $\mu_I=\mu_S\neq\mu_N$. For such junctions, normal and Andreev reflections take place separately at each side of the intermediate region. In the setup of \fref{fig:sketch}(a), when $|E|\le\Delta$, normal (Andreev) reflection is the only scattering process at $x=0$ ($x=d$). The reflection amplitudes are then
\begin{subequations}\label{eq:amplitudes}
\begin{align}
 a_1(E,\theta)={}&\frac{2\cos\theta_N\cos\theta_S}{t(E,\theta)}\Gamma_- \quad , \\
 b_1(E,\theta)={}&i\frac{\sin\theta_N-\sin\theta_S}{t(E,\theta)}\left(1-\Gamma_+\Gamma_-\right) \quad , \\
 t(E,\theta)={}& 1+\cos\left(\theta_N+\theta_S\right) \nonumber \\ & - \Gamma_+\Gamma_-\left[1-\cos\left(\theta_N - \theta_S\right)\right] \quad , \nonumber
\end{align}
\end{subequations}
with $\theta_{N,S}=\sin^{-1}(k_y/\mu_{N,S})$ and 
\begin{align*}
 \Gamma_{\pm}={}&\frac{\Delta^*}{\epsilon_{\pm}+\sqrt{(\epsilon_{\pm})^2-|\Delta|^2}}
 \quad , \\
 \epsilon_{\pm}={}& E\pm m_{\parallel}\cos(\theta_S-\alpha) \quad .
\end{align*}

Due to the spin-momentum locking at the surface of the 3DTI, the magnetic field affects differently quasiparticles moving in opposite directions \cite{Reinthaler_2015}. That is the reason why Andreev reflection is asymmetric with the magnetization, being proportional to $\Gamma_-$ while normal reflection is not, since it goes with $(1-\Gamma_+\Gamma_-)$. 
As a consequence, the conductance becomes more asymmetric as we increase $m_{\parallel}$, as shown in \fref{fig:NScond}(b). The energy range where Andreev reflection is enhanced also becomes wider. In the figure, this happens always for positive energies. 

Finally, the asymmetry in the conductance is also partly a consequence of having limited our analysis of transport to only one surface of the 3DTI. In \eref{eq:hamil}, the Dirac Hamiltonian chosen to describe electron-like quasiparticles is $\hat{h}(\mbf{k})=[v_F \mbf{k}-e\mbf{A}]\cdot\hat{\sigma}$. With this choice, quasi-particles have positive helicity, i.e., their spin points in the direction of their momentum. For the opposite surface state of the 3DTI, however, electron-like quasiparticles are described by $\hat{h}(\mbf{k})=[-v_F \mbf{k}-e\mbf{A}]\cdot\hat{\sigma}$. Therefore, the sign of their helicity changes and their spin points in the opposite direction of their momentum. As a consequence, the effect of the in-plane magnetization is reversed, with $\Gamma_{\pm}\rightarrow\Gamma_{\mp}$. The normal reflection amplitude is unchanged, but the Andreev reflection amplitude becomes proportional to $\Gamma_+$ in \eref{eq:amplitudes}. The asymmetry of the conductance is reversed, as it is shown by the black solid and dashed lines of \fref{fig:NScond}(d), corresponding to the $m_{\parallel}/\Delta=1.6$ case of \fref{fig:NScond}(b). If we consider a conductance measurement that includes opposite surfaces of the 3DTI, the symmetry with the energy is recovered. For the same parameters as \fref{fig:NScond}(b), we show such conductance spectroscopy in \fref{fig:NScond}(d). 


\section{Conclusions \label{sec:conc}}
We have analyzed the superconducting proximity effect at 3DTI-superconductor hybrids in the presence of an effective magnetization. We find that the interplay between the magnetization and the spin-momentum locking of the surface state gives rise to interesting odd-frequency triplet terms in the induced pair potential. Contrary to the usual even-frequency triplet terms, averaging over the angle of incidence does not cancel these components, and their effect can be observed in experimental transport observables such as the NIS conductance and LDOS. 

For an out-of-plane magnetization, the OTE terms increase with the strength of the field in the regime $\mu>m_z>0$. At the same time, the singlet term of the induced pair potential is suppressed. As a consequence, the conductance evolves from a conventional gapped profile for $m_z=0$ into a zero-bias conductance peak when $m_z\lesssim\mu$. 
Such strong magnetization can be induced in the system when a ferromagnetic insulator such as EuS is deposited on top of the NIS junction \cite{Moodera_2013}, in a setup like the one sketched in \fref{fig:sketch}(b). 

If we consider the effect of an orbital magnetic field on the junction, we find a similar gap to zero-energy peak transition on the conductance for in-plane magnetization parallel to the NIS interface. 
The main difference is that the critical value that determines the emergence of a zero-energy peak is now $m_{\parallel}=\Delta$, where the energy bands feature an indirect closing of the gap. 
For $m_{\parallel}\ge\Delta$, the OTE term is greater than the ESE one at the interface and similar in magnitude in the bulk superconducting region. As a result, we have connected the emergence of odd-frequency triplet in the induced pair potential on the surface of the 3DTI with a distinctive profile of the conductance. 

When the in-plane magnetization is not oriented parallel to the interface, the conductance becomes asymmetric with the energy. This is a manifestation of the breaking of the rotational symmetry of the band structure of the combined 3DTI-superconductor system. Only a transport measurement that probes one surface state of the 3DTI can detect the asymmetry of the conductance. Such setup, would be a perfect platform to probe the coexistence of topological order and superconductivity. 

Our results show that 3DTI-based NIS junctions feature a very rich induced superconducting pair potential that displays the elusive odd-frequency triplet component. 
Using very basic ingredients such as a conventional superconductor and an external magnetic field or ferromagnetic insulator, standard experimental techniques like conductance spectroscopy can be used to detect signatures of OTE superconductivity. 


\section*{Acknowledgments}
We thank D. Bercioux, F. S. Bergeret and F. Cr\'{e}pin for helpful discussions and comments. 
We acknowledge financial support by the DFG-JST research unit ``Topotronics'', the Helmholtz Foundation (VITI), the DFG Priority Program SPP1666, FOR1162, SFB 1170 ``ToCoTronics'', the ENB Graduate School on ``Topological Insulators'', and DFG Grant No. TK 60/1-1. 
This work was also supported by Topological Materials Science (TMS) (Grant No. 15H05853) and Grant No. 25287085 from the Ministry of Education, Culture, Sports, Science, and Technology, Japan (MEXT), and by the Core Research for Evolutional Science and Technology (CREST)
of the Japan Science and Technology Corporation (JST). 

