# Surface supercurrent diode effect

Noah F. Q. Yuan<EMAIL_ADDRESS>School of Science, Harbin Institute of
Technology, Shenzhen 518055, China

###### Abstract

We propose a new type of supercurrent diode effect on the surface of a
superconductor with surface states under in-plane magnetic fields. Surface
supercurrent diode effect can lead to a perfect supercurrent diode in a
considerably wide range of fields. For comparison, the conventional
supercurrent diode effect due to the spin-orbit coupling in a two-dimensional
superconductor cannot be perfect in usual cases. Candidates such as the (001)
surface of iron-based superconductors BaFe2-xCoxAs2 are discussed.
Calculations are performed under the Ginzburg-Landau formalism.

When both inversion and time-reversal symmetries are broken, the critical
current of a superconductor can be nonreciprocal, known as the supercurrent
diode effect (SDE). To theoretically elaborate the mechanism of SDE and
experimentally demonstrate SDE turn out to be not any easy task.

Experimentally, nonreicprocal transport in superconducting systems has been
found in the fluctuating region Waka ; Yasu ; Qin ; Hoshi , where resistance
is in general finite and nonlinear transport is significantly enhanced
compared to the normal conducting state. Recently, Ando et. al. Ando  realized
a superconducting diode that has zero resistance in one direction but finite
in the opposite direction, which inspired tremendous discoveries of SDEs in
various superconducting systems, such as two-dimensional electron gas (2DEG)
in quantum wells Chris , transition metal dichalcogenides Banabir ; Lorenz ,
twisted bilayer graphene Diez  and twisted trilayer graphene LinJ .

As far as we know, SDE was first theoretically proposed by Victor M. Edelstein
EdelEE ; EdelEE1 ; EdelME ; EdelSDE . After the discovery of the famous
Edelstein effect EdelEE ; EdelEE1 , Edelstein proposed the magnetoelectric
effect EdelME  and the supercurrent diode effect EdelSDE  for polar
superconductors. As will be elaborated Section I, Edelstein’s physical picture
is correct while the detailed calculation was in fact inadequate to show SDE
due to the gauge symmetry.

Recently theories of SDE were developed Yuan ; Akito ; James ; Harley ; Zhai ;
Ilic , including the one by Liang Fu and the author Yuan . Such theories apply
to two-dimensional (2D) superconductors and can be unified in the Ginzburg-
Landau formalism, where Cooper pairs are boosted to finite momentum by
external magnetic fields, so that critical currents parallel and anti-parallel
to the Cooper pair momentum can be unequal. In order to capture the critical
current nonreciprocity, the free energy is expanded to the quartic order in
Cooper pair momentum $q$. In this work, we refer to this theory as the
conventional theory of SDE, and this paper is devoted to address two issues
related to it.

Figure 1: Supercurrent diode coefficient $\eta$ under in-plane magnetic field
$H$ and temperature $T$ in the surface supercurrent diode effect, when field
and current are perpendicular. The solid blue line denotes the in-plane
critical field $H_{c3\parallel}$ for surface superconductivity, and the dashed
orange line denotes the field strength $H_{\rm M}$ beyond which $\eta=1$ (see
maintext).

First, the microscopic mechanism of the conventional SDE theory should be more
elaborated. As elaborated in Sec. I, we show that the Edelstein effect is the
microscopic driving force of the conventional supercurrent diode effect in
two-dimensional superconductors with spin-orbit coupling (SOC). Second,
mechanisms other than the conventional theory should be explored for SDE. In
Sec. II, we consider a three-dimensional (3D) superconductor with surface
states, and find that the Meissner effect of surface states can also lead to
SDE on its surface.

In this work, we find the conventional supercurrent diode effect (CSDE)
carried by conventional 2DEG with SOC and the surface supercurrent diode
effect (SSDE) carried by surface states of a 3D superconductor. In both types
of SDE, the critical currents $J_{c}^{\pm}$ along opposite directions can be
unequal, quantified by the dimensionless supercurrent diode coefficient

$\displaystyle\eta\equiv\frac{J_{c}^{+}-J_{c}^{-}}{J_{c}^{+}+J_{c}^{-}}\in[-1,1].$
(1)

In Fig. 1 and Fig. 2 the supercurrent diode coefficients of SSDE and CSDE are
plotted in the field-temperature plane respectively, whose mechanisms and
properties will be elaborated in the following two Sections.

##  I Conventional supercurrent diodes

###  I.1 Phenomenological theory

The free energy of a 2D superconductor per area is

$\displaystyle f_{2{\rm
D}}=\int\frac{d^{2}\bm{q}}{(2\pi)^{2}}\left\\{\sigma|\Delta_{\bm{q}}|^{2}+\frac{1}{2}\gamma|\Delta_{\bm{q}}|^{4}\right\\},$
(2)

up to the quartic order of the order parameter $\Delta_{\bm{q}}$ with Cooper
pair momentum $\bm{q}$, with quadratic and quaric coefficients $\sigma,\gamma$
respectively. In order to describe SDE, the quartic expansion of $\sigma$ is
carried out

$\displaystyle\sigma$ $\displaystyle=$ $\displaystyle
a+b(\hat{\bm{n}}\times\bm{H})\cdot\bm{q},$ (3) $\displaystyle a$
$\displaystyle=$ $\displaystyle t+a_{0}q^{2}-a_{1}q^{4},$ $\displaystyle b$
$\displaystyle=$ $\displaystyle b_{0}-b_{1}q^{2},$

where $\bm{H}$ is the in-plane external magnetic field, $\hat{\bm{n}}$ is the
normal vector of the 2D superconductor, $t=(T-T_{c})/T_{c}$ is the reduced
temperature with critical temperature $T_{c}$, and $a_{0,1},b_{0,1},\gamma$
are coefficients to be given later. The expansion of $\gamma$ in terms of $q$
may also be carried out Ilic , which leads to higher order corrections.

By minimizing $f$, equilibrium Cooper pairs are found to have the finite
momentum boosted by $\bm{H}$

$\bm{q}_{0}=\frac{b_{0}}{2a_{0}}\bm{H}\times\hat{\bm{n}}.$ (4)

In the presence of a gauge potential $\bm{A}$, the free energy is changed upon
minimal coupling $\bm{q}\to\bm{q}-2e\bm{A}$, and the supercurrent line density
can be computed $\bm{J}=-\partial f_{2{\rm D}}/\partial\bm{A}$. In the 2D
limit under in-plane fields we take the limit $\bm{A}\to\bm{0}$, the
supercurrent thus reads
$\bm{J}=\int\frac{d^{2}\bm{q}}{(2\pi)^{2}}\bm{J}_{\bm{q}}$, where
$\bm{J}_{\bm{q}}=e|\sigma|\partial_{\bm{q}}\sigma/\gamma$ is a function of
$\bm{q}$, defined in the range such that $\sigma\leq 0$.

Although $q$-linear term is sufficient to describe finite-momentum
superconductivity at zero current, higher order $q$-terms are needed for the
critical current to be nonreciprocal. When
$\sigma=t+a_{0}q^{2}+b_{0}(\hat{\bm{n}}\times\bm{H})\cdot\bm{q}$ is truncated
at $q^{2}$-term, although $\sigma(\bm{q})\neq\sigma(-\bm{q})$ is not an even
function of $\bm{q}$, it is found that
$\sigma(\bm{q}_{0}+{\bm{q}})=\sigma(\bm{q}_{0}-{\bm{q}})$ with respect to the
equilibrium momentum $\bm{q}_{0}$. As a result,
$\bm{J}_{\bm{q}_{0}+{\bm{q}}}=-\bm{J}_{\bm{q}_{0}-{\bm{q}}}$. Such symmetry
guarantees the critical current reciprocity $J_{c}^{+}=J_{c}^{-}$ and
$\eta=0$.

To describe SDE, we consider the following expansion near the equilibrium
momentum

$\displaystyle\sigma(\bm{q}+\bm{q}_{0})=a-b^{\prime}_{1}q^{2}(\hat{\bm{n}}\times\bm{H})\cdot\bm{q},\quad
b^{\prime}_{1}=b_{1}-2\frac{a_{1}}{a_{0}}b_{0},$ (5)

and find that the finite-momentum Cooper pairs boosted by external field would
change the odd order coefficients $b_{0}\to 0,b_{1}\to b^{\prime}_{1}$. Due to
this modification, up to the linear order in $\bm{H}$, the supercurrent diode
coefficient reads

$\eta=\left(\frac{b_{1}}{b_{0}}-2\frac{a_{1}}{a_{0}}\right)\sqrt{\frac{|t|}{3a_{0}^{3}}}b_{0}{(\bm{H}\times\hat{\bm{n}})\cdot\hat{\bm{i}}},$
(6)

where $\hat{\bm{i}}$ is the current direction. Details of the calculations can
be found in the Appendix.

Figure 2: Supercurrent diode coefficient $\eta$ under in-plane magnetic field
$B$ and temperature $T$ in the conventional supercurrent diode effect of a
Rashba superconductor, when the supercurrent and the magnetic field are
perpendicular.

###  I.2 Microscopic theory

Next we compute the Ginzburg-Landau coefficients $a_{0,1},b_{0,1},\gamma$ in
Eq. (3) from a microscopic model of the Rashba superconductor Yuan . Under an
in-plane magnetic field $\bm{H}$, the orbital effect is screened, and the
Zeeman effect dominates. Denote $c_{s}$ as the electron with spin
$s=\uparrow,\downarrow$, then on the Nambu basis
$(c_{\uparrow},c_{\downarrow},-c^{\dagger}_{\downarrow},c^{\dagger}_{\uparrow})$,
the Bogouliubov-de Gennes (BdG) Hamiltonian reads

$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}H_{\bm{k}+\frac{1}{2}\bm{q}}&\Delta_{\bm{q}}\\\
\Delta_{\bm{q}}^{*}&-\mathcal{T}H_{-\bm{k}+\frac{1}{2}\bm{q}}\mathcal{T}\end{pmatrix},$
(7) $\displaystyle H_{\bm{k}}$ $\displaystyle=$
$\displaystyle\frac{|\bm{k}|^{2}}{2m}-\varepsilon_{\rm F}+\alpha_{\rm
R}(\bm{k}\times\bm{\sigma})\cdot\hat{\bm{n}}+\mu_{\rm
m}\bm{H}\cdot\bm{\sigma},$ (8)

where $H_{\bm{k}}$ describes the normal state electrons with momentum
$\bm{k}$, spin $\bm{\sigma}$, mass $m$, Fermi energy $\varepsilon_{\rm F}$,
Rashba spin-orbit coupling (SOC) strength $\alpha_{\rm R}$ and the magnetic
moment $\mu_{\rm m}$. Electrons form Cooper pairs with momentum $\bm{q}$ under
the pairing potential $\Delta_{\bm{q}}$, and $\mathcal{T}=i\sigma_{y}K$ is the
antiunitary time-reversal operator with complex conjugation $K$.

To the linear order in $\bm{q}$, we have

$\mathcal{H}=H_{0}\tau_{z}+\Delta_{\bm{q}}\tau_{x}+\mu_{\rm
m}\bm{H}\cdot\bm{\sigma}+\frac{\bm{k}\cdot\bm{q}}{2m}+\frac{\alpha_{\rm
R}}{2}(\bm{q}\times\bm{\sigma})\cdot\hat{\bm{n}}.$ (9)

The zero-field normal Hamiltonian $H_{0}=H_{\bm{k}}|_{\bm{H}=\bm{0}}$ gives
rise to two energy bands $\xi_{\pm}={|\bm{k}|^{2}}/({2m})-\varepsilon_{\rm
F}\pm\alpha_{\rm R}|\bm{k}|$ split by SOC. Correspondingly the two Fermi
surfaces $\xi_{\pm}=0$ will have the same Fermi velocity $v_{\rm
F}=\sqrt{2\varepsilon_{\rm F}/m+\alpha_{\rm R}^{2}}$ but different density of
states (DOS) $N_{\pm}=N_{0}(1\mp{\alpha_{\rm R}}/{v_{\rm F}})$ with
$N_{0}=m/(2\pi)$ EdelME ; Samokhin ; Dimi1 ; Dimi2 .

There are three depairing terms, from the Zeeman effect $\mu_{\rm
m}\bm{H}\cdot\bm{\sigma}$, the Doppler effect $\bm{k}\cdot\bm{q}/(2m)$ of the
kinetic energy, and most importantly the Edelstein effect
$\frac{1}{2}{\alpha_{\rm R}}(\bm{q}\times\bm{\sigma})\cdot\hat{\bm{n}}$ of the
SOC EdelEE ; EdelEE1 , respectively.

The even order coefficients $a_{0,1}$ and the coefficient $\gamma$ can be
derived as in the conventional BCS theory

$\displaystyle a_{0}=\frac{1}{4}\frac{C_{0}v_{\rm F}^{2}}{(\pi
T_{c})^{2}},\;a_{1}=\frac{1}{8}\frac{C_{1}v_{\rm F}^{4}}{(\pi
T_{c})^{4}},\;\gamma=\frac{C_{0}}{(\pi T_{c})^{2}},$ (10)

where $C_{0}=7\zeta(3)/8$, $C_{1}=93\zeta(5)/2^{8}$ are numerical constants,
$v_{\rm F}=\sqrt{2\varepsilon_{\rm F}/m+\alpha_{\rm R}^{2}}$ is the Fermi
velocity. The odd order coefficients $b_{0,1}$ are determined by SOC, through
both band splitting effect and Edelstein effect:

$\displaystyle b_{0}=\frac{3}{2}\frac{C_{0}\alpha_{\rm R}}{(\pi
T_{c})^{2}}\mu_{\rm m},\;b_{1}=\frac{5}{2}\frac{C_{1}v_{\rm F}^{2}\alpha_{\rm
R}}{(\pi T_{c})^{4}}\mu_{\rm m}.$ (11)

From the coefficients above, we can work out the conventional supercurrent
diode coefficient

$\eta={0.64}\frac{\alpha_{\rm R}}{v_{\rm
F}}\frac{(\bm{H}\times\hat{\bm{n}})\cdot\hat{\bm{i}}}{H_{P}}\sqrt{1-\frac{T}{T_{c}}},$
(12)

