# Statistics of Energy in Isothermal Supersonic Turbulence

Branislav Rabatin Department of Physics, Florida State University,
Tallahassee, Florida  David C. Collins Department of Physics, Florida State
University, Tallahassee, Florida

###### Abstract

Turbulence is a key process in many astrophysical systems. In this work we
explore the statistics of thermal and kinetic energy of isothermal,
supersonic, turbulent gas. We develop analytic formulas for the PDF of thermal
and kinetic energies and their joint PDF. We compare these analytical models
with a suite of simulations with a fixed resolution of $1024^{3}$ cells across
4 different Mach numbers (1, 2, 4, 8) and three different driving patterns
(compressive, mixed, solenoidal). We discover an interesting jump
discontinuity in the thermal energy PDF, which carries onto the joint PDF.

star formation, ISM, turbulence

##  1 Introduction

Supersonic turbulence is a dynamical mechanism governing many aspects of
galaxy evolution. Statistical modelling of turbulent gas offers insight into
the dynamical evolution of the interstellar medium (ISM) and star formation.
Supersonic motion is characteristic for star-forming regions, wherein the
occurence of local collapse is triggered by large density fluctuations typical
of supersonic flows (Scalo & Elmegreen 2004, Mac Low & Klessen 2004) Turbulent
molecular clouds in astrophysical simulations are commonly treated as
isothermal. This expedient simplification is warranted due to the relatively
constant temperature of the molecular clouds maintainted by their effective
cooling rates.

In this study we set out to explore the statistical properties of thermal and
kinetic energy involved in the dynamics of supersonic isothermal turbulence by
describing the probability density function (PDF) of each component, as well
as the joint statistics of the two. Through this analysis, further
understanding of the statistical nature of isothermal supersonic turbulence is
obtained. This focus opens further possibilities to model the complex dynamics
of the ISM, extended to magnetized, self-gravitating turbulence with more
realistic thermal evolution.

This is the third paper in a series. The first paper (Rabatin & Collins,
2023b) explored the PDF of density. This is treated as a lognormal in many
works (e.g. Vazquez-Semadeni, 1994). A lognormal is reasonable, as the density
in supersonic isothermal turbulence can be treated as the result of a series
of shocks, and an infinite number of random multiplications results in a
lognormal distribution. If one instead considers that a parcel of gas
experiences only a small finite number of shocks, one obtains the Finite Shock
Model (FSM). In Rabatin & Collins (2023b) we show that FSM fits simulated data
quite well, and naturally explains the variation with forcing parameter. In
the second paper (Rabatin & Collins, 2023a), we examine the joint distribution
of density and velocity, and develop an analytic model for their correlation.
Often assumed to be uncorrelated, we find a correlation that depends on Mach
number, and can be described as a combination of differently weighted
marginalized distributions.

This paper is organized into several parts. In Sec. 2 we describe the
simulation code, parameters, and numerical data analysis. In Sec. 3 we review
the dynamical quantities describing isothermally turbulent medium - density
$\rho$, and velocity $\vec{v}$, as well as their statistical properties. In
Sec. 4 we show, that the thermal energy density per unit volume of isothermal
medium only depends on density. Therefore, its PDF, discussed in detail in
Sec. 4.4, can be built on our previously developed Finite Shock Model of
density. In Sec. 5 we focus on the kinetic energy density per unit volume.
Properties of the joint PDF of thermal and kinetic energy are discussed in
Sec. 6. We conclude in Sec. 7.

##  2 Methods

In this section we describe simulations that were used to verify our new model
and present the methods used for statistical analysis.

###  2.1 Simulations

To test our theoretical results a set of turbulent, supersonic, isothermal
simulations were performed using the hydrodynamic code Enzo (Bryan et al.,
2014). The piecewise parabolic method (Woodward & Colella, 1984) is used to
capture the shocks naturally occurring in a supersonically turbulent medium.
The simulations are performed within a cube of unit length $L$, periodic in
all directions, at a resolution of $1024^{3}$ cells, each having the same
volume $\delta V=1/1024^{3}$.

To maintain a steady, turbulent state, the Stochastic forcing module
(Federrath et al., 2008) is used to drive the medium by adding energy at the
large scale, at a rate proportional to the dissipation, $\dot{E}\sim M^{3}/L$,
where $M$ is the Mach number (Mac Low, 1999; Mac Low & Klessen, 2004). The
large-scale driving pattern is modelled after the Uhlenbeck & Ornstein (1930)
process. In addition to the 1D r.m.s. Mach number $\mathcal{M}_{\rm 1D}$, each
simulation is also characterized by the forcing mode $\xi$. The forcing mode
determines the proportion of the solenoidal components within the forcing
field; a purely solenoidal driving corresponds to $\xi=1$, while $\xi=0$
represents compressive driving.

For the purpose of extracting statistical properties from the simulations we
define the turnover time scale $\tau$. The turnover scale is roughly equal to
the turbulent crossing time $\tau_{\text{turb.}}=(L/2)/\mathcal{M}_{\rm 1D}$,
where $L/2$ is the size of the driving pattern and $\mathcal{M}_{\rm 1D}$ is
the 1D r.m.s. Mach number. Two frames, when taken at times separated by at
least $\tau$ are said to be statistically uncorrelated. All simulations are
run for 9 turnover scales, while the individual snapshots are taken with the
step of $\tau/10$. To capture the steady turbulent state, frames with
$t<2\tau$ are discarded, leaving 71 snapshots available for statistical
analysis. The 71 frames are analyzed together to provide a robust statistical
ensemble. A similar approach is common in other works (Porter et al. (1999),
Porter & Woodward (2000), Federrath et al. (2010), Federrath (2013), Federrath
et al. (2021)).

Numerical simulations were performed using three forcing modes, compressive,
mixed and solenoidal ($\xi=0,1/2,1$) and four target Mach numbers,
$\mathcal{M}_{\rm 1D}=1,2,4,8$. Simulation parameters and other extracted
statistical quantities are summarized in Table 1. The first column of the
table identifies each simulation in the form of $\xi-M$, with forcing $\xi$
and target Mach number $M$. The actual 1D r.m.s. Mach number $\mathcal{M}_{\rm
1D}$ is listed in the second column, followed by the ratio of volume-weighted
Mach number to mass-weighted Mach number squared, $\mathfrak{X}$. The third
and last columns show the average thermal and kinetic energy, $\langle
F\rangle=\mu_{M}$, $\langle K\rangle=\frac{3}{2}\mathcal{M}_{M}^{2}$.

Speed of sound $c_{s}=\sqrt{p/\rho}$ and mean density $\rho_{0}=M/V$ are
simulation constants which will be throughout this work, for the sake of
brevity, set to 1.

Sim. ($\xi-M$) | $\mathcal{M}_{\rm 1D}$ | $\mathfrak{X}$ | $\left\langle F\right\rangle$ | $\left\langle K\right\rangle$
---|---|---|---|---
$0-1$ | $0.977$ | $1.35$ | $0.910$ | $1.06$
$0-2$ | $1.99$ | $1.38$ | $1.75$ | $4.31$
$0-4$ | $3.92$ | $1.43$ | $2.49$ | $16.1$
$0-8$ | $7.77$ | $1.41$ | $3.14$ | $64.1$
$1/2-1$ | $0.999$ | $1.14$ | $0.191$ | $1.31$
$1/2-2$ | $2.00$ | $1.21$ | $0.588$ | $4.99$
$1/2-4$ | $3.98$ | $1.16$ | $1.03$ | $20.5$
$1/2-8$ | $7.89$ | $1.14$ | $1.37$ | $81.8$
$1-1$ | $0.993$ | $1.15$ | $0.160$ | $1.29$
$1-2$ | $1.98$ | $1.16$ | $0.489$ | $5.08$
$1-4$ | $3.71$ | $1.14$ | $0.840$ | $18.2$
$1-8$ | $8.04$ | $1.13$ | $1.12$ | $86.1$
Table 1: Simulation parameters. The first column denotes each simulation in
the form of $\xi-M$ where $\xi$ is the forcing mode and $M$ is the nominal 1D
r.m.s. Mach number. The second column lists the measured 1D r.m.s. Mach
number. The third column represents $\mathfrak{X}=\langle
v^{2}\rangle/\langle\rho v^{2}\rangle$, the ratio between the volume- and
mass-weighted Mach numbers, squared. Last two columns show the average thermal
and kinetic energy, $\langle F\rangle,\langle K\rangle$. The average thermal
energy is related to the mass-weighted mean of log density, $\mu_{M}$ (see
Sec. 4.3), whereas the average kinetic energy is related to the mass-weighted
Mach number, $\langle K\rangle=\frac{3}{2}\mathcal{M}_{M}^{2}$.

###  2.2 Statistical analysis

Any random quantity $Q(\vec{x})$ defined within a simulation domain $V$ gives
rise to its corresponding probability density function (PDF), which can be
calculated as

$f_{Q}(q)=\frac{1}{V}\int_{V}\mathrm{d}^{3}x\,\delta\big{(}q-Q(\vec{x})\big{)}$
(1)

where $q$ is the corresponding random variable and $f_{Q}(q)\,\mathrm{d}q$
represents the probability of $Q$ occurring on the interval between $q$ and
$q+\mathrm{d}q$.

Two random quantities $Q(\vec{x}),U(\vec{x})$, in addition to their individual
statistics, also define the joint PDF, which can be found as follows

$f_{(Q,U)}(q,u)=\frac{1}{V}\int_{V}\mathrm{d}^{3}x\,\delta\big{(}q-Q(\vec{x})\big{)}\,\delta\big{(}u-U(\vec{x})\big{)}$
(2)

with $q,u$ being the corresponding random variables. This approach is
generalizable to an arbitrary number of random quantities.

Given a known PDF of $q$, $f_{Q}(q)$, and known relation between the two
correpsonding random quantities $R=\phi(Q)$, the PDF of $R$, $f_{R}(r)$, can
be obtained by the following simple random variable transformation

$f_{R}(r)=\int\mathrm{d}q\,f_{Q}(q)\,\delta\big{(}r-\phi(q)\big{)}$ (3)

which naturally preserves the probability fraction,
${f_{R}(r)\,\mathrm{d}r=f_{Q}(q)\,\mathrm{d}q}$.

A known multivariate PDF of $q,u,\dots$ can be used to determine PDF of any
transformed random variables depending on $q,u,\dots$ Specifically, PDFs for
the cases of $2\to 1$ variables, with the relation between the quantities
given by $R=\phi(U,Q)$, and $2\to 2$ variables, given by $R=\phi_{1}(Q,U)$,
$V=\phi_{2}(Q,U)$, can be determined using the following formulas

$\displaystyle f_{R}(r)$
$\displaystyle=\int\mathrm{d}q\int\mathrm{d}u\,f_{(Q,U)}(q,u)\,\delta\big{(}r-\phi(q,u)\big{)}$
(4) $\displaystyle f_{(R,V)}(r,v)$
$\displaystyle=\int\mathrm{d}q\int\mathrm{d}u\,f_{(Q,U)}(q,u)\,\times$
$\displaystyle\hskip
65.0pt\delta\big{(}r-\phi_{1}(q,u)\big{)}\,\delta\big{(}v-\phi_{2}(q,u)\big{)}$
(5)

In this work, log density $s=\log\rho$ and speed $v$ will take the role of the
baseline random variables with known statistical properties. Their dependent
counterparts will be thermal energy $F(s)=se^{s}$ and kinetic energy
$K(s,v)=e^{s}v^{2}/2$. Equation (3) proves useful in determining the PDF of
thermal energy in Sec. 4, transformation (4) is used to obtain the PDF of
kinetic energy in Sec. 5 and the joint PDF of thermal and kinetic energy is
obtained in Sec. 6 using (5).

##  3 Isothermal turbulence: overview

The dynamics of isothermal compressible gas in the inertial subrange, where
both driving and dissipation can be neglected, is governed by the equation of
continuity and Euler’s equation

$\displaystyle\frac{\partial\rho}{\partial
t}+\vec{\nabla}\cdot\left(\rho\vec{v}\,\right)$ $\displaystyle=0$ (6)
$\displaystyle\rho\left(\frac{\partial\vec{v}}{\partial
t}+\vec{v}\cdot\vec{\nabla}\vec{v}\right)$ $\displaystyle=-\vec{\nabla}p$ (7)

where $\rho$ is density and $\vec{v}$ is velocity. The system of governing
equations is closed by the isothermal equation of state, in which pressure is
directly proportional to density

$p=c_{s}^{2}\rho$ (8)

where $c_{s}$ is the speed of sound, conventionally set to 1 in all subsequent
expressions.

We will find it convenient to work with the log of density,
$s=\log\rho/\rho_{0}$, where $\rho_{0}$ is the mean density.

Assuming given initial conditions at $t=0$ and boundary conditions at all
later times the equations determine the fields evolving over time and space.
Periodic boundary conditions in a box with a side of length $L$ are used.
Driving terms at largest scales and dissipation on the smallest scale given by
the grid resolution balance out.

###  3.1 Single-variable statistics

In this section we summarize the statistical properties of the basic dynamical
variables, density $\rho$ and speed $v$, as well as develop mathematical tools
that will be useful later.

####  3.1.1 Density Distribution: Finite Shock Model

We consider a cloud of supersonic turbulence as an ensemble of shocks. Each
shock adjusts the pre-shock density by a multiplicative factor $m^{2}$ where
$m=v/c_{s}$ is the local Mach number, drawn from the the velocity
distribution, $f_{v}(v)$. The perfect log-normal behavior emerges as a
consequence of the central limit theorem if a parcel of gas, on average,
experiences a large amount of shocks $n$, each adjusting its density by a
random multiplicative factor. In Rabatin & Collins (2023b), we improve this
model by assuming that $n$ is finite and perhaps small.

