# Fraunhofer response and supercurrent spin switching in Black Phosphorus with
strain and disorder

Mohammad Alidoust  Department of Physics, K.N. Toosi University of Technology,
Tehran 15875-4416, Iran  Morten Willatzen  Beijing Institute of Nanoenergy and
Nanosystems, Chinese Academy of Sciences, No. 30 Xueyuan Road, Haidian
District, Beijing 100083, China  Department of Photonics Engineering,
Technical University of Denmark, DK-2800 Kongens Lyngby, Denmark  Antti-Pekka
Jauho  Center for Nanostructured Graphene (CNG), Department of Micro- and
Nanotechnology, Technical University of Denmark, DK-2800 Kongens Lyngby,
Denmark

###### Abstract

We develop theory models for both ballistic and disordered superconducting
monolayer black phosphorus devices in the presence of magnetic exchange field
and strain. The ballistic case is studied through a microscopic Bogoliubov-de
Gennes formalism while for the disordered case we formulate a quasiclassical
model. Utilizing the two models, we theoretically study the response of
supercurrent to an externally applied magnetic field in two-dimensional black
phosphorus Josephson junctions. Our results demonstrate that the response of
the supercurrent to a perpendicular magnetic field in ballistic samples can
deviate from the standard Fraunhofer interference pattern when the Fermi level
and mechanical strain are varied. This finding suggests the combination of
chemical potential and strain is an efficient external knob to control the
current response in highly sensitive strain-effect transistors and
superconducting quantum interference devices. We also study the supercurrent
in a superconductor-ferromagnet-ferromagnet-superconductor junction where the
magnetizations of the two adjacent magnetized regions are uniform with
misaligned orientations. We show that the magnetization misalignment can
control the excitation of harmonics higher than the first harmonic
$\sin\varphi$ (in which $\varphi$ is the phase difference between the
superconductors) in supercurrent and constitutes a full spin switching current
element. Finally, we discuss possible experimental implementations of our
findings. We foresee our models and discussions could provide guidelines to
experimentalists in designing devices and future investigations.

###### pacs:

74.78.Na, 74.20.-z, 74.25.Ha

##  I introduction

Orthorhombic bulk black phosphorus (BP) under pressure experiences topological
phase transitions, making phosphorus allotropes an attractive new research
areaXiang ; Guo ; Zhao ; Gong . Among phosphorus allotropes, black phosphorus
is the most stable crystal structure Bridgman ; Keyes . Possessing a direct
bandgap that is tunable from $0.3$ eV to $2.0$ eV by applying strain; and
electric field, manipulating the number of layers involved, and in situ doping
with different atoms such as $\rm K$ and $\rm Rb$kim1 ; kim2 , BP offers an
excellent platform with great control over the density of charge carriers by
different means Carvalho ; wylsuperc_exp1 ; Rudenko ; Pereira ; wei ; Y.Ren1 ;
Hedayati . These features of BP are applicable to atomically thin electronics,
modern two-dimensional mechanical transistors, ultrasensitive sensors with
high on-off ratios, and optomechanic devices working under blue or ultraviolet
light Li1 ; Wang2 ; Wu1 ; Sangwan ; liu1 .

Recent research on phosphorus allotropes led to the experimental realization
of phosphoreneliu1 . Also, it has been found that black phosphorus in both
orthorhombic and rhombohedral structures can intrinsically develop
superconductivity Kawamura1 ; Karuzawa2 ; Flores ; Wang1 ; Shirotani ; Livas ;
Q.Huang ; Jun-JieZhang ; YanqingFeng ; R.Zhang ; Yuan . However, the dominant
pairing type in BP has brought up discussions and is still unclear. The same
issue arises in even monolayer graphene which is a much simpler two-
dimensional crystal. Particularly, several different and exotic pairing types
have been proposed, ranging from the usual $s$-wave to $p+ip$, $d+id$, and
$f$-wave pairings Ma ; Faye ; Nandkishore ; Kiesel ; Uchoa  and there is still
no agreement about the dominant pairing symmetry. In effect, theoretical
predictions of pairing symmetries in these systems strongly depend on the
strength of interactions in the model considered, pairing mechanisms, and
approximations made. Therefore, detailed experimental information on the
predominant interactions and pairing mechanism is critically important to
identify the pairing(s) that is (are) the most energetically favorable.

The presence of defects and impurities when synthesizing BP sheets is
inevitable and can influence the physical properties and characteristics of
devices made from them. Therefore, the study of the influence of disorder at a
microscopic level sheds light on the physical origins of experimental
observations and facilitates the optimization of devices made of nonideal BP
sheetsR.Zhang ; Koenraad ; Doganov ; Koenig ; Han ; Zou1 ; Liu2 ; qiu1 . In
this paper, we establish physical models for hybrid structures made of
superconducting and magnetized monolayer BP (from now on, by ‘black
phosphorus’ our mean is ‘monolayer black phosphorus’). We consider two-
dimensional BP Josephson junctions with the possibility of inclusion of
magnetism with arbitrary magnetization configurations, strain, and external
magnetic field as depicted in Fig. 1. We study both ballistic and disordered
devices containing nonmagnetic impurities. In ballistic systems, we develop a
Bogoliubov-de Gennes microscopic theory and employ realistic band parameters
for a BP sheet by means of density functional theory and symmetry
calculations. We study the response of supercurrent to an external magnetic
field and uniaxial/biaxial strain in a superconductor (S)-normal
(N)-superconductor (S) BP junction where the superconductors are identical and
possess different superconducting phases, creating a phase gradient $\varphi$
across the junction. Our results demonstrate that by manipulating the biaxial
strain and Fermi level, the response of critical supercurrent to an external
magnetic field deviates from the standard Fraunhofer interference pattern.
Also, the critical current at current reversal points displays a finite
nonvanishing supercurrent. In general, the current phase relation consists of
multiple harmonics $\sin n\varphi$, $n=\pm 1,\pm 2,...$. We plot the current
phase relation around the current reversal points and show that the appearance
of higher harmonics is responsible for the nonvanishing critical current.
Additionally, we examine the influence of the magnetization misalignment angle
on supercurrent flow in a superconductor (S)-ferromagnet ($\rm
F_{1}$)-ferromagnet ($\rm F_{2}$)-superconductor (S) $\rm SF_{1}F_{2}S$ BP
junction where F regions are of unequal thickness and have uniform
magnetizations with misaligned orientations. Our investigations reveal a full
spin switching of supercurrent upon increasing the magnetizations’ relative
misalignment angle $\theta$. For magnetizations smaller than $0.1$ eV,
increasing the relative misalignment angle removes higher harmonics, while for
stronger magnetizations ($\sim 0.2$ eV) it first induces higher harmonics
close to $0$-$\pi$ crossovers at $\theta\sim\pi/2$ and then fully switches the
supercurrent direction at antiparallel configuration, i.e., $\theta=\pi$. For
disordered systems in the presence of magnetism and superconductivity, we
formulate a quasiclassical Keldysh Green’s function model. To this end, we
utilize a low-energy model Hamiltonian and construct retarded, advanced, and
Keldysh propagators. Applying the quasiclassical approximation, we derive the
Eilenberger equation that is valid for ballistic and moderately disordered
systems where the mean free time of particles $\tau$ is either infinite or
sufficiently large. To further expand our theory, we consider the so-called
diffusive regime, where the impurities scatter quasiparticles in all
directions and the mean free time of particles goes to small values
($\tau\rightarrow 0$); expand the Green’s function in terms of zeroth and
first harmonics, and finally derive the Usadel equation. This approach was
recently generalized for the surface channels of topological insulators
(containing impurities and disorders) as well zyuzin ; bobkova ; alidoustphi0
. We derive charge current density in the diffusive regime and apply our
theory to a two-dimensional SNS BP Josephson junction subject to an external
magnetic field in the weak proximity limit, where the inducted superconducting
gap into the normal region is less than $10\%$ of the gap deep inside the
superconducting regions. Our results demonstrate that the delicate features
explored in the ballistic regime, using the microscopic Bogoliubov-de Gennes
theory, are washed out in the weak proximity limit of the diffusive regime,
leading to the standard Fraunhofer response.