where $H_{P}=1.25T_{c}/\mu_{\rm m}$ is the Pauli limit. This formula has been
numerically verified in Ref. Yuan  with the correct numerical factor.

There have been a series of works in deriving the correct formula for the SDE.
As far as we know, Victor M. Edelstein himself was one of the pioneers to
theoretically study the possible SDE in a superconductor without inversion
symmetry under an external magnetic field EdelSDE . However, the expansion of
$\sigma$ was truncated at $q^{2}$ term, and hence was not able to describe the
CSDE as elaborated above. Nevertheless, in Ref. EdelSDE , the Ginzburg-Landau
coefficient $b_{0}$ is the same as that in Eq. (11).

Liang Fu and the author studied the SDE of a Rashba superconductor (which we
dub as CSDE) in Ref. Yuan  by both numerical and analytical methods. The
numerical model includes the full effect of Zeeman field, the Doppler effect
and the Edelstein effect, hence numerical calculations lead to the
supercurrent diode coefficient consistent with Eq. (12). However, in the
analytical calculations, the Edelstein effect was absent and the the expansion
of $\sigma$ was truncated at $q^{3}$ term only. Interestingly, the
analytically derived supercurrent diode coefficient was still the same as Eq.
(12), in spite of the absence of two ingredients in the expansion of free
energy ($q^{4}$ term and Edelstein term).

Recently, S. Ilić and F. S. Bergeret further studied the SDE Ilic  by expand
$\sigma$ up to $q^{4}$ term. It was then found that the supercurrent diode
coefficient was zero to the linear order of the field without the Edelstein
effect.

The Edelstein effect is in fact the driving force of the CSDE. Due to the
Edelstein effect, Cooper pairs with momentum $\bm{q}$ will also carry spin
polarization $\bm{S}\propto\hat{\bm{n}}\times\bm{q}$ EdelME . The external
magnetic field $\bm{H}$ then couples to Cooper pairs via the Zeeman effect
$\bm{S}\cdot\bm{H}\propto(\hat{\bm{n}}\times\bm{H})\cdot\bm{q}$ and hence
favors a particular direction for $\bm{q}$ and for the supercurrent, leading
to the CSDE. To explicitly demonstrate this, one can calculate the odd order
coefficients without the Edelstein effect, which leads to smaller values
$b_{0}=C_{0}{\alpha_{\rm R}}\mu_{\rm m}/{(\pi T_{c})^{2}},\;b_{1}=C_{1}{v_{\rm
F}^{2}\alpha_{\rm R}}\mu_{\rm m}/{(\pi T_{c})^{4}}$ compared with those in Eq.
(11). Since $b_{1}/b_{0}=2a_{1}/a_{0}$ without the Edelstein effect, the
supercurrent diode coefficient vanishes $\eta\equiv 0$ according to Eq. (6).
Moreover, the critical field diverges at zero temperature when the Edelstein
effect is not included Samokhin ; Dimi1 ; Dimi2 . To have nonzero supercurrent
diode effect at finite temperatures and finite critical field at zero
temperature, it would be necessary to include the Edelstein effect in a
consistent microscopic theory of the two-dimensional Rashba superconductors.

##  II Surface supercurrent diodes

###  II.1 Ginzburg-Landau formalism

To improve the supercurrent diode coefficient, we can consider strong-SOC
materials. In particular, we can consider the strongest limit of SOC
$\alpha_{\rm R}\to v_{\rm F}$, which corresponds to the surface states
described by a single Dirac cone Hamiltonian $H_{\bm{k}}=v_{\rm
F}(\bm{k}\times\bm{\sigma})_{z}$, so that $\alpha_{\rm R}=v_{\rm F}$.

However, a single Dirac cone cannot exist by itself, but can live with the
whole bunch of bulk states. As a result, to describe the supercurrent diode
effect of superconducting surface states, one may include the bulk
superconductor as well.

We consider a semi-infinite superconductor in the lower half space $z\leq 0$,
which its normal conducting phase possesses surface states. In the Ginzburg-
Landau (GL) regime, the temperature and the field are near the superconducting
phase transitions and the length scales are greater than the zero-temperature
bulk coherence length. Since surface states usually extend to the bulk by a
few atomic layers TI , in the GL regime the surface states only live on the
surface $z=0$ of the superconductor. After washing out short-range details,
the GL free energy is composed of the bulk term and the surface term

$F=\int_{V}f_{3{\rm D}}d^{3}\bm{r}+\int_{\partial V}f_{2{\rm D}}d^{2}\bm{r},$
(13)

where the bulk energy density is

$f_{3{\rm
D}}=\alpha|\psi|^{2}+\frac{\beta}{2}|\psi|^{4}+\frac{|\bm{D}\psi|^{2}}{2m}+\frac{1}{2}(\bm{B}-\bm{H})^{2},$
(14)

and the surface energy density is

$f_{2{\rm D}}=\sigma|\psi|^{2}+\frac{1}{2}\gamma|\psi|^{4}.$ (15)

Here, $\psi(\bm{r})$ is the superconducting order parameter,
$\bm{D}=\nabla-2ie\bm{A}$ is the covariant derivative, and $\bm{A}(\bm{r})$ is
the vector potential. Due to the Meissner effect in 3D superconductors,
besides the internal field $\bm{B}=\nabla\times\bm{A}$, we also introduce the
external field $\bm{H}$.

The bulk free energy density $f_{3{\rm D}}$ is characterized by three
parameters $m,\alpha$ and $\beta$, where $m$ is the electron mass,
$\alpha=\alpha_{0}(T-T_{c0})$ is a function of temperature $T$, and $T_{c0}$
is the bulk critical temperature. The stability of the bulk free energy
requires $m,\alpha_{0},\beta$ to be all positive.

The surface free energy density $f_{2{\rm D}}$ is characterized by the
quadratic surface energy $\sigma$ Noah ; Simonin ; Andryushin ; deGennes ;
saint  and quartic term $\gamma$ as described in Section I. As the bulk free
energy describes the second-order phase transition between two phases (i.e.
superconducting versus normal), in the boundary free energy the quartic term
is irrelevant Diehl  and hence can be neglected in the rest of this
manuscript.

When surface states are present, the surface energy is negative $\sigma<0$,
suggesting the preferential accumulation of Cooper pairs on the surface. As a
result, the superconducting phase transition consists of two steps. The first
step is at temperatures above the bulk critical temperature, known as the
superconducting nucleation, where the electrons condense into Cooper pairs in
the surface layer, while the bulk electrons stay mostly normal. In the second
step, the temperature is lower than the bulk critical temperature, the bulk is
superconducting as well as the surface layer.

By minimizing the GL free energy with respect to $\psi$, we derive the GL
equation

$\displaystyle\left(-\frac{\bm{D}^{2}}{2m}+\alpha+\beta|\psi|^{2}\right)\psi=0,\
(\bm{r}\in V)$ (16)

and the de Gennes boundary condition Noah ; Simonin ; Andryushin ; deGennes ;
saint

$\displaystyle\hat{\bm{n}}\cdot\bm{D}\psi=\frac{\psi}{l_{s}},\
(\bm{r}\in\partial V)$ (17)

where $l_{s}=1/(2m|\sigma|)>0$ is the extrapolation length denoting the
thickness of surface superconductivity, and $\hat{\bm{n}}$ is the normal
vector of the surface pointing outward.

By minimizing the GL free energy with respect to $\bm{A}$, the Ampere’s law is
derived with supercurrent density $\bm{j}$

$\displaystyle\nabla\times\bm{B}=\bm{j}\equiv\frac{2e}{m}{\rm
Im}(\psi^{*}\bm{D}\psi)\ (\bm{r}\in V).$ (18)

According to the de Gennes boundary condition in Eq. (17), no bulk
supercurrent penetrates through the surface
$\hat{\bm{n}}\cdot\bm{j}|_{\partial V}=0$, since
$\psi^{*}\hat{\bm{n}}\cdot\bm{D}\psi$ is real on the surface.

###  II.2 Surface superconductivity

Near the superconducting phase transition, $|\psi|$ will be small and the GL
equation Eq. (16) can be linearized. We then find the solutions to the
linearized GL equation to obtain properties on the superconducting phase
transitions. Details on the critical temperature and critical field at zero
current can be found in Ref. Noah . We review the main results below.

At zero field and zero current, the surface superconductivity is confined
within the surface layer of thickness $\sim l_{s}$, and the onset critical
temperature for the surface superconductivity is
$T_{c}=T_{c0}+({2m\alpha_{0}l_{s}^{2}})^{-1}$, higher than the bulk critical
temperature. Near $T_{c}$, we find

$\psi=\phi\exp\left(-\frac{|z|}{l_{s}}\right),$ (19)

where $\phi=\sqrt{{2\alpha_{0}(T_{c}-T)}/{\beta}}$ is determined by the
nonlinear GL equation as shown in the Appendix.

Under finite fields and zero current, superconductivity is in general
suppressed by the fields. It is well known that the upper critical field for
the bulk is $H_{c2}=\Phi_{0}/(2\pi\xi_{0}^{2})\propto(T_{c0}-T)^{1}$, where
$\xi_{0}(T)\equiv\sqrt{1/(2m|\alpha|)}$ is the bulk coherence length and
$\Phi_{0}=h/(2e)$ is the flux quantum Abrikosov1 . When $H>H_{c2}$, bulk
superconductivity is killed, while surface superconductivity can still survive
until a new critical field, which is the critical field $H_{c3}$ of the
surface superconductivity. For out-of-plane fields
$H_{c3\perp}=\Phi_{0}/(2\pi\xi^{2})\propto(T_{c}-T)^{1}$, where the coherence
length is $\xi\equiv 1/\sqrt{2m\alpha_{0}(T_{c}-T)}$, while for in-plane
fields
$H_{c3\parallel}=0.525{\Phi_{0}}/{\left(\xi^{2}l_{s}\right)^{\frac{2}{3}}}\propto(T_{c}-T)^{2/3}$
Noah .

In this work, we consider weak in-plane magnetic fields. In conventional
superconductors we would expect Meissner effect, where magnetic fields are
repelled from the bulk. As a result, the magnetic field can only penetrate
into the surface layer, and supercurrent can only flow on the surface layer as
well. The thickness of this surface layer is another fundamental length scale
of the superconductor known as the London penetration depth
$\lambda_{L}\equiv(2e)^{-1}\sqrt{m\beta/|\alpha|}\propto|T_{c0}-T|^{-1/2}$.

However, in our case of superconductors with surface states, when
$T_{c0}<T<T_{c}$ superconductivity only survives in the surface layer of
thickness $\sim l_{s}$ already, and the Meissner effect is different. To see
this, we plug in the order parameter Eq. (19) of the surface superconductivity
into the Ampere’s law Eq. (18), then near $T_{c}$

$\bm{B}=\bm{H}\exp\left(-\frac{2|z|}{l_{s}}\right).$ (20)

It can be seen that the London penetration depth $\lambda_{L}$ does not appear
in the spatial distribution of the internal magnetic field $\bm{B}$ at least
to the leading order. Instead, the magnetic field fully penetrates into the
surface layer of thickness $\sim l_{s}$ where surface superconductivity lives.
Deep in the bulk, both magnetic field and superconductivity vanish. This is
the Meissner effect in the case of superconductors with surface states.

To be concrete, one can calculate that at zero external current but finite
field $\bm{H}$, surface Cooper pairs have zero momentum, while a persistent
current $\bm{J}_{0}=\int_{-\infty}^{0}\bm{j}dz=\hat{\bm{n}}\times\bm{H}$ will
be induced to screen magnetic field in the bulk. Notice that the supercurrent
line density $\bm{J}_{0}$ has the same dimension as the magnetic field
$\bm{H}$. At finite external current and field, surface Cooper pairs can have
finite momentum $\bm{q}$, and the supercurrent line density is

$\displaystyle\bm{J}=\int_{-\infty}^{0}\bm{j}dz=\left(I_{0}\bm{q}+\bm{J}_{0}\right)(1-\xi^{2}q^{2}),$
(21)

where $\bm{q}$ is confined in the range $q<1/\xi$ such that the system remains
superconducting, and

$I_{0}=\frac{l_{s}}{\kappa^{2}}H_{c3\perp},\quad\bm{J}_{0}=\hat{\bm{n}}\times\bm{H}.$
(22)

Notice that $H_{c3\perp}=\Phi_{0}/(2\pi\xi^{2})$ is the out-of-plane critical
field of the surface superconductivity with the coherence length $\xi\equiv
1/\sqrt{2m\alpha_{0}(T_{c}-T)}$ and the bulk GL parameter
$\kappa=\lambda_{L}/\xi_{0}$, and $\bm{J}_{0}$ is the field-induced current at
$\bm{q}=\bm{0}$.

We can then calculate the critical currents parallel and anti-parallel to a
given direction $\hat{\bm{i}}$ by maximizing and minimizing
$\bm{J}\cdot\hat{\bm{i}}$ over $\bm{q}$ in the range $q<1/\xi$ respectively.
Namely, $J_{c}^{+}={\rm max}(\bm{J}\cdot\hat{\bm{i}})$ and $J_{c}^{-}=-{\rm
min}(\bm{J}\cdot\hat{\bm{i}})$. At weak field, we find
$J_{c}^{\pm}=J_{c}(1\pm\eta)$, the average is
$J_{c}=\frac{2\sqrt{3}}{9}I_{0}/\xi$, and the surface supercurrent diode
coefficient reads

$\displaystyle\eta=\sqrt{3}\kappa^{2}\frac{(\bm{H}\times\hat{\bm{n}})\cdot\hat{\bm{i}}}{H_{c3\perp}}\frac{\xi}{l_{s}}.$
(23)