We model $s$ as the sum of $n$ independent, identical events ${Y\sim\log
m^{2}}$, normalized to zero mean and unit variance, as

$s=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}Y_{i},$ (9)

which gives a charachteristic function of

$\displaystyle\phi(\omega;n)\equiv\left\langle e^{i\omega
s}\right\rangle=\left\langle e^{i\omega Y/\sqrt{n}}\right\rangle^{n}.$ (10)

The distribution is then the inverse Fourier transform of the characteristic
function,

$f_{\text{sh.}}(s;n)=\int\limits_{-\infty}^{\infty}\frac{\mathrm{d}\omega}{2\pi}e^{-i\omega
s}\phi(\omega;n).$ (11)

The behavior can be captured to a great degree of accuracy by measuring two
parameters from the ensemble; mean $\mu\equiv\langle s\rangle$ and variance
$\sigma^{2}\equiv\langle s^{2}\rangle-\mu^{2}$. Values $\mu$, $\sigma$, $n$
and $\langle e^{s}\rangle$ are related as follows

$\langle e^{s}\rangle=e^{\mu}\phi(-i\sigma;n)$ (12)

Since $\langle e^{s}\rangle=\langle\rho/\rho_{0}\rangle=1$, mean $\mu$ and
standard deviation $\sigma$ of log density uniquely determine the value for
$n$ by way of (12). All other expected values involving density can also be
calculated using the characteristic function, further details can be found in
(Rabatin & Collins, 2023b).

####  3.1.2 Velocity Distribution: Stretched Maxwell

Due to the symmetries of the turbulent motion, the statistics of speed can be
approximated by the Maxwellian distribution

$\mathcal{M}(v;\mathcal{M}_{\rm 1D})=\frac{4\pi v^{2}}{(2\pi\mathcal{M}_{\rm
1D}^{2})^{3/2}}\exp\left(-\frac{v^{2}}{2\mathcal{M}_{\rm 1D}^{2}}\right),$
(13)

where $\mathcal{M}_{\rm 1D}=\sqrt{\langle v^{2}\rangle/3}$ is the 1D r.m.s.
Mach number. As described in (Rabatin & Collins, 2023a), the statistics of
speed in our simulations can be described much more accurately by introducing
a correction term of a higher power in $v$ within the exponential,

$\mathcal{M}_{t}(v;\mathcal{M}_{\rm 1D},b)\propto
v^{2}\exp\left[-\frac{v^{2}}{2a^{2}}\left(1-b+\frac{bv^{2}}{a^{2}}\right)\right]$
(14)

where the strength of the correction is adjusted by the sloping parameter
$b\in[0,1]$ and $a$ is tuned so that the r.m.s. Mach number matches the
desired value, $3\mathcal{M}_{\rm 1D}^{2}=\langle v^{2}\rangle$. Further
details can be found in (Rabatin & Collins, 2023a).

####  3.1.3 Velocity Distribution: Finite Shock Model

We can also approximate the velocity distribution with FSM. The density
fluctuations are generated by random multiplicative factors drawn from the
velocity distribution. Similarly, the statistics of $v$ can be described in
terms of FSM itself, with exactly one shock. Specifically, the PDF of the
random variable $w\equiv\log v^{2}$ is given by

$\displaystyle f_{w}(w;\mathcal{M}_{\rm 1D},b)$
$\displaystyle=\mathcal{M}_{t}^{(\text{log})}(w;\mathcal{M}_{\rm 1D},b)$
$\displaystyle\equiv e^{w/2}\mathcal{M}_{t}(e^{w/2};\mathcal{M}_{\rm 1D},b)$
$\displaystyle\approx f_{\text{sh.}}(w;\mu_{w},\sigma_{w},n_{w})$ (15)

where $\mu_{w},\sigma_{w},n_{w}=1$ are the mean, standard deviation and the
effective number of shocks of $w$. The mean and variance can be calculated
analytically for $b=0$

$\displaystyle\mu_{w}(\mathcal{M}_{\rm 1D},b=0)$
$\displaystyle=2\log\mathcal{M}_{\rm 1D}+\Upsilon_{0}$ (16)
$\displaystyle\sigma_{w}(b=0)$ $\displaystyle=\Sigma_{0}$ (17)

with $\Upsilon_{0}=2-\gamma-\log 2\approx 0.730$ and
$\Sigma_{0}=\sqrt{\pi^{2}/2-4}\approx 0.967$ are numerical parameters related
to the ideal Maxwellian distribution (here $\gamma$ refers to the Euler
constant). Note, that while the mean of $\log v^{2}$ depends on the Mach
number, the variance is scale free.

The stretched Maxwellian can also be approximated by the Finite Shock Model,
with modified parameters reflecting the change in mean and variance of $\log
v^{2}$. The most extreme case when $b=1$ can be again obtained analytically,

$\displaystyle\mu_{w}(\mathcal{M}_{\rm 1D},b=1)$
$\displaystyle=2\log\mathcal{M}_{\rm 1D}+\Upsilon_{1}$ (18)
$\displaystyle\sigma_{w}(b=1)$ $\displaystyle=\Sigma_{1}$ (19) $\displaystyle
n_{w}$ $\displaystyle\approx 1$ (20)

where $\Upsilon_{1}=\pi/4-\gamma/2+(1/2)\log
2+\log\big{[}\Gamma(7/4)/\Gamma(5/4)\big{]}\approx 0.857$,
$\Sigma_{1}=\sqrt{\pi^{2}/4-2C}\approx 0.797$, with $C$ being the Catalan
constant.

The approximation (15) can be used across the full range of parameters $b$ for
the following interpolation between the cases for $b=0$ and $b=1$

$\displaystyle\mu_{w}(\mathcal{M}_{\rm 1D},b)$ $\displaystyle\approx
2\log\mathcal{M}_{\rm 1D}+0.730+\frac{2.170b}{1+16.01b^{0.82}}$ (21)
$\displaystyle\sigma_{w}(b)$ $\displaystyle\approx
0.967-\frac{2.349b}{1+12.84b^{0.78}}$ (22) $\displaystyle n_{w}$
$\displaystyle=1$ (23)

####  3.1.4 Statistical weighting

We will find it useful to compute distribution with different weighting. Any
PDF calculated from an ensemble via eq. 1 can be weighted by a non-negative
quantity $W$,

$f^{(W)}_{Q}(q)=\frac{1}{W_{\text{net}}}\int_{V}d^{3}x\,W(\vec{x})\,\delta(q-Q(\vec{x})),$
(24)

where $W_{\text{net}}$ is the total of $W$ within the ensemble. Alternative
weighting of the same quantity can be useful in certain circumstances. In the
context of this work, weighting by volume ($V$), mass ($M$) and kinetic energy
($K$) is utilized.

In addition to the log-density PDF, by default weighted by volume, the PDF
weighted by mass and kinetic energy is required. In (Rabatin & Collins, 2023a)
we find, that the mass- and kinetic energy-weighted instances can be related
to the volume-weighted PDF, ${f^{(V)}_{s}(s)=f_{\text{sh.}}(s;\mu,\sigma,n)}$,
in the following way

$\displaystyle f^{(M)}_{s}(s)$
$\displaystyle=e^{s}f^{(V)}_{s}(s)=e^{s}f_{\text{sh.}}(s;\mu,\sigma,n)$ (25)
$\displaystyle f^{(K)}_{s}(s)$
$\displaystyle=f^{(M)}_{s}(s+\log\mathfrak{X})=\mathfrak{X}\,e^{s}f_{\text{sh.}}(s;\mu-\log\mathfrak{X},\sigma,n)$
(26)

where $\mathfrak{X}=(\mathcal{M}_{\rm 1D}/\mathcal{M}_{M})^{2}$ is the square
of the ratio between the 1D r.m.s. Mach number and its mass-weighted
counterpart.

The mass-weighted PDF of speed becomes relevant to build the joint PDF of
density and speed. It cannot be related to its volume-weighted counterpart,
but can still be simply described by the stretched Maxwellian using its own
set of parameters

$f^{(M)}_{v}(v)=\mathcal{M}_{t}(v;\mathcal{M}_{M},b_{M})$ (27)

####  3.1.5 Convolutions of Finite Shock Models

Two FSM distributions can be convolved together. This can be interpreted as a
random variable $s=s_{1}+s_{2}$ being the sum of two independent random
variables. The result of such convolution can be approximated by another FSM
formula

$f_{\text{conv}}(s;\mu_{1},\mu_{2},\sigma_{1},\sigma_{2},n_{1},n_{2})=\\\
\int\limits_{-\infty}^{\infty}\mathrm{d}s^{\prime}\,f_{\text{sh.}}(s;\mu_{1},\sigma_{1},n_{1})f_{\text{sh.}}(s-s^{\prime};\mu_{2},\sigma_{2},n_{2})\approx\\\
f_{\text{sh.}}(s;\mu_{\text{conv}},\sigma_{\text{conv}},n_{\text{conv}})$ (28)

with parameters $\mu_{\text{conv}}(\mu_{1},\mu_{2})=\mu_{1}+\mu_{2}$,
$\sigma^{2}_{\text{conv}}(\sigma_{1},\sigma_{2})=\sigma_{1}^{2}+\sigma_{2}^{2}$
and $n_{\text{conv}}(\sigma_{1},\sigma_{2},n_{1},n_{2})$, that can be
determined from further physical considerations, for example matching the mean
of $e^{s}$.

We will find the convolution of density and speed useful momentarily, which we
now see simplifies to another $f_{\text{sh.}}$. Using both (15) and (28), we
find

$\displaystyle\text{C}(x;M,b,\mu,\sigma)$
$\displaystyle=f_{\text{sh.}}(\mu,\sigma,n)*\mathcal{M}_{t}^{(\log)}(M,b)$
$\displaystyle\approx
f_{\text{sh.}}(\mu,\sigma,n)*f_{\text{sh.}}(w;\mu_{w},\sigma_{w},n_{w})$
$\displaystyle\approx
f_{\text{sh.}}(x;\mu_{\text{conv}},\sigma_{\text{conv}},n_{\text{conv}})$ (29)

where

$\displaystyle\mu_{\text{conv}}(\mathcal{M}_{\rm 1D},b)$
$\displaystyle=\mu+\mu_{w}(\mathcal{M}_{\rm 1D},b)-\log 2$ (30)
$\displaystyle\sigma_{\text{conv}}(b)$
$\displaystyle=\sqrt{\sigma^{2}+\sigma^{2}_{w}(b)}$ (31)

and $n_{\text{conv}}\in[1,\infty)$ is some yet undetermined effective
convoluted shock parameter.

###  3.2 Joint statistics of density and speed

In what follows we revisit the statistical considerations of an ensemble
emerging from snapshots of our isothermal, turbulent simulations. The primary
dynamical quantities, density $\rho$ and velocity $\vec{v}$, can be treated as
random variables with their own distributions. The statistical properties of
each isothermal, turbulent ensemble are fully captured in the joint PDF of log
density $s=\log\rho/\rho_{0}$ and magnitude of velocity, speed $v=|\vec{v}|$.
For the purposes of this work, we use the detailed model proposed in Rabatin &
Collins (2023a), that accurately captures the statistics along with the
correlations between the two quantities. The model gives an explicit formula
for the joint PDF,

$f_{(s,v)}^{(V)}(s,v)\approx f_{s}^{(V)}(s)f_{v}^{(M)}(v)\,+\\\
+\,(\mathfrak{X}-1)^{-1}\Big{(}e^{-s}f_{s}^{(K)}(s)-f_{s}^{(V)}(s)\Big{)}\Big{(}f_{v}^{(V)}(v)-f_{v}^{(M)}(v)\Big{)},$
(32)

given 6 parameters directly measured from a simulated ensemble: 1D r.m.s. Mach
number $\mathcal{M}_{\rm 1D}$, mass-weighted Mach number
$\mathcal{M}_{M}=\sqrt{\langle\rho v^{2}\langle/3}$, the ratio between
$\mathcal{M}_{\rm 1D}^{2}$ and its mass-weighted counterpart,
$\mathfrak{X}=(\mathcal{M}_{\rm 1D}/\mathcal{M}_{M})^{2}$, average speed
$u=\langle v\rangle$, average speed weighted by mass $u_{M}=\langle\rho
v\rangle$, mean of log density $\mu$ and standard deviation of density,
$\sigma$. These parameters, presented in Table 1, in combination with the
model presented in Rabatin & Collins (2023a) allow us to calculate any
statistic involving density and speed.

##  4 Thermal energy

Thermal energy, $F$, is derived from the equations of motion in Section 4.1 as
well as from thermodynamic considerations in Section 4.2 General properties of
the thermal energy are discussed in Section 4.3 and the PDF of thermal energy
is discussed in Section 4.4.

###  4.1 Conservation of energy

In this section, we manipulate equations (6, 7) in line with Banerjee &
Kritsuk (2018). The goal is to obtain a conserved scalar quantity, total
energy, along with its flux. We start by taking the dot product of $\vec{v}$
and Eq. (7), and subsequently using Eq. (6). These steps yield the following
equality

$\frac{\partial}{\partial t}\left(\frac{1}{2}\rho
v^{2}\right)+\vec{\nabla}\cdot\left(\frac{1}{2}\rho
v^{2}\vec{v}\,\right)=-\vec{v}\cdot\vec{\nabla}p$ (33)

Replacing $p$ with $\rho$ on the right-hand side and using the chain rule
gives us

$-\vec{v}\cdot\vec{\nabla}p=-\vec{\nabla}\cdot\left(\rho\vec{v}\left(1+\log\rho/\rho_{0}\right)\right)-\frac{\partial}{\partial
t}\left(\rho\log\rho/\rho_{0}\right)$ (34)

where $\rho_{0}$ is an arbitrary density scale, for example the ambient
density, or, in our case, the average density, that is conventionally set to 1
in our simulations.

