\section{model system}\label{s2}
We start from the Hamiltonian considering the coupling of a 2D TI and a $d$-wave superconductor, with $H =H_{\mathrm{TI}}+H_{\mathrm{SC}}+H_{\mathrm{I}}$.
$H_{TI}$ describes a 2D TI in the square lattice~\cite{re200,re201},
\begin{equation}
\begin{aligned}
H_{TI} &=\sum_{{\bf k}}C_{{\bf k}}^\dagger (h_\mathbf{k}\sigma_3s_0+2\lambda_{0}\sin k_x\sigma_1s_3\\&+2\lambda_0\sin k_y\sigma_2 s_0)C_{{\bf k}},
\end{aligned}
\end{equation}
with $h_\mathbf{k}=h_0-2t(\cos k_x+\cos k_y)$. $C^\dagger_{\bf k}$ is a forth order wave vector with $C^\dagger_{\bf k}=(c^\dagger_{{\bf k}1\uparrow},c^\dagger_{{\bf k}2\uparrow},c^\dagger_{{\bf k}1\downarrow},
c^\dagger_{{\bf k}2\downarrow})$. The subscripts $1,2$ and $\uparrow$, $\downarrow$ represent
the orbital and spin indices, respectively. $s_i$ and $\sigma_i$ are identity matrix $(i=0)$ or Pauli matrices $(i=1,2,3)$ in the spin and orbital spaces, respectively.


$H_{SC}$ describes a $d$-wave superconductor,
\begin{equation}
H_{SC}=\sum_{{\bf k}\sigma}\varepsilon_{\bf k}d^\dagger_{{\bf k}\sigma}d_{{\bf k}\sigma}+\sum_{{\bf k}}\Delta_{\bf k}(d^\dagger_{{\bf k}\uparrow}d^\dagger_{{-\bf k}\downarrow}+h.c.),
\end{equation}
with $\varepsilon_{\bf k}=-2t(\cos k_x+\cos k_y)-\mu$, and $\Delta_{\bf k}=2\Delta_0 (\cos k_x-\cos k_y)$.

$H_I$ is the interlayer single particle hopping term, with
\begin{equation}
H_{I}=-t_\perp \sum_{{\bf k}\alpha\sigma}(c_{{\bf k}\alpha\sigma}^\dagger d_{{\bf k}\sigma}+h.c.).
\end{equation}





To study the edge states and the Majorana corner states numerically,
let us transform
the Hamiltonian to the real space
by preforming a Fourier transformation. Then the Hamiltonian is rewritten as,
\begin{equation}
\begin{aligned}
H_{TI} &=-t\sum_{{\bf i}\alpha}(C_{{\bf i}}^\dagger \sigma_3s_0   C_{{\bf i}+\alpha}+h.c.)+h_0\sum_{{\bf i}\alpha}C_{{\bf i}}^\dagger \sigma_3s_0   C_{{\bf i}}
\\&-\lambda_0\sum_{{\bf i}}(iC_{{\bf i}}^\dagger\sigma_1s_3C_{{\bf i}+\hat{x}}-iC_{{\bf i}}^\dagger\sigma_2s_0C_{{\bf i}+\hat{y}}+h.c.),
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
H_{sc} &=-t\sum_{{\bf i}\alpha\sigma}(d_{{\bf i}\sigma}^\dagger d_{{\bf i}+\alpha,\sigma}+h.c.)-\mu\sum_{{\bf i}\sigma}d_{{\bf i}\sigma}^\dagger d_{{\bf i}\sigma}
\\&+\sum_{\langle{\bf ij}\rangle}(\Delta_{\bf ij}d_{{\bf i}\uparrow}^\dagger d_{{\bf j}\downarrow}^\dagger+h.c.),
\end{aligned}
\end{equation}
and
\begin{equation}
H_{I}=-t_\perp \sum_{{\bf i}\alpha\sigma}(c_{{\bf i}\alpha\sigma}^\dagger d_{{\bf i}\sigma}+h.c.),
\end{equation}
For the $d$-wave pairing, the site ${\bf j}$ is the nearest-neighbor site to the site ${\bf i}$, with
$\Delta_{\bf ij}=\pm\Delta_0$ ($\pm$ depends on $\langle{\bf ij}\rangle$ along the $x$ direction or $y$ direction).