This surface supercurrent diode coefficient $\eta$ is plotted in the phase
plane of magnetic field $H$ and temperature $T$, with supercurrent
perpendicular to the field, as shown in Fig. 1. At weak fields, $\eta\propto
H(T_{c}-T)^{-3/2}$ and $J_{c}\propto H(T_{c}-T)^{3/2}$, hence the critical
current difference $J_{c}^{+}-J_{c}^{-}\propto H$ is temperature independent.
At higher fields when $H\geq H_{\rm M}\equiv(l_{s}/\xi)H_{c3\perp}/\kappa^{2}$
we find $\eta=1$ until the critical field $H_{c3\parallel}$ where surface
superconductivity vanishes. A perfect supercurrent diode hence could be
realized in the field regime $H_{\rm M}<H<H_{c3\parallel}$ on the surface of a
superconductor with surface states, as shown in Fig. 1 where the Meissner
field $H_{\rm M}\propto(T_{c}-T)^{3/2}$ is denoted by the dashed orange line
and the critical field of surface superconductivity
$H_{c3\parallel}\propto(T_{c}-T)^{2/3}$ is denoted by the solid blue line.

For further lower temperatures $T<T_{c0}$, the surface together with the bulk
are all superconducting. Since the surface states contribution now becomes
negligible compared with the bulk, the supercurrent becomes reciprocal deep in
the bulk superconducting state. Without surface states, the surface energy is
zero $\sigma\to 0$, $l_{s}\to\infty$, and the critical current also becomes
reciprocal according to Eq. (23) Abrikosov2 . Thus to realize SSDE, surface
states are needed on the surface of a superconductor and the temperature range
is close to the onset critical temperature so that only surface
superconductivity exists.

###  II.3 Candidates

We have considered an ideal model of the semi-infinite superconductor with
surface states on its only one surface (the top surface). However, realistic
materials will at least have two surfaces, the top and the bottom surfaces. We
suppose surface states exist on both top and bottom surfaces, with
extrapolation lengths ${l_{s}^{\rm top}}$ and ${l_{s}^{\rm bot}}$
respectively, and states on the top surface are decoupled from the bottom
surface. Without loss of generality, we assume $l_{s}^{\rm top}<l_{s}^{\rm
bot}$, then the onset critical temperature for the surface superconductivity
is $T_{c}=T_{c0}+({2m\alpha_{0}|l_{s}^{\rm top}|^{2}})^{-1}$. In the
temperature regime $T_{c1}<T<T_{c}$ with
$T_{c1}=T_{c0}+({2m\alpha_{0}|l_{s}^{\rm bot}|^{2}})^{-1}$, only the top
surface is superconducting while the bulk and the bottom surface are not. Eq.
(23) then applies in this temperature regime. For lower temperatures
$T_{c0}<T<T_{c1}$, both the top and the bottom surfaces are superconducting
while the bulk is not, and the total surface supercurrent diode coefficient
from these two surfaces is $\eta=\eta(l_{s}^{\rm top})-\eta(l_{s}^{\rm bot})$.

The iron-based superconductors FeSe1-xTex ($x\sim 0.5$) Gang ; Zhang ; Wang ;
Kong  provide candidates of superconductors with surface states and hence of
SSDE. It has been theoretically calculated that Gang , the $3d$ electrons from
Fe atoms will exhibit band inversion in the bulk spectrum, and topological
surface states would emerge on the (001) surfaces (top and bottom) as observed
experimentally Zhang ; Wang ; Kong . However, due to the topological
protection, both the top and bottom surfaces will host surface states, and one
cannot eliminate states on one surface while remaining another by local
operations. Therefore, in usual samples of FeSe1-xTex, we expect $l_{s}^{\rm
top}\approx l_{s}^{\rm bot}$, and hence the SSDE would be weak, unless the top
and bottom surfaces can be made sharply different.

Notice that our theory of SSDE in fact does not specify the type and origin of
the surface states. Surface states in SSDE may have various origins, such as
topological surface states in FeSe1-xTex Gang ; Zhang ; Wang ; Kong  as
discussed above, surface reconstructions in cleaved or covered superconductors
Gao ; Heumen ; Par  as will be discussed below, and other crystal defect
states Matsu ; Khl .

We thus turn to superconductors with surface reconstruction, which fall into
another family of iron-based superconductors, the doped iron arsenides
BaFe2-xCoxAs2 Gao ; Heumen ; Par . The surface Ba atoms predominantly favor a
$\sqrt{2}\times\sqrt{2}$ order Gao , which leads to the lattice reconstruction
of the surface layer. As a result, surface related Fe 3$d$ states are present
in the electronic structure in the form of surface bands as both theoretically
calculated and experimentally observed Heumen . By changing the surface
condition upon cleavage or coverage, one can change the reconstruction and
hence the onset critical temperature as observed experimentally Par . Compared
with topological surface states in FeSe1-xTex, the surface reconstruction in
BaFe2-xCoxAs2 can be well adjusted on just one surface. Consequently, the SSDE
in BaFe2-xCoxAs2 is expected to be significant upon proper surface conditions.

Finally we address the dependence of supercurrent diode coefficient $\eta$ on
field $\bm{H}$ and current $\hat{\bm{i}}$ in the two types of SDEs in this
work, which turns out to be dictated by the underlying point group of the
superconducting system. Notice that $\eta$ is time-reversal even while
$\bm{H}$ and $\hat{\bm{i}}$ are time-reversal odd, thus we expect $\eta$ to be
expressed as a quadratic form of $\bm{H}$ and $\hat{\bm{i}}$, which is
invariant under the point group. On the surface of a three-dimensional
superconductor or in the plane of a two-dimensional Rashba superconductor, the
point groups are of the same type, namely $C_{nv}$ or its subgroup, where
$n=2,3,4,6,\infty$. When $n>2$, $\bm{H}$ and $\hat{\bm{i}}$ furnish the same
2D irreducible representation, while the normal vector $\hat{\bm{n}}$
furnishes the 1D irreducible representation $A_{2}$. In order to furnish the
trivial representation we find
$\eta\propto({\bm{H}}\times\hat{\bm{n}})\cdot\hat{\bm{i}}$ as shown in Eqs.
(12) and (23). In a two-dimensional superconductor without the inversion
center, however, other quadratic invariants will become possible, since the
point group is not limited to $C_{nv}$ Yuan ; Yuan1 . For example
$\eta\propto\bm{H}\cdot\hat{\bm{i}}$ when the point group is $D_{n}$
($n=2,3,4,6$), and
$\eta\propto\bm{H}\cdot\hat{\bm{i}}-2(\bm{H}\cdot\hat{\bm{m}})(\hat{\bm{m}}\cdot\hat{\bm{i}})$
when the point group is $D_{2d}$, where $\hat{\bm{m}}$ is the normal vector of
the mirror plane.

##  III Conclusions

In this work, we first show that in a two-dimensional superconductor with
spin-orbit coupling, the Edelstein effect makes supercurrent spin-polarized
and hence an external Zeeman field will favor the supercurrent along
particular directions, leading to the conventional supercurrent diode effect.
Second, we find that in three-dimensional superconductors with surface states,
close to the superconducting-normal phase transitions, superconductivity can
only survive in the surface layer. An external in-plane field can generate a
persistent current in the surface layer and hence affect the critical current
of the surface superconductivity, resulting in the surface supercurrent diode
effect. In the phase plane spanned by the temperature and the field, the
supercurrent diode coefficient of CSDE decreases when approaching the critical
field line (Fig. 2), while that of SSDE increases close to the
superconducting-normal phase transition (Fig. 1). Moreover, a perfect
supercurrent diode can be realized in SSDE near the critical field regime,
where currents along one direction can be supercurrent when smaller than the
critical value, while along the opposite direction no current can become
superconducting. Candidates for SSDE are discussed.

Acknowledgement—This work is supported by the National Natural Science
Foundation of China (Grant. No. 12174021).

##  Appendix A Surface supercurrent diode effect

###  A.1 Mean-field derivations of Ginzburg-Landau free enrgy

We consider the model Hamiltonian at zero field with effective on-site
attractive interaction among electrons

$H=\int
d^{3}\bm{r}d^{3}\bm{s}c^{\dagger}(\bm{r})\mathcal{H}(\bm{r},\bm{s})c(\bm{s})-\frac{g}{2}\int
d^{3}\bm{r}\left\\{c^{\dagger}(\bm{r})c(\bm{r})\right\\}^{2},$ (24)

where $c=\\{c_{\mu}\\}^{\rm T}$, $c_{\mu}(\bm{r})$ denotes an electron at site
$\bm{r}$ with band index $\mu$, $\mathcal{H}(\bm{r},\bm{s})$ is the normal
state Hamiltonian matrix of the topological material, and $g>0$ is the on-site
attraction strength between electrons. We would like to compute the $s$-wave
pairing correlation $\Delta(\bm{r})=g\langle
c_{\uparrow}(\bm{r})c_{\downarrow}(\bm{r})\rangle$ in order to obtain the
physics of superconductivity in this system, where $\langle O\rangle\equiv{\rm
Tr}(Oe^{-\frac{H}{k_{\rm B}T}})/Z$ denotes the thermodynamic average at
temperature $T$, and $Z\equiv{\rm Tr}e^{-\frac{H}{k_{\rm B}T}}$ is the
partition function. To capture the long-range physics of $\Delta(\bm{r})$, we
expand the free energy $F\equiv{-{k_{\rm B}T}}\log Z$ within mean field

$F=\int d^{3}\bm{r}\frac{|\Delta(\bm{r})|^{2}}{g}-\int
d^{3}\bm{r}d^{3}\bm{s}K(\bm{r},\bm{s})\Delta^{*}(\bm{r})\Delta(\bm{s}),$ (25)

where $K$ is the kernel in terms of sum over Mastubara frequency
$\omega=(2n+1)\pi k_{\rm B}T$ $(n\in\mathbb{Z})$:

$K(\bm{r},\bm{s})=k_{\rm B}T\sum_{\omega
ab}\frac{\phi_{a}(\bm{r})^{\dagger}\phi_{a}(\bm{s})}{\xi_{a}-i\omega}\left[\frac{\phi_{b}(\bm{r})^{\dagger}\phi_{b}(\bm{s})}{\xi_{b}-i\omega}\right]^{*},$
(26)

and eigenstates $\mathcal{H}\phi_{a}=\xi_{a}\phi_{a}$ of the normal
Hamiltonian. Deep in the bulk $|\bm{r}|\to\infty$, the kernel only involves
bulk states $\phi^{\rm B}$ and hence has full translation symmetry

$\displaystyle K^{\rm B}(\bm{r},\bm{s})$ $\displaystyle=$ $\displaystyle
k_{\rm B}T\sum_{\omega ab}\frac{\phi^{\rm B}_{a}(\bm{r})^{\dagger}\phi^{\rm
B}_{a}(\bm{s})}{\xi_{a}-i\omega}\left[\frac{\phi^{\rm
B}_{b}(\bm{r})^{\dagger}\phi^{\rm B}_{b}(\bm{s})}{\xi_{b}-i\omega}\right]^{*}$
(27) $\displaystyle=$ $\displaystyle K^{\rm B}(\bm{r}-\bm{s}).$

As a result, the free energy is composed of the bulk term and the surface term
$F=F_{0}+F_{s}$

$\displaystyle F_{0}$ $\displaystyle=$ $\displaystyle\int
d^{3}\bm{r}\frac{|\Delta(\bm{r})|^{2}}{g}-\int d^{3}\bm{r}d^{3}\bm{s}K^{\rm
B}(\bm{r}-\bm{s})\Delta^{*}(\bm{r})\Delta(\bm{s}),$ $\displaystyle F_{s}$
$\displaystyle=$ $\displaystyle\int
d^{3}\bm{r}d^{3}\bm{s}\Delta^{*}(\bm{r})G(\bm{r}),$ $\displaystyle G$
$\displaystyle=$ $\displaystyle\int d^{3}\bm{s}[K^{\rm
B}(\bm{r}-\bm{s})-K(\bm{r},\bm{s})]\Delta(\bm{s}).$

To rewrite the free energy within GL formalism, we introduce the order
parameter

$\psi(\bm{r})={\sqrt{mL}}\Delta(\bm{r}),\quad L=\frac{1}{3}\int K^{\rm
B}(\bm{r})|\bm{r}|^{2}d^{3}\bm{r},$ (28)

then the bulk free energy becomes the conventional GL free energy

$\displaystyle F_{0}=\int
d^{3}\bm{r}\left\\{\frac{|\nabla\psi|^{2}}{2m}+\alpha|\psi|^{2}\right\\},$
(29) $\displaystyle\alpha=\frac{1}{mL}\left\\{\frac{1}{g}-\int K^{\rm
B}(\bm{r})d^{3}\bm{r}\right\\}.$ (30)

From the bulk free energy above, the bulk critical temperature $T_{c0}$ and
the bulk coherence length $\xi_{00}$ at zero temperature and zero field are
determined by

$g\left.\int K^{\rm
B}(\bm{s})d^{3}\bm{s}\right|_{T=T_{c0}}=1,\quad\xi_{00}=0.18\frac{\hbar
v_{F}}{k_{\rm B}T_{c0}},$ (31)

where $v_{F}$ is the bulk Fermi velocity.

We consider the semi-infinite superconductor $\Sigma=\\{(x,y,z)|z\leq 0\\}$,
and the surface $\partial\Sigma=\\{(x,y,0)\\}$ contains surface states
$\phi^{\rm S}$ whose size of extension along $z$-axis is much smaller than
$\xi_{00}$. As a result, the translation symmetries along in-plane $x,y$
directions are preserved, while the translation symmetry along $z$-axis is
broken by the planar defect of the surface $z=0$. In this case, along $z$-axis
we need to perform the Laplace transform

$\displaystyle\hat{G}(\bm{\rho},p)=\int_{-\infty}^{0}G(\bm{\rho},z)e^{pz}dz$
(32)

In the following we suppress $\bm{\rho}$ and write
$\hat{G}(\bm{\rho},p)=\hat{G}(p),G(\bm{\rho},z)=G(z)$ for short. In the GL
regime, we focus on the long-range physics, which is $p\to 0$ in the Laplace
domain $\hat{G}(p)\approx\hat{G}(0)$. Then we perform the inverse Laplace
transform and get

$\displaystyle
G(z)=\mathcal{L}^{-1}\\{\hat{G}(p)\\}=\hat{G}(0)\mathcal{L}^{-1}\\{1\\}=\hat{G}(0)\delta(z).$
(33)

