text:
# Tunneling conductance of normal metal / 𝑑_{𝑥²-𝑦²}-wave superconductor junctions in the presence of broken time reversal symmetry states near interfaces

## I Introduction

Identifying of the pairing symmetry in high-$`T_c`$ superconductors is important to clarify the mechanism of the superconductivity. A great deal of experimental and theoretical studies have revealed that the superconducting pair potential has $`d_{x^2y^2}`$-wave symmetry in the bulk state . Since the pair potential of the $`d_{x^2y^2}`$-wave superconductor is anisotropic, the amplitude of the pair potential near the surface or interface is significantly reduced. This suppression of the pair potential causes a very interesting phenomena, $`i`$.$`e`$. the formation of Andreev bound state (ABS) at the Fermi energy (zero-energy) near a specularly reflecting surface  when the angle between the lobe direction of the $`d_{x^2y^2}`$-wave pair potential and the normal to the interface is nonzero. This state is originated from the interference effect in the effective pair potential of the $`d_{x^2y^2}`$-wave symmetry through the reflection at the surface or interface. The ABS manifests itself as a sharp peak in the middle of the tunneling conductance spectra, the so-called zero-bias conductance peak (ZBCP),  and the consistency between theory and experiments has been checked in details .

On the other hand, there still remains a controversial issue; formation of a broken time reversal symmetry state (BTRSS) at low temperature due to the mixing of subdominant $`s`$-, or $`d_{xy}`$-wave component  as the imaginary part of the pair potential to the predominant $`d_{x^2y^2}`$-wave component . Theoretical studies based on the quasiclassical approximation  and several lattice models  reported the presence of the induced subdominant pair potential near the surface which breaks the time reversal symmetry . The resulting surface density of states with the BTRSS shows the splitting of the zero-energy peak  and the corresponding tunneling conductance shows the ZBCP splitting. It has also been clarified that the magnitude of the splitting depends on the induced subdominant pair potential near the surface. Actually, Covington $`et`$ $`al.`$ , Krupke $`et`$ $`al.`$ , and Sharoni $`et`$ $`al.`$  reports the ZBCP splitting at low temperature and they ascribed the origin of the ZBCP splitting to the above BTRSS . However, at the same time, there are many experiments which do not show the ZBCP splitting, and in these experiments, ZBCP survives even at low temperatures . Although there are several pre-existing theories which discuss the BTRSS , almost all these theories treat the semi-infinite $`d_{x^2y^2}`$-wave superconductor and the relevance to the actual experiments of tunneling spectroscopy has not been fully clarified yet. At this stage, it is an important problem to clarify the stability and the possible observability of the formation of BTRSS, $`i.e.`$, the ZBCP splitting, for various condition of the junctions. There are several factors which determine the magnitude of the splitting of ZBCP; (i) inter-electron potential which induces the subdominant pair potential near the interface, (ii) transmission probability of the particles at the interface, (iii) orientation of the junctions, $`i.e.`$, the angle between the normal to the interface and the crystal axis of $`d_{x^2y^2}`$-wave superconductor, (iv) finite temperature effect, (v) roughness of the interface , and (vi) impurity concentration in the superconductor .

Although there are several factors which determine the magnitude of the splitting of ZBCP, in this paper, we concentrate on (i), (ii), (iii), and (iv). First, we calculate the spatial variation of the pair potential in the normal metal / $`d_{x^2y^2}`$ -wave superconductor junction (N/D junction) in the presence of subdominant pair potential using a quasiclassical formalism for various conditions of the junctions. Using these results, we calculate tunneling conductance.

The organization of this paper is as follows. In Sec. II, a formulation to calculate the spatial dependence of the pair potential, and the tunneling conductance is presented. In Sec. III, results of the numerical calculations for both $`d_{x^2y^2}`$+i$`s`$\- and $`d_{x^2y^2}`$+i$`d_{xy}`$-wave states are discussed in detail. Section IV is devoted to conclusions and future problems.

## II Theoretical formulation

We study N/D junctions in the clean limit, where the normal metal is located at $`x<0`$ and the $`d_{x^2y^2}`$-wave superconductor extends elsewhere. For the simplicity, two dimensional system is assumed and the $`x`$-axis is taken perpendicular to the flat interface located at $`x=0`$. When quasiparticles are in the $`xy`$-plane, a transmitted electron-like quasiparticle and hole-like quasiparticle feel different effective pair potentials $`\mathrm{\Delta }(\varphi _+)`$ and $`\mathrm{\Delta }(\varphi _{})`$, with $`\varphi _+=\varphi `$ and $`\varphi _{}=\pi \varphi `$. Here, $`\varphi `$ is the azimuthal angle in the $`xy`$-plane given by $`(k_x+\mathrm{i}k_y)/|𝒌|=\mathrm{e}^{\mathrm{i}\varphi }`$. The barrier potential at the interface has a $`\delta `$-functional form $`H\delta (x)`$, where $`\delta (x)`$ and $`H`$ are the $`\delta `$-function and its amplitude, respectively. A cylindrical Fermi surface is assumed and the magnitude of the Fermi momentum and the effective mass are chosen to be equal both in the normal metal and in the superconductor.

The quasiclassical Green’s function method  developed by Ashida $`et`$ $`al.`$  is used in order to determine the spatial variation of the pair potential self-consistently. In the following, we briefly summarize this quasiclassical method we used. We start with the Bogoliubov-de Gennes (BdG) equation for unconventional spin-singlet superconductors ,

$`E_n\stackrel{~}{u}_n(𝒓)`$ $`=`$ $`H_0\stackrel{~}{u}_n(𝒓)+{\displaystyle 𝑑𝒓^{}\mathrm{\Delta }(𝒓,𝒓^{})\stackrel{~}{v}_n(𝒓^{})},`$ (1)
$`E_n\stackrel{~}{v}_n(𝒓)`$ $`=`$ $`H_0\stackrel{~}{v}_n(𝒓)+{\displaystyle 𝑑𝒓^{}\mathrm{\Delta }^{}(𝒓,𝒓^{})\stackrel{~}{u}_n(𝒓^{})},`$ (2)
$`H_0`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}^2\mu ,`$ (3)
where $`\mu `$ is the chemical potential, while $`\stackrel{~}{u}_n(𝒓)`$ and $`\stackrel{~}{v}_n(𝒓)`$ denote the electron-like and hole-like components of the wave function,

$`\stackrel{~}{\mathrm{\Psi }}_n(𝒓)`$ $`=`$ $`\left(\begin{array}{c}\stackrel{~}{u}_n(𝒓)\\ \stackrel{~}{v}_n(𝒓)\end{array}\right),`$ (6)
$``$ $`\left(\begin{array}{c}u_n(\widehat{𝒌},𝒓)\\ v_n(\widehat{𝒌},𝒓)\end{array}\right)e^{\mathrm{i}k_\mathrm{F}r}=\mathrm{\Psi }_n(\widehat{𝒌},𝒓)e^{\mathrm{i}k_\mathrm{F}r}.`$ (9)
Here the quantities $`\widehat{𝒌}`$ and $`𝒓`$ stand for the unit vector of the wave number of the Cooper pair which is fixed on the Fermi surface $`(\widehat{𝒌}=𝒌_\mathrm{F}/|𝒌_\mathrm{F}|)`$, and the position of the center of mass of Cooper pair, respectively. After applying the quasiclassical approximation, the BdG equation is reduced to the Andreev equation ,

$`E_n\mathrm{\Psi }_n(\widehat{𝒌},𝒓)`$ $`=`$ $`\left[\mathrm{i}\mathrm{}v_\mathrm{F}𝒌+\widehat{\mathrm{\Delta }}(\widehat{𝒌},𝒓)\right]\widehat{\tau }_3\mathrm{\Psi }_n(\widehat{𝒌},𝒓),`$ (10)
$`\widehat{\mathrm{\Delta }}(\widehat{𝒌},𝒓)`$ $`=`$ $`\left(\begin{array}{cc}0& \mathrm{\Delta }(\widehat{𝒌},𝒓)\\ \mathrm{\Delta }^{}(\widehat{𝒌},𝒓)& 0\end{array}\right),`$ (13)
where $`v_\mathrm{F}`$ and $`\widehat{\tau }_i(i=1,2,3)`$ stand for Fermi velocity and Pauli matrices, respectively. The wave function $`\mathrm{\Psi }_n(\widehat{𝒌},𝒓)`$ is obtained by neglecting the rapidly oscillating plane-wave part following the quasiclassical approximation . The $`\widehat{𝒌}`$ dependence of $`\mathrm{\Delta }(\widehat{𝒌},𝒓)`$ represents the symmetry of the pair potential.

Now, we consider the case where a specularly reflecting surface or interface runs along the $`y`$-direction. In this case, the pair potential depends only on $`x`$ since the system is homogeneous along the $`y`$-direction. It is convenient to introduce the following directional notation,

$`\mathrm{\Psi }_n(\widehat{𝒌},𝒓)`$ $`=`$ $`\mathrm{\Phi }_n^{(+)}(\varphi _+,x)e^{\mathrm{i}|k_{\mathrm{F}x}|x}+\mathrm{\Phi }_n^{()}(\varphi _{},x)e^{\mathrm{i}|k_{\mathrm{F}x}|x},`$ (14)
$`\mathrm{\Phi }_n^{(\alpha )}(\varphi _\alpha ,x)`$ $`=`$ $`\left(\begin{array}{c}u_n^{(\alpha )}(\varphi _\alpha ,x)\\ v_n^{(\alpha )}(\varphi _\alpha ,x)\end{array}\right).`$ (17)
Here $`\pm `$ represents the sign of the $`x`$ component of the Fermi wave number $`k_{\mathrm{F}x}`$ and $`\alpha (\beta )=\pm `$. We define a Green’s function $`G_{\alpha \beta }(\varphi ,x,x^{})`$ and a quasiclassical Green’s function $`g_{\alpha \beta }(\varphi ,x)`$,

$`G_{\alpha \beta }(\varphi ,x,x^{})`$ $`=`$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{\mathrm{\Phi }_n^{(\alpha )}(\varphi _\alpha ,x)\mathrm{\Phi }_n^{(\beta )}(\varphi _\beta ,x)}{\mathrm{i}\omega _mE_n}},`$ (18)
$`g_{\alpha \beta }(\varphi ,x)\pm \mathrm{i}(\widehat{\gamma }_3)_{\alpha \beta }`$ $`=`$ $`2\mathrm{}|v_{\mathrm{F}x}|\widehat{\tau }_3G_{\alpha \beta }(\varphi ,x\pm 0,x).`$ (19)
In the above, $`\widehat{\gamma }_3`$ is the Pauli matrix in the directional space and $`v_{\mathrm{F}x}`$ is the $`x`$-component of the Fermi velocity. The quasiclassical Green’s function $`g_{\alpha \beta }(\varphi ,x)`$ obeys the Eilenberger equation,