This paper is organized as follows. In Sec. II, we study ballistic systems
using a microscopic Bogoliubov-de Gennes formalism. Specifically, In Sec.
II.1, we study the response of supercurrent to an external magnetic field,
perpendicular to the junction plane, with the incorporation of biaxial and
uniaxial strain and the variation of Fermi level. In Sec. II.2, we study the
supercurrent in a $\rm SF_{1}F_{2}S$ BP Josephson junction and how
magnetization rotation can alter the supercurrent. In Sec. III, we present the
generalized Eilenberger and Usadel theories for BP devices in the presence of
superconductivity and magnetism. We apply this model to a SNS BP Josephson
junction subject to an external magnetic field. Finally, we give concluding
remarks in Sec. IV.

figurec

Figure 1:  (Color online). Schematic of the monolayer black phosphorus
Josephson junctions. The two-dimensional system resides in the $xy$ plane, so
that the interfaces are along the $x$ direction, and the $y$ axis is normal to
the junctions. The uniaxial/biaxial strain is exerted into the plane of the
black phosphorus sheet, and an externally applied magnetic field ${\bm{H}}$ is
directed along $z$, perpendicular to the junction plane. The two $s$-wave
superconducting electrodes (S), with different macroscopic phases
$\varphi_{l,r}$, are proximity coupled to black phosphorus. The chemical
potentials of superconducting parts and nonmagnetized, and magnetized regions
are marked by $\mu_{S}$, $\mu_{N}$, and $\mu_{F1},\mu_{F2}$. (a) The junction
area is nonmagnetized BP of length $d$ and width $W$. (b) The junction area is
magnetized through proximity coupling to ferromagnets ($\rm F_{1}$ and $\rm
F_{2}$). The ferromagnets possess different lengths $d_{F1}\neq d_{F2}$ with
identical widths $W$. The magnetization orientation in each F region is
determined through its angle with respect to the $z$ axis, i.e.,
$\theta_{1,2}$.

##  II ballistic systems

The low-energy effective Hamiltonian of black phosphorus under strain
$\varepsilon_{ii}$ and additionally subject to an external magnetic field with
an associated vector potential $\text{A}_{j}$ can be expressed byVoon1 ; Voon2

$\begin{split}&H=\int\frac{d\textbf{k}}{(2\pi)^{2}}\hat{\psi}^{\dagger}({\textbf{k}})\Big{\\{}\big{[}u_{0}+\alpha_{i}\varepsilon_{ii}+(\eta_{j}+\beta_{ij}\varepsilon_{ii})(\text{k}_{j}-e\text{A}_{j})^{2}\big{]}\tau_{0}+\\\
&+\big{[}\delta_{0}+\mu_{i}\varepsilon_{ii}+(\gamma_{j}+\nu_{ij}\varepsilon_{ii})(\text{k}_{j}-e\text{A}_{j})^{2}\big{]}\tau_{x}-\chi_{y}(\text{k}_{y}-e\text{A}_{y})\tau_{y}\Big{\\}}\hat{\psi}({\textbf{k}}),\end{split}$
(1)

where the indices stand for coordinates, i.e., $i,j\equiv x,y$ , and summation
over repeated indices is assumed. The band parameters calculated through
density functional theory and symmetry computations are given in Table 1. When
$\varepsilon_{ii}$ is negative (positive), the strain in that direction (the
$x$ direction for $\varepsilon_{xx}$ and the $y$ direction for
$\varepsilon_{yy}$) is of the compression (stretch) type. The matrices
$\tau_{i}$ are the Pauli matrices in pseudo spin space, and an exchange field
${\bf h}$ invokes real-spin space indicated by
$\text{h}_{x}\sigma_{x}+\text{h}_{y}\sigma_{y}+\text{h}_{z}\sigma_{z}$ in
which $\sigma_{i}$ are the Pauli matrices in real-spin space. Therefore, in
the presence of magnetism, the associated field operator is given by
$\hat{\psi}^{\dagger}(\textbf{k})=(\psi_{A\uparrow}^{\dagger},\psi_{A\downarrow}^{\dagger},\psi_{B\uparrow}^{\dagger},\psi_{B\downarrow}^{\dagger})$.
Here the sublattices and real spins are labeled by $AB$ and
$\uparrow\downarrow$, respectively.

We describe the superconductivity of a BP sheet through the BCS picture where
particles with opposite spins are coupled by intervalley interactions. This
type of pairing can be expressed by

$\Delta^{AB}_{\uparrow\downarrow}\Big{\langle}\psi^{\dagger}_{A\uparrow}\psi^{\dagger}_{B\downarrow}\Big{\rangle}+\text{H.c.},$
(2)

in which $\Delta^{AB}_{\uparrow\downarrow}$ is the associated superconducting
gap. Note that other pairing types can also be considered, as described in
Ref. BP_Majorana, , and these will influence the final results. Nonetheless,
our investigations in Ref. BP_Majorana,  demonstrated that intravalley spin-
singlet $s$-wave pairing has insubstantial influence. Other symmetry
combinations such as $p$-wave, $d$-wave, $f$-wave, etc., may arise in BP under
strain. The consequences of these combinations will be considered in a future
work. There are some indications that a BP sheet under pressure and/or with
electron doping can be driven into the superconducting phase Guo ; Livas ;
Shirotani ; Jun-JieZhang ; YanqingFeng ; Q.Huang ; Xia ; Wang1 , although
still not conclusively R.Zhang ; Yuan . Nonetheless, we assume that singlet
superconductivity can be extrinsically induced into a BP sheet by means of the
proximity effect when BP is proximity coupled to an $s$-wave superconducting
electrode as shown in Fig. 1. In this case, the superconducting BP sheet can
be described by the following Bogoliubov-de Gennes Hamiltonian in the Nambu
space:

${\cal
H}(\textbf{k})=\left(\begin{array}[]{cc}H(\textbf{k})-\mu&\hat{\Delta}\\\
\hat{\Delta}^{\dagger}&-H^{\text{T}}(-\textbf{k})+\mu\end{array}\right),$ (3)

where $\Delta$ is the proximity-induced superconducting gap and $\mu$ is the
chemical potential. The associated $1\times 8$ vector field operator in k
space can now be expressed as
$\check{\psi}_{\text{BCS}}^{\dagger}(\textbf{k})=[\hat{\psi}^{\dagger}(\textbf{k}),\hat{\psi}(\textbf{-k})]$.
In the calculation of supercurrent below, the energies are given in units of
$|\Delta|$.

Table 1: Band parameters of a monolayer black phosphorus subject to externally applied strain Voon1 ; Voon2 .  $u_{0}$(eV) |  $\delta_{0}$(eV) |  $\alpha_{x}$(eV) |  $\alpha_{y}$(eV) |  $\mu_{x}$(eV)
---|---|---|---|---
-0.42 | +0.76 | +3.15 | -0.58 | +2.65
$\mu_{y}$(eV) |  $\eta_{x}$(eV$\cdot\textup{\AA}^{2}$) |  $\eta_{y}$(eV$\cdot\textup{\AA}^{2}$) |  $\gamma_{x}$(eV$\cdot\textup{\AA}^{2}$) |  $\gamma_{y}$(eV$\cdot\textup{\AA}^{2}$)
+2.16 | +0.58 | +1.01 | +3.93 | \+ 3.83
$\beta_{xx}$(eV$\cdot\textup{\AA}^{2}$) |  $\beta_{yx}$(eV$\cdot\textup{\AA}^{2}$) |  $\beta_{xy}$(eV$\cdot\textup{\AA}^{2}$) |  $\beta_{yy}$(eV$\cdot\textup{\AA}^{2}$) |
-3.48 | -0.57 | +0.80 | +2.39 |
$\nu_{xx}$(eV$\cdot\textup{\AA}^{2}$) |  $\nu_{yx}$(eV$\cdot\textup{\AA}^{2}$) |  $\nu_{xy}$(eV$\cdot\textup{\AA}^{2}$) |  $\nu_{yy}$(eV$\cdot\textup{\AA}^{2}$) |  $\chi_{y}$(eV$\cdot\textup{\AA}$)
-10.90 | -11.33 | -41.40 | -14.80 | +5.25

###  II.1 SNS black phosphorus Josephson junction subject to external magnetic
field