Thus the surface free energy reads

$\displaystyle F_{s}=\int
d^{3}\bm{r}\Delta^{*}(\bm{r})G(\bm{r})=\int_{z=0}d^{2}\bm{\rho}\Delta^{*}(\bm{\rho})\hat{G}(0),$
(34)

where the zeroth order Laplace transform is just the integral over the lower
$z$-axis

$\displaystyle\hat{G}(0)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{0}G(\bm{\rho},z)dz$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{0}dz\int d^{3}\bm{s}[K^{\rm
B}(\bm{\rho}-\bm{s})-K(\bm{\rho},\bm{s})]\Delta(\bm{s}).$

In terms of density of states $N,N_{0}$, bulk critical temperature $T_{c0}$,
and Debye frequency $\omega_{D}$, we have

$\displaystyle\int{d^{3}\bm{s}}K(\bm{r},\bm{s})$ $\displaystyle=$
$\displaystyle N(\bm{r})\log\frac{1.14\hbar\omega_{D}}{k_{\rm B}T_{c0}},$ (36)
$\displaystyle\int{d^{3}\bm{s}}K^{\rm B}(\bm{r}-\bm{s})$ $\displaystyle=$
$\displaystyle N_{0}\log\frac{1.14\hbar\omega_{D}}{k_{\rm
B}T_{c0}}=\frac{1}{g}.$

Here we define the norm squared $|\phi|^{2}\equiv{\phi^{\dagger}\phi}$, then
$N(\bm{r})=\sum_{a}|\phi_{a}(\bm{r})|^{2}\delta(\xi_{a})$ is the local density
of states, and $N_{0}=\underset{z\to-\infty}{\lim}N(\bm{r})$ is the bulk
density of states. Finally we arrive at the expression of the surface free
energy in the GL regime

$\displaystyle F_{s}$ $\displaystyle=$ $\displaystyle\int
d^{2}\bm{r}U(\bm{r})|\psi(\bm{r})|^{2},$ (37) $\displaystyle U(\bm{r})$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{0}\frac{\Delta(\bm{r},z)}{\Delta_{0}}\left\\{1-\frac{N(\bm{r},z)}{N_{0}}\right\\}\frac{dz}{mgL}.$

where $\Delta_{0}=\Delta(\bm{r})|_{z=0}$ is the pairing potential on the
surface. In the clean limit, $N(\bm{r},z)$ is a constant in the $xy$-plane on
the scale of bulk coherence length, and $U(\bm{r})=\sigma$ is spatially
constant.

By minimal coupling principle, we replace $\nabla$ by covariant derivative
$\bm{D}=\nabla-2ie\bm{A}$ and arrive at the GL free energy

$\displaystyle F$ $\displaystyle=$ $\displaystyle F_{0}+F_{s},$ (38)
$\displaystyle F_{0}$ $\displaystyle=$
$\displaystyle\int_{V}d^{3}\bm{r}\left\\{\frac{|\bm{D}\psi|^{2}}{2m}+\alpha|\psi|^{2}+\frac{\beta}{2}|\psi|^{4}+\frac{1}{2}(\bm{B}-\bm{H})^{2}\right\\},$
$\displaystyle F_{s}$ $\displaystyle=$ $\displaystyle{\sigma}\int_{\partial
V}|\psi|^{2}d^{2}\bm{r}.$

###  A.2 Order parameter spatial profiles

By minimizing the GL free energy with respect to $\psi$, we derive the GL
equation

$\displaystyle\left(-\frac{\bm{D}^{2}}{2m}+\alpha+\beta|\psi|^{2}\right)\psi=0,\
(\bm{r}\in V)$ (39)

and the de Gennes boundary condition

$\displaystyle\hat{\bm{n}}\cdot\bm{D}\psi=\frac{\psi}{l_{s}},\
(\bm{r}\in\partial V)$ (40)

where $l_{s}=1/(2m|\sigma|)>0$ is the extrapolation length denoting the
thickness of surface superconductivity.

Now we solve the nonlinear GL equation under the generalized de Gennes
boundary condition in the absence of external magnetic fields $\bm{A}=\bm{0}$.
We introduce the bulk order parameter $\psi_{0}\equiv\sqrt{|\alpha|/\beta}$,
the bulk coherence length $\xi_{0}\equiv\sqrt{1/(2m|\alpha|)}$, and the
dimensionless temperature $t=(T-T_{c0})/(T_{c}-T_{c0})$, then the integral to
the above equation reads

$-\xi_{0}^{2}\psi_{0}^{2}\left(\frac{d\psi}{dz}\right)^{2}+(\psi^{2}+\psi_{0}^{2}{\rm
sgn}\alpha)^{2}=C,$ (41)

where the integral constant $C$ depends on the temperature $T$.

In the regime $T_{c0}<T<T_{c}$, most of the superconductor could not be
superconducting except the surface layer. Deep in the bulk ($z\to-\infty$),
$\psi\to 0$, $\alpha>0$ and hence $C=\psi_{0}^{4}>0$. The order parameter is
solved as

$\psi(z)=\sqrt{2}\psi_{0}{{\rm
csch}\left(\frac{|z|}{\xi_{0}}+\theta\right)},\quad(T_{c0}<T<T_{c})$ (42)

where $\theta={\rm arctanh}\sqrt{|t|}$. At temperatures $T\lesssim T_{c}$
and/or in the long range $|z|\gg\xi$, superconductivity decays exponentially
$\psi(z)\sim\exp\left(-{|z|}/{\xi_{0}}\right)$ with bulk coherence length as
the decay length. In particular when $T=T_{c}$, $\xi_{0}=l_{s}$,
superconductivity is localized $\psi\propto\exp(z/{l_{s}})$ in the surface
layer with thickness $\sim l_{s}$. At $T=T_{c0}$, the surface
superconductivity reduces to power-law-like decay into the bulk
$\psi(z)\propto({{|z|}/{l_{s}}+1})^{-1}.$

In the regime $T<T_{c0}$, however, the bulk is already superconducting, and
$\psi$ is mostly uniform except the surface layer. Deep in the bulk
($z\to-\infty$), $\psi\to\psi_{0}$, $\alpha<0$ and hence $C=0$. The order
parameter is

$\psi(z)={\psi_{0}}{\coth\left(\frac{|z|}{\sqrt{2}\xi_{0}}+\varphi\right)},\quad(T<T_{c0})$
(43)

where $\varphi=\frac{1}{2}{\rm arcsinh}\sqrt{2|t|}$. Below $T_{c0}$ and in the
long range $|z|\gg\xi$, superconductivity is asymptotically uniform
$\psi(z)\approx\psi_{0}$.

When $T\to T_{c}$, $t\to 1$ and Eq. (42) becomes Eq. (19)

$\psi(z)=\phi\exp\left(-\frac{|z|}{l_{s}}\right),\quad\phi=\sqrt{\frac{2\alpha_{0}(T_{c}-T)}{\beta}}.$
(44)

###  A.3 Meissner effect of field

By minimizing the GL free energy with respect to $\bm{A}$, the Ampere’s law is
derived

$\displaystyle\nabla\times\bm{B}=\bm{j},\ (\bm{r}\in V)$ (45)

with the Ampere’s continuity boundary condition
$\hat{\bm{n}}\times\left(\bm{H}-\bm{B}\right)=\bm{0},\ (\bm{r}\in\partial V)$
and the supercurrent density

$\bm{j}=\frac{2e}{m}{\rm Im}(\psi^{*}\bm{D}\psi).$ (46)

Now we solve the Ampere’s law under the Ampere’s continuity boundary condition
in the presence of an external in-plane magnetic field $\bm{H}$. We choose the
gauge $\bm{A}=\bm{H}\times\hat{\bm{n}}f(z)$, then we have the Ampere’s law

$f^{\prime\prime}(z)=\frac{2(1-t)}{\kappa^{2}l_{s}^{2}}\exp\left(\frac{2z}{l_{s}}\right)f(z)$
(47)

with boundary condition $f^{\prime}(0)=1$. The solution is found

$f(z)=A_{0}^{-1}{I_{0}\left[2\delta\exp\left(\frac{z}{l_{s}}\right)\right]},\
A_{0}=\delta I_{1}(2\delta),$ (48)

where $\delta\equiv{\sqrt{2(1-t)}}/{\kappa}$ is a dimensionless quantity
depending on temperature.

The internal magnetic field is

$\displaystyle\bm{B}$ $\displaystyle=$
$\displaystyle\nabla\times\bm{A}=\bm{H}f^{\prime}(z),$ (49) $\displaystyle
f^{\prime}(z)$ $\displaystyle=$
$\displaystyle\frac{\delta}{A_{0}l_{s}}\exp\left(\frac{z}{l_{s}}\right){I_{1}\left[2\delta\exp\left(\frac{z}{l_{s}}\right)\right]}.$
(50)

Notice that $\delta\equiv{\sqrt{2(1-t)}}/{\kappa}$ is small near $T_{c}$ in
type-II superconductors, we have

$A_{0}f(z)=1+\delta^{2}\exp\left(\frac{2z}{l_{s}}\right)+O(\delta^{4}).$ (51)

Hence we find the internal magnetic field reads

$\bm{B}=\bm{H}\left\\{e^{2z/l_{s}}+\frac{\delta^{2}}{2}e^{2z/l_{s}}\left(e^{2z/l_{s}}-1\right)\right\\}+O(\delta^{4}).$
(52)

###  A.4 Surface supercurrent diode effect

We choose the ansatz $\psi=\Delta\exp(z/l_{s}+i\bm{q}\cdot\bm{r})$ with Cooper
pair momentum $\bm{q}=(q_{x},q_{y})$ to account for finite supercurrent. To
determine the parameter $\Delta$, we need to take into account the nonlinear
part of the GL equation, or equivalently we calculate the total free energy up
to $|\Delta|^{4}$. By minimizing the total free energy $F=F_{0}+F_{s}$, we
find $|\Delta|^{2}=-2\tilde{\alpha}_{\bm{q}}/\beta$, where
$\tilde{\alpha}_{\bm{q}}=\alpha_{0}(T-T_{c})+|\bm{q}|^{2}/(2m)$. Thus the
supercurrent line density becomes

$\displaystyle\bm{J}=\int_{-\infty}^{0}\bm{j}dz=I_{0}\left(\bm{q}-\bm{q}_{0}\right)(1-\xi^{2}q^{2}),$
(53)

where $\bm{q}_{0}=-\bm{J}_{0}/I_{0}=\bm{H}\times\hat{\bm{n}}/I_{0}$. The
dimensionless supercurrent diode coefficient is a function of the
dimensionless quantity $x\equiv\xi\bm{q}_{0}\cdot\hat{\bm{i}}$

$\eta=\begin{cases}{x(9-x^{2})}{(3+x^{2})^{-3/2}}&|x|\leq 1\\\ {\rm
sgn}(x)&|x|>1\end{cases}.$ (54)

At weak fields $|x|\ll 1$, $\eta\approx\sqrt{3}x$ and we recover Eq. (23).

Figure 3: Supercurrent diode coefficient $\eta$ as functions of the
dimensionless quantity. a) In the surface supercurrent diode effect (SSDE),
the dimensionless quantity is $\xi\bm{q}_{0}\cdot\hat{\bm{i}}$. At weak fields
$|x|\ll 1$, $\eta\approx\sqrt{3}x$ as shown by the green dashed line. b) In
the conventional supercurrent diode effect (CSDE), the dimensionless quantity
is $\xi_{0}^{2}/\xi\bm{q}_{0}\cdot\hat{\bm{i}}$. Green dashed lines are linear
approximation and orange dashed lines denote the perfect diode case $\eta=1$.

##  Appendix B Conventional supercurrent diode effect

###  B.1 Ginzburg-Landau coefficients

We compute the the function $\sigma$ microscopically within mean field theory.
The free energy is $F=\int f(\bm{q})d^{2}\bm{q}$, and the free energy density
can be expressed by the BdG Hamiltonian $f=|\Delta|^{2}/V-T\int
d^{2}\bm{k}{\rm tr}\log[1+e^{-\mathcal{H}/T}]$, where $V$ is the attractive
interaction strength and tr denotes trace in spin space.

By expansion in terms of $\Delta_{\bm{q}}$ we have

$\displaystyle\sigma=N_{0}\log\frac{T}{T_{c}}+\sum_{\lambda=\pm}\oint_{{\rm
FS}_{\lambda}}\frac{dk}{|\bm{v}^{\lambda}|}\left\\{\phi\left(\frac{Q+\lambda\mathcal{E}_{+}}{2\pi
T}\right)\cos^{2}\frac{\theta}{2}+\phi\left(\frac{Q+\lambda\mathcal{E}_{-}}{2\pi
T}\right)\sin^{2}\frac{\theta}{2}\right\\},$ (55)

where $\lambda=\pm$ denote contributions from inner ($\lambda=+$) and outer
($\lambda=-$) Fermi surfaces respectively, $\phi(x)={\rm
Re}\left[\psi\left(\frac{1+ix}{2}\right)\right]-\psi\left(\frac{1}{2}\right),$
and $\psi$ is the digamma function.

Here, $\bm{v}^{\lambda}=\partial_{\bm{k}}\xi_{\bm{k}}^{\lambda}$ is the
electron velocity, $Q=\bm{v}^{\lambda}\cdot\bm{q}$ is the depairing energy of
finite momentum pairing, $\mathcal{E}_{\pm}=|\bm{h}_{+}|\pm|\bm{h}_{-}|$ is
the depairing energy of Zeeman splitting for inter- ($+$) or intra-pocket
($-$) Cooper pairs with
$\bm{h}_{\pm}=\bm{B}+\bm{g}_{\frac{1}{2}\bm{q}\pm\bm{k}}$, and the angle
$\theta=\langle\bm{h}_{+},\bm{h}_{-}\rangle$ between $\bm{h}_{\pm}$ controls
the ratio between inter- or intra-pocket Cooper pairs. Supercurrent affects
$\bm{q}$ and hence depairing energy $Q$, while magnetic field $\bm{B}$
together with SOC affects depairing energies $\mathcal{E}_{\pm}$ and angle
$\theta$.