The total Hamiltonian can be rewritten as the matrix form.
In the momentum space, the matrix is a $12\times 12$ matrix $\hat{M}_{\bf k}$ with
$H=\sum_{\bf k}\Psi^{\dagger}_{\bf k}\hat{M_{\bf k}}\Psi_{\bf k}$. The vector $\Psi^{\dagger}_{\bf k}$ is expressed as,
\begin{equation}
\Psi^\dagger_{\bf k}=(C^\dagger_{\bf k},C_{-{\bf k}},d^\dagger_{{\bf k}\uparrow},d^\dagger_{{\bf k}\downarrow},d_{-{\bf k}\uparrow},d_{-{\bf k}\downarrow}).
\end{equation}
In the real space, the Hamiltonian matrix $\hat{M}$ is a $N\times N$ lattice ($H=\Psi^{\dagger}\hat{M}\Psi$),
with $N=8N_1+4N_2$ ($N_1$ and $N_2$ are the number of sites in the 2D TI layer and the $d$-wave
superconducting layer, respectively). The vector $\Psi^{\dagger}$ is expressed as,
\begin{eqnarray}
\Psi^\dagger =(C^\dagger_1,C_{1},\cdots,C^\dagger_{N_1},C_{N_1},d^\dagger_{1\uparrow},\cdots,d_{N_2\downarrow}).
\end{eqnarray}
Here the vector $C^\dagger_{i}$ is the Fourier transformation of the vector $C^\dagger_{\bf k}$

Diagonalizing the Hamiltonian matrix, we obtain the retarded Green's Function matrix $\hat G$, with the elements being expressed as,
\begin{equation}
G_{ij}(E)=\sum_n\frac{u_{in}u^{*}_{jn}}{E-E_n+i\Gamma}.
\end{equation}
Here $u_{in}$ and $E_n$ are the eigen-vectors and eigen-values of the Hamiltonian matrix, respectively.

In the momentum space,
the spectral function of the 2D TI layer is calculated through,
\begin{equation}
A({\bf k},E)=-\frac{1}{\pi}\sum^4_{p=1}\mathrm{Im} G_{pp}({\bf k},E).
\end{equation}

The proximity induced pairing term for the orbital $\tau$ can be studied through the mean-field pairing order parameter, expressed as,
\begin{equation}
\Delta_\tau({\bf k})=\langle c^\dagger_{{\bf k}\tau\uparrow}c^\dagger_{-{\bf k}\tau\downarrow}\rangle=\sum_n{u^{*}_{\tau,n}({\bf k})u_{\tau+6,n}({\bf k})}f(E_n),
\end{equation}
where $f(x)$ is Fermi distribution function.




In the real space,
the effective pairing order parameter of the sites $i$ and $j$ for the orbital $\tau$ can be expressed as,
\begin{equation}
\Delta^\tau_{ij}=\sum_nu^{*}_{h(i),n}u_{h(j)+6,n}f(E_n),
\end{equation}
with $h(i)=\tau+8(i-1)$.

We define the site-dependent $d$-wave pairing magnitude for the orbital $\tau$, with
\begin{equation}
\Delta^\tau_i=\mid \Delta^\tau_{i,i+\hat{x}}+\Delta^\tau_{i,i-\hat{x}}
-\Delta^\tau_{i,i+\hat{y}}-\Delta^\tau_{i,i-\hat{y}}\mid.
\end{equation}
At the system edge, $\Delta^\tau_i$ is expressed as
\begin{equation}
\Delta^\tau_i=\mid 2(\Delta^\tau_{i,i+\hat{\alpha}}+\Delta^\tau_{i,i-\hat{\alpha}}) \mid  \qquad (\alpha=x,y).
\end{equation}

  \begin{figure}
