Current address: ]Department of Materials Science and Engineering, Stanford
University, Stanford, CA, USA, 94305

Current address: ]Department of Physics and Astronomy, University of
Tennessee, Knoxville, TN, 37996, USA

# Evidence of $\phi_{0}$-Josephson junction from skewed diffraction patterns
in Sn-InSb nanowires

B. Zhang  Department of Physics and Astronomy, University of Pittsburgh,
Pittsburgh, PA, 15260, USA  Z. Li  Department of Physics and Astronomy,
University of Pittsburgh, Pittsburgh, PA, 15260, USA  V. Aguilar  Department
of Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA, 15260, USA
P. Zhang  Department of Physics and Astronomy, University of Pittsburgh,
Pittsburgh, PA, 15260, USA  M. Pendharkar  [  Electrical and Computer
Engineering, University of California, Santa Barbara, CA, 93106, USA  C.
Dempsey  Electrical and Computer Engineering, University of California, Santa
Barbara, CA, 93106, USA  J.S. Lee  [  California NanoSystems Institute,
University of California Santa Barbara, Santa Barbara, CA, 93106, USA  S.D.
Harrington  Materials Department, University of California Santa Barbara,
Santa Barbara, CA, 93106, USA  S. Tan  Department of Electrical and Computer
Engineering, University of Pittsburgh, Pittsburgh, PA, 15260, USA  Petersen
Institute of NanoScience and Engineering, University of Pittsburgh,
Pittsburgh, PA, 15260, USA  J.S. Meyer  Univ. Grenoble Alpes, CEA, Grenoble
INP, IRIG, Pheliqs, 38000, Grenoble, France.  M. Houzet  Univ. Grenoble Alpes,
CEA, Grenoble INP, IRIG, Pheliqs, 38000, Grenoble, France.  C.J. Palmstrøm
Electrical and Computer Engineering, University of California, Santa Barbara,
CA, 93106, USA  S.M. Frolov<EMAIL_ADDRESS>Department of Physics and
Astronomy, University of Pittsburgh, Pittsburgh, PA, 15260, USA

###### Abstract

We study Josephson junctions based on InSb nanowires with Sn shells. We
observe skewed critical current diffraction patterns: the maxima in forward
and reverse current bias are at different magnetic flux, with a displacement
of 20-40 mT. The skew is greatest when the external field is nearly
perpendicular to the nanowire, in the substrate plane. This orientation
suggests that spin-orbit interaction plays a role. We develop a
phenomenological model and perform tight-binding calculations, both methods
reproducing the essential features of the experiment. The effect modeled is
the $\phi_{0}$-Josephson junction with higher-order Josephson harmonics. The
system is of interest for Majorana studies: the effects are either precursor
to or concomitant with topological superconductivity. Current-phase relations
that lack inversion symmetry can also be used to design quantum circuits with
engineered nonlinearity.

##  I Introduction

Context. Interest in superconductor-semiconductor hybrid structures is along
two directions. On the one hand, they are explored as materials for quantum
technologies, such as superconducting qubits [1, 2]. On the other hand, they
are a platform with high potential for the discovery of topological
superconductivity [3].

Background: Josephson $\varphi_{0}$-junction. In semiconductor nanowires, a
combination of induced superconductivity, spin-orbit interaction and spin
splitting can famously induce Majorana modes and topological superconductivity
[4, 5]. The same ingredients can induce an anomalous Josephson effect, known
as $\varphi_{0}$-junction [6]. The primary characteristic of a Josephson
junction is the current phase relation (CPR) [7]. The most common CPR is a
sinusoidal function $I(\phi)=I_{c}\sin(\phi)$, where $I(\phi)$ is the
Josephson supercurrent, $\phi$ is the phase difference between superconducting
leads, and $I_{c}$ is the critical current. In a $\varphi_{0}$-junction,
$I(\phi=0)\not=0$, which is equivalent to a phase offset $\varphi_{0}$ in a
sinusoidal CPR [8, 9, 10, 11, 12, 13, 6, 14, 15, 16, 17, 18]. The
$\varphi_{0}$-junction state can be accompanied by bias direction-dependent
critical current [6, 14, 15, 16], which was dubbed the “supercurrent diode
effect” [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31].

A related effect is the $\pi$-junction effect where the additional phase shift
is equal to $\pi$ [32]. Note that the $\pi$-junction is not a special case of
a $\varphi_{0}$-junction. The $\pi$-junction can be realized under more basic
conditions, for instance without any spin-orbit interaction [33].

Challenge. Attempts were made to detect the phase shift $\varphi_{0}$ in
superconducting quantum interference devices (SQUIDs) [34, 35]. However, large
magnetic fields were used. These fields were in the SQUID plane in order not
to induce flux in the loop. However, at high fields of hundreds of milliTesla
to a few Tesla, and given the SQUID area in the tens of square microns,
multiple flux quanta thread the SQUID due to fringing fields and imperfect
alignment even when the utmoust effort is applied to ensure strictly in-plane
applied fields. Furthermore, Josephson junctions based on semiconductor
nanowires are gate-tunable. This is typically thought of as an advantage due
to an extra control knob. But the electric field from the gate changes the
path of supercurrent in the nanowire, and with it the enclosed flux in the
SQUID. A change by a fraction of a flux quantum is plausible for a 100 nm
nanowire in a field of hundreds of mT. This shifts the SQUID interference
pattern due to extra flux and not due to the intrinsic spin-orbit interaction
and the associated $\varphi_{0}$-junction effect.

Figure 1:  A phenomenological model for the skewed diffraction pattern: (a)
CPR of a $\phi_{0}$-Josephson junctions ($\delta_{12}=0$) at an external field
$B=0.5B_{c}$ with $\phi_{0}=B$. Local maximum and minimum of CPR are taken as
critical current flow through positive and negative bias $I_{c+}$ and
$I_{c-}$. In this configuration we have $|I_{c+}|=|I_{c-}|$. (b) CPR of first
(red dotted), second harmonic (blue dot-dashed), and their sum (black solid)
at external field $B=0.5B_{c}$. The second harmonic has a shift of ground-
state phase $\delta\phi$ from the first order with a magnitude of
$\delta\phi=\delta_{12}/2$. (c) Skewed diffraction pattern generated with our
phenomenological model. (d) Coefficient $\gamma$ as a function of field
$B_{x}$ when $\delta_{12}=0$ and $\delta_{12}\neq 0$, respectively.

Approach. We use supercurrent diffraction patterns as means of investigating
the $\varphi_{0}$-junction state [19]. The diffraction pattern is the
evolution of the critical current in magnetic field. It can reveal exotic
effects such as d-wave superconductivity in corner junctions [36], the
presence of edge states in planar junctions [37], higher-order Josephson
harmonics [38]. The fields at which the effect manifests are as low as 10-20
mT, smaller than in previous works [34, 35] but large enough to dominate over
possible self-field effects.

Results list. In InSb-Sn nanowire junctions, we observe skewed diffraction
patterns. When the magnetic field is perpendicular to the nanowire and in-
plane (along the $\hat{x}$-direction), the switching current $I_{c}$ is
inversion symmetric with respect to flux and current bias. The pattern is
nearly-symmetric along the out-of-plane direction ($\hat{y}$) and along the
current flow direction ($\hat{z}$) [Fig. 2(a)]. The effect is observed over
wide ranges of gate voltage, which works against a fine-tuning explanation. At
the same time, no consistent effect of gate voltage on the skew magnitude is
observed. To interpret our results, we develop two models. The first is a
phenomenological model that illustrates how a two-component CPR with
$\varphi_{0}$ can result in a skewed diffraction pattern. The second is a
numerical tight-binding model which yields the $\varphi_{0}$-junction. Using
realistic junction parameters, the numerical model is capable of reproducing
the key experimental observations.

Brief Methods. Junctions are prepared by coating the standing InSb nanowire
with a 15 nm layer of Sn [39]. In front of the nanowire, another nanowire
shadows the flux of Sn to create two disconnected Sn segments. Shadow junction
wires are transferred onto chips patterned with local gates, contacts to wires
are made using standard electron beam lithography and thin film deposition.
Measurements are done in a dilution refrigerator with a base temperature of
$\sim$50 mK equipped with a 3D vector magnet. Experimental panels ploting
switching current are extracted from current-voltage measurements.

##  II Figure 1: Phenomenological Model

We present a minimal model capable of reproducing skewed diffraction patterns
due to the $\varphi_{0}$-junction effect. This is not a microscopic model, but
we also perform microscopic tight-binding calculations further below. In the
phenomenological model, we postulate a CPR with the first and second
sinusoidal harmonics. In addition to a variable phase across the junction
$\phi$, we allow for two phase offsets: the global parameter $\varphi_{0}$ and
the relative phase offset between the first and the second harmonics,
$\delta_{12}$:

$\displaystyle I(\phi)=$ $\displaystyle
I_{1}\sin(\phi+\varphi_{0})+I_{2}\sin(2\phi+2\varphi_{0}+\delta_{12})$ (1)

where $I_{1}$ and $I_{2}$ are the amplitude of each harmonic at zero external
field.

We first explain how this CPR realizes the so-called supercurrent diode
effect. If $I_{2}=0$, the CPR exhibits $I(\phi=0)\not=0$ with the trace
shifted by $\varphi_{0}$ [Fig. 1(a)]. This offset can, in principle, be
detected in a SQUID, but not in a single junction measurement. This is because
$I_{c+}=I_{c-}$ for the maximum supercurrent in the positive and negative bias
direction. The same is true if $I_{2}\not=0$, but $\delta_{12}=0$. However, if
a phase offset between the first and the second harmonics $\delta_{12}\not=0$,
we get $I_{c+}\not=I_{c-}$, and the current-voltage characteristic of the
junction becomes shifted upwards or downwards [Fig. 1(b)].

This is the supercurrent diode effect. In a single junction, the phase $\phi$
is free to adjust until the maximum supercurrent is reached, detected by a
switch into the finite voltage state. If the CPR has at least two components
with a phase offset between them, the switching current is different for
positive and negative bias directions. Note that the diode effect is not to be
confused with hysteretic supercurrent in underdamped junctions.

Next we generate skewed critical $I_{c}$ diffraction patterns. For this we
assume a phenomenological magnetic field dependence for model parameters:
$I_{1},I_{2}\propto(1-{B/B_{c}}^{2})$, where $B_{c}$ is the critical field,
reflecting suppression of critical current by magnetic field, and
$\varphi_{0},\delta_{12}\propto B$ to model the effect of spin-orbit coupling
[6]. Fig. 1(c) shows the maxima in $I_{c+}$ and $I_{c-}$ located symmetrically
around $B=0$, the diffraction pattern is skewed and inversion-symmetric.

To quantify the skew, we introduce a coefficient $\gamma$, which share the
same function as ’supercurrent diode coefficient/quality-factor’ used in Ref
[27, 24]:

$\displaystyle\gamma=\frac{\Delta
I_{c}}{(|I_{c}|)}=\frac{|I_{c+}|-|I_{c-}|}{(|I_{c+}|+|I_{c-}|)/2}$ (2)

Here $\Delta I_{c}$ is the difference between the magnitudes of $I_{c+}$ and
$I_{c-}$, and $|I_{c}|$ is the average critical current. In Fig. 1(d),
$\gamma=0$ for any field when $\delta_{12}=0$. When $\delta_{12}\not=0$,
$\gamma$ peaks at a finite field.

Experimentally $\gamma$ can also change sign as the field is increased without
going through zero. The phenomenological model provides a simple explanation
for this. Within the model, $\gamma=-2I_{2}/I_{1}\sin\delta_{12}$ when
$I_{2}\ll I_{1}$. Hence, $\gamma$ changes sign when $\delta_{12}=\pi$. There
is no special significance to this situation, and in particular it does not
correspond to any topological transition.

##  III Figure 2: Device Description

A schematic of the nanowire device is depicted in Fig. 2(a). InSb nanowire
(blue) covered by Sn shell (silver) is placed on top of local gate electrodes
and contacted by Ti/Au (gold) contacts. The direction of supercurrent is along
$\hat{z}$. The direction of spin-orbit field $B_{so}$ is induced by the
breaking of inversion symmetry in the device geometry, and indicated along
$-\hat{x}$ based on previous experiments [40].

The scanning electron microscope image of device A is in Fig. 2(b). The Sn
shell covers 3 of the 6 facets of the hexagonal nanowire cross-section. For
device A, the Sn shell faces the bottom, meaning it is oriented opposite to
the schematic in Fig. 2(a). For devices B and C the shell is likely on the
side (see supplementary information for cross-sectional STEM and AFM imaging).
In principle, the shell orientation can influence the direction of $B_{so}$,
but we see no evidence of this from the measurements.

Figure 2:  (a) Cartoon of a shadow nanowire Josephson Junction device.
Effective spin-orbit magnetic field ($B_{so}$) is indicated. (b) Scanning
Electron Microscope image of Device A.  Figure 3:  Skewed critical current
diffraction pattern in device A. (a) At gate voltage $V_{gate}=-2\mathrm{V}$.
The normal resistance is 4 k$\Omega$. (b) Line-cuts from (a) labeled with
color squares. Critical current at positive and negative bias are labeled as
$I_{c+}$ and $I_{c-}$. (c) Upper panel: switching currents extracted from
panel (a), Lower panel: Coefficient $\gamma$ calculated from panel (c).

##  IV Figure 3: Skewed Diffraction Pattern

Fig. 3 shows a representative skewed diffraction pattern from device A. The
field is applied along $\hat{x}$, in-plane and perpendicular to the nanowire.
In Fig. 3(a) the switching current is where the differential resistance
$dV/dI$ changes from zero (dark blue) to a finite value. The pattern is
visibly inversion-symmetric. In data processing, we treat current source that
gives voltage drop across the device smaller than $10\mu$V as superconducting
regime and vice versa. The switching currents $I_{c+}$ and $I_{c-}$ extracted
from panel (a) exhibit maxima displaced to positive and negative fields
respectively [Figs. 3(b),3(c)]. By taking I-V traces at $B_{x}=50$ mT and
$B_{x}=-50$ mT, we get $\Delta I_{c}=5-8$ nA [Fig. 3(b)], which is about 30%
of $|I_{c}|$. The maxima of $I_{c\pm}$ are at $B_{x}=\pm 20$ mT. The
coefficient $\gamma$ peaks at a higher field of order 100 mT right before the
collapse of $I_{c}$ at higher fields.

Supercurrent is not hysteretic in this regime. A hysteresis would manifest in
a vertical shift in the $\gamma$ dependence which would not be inversion-
symmetric. We discuss data acquisition and processing in the underdamped
regime in supplementary information. A 10 mT hysteresis in magnet field is
found in the measurements (Supplementary Section.V) and considered in our
discussion.

Figure 4:  Field rotation in Device B. Coefficient $\gamma$ as a function of
angle when the field is rotating in three orthogonal planes: x-y(a), x-z(b),
and y-z(c) with fixed strength $|B|=50~{}\text{mT}$. In (a) the external field
is along the $\hat{y}$-axis when $\theta_{xy}=90^{\circ}$ and $270^{\circ}$.
In (b) the external field is along $\hat{x}$-axis when $\theta_{z-x}$ =
$90^{\circ}$ and $270^{\circ}$.In (C) the external field is along
$\hat{y}$-axis when $\theta_{z-y}$ = $90^{\circ}$ and $270^{\circ}$.

##  V Figure 4: Magnetic Field Anisotropy

Device B has the same geometry as device A and its SEM picture can be found in
Fig. S1(c). We measure device B in the external field $|B|=50$ mT and rotate
the field in three orthogonal planes. Critical current are traced from zero
bias and extracted at where current bias gives differential resistance larger
than 2 k$\Omega$. Coefficients $\gamma$ are calculated based on extracted
$I_{c+}$ and $I_{c-}$ [Fig.4(a)-(c)]. $\gamma$ in x-y and x-z planes reach
zero when the external field is aligned to the $\hat{y}$ or $\hat{z}$ axis,
and is large along $\hat{x}$. In the y-z plane, $\gamma$ is significantly
reduced. More critical current diffraction patterns in fields at a different
angle in the x-y plane can be found in Fig. S4. The junction is made with
nanowires that are half-covered by superconductor shells. Shell orientation is
studied in the supplementary materials. Device A has its Sn shell on the
bottom and device B has the Sn shell on the side [Fig. S2]. Nevertheless, the
same general behavior, with skew coefficient $\gamma$ being largest along
$\hat{x}$ is observed in devices A, B, C (see Supplementary III).

Figure 5:  Numerical simulation results: (a) Critical currents ($I_{c+}$ and
$I_{c-}$) as a function of magnetic field $B_{x}$. The chemical potential is
set to two transverse or four spin-full modes ($\mu=8\mathrm{meV}$). (b)
Coefficient $\gamma$ as a function of $B_{x}$ corresponding to different
combinations of terms in the Hamiltonian (see legend). (c) Coefficient
$\gamma$ as a function of angle $\theta$ when the external field is rotating
in three orthogonal planes with fixed strength $|B|=50$ mT. The Zeeman effect
(g = 50) is present in all the results. Other parameters used in the
simulation are $\alpha=200$ $\mathrm{nm}\cdot\mathrm{meV}$,
$m_{eff}=0.015m_{e}$, temperature $T=100$ mK. The lattice constant $a=8$ nm,
the nanowire diameter $d_{1}=120$ nm, the outer diameter (with Sn shell)
$d_{2}=140$ nm and the coverage angle $\phi=180^{\circ}$. How chemical
potentials were chosen in the simulation is discussed in Supplementary Part V.

##  VI Figure 5: Tight-Binding Model

We numerically study the microscopic properties of the system within a tight-
binding model that has the same geometry as experiments (see supplementary
information for the description of the 3D model). This model was developed to
study supercurrent interference in nanowires using KWANT [41, 42, 43]. In that
project, the field orientation along the nanowire was primarily investigated.
Here, we rotate the field. We can toggle on and off spin-orbit interaction
($\alpha$), the orbital vector-potential effect of external field (A), while
Zeeman effect and disorder always remain on. For details of this model, see
Supplementary materials Section VII.

In Fig. 5(a) we reproduce a skewed diffraction pattern within the tight-
binding model for fields oriented along $\hat{x}$. Fig. 5(b), we illustrate
the role of spin-orbit interaction. The characteristic peak-antipeak structure
is present whenever $\alpha\neq 0$. The coefficient $\gamma$ remains zero when
only Zeeman effect of magnetic field is included.

Orbital effect only ($\alpha=0$, $\textbf{{A}}\not=0$) yields a similar
structure (see $B_{x}=80\mathrm{mT}$), limited in magnitude and field. We
believe this is a simulation artefact that appears when the external field is
perpendicular to the line connecting the center of the wire and the center of
the shell. (See Supplementary Section. VII for simulation at other chemical
potential) With all contributions toggled on, we reproduce the magnetic field
anisotropy of coefficient $\gamma$, compare Fig. 4 and Fig. 5(c).

The tight-binding model allows us to generate current-phase relations, which
exhibit $\varphi_{0}$ shifts as well as higher order harmonics and an extra
$\delta_{12}$ emerge and increase when field is along $\hat{x}$ [Fig. S16].
This agrees with the phenomenological model of Fig. 1. The parameters used in
the numerical model are discussed in the supplementary. Some of the
parameters, such as spin-orbit strength, exceed those previously reported for
InSb nanowires [40], however we cannot claim that matching this model to data
is a reliable way of extracting spin-orbit interaction strength.

Figure 6:  Comparison between the experiment (a) and the numerical simulation
(b). Coefficient $\gamma$ versus external field $B_{x}$ and gate voltage
(chemical potential $\mu$) are plotted as 2D maps to study skew shape as
function of gate voltage, $V_{gate}$ ($\mu$ in simulation). The experimental
data are taken from device B. Parameters used for the simulation are the same
as in Fig 5. Based on the discussion in Supplementary Section V, the gate
voltage range in our experiment is corresponding to chemical potential $\mu$ =
15-30 meV in the simulation. Another 2D map derived from Device C can be found
in Fig. S18

##  VII Figure 6: Gate Voltage Dependence

For device B we take gate sweeps of the supercurrent in a series of external
fields that are along the $\hat{x}$-axis [Fig.6(a)]. We observe that for
$V_{gate}>0$, $\gamma$ has a characteristic double-peak shape in magnetic
field which describes the skewed diffraction pattern. The behavior is observed
over a significant range of gate voltages, from near pinch-off to near
saturation of normal state conductance. The less regular traces at negative
gate voltages are too close to pinch-off where supercurrent is small. See
supplementary information for unprocessed gate voltage data for this and other
devices.

The data demonstrate that skewed diffraction patterns do not appear only at
fine-tuned values of gate voltage. On the other hand, there is no clear gate
voltage dependence of the magnitude of magnetic field extent of $\gamma$. The
behavior looks qualitatively similar in the tight-binding simulation [Fig.
6(b)]. Given the gate voltage range, we do not expect being able to
significantly tune the bulk Rashba spin-orbit coefficient in these InSb
nanowires.

##  VIII Alternative Explanation: Self-field effects

Skewed diffraction patterns were observed for decades in Josephson junctions
due to self-field effects, where current through the junction generates flux
in the junction. This effect is more pronounced in wide junctions with high
critical current density. We estimate the magnetic field generated by current
flow through a nanowire with the Biot-Savart Law. With current I = 300 nA,
nanowire radius r = 60 nm and Junction width L = 120 nm, we get self induced
field B $\sim$ 1 $\mu$T at the surface. In Fig.2 (a) and (c), the maximum and
minimum of $I_{c}$ are achieved at $|B|=20$ mT, which is thousands of times
larger than the self-induced field. Therefore, the self-field effect is not
enough to explain the skewed pattern we observe. Furthermore, the skew due to
self-fields should generally manifest for external fields along both $\hat{x}$
and $\hat{y}$, and be sensitive to the shell orientation, while the skew
reported here is strongest along $\hat{x}$ for two different shell
orientations.

##  IX Conclusion

Skewed diffraction patterns, observed in our experiment, are reproduced using
two theoretical models. The phenomenological model, and the tight-binding
model, agree that the imbalance between $I_{c+}$ and $I_{c-}$ and the critical
current maxima displaced for zero field are related to the current-phase
relation with two Josephson harmonics, and a phase-shift between them. The
tight-binding model yields phase-shifts $\varphi_{0}$ and $\delta_{12}$ due to
strong Rashba spin-orbit interaction. In turn, experiment shows that skew in
diffraction pattern is largest when the field is oriented along $\hat{x}$, the
likely orientation of spin-orbit effective field in nanowires. The effects are
observed in multiple nanowires and require no fine-tuning with gate voltages.

##  X Acknowledgements

We thank G. Badawy, S. Gazibegovic, E. Bakkers for providing InSb nanowires.
We thank Y. Nazarov, R. Mong for discussions. We acknowledge the use of shared
facilities of the NSF Materials Research Science and Engineering Center
(MRSEC) at the University of California Santa Barbara (DMR 1720256) and the
Nanotech UCSB Nanofabrication Facility.

##  XI Funding

Work supported by the NSF PIRE:HYBRID OISE-1743717, NSF Quantum Foundry funded
via the Q-AMASE-i program under award DMR-1906325, U.S. ONR and ARO and France
ANR through Grant No. ANR-17-PIRE-0001 (HYBRID).

##  XII Data Availability

Curated library of data extending beyond what is presented in the paper, as
well as simulation and data processing code are available at [44].

##  XIII Duration and Volume of Study

This project was started in May 2021, when the skewed diffraction pattern was
first observed in nanowires Josephson junction (devices were studied longer).
The experiments measurement ended in May 2022 and the simulation analysis
ended in August 2022. Devices studied in this project are made of nanowires
that is reported in ref [39]. 28 devices on 3 chips are measured during 5
cooldowns in dilution refrigerators, producing about 5200 datasets.

##  XIV References

## References

  * Larsen  _et al._ [2015] T. W. Larsen, K. D. Petersson, F. Kuemmeth, T. S. Jespersen, P. Krogstrup, J. Nygård, and C. M. Marcus, Physical review letters 115, 127001 (2015).
  * De Lange  _et al._ [2015] G. De Lange, B. Van Heck, A. Bruno, D. Van Woerkom, A. Geresdi, S. Plissard, E. Bakkers, A. Akhmerov, and L. DiCarlo, Physical review letters 115, 127002 (2015).
  * Frolov  _et al._ [2020] S. Frolov, M. Manfra, and J. Sau, Nature Physics 16, 718 (2020).
  * Lutchyn  _et al._ [2010] R. M. Lutchyn, J. D. Sau, and S. D. Sarma, Physical review letters 105, 077001 (2010).
  * Oreg  _et al._ [2010] Y. Oreg, G. Refael, and F. Von Oppen, Physical review letters 105, 177002 (2010).
  * Yokoyama  _et al._ [2014] T. Yokoyama, M. Eto, and Y. V.Nazarov, Physical Review B 89, 195407 (2014).
  * Golubov  _et al._ [2004] A. A. Golubov, M. Y. Kupriyanov, and E. Il’Ichev, Reviews of modern physics 76, 411 (2004).
  * Buzdin [2008] A. Buzdin, Physical review letters 101, 107005 (2008).
  * Reynoso  _et al._ [2008] A. Reynoso, G. Usaj, C. Balseiro, D. Feinberg, and M. Avignon, Physical review letters 101, 107001 (2008).
  * Zazunov  _et al._ [2009] A. Zazunov, R. Egger, T. Jonckheere, and T. Martin, Physical review letters 103, 147004 (2009).
  * Margaris  _et al._ [2010] I. Margaris, V. Paltoglou, and N. Flytzanis, Journal of Physics: Condensed Matter 22, 445701 (2010).
  * Reynoso  _et al._ [2012] A. A. Reynoso, G. Usaj, C. B. an D. Feinberg, and M. Avignon, Physical Review B 86, 214519 (2012).
  * Brunetti  _et al._ [2013] A. Brunetti, A. Zazunov, A. Kundu, and R. Egger, Physical Review B 88, 144515 (2013).
  * Silaev  _et al._ [2014] M. Silaev, A. Y. Aladyshkin, M. Silaeva, and A. Aladyshkina, Journal of Physics: Condensed Matter 26, 095702 (2014).
  * Chen  _et al._ [2018] C.-Z. Chen, J. J. He, M. N. Ali, G.-H. Lee, K. C. Fong, and K. T. Law, Physical Review B 98, 075430 (2018).
  * Minutillo  _et al._ [2018] M. Minutillo, D. Giuliano, P. Lucignano, A. Tagliacozzo, and G. Campagnano, Physical Review B 98, 144510 (2018).
  * Campagnano  _et al._ [2015] G. Campagnano, P. Lucignano, D. Giuliano, and A. Tagliacozzo, Journal of Physics: Condensed Matter 27, 205301 (2015).
  * Dolcini  _et al._ [2015] F. Dolcini, M. Houzet, and J. S. Meyer, Physical Review B 92, 035428 (2015).
  * Sickinger  _et al._ [2012] H. Sickinger, A. Lipman, M. Weides, R. Mints, H. Kohlstedt, D. Koelle, R. Kleiner, and E. Goldobin, Physical review letters 109, 107002 (2012).
  * Ando  _et al._ [2020] F. Ando, Y. Miyasaka, T. Li, J. Ishizuka, T. Arakawa, Y. Shiota, T. Moriyama, Y. Yanase, and T. Ono, Nature 584, 373 (2020).
  * Lyu  _et al._ [2021] Y.-Y. Lyu, J. Jiang, Y.-L. Wang, Z.-L. Xiao, S. Dong, Q.-H. Chen, M. V. Milošević, H. Wang, R. Divan, J. E. Pearson, _et al._ , Nature communications 12, 1 (2021).
  * Baumgartner  _et al._ [2022] C. Baumgartner, L. Fuchs, A. Costa, S. Reinhardt, S. Gronin, G. C. Gardner, T. Lindemann, M. J. Manfra, P. E. Faria Junior, D. Kochan, _et al._ , Nature nanotechnology 17, 39 (2022).
  * Wu  _et al._ [2022] H. Wu, Y. Wang, Y. Xu, P. K. Sivakumar, C. Pasco, U. Filippozzi, S. S. Parkin, Y.-J. Zeng, T. McQueen, and M. N. Ali, Nature 604, 653 (2022).
  * He  _et al._ [2022] J. J. He, Y. Tanaka, and N. Nagaosa, New Journal of Physics 24, 053014 (2022).
  * Kopasov  _et al._ [2021] A. Kopasov, A. Kutlin, and A. Mel’nikov, Physical Review B 103, 144520 (2021).
  * Daido  _et al._ [2022] A. Daido, Y. Ikeda, and Y. Yanase, Physical Review Letters 128, 037001 (2022).
  * Yuan and Fu [2022] N. F. Yuan and L. Fu, Proceedings of the National Academy of Sciences 119, e2119548119 (2022).
  * Hou  _et al._ [2022] Y. Hou, F. Nichele, H. Chi, A. Lodesani, Y. Wu, M. F. Ritter, D. Z. Haxell, M. Davydova, S. Ilić, F. S. Bergeret, _et al._ , arXiv preprint arXiv:2205.09276 (2022).
  * Gupta  _et al._ [2022] M. Gupta, G. V. Graziano, M. Pendharkar, J. T. Dong, C. P. Dempsey, C. Palmstrøm, and V. S. Pribiag, arXiv preprint arXiv:2206.08471 (2022).
  * Bauriedl  _et al._ [2022] L. Bauriedl, C. Bäuml, L. Fuchs, C. Baumgartner, N. Paulik, J. M. Bauer, K.-Q. Lin, J. M. Lupton, T. Taniguchi, K. Watanabe, _et al._ , Nature communications 13, 1 (2022).
  * Pal  _et al._ [2022] B. Pal, A. Chakraborty, P. K. Sivakumar, M. Davydova, A. K. Gopi, A. K. Pandeya, J. A. Krieger, Y. Zhang, M. Date, S. Ju, _et al._ , Nature physics 18, 1228 (2022).
  * Ryazanov  _et al._ [2001] V. Ryazanov, V. Oboznov, A. Y. Rusanov, A. Veretennikov, A. A. Golubov, and J. Aarts, Physical review letters 86, 2427 (2001).
  * Heikkilä  _et al._ [2000] T. T. Heikkilä, F. K. Wilhelm, and G. Schön, EPL (Europhysics Letters) 51, 434 (2000).
  * Szombati  _et al._ [2016] D. Szombati, S. Nadj-Perge, D. Car, S. Plissard, E. Bakkers, and L. Kouwenhoven, Nature Physics 12, 568 (2016).
  * Pita-Vidal  _et al._ [2020] M. Pita-Vidal, A. Bargerbos, C.-K. Yang, D. J. Van Woerkom, W. Pfaff, N. Haider, P. Krogstrup, L. P. Kouwenhoven, G. De Lange, and A. Kou, Physical Review Applied 14, 064038 (2020).
  * Wollman  _et al._ [1995] D. Wollman, D. J. Van Harlingen, J. Giapintzakis, and D. Ginsberg, Physical review letters 74, 797 (1995).
  * Hart  _et al._ [2014] S. Hart, H. Ren, T. Wagner, P. Leubner, M. Mühlbauer, C. Brüne, H. Buhmann, L. W. Molenkamp, and A. Yacoby, Nature Physics 10, 638 (2014).
  * Stoutimore  _et al._ [2018] M. Stoutimore, A. Rossolenko, V. Bolginov, V. Oboznov, A. Rusanov, D. Baranov, N. Pugach, S. Frolov, V. Ryazanov, and D. Van Harlingen, Physical review letters 121, 177702 (2018).
  * Pendharkar  _et al._ [2021] M. Pendharkar, B. Zhang, H. Wu, A. Zarassi, P. Zhang, C. Dempsey, J. Lee, S. Harrington, G. Badawy, S. Gazibegovic, _et al._ , Science 372, 508 (2021).
  * Nadj-Perge  _et al._ [2012] S. Nadj-Perge, V. Pribiag, J. Van den Berg, K. Zuo, S. Plissard, E. Bakkers, S. Frolov, and L. Kouwenhoven, Physical review letters 108, 166801 (2012).
  * Groth  _et al._ [2014] C. W. Groth, M. Wimmer, A. R. Akhmerov, and X. Waintal, New Journal of Physics 16, 063065 (2014).
  * Nijholt and Akhmerov [2016] B. Nijholt and A. R. Akhmerov, Physical Review B 93, 235434 (2016).
  * Zuo  _et al._ [2017] K. Zuo, V. Mourik, D. B. Szombati, B. Nijholt, D. J. Van Woerkom, A. Geresdi, J. Chen, V. P. Ostroukh, A. R. Akhmerov, S. R. Plissard, _et al._ , Physical review letters 119, 187704 (2017).
  * [44] DOI: 10.5281/zenodo.7374094 .

Supplementary Materials

##  XV Fabrication and Measurements

Nanowires are placed onto a Si chip that has predefined local gates.
Electrostatic local gates are patterned by 100 keV Electron Beam Lithography
(EBL) on undoped Si substrates. Local gates have mixed widths of 80 and 200nm
and are separated with a distance of 40 nm. Electron beam evaporation of
1.5/6nm Ti/PdAu is used to metalize the gates which are covered by 10nm of ALD
HfOx that serves as a dielectric layer [Fig. S1(a)].

Figure S1:  Fabrication steps in making nanowires Josephson junction (a) Local
gate chips made on a Si substrate (purple). Ti/AuPd gates (gold) are covered
by HfOx dielectric (Green) (b) InSb shadow nanowires that are half-covered by
Sn shells are transferred onto the chips. (c) Metal leads made with Ti/Au are
connected to the shell to form a Josephson junction.

Based on the Scanning electron microscope images of all devices that are
studied in this report, the length of the JJs made with Sn-InSb nanowires is
120-150nm. The width of the junction is the same as the width of the
nanowires, which is $\approx$ 120nm. After placing nanowires onto gates [Fig.
S1(b)], we cover the whole chip with PMMA 950 A4 electron beam resist. Resist
is dried at room temperature by pumping with a mechanical pump in a metal
desiccator for 24 hours. Then we use EBL to define normal lead patterns. After
development, we clean the residue of resist in an oxygen plasma asher. In the
electron beam evaporator, we first use an in-situ ion mill to remove AlOx
capping layer from the nanowires in the contact area, after which we deposit
10nm/130nm Ti/Au on the chips [Fig. S1(c)].

Transport 2-point measurements with currents source and voltage measurement
model in parallel are used, with several stages of filtering placed at
different temperatures.

##  XVI Device list and shell orientation

3 chips are studied in this project. Each chip contains multiple devices.
Among these chips, 5 devices demonstrate sharp critical currents that are
tuned by the gates. Skewed diffraction patterns taken from these devices are
all plotted and posted in this report. In the main text, we present skewed
diffraction patterns from device A [Fig. 3]. In device B, we study field
direction dependence by rotating fields in different planes and study gate
voltage and field effect on the skewed diffraction pattern with a 2D map [Fig.
4 and 6].

Figure S2:  Scanning electron microscope (SEM) image of Devices measured in
report. (a) Right: SEM image of Device A (Chip: QPC2). Left: Scanning
transmission electron microscopy (STEM) based Energy-dispersive x-ray
spectroscopy (EDS) spectroscopic elemental image of Device A. (b) Right: SEM
image of Device A1 (Chip: QPC2). Left: STEM based EDS spectroscopic elemental
image of Device A1. (c) Right: SEM image of Device B (Chip: QPC4). Left:
Atomic force microscopy (AFM) image of Device B. (d) SEM image of Device
C(Chip: QPC3) (e) SEM image of Device C1 (Chip: QPC3)

Devices A and A1 are on Chip QPC2. The EDX results show InSb nanowires and Sn
shells. No ferromagnetic materials were found in the junction. Sn shells are
half-covering the nanowires from the bottom side and are in touch with the
HfOx dielectric layer.

Device B is on Chip QPC4. The shell orientation of the device is studied with
Atomic Force Microscopy [Fig. S2(c)]. Based on the AFM images, we conclude
that the Sn shell is on the side.

Devices C and C1 are on Chip QPC3. The shell orientations are not studied.

##  XVII Diffraction pattern at different gates in field along three axes

###  XVII.1 Device A

Figure S3:  Diffraction patterns at different gates in the field along three
axes, measured in Device A (QPC2). (a) dV/dI differential resistance as a
function of $\hat{x}$-direction field $B_{x}$ and current bias. $\gamma$ is
calculated with the extracted magnitude $\Delta I_{c}$. (b) Diffraction
pattern when the field is applied parallel to the nanowires, along
$\hat{z}$-direction. (c) Diffraction pattern when the field is applied out of
substrate plane, along $\hat{y}$-direction.

###  XVII.2 Device A1

Figure S4:  Diffraction patterns at different gates in the field along
$\hat{x}$ and $\hat{z}$ axes, measured in Device A1 (QPC2). (a) dV/dI
differential resistance as function of $\hat{x}$-direction field $B_{x}$ and
current bias. $\gamma$ is calculated with the extracted magnitude $\Delta
I_{c}$. (b) Diffraction pattern when the field is applied parallel to the
nanowires, along $\hat{z}$-direction.

###  XVII.3 Device B

Figure S5:  Diffraction patterns at different gates in the field along
$\hat{x}$ and $\hat{z}$ axes, measured in Device B (QPC4). (a) dV/dI
differential resistance as function of $\hat{x}$-direction field $B_{x}$ and
current bias. $\gamma$ is calculated with the extracted magnitude $\Delta
I_{c}$. (b) Diffraction pattern when the field is applied parallel to the
nanowires, along $\hat{z}$-direction.

###  XVII.4 Device C

Figure S6:  Diffraction patterns at different gates in the field along three
axes, measured in Device C (QPC3)(a) dV/dI differential resistance as function
of $\hat{x}$-direction field $B_{x}$ and current bias. $\gamma$ is calculated
with the extracted magnitude $\Delta I_{c}$. (b) Diffraction pattern when
field is applied parallel to the nanowires, along $\hat{z}$-direction. (c)
Diffraction pattern when the field is applied out of substrate plane, along
$\hat{y}$-direction.

###  XVII.5 Device C1

Figure S7:  Diffraction patterns at different gates in the field along three
axes, measured in Device C1 (QPC3) (a) dV/dI differential resistance as
function of $\hat{x}$-direction field $B_{x}$ and current bias. $\gamma$ is
calculated with the extracted magnitude $\Delta I_{c}$. There is a shift of
gate, so strength of $I_{c}$ may be different from scans with other two fields
directions. (b) Diffraction pattern when field is applied parallel to the
nanowires, along $\hat{z}$-direction. (c) Diffraction pattern when the field
is applied out of substrate plane, along $\hat{y}$-direction.

##  XVIII Diffraction pattern at x-y plane, field is applied along angle
$\theta$ between $0^{\circ}$ to $180^{\circ}$, Device C

Figure S8:  Diffraction pattern at the x-y plane, the field is applied along
angle $\theta_{x-y}$ between $0^{\circ}$ to $180^{\circ}$, measured in Device
C. Critical current difference $\Delta I_{c}$ is extracted from measurement
data using a peak finder Python script.

##  XIX Characteristics of the Junctions and hysteresis in the measurement
setup

Figure S9:  Gate dependence taken at external field $|B|$ = 0, measurement
data are taken from Device A (a-c), Device B (d-f) and Device C (g-i). (a,d,g)
Current through the JJs as function of gate voltage $V_{gate}$ when bias
voltage is set to be $V_{bias}$ = 10mV. (b,e,h) V-I characteristics taken at
different gate voltages. (c,f,i) Upper panel: dV/dI differential resistance as
a function of current source and $V_{gate}$. Lower panel: extracted critical
current $I_{c}$ (blue) and $I_{c}R_{N}$ product (red) as functions of
$V_{gate}$.(j) Electrical conductance of the junction in the unit of quantum
conductance ($2e^{2}/h$) as a function of chemical potential $\mu$ when
simulating with parameters used in other figures. Followed by normal state
conductance as a function of gate voltage $V_{gate}$ measured in device A, B,
and C, respectively. Differential resistance are extracted from (c,f,i) and
converted to conductance for plotting. (k) Critical current difference
calculated with $\Delta I_{c}$ as functions of gate voltage at zero field,
measured in device A, B, and C. $I_{c+}$ and $I_{c-}$ are extracted from
(c,f,i)

The gate effect is studied by applying a bias voltage across the device
$V_{bias}=10mV$. Conduction channels in Device A [Fig. S9(a)] can be fully
closed by the tunnel gate. While Devices B and C [Fig. S9(d),(g)] cannot be
fully closed. The Josephson effect at zero magnetic field is best studied in
the current-bias configuration [Fig. S9(c),(f),(i)]. The switching current
from superconducting to normal state regime is demonstrated with a read peak
in differential resistance (referred to in other Figures as $I_{c}$). The
magnitude of $I_{c}$ is calculated with $|I_{c}|=(|I_{c+}|+|I_{c-}|)/2$, where
$I_{c+}$ and $I_{c-}$ are extracted with a peak finder Python script by
finding the two of largest differential resistance at each gate voltage.
$I_{c}$ increases at more positive gate voltage, while the extracted products
$I_{c}R_{N}$ ($R_{N}$ is the normal state resistance) are in the range of
200-300 $\mu\cdot eV$.

###  XIX.1 Chemical potential used in simulation and corresponding gate
voltage in experimental measurements

In an attempt to establish additional correspondence between experiment and
theory, we study conductance as a function of chemical potential $\mu$ in the
simulation [Fig. S9(j)]. The normal state resistance read from skewed pattern
in Fig. 2(a) in main text is 4 kOhm, which is close to two transverse modes or
four spin-full modes. So we choose $\mu$ = 8 meV to study the skew shape in
Fig. 4. When studying the direction-dependent supercurrent transport as
function of field direction, we get the normal state resistance about 1.5-2
kOhm. Which is corresponding to six to eight transverse modes or twelve to
eighteen spin-full mode, and chemical potential $\mu\approx 20meV$ in the
simulation.

The skew map presented in Fig. 5 in main text is measured from Device B. By
comparing normal state conductance $G$ as a function of $V_{gate}$ in
simulation and measurement data [Fig. S9(j)], we conclude the range used in
Fig. 5, which is -0.5V to 2V, is corresponding to chemical potential
$\mu=15-60meV$. Another skew map measured with Device C in Fig. S18 (c) has
gate voltage range from -3V to 2V, which is corresponding to $\mu=30-60meV$.

###  XIX.2 Hysteresis in the Josephson junction

Hysteresis in the current-voltage characteristics is studied by extracting
$\Delta I_{c}=|I_{c+}|-|I_{c-}|$ from devices A, B, and C and plotting them as
a function of gate voltage at zero magnetic fields [Fig. S9(k)]. We find the
magnitude of $I_{c-}$ is larger than $I_{c+}$ at more positive gate voltage in
Devices A and C. When extracting $\Delta I_{c}$ from skewed diffraction
patterns [Figs. S3(b), S6(c)], it is also easy to see the critical current is
larger at negative bias, which results in the gamma as a function of the
$B_{x}$ field that is no longer inversion symmetric about zero fields, zero
current. While in Device B, we didn’t see the obvious difference between
$I_{c+}$ and $I_{c-}$.

There are two methods used in current configuration measurements: 1\.
Unidirectional current sweeps, either from positive to negative, or from
negative to positive bias. This method is used in Devices A and C. Which
results in the hysteresis in current bias.

2\. Sweep from zero bias. At a fixed gate voltage or field, scan from zero
current bias to the positive and negative side respectively (0 to $I_{+}$ and
0 to $I_{-}$). So one data set only records scans with $I$ larger than zero
and vice versa. Then we combine two datasets to get a full scan. By doing
this, we can get rid of the hysteresis because only switching currents are
measured. This method is used in device B and the results are plotted in [Fig
.4, 6].

In the main text, the skewed diffraction pattern from Device A are studied at
a smaller gate voltage, $V_{gate}=-2V$ [Fig.3]. At this gate voltage,
hysteresis in the junction is not observed.

###  XIX.3 Hysteresis in the superconducting magnet

Figure S10:  (a) Skewed critical current diffraction patterns, field is
scanned from positive to negative and from negative to positive direction in
adjacent panels. Device A is used in this scan and gate voltage is set to be
$V_{gate}$ = -1V for this scan. (b) Extracted critical current difference
$\Delta I_{c}$ is plotted as a function of $\hat{x}$-field from two scans in
(a).

##  XX Model used in simulation

To simulate the superconductor-nanowire-superconductor Josephson junction in
the presence of external magnetic field, we consider the following Hamiltonian
for a nanowire that is covered by superconductor lead at both ends.

$\displaystyle H=\left(\frac{\mathbf{p}^{2}}{2m^{*}}-\mu+\delta
U\right)\tau_{z}+\alpha(p_{z}\sigma_{x}-p_{x}\sigma_{z})\tau_{z}+g\mu_{B}\mathbf{B}\cdot\hat{\sigma}+\Delta\tau_{x},$
(S1)

where $\tau_{i}$ and $\sigma_{i}$ are Pauli matrices act on particle-hole and
spin space respectively. $\mathbf{p}=-i\hbar\nabla+e\mathbf{A}\tau_{z}$ is the
canonical momentum, and the magnetic potential $\mathbf{A}$ is chosen to be
$[0,B_{z}x-B_{x}z,-B_{y}x]$, so that it is invariant along the $x-$direction.
Further, $m^{*}$ is the effective mass, $\mu$ is the chemical potential and
$\delta U$ represent the onsite disorder inside the nanowires. The Zeeman
effect is given by $g\mu_{B}\mathbf{B}\cdot\hat{\sigma}$ and the Rashba spin-
orbit coupling is given by $\alpha(p_{z}\sigma_{x}-p_{x}\sigma_{z})$. Finally,
$\Delta$ is the superconducting pairing potential.

We first construct a tight-binding model based on the Hamiltonian (S1), then
the critical current under different parameter configurations can be obtained
from the imaginary part of the Green’s function. These Green’s functions can
be obtained by using the KWANT package. We explain the code in detail below.

The first step is to construct a system with the scattering region and leads.
Here we use the function kwant.continuum.discretize to convert the 3D
translational symmetric Hamiltonian (S2) into a tight-binding system (Figure
S11).

$\displaystyle
H=\left(\frac{\hbar^{2}(p_{x}^{2}+p_{y}^{2}+p_{z}^{2})}{2m^{*}}-\mu+\delta
U\right)\tau_{z}+\alpha(p_{z}\sigma_{x}-p_{x}\sigma_{z})\tau_{z}+g\mu_{B}\mathbf{B}\cdot\hat{\sigma}+\Delta\tau_{x},$
(S2)

Here the Hamiltonian does not contain the orbital effect because
kwant.continuum.discretize cannot handle the systems with lower symmetry. To
include the orbital effect, we need to apply the Peierls substitution to the
hopping term. The hopping between two sites $\vec{x}$ and $\vec{x}_{0}$
becomes $t\to te^{i\phi}$, where
$\phi=-e\mathbf{A}\cdot(\vec{x}-\vec{x}_{0})/\hbar.$

Figure S11: Tight-binding model generated by the function
kwant.continuum.discretize. Red dots represent the infinite leads.

In order to calculate the critical current, besides normal leads and
superconducting leads (red region in Fig. S11), we need to add a virtual self-
energy lead to this system. Here we attach the lead in the middle of the
nanowire (yellow region in Fig. S12). Notice this self-energy lead is not
connected to external devices, and is only used to calculate the Green’s
function. Usually, the Green’s function of an infinite system contains
infinite entries. But now we can divide this nanowire into two parts: self
energy lead (L) and the other region (R). Then, we have

$\displaystyle\begin{pmatrix}G_{L}^{r}&G_{LR}^{r}\\\
G_{RL}^{r}&G_{R}^{r}\end{pmatrix}=\begin{pmatrix}E+i\eta-H_{Lead}&H_{C}\\\
H_{C}^{\dagger}&E+i\eta-H_{R}\end{pmatrix}^{-1},$ (S3)

where $H_{C}$ and $H_{C}^{\dagger}$ are the hopping between the lead and the
rest of the nanowire. By solving this equation we get

$\displaystyle
G_{L}^{r}=\left(E-H_{Lead}-H_{C}^{\dagger}(E-H_{R})^{-1}H_{C}^{\dagger}\right)^{-1}.$
(S4)

Thus the Green’s function of the finite self energy lead contains the
information about the whole system.

In the KWANT package, the retarded Green’s function of the self-energy lead
can be obtained by using the function kwant.solvers.greens_function. We first
calculate the Green’s function $G^{r}_{L}(0)$ without the phase difference
between the two superconducting leads. Then the Green’s function with the
phase difference $\varphi$ can be obtained by modifying the Hamiltonian
$H_{Lead}$ in equation (S4). To be more precise, we change the hopping term
$t$ to $te^{i\varphi}$, where $t$ is the hopping from the left side of the
self-energy lead to the right side.

Critical current under finite temperature $T$ can be calculated by using the
imaginary Green’s function. Consider the self-energy lead as a subsystem. The
current equals the change in the number of electron on the left side of the
lead.

$\displaystyle I=ie\left\langle\sum_{i\in
L}\frac{\mathrm{d}n_{i}}{\mathrm{d}\tau}\right\rangle=\frac{ie}{\hbar}\left\langle\sum_{i\in
L}[c_{i}^{\dagger}(\tau)c_{i}(\tau),H_{Lead}]\right\rangle.$ (S5)

Here we consider the imaginary time evolution and $i$ runs through the
positions to the left of the self-energy lead. For any diagonal term
$c^{\dagger}_{j}c_{j}$ in $H_{Lead}$, we have

$\displaystyle[c_{i}^{\dagger}c_{i},c^{\dagger}_{j}c_{j}]=0.$ (S6)

For $j,k\neq i$, we have

$\displaystyle[c_{i}^{\dagger}c_{i},c_{j}c_{k}]=[c_{i}^{\dagger}c_{i},c^{\dagger}_{j}c_{k}]=[c_{i}^{\dagger}c_{i},c_{j}c_{k}^{\dagger}]=0.$
(S7)

Therefore, to have the non-zero commutator, we need at least one operator
$c_{j}$ or $c_{j}^{\dagger}$ such that $j\in L$. Suppose $j,k\in L$, then we
have

$\displaystyle[c_{j}^{\dagger}c_{j},c_{j}^{\dagger}c_{k}]=c_{j}^{\dagger}c_{j}c_{j}^{\dagger}c_{k}-c_{j}^{\dagger}c_{k}c_{j}^{\dagger}c_{j}=c_{j}c_{k}-0=c_{j}^{\dagger}c_{k}$
(S8)
$\displaystyle[c_{k}^{\dagger}c_{k},c_{j}^{\dagger}c_{k}]=c_{k}^{\dagger}c_{k}c_{j}^{\dagger}c_{k}-c_{j}^{\dagger}c_{k}c_{k}^{\dagger}c_{k}=0-c_{j}c_{k}=-c_{j}^{\dagger}c_{k}$
(S9)

These two term cancel each other, thus only hopping between the left side and
right side of the lead contribute to the critical current. Then the equation
(S5) simplifies to

$\displaystyle I$ $\displaystyle=\frac{ie}{\hbar}\left\langle\sum_{i\in L,j\in
R}[c_{i}^{\dagger}(\tau)c_{i}(\tau),t_{ji}c_{i}^{\dagger}(\tau)c_{j}(\tau)-t_{ij}c_{j}^{\dagger}(\tau)c_{i}(\tau)]\right\rangle$
(S10) $\displaystyle=\frac{ie}{\hbar}\sum_{i\in L,j\in R}(t_{ji}\langle
c_{i}^{\dagger}(\tau)c_{j}(\tau)\rangle-t_{ij}\langle
c_{j}^{\dagger}(\tau)c_{i}\rangle(\tau))$ (S11)

By using the definition $G(\tau,\tau^{\prime})_{ij}=\langle\langle
c_{i}^{\dagger}(\tau)c_{j}(\tau^{\prime})\rangle$ for $\tau>\tau^{\prime}$ we
have

$\displaystyle I=\frac{ie}{\hbar}\sum_{i\in L,j\in
R}(t_{ij}G(\tau,\tau^{\prime})_{ij}-t_{ji}G(\tau,\tau^{\prime})_{ji}).$ (S12)

Then by apply the inverse Fourier transformation on the right hand side of
this equation and take the limit $\tau-\tau^{\prime}\to 0^{+}$, we get

$\displaystyle I$ $\displaystyle=\frac{ie}{\hbar}\sum_{i\in L,j\in
R}\sum_{n\in\mathbb{Z}}k_{B}T(t_{ji}e^{i\omega_{n}(\tau-\tau^{\prime})}G(i\omega_{n})_{ij}-t_{ij}e^{i\omega_{n}(\tau-\tau^{\prime})}G(i\omega_{n})_{ij})$
(S13) $\displaystyle=\frac{iek_{B}T}{\hbar}\sum_{n\in\mathbb{Z}}\sum_{i\in
L,j\in R}(t_{ji}G(i\omega_{n})_{ij}-t_{ij}G(i\omega_{n})_{ij})$ (S14)
$\displaystyle=\frac{-4ek_{B}T}{\hbar}\sum_{n\in\mathbb{N}}\mathrm{Im}\\{\mathrm{Tr}\left(T_{RL}G(i\omega_{n})_{LR}-T_{LR}G(i\omega_{n})_{RL}\right)\\},$
(S15)

where $T_{LR}$ and $T_{RL}$ are the hopping matrices from left (right) to
right (left), $\omega_{n}=(2n+1)\pi k_{B}T$ is the $n$-th Matsubara frequency
for electron. Factor $4$ on the last line comes from positive-negative
symmetry when sum over all the integers and particle-hole symmetry of the
system.

Figure S12: Cross section of the nanowire. Here we add a two-layer self energy
lead in the middle of the wire, where $t$ is the hopping term from the left
side of self energy lead to the right side.

##  XXI Parameters used in the simulations

In this section, we discuss the strength of parameters used in the simulation
that give the best match to experimental data.

###  XXI.1 Temperature

In the measurements, lattice temperature can only be estimated by reading the
temperature from the sensor on the dilution refrigerator mixing chamber plates
and it is varying from 50 to 60 $\mathrm{mK}$ while scanning the external
field. However, electron temperature is the relevant parameter for the
simulation. In our case, electrons are cooled down from room temperature to
base temperature by several stages of filters, especially the cooper powder
filter. However, the temperature of electrons is usually a bit higher than the
device temperature. So we simulate the skewed shape in different temperatures
(Fig.S13) and choose 100$\mathrm{mK}$ for all the simulation results presented
in the main text.

Figure S13:  (a) Skewed diffraction pattern from simulation when temperature T
= 50$\mathrm{mK}$, 100$\mathrm{mK}$, 150$\mathrm{mK}$ and 250$\mathrm{mK}$,
respectively. (b) Coefficient $\gamma$ as function of $B_{x}$ at different
temperatures.

From the simulation results, we also find that the maximum and minimum value
of coefficient $\gamma$ become smaller at higher temperature. This is
consistent with our temperature dependence measurement results [Fig. S17].

###  XXI.2 Strength of Spin-orbit interaction

The strength of Spin-orbit coupling is estimated by studying the critical
current diffraction pattern and choosing the one which best reproduced the
experiment results. We set chemical potential $\mu$ = 8 meV, which is
corresponding to two transverse or four spin-full modes, same value is used in
the [FIg. 5] in the main text. The skew shape we observe experimentally has
two features that we aim to reproduce: 1) the largest critical current $I_{c}$
is not located at zero field; 2) the largest critical current difference is
around $B_{x}=\pm 50mT$. Based on the skewed shape simulated with different
strengths of spin-orbit coupling, we find the skew is best reproduced with
$\alpha=200nm\cdot meV$. So we choose $\alpha=200$ for all simulation results
in the main text. Note that the true strength of spin-orbit interaction in
nanowires may differ from the parameters in a tight-binding simulation.