Near $T_{c}$, the temperature dependence of $\sigma_{\bm{q}}$ can be captured
by the first term $N_{0}\log(T/T_{c})$, and we can set $T=T_{c}$ in the Fermi
surface integrals. To evaluate the integral, notice that the field is weak
$H\ll H_{P}$, one can expand the special function

$\displaystyle\phi(x)=2.10x^{2}-2.01x^{4}+2.00x^{6}+O(x^{8})$ (56)

and then integrate order by order to obtain GL coefficients.

###  B.2 Nonreciprocal critical current and polarity-dependent critical field

When the expansion of $\sigma$ near its minimum $\bm{q}_{0}$ reads

$\displaystyle\sigma=\sigma_{0}+a\delta q_{\parallel}^{2}-b\delta
q_{\parallel}^{3},$ (57)

where $a>0$ and $\delta
q_{\parallel}=(\bm{q}-\bm{q}_{0})\cdot\hat{\bm{q}}_{0}$. We find the
supercurrent is

$\displaystyle\gamma
J_{\parallel}/e=|\sigma_{\bm{q}}|\partial_{\parallel}\sigma_{\bm{q}}=2a\sigma_{0}\delta
q_{\parallel}-3\left(\sigma_{0}b\right)\delta q_{\parallel}^{2}$
$\displaystyle+2a^{2}\delta q_{\parallel}^{3}-5(ab)\delta
q_{\parallel}^{4}+3b^{2}\delta q_{\parallel}^{5}.$ (58)

Notice that to the leading order of $|\sigma_{0}|^{\frac{1}{2}}$, critical
currents $\pm J_{c}^{\pm}$ correspond to $\delta q_{\parallel}=\mp\delta
q_{c}$ respectively, where $\delta q_{c}=\sqrt{|\sigma_{0}|/3a}$. Then we have

$\gamma J_{c}^{\pm}/e=\frac{4|\sigma_{0}|^{3/2}}{9a}\left(\sqrt{3a^{3}}\pm
b\sqrt{|\sigma_{0}|}\right)+O(|\sigma_{0}|^{5/2}).$ (59)

To include higher order contributions, the supercurrent diode coefficient is

$\eta=g\left(x\right),\quad x=\sqrt{\frac{|\sigma_{0}|}{a^{3}}}{b}.$ (60)

The special function is

$g(x)=\frac{J(x)-J(-x)}{J(x)+J(-x)},$ (61)
$J(x)=[Q(x)-\frac{3}{2}xQ(x)][-1+Q(x)^{2}-xQ(x)^{3}]$ (62)
$Q(x)=\frac{1}{3x}-\frac{\sqrt{n}}{2x}-\frac{1}{2}\sqrt{{\frac{4}{5\sqrt{n}}}\left(1-\frac{2}{27x^{2}}\right)-\frac{n}{x^{2}}+\frac{8}{15x^{2}}}$
(63) $n=\frac{2}{15}\left(z+1/{z}+\frac{4}{3}\right),$ (64)
$z=\sqrt[3]{t+\sqrt{t^{2}-1}},\quad t=\frac{135}{4}x^{4}-5x^{2}+1.$ (65)

When $|x|<\frac{2}{3\sqrt{3}}$, $\frac{22}{27}<t<1$ and $z$ is complex. Denote
$t=\cos\theta$, then $z=e^{i\theta/3}$ and
$n=\frac{4}{15}(\cos\frac{\theta}{3}+\frac{2}{3})$ is real. Since $g(x)\approx
x/\sqrt{3}$ we get the leading order contribution
$\eta=\sqrt{\frac{|\sigma_{0}|}{3a^{3}}}{b}$.

The expansion of $\sigma_{\bm{q}}$ near its minimum $\bm{q}_{0}$ in general
can be anisotropic

$\displaystyle\sigma_{\bm{q}+\bm{q}_{0}}=\sigma_{0}+a(1+\epsilon)q_{x}^{2}+a(1-\epsilon)q_{y}^{2}+2a\eta
q_{x}q_{y}$
$\displaystyle-(b_{1}q_{x}^{3}+b_{2}q_{y}^{3}+b_{3}q_{x}q_{y}^{2}+b_{4}q_{x}^{2}q_{y}),$

where $a>0$ and $\epsilon^{2}+\eta^{2}<1$ for stability. The supercurrent
diode coefficient for supercurrent $\bm{J}=J(\cos\theta,\sin\theta)$ can be
worked out as

$\displaystyle\eta=\sqrt{\frac{|\sigma_{0}|}{3a^{3}}}\frac{\left(\frac{b_{1}+b_{3}}{2}+\frac{b_{1}-b_{3}}{2}\cos
2\theta\right)\cos\theta+\left(\frac{b_{2}+b_{4}}{2}-\frac{b_{2}-b_{4}}{2}\cos
2\theta\right)\sin\theta}{(1+\epsilon\cos 2\theta+\eta\sin 2\theta)^{3/2}}.$
(66)

## References

  * (1) R. Wakatsuki, Y. Saito, S. Hoshino, Y. M. Itahashi, T. Ideue, M. Ezawa, Y. Iwasa, and N. Nagaosa, Nonreciprocal charge transport in noncentrosymmetric superconductors, Sci. Adv. 3, e1602390 (2017).
  * (2) F. Qin, W. Shi, T. Ideue, M. Yoshida, A. Zak, R. Tenne, T. Kikitsu, D. Inoue, D. Hashizume and Y. Iwasa, Superconductivity in a chiral nanotube, Nat. Commun. 8, 14465 (2017).
  * (3) Kenji Yasuda, Hironori Yasuda, Tian Liang, Ryutaro Yoshimi, Atsushi Tsukazaki, Kei S. Takahashi, Naoto Nagaosa, Masashi Kawasaki and Yoshinori Tokura, Nonreciprocal charge transport at topological insulator/superconductor interface, Nat. Commun. 10, 2734 (2019).
  * (4) S. Hoshino, R. Wakatsuki, K. Hamamoto, and N. Nagaosa, Nonreciprocal charge transport in two-dimensional noncentrosymmetric superconductors, Phys. Rev. B 98, 054510 (2018).
  * (5) Fuyuki Ando, Yuta Miyasaka, Tian Li, Jun Ishizuka, Tomonori Arakawa, Yoichi Shiota, Takahiro Moriyama, Youichi Yanase and Teruo Ono, Observation of superconducting diode effect, Nature 584, 373 (2020).
  * (6) Christian Baumgartner, Lorenz Fuchs, Andreas Costa, Simon Reinhardt, Sergei Gronin, Geoffrey C. Gardner, Tyler Lindemann, Michael J. Manfra, Paulo E. Faria Junior, Denis Kochan, Jaroslav Fabian, Nicola Paradiso and Christoph Strunk, Supercurrent rectification and magnetochiral effects in symmetric Josephson junctions, Nature Nanotechnology 17, 39 (2022).
  * (7) Banabir Pal, Anirban Chakraborty, Pranava K. Sivakumar, Margarita Davydova, Ajesh K. Gopi, Avanindra K. Pandeya, Jonas A. Krieger, Yang Zhang, Mihir Date, Sailong Ju, Noah Yuan, Niels B. M. Schröter, Liang Fu and Stuart S. P. Parkin, Josephson diode effect from Cooper pair momentum in a topological semimetal, Nature Physics 18, 1228 (2022).
  * (8) Lorenz Bauriedl, Christian Bäuml, Lorenz Fuchs, Christian Baumgartner, Nicolas Paulik, Jonas M. Bauer, Kai-Qiang Lin, John M. Lupton, Takashi Taniguchi, Kenji Watanabe, Christoph Strunk and Nicola Paradiso, Supercurrent diode effect and magnetochiral anisotropy in few-layer NbSe2, Nat. Comm. 13, 4266 (2022).
  * (9) J. Diez-Merida, A. Diez-Carlon, S. Y. Yang, Y.-M. Xie, X.-J. Gao, K. Watanabe, T. Taniguchi, X. Lu, K. T. Law, and Dmitri K. Efetov, Magnetic Josephson Junctions and Superconducting Diodes in Magic Angle Twisted Bilayer Graphene, arXiv:2110.01067.
  * (10) Jiang-Xiazi Lin, Phum Siriviboon, Harley D. Scammell, Song Liu, Daniel Rhodes, K. Watanabe, T. Taniguchi, James Hone, Mathias S. Scheurer and J.I.A. Li, Zero-field superconducting diode effect in small-twist-angle trilayer graphene, Nat. Phys. 18, 1221 (2022).
  * (11) Victor M. Edelstein, Characteristics of the Cooper pairing in two-dimensional noncentrosymmetric electron systems, Sov. Phys. JETP 68 (6), (1989).
  * (12) Victor M. Edelstein, Spin polarization of conduction electrons induced by electric current in two-dimensional asymmetric electron systems, Solid State Commun. 73, 233 (1990).
  * (13) Victor M. Edelstein, Magnetoelectric Effect in Polar Superconductors, Phys. Rev. Lett. 75, 2004 (1995).
  * (14) Victor M. Edelstein, The Ginzburg-Landau equation for superconductors of polar symmetry, J. Phys.: Condens. Matter 8, 339 (1996).
  * (15) Noah F. Q. Yuan and Liang Fu, Supercurrent diode effect and finite-momentum superconductors, PNAS 119 (15) e2119548119 (2022).
  * (16) Akito Daido, Yuhei Ikeda, and Youichi Yanase, Intrinsic Superconducting Diode Effect, Phys. Rev. Lett. 128, 037001 (2022).
  * (17) James Jun He, Yukio Tanaka, and Naoto Nagaosa, A phenomenological theory of superconductor diodes, New J. Phys. 24, 053014 (2022).
  * (18) Harley D Scammell, J I A Li and Mathias S Scheurer, Theory of zero-field superconducting diode effect in twisted trilayer graphene, 2D Mater. 9, 025027 (2022).
  * (19) B. Zhai, B. Li, Y. Wen, F. Wu, and J. He, Prediction of ferroelectric superconductors with reversible superconducting diode effect, Physical Review B 106, L140505 (2022).
  * (20) S. Ilić and F. S. Bergeret, Theory of the Supercurrent Diode Effect in Rashba Superconductors with Arbitrary Disorder, Phys. Rev. Lett. 128, 177001 (2022).
  * (21) V. Dimitrova, M. V. Feigel’man, Phase diagram of a surface superconductor in parallel magnetic field, JETP Lett. 78, 637 (2003).
  * (22) V. Dimitrova, M. V. Feigel’man, Theory of a two-dimensional superconductor with broken inversion symmetry, Phys. Rev. B 76, 014522 (2007).
  * (23) K. V. Samokhin, Upper critical field in noncentrosymmetric superconductors, Phys. Rev. B 78, 224520 (2008).
  * (24) C.-X. Liu, X.-L. Qi, H. Zhang, X. Dai, Z. Fang, and Zhang S.-C. Model Hamiltonian for topological insulators, Phys. Rev. B 82, 045122 (2010).
  * (25) J. Simonin, Surface term in the superconductive Ginzburg-Landau free energy: Application to thin films, Phys. Rev. B 33, 7830(R) (1986).
  * (26) E. A. Andryushin, V. L. Ginzburg, and A. P. Silin, Boundary conditions in the macroscopic theory of superconductivity, Phys.-Usp. 36, 854 (1993).
  * (27) P. G. de Gennes, Superconductivity of metals and alloys (W. A. Benjamin, New York, 1966).
  * (28) D. Saint-James, and P. G. de Gennes, Onset of superconductivity in decreasing fields, Phys. Lett. 7, 306 (1963).
  * (29) Noah F. Q. Yuan and X.-J. Chen, Critical field measure for topological, arXiv:2209.00816.
  * (30) H. W. Diehl, The theory of boundary critical phenomena, Int. J. Mod. Phys. B11: 3503 (1997).
  * (31) A. A. Abrikosov, On the magnetic properties of superconductors of the second group, Sov. Phys. JETP 5, 1174 (1957).
  * (32) Gang Xu, Biao Lian, Peizhe Tang, Xiao-Liang Qi, and Shou-Cheng Zhang, Topological superconductivity on the surface of Fe-based superconductors, Phys. Rev. Lett. 117, 047001 (2016).
  * (33) Peng Zhang, Koichiro Yaji, Takahiro Hashimoto, Yuichi Ota, Takeshi Kondo, Kozo Okazaki, Zhijun Wang, Jinsheng Wen, G. D. Gu, Hong Ding, and Shik Shin, Observation of topological superconductivity on the surface of iron-based superconductor, Science 360, 182 (2018).
  * (34) Dongfei Wang, Lingyuan Kong, Peng Fan, Hui Chen, Shiyu Zhu, Wenyao Liu, Lu Cao, Yujie Sun, Shixuan Du, John Schneeloch, Ruidan Zhong, Genda Gu, Liang Fu, Hong Ding, and Hong-Jun Gao, Evidence for Majorana bound states in an iron-based superconductor, Science 362, 333 (2018).
  * (35) Lingyuan Kong, Shiyu Zhu, Michał Papaj, Lu Cao, Hiroki Isobe, Wenyao Liu, Dongfei Wang, Peng Fan, Hui Chen, Yujie Sun, Shixuan Du, John Schneeloch, Ruidan Zhong, Genda Gu, Liang Fu, Hong-Jun Gao, and Hong Ding, Half-integer level shift of vortex bound states in an iron-based superconductor, Nat. Phys. 15, 1181 (2019).
  * (36) Miao Gao, Fengjie Ma, Zhong-Yi Lu, and Tao Xiang, Surface structures of ternary iron arsenides AFe2As2 (A=Ba, Sr, or Ca), Phys. Rev. B 81, 193409 (2010).
  * (37) Erik van Heumen, Johannes Vuorinen, Klaus Koepernik, Freek Massee, Yingkai Huang, Ming Shi, Jesse Klei, Jeroen Goedkoop, Matti Lindroos, Jeroen van den Brink, and Mark S. Golden, Existence, character and origin of surface-related bands in the high temperature iron pnictide superconductor BaFe2-xCoxAs2, Phys. Rev. Lett. 106, 027002 (2011).
  * (38) Christopher T. Parzyck , Brendan D. Faeth , Gordon N. Tam, Gregory R. Stewart, and Kyle M. Shen, Enhanced surface superconductivity in Ba(Fe0.95Co0.05)2As2, Appl. Phys. Lett. 116, 062601 (2020).
  * (39) M. Matsumoto, C. Belardinelli, and M. Sigrist, Upper critical field of the 3 Kelvin phase in Sr2RuO4, J. Phys. Soc. Jpn. 72, 1623 (2003).
  * (40) I. N. Khlyustikov, and A. I. Buzdin, Twinning-plane superconductivity, Adv. Phys. 36, 271 (1987).
  * (41) A. A. Abrikosov, Concerning surface superconductivity in strong magnetic fields, Sov. Phys. JETP 20, 480 (1965).
  * (42) Noah F. Q. Yuan and Liang Fu, Topological metals and finite-momentum superconductors, PNAS 118 (3), e2019063118 (2021).