$`\mathrm{i}|v_{\mathrm{F}x}|{\displaystyle \frac{}{x}}g_{\alpha \beta }(\varphi ,x)=`$ $``$ $`\alpha \left[\mathrm{i}\omega _m\widehat{\tau }_3+\widehat{\mathrm{\Delta }}(\varphi _\alpha ,x)\right]g_{\alpha \beta }(\varphi ,x)`$ (20)
$`+`$ $`\beta `$ $`g_{\alpha \beta }(\varphi ,x)\left[\mathrm{i}\omega _m\widehat{\tau }_3+\widehat{\mathrm{\Delta }}(\varphi _\beta ,x)\right],`$ (21)
$`\widehat{\mathrm{\Delta }}(\varphi _\alpha ,x)`$ $`=`$ $`\left(\begin{array}{cc}0& \mathrm{\Delta }(\varphi _\alpha ,x)\\ \mathrm{\Delta }^{}(\varphi _\alpha ,x)& 0\end{array}\right),`$ (24)
where $`\omega _m`$ is the Matsubara frequency. The quasiclassical Green’s function can be written by the following evolution operator $`U_\alpha (\varphi _\alpha ,x,x^{})`$ as

$`g_{\alpha \beta }(\varphi ,x)=U_\alpha (\varphi _\alpha ,x,x^{})g_{\alpha \beta }(\varphi ,x^{})U_\beta ^1(\varphi _\beta ,x,x^{}),`$ (25)
where $`U_\alpha (\varphi _\alpha ,x,x^{})`$ satisfies the Andreev equation

$`\mathrm{i}\mathrm{}|v_{\mathrm{F}x}|{\displaystyle \frac{}{x}}U_\alpha (\varphi _\alpha ,x,x^{})=\alpha \left[\mathrm{i}\omega _m\widehat{\tau }_3+\widehat{\mathrm{\Delta }}(\varphi _\alpha ,x)\right]U_\alpha (\varphi _\alpha ,x,x^{}),`$ (26)
with $`U_\alpha (\varphi _\alpha ,x,x)=1`$.

Considering a semi-infinite N/D junction geometry, the pair potential in the superconductor side approaches to the bulk value $`\mathrm{\Delta }(\varphi _\alpha ,\mathrm{})`$ at sufficiently large $`x`$. Hence, the evolution operator can be divided into a growing part and a decaying part:

$`U_\alpha (\varphi _\alpha ,x,x^{})`$ $`=`$ $`\mathrm{\Lambda }_\alpha ^{(+)}(\varphi _\alpha ,x,x^{})e^{\kappa _\alpha (xx^{})}+\mathrm{\Lambda }_\alpha ^{()}(\varphi _\alpha ,x,x^{})e^{\kappa _\alpha (xx^{})},`$ (27)
$`\mathrm{\Lambda }_\alpha ^{(+)}(\varphi _\alpha ,x,x^{})`$ $`=`$ $`{\displaystyle \frac{1}{W_\alpha }}\mathrm{\Phi }_n^{(+)}(\varphi _\alpha ,x){}_{}{}^{\mathrm{T}}\mathrm{\Phi }_{n}^{()}(\varphi _\alpha ,x^{})\widehat{\tau }_2,`$ (28)
$`\mathrm{\Lambda }_\alpha ^{()}(\varphi _\alpha ,x,x^{})`$ $`=`$ $`{\displaystyle \frac{1}{W_\alpha }}\mathrm{\Phi }_n^{()}(\varphi _\alpha ,x){}_{}{}^{\mathrm{T}}\mathrm{\Phi }_{n}^{(+)}(\varphi _\alpha ,x^{})\widehat{\tau }_2,`$ (29)
where