Putting everything together, we get the conservation of energy

$\frac{\partial E_{\rm net}}{\partial
t}+\vec{\nabla}\cdot\vec{j}_{E,\text{net}}=0$ (35)

where $E_{\rm net}=K+F$ is the total energy density per unit volume consisting
of the sum of kinetic and thermal energy density per unit volume,

$\displaystyle K$ $\displaystyle=\frac{1}{2}\rho v^{2}$ (36) $\displaystyle F$
$\displaystyle=\rho\log\rho,$ (37)

and $\vec{j}_{E,\text{net}}$ is the energy density flux per unit area and
time,

$\vec{j}_{E,\text{net}}=\left[\frac{v^{2}}{2}+\left(1+\log\rho\right)\right]\rho\vec{v}.$
(38)

The conservation of energy as presented here holds in the inertial subrange of
the length scales of the turbulent motion, in which the energy cascades into
the smaller scales and both the driving and dissipation can be neglected.

While the conservation of energy provides a useful insight into the various
types of energy, in what follows we find the context for the thermal energy
$F$ within the framework of the thermodynamical laws.

###  4.2 Thermodynamical interpretation of $F$

To get a handle on the thermal energy found in the previous section, we start
from the first law of thermodynamics that holds for the LTE of an arbitrary,
small volume $\mathrm{d}V$ within the turbulent fluid,

$\mathrm{d}U=T\mathrm{d}S-p\mathrm{d}V$ (39)

where $U$ and $S$ are internal energy and entropy contained within the volume,
respectively.

The first law of thermodynamics can be re-expressed in terms of quantities
defined per unit mass; $U\to u$ (internal energy per unit mass), $S\to s$
(entropy per unit mass) and $V\to 1/\rho$ (volume per unit mass).
Additionally, in the isothermal framework, we can replace $p=\rho$. The first
law of thermodynamics then simplifies as follows

$\mathrm{d}u=T\mathrm{d}s+\frac{1}{\rho}\mathrm{d}\rho$ (40)

The specific internal energy is a function of entropy and density. However, to
utilize the fact, that the temperature is kept constant, a Legendre
transformation introducing a thermodynamical potential involving temperature
and density will prove useful. The Legendre transformation applies for the
specific quantities in the standard way,

$\mathcal{F}=u-Ts$ (41)

where $\mathcal{F}$ is the specific Helmholtz free energy (free energy per
unit mass). Replacing $u$ in Eq. (40) with $\mathcal{F}+Ts$ results in

$\mathrm{d}\mathcal{F}=-s\mathrm{d}T+\frac{1}{\rho}\mathrm{d}\rho$ (42)

which shows, that $\mathcal{F}$ is a function of $T,\rho$.

This equation allows us to directly use the isothermal assumption,
$T=\text{const.}$, and therefore $\mathrm{d}T=0$. Integrating over $\rho$
results in

$\mathcal{F}=\log\rho/\rho_{0}$ (43)

where $\rho_{0}$ again serves as an arbitrary constant with units of density
and is conventionally set to 1. Going back to the energy per unit volume,
$F=\rho\mathcal{F}$, gives us

$F=\rho\log\rho$ (44)

This derivation formalizes $F$ within the thermodynamical framework as the
Helmholtz energy density per unit volume.

Figure 1:  Function $F(s)=se^{s}$ together with it’s inverse,
${s(F)=W_{-1,0}(F)}$. The branch -1 is denoted with a dashed line, while the
branch 0 is plotted with solid line. Figure 2: Plots of the PDF (vertical
axis) of the shifted log of thermal energy, $\varphi=\log(1+eF)$ (horizontal
axis). Dots along with horizontal bars represent data points and error bars.
Solid curve depicts the exact formula (48). The PDF has been rescaled
vertically by the value $\displaystyle
f^{(+)}_{\text{max}}=\lim_{\varepsilon\to 0^{+}}f_{\varphi}(\varepsilon)$ for
the sake of clarity.

###  4.3 General properties, geometrical considerations

In this section we explore the properties of the Helmholtz free energy,
$F=\rho\log\rho=se^{s}$ more closely.

The dimensionless thermal energy $F=se^{s}$ is a non-injective function of
log-density $s$ that attains its minimum $F_{\text{min}}=-1/e$ at $s=-1$. For
$0<s<1$, $F<0$ and is multivalued, as ${F(\rho=0)=F(\rho=1)=0}$. The many-to-
one property can be observed from Fig. 1. Since the inverse function to
$se^{s}$ has two branches for $F<0$, an additional contribution to the PDF of
$F$ is expected for negative thermal energies. This additional weight is
expected to become important should the PDF of log-density, $f_{s}(s)$, be
broad enough to to reach both positive and negative tail in $s$. Therefore,
the weight of the negative thermal energy is expected to be more prominent in
high Mach numbers and/or compressive driving, as the deviation in $s$,
$\sigma$, is large in both cases. Additionally, the presence of a minimum at
$F=-1/e$ foreshadows a power law spike in the thermal energy PDF near the left
edge of the domain. This feature will be used in 4.4 to provide an analogous
random variable in a logarithmic scale, highlighting the power law behavior,
$\varphi=\log(1+eF)$. Additionally, this choice preserves the position of the
singularity at $\varphi=F=0$.

The mean value of $F$ can be calculated using the FSM (Rabatin & Collins,
2023b). Moreover, this mean value can be re-interpreted as the mass-weighted
average of log density,

$\langle F\rangle=\langle se^{s}\rangle=\langle s\rangle_{M}=\mu_{M}$ (45)

For ideally lognormal statistics of density, the mass-weighted mean of log-
density is simply related to its volume-weighted counterpart,
$\mu_{M}=-\mu_{V}=\sigma^{2}/2>0$. The positive sign of $\mu_{M}$ does not
change for a finite value of $n$, therefore, the mean thermal energy is always
positive $\langle F\rangle>0$.

###  4.4 Statistical properties

The Helmholtz free energy of the form $\rho\log\rho$, or, in terms of variable
$s$, $se^{s}$, is only a function of one of the primary variables - density.
As such, its PDF can be obtained solely using the PDF of density by a random
variable transformation and is given by equation (3) where $s,F$ take the role
of the baseline and transformed variables, respectively.

To explicitly integrate over the $\delta$-function, the roots of equation
$F=se^{s}$ must be found for all values of $F\in[-1/e,\infty)$. These roots
are captured in a special function, Lambert product logarithm $W$. Since the
function $se^{s}$ is not monotonic, its inverse, as depicted in Figure 1, is
multivalued. For the purposes of the thermal energy PDF, we recognize two real
branches of the product log, the principal branch $W_{0}$ for arguments
$\in[-1/e,\infty)$ and values ranging from $-1$ to $\infty$, and an additional
branch denoted by $W_{-1}$ for arguments $\in[-1/e,0)$ and values
$\in(-\infty,-1]$.

With these definitions, the $\delta$-function appearing in Equation (3) can be
rewritten as

$\delta(F-se^{s})=\sum_{k=-1,0}\frac{\delta(s-W_{k}(F))\,\theta_{k}(F)}{e^{s}\left|1+s\right|}$
(46)

where the sum is over the two real branches of the product log and each branch
only applies to its appropriate range as described earlier,

$\displaystyle\theta_{-1}(x)$ $\displaystyle=\begin{cases}1;\;-1/e\leq x<0\\\
0;\;\text{otherwise}\end{cases}$ $\displaystyle\theta_{0}(x)$
$\displaystyle=\begin{cases}1;\;x\geq 0\\\ 0;\;\text{otherwise}\end{cases}$

The PDF of $F$ is

$f_{F}(F)=\sum_{k=-1,0}\frac{f_{s}\Big{(}W_{k}(F)\Big{)}\,\theta_{k}(F)}{\exp\Big{(}W_{k}(F)\Big{)}\left|1+W_{k}(F)\right|}$
(47)

Since the PDF of thermal energy diverges when $F\to-1/e$ from above, we find
it advantageous to work with the PDF of its logarithmic counterpart,
$\varphi=\log(1+eF)$. The PDFs of $\varphi$ and $F$ are related by a simple
variable transformation

$f_{\varphi}(\varphi)=e^{\varphi-1}f_{F}(e^{\varphi-1}-e^{-1})$ (48)

Fig. 2 shows the histograms of log thermal energy density obtained from the
numerical simulations as solid points with error bars. Our model is depicted
as a solid black curve, showing a remarkable match. The power law behavior at
the left side, $F=-1/e$, $f_{F}\sim(F+1/e)^{\alpha}$ is apparent from the
linear tail as seen in plots. One curious feature is the discontinuity around
$\varphi=F=0$. This feature is more pronounced for the compressively driven
simulations, and for all driving modes the discontinuity increases with Mach
number. As expected, more high energy gas is created with increasing Mach
number, and this effect is more pronounced with compressive driving. Our model
perfectly captures these features.

##  5 Kinetic energy

Figure 3: Plots of the PDF (vertical axis) of logarithm of kinetic energy
density, $\kappa$ (horizontal axis). Dots and vertical lines represent data
points with error bars. Solid curve depicts the exact formula using Eqs. (51,
52), dashed line represents the approximation (53). For the sake of clarity,
the horizontal axis is shifted and rescaled by the corresponding mean and
standard deviation, while the vertical axis is rescaled by the maximum value.

In Section 4.1 we showed, that the kinetic energy is $K=\rho v^{2}/2$. Since
it depends on both density and speed, statistics of both, as well as the
correlations between the two, become important.

Some general properties of the kinetic energy are discussed in Section 5.1.
The exact formula for the PDF is discussed in Section 5.2, an approximation is
given in 5.3.

###  5.1 General properties

The average kinetic energy density $\langle K\rangle$ is directly related to
one of the parameters used to construct the joint statistics of density and
speed,

$\langle K\rangle=\frac{1}{2}\langle\rho v^{2}\rangle=\frac{1}{2}\langle
v^{2}\rangle_{M}=\frac{3}{2}\mathcal{M}_{M}^{2}$ (49)

where $\mathcal{M}_{M}$ is the mass-weighted Mach number. For specific values,
see Table 1. It should be noted, that since density and speed are correlated,
$\mathcal{M}_{\rm 1D}\neq\mathcal{M}_{M}$ and therefore, the value for
$\langle K\rangle$ differs from $(3/2)\mathcal{M}_{\rm 1D}^{2}$. This
difference is especially prominent in the transsonic regime.

Due to the additive nature of logarithms, we will later introduce
logarithmically transformed kinetic energy, $\kappa=\log K=s+2\log v-\log 2$.
This choice comes natural, as the result is a sum of near-Gaussian (albeit
correlated) random variables, $s$ and $\log v$. Log-density $s$ is near-normal
depending on the number of shocks, $n$, and $\log v$ is distributed according
to the FSM with effective $n=1$.

###  5.2 Kinetic PDF: exact formula

Following Eq. (4), the PDF of kinetic energy density, $K(s,v)=e^{s}v^{2}/2$,
can be formally written as the following transformation of random variables,
given a known joint PDF of the baseline variables $s,v$

$f_{K}(K)=\int\limits_{-\infty}^{\infty}\mathrm{d}s\int\limits_{0}^{\infty}\mathrm{d}v\,f_{(s,v)}(s,v)\,\delta\left(K-\frac{1}{2}e^{s}v^{2}\right)$
(50)

which can further be simplified by integrating over $v$ using the
$\delta$-function

$f_{K}(K)=\int\limits_{-\infty}^{\infty}\mathrm{d}s\,\frac{f_{(s,v)}\Big{(}s,\sqrt{2e^{-s}K}\Big{)}}{\sqrt{2e^{s}K}}$
(51)

This exact formula cannot be simplified any further without considering a
specific model for the joint statistics for density and speed, $f_{(s,v)}$.
Even for the simplest case of the lognormal PDF of density and Maxwellian PDF
of speed, the integral cannot be performed in a closed form.

Due to the lognormal nature of density, we examine $\kappa=\log K$, for which
the PDF is a simple variable transformation of $f_{K}$,

$f_{\kappa}(\kappa)=e^{\kappa}f_{K}(e^{\kappa})$ (52)

where $f_{K}$ is the baseline PDF of $K=\frac{1}{2}\rho v^{2}$.

###  5.3 Logarithmic approximation

In this section we will derive an approximate analytic formula directly for
PDF of $\kappa=\log K$ given a few basic assumptions.

We begin by integrating over the joint distribution, (32), by way of (4),
where $R$ now represents $\kappa=\log K=s+w-\log 2$. Since the terms in
$f_{(s,v)}$ are uncorrelated, this then becomes the convolution of five pairs
of density and speed distributions, which by (29) becomes

$\displaystyle f_{\kappa}(\kappa)=\,$
$\displaystyle\text{C}(\kappa;\mathcal{M}_{M},b_{M},\mu,\sigma)\,+$ (53)
$\displaystyle(\mathfrak{X}-1)^{-1}\Big{(}\mathfrak{X}\,$
$\displaystyle\text{C}(\kappa;\mathcal{M}_{\rm
1D},b,\mu-\log\mathfrak{X},\sigma)\,-$ $\displaystyle\mathfrak{X}\,$
$\displaystyle\text{C}(\kappa;\mathcal{M}_{M},b_{M},\mu-\log\mathfrak{X},\sigma)\,-$
$\displaystyle\text{C}(\kappa;\mathcal{M}_{\rm 1D},b,\mu,\sigma)\,+$
$\displaystyle\text{C}(\kappa;\mathcal{M}_{M},b_{M},\mu,\sigma)\Big{)}$

Since for the general shape of the FSM formula, $\mu,\sigma,n,\langle
e^{s}\rangle$ are tied together via eq. (12), we can re-express the mean
kinetic energy $\langle K\rangle=\langle
e^{\kappa}\rangle=\frac{3}{2}\mathcal{M}_{M}^{2}$, or, after multiplying by 2,