# Radar Sensing with OTFS: Embracing ISI and ICI to Surpass the Ambiguity
Barrier

Musa Furkan Keskin1, Henk Wymeersch1, Alex Alvarado2
1Chalmers University of Technology, Sweden 2Eindhoven University of
Technology, The Netherlands

###### Abstract

Orthogonal time frequency space (OTFS) is a promising alternative to
orthogonal frequency division multiplexing (OFDM) in high-mobility beyond 5G
communications. In this paper, we consider the problem of radar sensing with
OTFS joint radar-communications waveform and derive a novel OTFS radar signal
model by explicitly taking into account the inter-symbol interference (ISI)
and inter-carrier interference (ICI) effects. On the basis of the new model,
we show how ISI and ICI phenomena can be turned into an advantage to surpass
the maximum unambiguous detection limits in range and velocity, arising in
existing OFDM and OTFS radar systems. Moreover, we design a generalized
likelihood ratio test based detector/estimator that can _embrace_ ISI and ICI
effects. Simulation results illustrate the potential of embracing ISI/ICI and
demonstrate its superior detection and estimation performance over
conventional baselines.

Index Terms– OTFS, OFDM, joint radar-communications, inter-symbol
interference, inter-carrier interference.

##  I Introduction

As 5G systems are being rolled out, the time has come to conceive and develop
Beyond 5G (B5G) communication systems. There are now several initiatives in
Europa, the USA, and Asia to define what B5G will be in terms of use cases and
requirements [1, 2, 3]. As with all previous generations, one requirement will
be a 10-fold increase in peak data rate. Different from previous generations
is the increased emphasis on joint radar and communication (JRC) [4], driven
not only by localization use cases but also the inherent geometric nature of
the wireless propagation channel [5, 6].

In the pursuit of higher data rates, lower latency, and higher sensing
accuracies, we have no choice but to consider larger carrier frequencies,
above the 24 GHz band in 5G, as this is where larger bandwidths are available
[7]. At lower frequencies, (despite intense competition) OFDM has remained the
communication waveform of choice, due to its robustness to multipath, its
straightforward integration with multi-antenna systems, and its high
flexibility in terms of power and rate allocation [8]. In addition, OFDM is
suitable for JRC with standard FFT-processing, in both mono-static and bi-
static configurations [9]. However, at B5G frequencies, OFDM is challenged by
several effects, which force us to consider alternative modulation formats:
OFDM suffers from a high peak-to-average power ratio (PAPR), leading to
reduced power efficiency, which becomes a limiting factor at high carriers.
Secondly, OFDM require frequent adaptation due to mobility and fading, which
would lead to prohibitive overheads at high carriers due to the short
coherence times [8]. Moreover, the robustness to multipath comes at a cost of
inserting a cyclic prefix (CP) between OFDM symbols, resulting in a rate loss,
to combat inter-symbol interference (ISI). Finally, for radar, OFDM is
sensitive to inter-carrier-interference (ICI), resulting from Doppler shifts
under high target velocities [10].

These perceived drawbacks of OFDM have renewed the interest in alternative
modulation formats, in particular with favorable properties in terms of JRC
performance. Orthogonal time frequency space (OTFS) has become a promising
candidate in this respect, as it has lower PAPR [11], requires less frequent
adaptation [12], has a much lower cyclic overhead [13], and can cope with much
larger Doppler shifts [14]. In contrast to OFDM’s 2D modulation in the
time–frequency domain, OTFS relies on the delay–Doppler domain. This implies
that a time-varying channel with constant Doppler will appear time-invariant
to OTFS [15]. Since radar detections are of the form of range (delay) and
velocity (Doppler) tuples, OTFS is a natural candidate for JRC, as evidenced
by recent activity in this area [16, 17, 18, 19]. Target detection and
estimation of their range and velocity is done by converting the received
time-domain signal back to the delay-Doppler domain. This leads to high
complexity, e.g., requiring iterative interference cancellation-based
processing due to significant side-lobe levels [17]. In addition, the
expression of the signal in the delay-Doppler domain is complicated (see,
e.g., [19, Eq. (12)] and [17, Eq. (11)]), making it difficult to derive
insights into the structure of the OTFS signal in terms of ISI and ICI
effects. Finally, the existing approaches are limited in terms of ISI and ICI
due to standard ambiguity limits in [17, 18, 19].

Figure 1: OTFS JRC system, where the delay-Doppler data
$\mathbf{X}^{\rm{DD}}\in\mathbb{C}^{N\times M}$ is converted to the time-
frequency domain, and then to the time-domain, where pulse shaping occurs and
a single CP is added to the entire frame. The signal $y(t)$ is sampled in the
time domain at $t=\ell T/N$ for $\ell=0,\ldots,NM-1$ and the frequency domain
at $f=n\Delta f/M$ for $n=0,\ldots,NM-1$.

In this paper, we pursue a different and computationally more appealing
approach where the time-domain observations are directly utilized to estimate
target parameters without transformation from time to delay-Doppler domain.
This reveals the explicit dependency of the observation on target delays and
Dopplers, which facilitates detector/estimator design, leading to a low-
complexity multi-target generalized likelihood ratio test (GLRT) detector. In
addition, the proposed time-domain processing approach provides insights into
the OTFS waveform structure and thus enables exploitation of the ISI and ICI
effects to surpass the standard ambiguity limits. The main contributions of
this paper can be summarized as follows:

  * •

We derive a novel time-domain signal model for OTFS radar by explicitly taking
into account the ISI and ICI effects. This allows detection of targets
directly using the time-domain observations in a low-complexity manner without
transformation to delay-Doppler domain and enables turning ISI/ICI into an
advantage for sensing.

  * •

Based on the new model, we propose an ISI/ICI embracing approach to surpass
the range/velocity ambiguity barrier encountered in existing OFDM [20, 21] and
OTFS [17, 18, 19] radar algorithms, which enables detection of any practically
relevant range/velocity.

  * •

We design a generalized likelihood ratio test (GLRT) based detector/estimator
for OTFS radar that simultaneously compensates for the ISI/ICI and harnesses
their properties to improve performance.

##  II OTFS Radar Signal Model

We consider an OTFS JRC system consisting of an OTFS JRC transceiver and an
OTFS communications receiver, as shown in Fig. 1. The OTFS transceiver is
equipped with an OTFS transmitter that generates the JRC signal for
communicating with the OTFS communication receiver and an OTFS radar receiver
that processes the backscattered signals for target detection. In this
section, we derive OTFS radar signal model for the single-CP OTFS modulation
architecture from [22, 18, 23, 19]. As our contribution is not on OTFS
communication, the processing at the OTFS communication receiver will be
ignored.

The OTFS system has a total bandwidth $N\Delta f$ and total frame duration
$MT$ (excluding any CP), where $N$ and $M$ denote the number of subcarriers
and the number of symbols, respectively, $\Delta f$ is the subcarrier spacing
and $T=1/\Delta f$ represents the symbol duration. Let
$\mathbf{X}^{\rm{DD}}\in\mathbb{C}^{N\times M}$ denote the two-dimensional
(2-D) OTFS frame in the delay-Doppler domain consisting of $NM$ transmit data
symbols that reside on the delay-Doppler grid (see Fig. 1)

$\mathcal{G}=\left\\{\left(\frac{n}{N\Delta
f},\frac{m}{MT}\right)\leavevmode\nobreak\ \Big{\lvert}\leavevmode\nobreak\
0\leq n\leq N-1,0\leq m\leq M-1\right\\}\leavevmode\nobreak\ .$

Applying an inverse symplectic finite Fourier transform (ISFFT) (i.e., an
$N$-point FFT over the columns and an $M$-point IFFT over the rows of
$\mathbf{X}^{\rm{DD}}$), we transform the 2-D transmit data block from the
delay-Doppler domain to the frequency-time domain [12, 24, 22, 25, 26]

$\mathbf{X}=\mathbf{F}_{N}\mathbf{X}^{\rm{DD}}\mathbf{F}_{M}^{H}\leavevmode\nobreak\
,$ (1)

where $\mathbf{X}\in\mathbb{C}^{N\times M}$ is the frequency-time domain
signal and $\mathbf{F}_{N}\in\mathbb{C}^{N\times N}$ is the unitary DFT matrix
with $\left[\mathbf{F}_{N}\right]_{\ell,n}=\frac{1}{\sqrt{N}}e^{-j2\pi
n\frac{\ell}{N}}$.

To map the frequency-time domain 2-D sequence $\mathbf{X}$ to a time domain
signal transmitted over the wireless channel, we apply the Heisenberg
transform [12, 26], which entails an $N$-point IFFT together with a transmit
pulse-shaping waveform $g_{\rm{tx}}(t)$ (which is time limited to
$\left[0,\,T\right]$). The time domain signal corresponding to the
${m^{\rm{th}}}$ symbol after the Heisenberg transform can be written as

$s_{m}(t)=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}x_{n,m}\,e^{j2\pi n\Delta
ft}g_{\rm{tx}}(t)\,,\leavevmode\nobreak\ 0\leq t\leq T\leavevmode\nobreak\ ,$
(2)

where $x_{n,m}\triangleq\left[\mathbf{X}\right]_{n,m}$. Hence, the time domain
signal for the entire OTFS frame without CP is given by

$s(t)=\sum_{m=0}^{M-1}s_{m}(t-mT)\,,\leavevmode\nobreak\ 0\leq t\leq
MT\leavevmode\nobreak\ .$ (3)

Finally, the entire time domain signal with CP is given by

$s_{\rm{CP}}(t)=\begin{cases}s(t),&\leavevmode\nobreak\ \leavevmode\nobreak\
0\leq t\leq MT\\\ s(t+MT),&\leavevmode\nobreak\ \leavevmode\nobreak\
-T_{\rm{cp}}\leq t\leq 0\end{cases}\leavevmode\nobreak\ ,$ (4)

where $T_{\rm{cp}}$ denotes the CP duration.

We consider a $K$-tap doubly selective, narrowband radar channel model as
$h(t,\tau)=\sum_{k=0}^{K-1}\alpha_{k}\delta(\tau-\tau_{k})e^{j2\pi\nu_{k}t}$
[27, 19], where the ${k^{\rm{th}}}$ target is characterized by a complex
channel gain $\alpha_{k}$, a round-trip delay $\tau_{k}=2R_{k}/c$ and a
Doppler shift $\nu_{k}=2v_{k}/\lambda$, with $R_{k}$, $v_{k}$, $c$ and
$\lambda$ denoting the range, radial velocity, speed of propagation and
carrier wavelength, respectively.

After transmission through the channel, the backscattered signal at the OTFS
radar receiver can be expressed as

$\displaystyle y_{\rm{CP}}(t)$
$\displaystyle=\sum_{k=0}^{K-1}\alpha_{k}s_{\rm{CP}}(t-\tau_{k})e^{j2\pi\nu_{k}t}+w(t)$
(5)

for $-T_{\rm{cp}}\leq t\leq MT$, where $s_{\rm{CP}}(t)$ is given by (4) and
$w(t)$ is additive white Gaussian noise. We assume that the CP duration is
larger than the round-trip delay of the furthermost target, i.e.,
$T_{\rm{cp}}\geq\max_{k\in\\{0,\ldots,K-1\\}}\tau_{k}$.

##  III Novel OTFS Signal Model for ISI/ICI Exploitation

In this section, we derive a novel compact representation of the received
signal in (5) for OTFS radar and show how it can be exploited to improve radar
performance.

###  III-A Derivation of Received Signal for OTFS Radar

We first remove the CP in (5) (i.e., the interval
$\left[-T_{\rm{cp}},\,0\right]$) to obtain

$y(t)=\begin{cases}y_{\rm{CP}}(t),&\leavevmode\nobreak\ \leavevmode\nobreak\
0\leq t\leq MT\\\ 0,&\leavevmode\nobreak\ \leavevmode\nobreak\
-T_{\rm{cp}}\leq t\leq 0\end{cases}\leavevmode\nobreak\ .$ (6)

Under the assumption $T_{\rm{cp}}\geq\max_{k}\tau_{k}$, the signal in (6)
becomes a superposition of cyclically shifted copies of the transmit signal
$s(t)$ in (3) [22, eq. (6)]:

$y(t)=\sum_{k=0}^{K-1}\alpha_{k}s([t-\tau_{k}]_{MT})e^{j2\pi\nu_{k}t}+w(t)\leavevmode\nobreak\
,$ (7)

where $[\cdot]_{T}$ denotes modulo-$T$.

Let $S(f)\triangleq\mathcal{F}\\{s(t)\\}=\int_{0}^{MT}s(t)e^{-j2\pi
ft}\;\mathrm{d}t$ denote the Fourier transform of $s(t)$. Then, a cyclic shift
of $s(t)$ corresponds to a phase shift in $S(f)$, i.e.,

$s([t-\tau_{k}]_{MT})=\mathcal{F}^{-1}\big{\\{}S(f)e^{-j2\pi
f\tau_{k}}\big{\\}}\leavevmode\nobreak\ ,$ (8)

