\section{INTRODUCTION}
\label{intro}
Nearby galaxies (z $\la$ 0.1) are characterized by large angular scale and high apparent
brightness, which can be studied in more details and with higher accuracy than high-redshift
galaxies. For nearby galaxies, the ultraviolet (UV) imaging provides unique information. Massive,
young stars emit strong UV energy, which dominates the integrated UV light of star-forming
galaxies. Thus, the UV flux is widely used as an excellent and accurate measurement of
the current star formation rate (SFR; \citealt{Wilkins12,Kennicutt98,Lanz13}). In addition,
the interstellar dust can absorb the UV light and then re-emit at far-infrared wavelength.
The comparison between infrared and UV emission can effectively trace the dust attenuation
in galaxies. Radiative transfer models suggest that the ratio of far-infrared to UV
luminosity is a reliable estimator of the dust attenuation \citep{Buat96,Witt00,Panuzzo03},
depending weakly on the geometry of stars and dust, the extinction law, and the stellar
population.

The UV galaxy samples in the local universe are important for understanding the
evolution of galaxies with cosmic time \citep{Bouwens12,Reddy12,Ellis13,Martin05}. There
have been considerable attempts to explore the SFR, morphology, dust attenuation, and their
evolutions by constructing UV-selected local galaxies
\citep{Sullivan04,Calzetti94,Buat02,Kuchinski00,Marcum01,Lauger05}.

Based on \textit{GALEX} Public Release 2 and 3 (GR2/GR3), \citet{Gil07} (hereafter
\textit{GALEX} Atlas) presented a UV
catalog of 1034 nearby galaxies selected from the Third Reference Catalog of Bright
Galaxies (RC3; \citealt{de91}). They selected those galaxies whose optical diameters
of the $\mu_B = 25$ mag arcsec$^{-2}$ isophote (D25) are larger than 1{\arcmin}.
The FUV and NUV images were mainly taken
from \textit{GALEX} Nearby Galaxies Survey (NGS). From the surface brightness profiles,
they obtained asymptotic magnitudes, colors, and concentration indices. In
combination with archival optical and infrared data, they analyzed the color-magnitude and
color-color diagrams, relations between colors and morphologies, IRX$-\beta$ relation,
dust attenuation, and structural properties for different types of galaxies.

In February 2013, \textit{GALEX} GR7 data was available. It includes more than
45000 images, almost three times larger than GR2/GR3. With deeper observations (up to
250 kilo-second), it is now possible to create a relative complete UV sample of nearby
galaxies. In this work, we provide an updated UV catalog of nearby galaxies with
deepest FUV and NUV images from GR6/GR7, archival optical and infrared photometry, and
make some similar analysis to \citet{Gil07}.

In Section~\ref{data}, we present our sample selection of nearby galaxies. In Section~\ref{sbm}, the method
of measuring the surface brightness profile is described in detail. The catalog content including the
UV parameters and other archival data are described in Section \ref{cc}. The parameter comparison
of our catalog with \textit{GALEX} Altas is shown in Section \ref{compga}. Section \ref{aa} presents
some analyses and application examples of our catalog. The summary is in Section \ref{sum}.
All magnitudes presented in the paper are corrected for the Galactic extinction using
the Galactic Reddening Map of \citet{Schlegel98}. We use the reddening law of  \citet{Cardelli89} to
convert $E(B-V)$ to the UV extinction:  $A_{FUV}=7.9E(B-V)$ and $A_{NUV}=8.0E(B-V)$.