Figure S14:  (a)Critical current as function of $\hat{x}$-field when strength
of spin-orbit interaction $\alpha=0,60,120,180,240,300nm\cdot meV$. The
chemical potential $\mu=8meV$. (b) Coefficient $\gamma$ are extracted from (a)
and plotted as function of $B_{x}$ simulated with different $\alpha$. (c)At
same chemical potential, plot coefficient $\gamma$ versus perpendicular field
$B_{x}$ and spin-orbit strength $\alpha$ as a 2D map to study how alpha affect
skew shape.

###  XXI.3 Skew shape and field rotation simulation at another chemical
potential

Figure S15:  Numerical results from the KWANT simulation with same parameters
used in Fig. 5 but with another chemical potential $\mu=20meV$. This chemical
potential is corresponding to four transverse or eight spin-full modes. (a)
Critical current flow through each polarization as a function of magnetic
field $B_{x}$. (b) Coefficient $\gamma$ as a function of $B_{x}$, the
magnitude of $I_{c}$ is derived from (a). (c) Coefficient $\gamma$ as a
function of angle $\theta$ when the external field is rotating in three
orthogonal planes with fixed strength $|B|=50\mathrm{mT}$.

In the main text we present skewed diffraction pattern from Device A and
rotation from Device B (field rotation was not performed for device A). Thus
simulation at two separate chemical potential should be considered in
preparing this report. To maintain consistency of simulation, we choose
$\mu=8\mathrm{meV}$ for the simulation in the main text, Fig.4. Here we show
simulation results at $\mu=20\mathrm{meV}$ here.