$\langle 2e^{\kappa}\rangle=3\mathcal{M}_{M}^{2}=\\\
\exp\left(\mu_{\text{conv}}(\mathcal{M}_{M},b_{M})\right)\phi\left(-i\sigma_{\text{conv}}(b_{M}),n_{\text{conv}}\right)$
(54)

where the cancellation of last 4 terms happens under the assumption that the
dependence of $n_{\text{conv}}$ on $b$ can be neglected. Using relations (21,
22, 30, 31), constraint (54) can be expressed as

$3\mathcal{M}_{M}^{2}=\exp\left(\mu+\mu_{w}(\mathcal{M}_{M},b_{M})\right)\times\\\
\phi\left(-i\sqrt{\sigma^{2}+\sigma_{2}^{2}(b_{M})},n_{\text{conv}}\right)$
(55)

This equation can be used to find $n_{\text{conv}}$ given values of
$\mu,\sigma,\mathcal{M}_{M},b_{M}$.

Fig. 3 shows the log kinetic energy density histograms extracted from the
simulations as solid points with error bars. These are compared to our exact
model (Equation (51), solid line) as well as the approximate formula (Equation
(53), dashed line). It should be noted that our approximation holds even in
cases where the statistics of density strongly deviates from lognormal and
speed from ideal Maxwellian (see Rabatin & Collins (2023a).) It deviates from
the true distributions in datasets exhibiting larger correlation between
density and speed, such as the transsonic simulation with mixed or purely
solenoidal driving.

##  6 Joint PDF

Figure 4: The joint PDF of $\varphi$ and $\kappa$ with our detailed model.
Color and dashed contour shows data. Solid line shows contours of our model.
Each simulation is labeled with $\xi-M$ in the top left corner. The axes are
shifted and rescaled for the sake of clarity.

The joint PDF of thermal and kinetic energy can be obtained from a known joint
distribution of $s,v$ and a simple random variable transform $(s,v)\to(F,K)$
using the transformation (5).

$\displaystyle f_{(F,K)}(F,K)$ $\displaystyle=$
$\displaystyle\int\limits_{-\infty}^{\infty}\mathrm{d}s\int\limits_{0}^{\infty}$
$\displaystyle\mathrm{d}v\,f_{(s,v)}(s,v)\,\delta\left(F-se^{s}\right)\delta\left(K-\frac{1}{2}e^{s}v^{2}\right)$
(56)

Similar to the derivation of the PDF of thermal energy, we solve for the two
roots $k=-1,0$ of $s$ given a value of $F$ using the product logarithm,
$s_{k}=W_{k}(F)$. The other $\delta$-function can be used to determine the
sole root for $v=\sqrt{2e^{-s}K}$. Together, we consider both
$\delta$-functions in their root form

$\delta\left(F-se^{s}\right)\delta\left(K-\frac{1}{2}e^{s}v^{2}\right)=\\\
\sum_{k=-1,0}\frac{\delta(s-W_{k}(F))\delta(v-\sqrt{2e^{-s}K})\theta_{k}(F)}{e^{2s}|1+s|v}$
(57)

and the joint energy PDF can be written analytically as follows

$f_{(F,K)}(F,K)=\\\
\sum_{k=-1,0}\frac{f_{(s,v)}\left(W_{k}(F),\sqrt{2\exp(-W_{k}(F))K}\right)\theta_{k}(F)}{\exp\left(\frac{3}{2}W_{k}(F)\right)\left|1+W_{k}(F)\right|\sqrt{2K}}$
(58)

Since density is distributed roughly lognormally and the thermal energy
exhibits a power-law behavior near its minimum at $F=-1/e$ for $s=-1$, we
consider the logarithm of both components of energy in the following form:
$\varphi=\log(1+eF),\kappa=\log(K)$, for which the joint PDF can be written as

$f_{(\varphi,\kappa)}(\varphi,\kappa)=e^{\varphi+\kappa-1}f_{(F,K)}\left(e^{\varphi-1}-e^{-1},e^{\kappa}\right)$
(59)

The joint statistics of $\varphi$ (horizontal axis) and $\kappa$ (vertical
axis) are shown in Fig. 4. The histograms extracted from the simulations are
displayed as dashed contours, while the color represents the fractional
probability. Our result from equation 58 transformed via 59 is shown via solid
contours. The underlying statistical model for the density-speed PDF uses the
detailed basis discussed in Rabatin & Collins (2023a). The bimodal behavior of
the correlation between energies is apparent from the figure and perfectly
reproduced by our model. Each simulation has two lobes, one with high values
of kinetic and thermal energies, and one with low values of kinetic and
thermal energies. Gas near the mean density rarely attains high kinetic
energy. The discontinuity at $F=\varphi=0$ appears again, with the curious
result that gas with high thermal energy must have high kinetic energy.

##  7 Conclusions

In this work we provided a deep dive into the statistics of energy in
isothermal supersonic turbulence. Using a baseline model of the joint PDF of
density-speed we were able to derive the marginalized PDF of thermal and
kinetic energy, as well as their joint PDF.

The PDF of thermal energy carries the signature of the density-dependent
thermal energy, $F=se^{s}$. This manifests in two ways; first, a power law
behavior emerges near the left edge of the thermal energy PDF, at $F=-1/e$, as
a direct consequence of said minimum. Second, a jump discontinuity in the PDF
occurs at $F=0$ due to the region $F<0$ of the thermal PDF carrying additional
weight originating in the non-injective nature of the thermal energy function
of log-density. Specifically, the region $s\in(-\infty,-1]$ gives rise to the
same values of thermal energy as $s\in[-1,0)$. These features will show
gradually with increasing width in log-density. Our model matches the
numerically simulated data exceptionally well.

It is advantageous to discuss the log of the kinetic energy. The PDF of
$\kappa=\log K=s+w-\log 2$ uses the additive nature of logarithms and,
provided $s$ and $w$ are independent, its PDF can be evaluated via convolution
of the underlying PDFs for $s$, $w=\log v^{2}$. This was used with the
combination with the underlying model for the joint statistics of $s,v$,
consisting of 5 marginalized terms, providing an approximate formula for the
PDF of $\kappa$. This formula approximates the exact formula, given as a
numerical integral, fairly well in all datasets except in those that have a
higher degree of correlation between density and speed. The exact formula
matches the data reasonably well in the same datasets, which indicates, that
the correction to the joint PDF of density and speed, carrying the
correlations between the two, is lacking in precision.

Finally, the joint PDF of energy was expressed in closed form and showed a
remarkable match with the numerically obtained histograms. The joint PDF
exhibits an interesting feature, more prominent in more supersonic
simulations, that pronounces the weight of for $F<0$. For highly supersonic
motion, regardless of the driving pattern, most of the gas, by volume,
exhibits negative thermal energy.

By exploring the statistics of energy in this simplified model of interstellar
turbulence, we hope to lay groundwork for further modelling of the complex
dynamics of the ISM.

## Acknowledgements

Support for this work was provided in part by the National Science Foundation
under Grant AST-1616026 and AST-2009870. Simulations were performed on
_Stampede2_ , part of the Extreme Science and Engineering Discovery
Environment (XSEDE; Towns et al., 2014), which is supported by National
Science Foundation grant number ACI-1548562, under XSEDE allocation TG-
AST140008.

## References

  * Banerjee & Kritsuk (2018) Banerjee, S., & Kritsuk, A. G. 2018, Phys. Rev. E, 97, 023107
  * Bryan et al. (2014) Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19
  * Federrath (2013) Federrath, C. 2013, MNRAS, 436, 1245
  * Federrath et al. (2021) Federrath, C., Klessen, R. S., Iapichino, L., & Beattie, J. R. 2021, Nature Astronomy, 5, 365
  * Federrath et al. (2008) Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79
  * Federrath et al. (2010) Federrath, C., Roman-Duval, J., Klessen, R. S., Schmidt, W., & Mac Low, M. M. 2010, A&A, 512, A81
  * Mac Low (1999) Mac Low, M.-M. 1999, ApJ, 524, 169
  * Mac Low & Klessen (2004) Mac Low, M.-M., & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125
  * Porter et al. (1999) Porter, D., Pouquet, A., Sytine, I., & Woodward, P. 1999, Physica A: Statistical Mechanics and its Applications, 263, 263, proceedings of the 20th IUPAP International Conference on Statistical Physics
  * Porter & Woodward (2000) Porter, D. H., & Woodward, P. R. 2000, The Astrophysical Journal Supplement Series, 127, 159
  * Rabatin & Collins (2023a) Rabatin, B., & Collins, D. C. 2023a, MNRAS, 525, 297
  * Rabatin & Collins (2023b) —. 2023b, MNRAS, 521, L64
  * Scalo & Elmegreen (2004) Scalo, J., & Elmegreen, B. G. 2004, ARA&A, 42, 275
  * Towns et al. (2014) Towns, J., Cockerill, T., Dahan, M., et al. 2014, Computing in Science and Engineering, 16, 62
  * Uhlenbeck & Ornstein (1930) Uhlenbeck, G. E., & Ornstein, L. S. 1930, Physical Review, 36, 823
  * Vazquez-Semadeni (1994) Vazquez-Semadeni, E. 1994, ApJ, 423, 681
  * Woodward & Colella (1984) Woodward, P., & Colella, P. 1984, Journal of Computational Physics, 54, 115


# Andreev bound states and supercurrent
in disordered spin–orbit-coupled nanowire SNS-junctions

Jonas Lidal  Center for Quantum Spintronics, Department of Physics, Norwegian
University of Science and Technology, NO-7491 Trondheim, Norway  Norwegian
National Security Authority, Norway  Jeroen Danon  Center for Quantum
Spintronics, Department of Physics, Norwegian University of Science and
Technology, NO-7491 Trondheim, Norway

###### Abstract

We use a single-channel scattering formalism to derive general expressions for
the Andreev bound-state energies and resulting current–phase relationship in a
one-dimensional semiconductor-based SNS-junction, including arbitrarily
oriented effective spin–orbit and Zeeman fields and taking into account
disorder in the junction in by including a single scatterer with transmission
probability $0\leq T^{2}\leq 1$, arbitrarily located in the normal region. We
first corroborate our results by comparing them to the known analytic
limiting-case expressions. Then we simplify our general result in several
additional limits, including the case of the scatterer being at a specific
location and the cases of small and large spin–orbit fields compared to the
Zeeman splitting, assuming low transparency ($T\ll 1$). We believe that our
results could be helpful for disentangling the main spin-mixing processes in
experiments on low-dimensional semiconductor-based SNS-junctions.

##  I Introduction

Hybrid superconductor–semiconductor devices have attracted much attention in
recent years, mostly due to the exotic spin physics that can emerge from the
interplay of spin–orbit interaction (SOI), a strong Zeeman effect, and
proximity-induced superconducting correlations [1, 2]. One potential
application of such devices lies in the engineering of topological
superconductivity [3, 4, 5, 6], which is expected to come with localized
Majorana-like states that could be useful for fault-tolerant topological
quantum computation [7, 8]. The simultaneous breaking of inversion and time-
reversal symmetry associated with the SOI and Zeeman splitting makes
superconductor–semiconductor devices also natural candidates for exploring the
so-called superconducting diode effect [9, 10], which is potentially useful in
the context of low-dissipation electronics. Furthermore, advancing this hybrid
platform could result in a seamless integration of semiconductor- and
superconductor-based technologies, for instance, enabling long-range qubit-
qubit coupling of spin qubits via superconducting resonators [11].

One of the simplest of such hybrid devices one can make is the
superconductor–normal–superconductor (SNS) junction, i.e., two superconducting
regions that are connected via a normal (semiconducting) region. A finite
phase difference between the two superconductors typically results in a
supercurrent through the junction, the so-called Josephson effect. In the
presence of SOI and a Zeeman splitting in the normal region, the Josephson
effect can become asymmetric [12], yielding, e.g., a $\varphi_{0}$-junction
(where the free energy is minimal at a finite phase difference and a finite
supercurrent can run at zero phase difference) or supercurrent rectification
(when the maximal supported supercurrent in the two directions of the junction
is different). This makes the semiconductor-based SNS-junction a versatile
device that can be used for realizing superconducting diodes [13, 14, 15, 16,
17, 18, 19, 20, 21], to construct “phase batteries” [22, 23], or to realize
tunable $0$–$\pi$ junctions that could find use in superconducting qubit
platforms [24]. Moreover, they can potentially be tuned to a topological
regime where the junction hosts Majorana bound states [25, 26], making them a
natural building block for two-dimensional arrays of Majorana states [27].
Apart from all these envisioned applications, another potential use of SNS-
junctions is as a probe to study the spin physics in the semiconductor,
possibly revealing details about the interplay of the SOI and the Zeeman
effect [28] or about the occurrence of topological phase transitions [29, 30,
31, 32].

Because of their interesting physics and the potential use of SNS-junctions,
there is a considerable amount of theoretical literature investigating the
(critical) current through such junctions [33]. In the absence of SOI and a
Zeeman splitting, results have been derived for both short [34, 35, 36] and
long [37, 38, 39] junctions (compared to the superconducting coherence
length). When including spin mixing in the normal region, the problem becomes
more complex [40, 41, 42, 43, 44, 45, 46, 47], and for higher-dimensional
junctions it is not always possible to arrive at analytic results [48, 49, 50,
51].