Figure 2:  (Color online). The response of critical current to an external
magnetic field in a ballistic SNS BP Josephson junction. (a) The chemical
potentials of both superconducting sides are varied, $\mu_{S}=\mu$, while the
middle normal BP is undoped, $\mu_{N}=0$. The critical supercurrent is plotted
as a function of magnetic flux passing through the junction area $\Phi$
(normalized by the magnetic flux quantum $\Phi_{0}=h/2e$) for different values
of $\mu_{S}=\mu$ inside the S regions. The circles surround kinks where
supercurrent reversals occur in cases with $\mu=0.4$ eV and $0.8$ eV. (b) The
chemical potential throughout the junction is considered to be the same
$\mu_{S}=\mu_{N}=\mu$ and the critical current is plotted for various values
of $\mu$, $\mu=0.4,0.8,1.2,1.6,2.0$ eV. In (c), we show the current phase
relation at a current reversal point that the external magnetic field induces.
We consider a representative case, which is $\mu=0.4$ eV in panel (a). The
arrows indicate the change in maximum supercurrent when varying the external
magnetic field from $0$ to $\sim 3.0\Phi_{0}$. The compression in (a)-(c) is
kept fixed at $\varepsilon_{xx}=0.0,\varepsilon_{yy}=-0.2$. In panels
(d1)-(e2) we consider an undoped middle normal region $\mu_{N}=0$ and set
differing scenarios for an applied strain. In (d1) and (d2) we set
$\mu_{S}=\mu=1.5$ and $3.0$ eV in the superconducting regions and uniaxial
stress $\varepsilon_{xx}=0.0,\varepsilon_{yy}=-0.12,-0.16,-2.0,-2.4,-2.8$. In
(e1) and (e2), we have increased $\mu$ to $2.0,4.0$ eV and set equal biaxial
compression components in both the $x$ and $y$ directions:
$\varepsilon_{xx}=\varepsilon_{yy}=-0.12,-0.16,-2.0,-2.4,-2.8$.

To begin, we consider a two-dimensional SNS BP Josephson junction as depicted
in Fig. 1(a). The plane of BP is in the $xy$ plane, and the interfaces reside
in the $x$ direction, so that the $y$ axis is perpendicular to the interfaces.
The junction has a finite size in the $xy$ plane and has a length and width of
$d$ and $W$, respectively. An external magnetic field $\bm{H}$ is exerted
perpendicular to the junction plane and directed along the $z$ axis. To
account for the external magnetic field, we consider a vector potential
$\textbf{A}=(\text{A}_{x},\text{A}_{y},\text{A}_{z})=(0,xH_{z},0)$ so that A
satisfies the Lorentz condition $\bm{\nabla}\cdot\textbf{A}=0$ and produces
the external field $H_{z}=\bm{\nabla}\times\textbf{A}$ jap_ma . We also
consider the possibility of the inclusion of stress applied in the $x$ and $y$
directions. To evaluate the response of supercurrent to an external magnetic
field, we calculate charge current density using the quantum definition,
involving Hamiltonian (3), i.e.,

$\begin{split}\frac{\partial\rho}{\partial
t}=\lim\limits_{{\bm{r}}\rightarrow{\bm{r}}^{\prime}}\sum\limits_{\rho\tau\sigma\rho^{\prime}\tau^{\prime}\sigma^{\prime}}\frac{1}{i\hbar}\Big{[}\psi^{\dagger}_{\rho\tau\sigma}({\bm{r}}^{\prime}){\cal
H}_{\rho\tau\sigma\rho^{\prime}\tau^{\prime}\sigma^{\prime}}({\bm{r}})\psi_{\rho^{\prime}\tau^{\prime}\sigma^{\prime}}({\bm{r}})\\\
-\psi^{\dagger}_{\rho\tau\sigma}({\bm{r}}^{\prime}){\cal
H}_{\rho\tau\sigma\rho^{\prime}\tau^{\prime}\sigma^{\prime}}^{\dagger}({\bm{r}}^{\prime})\psi_{\rho^{\prime}\tau^{\prime}\sigma^{\prime}}({\bm{r}})\Big{]},\end{split}$
(4)

where the time variation of charge density $\rho$ can be attributed to charge
sources and sinks. Here ${\cal
H}_{\rho\tau\sigma\rho^{\prime}\tau^{\prime}\sigma^{\prime}}$ is the component
form of Eq. (3) with spin, valley, and particle-hole indices. Throughout the
paper, we consider a steady-state regime where $\partial\rho/\partial t=0$. To
compute the total current, we integrate the charge current density
perpendicular to the interface $J_{y}$ over the $x$ direction, i.e.,
$I(\varphi)=\int_{0}^{W}dxJ_{y}(x,y,\varphi)$. Note that the current is
independent of the $y$-coordinate as required by charge conservation. We
obtain appropriate spinors $\psi_{\rho\tau\sigma}$ within the normal BP region
by matching wavefunctions $\hat{\psi}_{l}=\hat{\psi}_{r}$ and applying the
continuity condition $\partial_{\textbf{k}}{\cal
H}_{l}(\textbf{k})\hat{\psi}_{l}=\partial_{\textbf{k}}{\cal
H}_{r}(\textbf{k})\hat{\psi}_{r}$, at the left and right superconductor
interfaces. At the left boundary, ${\cal H}_{l}$ and ${\cal H}_{r}$ stand for
the Hamiltonians of the superconducting and normal regions ($\hat{\psi}_{l}$
and $\hat{\psi}_{r}$ are their associated wavefunctions), respectively,
whereas at the right boundary they stand for the Hamiltonians of the normal
and superconducting regions, respectively. The resultant analytic expressions
are very cumbersome, constituting $1\times 8$ spinors. We do not give the
explicit expressions but employ them directly in evaluating observable
quantities numerically. Supercurrent and its response to externally
controllable agents such as magnetic field and stress are often measured in
experiments. Figure 2 illustrates the behavior of supercurrent passing through
the device shown in Fig. 1(a). The junction length and width are fixed at $7$
nm and $20$ nm, respectively. On the one hand, the device experiences
compression in the $y$ direction, i.e.,
$\varepsilon_{xx}=0,\varepsilon_{yy}=-0.2$. On the other hand, the junction is
exposed to an external magnetic field where the magnetic flux penetrating the
junction area is given by $\Phi=\pi dW{H}_{z}$. Theoretically, it has been
found that BP can support high strains as large as $40\%$ without any rupture
or dislocations Carvalho ; wei ; Li1 . In Fig. 2(a), the middle segment of the
SNS junction [Fig. 1(a)] is undoped, $\mu_{N}=0$, while the chemical potential
in the superconducting segments is varied: $\mu_{S}=\mu=0.4,0.8,1.2,1.6,2.0$
eV. In Fig. 2(b) we set the chemical potential of the middle normal BP area
equal to the superconducting parts, $\mu_{N}=\mu_{S}=\mu$. Increasing the
chemical potential $\mu$ in the first case, the supercurrent is suppressed,
while in the second case, the supercurrent first increases and then decreases.
Our investigations of the density of states through the retarded Green’s
function, $N(\omega_{n};\textbf{r})=-\pi^{-1}{\rm Im}\\{{\rm
Tr}[G^{R}(\omega_{n};\textbf{r})]\\}$, show that this behavior has a direct
link to the nucleation of subgap bound states that the supercurrent flows
through them (the so-called Andreev bound states). Increasing the number of
Andreev states, and placing them closer to the superconducting gap edge
increase the number of Cooper pairs that can pass across the junction from one
superconductor to another under the applied superconducting phase gradient
($\varphi=\varphi_{l}-\varphi_{r}$). In both Figs. 2(a) and 2(b) we see that
the critical supercurrent is nonvanishing even close to the supercurrent
reversal points surrounded by circles. To illustrate the origin of
nonvanishing critical current, we plot the charge current as a function of the
superconducting phase difference between the two superconductors $\varphi$ for
a gradually increasing $\Phi$ from $0$ to $\approx 3\Phi_{0}$ in Fig. 2(c).
The supercurrent at zero flux $\Phi$ is proportional to $\sin\varphi$. By
increasing $\Phi$, the supercurrent starts to deviate from the sinusoidal
form. figurec

t!

Figure 3:  (Color online). Current phase relation $I(\varphi)$ in a $\rm
SF_{1}F_{2}S$ BP Josephson junction. The current is plotted with increasing
magnetization misalignment angle $\theta_{1}=0.1\pi-1.0\pi$ with a step of
$0.1\pi$. The magnetization in $\rm F_{2}$ is kept fixed along the $z$
direction, $\theta_{2}=0$. Hence, $\theta_{1}=\pi/2$ is equivalent to
perpendicular magnetizations, while $\theta_{1}=\pi$ means antiparallel
alignment of magnetizations. In panels (a) and (b) we set equal magnetization
strengths in both $\rm F_{1}$ and $\rm F_{2}$, i.e.,
$|\textbf{h}_{1}|=|\textbf{h}_{2}|=\text{h}_{0}$ and change it as (a):
$\text{h}_{0}=0.1$ eV, (b): $\text{h}_{0}=0.2$ eV.