##  XXII Current phase relation derived from simulation results

Figure S16:  (a) CPR when external field is along three directions related to
the device at field strength equals to 0T, 0.05T and 0.1T. (b) Amplitude of
Sine and Cosine terms that are derived from Fourier expansion of CPR when
external field is along $\hat{x}$-direction. The combined sine term is plotted
as function of $B_{x}$ field and order of harmonics. (c) Ground state phase
$\phi_{n0}$ at the first and the second order of combined sin harmonics as
function of external field $B_{x}$. A constant $\pi$ was subtracted from all
second order harmonics but has no effect as second order harmonic has a period
of $\pi$.

The current phase relation (CPR) is studied within the model and its
parameters suggested by comparison with experiment. We plot the CPR curve when
the strength of the external field is 0T, 0.05T, 0.1T in $B_{x}$, $B_{y}$ and
$B_{z}$ directions [Fig. S16]. Parameters used in this simulation is same as
that in Fig.4 in main text. In Fig.S16(a), we find that only when the external
field is along $\hat{x}$-direction, there is a shift of the ground state phase
in CPR. The numerical CPR curves are similar to those postulated in the
phenomenological model (Eq. 1 and Fig. 1).

To study how the time-reversal symmetry is broken when external field is along
$\hat{x}$-direction. We first perform a Fourier expansion with the simulated
CPR [Fig.S16 (b)]. We find there is a significant second order sin term in the
simulation. What’s more interesting is the first and second order cos terms
are also large compared to the sin term and they are first increasing when
field is applied but decreasing at higher field. This explains why the phase
of ground state shift by such a large value. It is known that harmonics
functions can be combined into a sine using the following:

$\displaystyle\sum_{i=1}A_{i}\sin(i\phi)+\sum_{i=0}B_{i}\cos(i\phi)=\sum_{i=1}A^{\prime}_{i}\sin(i\phi+\phi{i0})+Constant$
(S16)

Here we find constant is contributing less than $2\%$ of the combined function
in any case, so we drop it. We find the amplitude of first and second order of
combined sin function has a ratio about 4:1 when field strength is near zero.
This ratio is used in the minimal model presented in the main text.

Another interesting result is when we study the ground state phase $\phi_{i0}$
in the first and second harmonics, we find they are increasing linearly with
$B_{x}$ field within the range 0-100mT. This was mentioned in [6] but here we
provide more details. Based on the simulation, we get $\phi_{20}>2\phi_{10}$
and the grey shadow region is indicating the $\delta_{12}$ increase in with
field. Hence we can confirm there is a $\delta\phi$ term in the CPR from the
simulation. How this $\delta\phi$ related to the strength of spin-orbit
interaction is worth a further discussion in future works.

##  XXIII Temperature dependence of skewed diffraction patterns

Figure S17:  Temperature dependence in Device A. (a) For $V_{gate}=-2V$ and
external field $B_{x}=50mT$. dV/dI differential resistance is plotted as
function of current source I and temperature T. (b) Coefficient $\gamma$ is
plotted as function of temperature T. (c) Diffraction patterns taken at
T=1.1K, $\gamma$ is extracted from 2D scan results and plotted as function of
$B_{x}$