where $\mathcal{F}^{-1}\\{\cdot\\}$ represents the inverse Fourier transform.
Sampling the time domain at $t=\ell T/N$ for $\ell=0,\ldots,NM-1$ and the
frequency domain at $f=n\Delta f/M$ for $n=0,\ldots,NM-1$, the equivalent
discrete-time representation of (7) can be written as

$\mathbf{y}=\sum_{k=0}^{K-1}\alpha_{k}\,\mathbf{F}^{H}\big{(}\mathbf{F}\mathbf{s}\odot\mathbf{b}(\tau_{k})\big{)}\odot\mathbf{c}(\nu_{k})+\mathbf{w}\leavevmode\nobreak\
,$ (9)

where $\odot$ is the Hadamard (element-wise) product,
$\mathbf{F}\in\mathbb{C}^{NM\times NM}$ is the unitary DFT matrix,
$\mathbf{s}\in\mathbb{C}^{NM\times 1}$ and $\mathbf{y}\in\mathbb{C}^{NM\times
1}$ denote the sampled versions of $s(t)$ and $y(t)$, respectively,
$\mathbf{w}\in\mathbb{C}^{NM\times 1}$ is the additive white Gaussian noise
with $\mathbf{w}\sim{\mathcal{CN}}({{\bm{0}}},\sigma^{2}{\bm{\mathrm{I}}})$,

$\mathbf{b}(\tau)=\mathbf{b}_{N}(\tau)\otimes\mathbf{b}_{M}(\tau)\in\mathbb{C}^{NM\times
1}$ (10)

is the frequency-domain steering vector with

$\displaystyle\mathbf{b}_{M}(\tau)$
$\displaystyle\triangleq{\left[1,e^{-j2\pi\frac{1}{M}\Delta
f\tau},\ldots,e^{-j2\pi\frac{M-1}{M}\Delta
f\tau}\right]}^{T}\in\mathbb{C}^{M\times 1}\leavevmode\nobreak\ ,$
$\displaystyle\mathbf{b}_{N}(\tau)$
$\displaystyle\triangleq{\left[1,e^{-j2\pi\Delta
f\tau},\ldots,e^{-j2\pi(N-1)\Delta f\tau}\right]}^{T}\in\mathbb{C}^{N\times
1}\leavevmode\nobreak\ ,$

and,

$\mathbf{c}(\nu)=\mathbf{c}_{M}(\nu)\otimes\mathbf{c}_{N}(\nu)\in\mathbb{C}^{NM\times
1}$ (11)

is the temporal steering vector with

$\displaystyle\mathbf{c}_{M}(\nu)$ $\displaystyle\triangleq{\left[1,e^{j2\pi
T\nu},\ldots,e^{j2\pi(M-1)T\nu}\right]}^{T}\in\mathbb{C}^{M\times
1}\leavevmode\nobreak\ ,$ $\displaystyle\mathbf{c}_{N}(\nu)$
$\displaystyle\triangleq{\left[1,e^{j2\pi\frac{T}{N}\nu},\ldots,e^{j2\pi\frac{T(N-1)}{N}\nu}\right]}^{T}\in\mathbb{C}^{N\times
1}\leavevmode\nobreak\ .$

Here, $\otimes$ denotes the Kronecker product. The vectors
$\mathbf{b}_{N}(\tau)$ and $\mathbf{c}_{M}(\nu)$ are commonly encountered in
OFDM radar, used for recovering the range and velocity, respectively. In
contrast, the vectors $\mathbf{b}_{M}(\tau)$ and $\mathbf{c}_{N}(\nu)$ are
commonly seen as disturbances that degrade the performance. In particular,
$\mathbf{b}_{N}(\tau)$ quantifies delay-dependent frequency-domain phase
rotations corresponding to the Fourier transform of fast-time (hereafter
called fast-frequency domain), while $\mathbf{b}_{M}(\tau)$ involves inter-
symbol (slow-time) delay-dependent phase rotations (hereafter called slow-
frequency domain), leading to inter-symbol interference (ISI). Similarly,
$\mathbf{c}_{N}(\nu)$ represents Doppler-induced fast-time phase rotations
causing inter-carrier interference (ICI), similar to the carrier frequency
offset (CFO) effect in OFDM communications [28], while $\mathbf{c}_{M}(\nu)$
captures Doppler-dependent slow-time phase progressions.

###  III-B Increasing Unambiguous Detection Intervals via ISI and ICI
Exploitation

The steering vector structures in (10) and (11) suggest that radar receivers
in OTFS suffer from both ISI and ICI effects [19], similar to the case of OTFS
communications [26]. However, unlike OTFS communications, these two effects
can be turned into an advantage for OTFS radar. In particular, ISI manifests
itself through the slow-frequency steering vector $\mathbf{b}_{M}(\tau)$ and
enables sampling the available bandwidth $N\Delta f$ at integer multiples of
$\Delta f/M$ as observed from the Kronecker structure in (10). In ISI-free
operation, the steering vector in (10) would only involve the fast-frequency
component $\mathbf{b}_{N}(\tau)$, which can sample the bandwidth with a
spacing of $\Delta f$, i.e., the subcarrier spacing. Hence, ISI can increase
the maximum detectable unambiguous delay by a factor of $M$ by sampling the
frequency domain $M$ times faster compared to a standard ISI-free radar
operation (e.g., in OFDM-based OTFS radar111Contrary to single-CP OTFS, OFDM-
based OTFS systems use separate CPs for each symbol in the OTFS/OFDM frame to
circumvent the ISI effect [24, 29].). More precisely, the unambiguous delays
with and without ISI are given, respectively, by

$\displaystyle\tau^{\rm{ISI}}_{\rm{max}}=\min\Bigg{\\{}\frac{M}{\Delta
f},T_{\rm{cp}}\Bigg{\\}},\leavevmode\nobreak\
\tau_{\rm{max}}=\min\Bigg{\\{}\frac{1}{\Delta
f},T_{\rm{cp}}\Bigg{\\}}\leavevmode\nobreak\ ,$ (12)

where $T_{\rm{cp}}$ is the upper limit to prevent inter-frame interference
(i.e., between consecutive OTFS frames) in OTFS radar and to prevent ISI
between consecutive symbols in OFDM-based OTFS radar.

Similarly to ISI, ICI can be exploited to increase the maximum detectable
unambiguous Doppler by a factor of $N$ [30]. In addition to the standard slow-
time steering vector $\mathbf{c}_{M}(\nu)$, Doppler-dependent phase rotations
can also be captured by the fast-time steering vector $\mathbf{c}_{N}(\nu)$,
which allows sampling the entire time window $MT$ with an interval of $T/N$.
Therefore, the unambiguous Doppler values with and without ICI can be
expressed, respectively, as

$\displaystyle\nu^{\rm{ICI}}_{\rm{max}}=\frac{N}{T},\leavevmode\nobreak\
\nu_{\rm{max}}=\frac{1}{T}\leavevmode\nobreak\ .$ (13)

##  IV Detection and Estimation with OTFS Radar

In this section, we provide the problem statement for OTFS radar sensing and
design a low-complexity detection and estimation scheme based on the signal
model in Sec. III-A.

###  IV-A Problem Formulation for OTFS Radar Sensing

Given the transmit signal $\mathbf{s}$, the problem of interest for OTFS radar
sensing is to detect the presence of multiple targets and estimate their
parameters, i.e., their gain-delay-Doppler triplets
$\\{(\alpha_{k},\tau_{k},\nu_{k})\\}_{k=0}^{K-1}$. This detection and
estimation is done from the time domain observations in (7), or, equivalently,
from its sampled version in (9). Unlike the existing works in the OTFS radar
literature, e.g., [17, 18, 19], where estimator design is based on the
received symbols in the delay-Doppler domain, we adopt a low-complexity
approach that performs detection/estimation directly using time-domain
observations without transforming them into delay-Doppler domain.

###  IV-B GLRT for Detection/Estimation in OTFS Radar

To design a detector, we first rewrite the signal model in (9) as

$\displaystyle\mathbf{y}=\sum_{k=0}^{K-1}\alpha_{k}\,\mathbf{C}(\nu_{k})\mathbf{F}^{H}\mathbf{B}(\tau_{k})\mathbf{F}\mathbf{s}+\mathbf{w}\leavevmode\nobreak\
,$ (14)

where
$\mathbf{B}(\tau)\triangleq{\rm{diag}}\left(\mathbf{b}(\tau)\right)\in\mathbb{C}^{NM\times
NM}$ and
$\mathbf{C}(\nu)\triangleq{\rm{diag}}\left(\mathbf{c}(\nu)\right)\in\mathbb{C}^{NM\times
NM}$. The hypothesis testing problem to test the presence of a single target
in (14) can be expressed as

$\displaystyle\mathbf{y}=\begin{cases}\mathbf{w},&\leavevmode\nobreak\
\leavevmode\nobreak\ {\rm{under\leavevmode\nobreak\ \mathcal{H}_{0}}}\\\
\alpha\,\mathbf{C}(\nu)\mathbf{F}^{H}\mathbf{B}(\tau)\mathbf{F}\mathbf{s}+\mathbf{w},&\leavevmode\nobreak\
\leavevmode\nobreak\ {\rm{under\leavevmode\nobreak\
\mathcal{H}_{1}}}\end{cases}\leavevmode\nobreak\ ,$ (15)

where the hypotheses $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$ represent the
absence and presence of a target, respectively. To solve (15), we treat
$\alpha$, $\tau$ and $\nu$ as deterministic unknown parameters and resort to
the GLRT

$\Lambda(\mathbf{y})=\frac{\max_{\alpha,\tau,\nu}p(\mathbf{y}\,\lvert\,\mathcal{H}_{1};\alpha,\tau,\nu)}{p(\mathbf{y}\,\lvert\,\mathcal{H}_{0})}\underset{\mathcal{H}_{0}}{\overset{\mathcal{H}_{1}}{\gtrless}}\eta\leavevmode\nobreak\
,$ (16)

where $\eta$ is the threshold. Under the assumption
$\mathbf{w}\sim{\mathcal{CN}}({{\bm{0}}},\sigma^{2}{\bm{\mathrm{I}}})$, the
GLRT becomes

$\displaystyle\Lambda(\mathbf{y})$ (17)
$\displaystyle=\frac{\exp\left(-\frac{1}{\sigma^{2}}\min_{\alpha,\tau,\nu}\left\lVert\mathbf{y}-\alpha\,\mathbf{C}(\nu)\mathbf{F}^{H}\mathbf{B}(\tau)\mathbf{F}\mathbf{s}\right\rVert^{2}\right)}{\exp\left(-\frac{1}{\sigma^{2}}\left\lVert\mathbf{y}\right\rVert^{2}\right)}\underset{\mathcal{H}_{0}}{\overset{\mathcal{H}_{1}}{\gtrless}}\eta\leavevmode\nobreak\
.$

For fixed $\tau$ and $\nu$, the optimal channel gain that maximizes the
numerator in (17) is given by

$\displaystyle\alpha=\frac{\mathbf{s}^{H}\mathbf{F}^{H}\mathbf{B}^{H}(\tau)\mathbf{F}\mathbf{C}^{H}(\nu)\mathbf{y}}{\left\lVert\mathbf{s}\right\rVert^{2}}\leavevmode\nobreak\
.$ (18)

Plugging (18) back into (17) and taking the logarithm, we have the detection
test

$\max_{\tau,\nu}\frac{\Big{\lvert}\mathbf{s}^{H}\mathbf{F}^{H}\mathbf{B}^{H}(\tau)\mathbf{F}\mathbf{C}^{H}(\nu)\mathbf{y}\Big{\rvert}^{2}}{\sigma^{2}\left\lVert\mathbf{s}\right\rVert^{2}}\underset{\mathcal{H}_{0}}{\overset{\mathcal{H}_{1}}{\gtrless}}\widetilde{\eta}\leavevmode\nobreak\
,$ (19)

where $\widetilde{\eta}=\log\eta$. To account for the presence of multiple
targets, the decision statistic in (19) can be computed over a discretized
delay-Doppler region and targets are declared at those locations where there
is a peak exceeding the threshold [31, Ch. 6.2.4].

##  V Numerical Results

In this section, we assess the performance of the proposed 2-D GLRT based
detector/estimator in (19). The OTFS radar observations are generated using
(5) instead of the proposed model in (9) to provide an implicit verification
of the transition from (5) to (9). As a benchmark, we consider a standard 2-D
FFT based processing traditionally employed in OFDM radar [20, 21]. To
evaluate detection performances for both GLRT and FFT methods, a cell-
averaging CFAR detector is employed with the probability of false alarm
$P_{\rm{fa}}=10^{-4}$ to declare targets in the delay-Doppler domain. For a
target with channel gain $\alpha$, we define the signal-to-noise ratio (SNR)
as ${\rm{SNR}}=\lvert\alpha\rvert^{2}/\sigma^{2}$. In addition, a rectangular
pulse shaping waveform is used for $g_{\rm{tx}}(t)$ in (2). We consider two
different parameter sets for OTFS, as shown in Table I, to illustrate the
results in both ISI-dominant (i.e., high $\Delta f$) and ICI-dominant (i.e.,
small $\Delta f$) operation regimes. In the ISI-dominant regime, maximum range
is the limiting factor for radar detection performance, while in the ICI-
dominant regime, radar performance is mainly limited by maximum velocity.