The long downward arrow indicates the path that the maximum supercurrent
traverses by increasing $\Phi$. Close to the first supercurrent reversal point
in Fig. 2(a), the first harmonic $\sin\varphi$ is highly suppressed and higher
harmonics $\sin 2\varphi,\sin 3\varphi,...$ dominate. Note that because of the
presence of these higher harmonics the supercurrent is always nonzero for all
values of $\Phi$. A further increase in $\Phi$ induces a supercurrent
reversal, as the small downward arrow indicates, and reverts the current phase
relation to the one proportional to $\sin\varphi$, which is similar to the
$\Phi=0$ case except with a minus sign. With a further increase in $\Phi$ the
supercurrent mimics and repeats the behavior described above. Next, in Figs.
2(d1)-(e2) we examine the influence of uniaxial stress. In these panels we
keep the middle BP undoped, $\mu_{F1}=\mu_{F2}=0$, and change $\mu_{S}=\mu$ in
the superconducting regions to $1.5,3.0,2.0,4.0$ eV from left to right.
Figures 2(d1) and 2(d2) exhibit the effect of varying stress,
$\varepsilon_{yy}=-0.12,-0.16,-0.20,-0.24,-0.28$, in the absence of any
component in the $x$ direction; $\varepsilon_{xx}=0$. Increasing the
compressive strain leads to a higher critical supercurrent. Interestingly, the
contribution of higher harmonics to the supercurrent is now small and the
critical supercurrent is vanishingly small at the supercurrent reversal
points, which is the opposite of the cases shown in Fig. 2(b). In panels 2(e1)
and 2(e2) we set an isotropic biaxial strain
$\varepsilon_{xx}=\varepsilon_{yy}=-0.12,-0.16,-0.20,-0.24,-0.28$ with
strengths identical to those in Figs. 2(d1) and 2(d2). As can be seen, the
compressive strains $\varepsilon_{xx}$ and $\varepsilon_{yy}$ highly enhance
the supercurrent when comparing Fig. 2(d1) to Fig.2(e1). We also see that the
contribution of higher harmonics close to the supercurrent reversal points is
recovered, specifically for high strains
$\varepsilon_{xx}=\varepsilon_{yy}=-0.24,-0.28$. The more prominent feature is
shown in Fig. 2(e2). We see that not only do the higher harmonics prohibit
nonzero supercurrent at the current reversal points, but also high strains
$\varepsilon_{xx}=\varepsilon_{yy}=-0.20,-0.24,-0.28$ now result in non-
Fraunhofer responses of supercurrent to the external magnetic field. At
$\varepsilon_{xx}=\varepsilon_{yy}=-0.20$ the peak of the critical
supercurrent at $\Phi=0$ is suppressed, and its amplitude is smaller than the
second and third peaks. A further increase in compressive strains, i.e.,
$\varepsilon_{xx}=\varepsilon_{yy}=-0.24,-0.28$, causes a dip at zero external
field, and the first peak of supercurrent appears near $\Phi\approx\Phi_{0}$.
Note that the junction has a finite width of $20$ nm, and it results in
nonideal Fraunhofer responses, as exhaustively studied in Ref. jap_ma,  for
standard metallic systems.

###  II.2 ${\rm SF_{1}F_{2}S}$ black phosphorus Josephson junction

In this section we expand our previous studies of ballistic systems to the
$\rm SF_{1}F_{2}S$ devices depicted in Fig. 1(b). The magnetization and
superconductivity can be induced in the BP sheet by means of the proximity
effect. The two superconductors are coupled through a bilayer of uniformly
magnetized ferromagnets, $\rm F_{1}F_{2}$, with different thicknesses
$d_{F1}\neq d_{F2}$ and misaligned magnetization orientations
$\theta_{1}\neq\theta_{2}$. The arrows on top of the F regions in Fig. 1(b)
stand for the magnetization direction in each F layer that make angles
$\theta_{1,2}$ with the $z$ axis perpendicular to the junction plane. Without
loss of generality, we consider in-plane magnetizations, i.e.,
${\mathbf{h}}=(\text{h}_{x},0,\text{h}_{z})$. The relative angle between the
magnetization of the two F segments can be controlled in experiment by
choosing different ferromagnetic materials for each F and applying an in-plane
external magnetic field in the $x$ direction. It is known that various
magnetic materials respond differently to an external magnetic field. The
magnetization of strongly magnetized materials, e.g, $\rm Co$ or $\rm LCMO$
compounds, rotates harder and slower with respect to weakly magnetized
materials, e.g., $\rm Py$ or $\rm NiFe$ in a given external magnetic field. An
external magnetic field perpendicular to the junction plane can induce
superconducting vortices that make analyses and the isolation of the pure
effect of magnetization rotation inconclusive jap_ma ; robin . The in-plane
external field in the $x$ direction, however, ensures that vortices are not
generated, and therefore, if we choose $\rm F_{1}/F_{2}\equiv LCMO/Py$, the
role of the external field is limited to only the induction of misalignment in
the relative magnetization orientation of $\rm F_{1}/F_{2}$ robin . Figure 3
illustrates the behavior of the supercurrent phase relation $I(\varphi)$ when
the magnetization misalignment angle increases. To obtain these results, we
fix $d_{F1}=10$ nm and $d_{F2}=5$ nm and change the magnetization strength
$|{\mathbf{h}}_{1}|=|{\mathbf{h}}_{2}|=0.1,0.2$ eV in Figs. 3(a) and 3(b),
respectively. We chose a representative strain set
$\varepsilon_{xx}=0,\varepsilon_{yy}=-0.2$ and chemical potential
$\mu_{N}=0,\mu_{S}=0.2$ eV. We set $\theta_{2}=0$ and vary $\theta_{1}$ from
zero, equivalent to parallel alignment, to $\pi$, equivalent to antiparallel
magnetization alignment, by a step of $0.1\pi$. In Fig. 3(a), when
$\theta_{1}=0$, the supercurrent exhibits a sign change at a superconducting
phase difference of $\varphi\approx 0.6\pi$. This current reversal occurs due
to strong contributions of higher harmonics, i.e., $\sin 2\varphi,\sin
3\varphi,...$. The magnetization misalignment weakens this contribution and
removes the higher harmonics, so that at $\theta_{1}=0.6\pi$ the current
reversal is fully removed, and when the magnetizations are antiparallel, we
find the usual first harmonic $I(\varphi)\propto\sin\varphi$. The
magnetization rotation from a parallel to antiparallel configuration causes a
full switching of the supercurrent flow. This finding is suggestive of
creating experimentally well controlled spin switching devices using the BP
sheets. In Fig. 3(b) when the magnetization strength is $0.2$ eV, we see that
the weak phase shift at $\varphi=0$, also seen in Fig. 3(a), is now
pronounced. The magnetization misalignment removes it, induces a strong $\sin
2\varphi$ component at $\theta_{1}=0.5\pi$, and then at $\theta_{1}=\pi$
switches the supercurrent flow direction, including the phase shift at zero
phase difference $\varphi=0$, which now reverts to a negative value. Other
parameters such as junction thickness and temperature can also induce current
reversals and higher harmonics in the current phase relationship.

##  III nonideal and disordered systems

To model nonmagnetic impurities, we formulate the quasiclassical model for a
BP with the inclusion of superconductivity and magnetism. To this end, we
express the Hamiltonian Eq. (1) by redefining parameters, which simplifies our
subsequent notation, as follows

$\begin{split}H=\int\frac{d\textbf{k}}{(2\pi)^{2}}\hat{\psi}^{\dagger}_{\textbf{k}}\Big{\\{}&\sum_{k=x,y}\eta_{k}(\text{k}_{k}-e\text{A}_{k})^{2}\tau_{0}+\gamma_{k}(\text{k}_{k}-e\text{A}_{k})^{2}\tau_{x}\\\
&-\chi_{y}(\text{k}_{y}-e\text{A}_{y})\tau_{y}+\mu_{0}\tau_{0}+\mu_{x}\tau_{x}\Big{\\}}\hat{\psi}_{\textbf{k}}.\end{split}$
(5)