A particularly simple approach that often yields explicit results is to
consider a (quasi-)one-dimensional junction, where all transport takes place
through a single or few channels. In that case, the supercurrent can often be
calculated analytically from the bound-state spectrum in the junction, which
follows from solving a few-channel scattering problem. Both for long [52] and
short [53] junctions, expressions have been derived that include the effect of
SOI and a Zeeman splitting. Taking also disorder into account is relatively
straightforward in this setup: Both disorder at the SN-interfaces and inside
the normal region can be modeled by including scaterrers with finite
reflection probability, which resuled in analytic expressions for the
supercurrent in certain limiting cases [54, 55, 56, 57]. These insights in the
single-channel physics are not only crucial as building blocks for a deeper
understanding of the multi-channel (two- or three-dimensional) case, but they
are also experimentally relevant for hybrid devices based on one-dimensional
semiconducting nanowires [58].

In this paper, we develop the single-channel scattering approach of Refs. [55,
56] further, generalizing the analytic results to describe any desired
combination of effective spin–orbit and Zeeman fields in the junction,
including a single point-like scaterrer at an arbitrary location in the normal
region to account for disorder. We show how our expressions reduce to the
known results in all limiting cases, and we demonstrate how additional limits
can be investigated easily, yielding compact and insightful expressions for
the supercurrent and critical current through the junction. We believe that
our results present a valuable contribution to the steadily growing
understanding of hybrid superconductor–semiconductor devices and that they
could be a useful asset for disentangling spin–orbit and Zeeman physics in
experiments on low-dimensional SNS-junctions.

The rest of the paper is structured as follows. In Sec. II we introduce the
setup we consider and the model we use to describe the system. We then show
how all ingredients of the model can be combined into a straightforward self-
consistency equation for the bound-state energies in the long-junction limit.
In Sec. III we present a general analytic solution of this equation, which can
be used to find the (critical) current through the junction. We first compare
our general result to the known limiting-case expressions, and then both
investigate new limits in the low-transparency case and present simple closed-
form general expressions for fixed positions of the scatterer. Finally, we
present a short conclusion in Sec. IV. The full explicit expression of our
analytic solution of the scattering problem is presented in an Appendix.

##  II Model

Setup—We consider a one-dimensional single-channel semiconductor nanowire
connected to two superconducting contacts that are separated by a distance
$L$, as sketched in Fig. 1(a). Due to the proximity effect, the two ends of
the wire will thus acquire a finite superconducting pairing potential, which
we assume to be of equal strength on the two sides. The wire can have
significant spin–orbit coupling and we also include an external magnetic field
that gives rise to a Zeeman splitting inside the wire. We model the disorder
in the junction in the simplest way possible, i.e., by including a single
scatterer, which is characterized by a transmission amplitude $T$ and is
located at the position $x=\alpha L$ with $0\leq\alpha\leq 1$.

Figure 1: (a) Sketch of the setup considered. A one-dimensional semiconductor
nanowire (gray) is connected to two superconductors (blue), separated by a
distance $L$, which have a phase difference $\phi=\phi_{l}-\phi_{r}$. Disorder
in the wire is modeled by a single scatterer with transmission probability
$T$, located at position $x=\alpha L$. (b) The Zeeman and effective spin–orbit
fields (black) are not necessarily parallel. The total fields the electrons
traveling in the two directions experience are their sum and difference fields
(green).

We describe the electrons in the wire using a Bogoliubov–de Gennes
Hamiltonian, $H=\frac{1}{2}\int dx\,\Psi^{\dagger}{\cal H}\Psi$, where
$\Psi=[\psi_{\uparrow}(x),\psi_{\downarrow}(x),\psi_{\downarrow}^{\dagger}(x),-\psi_{\uparrow}^{\dagger}(x)]^{T}$
with $\psi_{\uparrow(\downarrow)}(x)$ being the annihilation operator for an
electron with spin up(down) at position $x$. We further use

${\cal H}=\left(\begin{array}[]{c c}{{\hat{\cal H}_{0}}}&{{\hat{\Delta}}}\\\
{{{\hat{\Delta}}^{\dagger}}}&{{-\mathcal{T}\hat{\cal
H}_{0}\mathcal{T}^{\dagger}}}\end{array}\right),$ (1)

where $\hat{\cal H}_{0}$ and $\hat{\Delta}$ are $2\times 2$ matrices in spin
space and $\mathcal{T}$ is the time-reversal operator. The Hamiltonian
$\hat{\cal H}_{0}$ reads formally as

$\hat{\cal
H}_{0}=-\frac{\hbar^{2}\partial_{x}^{2}}{2m}-\mu+V_{0}\delta\left(x-\alpha
L\right)+({\bf{B}}_{\text{Z}}-i\partial_{x}\boldsymbol{\alpha}_{\rm
so})\cdot{\boldsymbol{\sigma}},$ (2)

where $V_{0}$ determines the strength of the scattering potential, ${\bf
B}_{\rm Z}$ is the Zeeman field, and the vector $\boldsymbol{\alpha}_{\rm so}$
characterizes the spin–orbit coupling. The vector
$\boldsymbol{\sigma}=\\{\sigma_{x},\sigma_{y},\sigma_{z}\\}$ consists of the
three Pauli matrices acting in spin space. In our Nambu basis, the off-
diagonal block in ${\cal H}$ becomes $\hat{\Delta}=\Delta(x)\mathbb{1}$, where
we assume the pairing potential to be

$\Delta(x)=\begin{cases}\Delta_{0}e^{i\phi_{l}},\quad&\text{for }x<0,\\\
\Delta_{0}e^{i\phi_{r}},\quad&\text{for }x>L,\\\ 0,\quad&\text{for }0\leq
x\leq L,\end{cases}$ (3)

which describes the induced superconductivity in the two ends of the wire,
$\phi_{l(r)}$ being the phase of the left(right) superconducting contact.

The goal is to find the supercurrent through the junction $I(\phi)$ as a
function of the phase difference $\phi=\phi_{l}-\phi_{r}$ between the two
superconductors. Assuming the short-junction limit, i.e. the length of the
junction being much smaller than the coherence length in the normal part,
$L\ll\xi_{N}\sim\hbar v_{\rm F}/\Delta_{0}$, with $v_{\rm F}$ the Fermi
velocity, and for simplicity focusing on the limit of zero temperature, the
supercurrent follows as [35]

$I(\phi)={\frac{2e}{\hbar}}{\sum_{n}}^{\prime}{\frac{dE_{n}}{d\phi}},$ (4)

where the $E_{n}$ are the energies of the subgap Andreev bound states in the
junction and the prime in the sum indicates that we include negative Andreev
levels only. This expression gives the supercurrent in the ground-state
configuration of the junction, thus assuming that it had enough time to
completely thermalize. When the phase difference $\phi$ varies over time, then
the supercurrent can be different depending on whether the variation is slow
enough to keep the system (approximately) in its ground state continuously or
not [55]; this also leads to different maximal (critical) supercurrents that
can be supported in thermalized and non-thermalized junctions. Here, we will
assume for simplicity that the system is always in its ground state, in which
case Eq. (4) gives the correct current–phase relationship.

We calculate the subgap Andreev levels, $E_{n}$, using a quantization
argument, where we describe the propagation of subgap electrons and holes in
the system in terms of a scattering matrix. We split our scattering
description into the following parts, as illustrated in Fig. 2(a): (i) free
propagation of electrons and holes on both sides of the scatterer, pictured
with the straight solid and dashed arrows, respectively (the small black
vertical arrows label the two different spin species), (ii) Andreev reflection
at the SN-interfaces on the left and right ends of the wire, illustrated by
the curves in the grey regions connecting incoming electrons to outgoing holes
and vice versa, and finally (iii) elastic scattering by the impurity (light
blue strip) at position $\alpha L$. Below we will consider these three
ingredients in detail and explain how we combine them to find the energies of
the Andreev bound states.

Figure 2: (a) Diagram summarizing all scattering processes. The vertical blue
shaded region represents the scatterer, where electrons and holes can be
transmitted or reflected. The horizontal arrows depict the free-propagation
amplitudes of the electrons (solid) and holes (dashed), see Eq. (7), where the
small vertical arrows indicate the two spin directions in the local
eigenbasis. Andreev reflection at the interfaces is represented by the bent
amplitudes in the gray shaded areas, connecting electron to hole amplitudes
and including spin mixing due to our choice of basis in the normal region. (b)
Combinations of scattering processes that are combined into the matrices
$M_{L,R}$ and $\tilde{M}_{L,R}$, see Eqs. (12,13).

Free propagation—We start by assuming translational invariance in the two non-
proximitized parts of the wire on both sides of the scatterer, i.e., for
$0\leq x<\alpha L$ and $\alpha L<x\leq L$, which allows us to write for these
regions

$\displaystyle\hat{\cal
H}_{0}^{L,R}=\frac{\hbar^{2}k^{2}}{2m}-\mu+({\bf{B}}_{\text{Z}}+k\,\boldsymbol{\alpha}_{\rm
so})\cdot{\boldsymbol{\sigma}},$ (5)

in terms of the wave vector $k$. Assuming that $\Delta_{0}\ll\mu$ and
anticipating that we will focus on the subgap physics of the junction, we
linearize this Hamiltonian around the Fermi wave vector,

$\displaystyle\hat{\cal H}_{0,\eta}^{L,R}\approx E_{\rm F}+\hbar v_{\rm
F}(\eta k-k_{\rm F})+({\bf{B}}_{\text{Z}}+\eta{\bf B}_{\rm
so})\cdot{\boldsymbol{\sigma}},$ (6)

where $\eta=\pm 1$ corresponds to electrons moving in the
$\pm\hat{x}$-direction and $v_{\rm F}$ is the Fermi velocity. We also
introduced the effective spin–orbit field ${\bf B}_{\rm so}=k_{\rm
F}\boldsymbol{\alpha}_{\rm so}$, dropping the small correction $(k-k_{\rm
F})\boldsymbol{\alpha}_{\rm so}$.

We will describe the free propagation always in a “local” eigenbasis of spin,
i.e., always setting the spin quantization axis along
${\bf{B}}_{\text{Z}}+{\bf B}_{\rm so}$ and ${\bf{B}}_{\text{Z}}-{\bf B}_{\rm
so}$ depending on the direction of propagation of the particle under
consideration [see Fig. 1(b)]. Using all the conventions mentioned above, we
then find that the transfer matrices that describe the propagation (see Fig.
2) can be written as

$t_{e,h}^{L\pm}=e^{i\sigma_{z}\chi_{\pm}^{e,h}\alpha}\quad\text{and}\quad
t_{e,h}^{R\pm}=e^{i\sigma_{z}\chi_{\pm}^{e,h}(1-\alpha)},$ (7)

using the phases $\chi^{e}_{\pm}=-V_{\pm}L/\hbar v_{\rm F}+EL/\hbar v_{\rm
F}+k_{\rm F}L$ and $\chi^{h}_{\pm}=-V_{\pm}L/\hbar v_{\rm F}-EL/\hbar v_{\rm
F}-k_{\rm F}L$ with $V_{\pm}=|{\bf B}_{\rm Z}\pm{\bf B}_{\rm so}|$, $E$ the
energy of the bound state, and $k_{\rm F}$ the Fermi wave number. Since we
work in the short-junction limit, we anticipate that the contributions
$EL/\hbar v_{\rm F}$ can be neglected. Furthermore, we see from Fig. 2 that
for the case of perfect Andreev reflection at the SN-interfaces the phase
factors $\pm k_{\rm F}L$ cancel since they always appear with opposite signs
in Andreev-reflected pairs of amplitudes. Under our assumptions, it thus
suffices to focus on the phases $V_{\pm}$, and we can drop the indices $e,h$
in Eq. (7).

Andreev reflection—The price to pay for aligning the spin quantization axis
with the local fields is that all reflection processes become effectively
spin-mixing. The Andreev reflection coefficients thus become matrices that
include a spin rotation that depends on the misalignment of the fields ${\bf
B}_{\rm Z}\pm{\bf B}_{\rm so}$. We find that the four matrices can be written
as

$\displaystyle r_{he}^{L}={}$
$\displaystyle{}e^{-i\alpha}R(-\theta)e^{-i\phi_{l}},$ $\displaystyle
r_{eh}^{L}={}$ $\displaystyle{}e^{-i\alpha}R(-\theta)e^{i\phi_{l}},$ (8)
$\displaystyle r_{he}^{R}={}$
$\displaystyle{}e^{-i\alpha}R(\theta)e^{-i\phi_{r}},$ $\displaystyle
r_{eh}^{R}={}$ $\displaystyle{}e^{-i\alpha}R(\theta)e^{i\phi_{r}},$ (9)

where $\alpha=\arccos(E/\Delta_{0})$ and

$R(\theta)=\begin{pmatrix}\cos(\frac{\theta}{2})&\sin(\frac{\theta}{2})\\\
-\sin(\frac{\theta}{2})&\cos(\frac{\theta}{2})\end{pmatrix},$ (10)

implements a spin rotation over the angle

$\theta=\arccos\left(\frac{B_{\rm Z}^{2}-B_{\rm so}^{2}}{V_{+}V_{-}}\right),$
(11)

i.e., the angle between ${\bf B}_{\rm Z}+{\bf B}_{\rm so}$ and ${\bf B}_{\rm
Z}-{\bf B}_{\rm so}$, as indicated in Fig. 1(b).

We can now construct the matrices

$\displaystyle M_{L}={}$ $\displaystyle{}t_{e}^{L+}r_{eh}^{L}t_{h}^{L-},$
$\displaystyle\tilde{M}_{L}={}$
$\displaystyle{}t_{h}^{L+}r_{he}^{L}t_{e}^{L-},$ (12) $\displaystyle M_{R}={}$
$\displaystyle{}t_{h}^{R-}r_{he}^{R}t_{e}^{R+},$
$\displaystyle\tilde{M}_{R}={}$
$\displaystyle{}t_{e}^{R-}r_{eh}^{R}t_{h}^{R+},$ (13)

which fully describe all propagation on the left and right side of the
scatterer, see also Fig. 2(b).