\centering
  \includegraphics[width=1.8in]{fig1.eps}
\caption{(Color online) (a)The energy bands in the momentum space.
(b) The spectral function as functions of the energy and the momentum in the momentum space.
}\label{fig1}
\end{figure}

The local density of states (LDOS) of the site $i$ in the 2D TI layer can be calculated through the Green's function in the real space, with,
\begin{equation}
\rho_i(E)=-\frac{1}{\pi}\sum^4_{p=1}\mathrm{Im} G_{m+p,m+p}(E).
\end{equation}
with $m=8(i-1)$.

In the following presented results, the parameters are set as $t=1$, $\lambda_0=0.5$, $h_0=3$, $\mu=-0.3$, $\Delta_0=0.2$, $t_\perp=0.8$, and $\Gamma=0.01$. Our main results are not sensitive to the parameters we considered.






We first study the energy spectrum in the system bulk. The energy bands obtained by diagonalizing the $12\times12$ Hamiltonian in the momentum space are presented in Fig.~1(a). One striking feature is
that the $\mathcal{C}_4$ rotational symmetry of the energy bands is broken, namely, the quasiparticle bands are significantly different along the $(k_x,0)$ direction and the $(0,k_y)$ direction.
However, such a rotational symmetry is preserved by adding the $d$-wave pairing term directly into the 2D TI Hamiltonian and neglecting the inter-layer hopping~\cite{re100,RN4,re34}.
Experimentally, the energy bands can be detected through the spectral function. Here for the 2D TI layer, we have checked numerically that the spectra are fully gapped in the whole Brillouin zone.
The spectral function of the 2D TI layer as functions of the momentum and the energy is presented in Fig.~1(b).
Also, here the $\mathcal{C}_4$ symmetry is broken and such asymmetrical behavior may be detected by later experiments.

  \begin{figure}
\centering
  \includegraphics[width=3in]{fig2.eps}
\caption{(Color online) Numerical results considering a cylinder geometry. (a) The eigenvalues
of the Hamiltonian with the open boundary condition along the $x$-direction being considered.
(b) The spectral function at the $i_x=1$ boundary. (c) The spectral function at the $i_x=N_x$ boundary.
(d-f) Similar to panels (a-c) but for the open boundary condition along the $y$-direction.
}\label{fig2}
\end{figure}


We now study the edge states through considering the cylinder geometry with the open boundary condition along one direction and the periodic boundary condition along the other.
The energy bands with the open boundary condition along the $x$ direction
are presented in Fig.~2(a). The whole energy bands should include both the bulk states and the edge states. Note that, at these two boundaries, the energy bands are still fully gapped.
The gapped edge states can be studied more clearly through the spectral functions at the two system boundaries, which are presented in Figs.~2(b) and 2(c), respectively. As is displayed, the spectra at these two boundaries are exactly the same. An energy gap
around 0.05 is revealed. This gapped behavior originated from the proximity induced effective $d$-wave pairing term~\cite{re100,RN4,re34}.


We turn to discuss the edge states when the open boundary rotates to the $y$ direction.
The corresponding energy bands and the spectral functions at the system boundaries ($i_y=1$ and $i_y=N_y$) are presented in Figs.~2(d)-2(f).
As presented, both the energy bands and the spectral functions are different from those shown in Figs.~2(a)-2(c), indicating the $\mathcal{C}_4$ symmetry breaking. Interestingly, here from the energy band spectra [Fig.~2(a)], the edge states are gapless, significantly different from previous theoretical results~\cite{re100,RN4,re34}.
The spectral functions at the two boundaries are also different, namely, at the $i_y=1$ boundary, the spectrum is still gapped, with the gap magnitude around 0.08.
This value is even larger than those at the $i_x=1$ and $i_x=N_x$ boundaries. At the $i_y=N_y$ boundary, the edge state is gapless.
Thus for the cylinder system, the four boundaries are nonequivalent and the zero energy states appear at the $i_y=N_y$ boundary. Such asymmetric behavior may be detected by later experiments.

  \begin{figure}