In the presence of spin and sublattices $A,B$, we define propagators antti

$\displaystyle
G_{\tau\sigma\tau^{\prime}\sigma^{\prime}}(t,t^{\prime};\mathbf{r},\mathbf{r}^{\prime})$
$\displaystyle=$ $\displaystyle-\Big{\langle}{\cal
T}\Psi_{\tau\sigma}(t,\mathbf{r}^{\prime})\Psi_{\tau^{\prime}\sigma^{\prime}}^{{\dagger}}(t^{\prime},\mathbf{r}^{\prime})\Big{\rangle},~{}~{}~{}~{}$
(6a)
$\displaystyle\bar{G}_{\tau\sigma\tau^{\prime}\sigma^{\prime}}(t,t^{\prime};\mathbf{r},\mathbf{r}^{\prime})$
$\displaystyle=$ $\displaystyle-\Big{\langle}{\cal
T}\Psi^{{\dagger}}_{\tau\sigma}(t,\mathbf{r})\Psi_{\tau^{\prime}\sigma^{\prime}}(t^{\prime},\mathbf{r}^{\prime})\Big{\rangle},~{}~{}~{}~{}$
(6b) $\displaystyle
F_{\tau\sigma\tau^{\prime}\sigma^{\prime}}(t,t^{\prime};\mathbf{r},\mathbf{r}^{\prime})$
$\displaystyle=$ $\displaystyle+\Big{\langle}{\cal
T}\Psi_{\tau\sigma}(t,\mathbf{r})\Psi_{\tau^{\prime}\sigma^{\prime}}(t^{\prime},\mathbf{r}^{\prime})\Big{\rangle},~{}~{}~{}~{}$
(6c) $\displaystyle
F^{{\dagger}}_{\tau\sigma\tau^{\prime}\sigma^{\prime}}(t,t^{\prime};\mathbf{r},\mathbf{r}^{\prime})$
$\displaystyle=$ $\displaystyle+\Big{\langle}{\cal
T}\Psi^{{\dagger}}_{\tau\sigma}(t,\mathbf{r})\Psi^{{\dagger}}_{\tau^{\prime}\sigma^{\prime}}(t^{\prime},\mathbf{r}^{\prime})\Big{\rangle},~{}~{}~{}~{}$
(6d)

where $\Psi_{\tau\sigma}$ are field operators, ${\cal T}$ is the time-ordering
operator, and $t,t^{\prime}$ are the imaginary times at
$\mathbf{r},\mathbf{r}^{\prime}$ locations, respectively. We consider the
elastic scattering potential $V(\mathbf{r})$ in a BP sheet by the self-energy
term

$\Sigma_{\text{imp}}(\mathbf{r}-\mathbf{r}^{\prime})=\Big{\langle}V(\mathbf{r})G(\mathbf{r},\mathbf{r}^{\prime})V(\mathbf{r}^{\prime})\Big{\rangle},$
(7)

where we average over the locations of impurities. To find the mean free time
of particles in the disordered BP, we neglect anisotropic terms. Therefore,
assuming isotropic scattering we find the mean free time as $\tau^{-1}=2\pi
n_{i}N_{0}\int d\Theta_{{\bm{n}}_{\text{F}}}(2\pi)^{-1}|v(\Theta)|^{2}$, in
which $v(\Theta)$ is the Fourier transform of the scattering potential that
depends on the relative angle $\Theta$ between the particles’ incidence
direction and the particles’ scattering direction, $n_{i}$ is the
concentration of nonmagnetic impurities, and $N_{0}$ is the density of states
per spin at the Fermi level of the system. In what follows, we assign the
Pauli matrices $\sigma_{0,x,y,z}$ to real spin, $\tau_{0,x,y,z}$ to pseudo
spin, and $\rho_{0,x,y,z}$ to the particle-hole in the presence of
superconductivity. In the particle-hole space we find the following equation
for the Green’s function:

$\displaystyle\left(\begin{array}[]{cc}-i\omega_{n}+\hat{H}(\mathbf{r})&-\hat{\Delta}(\textbf{r})\\\
\hat{\Delta}^{\dagger}(\textbf{r})&i\omega_{n}+\tau_{y}\sigma_{y}\hat{H}^{*}(\mathbf{r})\tau_{y}\sigma_{y}\end{array}\right)\check{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})$
(10) $\displaystyle=\delta(\mathbf{r}-\mathbf{r}^{\prime})+\frac{1}{2\pi
N_{0}\tau}\check{G}(\omega_{n};\mathbf{r},\mathbf{r})\check{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime}),$
(11)

in which $\omega_{n}=\pi(2n+1)k_{B}T$ is the Matsubara frequency,
$n\in\mathbb{Z}$, $k_{B}$ is the Boltzmann constant, $T$ is temperature, and
$\hat{\Delta}^{\dagger}({\bm{r}})$ is the Hermitian conjugation of
$\hat{\Delta}({\bm{r}})$: the proximity-induced superconducting gap. The
matrix form of the Green’s function can be expressed by

$\check{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})=\left(\begin{array}[]{cc}-\hat{G}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})&-i\hat{F}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})\tau_{y}\sigma_{y}\\\
-i\tau_{y}\sigma_{y}\hat{F}^{\dagger}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})&\tau_{y}\sigma_{y}\hat{\bar{G}}(\omega_{n};\mathbf{r},\mathbf{r}^{\prime})\tau_{y}\sigma_{y}\end{array}\right).$

Here we have denoted $4\times 4$ matrices by using a hat symbol,
$\hat{\square}$, and $8\times 8$ matrices by a check symbol,
$\check{\square}$. We now subtract from Eq. (10) its conjugate and Fourier
transform with respect to the relative coordinates:
$\mathbf{R}=(\mathbf{r}+\mathbf{r}^{\prime})/2$ and
$\delta\mathbf{r}=\mathbf{r}-\mathbf{r}^{\prime}$. In order to simplify our
calculations, we assume that the Green’s function is localized at the Fermi
level and define the quasiclassical Green’s function

$\displaystyle\check{g}(\omega_{n};\mathbf{R},\mathbf{n}_{\text{F}})=\frac{i}{\pi}\int
d\xi_{p}\check{G}(\omega_{n};\mathbf{R},\mathbf{p}),$ (12)

in which $d\xi_{p}=v_{F}dp$ and $v_{F}$ is the Fermi velocity. Incorporating
these assumptions, we finally arrive at the Eilenberger equation eiln :

$\displaystyle p_{\text{F}}^{k}\Big{\\{}{\cal
J}_{k},{\check{\tilde{\nabla}}_{k}\check{g}\Big{\\}}+\Big{[}\omega_{n}\rho_{z}-i\check{\Delta}-i\check{\text{M}}+i{\chi_{y}}p_{\text{F}}^{y}\tau_{y}+\frac{1}{2\tau}\langle\check{g}\rangle,\check{g}}\Big{]}=0,$
$\displaystyle{\cal J}_{k}=\eta_{k}+\gamma_{k}\tau_{z},$ (13)
$\displaystyle\check{\tilde{\nabla}}_{k}\check{\text{X}}\equiv\check{\nabla}_{k}\check{\text{X}}-\Big{[}ie\text{A}_{k}\rho_{z},\check{\text{X}}\Big{]},$

where the average over disorder is indicated by $\langle...\rangle$ and
$\check{\text{M}}$ stands for the magnetization. The momentum at the Fermi
surface in the $k$ direction is shown by $p_{F}^{k}$. Here, we have assumed
that the Fermi energy is large enough so that the Fermi wavelength,
$\lambda_{F}$, is negligible compared to the spatial variation $\xi_{\bm{r}}$
of observable quantities, i.e., $\xi_{\bm{r}}\gg\lambda_{F}$.

The Eilenberger equation can be further simplified in systems with a high
density of impurities so that $\tau^{-1}\gg|\omega_{n}|,|\Delta|$. In this
case, the quasiparticles move diffusively with random directions and
trajectories, which is the so-called diffusive regime usadel . In the
diffusive regime, we integrate the quasiclassical Green’s function, Eq. (12),
over all possible directions of the quasiparticles’ momentum:

$\displaystyle\Big{\langle}\check{g}(\omega_{n};\mathbf{R},\mathbf{n}_{\text{F}})\Big{\rangle}\equiv\int\frac{d\Omega_{\mathbf{n}_{\text{F}}}}{2\pi}\check{g}(\omega_{n};\mathbf{R},\mathbf{n}_{\text{F}}),\;\mathbf{n}_{\text{F}}=\frac{\mathbf{p}_{\text{F}}}{|\mathbf{p}_{\text{F}}|}.~{}~{}~{}~{}~{}\;\;$
(14)

In this regime, the Green’s function can be expanded through the first two
harmonics, $s$-wave and $p$-wave:

$\check{g}(\omega_{n};\mathbf{R},\mathbf{n}_{\text{F}})=\check{g}_{s}(\omega_{n};\mathbf{R})+{n}_{\text{F}}^{k}\check{{g}}_{p}^{k}(\omega_{n};\mathbf{R}),$
(15)

where the $s$-wave harmonic in the expansion is isotropic and much larger than
the $p$-wave harmonic:
$\check{g}_{s}\gg{n}_{\text{F}}^{k}\check{{g}}_{p}^{k}$. By substituting this
expanded Green’s function into Eq. (13) and performing an integration over
momentum directions we obtain

$\displaystyle\check{g}_{p}^{k}=-\tau
p_{\text{F}}^{k}\check{g}_{s}\Big{\\{}{\cal
J}_{k},\check{\tilde{\nabla}}_{k}\check{g}_{s}\Big{\\}}-\tau
p_{\text{F}}^{y}\check{{g}}_{s}\Big{[}i\chi_{y}\tau_{y},\check{g}_{s}\Big{]}.$
(16)

Next, we substitute Eq. (16) into Eq. (13), assume that
$\nabla_{k}\gamma_{k}=\nabla_{k}\eta_{k}=\nabla_{y}\chi_{y}=\nabla_{k}\text{A}_{k}=0$,
and find the generalized Usadel equation usadel :

$\frac{p_{\text{F}}^{k}}{2}\Big{\\{}{\cal
J}_{k},\check{\tilde{\nabla}}_{k}\check{{g}}_{p}^{k}\Big{\\}}+\frac{p_{\text{F}}^{y}}{2}\Big{[}i\chi_{y}\tau_{y},\check{{g}}_{p}^{k}\Big{]}+\Big{[}\omega_{n}\rho_{z}-i\check{\Delta}-i\check{\text{M}},\check{g}_{s}\Big{]}=0.$
(17)

We next apply the quasiclassical model (17) to the SNS hybrid structure
displayed in Fig. 1(a). We consider interfaces with a low transparency (the
so-called tunneling limit) between the superconducting and normal BP regions.
Therefore, the following boundary condition describes the coupling between the
superconductor and the normal parts boundary_c1 ; boundary_c2 :

$\zeta{n}_{k}\check{{g}}_{p}^{k}=\Big{[}\check{g}_{s},\check{g}_{\text{SC}}\Big{]},$
(18)

in which $\zeta$ is the ratio between the resistance of the barrier region and
the resistance in the normal region that controls the proximity effect at the
boundaries, ${n}_{k}$ is a unit vector perpendicular to the boundary, and
$\check{g}_{\text{SC}}$ is the Green’s function of the bulk superconducting
segment. To study charge transport, we derive an expression for the charge
supercurrent flow (due to the superconducting phase gradient across the
device) in the $y$ direction [see Fig. 1(a)]. The quantum definition of
current density, derived from Eq. (4), is expressed through the Hamiltonian,
Eqs. (5). As stated in the previous section, we consider a steady state regime
and therefore set $\partial\rho/\partial t=0$ in Eq. (4). After some
calculations using Eq. (4), we finally arrive at the following expression for
the current density in the diffusive regime:

${J}_{k}({\bm{r}})=\frac{ie\pi}{2}N_{0}\text{p}_{\text{F}}^{k}T\sum_{n}\mathrm{Tr}\Big{[}\rho_{z}\big{(}\eta_{k}+\gamma_{k}\tau_{z}\big{)}\check{{g}}_{p}^{k}\Big{]}.$
(19)

To obtain Eq. (19) we have assumed a sufficiently small $\chi_{y}$ and
neglected terms of the order of $\chi_{y}\text{p}_{\text{F}}^{-1}$ and once
more assumed that $\xi_{\bm{r}}\gg\lambda_{F}$. We note that the formulated
quasiclassical Eilenberger and Usadel formalisms can be extended to other two-
dimensional materials. To the best of our knowledge, our work is the first
that discusses the Eilenberger and Usadel approaches for BP and no counterpart
of graphene type is yet available, which could be a potential subject for our
future research zyuzin ; bobkova ; alidoustphi0 .

Figure 4:  (Color online). The critical current in a diffusive SNS BP
Josephson junction as a function of external magnetic flux $\Phi$. Here we set
$\eta=1$, $\gamma=0.5$ and vary $\chi_{y}=0.1,0.2,0.3$.

In order to find the charge current density flowing across the junction
depicted in Fig. 1(a), we solve Eq. (17) together with proper boundary
conditions, Eq. (18), and substitute the resultant Green’s function into Eq.
(19) for the diffusive regime. For the ballistic and/or moderately disordered
systems, one solves Eq. (13) and uses the ballistic counterpart of the current
density Eq. (19) alidoustphi0 . In the diffusive regime, the Usadel equation
(17) results in nonlinear boundary value differential equations that must be
evaluated numerically vrtx_ali ; jap_ma . To further simplify our
calculations, we have Taylor expanded the Green’s function around its bulk
solution $\check{g}_{0}$, i.e., $\check{g}_{s}\simeq\check{g}_{0}+\check{f}$.
This limit is accessible in systems with weak proximity coupling between the
superconductors and normal BP. Still, after performing these approximations,
the resultant expressions are large coupled partial differential equations.
Rather than presenting them explicitly, we evaluate them numerically jap_ma .
To calculate the supercurrent passing through the set up depicted in Fig.
1(a), we consider $\textbf{A}=(0,xH_{z},0)$. The junction length and width are
set to $d=\xi_{S}$ and $W=10\xi_{S}$, respectively, where $\xi_{S}$ is the
superconducting coherence length. Figure 4 shows the critical current as a
function of applied magnetic flux $\Phi$. Because this simplified model is
parametric and the quasiclassical approximations are applied, we have set
representative values for parameters: $\eta=1,\gamma=0.5$ and
$\chi_{y}=0.1,0.2,0.3$. The applied external magnetic field can induce a
Zeeman field that can be large, depending on the $g$-factor of BP. The
combination of Zeeman field and spin-orbit coupling can result in the
generation of superconducting triplet correlations niu ; Birge ; Konschelle ;
Bergeret ; ali_so1 ; ali_so2 ; H.Chakraborti . Hence, we have considered a
dominant spin-singlet superconductivity and added a small component of triplet
pairings that can occur when making hybrid interfaces of superconducting BP
subject to an external magnetic field. The critical current in the weak
proximity limit of the diffusive regime shows the standard Fraunhofer
patterns, and the variation of the parameters involved can simply change the
overall amplitude of maximum current flow across the junction. As shown, in
this weak proximity limit the delicate features explored in the ballistic
regime (Sec. II), such as nonzero supercurrent at the current reversal points
and suppression of the central peak, vanish. Note that the results for the
full proximity limit, however, can be closer to those explored in the
ballistic limit (Sec. II) and might retrieve the manifestation of the second
harmonic discussed earlier vrtx_ali . The methods developed in this work can
be applied to the other limits (i.e., the full proximity limit, but accounting
for impurities) and different geometries as well, as experimental structures
emerge.

##  IV conclusion