Impurity scattering—Finally we consider the scatterer, which we assume to be
spin-independent and particle–hole symmetric. Defining the vector of incoming
and outgoing electron amplitudes as ${\bf
a}_{e}=(a^{e}_{1},a^{e}_{2},a^{e}_{3},a^{e}_{4})^{T}$ and ${\bf
b}_{e}=(b^{e}_{1},b^{e}_{2},b^{e}_{3},b^{e}_{4})^{T}$, respectively, and
similarly for the incoming and outgoing hole amplitudes [see Fig. 2(a)], we
thus assume that the effect of the scatterer is to relate the amplitudes as
${\bf b}_{e,h}=S{\bf a}_{e,h}$ with

$S=\begin{pmatrix}T\,\mathbb{1}&-\sqrt{1-T^{2}}\,R(\theta)\\\
\sqrt{1-T^{2}}\,R(-\theta)&T\,\mathbb{1}\end{pmatrix},$ (14)

where $T^{2}$ is the transmission probability of the scatterer. The rotation
matrices included with the probability amplitudes for reflection again take
care of the misalignment of the spin bases for left- and right-moving
particles.

Bound state—We now see that, with the notation introduced before, we can also
relate

${\bf a}_{h}=\begin{pmatrix}0&M_{R}\\\ \tilde{M}_{L}&0\end{pmatrix}{\bf
b}_{e},\quad{\bf a}_{e}=\begin{pmatrix}0&\tilde{M}_{R}\\\
M_{L}&0\end{pmatrix}{\bf b}_{h}.$ (15)

Combining Eq. (15) with the relation ${\bf b}_{e,h}=S{\bf a}_{e,h}$ we find an
eigenvalue problem that defines the bound states in the junction. It thus
follows that the equation

$\text{Det}\left[1-S\begin{pmatrix}0&\tilde{M}_{R}\\\
M_{L}&0\end{pmatrix}S\begin{pmatrix}0&M_{R}\\\
\tilde{M}_{L}&0\end{pmatrix}\right]=0,$ (16)

in general constitutes a self-consistency equation for the energy $E$ of the
bound state. In our case, however, where $E$ only appears in the prefactor of
all Andreev reflection coefficients, it is possible to write down a closed-
form expression for the allowed energies.

##  III Results

Figure 3: Calculated bound-state energies (top) and supercurrent (bottom)
using Eqs. (17) and (20), respectively. The green dots present known results
from the literature, for comparison. We set $\alpha=1/\sqrt{2}$ throughout and
further we used: (a) $\chi_{\rm Z}=\chi_{\rm so}=0$ and $T=1$ (blue) or
$T=\sqrt{0.85}$ (green); (b) $\chi_{\rm Z}=0.2\,\pi$, $\chi_{\rm so}=0$, and
$T=1$ (blue) or $T=\sqrt{0.95}$ (green); (c) $\chi_{\rm Z}=0$, $\chi_{\rm
so}=0.2\,\pi$, and $T=1$ (blue) or $T=\sqrt{0.85}$ (green); (d) $\chi_{\rm
Z}=0.2\,\pi$, $\chi_{\rm so}=0.8\,\pi$, $T=\sqrt{0.95}$, and ${\bf B}_{\rm
Z}\perp{\bf B}_{\rm so}$ (blue), and $\chi_{\rm Z}=\chi_{\rm so}=0.2\,\pi$,
$T=\sqrt{0.85}$, and an angle of $\pi/4$ between the two fields.

Using the explicit expressions for the matrices given above, it is possible to
rewrite Eq. (16) in the form

$E(\phi)=\pm\frac{\Delta_{0}}{\sqrt{2}}\sqrt{v\pm\sqrt{w}},$ (17)

where

$\displaystyle v={}$ $\displaystyle{}1-q_{1}(T)+q_{2}(T)\cos(\phi),$ (18)
$\displaystyle w={}$
$\displaystyle{}\tfrac{1}{2}+[q_{1}(T)-q_{2}(T)\cos(\phi)]^{2}$
$\displaystyle{}\hskip 10.00002pt-
q_{3}(T)-q_{4}(T)\cos(\phi)-T^{4}\cos^{2}(\phi),$ (19)

revealing the functional dependence of the energies on the phase difference
$\phi$. The coefficients $q_{i}$ are cumbersome, but straightforward functions
of $T$ as well as $\theta$, $\alpha$, and $\chi_{\pm}$, and are given
explicitly in App. A. Although the lengthy expressions for the $q_{i}$ make
the general result for $E(\phi)$ slightly obscure at first sight, we believe
that it is in fact very useful in several regards: Firstly, the expression
(17) is written in such a way that the positive- and negative-energy solutions
can be distinguished easily, allowing for an explicit analytic evaluation of
the zero-temperature equilibrium supercurrent via (4),

$I(\phi)=-\frac{e\Delta_{0}}{\hbar}\frac{1}{\sqrt{2}}\sum_{\eta=\pm
1}\frac{(\partial_{\phi}v)+\frac{1}{2}\eta(\partial_{\phi}w)/\sqrt{w}}{\sqrt{v+\eta\sqrt{w}}},$
(20)

where the differentiation with respect to $\phi$ yields

$\displaystyle\partial_{\phi}v={}$ $\displaystyle{}-q_{2}(T)\sin(\phi),$ (21)
$\displaystyle\partial_{\phi}w={}$
$\displaystyle{}[2q_{1}(T)q_{2}(T)+q_{4}(T)]\sin(\phi)$ $\displaystyle{}\hskip
10.00002pt+[T^{4}-q_{2}(T)^{2}]\sin(2\phi).$ (22)

Secondly, using the explicit expressions given in the Appendix, it is
straightforward to investigate limiting cases, e.g., of small or large $|{\bf
B}_{\rm Z}|/|{\bf B}_{\rm so}|$ or $T$, or to simplify the general result for
a specific location for the scatterer, such as $\alpha=0$ or $\alpha=1/2$.

Below we will present our results, based on Eq. (17), in the following way.
First, we will explore the general expressions for the bound-state energies
and supercurrent, investigate their dependence on our model parameters, and
compare the results with the expected behavior. Then we will consider the
limits of $T\ll 1$ and $T=1$, where it is straightforward to write down
analytic expressions for the critical current $I_{c}={\rm max}_{\phi}I(\phi)$,
which can be connected to the magnitude and relative orientation of the fields
${\bf B}_{\rm Z}$ and ${\bf B}_{\rm so}$. Finally we will illustrate how
fixing $\alpha$ to a specific value—we will pick $\alpha=0$ and
$\alpha=1/2$—can allow to greatly simplify the general result, yielding
digestible expressions for the bound-state energies and the supercurrent
through the junction.

###  III.1 General case

In Fig. 3 we plot the four bound-state energies (top row) and supercurrent
through the junction (bottom row) as a function of the phase difference
$\phi$, as given by Eqs. (17) and (20), respectively. In Fig. 3(a) we set both
phases $\chi_{\rm Z}\equiv|{\bf B}_{\rm Z}|L/\hbar v_{\rm F}$ and $\chi_{\rm
so}\equiv|{\bf B}_{\rm so}|L/\hbar v_{\rm F}$ to zero, to see if our
expressions reproduce the usual results. The blue curves are for $T=1$ and the
green curves for $T=\sqrt{0.85}$, and in both cases we set $\alpha=1/\sqrt{2}$
(which we keep at this value throughout the whole Figure). We see that the
energies indeed follow the expected $E(\phi)=\pm\Delta_{0}\cos(\phi/2)$ for
the fully transparent junction and open up an anticrossing around $\phi=\pi$
for lower transparency. The green dots show the energies given by Beenakker’s
well-known expression $E(\phi)=\pm\Delta_{0}\sqrt{1-T^{2}\sin(\phi/2)^{2}}$
for a short single-channel Josephson junction [35]. In Fig. 3(b) we add a
finite Zeeman field, setting $\chi_{\rm Z}=0.2\,\pi$, while keeping $\chi_{\rm
so}=0$. As expected, the doubly (spin) degenerate levels of Fig. 3(a) now
split. For $T=1$ (blue curves), the effect of the Zeeman field is simply an
extra phase shift that is acquired by electrons and holes crossing the
junction, opposite for the two spin directions [53]. The bound-state energies
thus follow $E(\phi)=\pm\Delta_{0}\cos(\phi/2\pm\chi_{\rm Z})$, and the system
can undergo a $0$–$\pi$ transition as the Zeeman field is increased [40, 59,
60, 61, 41, 62, 2]. In the presence of a finite scattering probability, level
anticrossings will emerge. This is illustrated by the green curve, where we
set $T=\sqrt{0.95}$. The green dots show the energies as given by the simple
expression derived in Ref. [54] [Eqs. (23,24), valid in the limit $\chi_{\rm
so}=0$], which agree with our more general expression. In Fig. 3(c) we show
the opposite case of finite spin–orbit coupling, $\chi_{\rm so}=0.2\,\pi$ and
zero Zeeman splitting, where the transparency is set to $T=1$ (blue) and
$T=\sqrt{0.85}$ (green). These results are identical to the effectively
spinless case shown in Fig. 3(a), where $\chi_{\rm Z}=\chi_{\rm so}=0$.
Indeed, in the short-junction limit, phase shifts purely due to spin–orbit
coupling cancel, and deviations from the spinless case appear only in higher
order in $\Delta_{0}\tau_{\rm dwell}/\hbar$, where $\tau_{\rm dwell}$ is the
dwell time of the current-carrying particles in the junction [63, 64].
Finally, in Fig. 3(d) we explore more general parameter settings: In blue we
show the case $\chi_{\rm Z}=0.2\,\pi$, $\chi_{\rm so}=0.8\,\pi$, and
$T=\sqrt{0.95}$, where the fields ${\bf B}_{\rm Z}$ and ${\bf B}_{\rm so}$
were assumed perpendicular, and in green the case $\chi_{\rm Z}=\chi_{\rm
so}=0.2\,\pi$, and $T=\sqrt{0.85}$, now using an angle of $\pi/4$ between the
two fields. In all columns, the lower panel shows the corresponding
current–phase relationship at zero temperature, which in principle follows
from Eq. (4) as the phase derivative of the sum of all negative bound-state
energies. We emphasize again that, due to the form of Eq. (17), our expression
for the supercurrent (20) is written in a closed form and does not require a
separate determination of which of the four bound-state energies are negative.
We see that we never find an anomalous Josephson current (a finite
supercurrent at $\phi=0$), which is already clear from the observation that
the bound-state spectrum is symmetric in $\phi$, i.e., $E(\phi)=E(-\phi)$
[54].

The results presented in this Section show how Eqs. (17) and (20) can be
straightforwardly used to calculate the bound-state spectrum as well as the
current–phase relationship in a single-channel SNS-junction, for any
combination of effective Zeeman and spin–orbit fields and across the whole
range of $T$, from a fully transparent to a fully opaque junction. We found
that our results exhibit the expected behavior as a function of all
parameters, and that our expressions reproduce the known results in limiting
cases. The fact that our analytic expressions are most general presents the
great advantage that they can be used to explore additional limits as well, by
straightforward simplification. Below, we will first consider the fully
transparent limit $T=1$, which has been investigated before in the literature
[55], and we will show how our expressions reproduce the expected results.
Then we will explore the experimentally relevant limit of $T\ll 1$, deriving
relatively compact expressions for the spectrum and the supercurrent, which
also allow us to find an analytic expression for the critical current through
the junction. We will derive compact analytic expressions for the current in
the limiting cases of strong spin–orbit coupling and vanishing spin–orbit
coupling. Finally, we will show how setting $\alpha$ to a specific value
highly reduces the complexity of the expressions, yielding compact results
that are valid for any $0<T\leq 1$. As an example, we give the results for
$\alpha=0$ (most scattering takes place at one of the interfaces) and
$\alpha=1/2$ (the scatterer is placed in the center of the wire).

###  III.2 High-transparency limit, $T=1$

In the case of a fully transparent wire with $T=1$, the location of the
scatterer becomes irrelevant and Eq. (17) can be simplified to

$E(\phi)=\pm\frac{\Delta_{0}}{\sqrt{2}}\sqrt{1+\cos(\phi\pm\zeta)},$ (23)

with the angle $0\leq\zeta\leq\pi$ defined through

$\cos(\zeta)=\cos(\chi_{+})\cos(\chi_{-})-\cos(\theta)\sin(\chi_{+})\sin(\chi_{-}).$
(24)

We decided to express the energies in the same form as in Eq. (17), but it is
clear that this result is equivalent to
$E(\phi)=\pm\cos[\frac{1}{2}(\phi\pm\zeta)]$, as, e.g., derived in Ref. [55].
In this fully-transparent limit the levels show the regular
$\cos(\phi/2)$-behavior, with a phase shift equal to $\pm\zeta$, which is the
result of the interplay of the effects of the Zeeman and spin–orbit fields on
the phase of the propagating particles. This is illustrated in Fig. 4(a),
where we plot the bound-state spectrum as given by Eq. (23) with $\chi_{\rm
Z}=0.4\,\pi$ and $\chi_{\rm so}=0.2\,\pi$ (blue curves) and $\chi_{\rm
Z}=0.2\,\pi$ and $\chi_{\rm so}=0.4\,\pi$ (green curves), setting the angle
between the two fields to $\pi/4$ in both cases.

Figure 4: Calculated (a) bound-state energies and (b) supercurrent using Eqs.
(23) and (25). Blue: $\chi_{\rm Z}=0.4\,\pi$ and $\chi_{\rm so}=0.2\,\pi$.
Green: $\chi_{\rm Z}=0.2\,\pi$ and $\chi_{\rm so}=0.4\,\pi$. In both cases we
set the angle between the fields to $\pi/4$.

The equilibrium supercurrent then follows as

$I(\phi)=\frac{e\Delta_{0}}{2\hbar}\sum_{\eta=\pm
1}\frac{\sin(\phi+\eta\,\zeta)}{|\cos[\frac{1}{2}(\phi+\eta\,\zeta)]|},$ (25)