##  XXIV $\gamma$ trace in different field

Figure S18:  Coefficient $\gamma$ as function of $V_{gate}$ in different
external field $B$. (a) External field with fixed strength $|B|=50mT$ and
along $\hat{x}$, $\hat{y}$, and $\hat{z}$ axis. (b) External field with fixed
strength $|B|=100mT$ and along $\hat{x}$ and $\hat{z}$ axis. (c) The 2D
$B_{x}$ versus $V_{gate}$ map of coefficient gamma of Device C. The gate
voltage range is corresponding to the chemical potential $\mu=6-30meV$


$(G-X,\mathbf{v})\models\varphi[\mathbf{x}]$ and $Z\subseteq X$, there is
$j\in\\{1,\ldots,s\\}$ such that $u_{j}\in X\setminus X$. Because $X$ is a
colorful solution, $u_{j}\in B$ . Therefore, the subroutine calls
$\textsc{FindX}(Z^{\prime},h^{\prime})$ for $Z^{\prime}=Y\setminus\\{u_{j}\\}$
and $h^{\prime}=h-1$. It is easy to see that (a′) and (b′) are fulfilled.

Assume that $(F,\mathbf{v})\models\varphi[\mathbf{x}]$. Then ${\sf
depth}(Z)\geq k$. Because ${\sf depth}(X)\leq k-1$, $X$ has a representation
$(T,\alpha)$ with ${\sf depth}(T)\leq k-1$. Because ${\sf depth}(T)\leq k-1$
and ${\sf depth}(Z)\geq k$, there are vertices $x,y\in Z$ such that the nodes
$x^{\prime}=\alpha^{-1}(x)$ and $y^{\prime}=\alpha^{-1}(y)$ have the lowest
common ancestor $z$ in $T$ such that $z\neq x,y$ and it holds that
$\alpha(A_{T}(s))\setminus Z\neq\emptyset$ and $\alpha(A_{T}(z))$ is an
$(x,y)$-separator in $G$. In particular, $x$ and $y$ cannot be both in
$N_{G}[V(C)]$. By symmetry, we can assume that $x\notin N_{G}[V(C)]$. This
means that there is a red-yellow component $H$ of $F$ such that properties
(i)–(iv) of step 5 are fulfilled. Because $X$ is a colorful solution, we have
that $S=N_{F}[V(H)]\cap(B\cup Y)\subseteq X$. Thus, (a′) and (b′) are
fulfilled for $Z^{\prime}=Z\cup S$ and $h^{\prime}=h-|S|$. As the subroutine
calls $\textsc{FindX}(Z^{\prime},h^{\prime})$, we conclude that the claim if
fulfilled.