\section{Sample Selection}
\label{data}
The nearby galaxies in this paper are extracted from HyperLeda \footnote{\url{http://leda.univ-lyon1.fr, revised in 2013 July 1}}
\citep{Paturel03,Makarov14}. HyperLeda is a catalog of
extragalactic sources, which gives accurate coordinates of objects (typical accuracy
better than 2$\arcsec$), morphological parameters (e.g. size, axis ratio, position
angle, morphogical type, etc.), and astrophysical parameters (optical magnitudes and colors,
surface brightness, distance, kinematic quantities, etc.). It includes 3.7 million galaxies,
much more than RC3 as used by \citet{Gil07}.
The catalog provides Doppler distances for 20314 galaxies and redshift-independent distances
for 2013 galaxies.
Figure \ref{fig1}a shows the distributions of major-axis diameters in both HyperLeda
and RC3 catalogs, while Figure \ref{fig1}b presents the same distributions but
with D25 $>1 \arcmin$. The majority of HyperLeda galaxies have diameters less than 1{\arcmin},
and over 99\% RC3 galaxies are included in HyperLeda.  The HyperLeda galaxies with diameters
above 1{\arcmin} are more than those of RC3 by a factor of 3.

We adopt the same criterion of optical diameters as \citet{Gil07} used, that is, D25 is larger than
1$\arcmin$. There are 22948 galaxies satisfying this constraint. We retrieve FUV and NUV images
for these galaxies from \textit{GALEX} GR6/GR7 to obtain their intensity maps. Three imaging
surveys are involved: Medium Imaging Survey (MIS), Nearby Galaxy Survey (NGS), and Deep Imaging
Survey (DIS). We require the galaxies (1.5 $\times$ D25) to be fully covered by
the \textit{GALEX} field of view (FoV; 1.2{\arcdeg} in diameter). Those galaxies with low
signal-to-noise ratio or contaminated by surrounding galaxies and bright stars are excluded.
The optical diameter D25 of M31 is about 3.0{\arcdeg}, far beyond the \textit{GALEX} FoV.
We obtain 26 tiles of M31 in FUV and NUV from NGS and DIS, which are stacked to
form very deep mosaics by Montage\footnote{\url{http://montage.ipac.caltech.edu/}}.
The final sample contains 4138 galaxies, 2321 of which have both FUV and NUV observations.
There are 14 and 1803 galaxies having only FUV and NUV observations, respectively.

We present their basic parameters given by HypeLeda in Table 1, including the
positions, sizes, morphological types and distances.
For galaxies whose distances are not available, we get their distances from NED (the
NASA/IPAC Extragalactic Database). There are 4047 galaxies with either luminosity distances
or redshift-dependent distances.

Figure \ref{fig2} shows the normalized distributions of D25 diameter, morphological type, and
distance in \textit{GALEX} Atlas and our sample. By contrast, our sample gives more galaxies with
smaller diameters, earlier types, and larger distances.

\section{Surface Brightness Measurement} \label{sbm}
\textit{GALEX} provides sky background subtracted images, but the background
near the galactic center is overestimated, especially for galaxies with large apparent sizes.
We introduce a new procedure to better estimate the sky background, which is quite
similar to the method in \citet{Zou11}. SExtractor \citep{Bertin96} is used to separate
source pixels from background pixels. The source pixels, together with the region with the galactocentric
distance $r < 2 \times $D25, are masked. Then, the remaining pixels are fitted with
a two dimensional polynomial function to generate the background map. Finally, the fitted
map is subtracted from the intensity image.

Foreground stars contaminate the galaxy flux and need to be masked in the intensity map.
We get the point sources (mostly stars) within the image from \textit{GALEX} photometric
catalogs. The mask aperture is determined by the growth curve, where the corresponding
aperture magnitude is close to the total magnitude presented in the catalog. However, the
catalogued point sources might include some star formation regions, whose UV color is much bluer
than normal stars. Since the number of star formation regions is strongly dependent on the galaxy
type, we separate our entire sample into two subsamples according to the morphological
type: one with T $\geq$ $-$0.5 identified as spiral/irregular galaxies (types from $Sa$ to $Irr$)
and the other with T $<$ $-$0.5 identified as elliptical/lenticular galaxies (types from $E$
to $S0a$).
For early-type galaxies, all point sources detected by \textit{GALEX} are regarded
as foreground stars and masked.
For late-type galaxies, we only mask the point sources with (FUV $-$ NUV) $>$ 1,
and keep bluer point sources with (FUV $-$ NUV) $<$ 1 as star forming regions.
The nuclei of galaxies, also point-like, are not masked.
For galaxies with single UV-band observations, we define foreground stars by visually
inspecting both the UV image and the color image of the Digitized Sky Survey.

We then use the MATLAB package of astronomy and astrophysics for the measurement \citep{Ofek14}.
The radial surface brightness profile is computed with a series of elliptical annuli, whose
major radii range from 3$\arcsec$ to at least 1.5 times the D25 radius. The brightness
errors are estimated with the method in \citet{Gil05}. Figure \ref{fig3}
shows some examples of the false color images and surface brightness profiles of galaxies
in our sample. The radial (FUV $-$ NUV) color profile is also plotted.

\section{The Catalog Content} \label{cc}
\subsection{Asymptotic and D25 Magnitudes}
We adopt the technique introduced by \citet{Cairos01} to calculate the asymptotic
magnitudes. First, the accumulated flux and the gradient of the growth curve
at each radius are computed. Then, an appropriate radial range away from the
galaxy center is visually chosen, where an error-weighted linear fit to the relation of
accumulated flux vs. the gradient of the growth curve is performed.
Finally, the intercept of this linear
fit is regarded as the asymptotic magnitude, and corresponding magnitude error
is derived from the error of the linear fit.

We also calculate the aperture magnitudes inside the D25 ellipses. The errors
are estimated from the random noise inside the apertures and the uncertainties
of the sky background subtraction (see \citealt{Gil05}). All the errors do not
include the uncertainty of NUV and FUV photometric zeropoints, which are about 0.15 and
0.09 mag, respectively.

\subsection{Effective surface brightness and concentration indices}
Based on the growth curve, we derive the radii containing 50\%, 20\%, 25\%,
75\%, 80\% of the total luminosity ($r_{\mathrm{eff}}$,$r_{20}$, $r_{25}$,
$r_{75}$, and $r_{80}$, resepectively). The effective surface brightness is
calculated as $I = 0.5L_{\mathrm{FUV}}/\pi r^{2}_\textrm{eff}$ and the
concentration indices are calucated as C31=$r_{75}$/$r_{25}$ and
C42 = 5log($r_{80}$/$r_{20}$) \citep{de77,Kent85}.

\subsection{Archival optical and infrared data}
Multi-wavelength data of nearby galaxies could improve our understanding of their
physical nature.
We add some corollary data from other optical and infrared surveys
into our catalog. HyperLeda provides the total $B$ magnitude ( $B_{T}$) and the total
$I$ magnitude ($I_{T}$). A total of 4064 galaxies have $B_{T}$ magnitudes and 3265
galaxies have $I_{T}$ magnitudes. There are 640 and 813 galaxies in our sample with
$(U-B)_{T}$ and $(B-V)_{T}$ colors, respectively. In near infrared (NIR), we compile
$JHK$-band photometry from the 2MASS Large Galaxy Atlas (LGA) \citep{Jarrett03}. If
galaxies are unavailable in LGA, $JHK$ photometry in 2MASS Extended Source
Catalog (XSC) are used. A total of 3708 galaxies have $JHK$ magnitudes. In mid and far
infrared (MIR and FIR), the photometry in 12$\mu$m, 25$\mu$m, 60$\mu$m and 100$\mu$m
are obtained from Infrared Astronomical Satellite (IRAS). We follow the priority
given by \citet{Gil07} to compile the photometry (\citealt{Rice88,Moshir90} and the
$IRAS$ Point Source Catalog).
The 60$\mu$m and 100$\mu$m are used to estimate the far infrared emission. There are 1570 galaxies
with infrared fluxes in both 60$\mu$m and 100$\mu$m. The UV properties and archival
optical and IR data are presented in presented in Table \ref{tab2} and Table \ref{tab3}, respectively.

\section{Comparison with \textit{GALEX} Atlas} \label{compga}
We show the photometric comparisons between \textit{GALEX} Atlas and our sample in Figure \ref{fig4}.
Figure \ref{fig4}a--d display the differences of asymptotic and aperture magnitudes between our catalog
and \textit{GALEX} Atlas. There are 94\% and 93\% galaxies in the FUV and NUV  with the difference of
asymptotic magnitudes within 0.5 mag. We check the galaxies with differences larger than 0.5 mag and
find that these differences might be caused by the different sky-background estimation and elliptical
parameters inherited from their parent catalogues. We fit a sky
background map for each galaxy, while \textit{GALEX} Atlas only use a single background value.
For example, the average background of PGC 4085 (ESO 243-G041) given by \textit{GALEX} Atlas
is 5.59 $\times$ 10$^{-4}$ counts s$^{-1}$ in the NUV, while the average value of our computed background
is 2.24 $\times$ 10$^{-3}$ counts s$^{-1}$, which leads to the magnitude difference of 0.59.
After checking the intensity image, our fitted background is more reasonable, since
the flux of no-signal area around the galaxies is about 2 $\times$ 10$^{-3}$ counts s$^{-1}$.
There are many galaxies in \textit{GALEX} Atlas and our sample with different axis ratios and D25 diameters.
When galaxies become fainter, the difference of the axis ratios can reach a factor of 9,
which leads to different growth curves and results in different magnitudes.
For example, the axis ratio of PGC 4190 (NGC 0407) given by HyperLeda is 6.9,
while \textit{GALEX} Atlas gives 4.3.
These differences lead to the different magnitudes between \textit{GALEX} Atlas and our sample.

Figure \ref{fig4}e and f show the differences between the (FUV-NUV) colors of
asymptotic and aperture magnitudes. The colors of our sample is consistent with
those in the catalogues of \textit{GALEX} Atlas with standard deviations of 0.19 and
0.13 for asymptotic and aperture colors, respectively.

The differences of the concentration indices for the two samples are displayed in Figure \ref{fig4}g and h.
The concentration indices in both FUV and NUV of our samples are on average smaller than those in \textit{GALEX}
Atlas, probably due to the differences of the growth curves.

\section{Analyses and Some Applications} \label{aa}
\subsection{UV Properties of the galaxies}
Figure \ref{fig5} shows the distributions of UV properties of the galaxies
in our catalog. Figure \ref{fig5}a-g presents the histograms of asymptotic magnitudes,
morphological type, luminosities, (FUV $-$ NUV) colors, effective radii,
and concentration indices. The median asymptotic magnitudes in FUV and NUV are
17.55 and 17.20 mag, respectively.
In our sample, about 5$\%$ galaxies show asymptotic magnitudes larger
than D25 magnitudes. For these galaxies, the diameters corresponding to
asymptotic magnitudes are smaller than the D25 diameters due to the lack of
diffuse emission. The residual flux after background subtraction within
D25 ellipses leads to additional UV flux in the D25 magnitudes.

There are 1030 (25\%) and 2899 (70\%) early and late-type galaxies (Figure \ref{fig5}b),
and 209 (5\%) galaxies without morphological types given by HyperLeda.
Figure \ref{fig5}c presents the luminosity distributions . The median FUV
and NUV luminosities are 7.8 $\times$ 10$^8$ and 9.1 $\times$ 10$^8$ $L_{\odot}$
(3.0 $\times$ 10$^{35}$ and 3.5 $\times$ 10$^{35}$ W), respectively.
Ten galaxies in our catalog are ultraviolet-luminous galaxies (UVLGs; \citealt{Heckman05})
defined as the FUV luminosity larger than 2$\times$10$^{10}$ $L_{\odot}$ ($7.6\times10^{36}$ W),
which are extremely rare in the local universe. Among the ten galaxies, PGC 2248, 36466,
59214 have been identified by \textit{GALEX} Atlas. PGC 53898 and 71035 are Seyfert
galaxies \citep{Veron06}. PGC 4007 and 17625 are a luminous infrared galaxies (LIRG; \citealt{Haan11,Sanders03}).
The other three galaxies, PGC 23064, 30400, and 51865, are star-forming
spiral galaxies with very UV-luminous arms. PGC 70348 (NGC 7469) is defined as a UVLG in \textit{GALEX} Atlas but not
qualified in our sample. Our FUV luminosity of this galaxy is about 7.3$\times$10$^{36}$ W
, close to the threshold.

Figure \ref{fig5}d presents the color distributions of asymptotic and D25 magnitudes.
The majority of galaxies has color bluer than (FUV $-$ NUV) = 1, which has been used to separate
foreground stars. The red tail of the distributions is populated by early-type galaxies.
The distributions of effective radii in kiloparsecs are shown in Figure \ref{fig5}e.
The median effective radii for FUV and NUV are 4.3 and 5.0 kpc, respectively.
The distributions and comparisons of concentration indices are presented in
Figure \ref{fig5}f-i. The galaxies are on average slightly
more concentrated in the NUV than in the FUV by differences of 0.03 in C31 and
0.10 in C42,
which may be due to the different fractions of the bulge component
in spiral galaxies in these two bands \citep{Gil07}. Figure \ref{fig5}j shows the color
as a function of the FUV magnitude. The late-type galaxies are brighter and located in a
narrow color range, while the early-type galaxies are fainter and distributed with
broader and redder color span.

\subsection{Effective Surface Brightness}
The plot of effective surface brightness versus the
FUV luminosity is presented in Figure \ref{fig6}. This figure shows a trend of slight increasing
effective surface brightness with increasing luminosity, indicating that
galaxies with large effective radii turn to be more luminous \citep{Hoopes07}.
\textbf{The effective surface brightness distribution of our sample is more concentrated
than \citet{Hoopes07} sample, since the contours in Figure \ref{fig6} enclose 96\%
of the galaxies in our sample but only 84\% in their sample. }

The compact UVLGs are UVLGs defined with a surface brightness $I_{FUV} >$ 10$^{9}$ $L_{\odot}$ kpc$^{-2}$
\citep{Hoopes07}. They have characteristics that are remarkably similar to the Lyman break galaxies
(LBGs), such as SFRs, metallicities, morphologies, kinematics, and attenuations
\citep{Hoopes07,Overzier08,Basu-Zych07}.
Only 12 galaxies in our sample have the surface brightness $I_\mathrm{FUV} >$ 10$^{9}$ $L_{\odot}$ kpc$^{-2}$.
Two of them, PGC 59214 and PGC 71035, have both high luminosity
and high surface brightness, which can be classified as compact UVLGs.
PGC 59214 is a BL Lac object and PGC 71035 is a Seyfert 1.5 \citep{Veron06}.
Both of them have very compact nuclei in the FUV, and the optically thick accretion
disks around super-mass black holes probably dominates their FUV luminosity
\citep{Malkan82,Ward87,Sanders89}.

\subsection{(FUV $- K$) tracing galactic morphology}
Compared to the optical colors \citep{Weiner05,Pozzetti10,Talia14}, the combination
of UV and IR should be more efficient in separating the early and late-type galaxies,
given their different SFRs and dust attenuation \citep{Gil07}.
Since the UV emission is very sensitive to the presence of recent star
formation activity and the $K$ band emission is sensitive to the accumulated star
formation \citep{Munoz07}, the (FUV $- K$) color provides a robust discrimination
between early and late-type galaxies as analyzed by \citet{Gil07}.

Figure \ref{fig7}a shows the relation between (FUV $- K$) color and galactic morphology.
Here we convert the 2MASS $K$ band Vega magnitudes to the AB system by adding 1.84 mag \citep{Munoz09,Cohen03}.
We calculate the dividing (FUV $- K$) color to separate early and late-type galaxies
by minimizing the number of galaxies that are misclassified by the dividing magnitude,
which leads to (FUV $- K$) = 6.84. About 68\% early-type galaxies and 86\% late-type
galaxies are correctly classified. The (FUV $- K$) is linearly correlated to the morphology of
late-type galaxies. The best linear fit to the relation (green line in Figures \ref{fig7})
is
\begin{equation}
  \mathrm{(FUV }- K) = \mathrm{(5.8\pm0.1) - (0.25\pm0.02) \times T}.
\end{equation}
Here T is the morphological type given by HyperLeda and the errors of the
parameters are derived from the linear fitting.

The concentration indices have also been used as a classification tool
\citep{Abraham96,Bershady00}. The (FUV$ - K$) color together with concentration indices are expected
to show better separation of different type galaxies, which is discussed in \citet{Gil07}.
We present the FUV concentration index C42 as a function of the  (FUV $ - K$) color in
Figure \ref{fig7}b. Early-type galaxies have larger C42 and redder (FUV $ - K$). The best
separation line in the plane of this figure is computed as
\begin{equation}
	(\mathrm{FUV }- K) = 14.4\pm0.3 - (2.5\pm0.1) \times \mathrm{C42_\mathrm{FUV}}.
  \end{equation}
With this, 76\% and 94\% of early and late-type galaxies are correctly classified.

\subsection{IRX-$\beta$ relation}
\label{damr}
The IR excess (IRX) is widely used as a good tracer of the dust attenuation in galaxies,
and the slope of the UV continuum ($\beta$) is considered as a proxy to estimate the IRX when
IR data are not available, which is known as IRX-$\beta$ relation \citep{Daddi07,Reddy09,Takeuchi10}.
IRX is defined as log($f_\mathrm{TIR}/f_\mathrm{FUV}$), where $f_\mathrm{TIR}$ is the total infrared flux,
and $\beta$ can be estimated by (FUV $-$ NUV).
Here, the IRX-$\beta$ relation is not valid for early-type galaxies, since a substantial
part of dust heating is due to photons emitted by old stars which do not emit primarily
in UV and the dust attenuation would be overestimated by IRX \citep{Buat11}.

Figure \ref{fig8} shows the (FUV $-$ NUV) color against IRX for our late-type galaxies together
with the local starburst relation derived from \citet{Meurer99}. Here we define
$f_\mathrm{FUV} = \lambda f_{\lambda}$ and adopt the formula proposed by \citet{Sanders96} to derive
the $f_\mathrm{TIR}$, which is estimated using the 12, 25, 60 and 100 $\mu$m from IRAS data. We separate
our galaxies into two groups according to their FUV luminosities.  In Figure \ref{fig8}, we can see
that most of our galaxies are located below the relation of \citet{Meurer99}. The galaxies with brighter
FUV luminosities ($L_{\mathrm{FUV}} \geq 10^{35} W$) have a tigher relation than that of fainter galaxies.
Both groups have steeper IRX$-\beta$ relations than that of \citet{Meurer99}. We linearly fit the relations
of these two groups, which are formulated as
\begin{equation}
	\mathrm{IRX} = (6.54\pm0.01) \times (\mathrm{FUV} - \mathrm{NUV}) - (2.10\pm0.01), L_\mathrm{FUV} \geq 10^{35} \mathrm{W},
\end{equation}
\begin{equation}
	\mathrm{IRX} = (2.76\pm0.05) \times (\mathrm{FUV} - \mathrm{NUV}) - (0.84\pm0.07), L_\mathrm{FUV} < 10^{35} \mathrm{W}.
\end{equation}

\section{SUMMARY}
\label{sum}
In this paper, we provide an updated catalog of 4138 nearby galaxies based on the
latest (GR6/GR7) of \textit{GALEX}, which is more than 3 times the number of galaxies
in the original \textit{GALEX} Atlas.
Our samples are selected from the extragalactic catalog of HyperLinda.
The D25 diameter is set to be larger than 1{\arcmin}. Compared with \textit{GALEX} Atlas, we
apply a more precise procedure to estimate the sky background in the \textit{GALEX} images.
Radial FUV and NUV surface brightness profiles are obtained. From these profiles we
calculate asymptotic magnitudes, aperture (D25) magnitudes, UV colors, structural
parameters (effective radii and concentration indices), luminosities, and effective
surface brightness. We also augment our data set with archival optical and infrared
photometry from HyperLeda, 2MASS and IRAS. With this updated catalog, we confirm
that the (FUV $- K$) color provides a good criterion to distinguish early and late-type galaxies,
which can be improved further with the concentration indices.
The IRX-$\beta$ relation is reformulated with our UV-selected nearby galaxies.

The GALEX images and catalogs of our nearby galaxies can be accessed via the website
\url{http://batc.bao.ac.cn/$\sim$zouhu/doku.php?id=projects:galax:start}

\begin{acknowledgements}
The authors acknowledge support from the National Science Foundation of China under
grants NSFC-11273028, NSFC-11333004 and NSFC-11203031, and support from the National
Astronomical Observatories, Chinese Academy of Sciences under the Young Researcher Grant.
We thank Mark Seibert and Scott Fleming for
providing subfunctions to generate false color images. We acknowledge the usage of the
HyperLeda database\footnote{\url{http://leda.univ-lyon1.fr}}.
This research has made use of the NASA/IPAC Extragalactic Database (NED), which is
operated by the Jet Propulsion Laboratory, California Institute of Technology,
under contract with the National Aeronautics and Space Administration.
Some of the data presented in this paper were obtained from the Mikulski Archive
for Space Telescopes (MAST).
This research has made use of the NASA/IPAC Infrared Science Archive,
which is operated by the Jet Propulsion Laboratory, California Institute of Technology,
under contract with the National Aeronautics and Space Administration.
This research made use of Montage, funded by the National Aeronautics and Space
Administration's Earth Science Technology Office, Computation Technologies Project,
under Cooperative Agreement Number NCC5-626 between NASA and the California Institute
of Technology. Montage is maintained by the NASA/IPAC Infrared Science Archive.

\end{acknowledgements}


\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).


\section{Introduction and Purpose}

Dose verification Quality Assurance (QA) measurements are critical in proton therapy for the validation of patient treatment plans \cite{ICRU78ch6}. These measurements are especially important in Intensity Modulated Proton Therapy (IMPT), which scans a focused proton beam spot over the target region many times at different energies and intensities to improve dose conformity. IMPT carries increased treatment planning complexity relative to passively scattered treatment modalities, making the need for dose verification over all beam energy layers especially critical.

An emerging concept for proton therapy dose verification systems is the use of liquid scintillator phantoms with position sensitive photon sensors to image the 3D dose distribution with the phantom volume. These novel detectors present several advantages over the current QA measurement systems that rely upon the physical movement of a two-dimensional array of gaseous detectors within a water phantom \cite{Beddar2009etal, Kroll2013etal}. In particular, liquid scintillator-based systems permit sub-millimeter spatial resolution \cite{Archambault2012etal}, and have the potential to perform validation measurements over all energy layers in no more time than is required for the patient's exposure to the beam. This real-time dose verification is an especially attractive feature of the active phantom approach, as it facilitates a high patient throughput for clinical operations.

Water-based Liquid Scintillator (WbLS) is a newly-developed water-scintillator emulsion \cite{Yeh2011etal} that is a promising candidate material for next-generation particle physics detectors \cite{Alonzo2014etal}.
As WbLS is predominantly water, it offers improved water-equivalence for a scintillator-filled dose-verification phantom, thereby offering superior dosimetry performance relative to traditional liquid scintillator. WbLS also largely mitigates the flammability and toxicity hazards that are traditionally associated with liquid scintillators. We are currently exploring the possibility of using WbLS for dose verification QA in proton therapy.

One aspect of the use of scintillating materials for dose-verification systems that has yet to be comprehensively investigated is the potential for radiation damage to affect the scintillating performance of the phantom volume. Such damage could lead to systematic measurement errors due to a variable scintillation light yield as a function of dose deposited. In this study, we have quantified the relationship between accumulated dose and scintillating light yield for pure liquid scintillator and WbLS samples. Our measurements of changes to the light yield have been assessed in terms of the
expected effects upon QA detector operation.

\setlength{\tabcolsep}{0pt}
\begin{figure*}[!t]
\centering
\begin{tabular}{cc}
    \subfloat[ ] {\label{subfig:BeamSetup} \includegraphics[trim={2.5cm 0 2cm 0},clip,width=0.49\textwidth]{BeamSetup_withDimensions.pdf}}
    & \subfloat[ ] {\label{subfig:BeamProfile} \includegraphics[width=0.49\textwidth]{BeamProfileWithLiquidVolumesOutlined.pdf}}   \\
\end{tabular}
\label{fig:IrradDetails}
\caption{Details of the proton beam irradiation geometry; (a)
a schematic drawing of the sample holders used in the irradiations, (b)
the proton beam intensity profile incident upon the scintillator volumes (outlined areas), taken using a phosphor screen coupled to a CCD. The fiducial marks in (b)
define the 10 cm by 10 cm region about the beam centre.}
\end{figure*}

\setlength{\tabcolsep}{11pt}

\section{Experiment}

\subsection{Materials}
For the pure liquid scintillator samples we used a Linear Alkyl Benzene (LAB) solvent with 3 g/L 2,5-Diphenyloxazole (PPO) as primary fluorophore and 15 mg/L 1,4-Bis(2-methylstyryl)benzene as a wavelength shifting fluorophore. Our WbLS was composed of 5\% LAB by mass with 2 g/L PPO in the organic phase as the primary fluorophore.

\subsection{Irradiation A}
Proton irradiations were carried out at the NASA Space Radiation Laboratory (NSRL) at Brookhaven National Laboratory. Duplicate samples in polypropylene vials were mounted in two HDPE mounts and exposed to a 10 cm x 20 cm beam of 201 MeV protons (figure \ref{subfig:BeamSetup}).
The uniformity of the beam intensity across the target area was measured to be approximately $\pm 13 \%$ (figure \ref{subfig:BeamProfile}).
Two consecutive beam exposures were used to deliver approximately 300 Gy and a further 500 Gy to samples of pure liquid scintillator and WbLS. Samples from the rear mount were removed following the first exposure, and samples in the front mount were irradiated by both exposures.
The dose accumulated in the rear mount was verified using a Fricke dosimetry solution that was placed in the sample holder along with the scintillating samples. The Fricke measurements suggest that the first exposure delivered a dose of 328 $\pm$ 10 Gy to the Fricke vial. The 800 Gy accumulated by the front holder's samples is twice the maximum dose measurable with Fricke dosimetry \cite{ASTMFricke}, so the NSRL beam instrumentation -- a calibrated ionisation chamber placed upstream from the samples -- was relied upon for the dose estimate.

During the irradiations, online measurements of the scintillation response of liquid scintillator and WbLS samples were taken to search for short term radiation damage effects. Liquid Scintillator and WbLS samples were placed 50 cm behind the HDPE holders and measurements were carried out using photomultipliers and a digital oscilloscope. A gas ionisation chamber was used as a reference to correct for variations in the beam intensity. No short term damage to the light yield of either sample was observed, to a precision of 3.3\% and 2.4\% for the WbLS and liquid scintillator samples, respectively. The dose delivered to the samples is estimated to be 556 Gy, using the simulation model described below.

Empty vials were also placed in the beam area and later filled with scintillator and water-based liquid scintillator to check for radiation-induced leaching from these vials. Measurements of the scintillating materials that resided in these vials for approximately 1 month revealed no change in the light yield, suggesting that vial leaching does not affect our results.

\begin{table*}[bht]
\renewcommand{\arraystretch}{1.3}
\caption{The simulated mean energy deposit per interacting proton, relative beam intensity, and relative dose in the irradiated samples. All results refer to a vial placed at the beam centre.}
\label{tab:DoseInfo}
\centering
\begin{tabular}{cccc}
\hline
\hline
Sample Location & Mean Energy Deposit   & Relative beam intensity   & Relative dose          \\
\hline
Front HDPE holder       & 9.37 MeV              & 1                         & 1                      \\
Rear HDPE holder        & 9.75 MeV              & 0.957                     & 0.996                  \\
Online measurement      & 9.99 MeV              & 0.652                     & 0.695                 \\
\hline
\hline
\end{tabular}
\end{table*}

\subsubsection{Simulation Model}
We have modeled our irradiation experiments using a Geant4 v10.0 \cite{Agostinelli2003etal} simulation of the measurements. The calculated dose information for the irradiated samples is presented in table \ref{tab:DoseInfo}. The doses delivered to the front and rear vials were identical to within 0.5\%. Figure \ref{fig:Verticies} plots the original beam position of protons that interact in the rear vial. The relatively sharp image of the liquid volume suggests the beam divergence upon entry into the rear vial is on the order of 1 mm. Since the beam spreading was small compared to the length scale of the variations in the incident beam intensity (figure \ref{subfig:BeamProfile}), the incident beam profile was used to estimate the relative dose at the sample locations for the rear sample holder.  We attribute the excess events above the scintillating volume in figure \ref{fig:Verticies} to delta rays that travel through the air above the scintillator.

\begin{figure}[tbh]
\centering
\includegraphics[width=\columnwidth]{RearVialVerticies_fine_withLiquidOutline.pdf}
\caption{The initial beam location of simulated protons that deposit energy in the scintillator encapsulated by the central rear vial. The shaded region denotes the location of the scintillator. }
\label{fig:Verticies}
\end{figure}

\subsection{Irradiation B}
A second irradiation was carried out to investigate the effect of low doses upon the scintillator light yield. A beam of 201 MeV protons was used to deliver a dose of approximately 20 Gy and 50 Gy in two irradiations to an irradiation geometry similar to figure \ref{subfig:BeamSetup}. A 21.5 cm thick HDPE block was used to slow the incident proton beam, allowing the vials in the rear holder to receive a greater dose from stopping protons. Fricke dosimetry measurements of the front and rear vials suggest that the dose ratio in the rear relative to the front was 1.6. A beam profile similar to figure \ref{subfig:BeamProfile} was obtained.

\subsection{Light yield measurements}
Scintillation measurements were carried out using a Beckman LS6500 liquid scintillation spectrometer. A
$^{137}$Cs source was placed inside the instrument's lead shielding near the sample chamber to excite the scintillator. The Compton electron spectrum due to this source was measured repeatedly for all samples both prior to and after the irradiation. All irradiated scintillating samples were dispensed into unirradiated vials and compared with an unirradiated control.
Some instrumental artifacts were observed in the measured scintillation spectrum. These fluctuations in the histogrammed pulse height data exhibited a bin-to-bin correlation in a quasi-periodic fashion (figure \ref{fig:Spec}). The amplitude of these systematic variations did not vary strongly with the number of events in the histogram, so that acquiring for longer periods could minimise this effect. The precise cause of these artifacts is not known. The effect of the instrumental artifacts is taken into account in the assigned systematic uncertainty.

The change in the scintillation light yield upon irradiation was assessed relative to the control sample using an algorithm that compared the sample's Compton electron spectrum with that of the control.
A multiplicative factor was used to scale the sample spectrum along the channel axis. The histogram binning was correctly accounted for in the scaling operation. The change in the measured light yield was determined as the scaling factor that minimised the $\chi ^{2}$ between the measured and control spectra in a region above a channel threshold and below the bin in the control spectrum where the number of counts fell below 20. The optimal channel threshold was determined empirically for each sample. The operation of this algorithm is illustrated in figure \ref{fig:Spec}.

\begin{figure}[!t]
\centering
\includegraphics[width=\columnwidth]{Spec_PureLS_DemonstrateAlgo_Legend_rev1p1.pdf}
\caption{An example of the measured Compton electron spectrum for the control (black) and sample of liquid scintillator irradiated by $\approx$800 Gy of 201 MeV protons (red). The green trace shows the irradiated spectrum scaled by the amount of light yield loss as determined by the algorithm described in the text. The non-statistical fluctuations are an artifact of the liquid scintillation spectrometer readout.}
\label{fig:Spec}
\end{figure}

To account for the possibility that our algorithm's multiplicative factor is a biased estimator of the light yield, a Geant4 v10.0 simulation was developed to synthesize a similar Compton electron spectrum to that measured in our detector.
An isotropic source of 662 keV gamma rays was placed 2.5 cm from a scintillating volume identical in size to that used in our measurements. A 2.5 cm thick cylindrical lead shield  with an inner diameter of 7.5 cm surrounded the scintillating volume and source.
The simulated spectrum of energy deposits in the scintillator was convolved with a Gaussian of sufficient width to give an empirical match for the measured liquid scintillator and WbLS spectra.
Systematic errors in our light yield algorithm were investigated by comparing the algorithm's estimate of the scaling factor with the known scaling factor. This analysis shows that the algorithm underestimates the amount of light yield loss for both the pure liquid scintillator and the WbLS, with the lower amplitude WbLS spectrum having a larger error. The magnitude of this error was also found to increase as the reduction in the light yield increased. The required correction was less than 0.5\% of the light yield for all measured light yield losses, and we have corrected our reported results for this effect.

\setlength{\tabcolsep}{0pt}
\begin{figure*}[bt]
\centering
\begin{tabular}{cc}
    \subfloat[Liquid Scintillator] {\label{subfig:PureLSAbs} \includegraphics[width=0.49\textwidth,clip=true,trim=0.5cm 6cm 1.15cm 6cm]{Absorbance_LS_grid_BetterFont_BoldLines.pdf}}
    & \subfloat[Water-based Liquid Scintillator] {\label{subfig:WbLSAbs} \includegraphics[width=0.49\textwidth,clip=true,trim=0.5cm 6cm 1.15cm 6cm]{Absorbance_5pcWbLS_grid_BetterFont_boldLines.pdf}}   \\
\end{tabular}
\caption{The optical absorbance spectra for unirradiated samples (black), and samples irradiated with protons (blue and red). The average doses received by samples in the same irradiation holder are given.}
\label{fig:AbsSpec}
\end{figure*}
\setlength{\tabcolsep}{11pt}

\subsection{Optical Characterisation}
Optical absorbance measurements of the irradiated samples were carried out using a Shimadzu UV1800 spectrophotometer. We observed an increase in the optical attenuation coefficient of the irradiated scintillators that increased with increasing dose delivered to the scintillator (figure \ref{fig:AbsSpec}).

\section{Results}
The reduction of the light yield at the proton doses received by the samples is given in figure \ref{fig:LightYields}. The samples were measured approximately 24 hours following irradiation to allow short lived activation products to decay. No increase in the count rate was observed relative to the pre-irradiation measurements, suggesting that long lived activation products were created at negligible levels. All samples from Irradiation B exhibited no significant change in the light yield relative to the control samples. Their average dose and relative light yield is given in the figure, with the dose uncertainty covering the range of doses received by these samples (20 Gy to 90 Gy). The light yield uncertainty for all samples is the combination of the measurement standard deviation
and an estimated relative uncertainty on the bias correction of 50\%. A small loss of light yield was observed at approximately 320 Gy for both WbLS and liquid scintillator samples. For doses of approximately 800 Gy, a clear degradation of the scintillation yield is evident in both samples. This damage level agrees with the online measurements within their uncertainty, suggesting that there was no measurable recovery of scintillator in the time between the irradiation and the measurement. The scintillation light yield was measured repeatedly for approximately one week for all samples following irradiation, and no recovery or further deterioration was observed.

\begin{figure}[tbh]
\centering
\includegraphics[width=\columnwidth,clip=true,trim=0.5cm 6cm 0.5cm 6cm]{RelativeLY_All_grid_reanalysed_withLegend.pdf}
\caption{The measured light yield in the liquid scintillator (red) and Water-based Liquid Scintillator (blue) following proton irradiation, relative to an unirradiated control.}
\label{fig:LightYields}
\end{figure}

The optical absorbance spectra (figure \ref{fig:AbsSpec}) of both samples exhibit an increase with dose at wavelength regions relevant for the detection of scintillation photons. The reduction in the measured light yield of this increased optical absorption can be quantified using an effective absorbance due to the radiation damage, $\alpha_{eff}$:
\begin{equation}
f_{LY} = e^{-\bar{x}\alpha_{eff}}
\end{equation}
where $f_{LY}$ is the factor by which the light yield is reduced, and $\bar{x}$ is the mean optical path length of the emitted photons within the scintillator. The effective absorbance is calculated by weighting the excess damage-induced absorbance ($\alpha(\lambda)$, the difference between the irradiated and unirradiated traces in figure \ref{fig:AbsSpec}) with the normalised scintillator emission spectrum ($S(\lambda)$) and photomultiplier quantum efficiency ($QE(\lambda)$):
\begin{equation}
\alpha_{eff} = \int \alpha(\lambda)S(\lambda)QE(\lambda) d\lambda
\end{equation}
The amount of the observed light yield change that is due to the damage-induced optical absorbance is given in table \ref{tab:OptAbsEffect}, approximating $\bar{x}$ to the radius of the sample vial. The additional optical absorbance induced in the samples is insufficient to fully account for the observed change in the scintillation light yield, which suggests that the damage is also evident as a reduction in the scintillation photon production. This may be due to the radiation-induced formation of chemical quenching agents, or modification of the solvent/fluorophore molecules involved in the scintillation process \cite{Chong1997, Mesquita2001}. We performed Fourier transform infrared spectroscopy to identify chemical changes in the samples. No significant changes were observed. Given that the doses in Refs \cite{Chong1997, Mesquita2001} were orders of magnitude larger than the maximum dose of this work, this is perhaps unsurprising.

\begin{table*}[tbh]
\renewcommand{\arraystretch}{1.3}
\caption{The effect of damage-induced optical absorbance upon the measured light yield.}
\label{tab:OptAbsEffect}
\centering
\begin{tabular}{x{3.4cm}cx{2.5cm}x{4.5cm}}
\hline
\hline
Material                        & Mean Dose          & Measured light yield loss & Light yield loss expected from optical absorbance \\
\hline
Liquid Scintillator             & 327 $\pm$ 10 Gy    & 0.7 $\pm$ 0.4 \%        & 0.1 $\pm$ 0.2\%                                            \\
Liquid Scintillator             & 797 $\pm$ 24 Gy    & 1.7 $\pm$ 0.3 \%        & 0.3 $\pm$ 0.2\%                                            \\
Water-based Liquid Scintillator & 342 $\pm$ 10 Gy    & 0.8 $\pm$ 0.5 \%        & 0.0 $\pm$ 0.2\%                                           \\
Water-based Liquid Scintillator & 834 $\pm$ 25 Gy    & 1.3 $\pm$ 0.5 \%        & 0.7 $\pm$ 0.2\%                                            \\
\hline
\hline
\end{tabular}
\end{table*}

\section{Discussion and Conclusion}
No damage to the scintillation light yield was evident for proton doses up to 90 Gy. Given that a typical patient dose in a proton therapy treatment fraction is 2 Gy, this suggests that a liquid scintillator-based proton therapy dose-verification system could service at least 45 patients with no change in performance.

We will now consider the performance of such a QA instrument in more detail, given the observed changes at 800 Gy of proton dose. The typical annual number of patients served by a single proton therapy treatment room is approximately 300. Therefore, a QA instrument that validates the treatment plan of every patient may receive approximately 600 Gy in a year of operation, assuming that the instrument is used to validate a single 2 Gy treatment fraction for each patient. While our study has shown that the scintillator is damaged by this level of dose; an important consideration is that the entire volume of the scintillating phantom is not irradiated during the QA measurement. Indeed, for a phantom volume of 30$^{3}$ cm$^{3}$, and a typical lateral target size of 9 cm$^{2}$ and target depth of 20 cm, the irradiated volume is approximately 1\% of the total phantom volume. We assume a flat dose profile delivering uniform dose to this volume, which is a conservative approximation. Therefore, for complete mixing of the irradiated scintillator with the larger volume, the light yield damage is diluted by a factor of 100 relative to our measurements. In our measurement geometry, 600 Gy of proton dose to a scintillator results in an apparent light yield loss of approximately 1.3\%. Half of this loss is assumed to be due to damage-induced optical attenuation increases, which corresponds to $\alpha_{eff} = 0.0083$ cm$^{-1}$. Therefore, the irradiated phantom volume would exhibit a scintillation light production reduction of $ 1.3\% / (2 \times 100) = 0.007\%$. For a mean optical path length of 15 cm, the light absorption in the scintillator volume would decrease the detected number of photons by $1 - e^{-15*0.0083/100} = 0.1\%$. Therefore, the amount of light yield reduction due to radiation damage in a QA phantom will be approximately 0.1\% over a year of full-time clinical use. It is probable that this level of systematic variation of a dose verification instrument is acceptable for clinical operation given the magnitude of other uncertainties in the dose measurement \cite{ICRU78ch8}.

\begin{acknowledgments}
The authors would like to thank Mike Sivertz, Adam Rusek, and Chiara La Tessa at the NASA Space Radiation Laboratory, as well as Russ Burns, Rich Sautkulis, and Jim Jardine for their assistance with this study.
This research was funded by a Technology Maturation Award from the Office of Technology Commercialization and Partnerships at Brookhaven National Laboratory.
\end{acknowledgments}

\bibliographystyle{unsrt}