cf. the expressions given in Ref. [55]. Fig. 4(b) shows the supercurrent
following from this expression, for the same two sets of parameters as used in
the left pane. A straightforward inspection of this result yields that the
corresponding critical current in the fully transparent case reads as

$I_{c}=\frac{e\Delta_{0}}{\hbar}[1+|\cos(\zeta)|],$ (26)

in line with the results presented in Ref. [55]. For comparison, we added the
values of $\pm I_{c}$ predicted by Eq. (26) to the bottom pane of Fig. 4 as
horizontal dashed lines.

We thus see that our expressions reproduce the known results for the
transparent case. The unorthodox form of Eq. (23), based on that of Eq. (17),
allows for a more straightforward evaluation of the ground-state energy and
thus the equilibrium supercurrent since the sign of each of the energies is
obvious from the expression.

###  III.3 Low-transparency limit, $T\ll 1$

Figure 5: Critical current in the opaque regime. (a) Correction to the
critical current in the opaque regime in the presence of strong spin–orbit
coupling. We plot $i_{c}$ to leading order in $\chi_{\rm Z}/\chi_{\rm so}$ in
the limit of $\chi_{\rm Z}/\chi_{\rm so}\ll 1$, as given by Eq. (29), as a
function of the angle $\theta$ between ${\bf B}_{\rm Z}$ and ${\bf B}_{\rm
so}$ (polar angle) and $\chi_{\rm Z}$ (radial coordinate) for $\chi_{\rm
so}=0.4\,\pi$ and $\alpha=0.4$. (b) Plot of the same correction $i^{(2)}_{c}$
for the same $\chi_{\rm so}$ and $\alpha$, at $\chi_{\rm Z}=0.2$, i.e., at the
outer circle of the plot shown in (a). (c) The critical current $i_{c}$ as a
function of $\chi_{\rm Z}$ in the opposite limit of vanishing spin–orbit
coupling, where we used $\alpha=0.1$. (d) The current-phase relationship for
$\chi_{\rm Z}=0.9\,\pi$ [vertical dotted line in (c)], $\chi_{\rm so}=0$,
$T=0.1$, and $\alpha=0.1$.

For the case of a relatively opaque junction, we take $T\ll 1$ and thus expand
Eqs. (17) and (20) to leading order in $T$. This yields, as expected for low
transparency, a sinusoidal current–phase relationship,

$I(\phi)=\frac{e\Delta_{0}}{\hbar}i_{c}T^{2}\sin(\phi),$ (27)

where the dependence of the current on the Zeeman and spin–orbit fields is
captured by

$i_{c}=\sum_{\eta=\pm
1}\frac{q_{2}(1)+\eta\frac{2q_{4}(\frac{1}{2}\sqrt{2})+4q_{1}(\frac{1}{2}\sqrt{2})q_{2}(\frac{1}{2}\sqrt{2})}{\sqrt{\frac{1}{2}+q_{1}(0)^{2}-q_{3}(0)}}}{\sqrt{2-2q_{1}(0)-2\eta\sqrt{\frac{1}{2}+q_{1}(0)^{2}-q_{3}(0)}}},$
(28)

written again in terms of the coefficients $q_{i}(T)$ as given in App. A.
Although thus cumbersome, this analytic expression for the critical current
can be used to reveal useful insights in the dependence of the current on the
Zeeman and spin–orbit fields and their relative orientation.

For example, for the case of a junction with strong spin–orbit coupling we can
investigate perturbatively the effect of the application of a small Zeeman
field in addition. We thus expand Eq. (28) to leading order in $\chi_{\rm
Z}/\chi_{\rm so}$, assuming also $\chi_{\rm Z}\ll 1$, yielding $i_{c}\approx
1+i_{c}^{(2)}$ with

$\displaystyle i_{c}^{(2)}={}$ $\displaystyle{}-2\frac{\chi_{\rm
Z}^{2}}{\chi_{\rm so}^{2}}\cos(\chi_{\rm so})\sin(\alpha\chi_{\rm
so})\sin(\tilde{\alpha}\chi_{\rm so})\sin^{2}(\theta)$
$\displaystyle{}-2\chi_{\rm Z}^{2}\tilde{\alpha}\alpha\cos^{2}(\theta),$ (29)

where $\tilde{\alpha}=1-\alpha$. In Fig. 5 we show the resulting second-order
correction for $\chi_{\rm so}=0.4\,\pi$ and $\alpha=0.4$, (a) as a function of
both $\theta$ and $\chi_{\rm Z}$ and (b) as a function of $\theta$ for
$\chi_{\rm Z}=0.2$. We see that the (negative) correction is maximal for ${\bf
B}_{\rm Z}$ and ${\bf B}_{\rm so}$ (anti)aligned and minimal when the two
fields are perpendicular. From Eq. (29) it is indeed easy to see that the
maximum correction $-2\chi_{\rm Z}^{2}\tilde{\alpha}\alpha$ is reached at
$\theta=0,\pi$ and the minimum correction $-2(\chi_{\rm Z}/\chi_{\rm
so})^{2}\cos(\chi_{\rm so})\sin(\alpha\chi_{\rm
so})\sin(\tilde{\alpha}\chi_{\rm so})$ at
$\theta=\frac{\pi}{2},\frac{3\pi}{2}$. In an experiment, one could thus not
only relate the directions of the applied ${\bf B}_{\rm Z}$ for which maxima
and minima in the critical current are observed to the direction of the
effective spin–orbit field, but one could also obtain information about the
(probably phenomenological) parameter $\alpha$ and about $\chi_{\rm so}$, by
following the development of the maxima and minima as a function of $\chi_{\rm
Z}$.

We can also consider the opposite case of vanishing spin–orbit interaction,
where the critical current $i_{c}$ becomes a function of $\alpha$ and
$\chi_{\rm Z}$ alone,

$\displaystyle i_{c}=\frac{\tan(2\alpha\chi_{\rm Z})|\cos(2\alpha\chi_{\rm
Z})|-\tan(2\tilde{\alpha}\chi_{\rm Z})|\cos(2\tilde{\alpha}\chi_{\rm
Z})|}{\sin(2[\alpha-\tilde{\alpha}]\chi_{\rm Z})}.$ (30)

In Fig. 5(c) we show the critical current following from Eq. (30) for the case
$\alpha=0.1$. As an example, we plot in Fig. 5(d) the current–phase
relationship for $\chi_{\rm Z}=0.9\,\pi$, $\chi_{\rm so}=0$, $T=0.1$, and
$\alpha=0.1$, as calculated from Eq. (20), where the blue dashed lines
represent the critical current as given by Eq. (30). The vertical dotted line
in Fig. 5(c) indicates the $\chi_{\rm Z}$ that was used in Fig. 5(d). We note
that the sign changes in $i_{c}$ correspond to $0$–$\pi$ transitions in the
junction.

###  III.4 Fixed position of the scatterer

Finally, we show how the general expression given for the bound-state spectrum
(17) simplifies when the location of the scatterer is fixed. Here we choose
the particularly simple values of $\alpha=0$ and $\alpha=1/2$, but similar
simplifications can be made for any value of $\alpha$.

When $\alpha=0$, i.e., the scattering takes places mostly at one end of the
wire, then the coefficients $q_{i}$ as given in App. A simplify considerably,
yielding

$E(\phi)=\pm\frac{\Delta_{0}}{\sqrt{2}}\sqrt{v_{0}\pm\sqrt{w_{0}}},$ (31)

with now

$\displaystyle v_{0}={}$
$\displaystyle{}1+\cos(\zeta)[\cos(\zeta)R^{2}+T^{2}\cos(\phi)],$ (32)
$\displaystyle w_{0}={}$
$\displaystyle{}\sin(\zeta)^{2}\big{(}1-[\cos(\zeta)R^{2}+T^{2}\cos(\phi)]^{2}\big{)},$
(33)

using the same angle $\zeta$ as before [see Eq. (24)] and the notation
$R^{2}=1-T^{2}$. Defining

$\cos(\kappa)=\cos(\zeta)R^{2}+T^{2}\cos(\phi),$ (34)

with $0\leq\kappa\leq\pi$, we see that the energies become

$E(\phi)=\pm\Delta_{0}\cos[\tfrac{1}{2}(\kappa\pm\zeta)],$ (35)

where the phase-dependence is implicit via the angle $\kappa$. This allows us
to write a compact expression for the equilibrium supercurrent as well,

$I(\phi)=\frac{e\Delta_{0}}{2\hbar}\sum_{\eta=\pm
1}\frac{T^{2}\sin(\kappa+\eta\,\zeta)\sin(\phi)}{|\cos[\frac{1}{2}(\kappa+\eta\,\zeta)]\sin(\kappa)|},$
(36)

Another simple choice for $\alpha$ that leads to relatively compact
expressions is $\alpha=1/2$, i.e., the scatterer placed in the center of the
wire. In this case, the energy reads again as

$E(\phi)=\pm\frac{\Delta_{0}}{\sqrt{2}}\sqrt{v_{\frac{1}{2}}\pm\sqrt{w_{\frac{1}{2}}}},$
(37)

with now

$\displaystyle v_{\frac{1}{2}}={}$
$\displaystyle{}1+\cos(\zeta)[R^{2}+T^{2}\cos(\phi)]-2\sin^{2}(\lambda)R^{2},$
(38) $\displaystyle w_{\frac{1}{2}}={}$
$\displaystyle{}16\sin^{2}(\zeta/2)T^{2}\sin^{2}(\phi/2)$ (39)
$\displaystyle{}\hskip
10.00002pt\times\left\\{\cos^{2}(\zeta/2)[1-T^{2}\sin^{2}(\phi/2)]-\sin^{2}(\lambda)R^{2}\right\\},$

using the same definitions as before and

$\displaystyle\sin(\lambda)=\sin(\chi_{+}/2)\sin(\chi_{-}/2)\sin(\theta).$
(40)

In this case, the current–phase relationship can still be straightforwardly
calculated, but does not reduce to an expression as elegant as Eq. (36).

Other values of $\alpha$ can also straightforwardly be explored using the
explicit results given in App. A; the resulting expressions will most likely
be of similar complexity as the ones for $\alpha=0,1/2$.

##  IV Conclusion

To conclude, we presented analytic expressions for the Andreev bound-state
energies and supercurrent in a single-channel semiconductor-based SNS-
junction. We took into account the effects of finite spin–orbit interaction
and a Zeeman splitting, allowing the two associated fields to have any desired
orientation, and we included a single scatterer that can be arbitrarily
located in the normal region. We first corroborated our results by comparing
them to all limiting cases for which analytic expressions are known. Then we
demonstrated how our general expressions can be simplified in several
additional limits, such as that of small and large spin–orbit coupling in the
low-transparency case, resulting in insightful expressions that could be
helpful for disentangling spin–orbit and Zeeman effects in experiment.

We gratefully acknowledge financial support via NTNU’s Onsager Fellow program.

##  Appendix A Explicit expressions for $q_{i}$

The four coefficients $q_{i}$ read explicitly as

$\displaystyle q_{1}(T)={}$
$\displaystyle{}R^{2}\sin(\theta)^{2}[\sin(\alpha\chi_{+})^{2}\sin(\alpha\chi_{-})^{2}+\sin(\tilde{\alpha}\chi_{+})^{2}\sin(\tilde{\alpha}\chi_{-})^{2}]-\tfrac{1}{2}R^{2}[f_{2}(2\alpha\chi_{+},2\alpha\chi_{-})+f_{2}(2\tilde{\alpha}\chi_{+},2\tilde{\alpha}\chi_{-})],$
(41) $\displaystyle q_{2}(T)={}$
$\displaystyle{}T^{2}f_{2}(\chi_{+},\chi_{-}),$ (42) $\displaystyle
q_{3}(T)={}$
$\displaystyle{}\tfrac{1}{4}R^{2}T^{2}\sin(\theta)^{2}[f_{2}(2\alpha\chi_{+},-2\tilde{\alpha}\chi_{-})+f_{2}(2\alpha\chi_{-},-2\tilde{\alpha}\chi_{+})]-\tfrac{1}{8}R^{4}\sin(\theta)^{4}-4T^{4}\sin(\theta)^{2}\sin(\chi_{+})^{2}\sin(\chi_{-})^{2}$
$\displaystyle{}+\tfrac{1}{4}R^{4}\sin(\theta)^{2}[2+f_{2}(2\alpha\chi_{+},-2\beta\chi_{-})+f_{2}(2\tilde{\alpha}\chi_{+},2\beta\chi_{-})+f_{2}(2\alpha\chi_{-},-2\beta\chi_{+})+f_{2}(2\tilde{\alpha}\chi_{-},2\beta\chi_{+})$
$\displaystyle{}\hskip 60.00009pt-
f_{2}(2\alpha\chi_{+},2\alpha\chi_{-})-f_{2}(2\tilde{\alpha}\chi_{+},2\tilde{\alpha}\chi_{-})]$
$\displaystyle{}+\tfrac{1}{4}R^{2}\sin(\theta)^{2}[f_{2}(2\chi_{+},2\alpha\chi_{-})+f_{2}(2\chi_{+},2\tilde{\alpha}\chi_{-})+f_{2}(2\chi_{-},2\alpha\chi_{+})+f_{2}(2\chi_{-},2\tilde{\alpha}\chi_{+})$
$\displaystyle{}\hskip 60.00009pt-
f_{2}(2\alpha\chi_{+},2\alpha\chi_{-})-f_{2}(2\tilde{\alpha}\chi_{+},2\tilde{\alpha}\chi_{-})]$
$\displaystyle{}+R^{2}T^{2}f_{6}(2\chi_{+},2\chi_{-})+\tfrac{1}{2}R^{4}[f_{8}(2\beta\chi_{+},2\beta\chi_{-})+f_{8}(2\chi_{+},2\chi_{-})+f_{8}(0,0)]+\tfrac{1}{2}T^{4}f_{2}(2\chi_{+},2\chi_{-})$
$\displaystyle{}+\tfrac{1}{4}\left(T^{2}+R^{2}\sin(\theta)^{2}[2\sin(\alpha\chi_{+})^{2}+2\sin(\tilde{\alpha}\chi_{+})^{2}-\sin(\chi_{+})^{2}-\sin(\beta\chi_{+})^{2}]\right)$
$\displaystyle{}\hskip
30.00005pt\times\left(T^{2}+R^{2}\sin(\theta)^{2}[2\sin(\alpha\chi_{-})^{2}+2\sin(\tilde{\alpha}\chi_{-})^{2}-\sin(\chi_{-})^{2}-\sin(\beta\chi_{-})^{2}]\right)$
$\displaystyle{}-\tfrac{1}{4}\left[T^{2}+\tfrac{1}{2}R^{2}\sin(\theta)^{2}\cos(2\beta\chi_{+})\right]\left[T^{2}+\tfrac{1}{2}R^{2}\sin(\theta)^{2}\cos(2\beta\chi_{-})\right]$
$\displaystyle{}-\tfrac{1}{4}\left[T^{2}-\tfrac{1}{2}R^{2}\sin(\theta)^{2}\cos(2\chi_{+})\right]\left[T^{2}-\tfrac{1}{2}R^{2}\sin(\theta)^{2}\cos(2\chi_{-})\right]+\tfrac{1}{4}\left[T^{2}-\tfrac{1}{2}R^{2}\sin(\theta)^{2}\right]^{2},$
(43) $\displaystyle q_{4}(T)={}$
$\displaystyle{}T^{2}R^{2}\left(2f_{2}(\beta\chi_{-},\beta\chi_{+})-\sin(\theta)^{2}[\cos(\chi_{+})-\cos(\beta\chi_{+})][\cos(\chi_{-})-\cos(\beta\chi_{-})]\right),$
(44)