\centering
  \includegraphics[width=3in]{fig3.eps}
\caption{(Color online) (a) Schematic illustration of a 2D TI being grown on a d-wave high-T$_c$ superconductor. (b) The eigenvalues of the Hamiltonian in the real space. (c) The intensity plot of the zero energy LDOS in the real space.
}\label{fig3}
\end{figure}

Now let us discuss the possible Majorana corner states.
We consider a 2D TI (with a finite-size lattice $40\times40$) being placed on a larger high-T$_c$ superconductor, as sketched in Fig.~3(a).
Considering open boundaries along both $x$ and $y$ directions, the Hamiltonian matrix is diagonalized. The corresponding eigenvalues are presented in Fig.~3(b).
Four zero energy eigenvalues are revealed. Note that there is no zero energy eigenvalue when periodic boundary conditions are taken.
Therefore, these four eigenvalues should locate at the system edges, relating to Majorana zero modes. For the topological superconducting system,
four zero energy eigenvalues should come from two zero energy physical quasiparticles, corresponding to two pairs of Majorana zero modes at the system boundaries.
The distributions of the two pairs of Majorana zero modes can be seen from the zero energy LDOS, as displayed in Fig.~3(c).
These two pairs locate at two lower corners of the system (one pair at each corner). As for the upper corners, no zero energy state exists.










The asymmetric results presented above can
be understood through exploring the pairing order parameter in the 2D TI layer. The order parameter for the orbital 1
 as a function of the momentum ${\bf k}$ is presented in Fig.~4(a). As is seen, here both the $\mathcal{C}_4$ symmetry and the inversion symmetry for the magnitude of the order parameter are broken. We can separate the whole pairing order parameter as the singlet channel $\Delta_e$ and the triplet channel $\Delta_o$ [$\Delta_1({\bf k})=\Delta_e({\bf k})+\Delta_o({\bf k})]$, with $\Delta_e({\bf k})=1/2[\Delta_1({\bf k})+\Delta_1(-{\bf k})]$ and $\Delta_o({\bf k})=1/2[\Delta_1({\bf k})-\Delta_1(-{\bf k})]$. The corresponding numerical results for these two channels are displayed in Figs.~4(b) and 4(c), respectively. For the singlet channel, the result is consistent with the $d_{x^2-y^2}$ pairing symmetry. For the triplet channel, a $p$-wave symmetry is revealed. The whole effective pairing symmetry should be $d$+$p$ wave. As a result, both the $\mathcal{C}_4$ symmetry and the inversion symmetry are broken for the spin-dependent spectral function for the orbital 1. As to the whole spectral function, the $\mathcal{C}_4$ symmetry is broken, while the inversion symmetry preserves due to the time reversal symmetry of the system.


  \begin{figure}
\centering
  \includegraphics[width=3.3in]{fig4.eps}
\caption{(Color online) Intensity plots of the pairing order parameter of the 2D TI layer.
(a) The whole pairing order parameter. (b) The pairing order parameter in the singlet channel part. (c) The pairing order parameter in the triplet channel part.
}\label{fig4}
\end{figure}

In the real space with open boundaries along both the $x$ direction and the $y$ direction, the site dependent $d$-wave pairing order parameters $\Delta_i$ for the orbital 1 are displayed in Fig.~5(a). The two dimensional cuts of $\Delta_i$ along the four boundaries are plotted in Fig.~5(b).
As is seen, at the $i_y=1$ boundary, the magnitude is relatively large.
For the $i_x=1$ and $i_x=N_x$ boundaries , the same magnitudes are revealed while they are smaller than that at the $i_y=1$ boundary.
As for the $i_y=N_y$ boundary, the magnitude of the induced pairing order parameter is rather small.
Therefore, with a cylinder geometry, the topological protected gapless edge states may survive at this boundary.
 These results are well consistent with the numerical results of the spectral function at the system boundaries [Fig.~2]. In the mean time, for the finite-size system with open boundaries, the Majorana corner states will generally emerge at the corners where the sign of the $d$-wave pairing term changes sign~\cite{re100,RN4,re34}. However,