Recall that we call $\textsc{FindX}(\emptyset,p+k)$ and note that conditions
(a) and (b) are trivially fulfilled for $Z=\emptyset$ and $h=t$. Observe also
that in each recursive call of the subroutine the parameter $h$ strictly
decreases. Thus, we conclude that we output $X$ is some recursive call of
$\textsc{FindX}(Z,h)$. ∎

Lemma 16 concludes the description of the algorithm and its correctness proof.
We summarize and evaluate the running time in the following lemma.

######  Lemma 17.

Elimination Distance–(depth) to $\varphi$ on $(p,k)$-unbreakable graphs for
$\varphi\in\Sigma_{3}$ can be solved in
$2^{\mathcal{O}((p+k)(\log(p+k)+p))}\cdot n^{\mathcal{O}(|\varphi|)}$ time.

###### Proof.

If $|V(G)|\leq(3p+2k)(p+1)$, the problem is solved by brute force in
$2^{(3p+2k)(p+1)}\cdot((3p+2k)(p+1))^{\mathcal{O}(k+|\varphi|)}$ time. Assume
that $|V(G)|>(3p+2k)(p+1)$. Then we construct $\mathcal{F}$ in
$2^{\mathcal{O}((p+k)\log(p+k))}\cdot n\log n$ time. The size of $\mathcal{F}$
is at most $2^{\mathcal{O}((p+k)\log(p+k))}\cdot\log n$, and for every
$R\in\mathcal{F}$, we construct $\mathcal{F}^{\prime}$ in
$2^{\mathcal{O}(\min\\{p,k\\}\log(p+k))}\cdot n\log n$ time. Recall that the
size of $\mathcal{F}^{\prime}$ is at most
$2^{\mathcal{O}(\min\\{p,k\\}\log(p+k))}\cdot\log n$. Then we consider at most
$n^{r}$ $r$-tuples of vertices $\mathbf{v}$ that can be listed in
$n^{\mathcal{O}(|\varphi|)}$ time. Finally, for every $R\in\mathcal{F}$,
$Y\in\mathcal{F}^{\prime}$, and every $\mathbf{v}$, we call
$\textsc{FindX}(\emptyset,p+k)$.

Thus, it remains to evaluate the running time of
$\textsc{FindX}(\emptyset,p+k)$. Notice that in each call, $|Z|\leq p+k$. Then
we can verify in $(p+k)^{\mathcal{O}(k)}\cdot n^{\mathcal{O}(1)}$ time whether
${\sf depth}(Z)\leq k-1$ using Lemma 7. Also we can check whether
$(F,\mathbf{v})\models\varphi[\mathbf{x}]$ in $n^{\mathcal{O}(|\varphi|)}$ by
Observation 1. Simultaneously, we find an $s$-tuple $\mathbf{u}$ of vertices
of $F$ such that $(F,\mathbf{vu})\not\models\varphi[\mathbf{xy}]$ if this is
not the case. In step 4, we perform at most $s$ recursive calls. In step 5,
finding $H$ can be done in polynomial time. Notice that we have at most
$|Z|\leq p+k$ recursive calls in this step. The depth if the recursion is
upper bounded by $k+p$. This implies that the running time of
$\textsc{FindX}(\emptyset,p+k)$ is $(p+k)^{\mathcal{O}(p+k)}\cdot
n^{|\varphi|}$.

Summarizing, we obtain that the total running time is
$2^{\mathcal{O}((p+k)(\log(p+k)+p))}\cdot n^{\mathcal{O}(|\varphi|)}$. ∎

##  5 Lower bound for $\Pi_{3}$-formulas

In this section, we complement Theorem 1 by proving that there are formulas in
$\Pi_{3}$ for which Elimination Distance–($\star$) to $\varphi$ is
$\operatorname{{\sf W}}[2]$-hard. We state now Theorem 2 formally.

Theorem 2. _For every $\star\in\\{\text{\sf conn},\text{\sf prop},\text{\sf
depth}\\}$, there is $\varphi\in\Pi_{3}$ such that Elimination
Distance–($\star$) to $\varphi$ is $\operatorname{{\sf W}}[2]$-hard._

###### Proof.

We show that the problems are $\operatorname{{\sf W}}[2]$-hard for the formula
$\varphi$ expressing the property that for every vertex $u$ of a graph, there
is a vertex $v$ at distance at most two from $u$ of degree at most one. Notice
that the property that a vertex $v$ of a given graph $G$ has degree at most
one can be written as follows: for every $z_{1},z_{2}\in V(G)$, if $v$ is
adjacent to $z_{1}$ and $z_{2}$, then $z_{1}=z_{2}$. Thus, we define the
formula

$\psi(v,z_{1},z_{2})=[((v\sim z_{1})\wedge(v\sim
z_{2}))\rightarrow(z_{1}=z_{2})]$

with three free variables and set

$\displaystyle\varphi=\forall x\exists y_{1}\exists y_{2}\forall z_{1}\forall
z_{2}~{}$ $\displaystyle[\psi(x,z_{1},z_{2})\vee((x\sim
y_{1})\wedge\psi(y_{1},z_{1},z_{2}))$ $\displaystyle\vee((x\sim
y_{1})\wedge(y_{1}\sim y_{2})\wedge\psi(y_{2},z_{1},z_{2}))].$

Clearly, $\varphi\in\Pi_{3}$.

To show $\operatorname{{\sf W}}[2]$-hardness, we reduce from the Set Cover
problem. The problem asks, given a universe $U$, a family $\mathcal{S}$ of
subsets of $\mathcal{S}$, and a positive integer $k$, whether there is
$\mathcal{S}^{\prime}\subset\mathcal{S}$ of size at most $k$ that covers $U$,
that is, for every $u\in U$, there is $S\in\mathcal{S}^{\prime}$ such that
$u\in S$. It is well-known that Set Cover is $\operatorname{{\sf
W}}[2]$-complete when parameterized by $k$ [13].

$w_{2}$$u_{2}^{(p)}$$u_{3}^{(p)}$$u_{1}^{(p)}$$u_{4}^{(p)}$$s_{1}$$v_{1}$$w_{1}$$s_{2}$$v_{2}$

Figure 2: Construction of $G$ for $n=2$ and $m=2$ with
$S_{1}=\\{u_{1},u_{2},u_{3}\\}$ and $S_{2}=\\{u_{2},u_{3},u_{4}\\}$; for
simplicity, just one copy of each $u_{i}^{(p)}$ for $p\in\\{0,\ldots,k\\}$ is
shown.

Let $(U,\mathcal{S},k)$ be an instance of Set Cover with
$U=\\{u_{1},\ldots,u_{n}\\}$, $\mathcal{S}=\\{S_{1},\ldots,S_{m}\\}$. We also
assume that $n\geq 2$ and $k\leq m$. We construct the following graph $G$ (see
Figure 2).

  * •

For every $i\in\\{1,\ldots,n\\}$, construct $k+2$ vertices
$u_{i}^{(1)},\ldots,u_{i}^{(k+2)}$, and then for every
$i,j\in\\{1,\ldots,n\\}$ and all $p,q\in\\{1,\ldots,k+2\\}$ such that
$(i,p)\neq(j,q)$, make $u_{i}^{(p)}$ and $u_{j}^{(q)}$ adjacent.

  * •

For every $j\in\\{1,\dots,m\\}$, construct three vertices $s_{j},v_{j},w_{j}$
and edges $s_{j}v_{j}$ and $v_{j}w_{j}$.

  * •

For every $i\in\\{1,\ldots,n\\}$ and every $j\in\\{1,\ldots,m\\}$, make
$s_{j}$ adjacent to $u_{i}^{(1)},\ldots,u_{i}^{(k+2)}$ if $u_{i}\in S_{j}$.