In summary, utilizing the low-energy effective Hamiltonian for black
phosphorus, we formulated the Bogoliubov-de Gennes microscopic theory and
quasiclassical model to describe ballistic and disordered black phosphorus
samples in the presence of superconductivity and magnetism. We also
incorporated an external magnetic field and biaxial/uniaxial stress in the
models. Particularly, we studied the responses of supercurrent in a two-
dimensional superconductor-normal-superconductor black phosphorus Josephson
junction to an external magnetic field perpendicular to the junction plane.
Our results demonstrate that by properly tuning the stress and Fermi energy,
the critical supercurrent deviates significantly from the standard Fraunhofer
interference pattern, so that the central peak is suppressed, while the
amplitudes of the next peaks are larger than the central peak with nonzero
current at supercurrent reversal points. We showed that the nonzero critical
current at current reversal points is a direct consequence of the appearance
of harmonics higher than the first usual harmonic in supercurrent around the
current reversal points that prohibits fully vanishing current. Furthermore,
we studied the influence of magnetization misalignment in a superconductor-
ferromagnet-ferromagnet-superconductor junction where the ferromagnetic
regions possess unequal thicknesses and magnetization orientations. Our
results show a full supercurrent switching effect through increasing the
magnetization misalignment angle. To complement our theory, we considered
nonmagnetic impurities in the system. Employing the quasiclassical
approximations, we derived Eilenberger and Usadel equations. The former is
applicable to systems with a moderate density of impurities, while the latter
describes systems containing a high density of impurities so that the
quasiparticles have diffusive motions. We applied this model to a two-
dimensional superconductor-normal-superconductor black phosphorus junction
subject to an external magnetic field. Our investigation in the weak proximity
limit of the diffusive regime showed the standard Fraunhofer response of
supercurrent to an external magnetic field. Nevertheless, a full proximity
limit can modify the results as this limit is closer to the ballistic regime.
Our models and results provide a detailed explanation of the supercurrent
behavior in BP Josephson junctions in the presence of impurities, external
magnetic field, magnetic exchange field, and strain. Furthermore, our findings
by the application of these models to specific configurations demonstrate that
both strain and the orientation of magnetic exchange fields offer effective
tools to control the behavior of current in strain-effect BP transistors and
highly sensitive BP-based devices such as superconducting quantum interference
devices.

###### Acknowledgements.

M.A. is supported by Iran’s National Elites Foundation (INEF). M.A. would like
to thank A. Zyuzin and G. Sewell for useful discussions. Support from the
Talent 1000 Foreign Expert program of China is acknowledged. Center for
Nanostructured Graphene is supported by the Danish National Research
Foundation (Project No. DNRF103).