since the effective $d$-wave order at the $i_y=N_y$ boundary is too small, the upper Majorana corner states will not come into being. Therefore, here the Majorana corner states only emerge at lower boundaries, as presented in Fig.~3.

\begin{figure}
\centering
\includegraphics[width=3in]{fig5.eps}
\caption{(Color online) (a) The intensity plot of the site-dependent $d$-wave pairing order parameter of the 2D TI layer.
(b) The replot of the $d$-wave pairing order parameter at the four boundaries.
}\label{fig5}
\end{figure}

The effective quasiparticle pairing in the 2D TI layer can also be explored through the anomalous Green's function~\cite{supp}.
In the momentum space, our analytical results of the anomalous Green's function verify that the singlet channel and the triplet channel pairing indeed coexist and leading to the $\mathcal{C}_4$ symmetry breaking. Numerically, the imaginary parts of the anomalous Green's function are consistent with the pairing order parameter shown in Fig.~4.




The asymmetric behavior at the system boundaries can also be well understood through the anomalous Green's function and the effective Hamiltonian at the system boundaries~\cite{supp}.
 Without coupling to the superconductor, the 2D TI is gapless at the boundaries, with the linear quasiparticle dispersion crossing the Fermi energy.
 The effective Hamiltonian at boundaries should include the intraorbital hopping term and the interorbital hopping term.
 When coupling to a superconductor, both terms contribute to the effective pairing at the boundaries.
 For the $i_y=N_y$ boundary, the interorbital hopping constant is opposite to the intraorbital one. These two contributions cancel out. Therefore, the effective pairing magnitude at this boundary is rather small.
For the $i_y=1$ boundary, these two contributions add up directly. The pairing magnitude at this boundary is relatively large. For the $i_x=1$ or $i_x=N_x$ boundary,
the interorbital hopping constant is imaginary.
An effective $d$-wave pairing term is also induced, while its magnitude is relatively smaller
than that of the $i_y=1$ boundary. The analytical results of the anomalous Green's function at these four boundaries are consistent with the above discussions.

At last, we would like to make several remarks. First, previously the proximity effect has been widely used to artificially create topological superconductors. Very recently, based on the proximity effect, several proposals have also been proposed to realize the second order topological superconductor in various hybrid systems. At this stage, studying the proximity effect in a more strict way, especially for the possible higher order topological superconducting system, is timely and of broad interest. Secondly, for the finite-size system, the Majorana bound state only emerge at lower corners. The zero energy degeneracy is reduced to half of those obtained from previous phenomenological theoretical results. The reduced degeneracy may make the Majorana bound states be more controllable. Thirdly, Our main results can be well understood. In the momentum space, the proximity induced pairings include both the $d$-wave component and the $p$-wave component, leading to the breaking of $\mathcal{C}_4$ symmetry.
At system boundaries in the real space, the proximity induced pairing terms are mainly contributed by the edge states of the TI.
The edge states are contributed by the intraorbital channel and the interorbital channel. The coexistence of these two channels may strengthen or weaken the induced pairing gap, leading to the interesting asymmetric behaviours.
At last, it is needed to pinpoint that the present system is different from both the first order topological superconductor and the higher order one, namely,
the system can be gapless for a certain one-dimensional boundary while it is fully gapped for other boundaries. In the meantime, the Majorana bound states emerge at partial corners of the finite system. It is also interesting that the zero energy states may shift from the upper boundary to the lower corners when the boundary condition along the $x$-direction changes.