$`\kappa _\alpha `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_\alpha }{|v_{\mathrm{F}x}|}},\mathrm{\Omega }_\alpha =\sqrt{\omega _m^2+|\mathrm{\Delta }(\varphi _\alpha ,\mathrm{})|^2},`$ (30)
$`W_\alpha `$ $`=`$ $`{}_{}{}^{\mathrm{T}}\mathrm{\Phi }_{n}^{(+)}(\varphi _\alpha ,x)\widehat{\tau }_2\mathrm{\Phi }_n^{()}(\varphi _\alpha ,x)={}_{}{}^{\mathrm{T}}\mathrm{\Phi }_{n}^{()}(\varphi _\alpha ,x)\widehat{\tau }_2\mathrm{\Phi }_n^{(+)}(\varphi _\alpha ,x)`$ (31)
$`=`$ $`\text{constant}.`$ (32)
In the above, $`{}_{}{}^{\mathrm{T}}\mathrm{\Phi }_{n}^{(\alpha )}(\varphi _\alpha ,x)`$ denotes the transposition of $`\mathrm{\Phi }_n^{(\alpha )}(\varphi _\alpha ,x)`$.

Retaining the most divergent term in semi-infinite limit, we find the quasiclassical Green’s function $`\widehat{g}_{\alpha \alpha }(\varphi ,x)`$ in the superconductor side given by ,

$`\widehat{g}_{\alpha \alpha }(\varphi _\alpha ,x)=\mathrm{i}\left({\displaystyle \frac{2\widehat{A}_{\mathrm{S}\alpha }(x)}{\mathrm{Tr}[\widehat{A}_{\mathrm{S}\alpha }(x)]}}1\right),`$ (33)
where

$`\widehat{A}_{\mathrm{S}+}(x)`$ $`=`$ $`\mathrm{\Lambda }_+^{()}(\varphi _+,x,L)\mathrm{\Lambda }_{}^{(+)}(\varphi _{},L,0)\widehat{R}_\mathrm{N}U_+(\varphi _+,0,x)e^{\kappa x}`$ (34)
$`=`$ $`\widehat{\lambda }_{\mathrm{S}+}(x,0)\widehat{R}_\mathrm{N}\stackrel{~}{U}_+(\varphi _+,0,x),`$ (35)
$`\widehat{\lambda }_{\mathrm{S}+}(x,0)`$ $``$ $`\left(\begin{array}{c}u_n^{()}(\varphi _+,x)\\ v_n^{()}(\varphi _+,x)\end{array}\right)\left(\begin{array}{cc}u_n^{()}(\varphi _{},0)& v_n^{()}(\varphi _{},0)\end{array}\right)\widehat{\tau _2}.`$ (39)
In the above, the matrix $`\widehat{R}_\mathrm{N}`$ represents resistance at the interface which is given by

$`\widehat{R}_\mathrm{N}\left(\begin{array}{cc}1& 0\\ 0& 1\sigma _\mathrm{N}(\varphi )\end{array}\right),\sigma _\mathrm{N}(\varphi )={\displaystyle \frac{4\mathrm{cos}\varphi ^2}{4\mathrm{cos}\varphi ^2+Z^2}},`$ (42)
where $`\sigma _\mathrm{N}(\varphi )`$ stands for tunneling conductance when the system is in the normal state  and $`Z`$ is the effective barrier height at the interface with $`Z=2mH/(\mathrm{}^2k_\mathrm{F})`$. In order to obtain the quantity $`\stackrel{~}{U}_+(\varphi _+,0,x)`$ in eq.(34), we rewrite $`\widehat{A}_{\mathrm{S}+}(x)`$ as

$`\widehat{A}_{\mathrm{S}+}(x)=\left(\begin{array}{c}u_n^{()}(\varphi _+,x)\\ v_n^{()}(\varphi _+,x)\end{array}\right)\left(\begin{array}{cc}X_+(x)& Y_+(x)\end{array}\right)\widehat{\tau _2},`$ (46)
where $`𝒟_\alpha (x)=(\mathrm{i})v_n^{()}(\varphi _\alpha ,x)/u_n^{()}(\varphi _\alpha ,x)`$ and $`_+(x)=\mathrm{i}X_+(x)/Y_+(x)`$ obey the following Riccati type equations

$`\mathrm{}|v_{\mathrm{F}x}|{\displaystyle \frac{}{x}}𝒟_\alpha (x)`$ $`=`$ $`\alpha \left[2\omega _m𝒟_\alpha (x)+\mathrm{\Delta }(\varphi _\alpha ,x)𝒟_\alpha ^2(x)\mathrm{\Delta }^{}(\varphi _\alpha ,x)\right],`$ (47)
$`\mathrm{}|v_{\mathrm{F}x}|{\displaystyle \frac{}{x}}_+(x)`$ $`=`$ $`2\omega _m_+(x)+\mathrm{\Delta }^{}(\varphi _+,x)_+^2(x)\mathrm{\Delta }(\varphi _+,x).`$ (48)
We can write the quasiclassical Green’s function in a compact form ,

$`\widehat{g}_{++}(\varphi _+,x)`$ $`=\mathrm{i}\left[{\displaystyle \frac{2}{1𝒟_+(x)_+(x)}}\left(\begin{array}{cc}1& \mathrm{i}_+(x)\\ \mathrm{i}𝒟_+(x)& 𝒟_+(x)_+(x)\end{array}\right)1\right].`$ (51)
Initial conditions of these equations are

$`𝒟_\alpha (\mathrm{})={\displaystyle \frac{\mathrm{\Delta }^{}(\varphi _\alpha ,\mathrm{})}{\omega _m+\alpha \mathrm{\Omega }_\alpha }},_+(0)={\displaystyle \frac{1\sigma _\mathrm{N}(\varphi )}{D_{}(0)}}.`$ (52)
The pair potential is given by

$`\mathrm{\Delta }(\varphi ,x)={\displaystyle \underset{0m<\omega _c/2\pi T}{}}{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\pi /2}^{\pi /2}}𝑑\varphi ^{}{\displaystyle \underset{\alpha }{}}V(\varphi ,\varphi _\alpha ^{})[\widehat{g}_{\alpha \alpha }(\varphi _\alpha ^{},x)]_{12},`$ (53)
with $`\widehat{g}_{}(\varphi _{},x)=\widehat{g}_{++}(\varphi _+,x)`$, where $`\omega _c`$ is the cutoff energy and $`[\widehat{g}_{\alpha \alpha }(\varphi _\alpha ,x)]_{12}`$ means the 12 element of $`\widehat{g}_{\alpha \alpha }(\varphi _\alpha ,x)`$. Here $`V(\varphi ,\varphi _\alpha )`$ is the effective inter-electron potential of the Cooper pair. In our numerical calculations, new $`\mathrm{\Delta }(\varphi _\alpha ,x)`$ and $`\widehat{g}_{\alpha \alpha }(\varphi _\alpha ,x)`$ are obtained using eqs.(47)-(51) and eq. (53). We repeat this iteration process until the sufficient convergence is obtained.

Next, we calculate the tunneling conductance spectra based on the self-consistently determined pair potential. The resulting normalized tunneling conductance $`\sigma _\mathrm{T}(eV)`$ with the bias voltage $`V`$ is given by

$`\sigma _\mathrm{T}(eV)`$ $`=`$ $`{\displaystyle \frac{{\displaystyle _{\pi /2}^{\pi /2}}𝑑\varphi {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑E\sigma _\mathrm{N}(\varphi )\sigma _\mathrm{S}(E,\varphi )\mathrm{sech}^2\left({\displaystyle \frac{EeV}{2k_\mathrm{B}T}}\right)\mathrm{cos}\varphi }{{\displaystyle _{\pi /2}^{\pi /2}}𝑑\varphi {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑E\sigma _\mathrm{N}(\varphi )\mathrm{sech}^2\left({\displaystyle \frac{EeV}{2k_\mathrm{B}T}}\right)\mathrm{cos}\varphi }},`$ (54)
$`\sigma _\mathrm{S}(E,\varphi )`$ $`=`$ $`{\displaystyle \frac{1+\sigma _\mathrm{N}(\varphi )|\mathrm{\Gamma }_{\mathrm{S}+}(E,\varphi _+,0)|^2+[\sigma _\mathrm{N}(\varphi )1]|\mathrm{\Gamma }_{\mathrm{S}+}(E,\varphi _+,0)|^2|\mathrm{\Gamma }_\mathrm{S}(E,\varphi _{},0)|^2}{|1+[\sigma _\mathrm{N}(\varphi )1]\mathrm{\Gamma }_{\mathrm{S}+}(E,\varphi _+,0)\mathrm{\Gamma }_\mathrm{S}(E,\varphi _{},0)|^2}}.`$ (55)
In the above, $`\mathrm{\Gamma }_{\mathrm{S}\alpha }(E,\varphi _\alpha ,x)`$ is obtained by solving following equations,

$`\mathrm{i}\mathrm{}|v_{\mathrm{F}x}|{\displaystyle \frac{}{x}}\mathrm{\Gamma }_{\mathrm{S}+}(E,\varphi _+,x)=2E\mathrm{\Gamma }_{\mathrm{S}+}(E,\varphi _+,x)\mathrm{\Delta }(\varphi _+,x)\mathrm{\Gamma }_{\mathrm{S}+}^2(E,\varphi _+,x)\mathrm{\Delta }^{}(\varphi _+,x),`$ (56)
$`\mathrm{i}\mathrm{}|v_{\mathrm{F}x}|{\displaystyle \frac{}{x}}\mathrm{\Gamma }_\mathrm{S}(E,\varphi _{},x)=2E\mathrm{\Gamma }_\mathrm{S}(E,\varphi _{},x)\mathrm{\Delta }^{}(\varphi _{},x)\mathrm{\Gamma }_\mathrm{S}^2(E,\varphi _{},x)\mathrm{\Delta }(\varphi _{},x).`$ (57)
In the following, almost calculations are performed on the temperature $`T/T_c=0.05`$, where $`T_c`$ is the critical temperature of the bulk $`d_{x^2y^2}`$-wave superconductor.

## III Broken time reversal symmetry state near an interface of a $`d_{x^2y^2}`$-wave superconductor

In this section, the spatial dependence of the self-consistently determined pair potential and the corresponding tunneling conductance are presented for the $`d_{x^2y^2}`$-wave superconducting state. In the middle of the $`d_{x^2y^2}`$-wave superconductor, the pair potential is given by $`\mathrm{\Delta }(\varphi _\alpha ,\mathrm{})=\mathrm{\Delta }_0\mathrm{cos}[2(\varphi \theta )]`$, where $`\theta `$ is the angle between normal to the interface and the lobe direction of the $`d_{x^2y^2}`$-wave pair potential, $`i.e.`$, the angle between the $`x`$-axis and the crystal $`a`$-axis of the $`d_{x^2y^2}`$-wave. In this paper, we choose various $`\theta `$ $`(0\theta \pi /4)`$ by changing the magnitude of $`Z`$ and $`T_s`$ ($`T_{d_{xy}}`$).

### A $`d_{x^2y^2}`$\+ i$`s`$-wave state

In this subsection, we show the spatial dependence of the pair potential and the resulting tunneling conductance of the $`d_{x^2y^2}`$+i$`s`$-wave state realized near the interface of the N/D junction. The spatial dependence of the pair potential is expressed as

$`\mathrm{\Delta }(\varphi ,x)=\mathrm{\Delta }_d(x)\mathrm{cos}[2(\varphi \theta )]+\mathrm{\Delta }_s(x),`$ (58)
where $`\mathrm{\Delta }_d(x)`$ and $`\mathrm{\Delta }_s(x)`$ correspond to the amplitude of the $`d_{x^2y^2}`$-wave and $`s`$-wave superconducting states, respectively. The attractive potential $`V(\varphi ,\varphi ^{})`$ is given by

$`V(\varphi ,\varphi ^{})=2V_d\mathrm{cos}[2(\varphi \theta )]\mathrm{cos}[2(\varphi ^{}\theta )]+V_s`$ (59)
where $`V_d`$ and $`V_s`$ denote the attractive potential of predominant $`d_{x^2y^2}`$-wave and subdominant $`s`$-wave, respectively, and they are given as

$`V_d`$ $`=`$ $`{\displaystyle \frac{2\pi k_\mathrm{B}T}{\mathrm{log}{\displaystyle \frac{T}{T_c}}+{\displaystyle \underset{0m<\omega _c/2\pi T}{}}{\displaystyle \frac{1}{m+1/2}}}},`$ (60)
$`V_s`$ $`=`$ $`{\displaystyle \frac{2\pi k_\mathrm{B}T}{\mathrm{log}{\displaystyle \frac{T}{T_s}}+{\displaystyle \underset{0m<\omega _c/2\pi T}{}}{\displaystyle \frac{1}{m+1/2}}}}.`$ (61)
Here, $`T_s`$ denotes the transition temperature of $`s`$-wave component of the pair potential without predominant $`d_{x^2y^2}`$-wave component. For $`\theta =0`$ or $`\theta =\pi /4`$, only $`\text{Re}[\mathrm{\Delta }_d(x)]`$ and $`\text{Im}[\mathrm{\Delta }_s(x)]`$ are nonzero. The spatial dependence of the pair potentials $`\text{Re}[\mathrm{\Delta }_d(x)]`$ and $`\text{Im}[\mathrm{\Delta }_s(x)]`$ is plotted in Fig. 1(a) for various $`\theta `$ with $`T_s/T_d=0.2`$ and $`Z=3`$. The $`x`$-axis of Fig. 1(a) is normalized by $`\xi _0=\mathrm{}v_\mathrm{F}/\mathrm{\Delta }_0`$ which is the coherence length of the superconductor. For $`\theta =0`$ \[(100) surface\], since $`\text{Im}[\mathrm{\Delta }_s(x)]=0`$ is satisfied, the time reversal symmetry is not broken and the amplitude of $`\text{Re}[\mathrm{\Delta }_d(x)]`$ is not suppressed at the interface. By changing the angle $`\theta `$ from zero, the magnitude of $`\text{Re}[\mathrm{\Delta }_d(x)]`$ is reduced near the interface, while $`\text{Im}[\mathrm{\Delta }_s(x)]`$ is induced at the interface . The suppression of $`\text{Re}[\mathrm{\Delta }_d(x)]`$ is originated from a depairing effect that the effective pair potentials $`\mathrm{\Delta }(\varphi _+,0)`$ and $`\mathrm{\Delta }(\varphi _{},0)`$ have reversed contribution to the pairing interaction for certain range of $`\varphi `$ for $`\theta 0`$. When $`\text{Re}[\mathrm{\Delta }_d(x)]`$ is suppressed at the interface, the quasiparticle forms the ABS near the interface at zero-energy . The ABS is unstable with the introduction of $`s`$-wave attractive potential, then $`\text{Im}[\mathrm{\Delta }_s(x)]`$ is induced at the interface . The magnitude of $`\text{Im}[\mathrm{\Delta }_s(x)]`$ becomes maximum at $`\theta =\pi /4`$, where the above suppression effect is most significant. In Fig. 1(b), $`\text{Re}[\mathrm{\Delta }_d(x)]`$ and $`\text{Im}[\mathrm{\Delta }_s(x)]`$ are plotted for various $`Z`$ with $`T_s/T_d=0.2`$ and $`\theta =\pi /4`$. Even if the BTRSS becomes to be most stable at $`\theta =\pi /4`$, when the height of barrier is small, the magnitude of the subdominant imaginary component of $`\mathrm{\Delta }_s(x)`$ is not induced at all. The induced imaginary component of $`\mathrm{\Delta }_s(x)`$ is enhanced with the increase of $`Z`$.

The spatial dependence of the pair potentials near the interface with the intermediate angle ($`\theta =\pi /6`$) is shown in Fig. 2 for various height of barrier. In such a case, both $`\text{Im}[\mathrm{\Delta }_d(x)]`$ and $`\text{Re}[\mathrm{\Delta }_s(x)]`$ becomes nonzero and the spatial dependence is much more complex as compared to that for $`\theta =0`$ or $`\theta =\pi /4`$. The amplitudes of $`\text{Im}[\mathrm{\Delta }_d(x)]`$, $`\text{Re}[\mathrm{\Delta }_s(x)]`$, and $`\text{Im}[\mathrm{\Delta }_s(x)]`$ are enhanced for larger magnitude of $`Z`$, where the suppression of $`\text{Re}[\mathrm{\Delta }_d(x)]`$ is significant. However, the amplitudes of $`\text{Im}[\mathrm{\Delta }_d(x)]`$ and $`\text{Re}[\mathrm{\Delta }_s(x)]`$ are one order smaller than that of $`\text{Im}[\mathrm{\Delta }_s(x)]`$.

Next, we look at the magnitude of subdominant components of the pair potential at the interface, $`\text{Im}[\mathrm{\Delta }_s(0)]`$, $`\text{Im}[\mathrm{\Delta }_d(0)]`$, and $`\text{Re}[\mathrm{\Delta }_s(0)]`$, for various $`T_s`$, $`Z`$, and $`\theta `$. As shown in Fig. 3(a), the magnitude of $`\text{Im}[\mathrm{\Delta }_s(0)]`$ increases monotonically with $`T_s`$ for fixed $`\theta `$ and $`Z`$, and it is enhanced for larger magnitude of $`Z`$. In other words, the amplitude of $`\text{Im}[\mathrm{\Delta }_s(0)]`$ is sensitive to the transmission probability of the junctions. In Fig. 3(b), $`\text{Im}[\mathrm{\Delta }_s(0)]`$ is plotted as a function of $`\theta `$ for sufficiently larger magnitude of $`Z(=5.0)`$. For $`\theta =0`$, $`i.e.`$, junction with (100) interface, the magnitude of $`\text{Im}[\mathrm{\Delta }_s(0)]`$ is negligibly small near the interface even at the larger magnitude of $`T_s`$. The magnitude of $`\text{Im}[\mathrm{\Delta }_s(0)]`$ is a monotonically increasing function with the increase of $`\theta `$ and has a maximum at $`\theta =\pi /4`$. As seen from Figs. 3(c) and 3(d), both the magnitude of $`\text{Im}[\mathrm{\Delta }_d(0)]`$ and $`\text{Re}[\mathrm{\Delta }_s(0)]`$ is enhanced and has a maximum at a certain $`\theta `$. In the intermediate $`\theta `$, $`i.e.`$, $`\theta 0`$ or $`\theta \pi /4`$, the magnitude of $`\mathrm{\Delta }(\varphi _+,x)`$ and that of $`\mathrm{\Delta }(\varphi _{},x)`$ does not coincide any more, the interference with the quasiparticle and the pair potential becomes complex. Then, not only $`\text{Im}[\mathrm{\Delta }_s(0)]`$ but also $`\text{Re}[\mathrm{\Delta }_s(0)]`$ and $`\text{Im}[\mathrm{\Delta }_d(0)]`$ become nonzero.

Using self-consistently determined pair potentials, let us look at the normalized tunneling conductance $`\sigma _\mathrm{T}(eV)`$. In order to clarify the temperature $`T`$ dependence of $`\sigma _\mathrm{T}(eV)`$, we choose $`T=0`$ in the left panels of Fig. 4 and $`T=0.05T_c`$ in the right panels. Only for $`\theta =0`$, line shape of $`\sigma _\mathrm{T}(eV)`$ is similar to that of the bulk density of states of $`d_{x^2y^2}`$-wave superconducting state. In other cases, $`\sigma _\mathrm{T}(eV)`$ has a zero bias enhanced line shape. As clarified in previous literatures , when $`\theta `$ deviates from zero, since the injected and reflected quasiparticles have a chance to feel the sign change of the pair potentials, zero-energy ABS is formed at the interface. This zero-energy ABS causes the ZBCP when the magnitude of $`T_s`$ is small. With the increase of the magnitude of $`T_s`$, the zero energy ABS is unstable and $`s`$-wave subdominant component is induced which breaks time reversal symmetry and it blocks the motion of the quasiparticles. Then, the energy levels of bound state shift from zero and the local density of states has a zero-energy peak splitting. The resulting $`\sigma _\mathrm{T}(eV)`$ has a ZBCP splitting as shown in Fig. 4(b). However, with the increase of $`T`$, the slight splitting of ZBCP fades out due to smearing effect by finite temperature and the resulting $`\sigma _\mathrm{T}(eV)`$ has a rather broad ZBCP \[see Fig. 4(c) and Fig. 4(d)\].

Next, we concentrate on how $`\sigma _\mathrm{T}(eV)`$ is influenced by the transmission probability of the junctions, $`i.e.`$, the magnitude of $`Z`$. In Fig. 5, $`\sigma _\mathrm{T}(eV)`$ with $`\theta =\pi /4`$ is plotted for $`T=0`$ (left panels) and $`T=0.05T_c`$ (right panels). For the junctions with high transmissivity, $`\sigma _\mathrm{T}(eV)`$ has a ZBCP, \[see Fig. 5(a)\] and the magnitude of $`\sigma _\mathrm{T}(0)`$ is firstly enhanced with the increase of $`Z`$. In this case, the predominant $`d_{x^2y^2}`$-wave component only exists near the interface as shown in Fig. 1(b). However, with the increase of $`Z`$, $`\sigma _\mathrm{T}(eV)`$ starts to have a ZBCP splitting at a certain value of $`Z`$, where the magnitude of subdominant component $`\text{Im}[\mathrm{\Delta }_s(x)]`$ at the interface becomes the same order as that of the predominant component $`\text{Re}[\mathrm{\Delta }_d(x)]`$. For sufficiently larger magnitude of $`Z`$, $`\sigma _\mathrm{T}(eV)`$ has a ZBCP splitting, \[see Fig. 5(b) and 5(c),\] and the magnitude of $`\sigma _\mathrm{T}(0)`$ decreases with the increase of $`Z`$. However, the above obtained results are influenced by finite temperature effect. The right panels of Fig. 5 is shown for the tunneling conductance in $`T/T_c=0.05`$. The slight enhanced structure of $`\sigma _T(eV)`$ at $`eV=\pm \mathrm{\Delta }_0`$ in Fig. 5(a), 5(b), and 5(c) is invisible due to the smearing effect by finite temperature \[see Fig. 5(c), 5(d), and 5(e)\]. With the increase of the magnitude of $`Z`$, $`\sigma _\mathrm{T}(eV)`$ has a ZBCP with tiny dip even at $`Z=2.5`$, where the order of the amplitude of $`\text{Im}[\mathrm{\Delta }_s(0)]`$ is 0.2$`\mathrm{\Delta }_0`$. With the further increase of $`Z`$, $`\sigma _\mathrm{T}(eV)`$ has a ZBCP splitting \[see Fig. 5(f)\], however the degree of the splitting is significantly weakened as compared to the corresponding curves in Fig. 5(c).

Finally, we look at the relation between the position of the splitted peak and the magnitude of $`\text{Im}[\mathrm{\Delta }_s(0)]`$. In Fig. 6(a), $`\sigma _\mathrm{T}(eV)`$ is plotted for various $`T_s`$ with $`Z=5`$ and $`T=0`$. As shown in Fig. 3(b), the magnitude of the induced subdominant imaginary component of $`\text{Im}[\mathrm{\Delta }_s(0)]`$ is about $`0.16\mathrm{\Delta }_0`$, $`0.3\mathrm{\Delta }_0`$, and $`0.42\mathrm{\Delta }_0`$ for $`T_s/T_d=0.1`$, $`T_s/T_d=0.3`$ and $`T_s/T_d=0.5`$, respectively. The corresponding $`\sigma _\mathrm{T}(eV)`$ has a splitted peak locating at $`\pm 0.16\mathrm{\Delta }_0`$, $`\pm 0.3\mathrm{\Delta }_0`$, and $`\pm 0.42\mathrm{\Delta }_0`$, respectively. With the increase of $`T`$, the height of these peaks are drastically suppressed as shown in Fig. 6(b). Finally, we show how the line shape of $`\sigma _\mathrm{T}(eV)`$ changes at the temperature $`T=T_{\overline{s}}`$ ($`T_{\overline{s}}=0.12T_c`$) where $`\text{Im}[\mathrm{\Delta }_s(x)]`$ becomes nonzero. As seen from Fig. 6(c), the magnitude of $`\sigma _T(0)`$ is reduced with the introduction of $`\text{Im}[\mathrm{\Delta }_s(x)]`$.

At the end of this subsection, we can summarize that even in the presence of BTRSS, the resulting $`\sigma _\mathrm{T}(eV)`$ does not always have a clear ZBCP splitting due to the finite temperature effect when the magnitude of the transmission probability of the junctions is not low.

### B $`d_{x^2y^2}`$+i$`d_{xy}`$-wave state

In this subsection, we study spatial dependence of the pair potentials of the $`d_{x^2y^2}`$+i$`d_{xy}`$-wave state and the resulting tunneling conductance in N/D junctions. The pair potential is given by

$`\mathrm{\Delta }(\varphi ,x)=\mathrm{\Delta }_d(x)\mathrm{cos}[2(\varphi \theta )]+\mathrm{\Delta }_{d_{xy}}(x)\mathrm{sin}[2(\varphi \theta )],`$ (62)
where $`\mathrm{\Delta }_{d_{xy}}(x)`$ is an amplitude of the $`d_{xy}`$-wave superconducting state and a complex number. The attractive inter-electron potential $`V(\varphi ,\varphi ^{})`$ is given by

$`V(\varphi ,\varphi ^{})=2V_d\mathrm{cos}[2(\varphi \theta )]\mathrm{cos}[2(\varphi ^{}\theta )]+2V_{d_{xy}}\mathrm{sin}[2(\varphi \theta )]\mathrm{sin}[2(\varphi ^{}\theta )],`$ (63)
where $`V_d`$ and $`V_{d_{xy}}`$ stand for the attractive potential of predominant $`d_{x^2y^2}`$-wave and subdominant $`d_{xy}`$-wave, respectively, and they are given as

$`V_d`$ $`=`$ $`{\displaystyle \frac{2\pi k_\mathrm{B}T}{\mathrm{log}{\displaystyle \frac{T}{T_c}}+{\displaystyle \underset{0m<\omega _c/2\pi T}{}}{\displaystyle \frac{1}{m+1/2}}}},`$ (64)
$`V_{d_{xy}}`$ $`=`$ $`{\displaystyle \frac{2\pi k_\mathrm{B}T}{\mathrm{log}{\displaystyle \frac{T}{T_{d_{xy}}}}+{\displaystyle \underset{0m<\omega _c/2\pi T}{}}{\displaystyle \frac{1}{m+1/2}}}}.`$ (65)
The spatial dependence of the pair potentials with a finite transmissivity for $`Z=3.0`$ and $`\theta =\pi /4`$ is shown in Fig. 7(a). As in the case for $`d_{x^2y^2}`$+i$`s`$-wave state, the amplitude of $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(x)]`$ vanishes at $`\theta =0`$. For $`\theta =\pi /4`$, the suppression of the magnitude of $`\text{Re}[\mathrm{\Delta }_d(x)]`$ is most significant, while $`\mathrm{\Delta }_{d_{xy}}(x)`$ is induced at the interface. At this $`\theta `$, the $`d_{xy}`$-wave component is not affected by depairing effect seriously since the lobe direction of $`d_{xy}`$-wave pair potential is parallel or perpendicular to the interface as in the case for predominant $`d_{x^2y^2}`$-wave with $`\theta =0`$. As shown in Fig. 7(b), the amplitude of $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(x)]`$ is enhanced with the increase of $`Z`$.

The spatial dependencies of $`\text{Re}[\mathrm{\Delta }_d(x)]`$, $`\text{Im}[\mathrm{\Delta }_d(x)]`$, $`\text{Re}[\mathrm{\Delta }_s(x)]`$, and $`\text{Im}[\mathrm{\Delta }_s(x)]`$ are plotted in Figs. 8(a), 8(b), 8(c), and 8(d), respectively, with the intermediate $`\theta `$ ($`\theta =\pi /6`$) for various $`Z`$. In such a case, the spatial dependence is much more complex as compared to that for $`\theta =0`$ or $`\theta =\pi /4`$, and both $`\text{Im}[\mathrm{\Delta }_d(x)]`$ and $`\text{Re}[\mathrm{\Delta }_{d_{xy}}(x)]`$ become nonzero. The amplitudes of $`\text{Im}[\mathrm{\Delta }_d(x)]`$, $`\text{Re}[\mathrm{\Delta }_{d_{xy}}(x)]`$, and $`\text{Im}[\mathrm{\Delta }_s(x)]`$ are enhanced for larger magnitude of $`Z`$, where the suppression of $`\text{Re}[\mathrm{\Delta }_d(x)]`$ is significant. The remarkable feature is that the amplitude of $`\text{Im}[\mathrm{\Delta }_d(x)]`$ and $`\text{Re}[\mathrm{\Delta }_{d_{xy}}(x)]`$ can become the same order of $`\text{Im}[\mathrm{\Delta }_s(x)]`$ for larger $`Z`$.

Next, we look at the magnitude of subdominant components of the pair potential at the interface, $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(0)]`$, $`\text{Im}[\mathrm{\Delta }_d(0)]`$, and $`\text{Re}[\mathrm{\Delta }_{d_{xy}}(0)]`$ for various parameters $`T_{d_{xy}}`$, $`Z`$ and $`\theta `$ shown in Fig. 9. The magnitude of $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(0)]`$ increases monotonically with $`T_{d_{xy}}`$ for fixed $`\theta `$ and $`Z`$. It is also enhanced for larger magnitude of $`Z`$. Comparing the corresponding situation of $`d_{x^2y^2}+\mathrm{i}s`$-wave shown in Fig. 3(a), the magnitude of $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(0)]`$ is suppressed for the small magnitude of $`T_{d_{xy}}`$. In Fig. 9(b), $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(0)]`$ is plotted as a function of $`\theta `$ for sufficiently larger magnitude of $`Z(=5.0)`$. The magnitude of $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(0)]`$ is a monotonically increasing function with the increase of $`\theta `$ and has a maximum at $`\theta =\pi /4`$. The amplitudes of $`\text{Im}[\mathrm{\Delta }_d(0)]`$ and $`\text{Re}[\mathrm{\Delta }_{d_{xy}}(0)]`$ are negligibly small for $`\theta =0`$ and $`\theta =\pi /4`$, and has a maximum at a certain $`\theta `$ as shown in Figs. 9(c) and 9(d), respectively.