## References

  * (1) Z. J. Xiang, G. J. Ye, C. Shang, B. Lei, N. Z. Wang, K. S. Yang, D. Y. Liu, F. B. Meng, X. G. Luo, L. J. Zou, Z. Sun, Y. Zhang, and X. H. Chen, Pressure-induced electronic transition in black phosphorus, Phys. Rev. Lett. 115, 186403 (2015).
  * (2) J. Zhao, R. Yu, H. Weng, and Z. Fang, Topological node-line semimetal in compressed black phosphorus, Phys. Rev. B 94, 195104 (2016).
  * (3) P. L. Gong, D. Y. Liu, K. S. Yang, Z. J. Xiang, X. H. Chen, Z. Zeng, S. Q. Shen, and L. J. Zou, Hydrostatic pressure induced three-dimensional Dirac semimetal in black phosphorus, Phys. Rev. B 93, 195434 (2016).
  * (4) J. Guo, H. Wang, F. von Rohr, W. Yi, Y. Zhou, Z. Wang, S. Cai, S. Zhang, X. Li, Y. Li, J. Liu, K. Yang, A. Li, S. Jiang, Q. Wu, T. Xiang, R. J. Cava, and L. Sun, Electron-hole balance and the anomalous pressure-dependent superconductivity in black phosphorus, Phys. Rev. B 96, 224513 (2017).
  * (5) P. Bridgman, Two new modifications of phosphorus,  J. Am. Chem. Soc. 36, 1344 (1914).
  * (6) R. W. Keyes, The electrical properties of black phosphorus,  hy. Rev. 92, 580 (1953).
  * (7) J. Kim, S. S. Baik, S. H. Ryu, Y. Sohn, S. Park, B.-G. Park, J. Denlinger, Y. Yi, H. J. Choi, and K. S. Kim, Observation of tunable band gap and anisotropic Dirac semimetal state in black phosphorus, Science 349, 723 (2015).
  * (8) J. Kim, S. Baik, S. Jung, Y. Sohn, S. Ryu, H. Choi, B. Yang, and K. Kim, Two-Dimensional Dirac Fermions Protected by Space-Time Inversion Symmetry in Black Phosphorus, Phys. Rev. Lett. 119, 226801 (2017).
  * (9) A. Carvalho, M. Wang, X. Zhu, A. Rodin, H. Su, and A. C. Neto, Phosphorene: from theory to applications, Nat. Rev. Materials 1, Article number: 16061 (2016).
  * (10) A. N. Rudenko and M. I. Katsnelson, Quasiparticle band structure and tight-binding model for single-and bilayer black phosphorus, Phys. Rev. B 89, 201408(R) (2014).
  * (11) J. M. Pereira Jr. and M. I. Katsnelson, Landau levels of single layer and bilayer phosphorene, Phys. Rev. B 92, 075437 (2015).
  * (12) Q. Wei and X. Peng, Superior mechanical flexibility of phosphorene and few-layer black phosphorus, Appl. Phys. Lett. 104, 251915 (2014).
  * (13) Y. Ren, P. Liu, F. Cheng and G. Zhou, Strain-induced effects in zigzag-edged blue phosphorene nanoribbons with edge sulfur passivation, J. Phys.: Condens. Matter 30, 395303 (2018).
  * (14) H. Hedayati and E. Faizabadi, Dwell time, Hartman effect and transport properties in a ferromagnetic phosphorene monolayer, J. Phys.: Condens. Matter 30, 085303 (2018).
  * (15) C. Heikes, I.-L. Liu, T. Metz, C. Eckberg, P. Neves, Y. Wu, L. Hung, P. Piccoli, H. Cao, J. Leao, J. Paglione, T. Yildirim, N. P. Butch, W. Ratcliff, Mechanical control of crystal symmetry and superconductivity in Weyl semimetal $\rm MoTe_{2}$, arXiv:1804.09093.
  * (16) Li, L. et al. Black phosphorus field-effect transistors, Nat. Nanotechnol. 9, 372 (2014).
  * (17) H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tomanek, and P. D. Ye, Phosphorene: an unexplored 2D semiconductor with a high hole mobility, ACS NANO 8, 4033 (2014).
  * (18) Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, Electronics and optoelectronics of two-dimensional transition metal dichalcogenides, Nat. Nanotech. 7, 699 (2012).
  * (19) W. Wu, et al. Piezoelectricity of single-atomic-layer $\rm MoS_{2}$ for energy conversion and piezotronics,  Nat. 514, 470 (2014).
  * (20) V. K. Sangwan, et al. Gate-tunable memristive phenomena mediated by grain boundaries in single-layer $\rm MoS_{2}$ ,  Nature Nanotech. 10, 403–406 (2015).
  * (21) H. Kawamura, I. Shirotani, and K. Tachikawa, Anomalous superconductivity and pressure induced phase transitions in black phosphorus, Solid State Commun. 54, 775 (1985).
  * (22) M. Karuzawa, M. Ishizuka, and S. Endo, The pressure effect on the superconducting transition temperature of black phosphorus  J. Phys.: Condens. Matter 14, 10759 (2002).
  * (23) J.A. Flores-Livas, A. Sanna, A. P. Drozdov, L. Boeri, G. Profeta, M. Eremets, and S. Goedecker, Interplay between structure and superconductivity: Metastable phases of phosphorus under pressure,  Phys. Rev. Mater. 1, 024802 (2017).
  * (24) B.-T. Wang, P.-F. Liu, T. Bo, W. Yin, F. Wang, Superconductivity in two-dimensional phosphorus carbide ($\rm\beta_{0}$-PC), arXiv:1803.03399.
  * (25) I. Shirotani, J. Mikami, T. Adachi, Y. Katayama, K. Tsuji, H. Kawamura, O. Shimomura, and T. Nakajima, Phase transitions and superconductivity of black phosphorus and phosphorus-arsenic alloys at low temperatures and high pressures, Phys. Rev. B 50, 16274 (1994).
  * (26) J. Flores-Livas, A. Sanna, A. Drozdov, L. Boeri, G. Profeta, M. Eremets, and S. Goedecker, Interplay between structure and superconductivity: Metastable phases of phosphorus under pressure, Phys. Rev. Materials 1, 024802 (2017).
  * (27) G. Q. Huang, Z. W. Xing, and D. Y. Xing, Prediction of superconductivity in Li-intercalated bilayer phosphorene, Appl. Phys. Lett. 106, 113107 (2015).
  * (28) J. Zhang and S. Dong, Prediction of above 20 K superconductivity of blue phosphorus bilayer with metal intercalations, 2D Mater. 3, 035006 (2016).
  * (29) Y. Feng, H. Sun, J. Sun, Z. Lu and Y. You, Prediction of phonon-mediated superconductivity in hole-doped black phosphorus, J. Phys.: Condens. Matter 30 015601 (2018).
  * (30) R. Zhang, J. Waters, A.K. Geim and I.V. Grigorieva, Intercalant-independent transition temperature in superconducting black phosphorus, Nat. Comm. 8, Article number: 15036 (2017).
  * (31) H. Yuan, L. Deng, B. Lv, Z. Wu, Z. Yang, S. Li, S. Huyan, Y. Ni, J. Sun, F. Tian, D. Wang, H. Wang, S. Chen, Z. Ren, C. Chu, Investigation on the reported superconductivity in intercalated black phosphorus, arXiv:1712.06716.
  * (32) T. Ma, F. Yang, H. Yao, and H.-Q. Lin, Possible triplet $p+ip$ superconductivity in graphene at low filling, Phys. Rev. B 90, 245114 (2014).
  * (33) J. P. L. Faye, P. Sahebsara, and D. Sénéchal, Chiral triplet superconductivity on the graphene lattice, Phys. Rev. B 92, 085121 (2015).
  * (34) R. Nandkishore, L. S. Levitov and A. V. Chubukov, Chiral superconductivity from repulsive interactions in doped graphene, Nat. Phys. 8, 158 (2012).
  * (35) M. L. Kiesel, C. Platt, W. Hanke, D. A. Abanin, and R. Thomale, Competing many-body instabilities and unconventional superconductivity in graphene, Phys. Rev. B 86, 020507(R) (2012).
  * (36) B. Uchoa and A. H. Castro Neto, Superconducting states of pure and doped graphene,  Phys. Rev. Lett. 98, 146801 (2007).
  * (37) P. M. Koenraad, M. E. Flatte, Single dopants in semiconductors,  Nat. Mater., 10, 91 (2011).
  * (38) R. A. Doganov, E. C. O’Farrell, S. P. Koenig, Y. Yeo, A. Ziletti, A. Carvalho, D. K. Campbell, D. F. Coker, K. Watanabe, T. Taniguchi, Transport properties of pristine few-layer black phosphorus by van der Waals passivation in an inert atmosphere, Nat. Comm. 6, 6647 (2015).
  * (39) C. Han, Z. Hu, A. Carvalho, N. Guo, J. Zhang, F. Hu, D. Xiang, J. Wu, B. Lei, L. Wang, Oxygen induced strong mobility modulation in few-layer black phosphorus, 2D Mater. 4, 021007 (2017).
  * (40) S. P. Koenig, R. A. Doganov, H. Schmidt, A. H. C. Neto, B. Olzyilmaz, Electric field effect in ultrathin black phosphorus, Appl. Phys. Lett. 104, 103106 (2014).
  * (41) Y.-L. Zou, J. Song, C. Bai, K. Chang, Impurity scattering and Friedel oscillations in monolayer black phosphorus, Phys. Rev. B 94, 035431 (2016).
  * (42) Y. Liu, T. Low, P. P. Ruden Mobility anisotropy in monolayer black phosphorus due to scattering by charged impurities, Phys. Rev. B 93, 165402 (2016).
  * (43) Z. Qiu, H. Fang, A. Carvalho, A. S. Rodin, Y. Liu, S. Jun R. Tan, M. Telychko, P. Lv, J. Su, Y. Wang, A. H. Castro Neto, and J. Lu, Resolving the spatial structures of bound hole states in black phosphorus, Nano Lett. 17, 6935 (2017).
  * (44) A. Zyuzin, M. Alidoust, D. Loss, Josephson junction through a disordered topological insulator with helical magnetization, Phys. Rev. B 93, 214502 (2016).
  * (45) I. V. Bobkova, A. M. Bobkov, A. A. Zyuzin, M. Alidoust, Magnetoelectrics in disordered topological insulator Josephson junctions, Phys. Rev. B 94, 134506 (2016).
  * (46) M. Alidoust and H. Hamzehpour, Spontaneous supercurrent and $\varphi_{0}$ phase shift parallel to magnetized topological insulator interfaces, Phys. Rev. B 96, 165422 (2017).
  * (47) F. N. Xia, H. Wang, and Y. C. Jia, Rediscovering black phosphorus as an anisotropic layered material for optoelectronics and electronics, Nat. Commun. 5, 4458 (2014).
  * (48) L. Voon, A. Lopez-Bezanilla, J. Wang, Y. Zhang, M. Willatzen, Effective Hamiltonians for phosphorene and silicene, New J. Phys. 17, 025004 (2015).
  * (49) L. Voon, J. Wang, Y. Zhang, and M. Willatzen, Band parameters of phosphorene, Journal of Physics: Conference Series 633, 012042 (2015).
  * (50) M. Alidoust, M. Willatzen, and A.-P. Jauho, Strain-engineered Majorana zero energy modes and $\varphi_{0}$ Josephson state in black phosphorus, Phys. Rev. B 98, 085414 (2018).
  * (51) M. Alidoust and K. Halterman, Proximity induced vortices and long-range triplet supercurrents in ferromagnetic Josephson junctions and spin valves, J. Appl. Phys. 117, 123906 (2015).
  * (52) A. Srivastava, L. A. B. O. Olthof, A. Di Bernardo, S. Komori, M. Amado, C. Palomares-Garcia, M. Alidoust, K. Halterman, M. G. Blamire, and J. W. A. Robinson, Magnetization-control and transfer of spin-polarized Cooper pairs into a half-metal manganite, Phys. Rev. Applied 8, 044008 (2017).
  * (53) H. Haug and A.-P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors, Springer-Verlag Berlin Heidelberg (2008).
  * (54) G. Eilenberger, Transformation of Gorkov’s equation for type II superconductors into transport-like equations, Z. Phys. 214, 195 (1968).
  * (55) K. D. Usadel, Generalized diffusion equation for superconducting alloys, Phys. Rev. Lett. 25, 507 (1977).
  * (56) A. V. Zaitsev, Sov. Phys. JETP 59, 1015 (1984).
  * (57) M. Y. Kuprianov and V. F. Lukichev, Sov. Phys. JETP 67, 1163 (1988).
  * (58) M. Alidoust, A. Zyuzin, K. Halterman, Pure Odd Frequency Superconductivity at the Cores of Proximity Vortices, Phys. Rev. B 95, 045115 (2017).
  * (59) Z. P. Niu, A spin triplet supercurrent in half metal ferromagnet/superconductor junctions with the interfacial Rashba spin-orbit coupling, Appl. Phys. Lett. 101 062601 (2012).
  * (60) F. Konschelle, Transport equations for superconductors in the presence of spin interaction, Euro. Phys. J. B 87, 119 (2014).
  * (61) F. S. Bergeret and I. V. Tokatly, Spin-orbit coupling as a source of long-range triplet proximity effect in superconductor-ferromagnet hybrid structures, Phys. Rev. B 89, 134517 (2014).
  * (62) M. Alidoust, K. Halterman, Spontaneous edge accumulation of spin currents in finite-size two-dimensional diffusive spin-orbit coupled SFS heterostructures, New J. Phys. 17, 033001 (2015).
  * (63) M. Alidoust, K. Halterman, Long-range spin-triplet correlations and edge spin currents in diffusive spin-orbit coupled SNS hybrids with a single spin-active interface, J. Phys: Cond. Matt. 27, 235301 (2015).
  * (64) M. Minutillo, D. Giuliano, P. Lucignano, A. Tagliacozzo, and G. Campagnano, Anomalous Josephson effect in S/SO/F/S heterostructures, Phys. Rev. B 98, 144510 (2018).
  * (65) N. Satchell and N. O. Birge, Supercurrent in ferromagnetic Josephson junctions with heavy metal interlayers, Phys. Rev. B 97, 214509 (2018).
  * (66) H. Chakraborti, S. Deb, R. Schott, V. Thakur, A. Chatterjee, S. Yadav, R. K. Saroj, A. Wieck, S. M. Shivaprasad, K. Das Gupta and S. Dhar, Coherent transmission of superconducting carriers through a $\sim 2\mu m$ polar semiconductor, Supercond. Sci. Technol. 31 085007 (2018).