Table I: OTFS Parameter Sets for Simulations Parameter | ISI-dominant | ICI-dominant
---|---|---
| Regime | Regime
Carrier Frequency, $f_{c}$ | $60\,\rm{GHz}$ | $60\,\rm{GHz}$
Subcarrier Spacing, $\Delta f$ | $781.3\,\rm{kHz}$ | $48.8\,\rm{kHz}$
Number of Subcarriers, $N$ | $64$ | $1024$
Total Bandwidth, $B$ | $50\,\rm{MHz}$ | $50\,\rm{MHz}$
Symbol Duration, $T$ | $1.28\,\rm{\mu s}$ | $20.48\,\rm{\mu s}$
Cyclic Prefix Duration, $T_{\rm{cp}}$ | $7.68\,\rm{\mu s}$ | $20.48\,\rm{\mu s}$
Number of Symbols, $M$ | $64$ | $8$
Frame Duration, $MT+T_{\rm{cp}}$ | $89.6\,\rm{\mu s}$ | $184.3\,\rm{\mu s}$
Range Resolution, $\Delta R$ | $3\,\rm{m}$ | $3\,\rm{m}$
Maximum Range, $R_{\rm{max}}$ | $192\,\rm{m}$ | $3072\,\rm{m}$
(Standard) |  |
Maximum Range, $R^{\rm{ISI}}_{\rm{max}}$ | $1152\,\rm{m}$ | $3072\,\rm{m}$
(ISI Embracing) |  |
Velocity Resolution, $\Delta v$ | $30.5\,\rm{m/s}$ | $15.3\,\rm{m/s}$
Maximum Velocity, $v_{\rm{max}}$ | $\pm 976.6\,\rm{m/s}$ | $\pm 61\,\rm{m/s}$
(Standard) |  |
Maximum Velocity, $v^{\rm{ICI}}_{\rm{max}}$ | no practical limit | no practical limit
(ICI Embracing) |  |

###  V-A Example 1: ISI-dominant Regime

In the ISI-dominant regime, we consider a scenario with four targets located
at the same velocity ($20$ m/s), but with different ranges, as shown in Fig. 2
(left). As shown in Fig. 2 (right), the standard FFT processing [20, 21] or
the existing OTFS detectors [17, 18, 19] can detect at most two targets as
Target 3 and Target 4 fall into the same range-velocity bin as Target 1 and
Target 2, respectively.

Figure 2: Scenario for OTFS sensing in the ISI-dominant regime, where target
SNRs are given by $\\{25,10,20,10\\}\,\rm{dB}$, respectively.

Figure 3: ISI-dominant regime: Range profile at $v=20\,\rm{m/s}$ obtained by
the different methods 3 for the scenario in Fig. 2, and 3 for a modified
version of the scenario in Fig. 2, where Target 3 moved $30\,\rm{m}$ further
and Target 4 moved $30\,\rm{m}$ closer.

Fig. 3 shows the range profiles222The range profile of the 2-D GLRT in (19) is
obtained by plotting the decision statistic in (19) for a fixed Doppler $\nu$
over an interval of delay values $\tau$. For the 2-D FFT method, range profile
corresponds to the range slice of the 2-D FFT output, taken from a certain
Doppler $\nu$. of the considered methods in two different scenarios for a
single noise realization. It is observed that by virtue of the ISI embracing,
the proposed GLRT detector in (19) can detect four targets separately by
increasing the maximum range by a factor of $6$ (see (12) and Table I),
whereas the 2-D FFT yields a peak only at the location of Target 1. In
addition, the GLRT detector achieves lower side-lobe levels than the FFT
method by taking into account the ISI in detector design, which enables
compensating for its effect on the range profile. Moreover, even when Target 3
and Target 4 are displaced in Fig. 3 so that four targets are resolvable in
the ambiguity region, the FFT can only detect Target 1 due to the strong ISI
effect, while all the target peaks are clearly visible in the range profile of
the GLRT detector. Therefore, the proposed approach can simultaneously
mitigate ISI to have low side-lobes and embrace the information conveyed by
ISI to detect targets beyond the standard maximum range limit $R_{\rm{max}}$.

To illustrate the detection and estimation performance, we simulate $100$
independent Monte Carlo noise realizations and choose Target 2 as the
reference target. Fig. 4 shows the probability of detection and root mean-
squared error (RMSE) of range estimates of the reference target, and the
average number of false alarms (detections that do not belong to any of the
targets) with respect to the SNR of the reference target. It is seen that the
proposed GLRT detector/estimator significantly outperforms the standard FFT
method in terms of both detection and estimation performances. This is
accomplished through the ISI-aware modeling in (9) and the corresponding
detector design in (19), which performs ISI compensation via the term
$\mathbf{B}^{H}(\tau)$. We note that the GLRT detector/estimator in (19)
performs block-wise processing of the entire OTFS frame ($NM$ symbols), while
the FFT method applies separate $N$\- and $M$-point FFTs over frequency and
time domains, respectively, which provides computational simplicity, but leads
to poor radar performance.

Figure 4: Detection and estimation performance of the considered OTFS radar
algorithms with respect to SNR in the ISI-dominant regime. 4 Probability of
detection, 4 average number of false alarms, and 4 range RMSE.

###  V-B Example 2: ICI-dominant Regime

For the ICI-dominant regime, the scenario in Fig. 5 is considered, where four
targets are located at the same range bin, but with different velocities.
Since the maximum velocity is small, Target 3 and Target 4 lie in the same
range-velocity bin as Target 1 and Target 2, respectively. Hence, standard
OFDM [20, 21] and OTFS [17, 18, 19] radar algorithms cannot distinguish Target
1-3 and Target 2-4 as two separate echoes.

Figure 5: Scenario for OTFS sensing in the ICI-dominant regime, where target
SNRs are given by $\\{25,10,20,10\\}\,\rm{dB}$, respectively.

Figure 6: ICI-dominant regime: Velocity profile at $R=120\,\rm{m}$ obtained by
the different methods 6 for the scenario in Fig. 5, and 6 for a modified
version of the scenario in Fig. 5, where Target 3 and Target 4 move
$10\,\rm{m/s}$ faster.

We plot the velocity profiles obtained by the considered methods in Fig. 6 for
a single noise realization. As indicated in (13) and Table I, the proposed
approach increases the maximum velocity by a factor of $N=1024$, which allows
detection of targets beyond the standard velocity limit $v_{\rm{max}}$. Hence,
the proposed approach can resolve four targets and detect their true
(unambiguous) velocities, while the FFT method can only detect two targets.
This is observed also in Fig. 6, where four targets have different ambiguous
velocities; the targets cannot be resolved by the FFT method due to poor
velocity resolution. Through the ICI embracing capability of the proposed
detector, all targets can be resolved.

Figure 7: Detection and estimation performance of the considered OTFS radar
algorithms with respect to SNR in the ICI-dominant regime. 7 Probability of
detection, 7 average number of false alarms, and 7 velocity RMSE.

Fig. 7 shows the detection and estimation performance curves with respect to
the SNR of the reference target, Target 2, averaged over $100$ realizations.
Similar to the ISI-dominant case, the proposed approach achieves significant
performance gains over the conventional FFT method in terms of both the
probability of detection and velocity RMSE by explicitly accounting for the
ICI effect in detector/estimator design.

##  VI Concluding Remarks

In this paper, we have considered the problem of radar sensing with OTFS
waveform and derived a low-complexity GLRT based detector/estimator by taking
into account the ISI and ICI effects caused by the lack of guard intervals in
time and frequency, respectively. The proposed radar processing approach not
only compensates for the ISI and ICI effects, but also turns them into an
advantage by significantly increasing maximum detectable range and velocity.
Simulation results have verified the ISI and ICI embracing capability of the
proposed radar receiver and demonstrated improved detection and estimation
performance with respect to the standard FFT method employed in OFDM radar.

## Acknowledgments

This work is supported by the EU Marie Curie Actions Individual Fellowship
project OTFS-RADCOM, grant nr. 888913, Vinnova grant 2018-01929 and the
European Commission through the H2020 project Hexa-X (Grant Agreement no.
101015956). The work of A. Alvarado is supported by the Dutch Technology
Foundation TTW, which is part of the Netherlands Organisation for Scientific
Research (NWO), under the project Integrated Cooperative Automated Vehicles
(i-CAVE).

## References

  * [1] W. Saad, M. Bennis _et al._ , “A vision of 6G wireless systems: Applications, trends, technologies, and open research problems,” _IEEE network_ , vol. 34, no. 3, pp. 134–142, 2019.
  * [2] N. Rajatheva, I. Atzeni _et al._ , “White paper on broadband connectivity in 6G,” _arXiv preprint arXiv:2004.14247_ , 2020.
  * [3] P. Yang, Y. Xiao _et al._ , “6G wireless communications: Vision and potential techniques,” _IEEE Network_ , vol. 33, no. 4, pp. 70–75, 2019\.
  * [4] F. Liu, C. Masouros _et al._ , “Joint radar and communication design: Applications, state-of-the-art, and the road ahead,” _IEEE Transactions on Communications_ , vol. 68, no. 6, pp. 3834–3862, 2020.
  * [5] C. De Lima, D. Belot _et al._ , “Convergent communication, sensing and localization in 6G systems: An overview of technologies, opportunities and challenges,” _IEEE Access_ , pp. 1–1, 2021.
  * [6] H. Wymeersch, J. He _et al._ , “Radio localization and mapping with reconfigurable intelligent surfaces: Challenges, opportunities, and research directions,” _IEEE Vehicular Technology Magazine_ , vol. 15, no. 4, pp. 52–61, 2020.
  * [7] H. Elayan, O. Amin _et al._ , “Terahertz communication: The opportunities of wireless technology beyond 5G,” in _2018 International Conference on Advanced Communication Technologies and Networking (CommNet)_. IEEE, 2018, pp. 1–5.
  * [8] P. Banelli, S. Buzzi _et al._ , “Modulation formats and waveforms for 5G networks: Who will be the heir of OFDM?: An overview of alternative modulation schemes for improved spectral efficiency,” _IEEE Signal Processing Magazine_ , vol. 31, no. 6, pp. 80–93, 2014.
  * [9] M. Braun, “OFDM radar algorithms in mobile communication networks,” _Ph. D. dissertation_ , 2014.
  * [10] T. Wang, J. G. Proakis _et al._ , “Performance degradation of ofdm systems due to Doppler spreading,” _IEEE Transactions on Wireless Communications_ , vol. 5, no. 6, pp. 1422–1432, 2006.
  * [11] G. Surabhi, R. M. Augustine _et al._ , “Peak-to-average power ratio of OTFS modulation,” _IEEE Communications Letters_ , vol. 23, no. 6, pp. 999–1002, 2019.
  * [12] R. Hadani, S. Rakib _et al._ , “Orthogonal time frequency space modulation,” in _2017 IEEE Wireless Communications and Networking Conference (WCNC)_. IEEE, 2017, pp. 1–6.
  * [13] P. Raviteja, Y. Hong _et al._ , “Practical pulse-shaping waveforms for reduced-cyclic-prefix OTFS,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 1, pp. 957–961, 2018.
  * [14] K. Murali and A. Chockalingam, “On OTFS modulation for high-doppler fading channels,” in _2018 Information Theory and Applications Workshop (ITA)_. IEEE, 2018, pp. 1–10.
  * [15] F. Wiffen, L. Sayer _et al._ , “Comparison of OTFS and OFDM in ray launched sub-6 ghz and mmwave line-of-sight mobility channels,” in _2018 IEEE 29th Annual International Symposium on Personal, Indoor and Mobile Radio Communications (PIMRC)_. IEEE, 2018, pp. 73–79.
  * [16] L. Gaudio, M. Kobayashi _et al._ , “Performance analysis of joint radar and communication using OFDM and OTFS,” _arXiv e-prints_ , p. arXiv:1902.01184, Feb 2019.
  * [17] L. Gaudio, M. Kobayashi _et al._ , “Hybrid digital-analog beamforming and MIMO radar with OTFS modulation,” _arXiv preprint arXiv:2009.08785_ , 2020.
  * [18] P. Raviteja, K. T. Phan _et al._ , “Orthogonal time frequency space (OTFS) modulation based radar system,” in _2019 IEEE Radar Conference (RadarConf)_. IEEE, 2019, pp. 1–6.
  * [19] L. Gaudio, M. Kobayashi _et al._ , “On the effectiveness of OTFS for joint radar parameter estimation and communication,” _IEEE Transactions on Wireless Communications_ , 2020.
  * [20] C. Sturm and W. Wiesbeck, “Waveform design and signal processing aspects for fusion of wireless communications and radar sensing,” _Proceedings of the IEEE_ , vol. 99, no. 7, pp. 1236–1259, July 2011.
  * [21] M. Braun, “OFDM radar algorithms in mobile communication networks,” _Karlsruher Institutes für Technologie_ , 2014.
  * [22] P. Raviteja, Y. Hong _et al._ , “Practical pulse-shaping waveforms for reduced-cyclic-prefix OTFS,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 1, pp. 957–961, 2018.
  * [23] Y. Ge, Q. Deng _et al._ , “Receiver design for OTFS with fractionally spaced sampling approach,” _arXiv preprint arXiv:2009.00806_ , 2020.
  * [24] W. Shen, L. Dai _et al._ , “Channel estimation for orthogonal time frequency space (OTFS) massive MIMO,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 16, pp. 4204–4217, 2019.
  * [25] A. Farhang, A. RezazadehReyhani _et al._ , “Low complexity modem structure for OFDM-based orthogonal time frequency space modulation,” _IEEE Wireless Communications Letters_ , vol. 7, no. 3, pp. 344–347, 2017.
  * [26] P. Raviteja, K. T. Phan _et al._ , “Interference cancellation and iterative detection for orthogonal time frequency space modulation,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 10, pp. 6501–6515, 2018.
  * [27] P. Kumari, J. Choi _et al._ , “IEEE 802.11ad-based radar: An approach to joint vehicular communication-radar system,” _IEEE Transactions on Vehicular Technology_ , vol. 67, no. 4, pp. 3012–3027, April 2018\.
  * [28] T. Roman, S. Visuri _et al._ , “Blind frequency synchronization in OFDM via diagonality criterion,” _IEEE Transactions on Signal Processing_ , vol. 54, no. 8, pp. 3125–3135, 2006.
  * [29] A. Farhang, A. RezazadehReyhani _et al._ , “Low complexity modem structure for OFDM-based orthogonal time frequency space modulation,” _IEEE Wireless Communications Letters_ , vol. 7, no. 3, pp. 344–347, 2017.
  * [30] M. F. Keskin, H. Wymeersch _et al._ , “MIMO-OFDM joint radar-communications: Is ICI friend or foe?” _arXiv preprint arXiv:2103.15694_ , 2021.
  * [31] M. A. Richards, _Fundamentals of Radar Signal Processing_. Tata McGraw-Hill Education, 2005.