where we used the function

$f_{n}(x,y)=\cos\left(\tfrac{1}{2}\theta\right)^{n}\cos(x+y)+\sin\left(\tfrac{1}{2}\theta\right)^{n}\cos(x-y),$
(45)

and the notation $\tilde{\alpha}=1-\alpha$, $\beta=\tilde{\alpha}-\alpha$, and
$R=\sqrt{1-T^{2}}$.

## References

  * Bercioux and Lucignano [2015] D. Bercioux and P. Lucignano, Rep. Prog. Phys. 78, 106001 (2015).
  * Cai  _et al._ [2023] R. Cai, I. Žutić, and W. Han, Adv. Quant. Tech. 6, 2200080 (2023).
  * Sau  _et al._ [2010] J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma, Phys. Rev. Lett. 104, 040502 (2010).
  * Oreg  _et al._ [2010] Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. 105, 177002 (2010).
  * Lutchyn  _et al._ [2010] R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. 105, 077001 (2010).
  * Lutchyn  _et al._ [2018] R. M. Lutchyn, E. P. A. M. Bakkers, L. P. Kouwenhoven, P. Krogstrup, C. M. Marcus, and Y. Oreg, Nat. Rev. Mater. 3, 52 (2018).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * Das Sarma  _et al._ [2015] S. Das Sarma, M. Freedman, and C. Nayak, npj Quant. Inf. 1, 15001 (2015).
  * 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).
  * Nadeem  _et al._ [2023] M. Nadeem, M. S. Fuhrer, and X. Wang, Nat. Rev. Phys. 5, 558 (2023).
  * Burkard  _et al._ [2020] G. Burkard, M. J. Gullans, X. Mi, and J. R. Petta, Nat. Rev. Phys. 2, 129 (2020).
  * Rasmussen  _et al._ [2016] A. Rasmussen, J. Danon, H. Suominen, F. Nichele, M. Kjaergaard, and K. Flensberg, Phys. Rev. B 93, 155406 (2016).
  * Mazur  _et al._ [2022] G. P. Mazur, N. van Loo, D. van Driel, J.-Y. Wang, G. Badawy, S. Gazibegovic, E. P. A. M. Bakkers, and L. P. Kouwenhoven, arXiv:2211.14283 (2022).
  * Baumgartner  _et al._ [2022a] C. Baumgartner, L. Fuchs, A. Costa, S. Reinhardt, S. Gronin, G. C. Gardner, T. Lindemann, M. J. Manfra, P. E. Faria Junior, D. Kochan, J. Fabian, N. Paradiso, and C. Strunk, Nat. Nano. 17, 39 (2022a).
  * Baumgartner  _et al._ [2022b] C. Baumgartner, L. Fuchs, A. Costa, J. Picó-Cortés, S. Reinhardt, S. Gronin, G. C. Gardner, T. Lindemann, M. J. Manfra, P. E. F. Junior, D. Kochan, J. Fabian, N. Paradiso, and C. Strunk, J. Phys.: Cond. Matt. 34, 154005 (2022b).
  * Davydova  _et al._ [2022] M. Davydova, S. Prembabu, and L. Fu, Science Adv. 8, eabo0309 (2022).
  * Souto  _et al._ [2022] R. S. Souto, M. Leijnse, and C. Schrade, Phys. Rev. Lett. 129, 267702 (2022).
  * Fominov and Mikhailov [2022] Y. V. Fominov and D. S. Mikhailov, Phys. Rev. B 106, 134514 (2022).
  * Valentini  _et al._ [2023] M. Valentini, O. Sagi, L. Baghumyan, T. de Gijsel, J. Jung, S. Calcaterra, A. Ballabio, J. A. Servin, K. Aggarwal, M. Janik, T. Adletzberger, R. S. Souto, M. Leijnse, J. Danon, C. Schrade, E. Bakkers, D. Chrastina, G. Isella, and G. Katsaros, arXiv:2306.07109 (2023).
  * Ciaccia  _et al._ [2023] C. Ciaccia, R. Haller, A. C. C. Drachmann, T. Lindemann, M. J. Manfra, C. Schrade, and C. Schönenberger, Phys. Rev. Res. 5, 033131 (2023).
  * Maiani  _et al._ [2023] A. Maiani, K. Flensberg, M. Leijnse, C. Schrade, S. Vaitiekėnas, and R. Seoane Souto, Phys. Rev. B 107, 245415 (2023).
  * Strambini  _et al._ [2020] E. Strambini, A. Iorio, O. Durante, R. Citro, C. Sanz-Fernández, C. Guarcello, I. V. Tokatly, A. Braggio, M. Rocci, N. Ligato, V. Zannier, L. Sorba, F. S. Bergeret, and F. Giazotto, Nat. Nano. 15, 656 (2020).
  * Mayer  _et al._ [2020] W. Mayer, M. C. Dartiailh, J. Yuan, K. S. Wickramasinghe, E. Rossi, and J. Shabani, Nat. Commun. 11, 212 (2020).
  * Kjaergaard  _et al._ [2020] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, Ann. Rev. Cond. Matt. Phys. 11, 369 (2020).
  * Pientka  _et al._ [2017] F. Pientka, A. Keselman, E. Berg, A. Yacoby, A. Stern, and B. I. Halperin, Phys. Rev. X 7, 021032 (2017).
  * Fornieri  _et al._ [2019] A. Fornieri, A. M. Whiticar, F. Setiawan, E. Portolés, A. C. C. Drachmann, A. Keselman, S. Gronin, C. Thomas, T. Wang, R. Kallaher, G. C. Gardner, E. Berg, M. J. Manfra, A. Stern, C. M. Marcus, and F. Nichele, Nature 569, 89 (2019).
  * Hell  _et al._ [2017] M. Hell, M. Leijnse, and K. Flensberg, Phys. Rev. Lett. 118, 107701 (2017).
  * Lidal and Danon [2023] J. Lidal and J. Danon, Phys. Rev. B 107, 085303 (2023).
  * San-Jose  _et al._ [2013] P. San-Jose, J. Cayao, E. Prada, and R. Aguado, New J. Phys. 15, 075019 (2013).
  * San-Jose  _et al._ [2014] P. San-Jose, E. Prada, and R. Aguado, Phys. Rev. Lett. 112, 137001 (2014).
  * Murthy  _et al._ [2020] C. Murthy, V. D. Kurilovich, P. D. Kurilovich, B. van Heck, L. I. Glazman, and C. Nayak, Phys. Rev. B 101, 224501 (2020).
  * Dominguez  _et al._ [2022] F. Dominguez, E. G. Novik, and P. Recher, arXiv:2210.02065 (2022).
  * Golubov  _et al._ [2004] A. A. Golubov, M. Y. Kupriyanov, and E. Il’ichev, Rev. Mod. Phys. 76, 411 (2004).
  * Kulik and Omelyanchuk [1975] I. O. Kulik and A. N. Omelyanchuk, Sov. Phys. JETP 21, 96 (1975).
  * Beenakker [1991] C. W. J. Beenakker, Phys. Rev. Lett. 67, 3836 (1991).
  * Beenakker and van Houten [1991] C. W. J. Beenakker and H. van Houten, Phys. Rev. Lett. 66, 3056 (1991).
  * Kresin [1986] V. Z. Kresin, Phys. Rev. B 34, 7587 (1986).
  * Altshuler and Spivak [1987] B. L. Altshuler and B. Z. Spivak, Sov. Phys. JETP 65, 343 (1987).
  * Furusaki  _et al._ [1991] A. Furusaki, H. Takayanagi, and M. Tsukada, Phys. Rev. Lett. 67, 132 (1991).
  * Buzdin [2005] A. I. Buzdin, Rev. Mod. Phys. 77, 935 (2005).
  * Buzdin  _et al._ [1992] A. I. Buzdin, B. Bujicic, and B. M. Y. Kupriyanov, Sov. Phys. JETP 74, 124 (1992).
  * Buzdin [2008] A. Buzdin, Phys. Rev. Lett. 101, 107005 (2008).
  * Bergeret and Tokatly [2015] F. S. Bergeret and I. V. Tokatly, Europhys. Lett. 110, 57005 (2015).
  * Bezuglyi  _et al._ [2002] E. V. Bezuglyi, A. S. Rozhavsky, I. D. Vagner, and P. Wyder, Phys. Rev. B 66, 052508 (2002).
  * Barash and Bobkova [2002] Y. S. Barash and I. V. Bobkova, Phys. Rev. B 65, 144502 (2002).
  * Dimitrova and Feigel’man [2006] O. V. Dimitrova and M. V. Feigel’man, J. Exp. Theor. Phys. 102, 652 (2006).
  * Gogin  _et al._ [2022] L. Gogin, F. Rossi, and F. Dolcini, New J. Phys. 24, 093025 (2022).
  * Liu and Chan [2010] J.-F. Liu and K. S. Chan, Phys. Rev. B 82, 125305 (2010).
  * Reynoso  _et al._ [2008] A. A. Reynoso, G. Usaj, C. A. Balseiro, D. Feinberg, and M. Avignon, Phys. Rev. Lett. 101, 107001 (2008).
  * Reynoso  _et al._ [2012] A. A. Reynoso, G. Usaj, C. A. Balseiro, D. Feinberg, and M. Avignon, Phys. Rev. B 86, 214519 (2012).
  * Ness  _et al._ [2022] H. Ness, I. A. Sadovskyy, A. E. Antipov, M. van Schilfgaarde, and R. M. Lutchyn, npj Comp. Mater. 8, 23 (2022).
  * Krive  _et al._ [2004] I. V. Krive, L. Y. Gorelik, R. I. Shekhter, and M. Jonson, Low Temp. Phys. 30, 398 (2004).
  * Cheng and Lutchyn [2012] M. Cheng and R. M. Lutchyn, Phys. Rev. B 86, 134522 (2012).
  * Yokoyama  _et al._ [2013] T. Yokoyama, M. Eto, and Y. V. Nazarov, J. Phys. Soc. Japan 82, 054703 (2013).
  * Yokoyama and Nazarov [2014] T. Yokoyama and Y. V. Nazarov, Europhys. Lett. 108, 47009 (2014).
  * Yokoyama  _et al._ [2014] T. Yokoyama, M. Eto, and Y. V. Nazarov, Phys. Rev. B 89, 195407 (2014).
  * Campagnano  _et al._ [2015] G. Campagnano, P. Lucignano, D. Giuliano, and A. Tagliacozzo, J. Phys: Cond. Matt. 27, 205301 (2015).
  * Doh  _et al._ [2005] Y.-J. Doh, J. A. van Dam, A. L. Roest, E. P. A. M. Bakkers, L. P. Kouwenhoven, and S. D. Franceschi, Science 309, 272 (2005).
  * Fulde and Ferrell [1964] P. Fulde and R. A. Ferrell, Phys. Rev. 135, A550 (1964).
  * [60] A. I. Larkin and Y. N. Ovchinnikov, Zh. Eksp. Teor. Fiz. 47 1136 (1964), Sov. Phys. JETP 20 762 (1965) .
  * Buzdin  _et al._ [1982] A. I. Buzdin, L. N. Bulaevskii, and S. V. Panjukov, JETP Lett. 35, 178 (1982).
  * Ryazanov  _et al._ [2001] V. V. Ryazanov, V. A. Oboznov, A. Y. Rusanov, A. V. Veretennikov, A. A. Golubov, and J. Aarts, Phys. Rev. Lett. 86, 2427 (2001).
  * Chtchelkatchev and Nazarov [2003] N. M. Chtchelkatchev and Y. V. Nazarov, Phys. Rev. Lett. 90, 226806 (2003).
  * Béri  _et al._ [2008] B. Béri, J. H. Bardarson, and C. W. J. Beenakker, Phys. Rev. B 77, 045311 (2008).