The resulting normalized tunneling conductance $`\sigma _\mathrm{T}(eV)`$ is plotted in Figs. 10(a) and 10(b) with $`T=0`$ and the corresponding quantities with $`T=0.05T_c`$ is plotted in Figs. 10(c) and 10(d) for various $`\theta `$. Only for $`\theta =0`$, line shape of $`\sigma _\mathrm{T}(eV)`$ is similar to that of the bulk density of states of $`d_{x^2y^2}`$-wave superconducting state. In other cases, $`\sigma _\mathrm{T}(eV)`$ has a zero-bias enhanced line shape due to the formation of ABS \[see Fig. 4(a)\]. The ABS causes the ZBCP when the magnitude of $`T_{d_{xy}}`$ is small. However, with the increase of the magnitude of $`T_{d_{xy}}`$, the zero-energy ABS is unstable and subdominant $`d_{xy}`$-wave component is induced. The energy levels of bound state shifts from zero and the local density of states has a zero-energy peak splitting. The resulting $`\sigma _\mathrm{T}(eV)`$ has a ZBCP splitting as shown in Fig. 10(b). As compared to the case of $`d_{x^2y^2}`$+i$`s`$ state, the magnitude of $`\sigma _\mathrm{T}(0)`$ is not reduced seriously since the effective $`d_{xy}`$-wave pair potential felt by quasiparticles is distributed from $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(0)]`$ to $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(0)]`$. With the increase of temperature, the ZBCP splitting is not visible any more. \[see Fig. 10(c) and Fig. 10(d)\].