We claim that $G$ has a set cover of size at most $k$ if and only if ${\sf
ed}_{\varphi}^{\star}(G)\leq k$ for $\star\in\\{\text{\sf conn},\text{\sf
prop},\text{\sf depth}\\}$. Notice that by the definition of $\varphi$,
$H\models\varphi$ if and only if $C\models\varphi$ for every component $C$ of
$H$. Therefore, ${\sf ed}_{\varphi}^{\text{\sf conn}}(G)={\sf
ed}_{\varphi}^{\text{\sf prop}}(G)={\sf ed}_{\varphi}^{\text{\sf depth}}(G)$,
and it is sufficient to prove that $G$ has a set cover of size at most $k$ if
and only if there is an elimination set $X$ of $G$ with ${\sf depth}(G)\leq
k-1$ such that $C\models\varphi$ for every component $C$ of $G-X$.

Suppose that sets $S_{j_{1}},\ldots,S_{j_{k}}\in\mathcal{S}$ form a set cover.
We define $X=\\{w_{j_{1}},\ldots,w_{j_{k}}\\}$. Since $|X|=k$, ${\sf
depth}(X)\leq k-1$. Notice that $H=G-X$ is connected. We claim that
$H\models\varphi$. Recall that $H\models\varphi$ if and only if for every
vertex $x\in V(H)$ there is a vertex $y\in V(H)$ at distance at most two such
that $d_{H}(y)\leq 1$. This property trivially holds if
$x\in\\{s_{j},v_{j},w_{j}\\}\setminus X$ for $j\in\\{1,\ldots,m\\}$. Consider
a vertex $u_{i}^{(p)}$ for some $i\in\\{1,\ldots,n\\}$ and
$p\in\\{1,\ldots,k+2\\}$. We have that there is $h\in\\{1,\ldots,k\\}$ such
that $u_{i}\in S_{j_{h}}$. Then $u_{i}s_{j_{h}}\in E(H)$. Because
$w_{j_{h}}\in X$, we obtain that $d_{H}(v_{j_{h}})=1$. Since
$s_{j_{h}}v_{j_{h}}\in E(H)$, $v_{j_{h}}$ is at distance two from
$u_{i}^{(p)}$ as required. We conclude that $H\models\varphi$.

For the opposite direction, assume that there is an elimination set $X$ of $G$
with ${\sf depth}(X)\leq k-1$ such that $C\models\varphi$ for every component
$C$ of $G-X$. Consider $Z=\\{u_{i}^{(p)}\mid 1\leq i\leq n,~{}1\leq p\leq
k+2\\}$. Because $Z$ is a clique, we have that $|Z\cap X|\leq k$. To see this,
consider a representation $(T,\alpha)$ of $X$ with ${\sf depth}(T)\leq k-1$.
Then there is a leaf $x$ of $T$ such that $\alpha^{-1}(X\cap Z)\subseteq
A_{T}(x)$. Since ${\sf depth}(T)\leq k-1$, we conclude that $|Z\cap X|\leq k$.
Note that the vertices of $Z\setminus X$ are in the same component $H$ of
$G-X$. Let $W=\\{w_{1},\ldots,w_{m}\\}$. By Observation 3, $|N_{G}(V(H))\cap
X|\leq k$. Hence, $|N_{G}(V(H))\cap(X\cap W)|\leq k$ as well. Let
$\\{w_{j_{1}},\ldots,w_{j_{\ell}}\\}=N_{G}(V(H))\cap(X\cap W)$. We claim that
the sets $S_{j_{1}},\ldots,S_{j_{k}}$ cover $U$.

Consider an arbitrary $i\in\\{1,\ldots,n\\}$. Because $|Z\cap W|\leq k$, there
are two distinct $p,q\in\\{1,\ldots,k+2\\}$ such that
$u_{i}^{(p)},u_{i}^{(q)}\in V(H)$. Since $H\models\varphi$, there is a vertex
$z\in V(H)$ at distance at most two from $u_{i}^{(p)}$ such that $d_{H}(z)\leq
1$. Since $n\geq 2$, we have that $|Z\setminus X|\geq 3$ and, therefore,
$d_{H}(u_{i}^{(p)})\geq 2$. Moreover, for every $h\in\\{1,\ldots,n\\}$ and
$r\in\\{1,\ldots,k+2\\}$, if $u_{h}^{(r)}\in V(H)$, then
$d_{H}(u_{h}^{(r)})\geq 2$. Then the construction of $G$ implies that there is
$h\in\\{1,\ldots,m\\}$ such that $s_{h}\in V(H)$ and $u_{i}^{(p)}s_{h}\in
E(G)$ with the property that either $d_{H}(s_{j})\leq 1$ or $s_{j}$ has a
neighboor in $H$ of degree at most one. As $s_{h}$ is adjacent to
$u_{i}^{(p)}$, this vertex is adjacent to $u_{i}^{(q)}$. Hence,
$d_{H}(s_{h})\geq 2$. We obtain that $v_{h}\in V(H)$ and $d_{H}(v_{h})\leq 1$
by the construction of $G$. This means that $w_{h}\notin H$, that is,
$w_{h}\in N_{G}(V(H))\cap(X\cap W)$. We conclude that there is
$t\in\\{1,\ldots,\ell\\}$ such that $j_{t}=h$. Finally, because $u_{i}^{(p)}$
is adjacent to $s_{j_{t}}$, $u_{i}\in S_{j_{t}}$ and this concludes the proof.
∎

##  6 Discussion

We established a parameterized complexity dichotomy for the elimination
problems whose aim is to satisfy a FOL formula $\varphi$ with respect to the
quantification structure of the prefix. For this, we considered three variants
of the elimination distance to the class of graphs modelling $\varphi$ and
defined the Elimination Distance–($\star$) to $\varphi$ for
$\star\in\\{\text{\sf conn},\text{\sf prop},\text{\sf depth}\\}$ corresponding
to the considered type of distance. In Theorem 1, we proved that for every FOL
formula $\varphi\in\Sigma_{3}$, Elimination Distance–($\star$) to $\varphi$ is
$\operatorname{{\sf FPT}}$ for $\star\in\\{\text{\sf conn},\text{\sf
prop},\text{\sf depth}\\}$. In Theorem 2, we showed that this result is tight
in the sense that there are FOL formulas $\varphi\in\Pi_{3}$ such that these
problems are $\operatorname{{\sf W}}[2]$-hard.

Notice that the above dichotomy is the same for all the considered variants of
the elimination problems. Moreover, it coincides with the structural dichotomy
obtained by for Deletion to $\varphi$ by Fomin, Golovach, and Thilikos in
[15]. This leads to the following natural question: is there a FOL formula
$\varphi$ such that the parameterized complexity of Elimination
Distance–($\star$) to $\varphi$ for $\star\in\\{\text{\sf conn},\text{\sf
prop},\text{\sf depth}\\}$ and Deletion to $\varphi$ differs? In particular,
is there a formula $\varphi$ such that Deletion to $\varphi$ is
$\operatorname{{\sf FPT}}$ but one of the problem Elimination
Distance–($\star$) to $\varphi$ for $\star\in\\{\text{\sf conn},\text{\sf
prop},\text{\sf depth}\\}$ turns to be, say $\operatorname{{\sf W}}[1]$ or
$\operatorname{{\sf W}}[2]$-hard? Note that Lemma 5 holds for every FOL
formula $\varphi$. Thus, solving Elimination Distance–($\star$) to $\varphi$
for $\star\in\\{\text{\sf conn},\text{\sf prop},\text{\sf depth}\\}$ can be
reduced to solving these problems on unbreakable graphs by Theorem 3. Since
Elimination Distance–($\star$) to $\varphi$ for $\star\in\\{\text{\sf
conn},\text{\sf prop},\text{\sf depth}\\}$ is somehow similar to Deletion to
$\varphi$ on unbreakable graphs, it may happen that Elimination
Distance–($\star$) to $\varphi$ for $\star\in\\{\text{\sf conn},\text{\sf
prop},\text{\sf depth}\\}$ are $\operatorname{{\sf FPT}}$ whenever Deletion to
$\varphi$ is $\operatorname{{\sf FPT}}$. However proving this would demand
applying different algorithmic tools as our techniques are tailored for
$\varphi\in\Sigma_{3}$. Also it would be interesting to know whether there are
FOL formulas such that Elimination Distance–($\star$) to $\varphi$ for
$\star\in\\{\text{\sf conn},\text{\sf prop},\text{\sf depth}\\}$ differ from
the parameterized complexity viewpoint.

In contrast with the same behaviour of the elimination and deletion problems
with respect to the inclusion in $\operatorname{{\sf FPT}}$, we would like to
point that they behave differently with respect to _kernelization_ (we refer
to the books [10, 16] for the definition of the notion). It was shown in [15]
that Deletion to $\varphi$ admits a polynomial kernel for
$\varphi\in\Sigma_{1}\cup\Pi_{1}$ (in fact, Deletion to $\varphi$ is
polynomial for $\varphi\in\Sigma_{1}$) and there are formulas
$\varphi\in\Pi_{2}$ and $\Sigma_{2}$ such that Deletion to $\varphi$ has no
polynomial kernel unless $\operatorname{{\sf NP}}\subseteq\operatorname{{\sf
coNP}}/{\sf poly}$. For the elimination problems, we can show the following
lower bound.

######  Proposition 1.

There are formulas $\varphi\in\Pi_{1}$ such that Elimination
Distance–($\star$) to $\varphi$ do not admit polynomial kernels unless
$\operatorname{{\sf NP}}\subseteq\operatorname{{\sf coNP}}/{\sf poly}$ for
$\star\in\\{\text{\sf conn},\text{\sf prop},\text{\sf depth}\\}$.

###### Proof.

We show the claim for the formula $\varphi$ expressing the property that a
graph has no triangles, that is, cycles of length three:

$\displaystyle\varphi=\forall x\forall y\forall
z~{}[(x=y)\vee(y=z)\vee(x=z)\vee\neg(x\sim y)\vee\neg(y\sim z)\vee\neg(x\sim
z)].$

It is straightforward to see that $G\models\varphi$ if and only if $G$ has no
triangles.

By the classical results of Lewis and Yannakakis [22], Deletion to $\varphi$
is $\operatorname{{\sf NP}}$-complete. Then it is easy to observe that the
problem remains $\operatorname{{\sf NP}}$ on instances $(G,k)$, where $G$ is a
$(k+1)$-connected graph. For example, we can reduce from Deletion to $\varphi$
on general graphs. Let $G$ be an $n$-vertex graph. We assume that $k<n-1$ as
otherwise the problem is trivial. We construct the graph $G^{\prime}$ from $G$
by adding $k+1$ copies of the complete bipartite graph $K_{n,n}$ and making
each vertex of one part of the vertex partition to a unique vertex of $G$.
Clearly, $G^{\prime}$ is $(k+1)$-connected and it is easy to see that $G-X$ is
triangle-free if and only if $G^{\prime}$ has no triangles for every
$X\subseteq V(G^{\prime})$. This proves the $\operatorname{{\sf NP}}$-hardness
for Deletion to $\varphi$ on $(k+1)$-connected graphs. Then Observation 4
implies that Elimination Distance–($\star$) to $\varphi$ is
$\operatorname{{\sf NP}}$-complete for every $\star\in\\{\text{\sf
conn},\text{\sf prop},\text{\sf depth}\\}$.

Let $(G_{1},k),\ldots(G_{t},k)$ be instances of Elimination Distance–($\star$)
to $\varphi$ for some $\star\in\\{\text{\sf conn},\text{\sf prop},\text{\sf
depth}\\}$. Let $G$ be the disjoint union of $G_{1},\ldots,G_{t}$. Then for
$\star\in\\{\text{\sf conn},\text{\sf prop}\\}$, we have that $(G,k)$ is a
yes-instance of Elimination Distance–($\star$) to $\varphi$ if and only if
$(G_{j},k)$ is a yes-instance of Elimination Distance–($\star$) to $\varphi$
for every $j\in\\{1,\ldots,t\\}$. Then by the result of Bodlaender, Jansen,
and Kratsch [3] (see also [16, Part III] for the introduction to the
technique), Elimination Distance–($\star$) to $\varphi$ does not admit a
polynomial kernel unless $\operatorname{{\sf NP}}\subseteq\operatorname{{\sf
coNP}}/{\sf poly}$. For Elimination Distance–(depth) to $\varphi$, consider
$G^{\prime}$ that is the disjoint union of $G$ and $K_{k+3}$. Clearly, ${\sf
ed}_{\varphi}^{\text{\sf depth}}(K_{k+2})=k+1$. Then by Lemma 2,
$(G^{\prime},k+1)$ is a yes-instance of Elimination Distance–(depth) to
$\varphi$ if and only if $(G_{j},k)$ is a yes-instance of Elimination
Distance–(depth) to $\varphi$ for every $j\in\\{1,\ldots,t\\}$. This implies
that Elimination Distance–(depth) to $\varphi$ has no polynomial kernel unless
$\operatorname{{\sf NP}}\subseteq\operatorname{{\sf coNP}}/{\sf poly}$. ∎

Notice that Proposition 1 does no exclude existence of _Turing kernels_ (we
again refer to [10, 16] for the definition of the notion). This makes it
natural to ask whether Elimination Distance–($\star$) to $\varphi$ admit
polynomial Turing kernels for $\varphi\in\Sigma_{3}$ for $\star\in\\{\text{\sf
conn},\text{\sf prop},\text{\sf depth}\\}$.

We defined the depth of a set $X\subseteq V(G)$ using a representation.
However, there is an equivalent definition that uses the notion of _tree-
depth_ (see, e.g., [5] for the definition). Let $G$ be a graph and let
$X\subseteq V(G)$. We define the _torso_ of $X$ as the graph $H$ obtained from
$G[X]$ by making every two vertices $u,v\in X$ adjacent if three is a
component $C$ of $G-X$ such that $u,v\in N_{G}(V(C))$. Then the following
property can be shown by the definition of tree-depth.

######  Observation 5.

For a set $X\subseteq V(G)$ and an integer $k$, ${\sf depth}(X)\leq k-1$ if
and only if the tree-depth of the torso of $X$ is at most $k$.

Thus, Elimination Distance–(depth) to $\varphi$ can be stated as follows:
given a graph $G$ and a nonnegative integer $k$, is there $X\subseteq V(G)$
whose torso has the tree-depth at most $k$ such that $G-X\models\varphi$? In
other words, we ask whether there is a set of vertices whose torso has bounded
tree-depth such that the graph obtained by the deletion of this set models our
formula. Then we can consider the variants of Elimination Distance–(depth) to
$\varphi$ for other “width-measures”. For example, what can be said about
parameterized complexity of the variant of Elimination Distance–(depth) to
$\varphi$, where the _tree-width_ (see, e.g, [10] for the defintion) of the
torso of $X$ should be at most $k-1$?

Finally, we believe that it could be interesting to consider yet another
variant of the elimination distance. Recall that in the definitions of ${\sf
ed}_{\varphi}^{\star}$ for $\star\in\\{\text{\sf conn},\text{\sf
prop},\text{\sf depth}\\}$, we considered properties of the components. In
particular,

${\sf ed}_{\varphi}^{\text{\sf prop}}(G)=\begin{cases}0,&\mbox{if
}G\models\varphi,\\\ 1+\min_{v\in V(G)}{\sf ed}_{\varphi}^{\text{\sf
prop}}(G-v),&\mbox{if\\! }G\not\models\varphi\text{\\! and\\! }G\text{\\!
is\\! connected},\\\ \max\\{1,\max\\{{\sf ed}_{\varphi}^{\text{\sf
prop}}(C)\mid C\text{ is a component of
}G\\}\\},&\mbox{otherwise}.\end{cases}$

However, we can consider _unions_ of components instead. We say that graphs
$G_{1},\ldots,G_{s}$ form a _component-partition_ of $G$ if every component of
$G$ is a component of $G_{i}$ for some $i\in\\{1,\ldots,s\\}$ and $G$ is the
disjoint union of $G_{1},\ldots,G_{s}$. Then, we can define

${\sf ed}_{\varphi}^{\text{\sf part}}(G)=\begin{cases}0,&\mbox{if\\!
}G\models\varphi,\\\ 1+\min_{v\in V(G)}{\sf ed}_{\varphi}^{\text{\sf
part}}(G-v),&\mbox{if\\! }G\not\models\varphi\text{\\! and\\! }G\text{\\!
is\\! connected},\\\ \min\\{\max\\{1,{\sf ed}_{\varphi}^{\text{\sf
part}}(G_{1}),\ldots,{\sf ed}_{\varphi}^{\text{\sf part}}(G_{s})\\}\mid
G_{1},\ldots,G_{s}&\\\ \hskip 48.36967pt\text{is a component partition of
}G\\}&\mbox{otherwise}.\end{cases}$

Then we can define the respective Elimination Distance–(part) to $\varphi$ and
investigate its parameterized complexity depending of $\varphi$. Note that our
approach for solving the elimination problems fails in this case. In
particular, we cannot express the problem using MSOL.

## References

  * [1] A. Agrawal, L. Kanesh, F. Panolan, M. S. Ramanujan, and S. Saurabh, An FPT algorithm for elimination distance to bounded degree graphs, in 38th International Symposium on Theoretical Aspects of Computer Science, (STACS), vol. 187 of LIPIcs, Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2021, pp. 5:1–5:11.
  * [2] A. Agrawal and M. S. Ramanujan, On the parameterized complexity of clique elimination distance, in 15th International Symposium on Parameterized and Exact Computation (IPEC), vol. 180 of Leibniz International Proceedings in Informatics (LIPIcs), Dagstuhl, Germany, 2020, Schloss Dagstuhl–Leibniz-Zentrum für Informatik, pp. 1:1–1:13.
  * [3] H. L. Bodlaender, B. M. P. Jansen, and S. Kratsch, Kernelization lower bounds by cross-composition, SIAM J. Discrete Math., 28 (2014), pp. 277–305.
  * [4] E. Börger, E. Grädel, and Y. Gurevich, The classical decision problem, Springer Science & Business Media, 2001.
  * [5] J. Bulian and A. Dawar, Graph isomorphism parameterized by elimination distance to bounded degree, Algorithmica, 75 (2016), pp. 363–382.
  * [6] , Fixed-parameter tractable distances to sparse graph classes, Algorithmica, 79 (2017), pp. 139–158.
  * [7] L. Cai, S. M. Chan, and S. O. Chan, Random separation: A new method for solving fixed-cardinality optimization problems, in Proceedings of the 2nd International Workshop Parameterized and Exact Computation (IWPEC), vol. 4169 of Lecture Notes in Computer Science, Springer, 2006, pp. 239–250.
  * [8] R. Chitnis, M. Cygan, M. Hajiaghayi, M. Pilipczuk, and M. Pilipczuk, Designing FPT algorithms for cut problems using randomized contractions, SIAM J. Comput., 45 (2016), pp. 1171–1229.
  * [9] B. Courcelle and J. Engelfriet, Graph Structure and Monadic Second-Order Logic - A Language-Theoretic Approach, vol. 138 of Encyclopedia of mathematics and its applications, Cambridge University Press, 2012\.
  * [10] M. Cygan, F. V. Fomin, L. Kowalik, D. Lokshtanov, D. Marx, M. Pilipczuk, M. Pilipczuk, and S. Saurabh, Parameterized Algorithms, Springer, 2015\.
  * [11] M. Cygan, D. Lokshtanov, M. Pilipczuk, M. Pilipczuk, and S. Saurabh, Minimum bisection is fixed-parameter tractable, SIAM J. Comput., 48 (2019), pp. 417–450.
  * [12] R. Diestel, Graph Theory, 4th Edition, vol. 173 of Graduate texts in mathematics, Springer, 2012.
  * [13] R. G. Downey and M. R. Fellows, Fundamentals of Parameterized Complexity, Texts in Computer Science, Springer, 2013.
  * [14] J. Flum and M. Grohe, Parameterized Complexity Theory, Texts in Theoretical Computer Science. An EATCS Series, Springer, 2006.
  * [15] F. V. Fomin, P. A. Golovach, and D. M. Thilikos, On the parameterized complexity of graph modification to first-order logic properties, Theory Comput. Syst., 64 (2020), pp. 251–271.
  * [16] F. V. Fomin, D. Lokshtanov, S. Saurabh, and M. Zehavi, Kernelization, Cambridge University Press, Cambridge, 2019.  Theory of parameterized preprocessing.
  * [17] M. Frick and M. Grohe, The complexity of first-order and monadic second-order logic revisited, Ann. Pure Appl. Logic, 130 (2004), pp. 3–31.
  * [18] G. Gottlob, P. G. Kolaitis, and T. Schwentick, Existential second-order logic over graphs: Charting the tractability frontier, J. ACM, 51 (2004), pp. 312–362.
  * [19] , Existential second-order logic over graphs: Charting the tractability frontier, J. ACM, 51 (2004), pp. 312–362.
  * [20] J. Guo, F. Hüffner, and R. Niedermeier, A structural view on parameterizing problems: Distance from triviality, in Proceedings of the 1st International Workshop Parameterized and Exact Computation (IWPEC), vol. 3162 of Lecture Notes in Computer Science, Springer, 2004, pp. 162–173.
  * [21] E.-M. C. Hols, S. Kratsch, and A. Pieterse, Elimination Distances, Blocking Sets, and Kernels for Vertex Cover, in 37th International Symposium on Theoretical Aspects of Computer Science (STACS), vol. 154 of Leibniz International Proceedings in Informatics (LIPIcs), Dagstuhl, Germany, 2020, Schloss Dagstuhl–Leibniz-Zentrum für Informatik, pp. 36:1–36:14.
  * [22] J. M. Lewis and M. Yannakakis, The node-deletion problem for hereditary properties is NP-complete, J. Comput. Syst. Sci., 20 (1980), pp. 219–230.
  * [23] A. Lindermayr, S. Siebertz, and A. Vigny, Elimination distance to bounded degree on planar graphs, in Proceeding of the 45th International Symposium on Mathematical Foundations of Computer Science (MFCS), vol. 170 of LIPIcs, Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2020, pp. 65:1–65:12.
  * [24] D. Lokshtanov, M. S. Ramanujan, S. Saurabh, and M. Zehavi, Reducing CMSO model checking to highly connected graphs, in 45th International Colloquium on Automata, Languages, and Programming (ICALP), vol. 107 of LIPIcs, Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2018, pp. 135:1–135:14.
  * [25] , Reducing CMSO model checking to highly connected graphs, CoRR, abs/1802.01453 (2018).
  * [26] J. Nešetřil and P. Ossona de Mendez, Tree-depth, subgraph coloring and homomorphism bounds, European J. Combin., 27 (2006), pp. 1022–1041.
  * [27] R. Niedermeier, Invitation to fixed-parameter algorithms, vol. 31 of Oxford Lecture Series in Mathematics and its Applications, Oxford University Press, 2006.
  * [28] C. Smorynski, The incompleteness theorems, in Handbook of mathematical logic, vol. 90 of Stud. Logic Found. Math., North-Holland, Amsterdam, 1977, pp. 821–865.
  * [29] M. Y. Vardi, The complexity of relational query languages (extended abstract), in Proceedings of the 14th Annual ACM Symposium on Theory of Computing, May 5-7, 1982, San Francisco, California, USA, ACM, 1982, pp. 137–146.
  * [30] R. Williams, Faster decision of first-order graph properties, in Joint Meeting of the Twenty-Third EACSL Annual Conference on Computer Science Logic (CSL) and the Twenty-Ninth Annual ACM/IEEE Symposium on Logic in Computer Science (LICS), ACM, 2014, pp. 80:1–80:6.