Next, we concentrate on how $`\sigma _\mathrm{T}(eV)`$ is influenced by the transmissivity of the junctions, $`i.e.`$, the magnitude of $`Z`$. In Fig. 11, $`\sigma _\mathrm{T}(eV)`$ with $`\theta =\pi /4`$ is plotted for $`T=0`$ (left panels) and $`T=0.05T_c`$ (right panels). As shown in Fig. 11, $`\sigma _\mathrm{T}(eV)`$ with $`\theta =\pi /4`$ depends on the transmissivity of the junction crucially. For the junctions with high transmissivity, $`\sigma _\mathrm{T}(eV)`$ has a ZBCP \[see Fig. 11(a)\] and the magnitude of $`\sigma _\mathrm{T}(0)`$ is firstly enhanced with the increase of $`Z`$. In this case, only the predominant $`d_{x^2y^2}`$-wave component exists near the interface as shown in Fig. 1(b). However, with the increase of $`Z`$, $`\sigma _\mathrm{T}(eV)`$ starts to have a peak splitting at a certain value of $`Z`$, where the magnitude of the induced subdominant component at the interface $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(0)]`$ overlaps that of the predominant one $`\text{Re}[\mathrm{\Delta }_d(x)]`$. For sufficiently larger magnitude of $`Z`$, $`\sigma _\mathrm{T}(eV)`$ has a ZBCP splitting \[see Figs. 11(b) and 11(c)\] and the magnitude of $`\sigma _\mathrm{T}(0)`$ decreases with the increase of $`Z`$. The right panels of Fig. 11 are shown for the $`\sigma _\mathrm{T}(eV)`$ with finite temperature. Only for the junctions with the lowest transmissivity ($`Z=5`$), the subtle splitting of ZBCP appears in $`\sigma _\mathrm{T}(eV)`$. For $`T=0.05T_C`$, similar line shapes of $`\sigma _T(eV)`$ with smeared ZBCP splitting are obtained \[see Fig. 11(d), Fig. 11(e) and Fig. 11(f)\]. However, the slight enhanced structure of $`\sigma _T(eV)`$ at $`eV=\pm \mathrm{\Delta }_0`$ vanishes.

Finally, we look at the relation between the position of the splitted peak and the magnitude of $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(0)]`$. In Fig. 12(a), we show $`T_{d_{xy}}`$ dependence of the tunneling conductance spectra $`\sigma _\mathrm{T}(eV)`$ for $`\theta =0`$ with $`Z=5`$. As shown in Fig. 9(b), the magnitude of the induced subdominant imaginary component of $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(0)]`$ is about $`0.15\mathrm{\Delta }_0`$, $`0.35\mathrm{\Delta }_0`$, and $`0.52\mathrm{\Delta }_0`$ for $`T_{d_{xy}}/T_d=0.1`$, $`T_{d_{xy}}/T_d=0.3`$ and $`T_{d_{xy}}/T_d=0.5`$, respectively. The corresponding $`\sigma _\mathrm{T}(eV)`$ has a splitted peak locating at $`\pm 0.16\mathrm{\Delta }_0`$, $`\pm 0.3\mathrm{\Delta }_0`$, and $`\pm 0.42\mathrm{\Delta }_0`$, respectively. The width of the two splitted peaks depends on the temperature as shown in Fig. 12(b). Next, we show how the line shape of $`\sigma _\mathrm{T}(eV)`$ changes at the temperature $`T=T_{\overline{d_{xy}}}`$ \[$`T_{\overline{d_{xy}}}=0.07T_c`$\] where $`\text{Im}[\mathrm{\Delta }_{d_{xy}}(x)]`$ becomes nonzero. As seen from Fig. 12(c), the magnitude of $`\sigma _\mathrm{T}(0)`$ is reduced with the introduction of $`\text{Im}[\mathrm{\Delta }_s(x)]`$.

At the end of this section, we can summarize that even in the presence of BTRSS, the resulting $`\sigma _\mathrm{T}(eV)`$ does not always have a clear ZBCP splitting due to the finite temperature effect when the transmission probability of the particles at the interface is not low. In order to observe the ZBCP splitting which is one of the evidence supporting BTRSS, we must measure $`\sigma _\mathrm{T}(eV)`$ for the junctions with low transmissivity with (110) oriented interface ($`\theta =\pi /4`$) at low temperatures. In such a case, we can classify the the $`d_{x^2y^2}`$+i$`s`$\- and $`d_{x^2y^2}`$+i$`d_{xy}`$-wave state through the magnitude of $`\sigma _\mathrm{T}(0)`$ and the width of the ZBCP splitting for the junctions with low transmission probability.

## IV Conclusion

In this paper, spatial dependence of the pair potential in the N/D junctions is determined on the basis of the quasiclassical theory in the presence of subdominant component of the pair potential near the interface. We clarified the influence of the spatial variation of the pair potentials on the tunneling conductance spectra for various conditions of the junctions. We selected two kinds of subdominant components $`s`$\- and $`d_{xy}`$-wave which are induced as $`d_{x^2y^2}`$+i$`s`$\- and $`d_{x^2y^2}`$+i$`d_{xy}`$-wave state, respectively. The amplitude of $`\mathrm{\Delta }_s(x)`$ \[$`\mathrm{\Delta }_{d_{xy}}(x)`$\] is enhanced with the increase of the magnitude of $`Z`$, $`T_s[T_{d_{xy}}]`$, and $`\theta `$ ($`0<\theta <\pi /4`$). In the intermediate $`\theta `$, both $`\text{Im}[\mathrm{\Delta }_d(x)]`$ and $`\text{Re}[\mathrm{\Delta }_s(x)]`$ $`\{\text{Re}[\mathrm{\Delta }_{d_{xy}}]\}`$ becomes nonzero, although these magnitudes are small as compared to those of $`\text{Re}[\mathrm{\Delta }_d(x)]`$ and $`\text{Im}[\mathrm{\Delta }_s(x)]`$ $`\{\text{Im}[\mathrm{\Delta }_{d_{xy}}]\}`$. The resulting $`\sigma _\mathrm{T}(eV)`$ depends on $`Z`$, $`T_s`$ (or $`T_{d_{xy}}`$), and $`\theta `$. For fixed $`Z`$, the magnitude of $`\sigma _\mathrm{T}(0)`$ increases with the increase of $`\theta `$ $`(0<\theta <\pi /4)`$ for small magnitude of $`T_s/T_d`$ \[$`T_{d_{xy}}/T_d`$\]. While $`\sigma _\mathrm{T}(0)`$ increases and decreases again with the increase of $`\theta `$ $`(0<\theta <\pi /4)`$ due to the formation of $`d_{x^2y^2}`$+i$`s`$\- \[$`d_{x^2y^2}`$+i$`d_{xy}`$\]-wave state for sufficiently larger magnitude of $`T_s/T_d`$ \[$`T_{d_{xy}}/T_d`$\]. For fixed $`\theta `$ and $`T_s/T_d`$ \[$`T_{d_{xy}}/T_d`$\], $`\sigma _\mathrm{T}(eV)`$ has a ZBCP for small magnitude of $`Z`$, while it has a ZBCP splitting for sufficiently larger magnitude of $`Z`$. Using junctions with small transmissivity, we can distinguish $`d_{x^2y^2}`$+i$`s`$-wave state from $`d_{x^2y^2}`$+i$`d_{xy}`$-wave state since the magnitude of $`\sigma _\mathrm{T}(0)`$ for $`d_{x^2y^2}`$+i$`s`$-wave state is much more reduced as compared to that for $`d_{x^2y^2}`$+i$`d_{xy}`$-wave state as seen from Fig. 4 \[Fig. 6\] to Fig. 10 \[Fig. 12\]. By taking into account finite temperature effect, the degree of the ZBCP splitting is suppressed and the fine structure at $`eV\mathrm{\Delta }_0`$ in $`\sigma _\mathrm{T}(eV)`$ becomes invisible.

In the light of our theory, one of the reason for the absence of ZBCP splitting for many tunneling experiments may be due to their high transmissivity of the junctions (small magnitude of $`Z`$) and high temperatures. If we choose $`T_c`$ as 90K, unless the transmissivity of the junction is sufficiently low, it is difficult to observe ZBCP splitting clearly at 4.5K. In order to see the ZBCP splitting clearly, we must observe much more low temperature ($`T<4.5K`$) using a clean junctions with low transmissivity for $`\theta =\pi /4`$.

In this paper, we have chosen a free electron model with cylindrical Fermi surfaces. In order to compare actual tunneling conductance, we must calculate spatial dependence of the pair potential and the tunneling conductance by taking into account of the actual shape of Fermi surface. For this purpose, it is a promising way to perform the calculation based on the $`t`$-$`J`$ model with Gutzwiller approximation , where the doping dependence is naturally taken into account. Using this model, we succeeded to explain detailed line shape of the scanning tunneling conductance spectra around Zn impurity . As regards quasiparticle states near the interface, the pre-existing calculations are performed for infinite barrier limit . In order to compare actual tunneling experiments, we must calculate the N/D junctions based on the $`t`$-$`J`$ model using Gutzwiller approximation. In the $`t`$-$`J`$ model, the subdominant $`d_{xy}`$-wave component is hard to be realized  since a nearest neighbor attractive potential is taken into account. In order to discuss the $`d_{xy}`$-wave component, we must consider much more long range interaction including three site hopping term.

It is also revealed by recent theory that the BTRSS influences significantly on the Josephson effect . The previous theory assumes a sufficient larger magnitude of barrier . In order to compare the existing experiments, we must calculate Josephson current in the presence of BTRSS for arbitrary transmissivity of the junctions.

## Acknowledgments

First author (Y.T.) would like to thank K. Machida and M. Ichioka for useful and fruitful discussions. He acknowledges the financial support of Research Fellowships of Japan Society for the Promotion of Science (JSPS) for Young Scientists. This work is supported by the Core Research for Evolutional Science and Technology (CREST) of the Japan Science and Technology Corporation (JST). The computational aspect of this work has been performed at the facilities of the Supercomputer Center, Institute for Solid State Physics, University of Tokyo and the Computer Center.


format: text

text:
# Disk Accretion at 10 Myr: Results from the TW Hydrae Association

## 1. Introduction

The TW Hydrae association (TWA) has come to be identified as one of the nearest associations of young stars (de la Reza et al. 1989; Webb et al. 1999; see also relevant articles in this volume), and hence has garnered considerable attention in the past few years with the characterization of its population. Moreover, its age of $``$ 10 Myr (Webb et al. 1999) makes it an extremely interesting region with which to study the processes of planet formation, since initial studies have shown that 10 Myr should be the approximate timescale for the dissipation of optically thick circumstellar disks (Strom, Edwards, & Skrutskie 1993).

An important question related to the dissipation of disks is whether gas accretion onto the central star, and hence a reservoir of gas in the disk, can persist for up to 10 Myr, and whether any such older disks show signs of evolution from their younger cousins. The timescale for the dissipation of gas disks is crucial for the formation of giant planets; core accretion models require up to 10 Myr or more to form giant planets (e.g. Bodenheimer, Hubickyj, & Lissauer 2000). The presence of active accretion in any of the TWA members would thus indicate that indeed in some cases, optically thick gas disks do survive for up to 10 Myr. TW Hya itself is already a well-known classical T Tauri star, with attendant accretion disk (Rucinski & Krautter 1983). Extending the initial results presented in Muzerolle et al. (2000), we investigate whether other members of the association also exhibit accretion signatures, and use these signatures to examine evidence of disk evolution.

## 2. Observations

Optical echelle spectra of 11 TWA members were obtained with HIRES on Keck. The spectral resolution was $`8`$ km s-1, with a wavelength coverage of about 6310 - 8750 Å. This dataset is supplemented with Keck ESI spectra of 4 of the 5 remaining members listed in Webb et al. (1999), the 3 candidates from Sterzik et al. (1999), and 5 of the 8 new candidates identified by Zuckerman et al. (2001). The ESI spectra are of lower resolution (37.5 km s-1), but have a larger wavelength coverage of about 3900 - 10750 Å. Though the resolution is much lower than the echelle spectra, it is still sufficient to resolve line profiles formed in magnetically-mediated accretion flows (FWHM $`200`$ km s-1), where the gas is undergoing ballistic infall onto the star (e.g. Muzerolle, Calvet, & Hartmann 2001). Table 1 lists the objects observed, along with the H$`\alpha `$ equivalent widths and FWHM line widths (multiple values indicate measurements from both ESI and HIRES, respectively), and whether the \[OI\] $`\lambda `$6300 line appears in emission.

## 3. The Accretors of TWA

Only 3 of the 23 TWA members observed, TW Hya, Hen 3-600A, and TWA 14, show broad, asymmetric Balmer emission line profiles indicative of active accretion. All other objects exhibit narrow (FWHM $`<100`$ km s-1) and symmetric Balmer line profiles, typical of the chromospheric emission of weak T Tauri stars (a few of these objects are probable spectroscopic binaries, where the photospheric absorption lines are broadened or doubled, in which case the Balmer lines are slightly broader at FWHM $`100140`$ km s-1). Two of the three objects with Balmer profiles showing gas infall, TW Hya and Hen 3-600A, also show \[OI\] $`\lambda `$6300 emission, an indicator of mass outflow seen only in actively accreting young stars; none of the remaining stars observed with ESI show this emission. Furthermore, these two objects exhibit infrared excess emission indicative of circumstellar disks (Jayawardhana et al. 1999a,b)111No published infrared observations yet exist for TWA 14.. Finally, TW Hya and TWA 14 exhibit continuum veiling of their photospheric absorption features, with values of 2.4 and 0.2, respectively, at 4400 Å(because Hen 3-600A is a double-lined binary, reliable veiling determinations cannot be made without knowing the stellar properties of each component).

In order to quantify the accretion activity, we have modeled the H$`\alpha `$ line profiles using radiative transfer models of magnetospheric accretion. These models successfully reproduce the observed line emission in typical 1 Myr-old classical T Tauri stars in Taurus (Muzerolle, Calvet, & Hartmann 2001). Muzerolle et al. (2000) have already presented results for TW Hya and Hen 3-600A; see that paper for details. Model comparisons from that paper are shown in Figure 1, along with new results for TWA 14. Note that the stellar parameters for this object are not well-determined; we assume the same distance and age as TW Hya, and, given the similar spectral types, also the same mass and radius. In each of the three accreting TWA members, the model accounts very well for the observed profile shape, width, and strength, with surprisingly small mass accretion rates ($`\dot{M}=5\times 10^{10}M_{}yr^1`$ for TW Hya, and $`\dot{M}=5\times 10^{11}`$ for Hen 3-600A and TWA 14) compared to typical values measured in 1 Myr-old accretors ($`\dot{M}10^8`$). We note that the two HIRES and ESI observations of TW Hya and Hen 3-600A show relatively similar profile shapes and line strengths, with TW Hya showing somewhat stronger blueshifted absorption in the ESI spectrum.

The H$`\alpha `$ models can only provide constraints on the mass accretion rate to within about a factor of 5, however, because of uncertainties in the parameterized gas temperature and magnetospheric size (see Muzerolle, Calvet, & Hartmann 2001). A much more definitive measurement of $`\dot{M}`$ can be made from the luminosity of the UV continuum excess, which is produced in the accretion shock formed as accreting material falls onto the surface of the star. Muzerolle et al. (2000) used the detailed accretion shock calculations of Calvet & Gullbring (1998) to model the UV and optical spectral energy distribution of TW Hya, the only one of the TWA accretors with sufficient UV observations. The resulting value of the accretion rate, $`4\times 10^{10}M_{}yr^1`$, is in good agreement with the H$`\alpha `$ models. 222Batalha and Alencar, this volume, present extensive observations of the blue excess and line profile variations, deriving larger accretion rates for TW Hya. The reason for the discrepancy is unclear, but possibly due in part to intrinsic variability.

## 4. Accretion at 10 Myr

Due to the similarity of the Balmer emission line profiles of the TWA accretors to those in the well-studied younger objects in Taurus, and the success the accretion models have in explaining these observations, it is clear that the standard disk accretion scenario for low-mass young stars still applies in these older objects. Thus, some stars are able to maintain magnetically-mediated accretion flows, and associated optically thick accretion disks, for at least 10 Myr. The mass accretion rates in the 10 Myr-old disks are significantly lower than in typical 1 Myr-old objects, which indicates that the gas disk surface densities have decreased substantially over 10 Myr. Such a decrease is in fact predicted by models of viscous disk evolution (Hartmann et al. 1998).

Moreover, the fraction of TWA members which are still accreting, and hence possess optically thick gas disks, is only 13%. This is significantly lower than that found in the 1 Myr-old star forming regions such as Taurus and Orion, where disk fractions as determined from near-infrared excesses range from 60-80% (Kenyon & Hartmann 1995; Lada et al. 2000). Both the decreased accretion disk fraction and the very small mass accretion rates in the three TWA accretors lend substantial evidence for significant disk evolution occurring over 10 Myr. Recent results using near- and mid-infrared excesses as signatures of optically thick inner and outer disk regions also point to similar disk dissipation timescales (e.g., Jayawardhana, this volume).

The presence of the three accreting 10 Myr-old T Tauri stars further begs the question of why their disks have survived for so long when optically thick gas disks in the remaining TWA members have since dissipated. The presence of close companions does not seem to preclude long-lived accretion disks, since Hen 3-600A is a spectroscopic binary with accretion signatures, and also has a third companion (Hen 3-600B, not accreting) $`50`$ AU distant. As noted by Muzerolle et al. (2000), a circumbinary disk around the two components of Hen 3-600A should be truncated to $`20`$ AU; for the gas to have lasted for 10 Myr, the initial mass of the circumbinary disk must have been more massive than typical ($`>0.1M_{}`$). It is possible that planet formation has led to the dissipation of gas disks in the other TWA stars, some of which do exhibit mid-infrared emission from debris disks (Jaywardhana et al. 1999b), which may represent the detritus of the planet formation process. Planet formation may be ongoing in the disks of the three accretors; the spectral energy distribution of TW Hya indicates that dust coagulation or clearing has occurred within the innermost $`1`$ AU of its disk (see D’Alessio, this volume).

## References

Bodenheimer, P., Hubickyj, O., & Lissauer, J. J. 2000, Icarus, 143, 2

Calvet, N. & Gullbring, E. 1998, ApJ, 509, 802

de la Reza, R., Torres, C. A. O., Quast, G., Castilho, B. V., & Vieira, G. L. 1989, ApJ, 343, L61

Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385

Jayawardhana, R. et al. 1999a, ApJ, 520, L41

Jayawardhana, R. et al. 1999b, ApJ, 521, L129

Kenyon, S. J. & Hartmann, L. 1995, ApJS, 101, 117

Lada, C. J. et al. 2000, AJ, 120, 3162

Muzerolle, J., Calvet, N., Briceño, C., Hartmann, L., & Hillenbrand, L. 2000, ApJ, 535, L47

Muzerolle, J., Calvet, N., & Hartmann, L. 2001, ApJ, 550, 944

Rucinski, S. M., & Krautter, J. 1983, A&A, 121, 217

Sterzik, M. F., Alcalá, J. M., Covino, E., & Petr, M. G. 1999, A&A, 346, L41

Webb, R. A. et al. 1999, ApJ, 512, L63

Zuckerman, B., Webb, R. A., Schwartz, M., & Becklin, E. E. 2001, ApJ, 549, L233


text:
# Exchange in silicon-based quantum computer architectures

\[

## Abstract

The silicon-based quantum computer proposal has been one of the actively pursued ideas during the past three years. Here we calculate the donor electron exchange in silicon and germanium, and demonstrate an atomic-scale challenge for quantum computing in Si (and Ge), as the six (four) conduction band minima in Si (Ge) lead to inter-valley electronic interference, generating strong oscillations in the exchange splitting of two-donor two-electron states. Donor positioning with atomic scale precision within the unit cell thus becomes a decisive factor in determining the strength of the exchange coupling—a fundamental ingredient for two-qubit operations in a silicon-based quantum computer.

\]

Following the seminal proposal  by Kane there has been a great deal of activity  in efforts to develop a silicon-based quantum computer (QC) architecture. The basic ideas of the Kane proposal are simple and attractive: to use donor nuclear spins as quantum bits (qubits), and to utilize the vast infrastructure and technology associated with the Si industry to fabricate precisely controlled Si nanostructures, where exchange effects between electrons and nuclei in neighboring donor impurities (e.g. 31P in Si) could serve as the two-qubit gates, similar to the electron-spin-based QC proposal by Loss and DiVincenzo . The motivation for a Si quantum computer is obvious: Once the basic one-qubit and two-qubit operations have been demonstrated using donor impurities in Si nanostructures, computer chip fabrication technology associated with the existing and dominant Si industry will easily enable the scale-up of information processing involving large number of donor nuclear spin qubits. Indeed, one of the formidable stumbling blocks in developing working quantum computer hardwares has been the scale-up problem, as the demonstrated qubits in trapped ion and liquid state NMR techniques are not readily scalable in any significant manner .

A great deal of experimental work is currently being aimed at developing suitable qubits in Si nanostructures with precisely introduced dopant impurities, using both a “top-down” approach with ion-implantation, and a “bottom-up” approach with MBE growth and scanning tunneling microscopy . In the Si QC model , donor electrons act as shuttles between different nuclear spins. For two-qubit operations, which are required for a universal QC, both electron-electron exchange and electron-nucleus hyperfine interaction need to be precisely controlled. These are unquestionably formidable experimental problems. In the original proposal, Kane used the Herring-Flicker exchange formula  for two hydrogenic centers to obtain an order of magnitude estimate of the electron exchange among donors in Si . However, as he also pointed out, donor exchange in Si is not hydrogenic.

In this Letter we show that exchange effects in proposed donor nuclei based Si QC architectures are actually very subtle due to quantum interference effects inherent in the complicated Si band structure. In particular, special care, going far beyond what is being currently attempted in the fabrication of Si QC, will be required in controlling the surface gates crucial to QC operations. We also provide results for Ge which seems to have certain advantages over Si as a candidate for dopant exchange based quantum computation. Unfortunately, Ge is not such an attractive practical alternative since there is no well-established Ge technology to take advantage of.

The study of shallow impurities in Si and Ge is a quite mature field . However, complications arising from the anisotropy of the electron effective mass and from interference among the Bloch wavefunctions at the degenerate conduction-band edges were never fully explored in detail. Both effects were discussed by Andres et al  statistically, in a study of magnetic susceptibility of Si:P system performed in the context of localization and magnetic phase transitions. Silicon conduction band has six minima located along the $`<`$100$`>`$ directions, at about 85% between the center ($`\mathrm{\Gamma }`$ point) and the boundary (X points) of the Brillouin zone (BZ): $`𝐤_{\pm z}=0.85(0,0,\pm 2\pi /\mathrm{a})`$, and the equivalent $`x,y`$ directions. The interference between these valleys causes fast oscillations (on the scale of the inter-atomic spacing, quite different from the slowly varying exchange interaction in the spin-based quantum dot QC model ) in the exchange interaction. Such oscillations cannot be accounted for by a simple calculation using hydrogenic centers. Here we perform a Heitler-London calculation of the exchange coupling  between two substitutional donors in bulk Si, incorporating the effects of both the anisotropic effective mass and the valley interference. Within the envelope function approach, the ground single valence donor electron state (with $`A_1`$ symmetry) is written as a symmetric superposition of the six conduction band minima :

$$\psi (𝐫)=\frac{1}{\sqrt{6}}\underset{\mu }{\overset{6}{}}F_\mu (𝐫)\varphi _\mu (𝐫),$$
(1)
where $`\varphi _\mu (𝐫)=u_\mu (𝐫)e^{i𝐤_\mu 𝐫}`$ are Bloch wavefunctions. Here $`𝐤_\mu `$ refer to the six conduction band minima, and $`F_\mu (𝐫)`$ are the corresponding envelope functions. We use the Kohn-Luttinger variational form for these anisotropic envelopes , e.g.,

$$F_{\pm z}(𝐫)=\frac{1}{\sqrt{\pi a^2b}}e^{\sqrt{(x^2+y^2)/a^2+z^2/b^2}}.$$
(2)
To determine the effective Bohr radii $`a`$ and $`b`$, we have used the most recently measured Si and Ge parameters  (especially longitudinal and transverse effective masses) to perform a variational calculation . The calculated wavefunction widths, as shown in Table I, are quite close to those obtained over 40 years ago . The large difference in transverse mass has only minimal effect on our results.

To calculate the exchange splitting between the ground singlet and triplet states for an impurity pair in Si, we use an approximate form of the Heitler-London approach , which leads to

$`J(𝐑)`$ $`=`$ $`{\displaystyle d^3𝐫_1d^3𝐫_2\psi ^{}(𝐫_1)\psi ^{}(𝐫_2𝐑)\frac{e^2}{ϵ|𝐫_1𝐫_2|}}`$ (4)
$`\times \psi (𝐫_1𝐑)\psi (𝐫_2)`$
$`=`$ $`{\displaystyle \frac{1}{36}}{\displaystyle \underset{\mu ,\nu }{}}\left[{\displaystyle \underset{𝐊,𝐊^{}}{}}|c_𝐊^\mu |^2|c_𝐊^{}^\nu |^2e^{i(𝐊𝐊^{})𝐑}\right]`$ (6)
$`\times j_{\mu \nu }(𝐑)\mathrm{cos}(𝐤_\mu 𝐑)\mathrm{cos}(𝐤_\nu 𝐑),`$
where $`𝐑`$ is the relative position of the impurity nuclei pair, and

$`j_{\mu \nu }(𝐑)`$ $`=`$ $`{\displaystyle d^3𝐫_1d^3𝐫_2F_\mu ^{}(𝐫_1)F_\nu ^{}(𝐫_2𝐑)\frac{e^2}{ϵ|𝐫_1𝐫_2|}}`$ (8)
$`\times F_\mu (𝐫_1𝐑)F_\nu (𝐫_2).`$
The second summation in (6) comes from the reciprocal lattice expansion of the periodic part of the Bloch function, $`u_\mu (𝐫)=_𝐊c_𝐊^\mu e^{i𝐊𝐫}`$, and is identically 1 when $`𝐑`$ is an fcc crystal lattice vector. Integrals in Eq. (4) with rapidly oscillating integrands are neglected in Eq. (6). For $`\mu =\nu `$, we calculate the integrals in (8) by replacing the denominator $`|𝐫_1𝐫_2|`$ by its value for the line along

the impurities , and assume $`j_{\mu \nu }=\sqrt{j_{\mu \mu }j_{\nu \nu }}`$ when $`\mu \nu `$. All are excellent approximations for large separations ($`|𝐑|a,b\mathrm{a}`$).

Germanium conduction band has four minima along the $`<`$111$`>`$ directions, at the zone boundary L points. Straightforward changes in (1) lead to expressions for the exchange coupling in Ge similar to (6), with $`\mu `$ and $`\nu `$ along the $`<`$111$`>`$ directions. Because the conduction band minima are located at the zone boundary, the oscillation in the exchange coupling should display a simpler behavior than in Si when both donors are on the same fcc sublattice, thus making Ge an intriguing candidate for the purpose of quantum computing.

We calculate the exchange energy for a pair of donors at $`𝐑_A`$ and $`𝐑_B`$ as a function of their relative position $`𝐑=𝐑_A𝐑_B`$ from Eq.(6), taking the expression in the square bracket equal to unity. This means that the results are appropriate for any lattice vector $`𝐑`$, while for general inter-donor distances our values should be taken as an estimate. For definiteness, we consider $`𝐑`$ along high-symmetry directions in the crystal lattice. Fig. 1 shows the calculated values of $`J(𝐑)`$ for donors in Si and Ge with $`𝐑`$ along the ,  and  directions (frames (a), (b) and (c) respectively). The solid lines give the results for Si, and exhibit the expected decay of $`J`$ with increasing $`|𝐑|`$ due to the decrease of the donor wavefunction overlap. Other general features of Si and Ge exchange are clearly illustrated in this figure, namely the oscillatory behavior of $`J`$ superimposed on its overall decay with distance, and the strong anisotropy of $`J(𝐑)`$ which is apparent by comparing different frames. Both features are consequences of the host material band structure and have not been considered in detail in previous studies , either for simplicity , or because such effects are averaged out for a random donor distribution . The filled circles in Fig. 1 indicate the accessible values of $`J(𝐑)`$ when the impurity pair in Si is located at lattice sites along the considered directions. Since for Si the conduction band minima correspond to points inside the BZ, the period of oscillation in $`J`$ and the lattice periodicity are not commensurate. We have also investigated the effect of small perturbations in the atomic positions: The open squares give the resulting exchange values when one of the impurities is slightly displaced to off-lattice positions. The set of squares around each circle in Fig. 1 corresponds to displacements along different directions, with the distance from the original lattice position arbitrarily taken as $`\delta =0.235`$ Å (i.e., 10% of the nearest-neighbor distance in Si). The squares follow to a very good approximation the behavior of the calculated $`J(𝐑)`$, with $`𝐑`$ along the unperturbed crystal direction, regardless of the direction of $`\delta `$. The upper bound of $`\mathrm{\Delta }J/J`$ calculated from the small-displacement data is about 2.5 to 5% along the  direction, which means that the requisite control over donor positioning should be much better than 10% of the nearest neighbor distance, i.e. better than 0.235 $`\AA `$, a rather difficult task.

The corresponding results for donor pairs in Ge are given by the dashed lines, diamonds and crosses in Fig. 1. In these calculations we have also assumed nearly free electrons, which means $`c_0^\mu c_{𝐊0}^\mu `$. The main qualitative difference between the calculated exchange coupling for donors in Ge and Si arises from the different locations of the conduction band minima in the BZ. Since for Ge the minima occur at the zone boundary L points, the oscillations in $`J(𝐑)`$ are commensurate with the lattice periodicity. For donors at lattice sites along the  direction, the accessible values of $`J(𝐑)`$ correspond exactly to successive local maxima of Eq.(4). As a consequence, off-lattice displacements from the original lattice sites by 10% of the nearest-neighbors distance in Ge ($`\delta =0.245`$ Å) have negligible effect on the exchange coupling for donors along this direction. Quantitatively, the longer range for the interactions in Ge as compared to Si is due to the larger values of the effective Bohr radii for Ge (see Table I).

Our results indicate that a 31P donor array along the  direction in either of the host materials meets the requirements for quantum computer implementation. The exchange energy for two 31P substitutional impurities 100 to 200 Åapart along this direction ranges between a few tenth of one meV to one $`\mu `$eV, corresponding to less than 100 picosecond to 100 nanosecond gate operation time. These values are robust with respect to small off-lattice displacements.

Displacements of the 31P atoms into neighboring lattice sites, which are bound to occur in the fabrication process due to either uncontrolled implantation or surface diffusion during MBE growth, also deserve careful theoretical investigation. In Fig. 2 we present the calculated $`J(𝐑)`$ for a specific relative position $`𝐑_0`$ along the  direction (black symbols), as well as for values of $`𝐑=𝐑_0+𝚫`$, with $`𝚫`$ ranging over the 4 nearest neighbors, 12 second nearest neighbors and 6 third nearest neighbors in the diamond structure (see inset of Fig. 2). These displacements cause a relatively small change in the inter-donor distance (see horizontal scale in the figure). For third neighbors (blue symbols) the trends in Fig. 1 are reasonably preserved. Surprisingly, most first and second neighbors displacements essentially destroy the exchange coupling. This effect is entirely due to the $`\mathrm{cos}(𝐤_\mu 𝐑)`$ factors in Eq. (6), coming from the plane-wave part of the Bloch functions, and is not directly related to the numerical values or approximations involved in obtaining the $`j_{\mu \nu }`$ coefficients. Therefore we do not attempt to refine the envelope function calculation at this stage, e.g. by including central cell corrections  and inter-valley coupling . The results obtained here are inconsistent with the commonly accepted idea that when the Bohr radii of the envelope functions are much bigger than the lattice parameter, the substitutional positioning of the impurity atoms is not of much importance in determining the exchange coupling.

In a Si QC architecture there are surface gates (so called A and J gates) that are used to control single and two-qubit operations . By adjusting the gate potential(s), one can shift the center of a donor electron wavefunction, thus quickly traversing the fast oscillating terrain in the exchange coupling. The results obtained in this paper also indicate that donor electron exchange depends very sensitively on the applied gate potential(s). Furthermore, whether the electrons are trapped by donor ionic potentials or by the applied gate potentials, the fast oscillation in exchange will persist as long as the electrons are in the bulk of a Si crystal. For a 2-qubit gate of a QC, it is the time integral of the exchange constant $`J`$ that determines the gating time . The implication of the oscillatory exchange is that the A and J gate voltages corresponding to the peak exchange coupling have to be well-controlled, optimally close to a local maximum where the exchange is least sensitive to the gate voltage. Since the oscillatory exchange period is close to lattice spacing, the positioning of the donor electrons by the A and J surface gates must be controlled at least to that precision. Furthermore, if one intends to achieve adiabaticity by increasing the switch-on time of the gate, so that $`J𝑑t=(2n+1)\pi \mathrm{}`$ with $`n1`$ , the error in $`J`$ due to gate inaccuracies would accumulate in the integral. It is thus crucial to precisely control the surface gates in order to minimize the possible errors.

Placing the donors and therefore their electrons in the middle of a symmetric Si quantum well will help reduce the fast phase oscillations because only two of the six bulk Si conduction band valleys would contribute to the donor electron ground state , suppressing the magnitude of the interference effect. In a symmetric quantum well the splitting between the ground state and the next excited state (in analogy to the $`A_1`$ and $`T`$ donor electron states in bulk) is much smaller than in the bulk (crudely estimated to be about 3 meV from symmetry arguments, in contrast to nearly 12 meV in bulk Si; the smaller splitting implies that the adiabatic condition is somewhat harder to satisfy in the quantum well). Strains, electric fields, and/or asymmetry in the quantum well might help eliminate the remaining valley degeneracy  so that all the fast oscillating factors in the exchange constant $`J`$ may be removed, thus leading to a slowly varying $`J`$ and a much easier control of the exchange gate. However, strains themselves are hardly controllable, while asymmetry in a quantum well may introduce additional complications that lead to decoherence. Further investigations of these factors are still ongoing and the results will be reported elsewhere .

As we have demonstrated above, moving one of the donors to its second nearest neighbor position causes strong suppression in exchange coupling between the two donor electrons. This cancelation of exchange due to valley interference might be useful for isolating the neighboring qubits. In a Si QC, for most of the time the qubits should not “talk” to each other. Unwanted interaction would lead to leak of information and decoherence. Therefore, the positions with vanishing exchange interaction might provide a quieter environment compared to an arbitrary pair of positions for qubits.

We thank Ravin Bhatt, Bruce Kane, Ellen Williams, Bob Clark, Rodrigo Capaz, and Luiz Davidovich for helpful discussions, and ARDA and NSA for financial support. BK acknowledges financial support from CNPq (Brazil).