$\displaystyle=\frac{1}{\pi}\sqrt{2\,\frac{N}{d_{s}(0)}}+\frac{1}{4}.$

###  .8 Growth of perturbations in rate networks

On small timescales the distance evolves in proportion to the Lyapunov
exponent

$\displaystyle d(t)$ $\displaystyle=d(0)\,\exp(\lambda_{\mathrm{max}}\,t).$

So
$d_{s}(t)-d_{n}(t)=\big{(}d_{s}(0)-d_{n}(0)\big{)}\,\exp\big{(}\lambda_{\mathrm{max}}\,t\big{)}$.
In particular, there is no maximum expected on a timescale of the neuronal
dynamics. The typical timescale instead is determined by the maximal Lyapunov
exponent.

To derive an equation for the time evolution of the correlation between
replicas in the rate network, we use the pair of equations obtained from the
mean-field description

$\displaystyle\tau\,\partial_{t}x^{\alpha}(t)$
$\displaystyle=-x^{\alpha}(t)+h^{\alpha}(t)$ (59) $\displaystyle\langle
h^{\alpha}(t)h^{\beta}(s)\rangle$ $\displaystyle=Q^{\alpha\beta}(t,s)\equiv
g^{2}\,\langle T(x^{\alpha}(t))T(x^{\beta}(s))\rangle.$

We may thus write

$\displaystyle x^{\alpha}(t)$
$\displaystyle=\frac{1}{\tau}\,\int_{-\infty}^{t}\,e^{-\frac{t-t^{\prime}}{\tau}}\,h^{\alpha}(t^{\prime})\,dt^{\prime}.$

So the correlation function $c^{\alpha\beta}(t,s):=\langle
x^{\alpha}(t)x^{\beta}(s)\rangle$ obeys

$\displaystyle c^{\alpha\beta}(t,s)$
$\displaystyle=\frac{1}{\tau}\,\int_{-\infty}^{t}\,e^{-\frac{t-t^{\prime}}{\tau}}\,\frac{1}{\tau}\,\int_{-\infty}^{s}\,e^{-\frac{s-s^{\prime}}{\tau}}\,\langle
h^{\alpha}(t^{\prime})h^{\beta}(s^{\prime})\rangle\,dt^{\prime}\,ds^{\prime},$

which becomes in differential form

$\displaystyle\tau\partial_{t}c^{\alpha\beta}(t,s)$
$\displaystyle=-c^{\alpha\beta}(t,s)+\frac{1}{\tau}\,\int_{-\infty}^{s}\,e^{-\frac{s-s^{\prime}}{\tau}}\,\langle
h^{\alpha}(t)h^{\beta}(s^{\prime})\rangle\,ds^{\prime}.$

This differential equation allows the integration along the $t$ direction by
one time-step $\delta$

$\displaystyle c^{\alpha\beta}(t+\delta,s)$
$\displaystyle=(1-\frac{\delta}{\tau})\,c^{\alpha\beta}(t,s)+\frac{\delta}{\tau}\,\int_{-\infty}^{s}\,e^{-\frac{s-s^{\prime}}{\tau}}\,\langle
h^{\alpha}(t)h^{\beta}(s^{\prime})\rangle\,ds^{\prime},$

which requires only $c^{\alpha\beta}$ and $\langle hh\rangle$ in the $t-$past
and in the $s-$past. The integration can be done for $t\in[0,T]$, where $T$ is
a desired final point.

As a result, one has $c^{\alpha\beta}(t,s)$ for $t\in[0,T]$. In the next
update step we move into the $s-$direction by

$\displaystyle c^{\alpha\beta}(t,s+\delta)$
$\displaystyle=(1-\frac{\delta}{\tau})\,c^{\alpha\beta}(t,s)+\frac{\delta}{\tau}\,\int_{-\infty}^{t}\,e^{-\frac{t-t^{\prime}}{\tau}}\,\langle
h^{\alpha}(t^{\prime})h^{\beta}(s)\rangle\,dt^{\prime},$

where the latter integral can be computed because it requires only

$\displaystyle\langle h^{\alpha}(t)h^{\beta}(s)\rangle$
$\displaystyle=g^{2}\langle
T(x_{1})T(x_{2})\rangle_{(x_{1},x_{2})\sim\mathcal{N}(0,c^{\alpha\beta}(t,s))}$
(60)

computed in the previous step.

It makes sense to introduce as an auxiliary variable

$\displaystyle q^{\alpha\beta}(t,s)$
$\displaystyle:=\frac{1}{\tau}\int_{-\infty}^{t}\,e^{-\frac{t-t^{\prime}}{\tau}}\,\langle
h^{\alpha}(t^{\prime})h^{\beta}(s)\rangle\,dt^{\prime},$

for which we can assume the symmetry
$q^{\alpha\beta}(t,s)=q^{\beta\alpha}(t,s)$ to write the updates

$\displaystyle c^{\alpha\beta}(t+\delta,s)$
$\displaystyle=(1-\frac{\delta}{\tau})\,c^{\alpha\beta}(t,s)+\frac{\delta}{\tau}\,q^{\alpha\beta}(s,t),$
(61) $\displaystyle c^{\alpha\beta}(t,s+\delta)$
$\displaystyle=(1-\frac{\delta}{\tau})\,c^{\alpha\beta}(t,s)+\frac{\delta}{\tau}\,q^{\alpha\beta}(t,s).$
(62)

The auxiliary variable $q^{\alpha\beta}(t,s)$ obeys the differential equation
$(\tau\partial_{t}+1)\,q^{\alpha\beta}(t,s)=\langle
h^{\alpha}(t)h^{\beta}(s)\rangle$, which yields the update equation

$\displaystyle q^{\alpha\beta}(t+\delta,s)$
$\displaystyle=(1-\frac{\delta}{\tau})\,q^{\alpha\beta}(t,s)+\frac{\delta}{\tau}\,\langle
h^{\alpha}(t)h^{\beta}(s)\rangle.$ (63)

So the required sequence of updates is:

  1. 1.

Start at $t=0$.

  2. 2.

Assume we have computed $c^{\alpha\beta}(t^{\prime},s^{\prime})$ and
$q^{\alpha\beta}(t^{\prime},s^{\prime})$ until this point $t$ for all
$(t^{\prime}<t,s^{\prime}<t)$.

  3. 3.

Compute $c^{\alpha\beta}(t+\delta,\forall s^{\prime}\leq t)$ using (61).

  4. 4.

Compute $c^{\alpha\beta}(\forall t^{\prime}\leq t,t+\delta)$ using (62).

  5. 5.

Compute $\langle h^{\alpha}(t+\delta)h^{\beta}(s^{\prime})\rangle$ $\forall
s^{\prime}\leq t$ and $\langle
h^{\alpha}(t^{\prime})h^{\beta}(t+\delta)\rangle$ $\forall t^{\prime}\leq t$
using (60) and the result from the previous step.

  6. 6.

Compute $q^{\alpha\beta}(t^{\prime}\leq t+\delta,t+\delta)$ by iterating (63)
with zero initial condition and starting with $t^{\prime}$ sufficiently far
back in the past.

  7. 7.

Compute $c^{\alpha\beta}(t+\delta,t+\delta)$ by average of (61) and (62) using
new value $q^{\alpha\beta}(t+\delta,t)$ and $q^{\alpha\beta}(t,t+\delta)$,
respectively.

  8. 8.

Go to the next time slice $t\to t+\delta$, return to step 3.

###  .9 Flux tubes in binary networks

It has been shown by Puelma Touzel and Wolf (2019) that the borders of flux
tubes in spiking networks of inhibitory LIF neurons are related to changes in
the global order of spikes. In particular, if a perturbation creates an
additional spike or causes the omission of an expected one, the mean firing
rate will stay constant but the order of future spikes is very likely to be
irrevocably changed. The divergence rate of two trajectories can be assessed
by calculating the mean number of unexpected spike order changes caused by a
single such perturbation, resulting in a branching process. In the context of
binary networks, we can ask the equivalent question: Given a flip of a single
neuron’s activity variable, how many “wrong” update results will occur on
average in the following time $\tau$ ?

The flip of one neuron $x_{j}\to-x_{j}$ causes a change $\Delta
h_{i}=-2J_{ij}x_{j}$ in the input of neurons it is connected to. Across
different target neurons $\Delta h$ is therefore distributed as

$\rho(\Delta h)=\mathcal{N}\left(\pm
2\frac{\bar{g}}{N},4\frac{g^{2}}{N}\right)$

and the probability of a neuron to be updated into the wrong state due to the
perturbation in the input is

$\displaystyle p\left(\text{flip}(x)|\Delta h\right)$
$\displaystyle=\left\langle\left|\mathrm{T}_{\mathrm{p}}(h+\Delta
h)-\mathrm{T}_{p}(h)\right|\right\rangle_{h\sim\mathcal{N}(R,Q_{0})}$
$\displaystyle\overset{\left|\Delta h\right|\ll
1,\,\mathrm{T}_{\mathrm{p}}^{\prime}\geq 0}{\approx}\left|\Delta
h\right|\left\langle\mathrm{T}_{p}^{\prime}(h)\right\rangle_{h\sim\mathcal{N}(R,Q_{0})},$

where the absolute value enters because both directions of perturbation cause
a positive probability of “wrong” updating, and we assume $T^{\prime}\geq 0$
for simplicity. Now we ask the following question: How many downstream flips
$n_{\text{spawns}}$ will, on average, be triggered in the network during one
time constant, given a single original flip? This quantity controls whether
the decorrelating flips will proliferate or not, because since every neuron is
updated on average once per time constant, if $n_{\text{spawns}}<1$ and the
neuron carrying the original flip is updated again, it is most likely updated
“correctly” again and the average number of flips in the network has
decreased. If $n_{\text{spawns}}>1$ on the other hand, the average number of
flips increases.

Being interested in the transition point, we can assume
$n_{\text{spawns}}\approx 1$ so that we do not need to take the interaction of
several flips into account. Then

$\displaystyle\bar{n}_{\text{spawns}}$ $\displaystyle=N\left\langle
p\left(\text{flip}(x)|\Delta h\right)\right\rangle_{\Delta h}$ (64)

and while we take the mean input $R$ into account, we neglect the perturbation
of the mean input $\left\langle\Delta h\right\rangle=\pm
2\bar{g}/N=\text{$\mathcal{O}(N^{-1})$}\approx 0$ as it is small compared to
the standard deviation $\sigma_{\Delta
h}=2g/\sqrt{N}=\mathcal{O}(N^{-\frac{1}{2}})$, allowing the simple calculation

$\displaystyle\bar{n}_{\text{spawns}}$
$\displaystyle=N\left\langle\mathrm{T}_{p}^{\prime}(h)\right\rangle_{h\sim\mathcal{N}(R,Q_{0})}\left\langle\left|\Delta
h\right|\right\rangle_{\Delta h\sim\mathcal{N}(0,\sigma_{\Delta h}^{2})}$
$\displaystyle=N\left\langle\mathrm{T}_{p}^{\prime}(h)\right\rangle_{h\sim\mathcal{N}(R,Q_{0})}\sqrt{\frac{2}{\pi}}\sigma_{\Delta
h}$
$\displaystyle=2\sqrt{\frac{2N}{\pi}}g\left\langle\mathrm{T}_{p}^{\prime}(h)\right\rangle_{h\sim\mathcal{N}(R,Q_{0})}.$

Finally, accounting for $\mathrm{T}_{p}^{\prime}(h)=\mathrm{T}^{\prime}(h)/2$,
given by (6), the chaos transition is expected at

$\displaystyle 1\overset{!}{=}\bar{n}_{\text{spawns}}$
$\displaystyle=\frac{\sqrt{2N}}{\sqrt{\pi}}g\left\langle\mathrm{T}^{\prime}(h)\right\rangle_{h\sim\mathcal{N}(R,Q_{0})},$
(65)

which is exactly the result (15), derived via the completely different route
of the replica calculation. While the derivation here is nicely and
intuitively interpretable, the derivation via field theory and replica
calculation allows for systematic generalizations. For example, it is not
clear how to obtain the residual correlation (13) in the ad hoc approach.

#### Fluxtube size

The flux-tube diameter is not a very informative measure for a binary network,
since the system trajectory in phase space is typically not in the middle of a
“tube” but close to some of its boundaries (given by the thresholds).
Therefore, the distance to a boundary strongly depends on the direction of
perturbation. As a relatively informative measure, we consider smearing the
trajectory in all directions with some variance $\text{Var}(\Delta
h)=\sigma_{\text{fl}}^{2}$, which is chosen such that on average, one flux-
tube boundary is crossed. This procedure makes sense insofar, as it is similar
to adding noise onto the input. It is important to be aware that
$\sigma_{\text{fl}}$ is not strictly the average distance to the closest
boundary, although the two quantities should covary.

The situation is analogous to the above calculation, because we again need to
consider the flips occurring during an update in (on average) all $N$ neurons,
which is given by (64) only with $\sigma_{\Delta h}$ replaced by
$\sigma_{\text{fl}}$. Demanding $\bar{n}_{\text{spawns}}\overset{!}{=}1$ then
yields

$\displaystyle 1$ $\displaystyle\overset{!}{=}N\left\langle
p\left(\text{flip}(x)|\Delta h\right)\right\rangle_{\Delta
h\sim\mathcal{N}(0,\sigma_{\text{fl}}^{2})}$
$\displaystyle\Rightarrow\quad\sigma_{\text{fl}}$
$\displaystyle=\frac{\sqrt{2\pi}}{N\left\langle\mathrm{T}^{\prime}(h)\right\rangle_{h\sim\mathcal{N}(R,Q_{0})}}.$

Of course, the $1/N$ scaling needs to be taken with caution, since our
perturbation goes into all $N$ phase-space directions, resulting in a total
length scaling as $1/\sqrt{N}$ .

###  .10 Equivalence of dynamical mean-field theories of binary and rate
networks

The dynamics (21) can equivalently be written as

$\displaystyle\tau\partial_{t}\boldsymbol{x}$
$\displaystyle=-\boldsymbol{x}+\mathrm{T}\Bigl{(}\boldsymbol{h}\Bigr{)},$ (66)
$\displaystyle\boldsymbol{h}$
$\displaystyle=\boldsymbol{Jx}+\sqrt{\tau}\boldsymbol{\eta},$

where the noise $\eta_{i}$ is an Ornstein-Uhlenbeck process (Uhlenbeck and
Ornstein, 1930), $(\tau\partial_{t}+1)\,\eta_{i}=\xi_{i}$. This form allows
the application of the model-independent field theory. The single-neuron,
single-time-slice probability functional is
$\rho[x(t)|h]=\delta[x(t)-\int_{-\infty}^{t}\,e^{-\frac{t-t^{\prime}}{\tau}}\,\mathrm{T}\left(h_{i}\left(t^{\prime}\right)\right)\,\frac{dt^{\prime}}{\tau}]$,
and the noise term is taken into account in
$\rho[\boldsymbol{h}]=\delta[\boldsymbol{h}-\boldsymbol{J}\boldsymbol{x}-\sqrt{\tau}\boldsymbol{\eta}]$.
Plugging this expression into (2) we obtain the same equation (5) for the mean
activity as for the binary neuron, if we choose the strength of the noise
$\sigma_{\xi}$ such that $Q(t,t)=g^{2}$ as well. The reason for the
equivalence is that the exponential function appearing in the convolution
equation is the Green’s function of $\tau\partial_{t}+1$. In a stationary
state, the saddle-point solution for $Q(\Delta t)\coloneqq Q(t,t+\Delta
t)=g^{2}\,\langle x(t)x(t+\tau)\rangle$, moreover, follows the same Newtonian
equation of motion (9) as for the binary model (Sompolinsky  _et al._ , 1988,
eq. 7).

#### Matching initial conditions

Knowing that the differential equations for the time-lagged autocorrelations
are the same, we have to adjust their respective initial conditions to
establish full equivalence. Here we use the subscripts $b$ and $r$ to refer to
the quantities of the binary and rate model, respectively. Two initial
conditions are needed for a unique solution. One is to require that
$\lim_{t\rightarrow\infty}\dot{Q}\left(t\right)=0$, which is the same in both
cases. So the autocorrelation for infinite time-lags is described by a single
value $Q_{\infty}$, which vanishes for point-symmetric activation functions,
but is in general nonzero and self-consistently determined by the static
variability across neurons, caused by the disorder (compare Figure 5). In the
binary case, the second condition is $Q_{b}(0)=g^{2}$ because the zero-lag
autocorrelation of a single spin is always one. In a rate network, however,
the input noise strength determines how quickly the autocorrelation decays,
resulting in the condition on the derivative
$\dot{Q}_{r}(0+)=-\sigma_{\xi}^{2}/2$ (Schuecker  _et al._ , 2018)111(note
that our notational convention differs by the factor $1/2$).

The idea is to choose the variance of the noise $\sigma_{\xi}^{2}$ in the rate
network such that $Q_{r}(0)=Q_{b}(0)$, so that the time-lagged solutions for
the variance $Q(\Delta t)$ match.

To do so, using $\dot{Q}(\infty)=0$ and conservation of total “energy”
$V_{Q_{0}}+\tau^{2}\dot{Q}^{2}/2$ implied by the Newtonian form of (10), the
condition $Q_{0}:=Q_{r}(0)\stackrel{{\scriptstyle!}}{{=}}Q_{b}(0)=g^{2}$ can
be expressed as a condition for the derivative and thus the noise amplitude

$\displaystyle\frac{1}{2}\tau^{2}\dot{Q}_{0}^{2}+V_{Q_{0}}(Q_{0})$
$\displaystyle=\left.V_{Q_{0}}(Q_{\infty})\right|_{Q_{0}=g^{2}}.$ (67)

Plugging in $\dot{Q}_{0}=-\sigma_{\xi}^{2}/2$ and solving for
$\sigma_{\xi}^{2}$ yields the condition (22) in the main text. This proves
that the binary and rate model with appropriate noise have equivalent mean
activities and time-lagged autocorrelations in dynamical mean-field
approximation.

#### Explanation of the result

Figure 14: Equivalence of binary (left) and rate models in the two different
forms (66) (middle) and (21) (right): Mapping from output to input by
identical matrix $J$; asynchronous update process $U$ with rate $\tau^{-1}$
implies exponential convolution kernel, leading to leaky-integration (cf.
(5)), identical to operator $L=(\tau\partial_{t}+1)^{-1}$ present explicitly
in (66) and (21). Transitions between discrete binary states effect red noise
$\eta$ in input $h$ (middle), which corresponds to white noise $\xi$ that is
low-pass filtered by $L$ (right). Rate models differ in the order of
application of this kernel and the connectivity, which yields equivalent
dynamics because the two operators commute.

Taking a step back, what is the intuition behind this result? When the binary
neurons are averaged over realizations of the update time disorder, the
Poisson update process with rate $\tau^{-1}$ becomes an exponential kernel
corresponding to that of the rate network. The discrete jumps of the binary
neurons around their mean become red noise (Lindner, 2009; Frey and Weber,
2017), corresponding to the low-pass-filtered noise $\eta$ of the rate network
(66). By nice conspiracy, this noise corresponds to simple white noise $\xi$
in (21), which is also the version treated in most works on rate networks with
noise, such as (Rajan  _et al._ , 2010; Aljadeff  _et al._ , 2015; Kadmon and
Sompolinsky, 2015; Martí  _et al._ , 2018; Schuecker  _et al._ , 2018;
Crisanti and Sompolinsky, 2018). This tight relation between the binary and
rate models is summarized conceptually in Figure 14.

###  .11 Slope of correlation transmission in binary and rate neurons

Here we show that the difference between discrete signaling and continuous
signaling leads to a qualitative difference in the slope of the correlation-
transmission curve and thus the transition to chaos.

Assume, as an approximation, that two neurons receive inputs that are jointly
Gaussian distributed as

$\displaystyle(h_{1},h_{2})$ $\displaystyle\sim\mathcal{N}(0,K),$

where the covariance matrix is given by

$\displaystyle K(c_{\mathrm{in}})$
$\displaystyle:=q\,\left(\begin{array}[]{cc}1&c_{\mathrm{in}}\\\
c_{\mathrm{in}}&1\end{array}\right).$

Here $c_{\mathrm{in}}\in[-1,1]$ controls the correlation between the inputs.

#### Continuous signaling

A neuron with continuous signaling has the output

$\displaystyle y_{i}$ $\displaystyle=\mathrm{T}(h_{i}),$

where $\mathrm{T}\in[-1,1]$ is an activation function. The mean output is thus

$\displaystyle\langle y_{i}\rangle$
$\displaystyle=\langle\mathrm{T}(h)\rangle_{h\sim\mathcal{N}(0,q)}.$ (68)

For a point-symmetric gain function that we assume in the following the mean
vanishes so that the variance of the outputs is

$\displaystyle a$
$\displaystyle=\langle\mathrm{T}^{2}(h)\rangle_{h\sim\mathcal{N}(0,q)}.$

The correlation coefficient between the outputs of a pair of neurons is

$\displaystyle c_{\mathrm{out}}^{\text{cont.}}(c_{\mathrm{in}})$
$\displaystyle:=a^{-1}\,\langle y_{1}y_{2}\rangle$ (69)
$\displaystyle=a^{-1}\,\langle\mathrm{T}(h_{1})\mathrm{T}(h_{2})\rangle_{(h_{1},h_{2})\sim\mathcal{N}(0,K(c_{\mathrm{in}}))}$

which has the slope

$\displaystyle\frac{dc_{\mathrm{out}}^{\mathrm{cont.}}(c_{\mathrm{in}})}{dc_{\mathrm{in}}}$
$\displaystyle=a^{-1}\,\langle\mathrm{T}^{\prime}(h_{1})\mathrm{T}^{\prime}(h_{2})\rangle_{(h_{1},h_{2})\sim\mathcal{N}(0,K(c_{\mathrm{in}}))},$

by Price’s theorem (Papoulis and Pillai, 2002). Evaluated at
$c_{\mathrm{in}}=1$ this is

$\displaystyle\frac{dc_{\mathrm{out}}^{\mathrm{cont.}}(1)}{dc_{\mathrm{in}}}$
$\displaystyle=a^{-1}\,\langle(\mathrm{T}^{\prime}(h))^{2}\rangle_{h\sim\mathcal{N}(0,q)}\stackrel{{\scriptstyle\mathrm{T}^{\prime}<\infty}}{{<}}\infty.$
(70)

For activation functions $\mathrm{T}$ with finite slope
$\mathrm{T}^{\prime}<\infty$ this slope is thus finite. For the signum
function $\mathrm{T}(x)=2H(x)-1$ we get $a=1$ and

$\displaystyle\frac{dc_{\mathrm{out}}^{\mathrm{cont.}}(c_{\mathrm{in}})}{dc_{\mathrm{in}}}$
$\displaystyle=2\,\langle\delta(h_{1})\delta(h_{2})\rangle_{(h_{1},h_{2})\sim\mathcal{N}(0,K(c_{\mathrm{in}}))}$
$\displaystyle=\frac{1}{\pi\,\sqrt{\det(K(c_{\text{in}}))}}=(q\pi)^{-1}\,(1-c_{\mathrm{in}}^{2})^{-\frac{1}{2}},$
(71)

where the latter line comes from the normalization condition of the two-
dimensional Gaussian distribution. Thus, the slope diverges if and only if the
output of the neuron becomes discrete.

#### Discrete signaling

Now consider a neuron with discrete output, but smooth activation function
$T\in[-1,1]$; a smooth function here corresponds to a probabilistic activation

$\displaystyle y_{i}$ $\displaystyle=\begin{cases}1&\text{with prob.
}\big{(}\mathrm{T}(h)+1\big{)}/2\\\ -1&\text{with prob.
}1-\big{(}\mathrm{T}(h)+1\big{)}/2\end{cases}.$

The mean output is thus

$\displaystyle\langle y_{i}\rangle$ $\displaystyle=\langle
1\cdot\big{(}\mathrm{T}(h)+1\big{)}/2-1\cdot(1-\big{(}\mathrm{T}(h)+1\big{)}/2)\rangle_{h\sim\mathcal{N}(0,q)}$
$\displaystyle=\langle\mathrm{T}(h)\rangle_{h\sim\mathcal{N}(0,q)},$

the same as for the continuous signaling (68). For a point-symmetric gain
function that we assume in the following the mean vanishes so that the
variance of the outputs is $a=1$. The correlation coefficient of the output is
then identical to the second moment between the outputs of a pair of neurons

$\displaystyle c_{\mathrm{out}}^{\text{disc}}(c_{\mathrm{in}})$
$\displaystyle:=\langle y_{1}y_{2}\rangle$ (72)
$\displaystyle=1\cdot(1-p_{\mathrm{diff}})-1\cdot
p_{\mathrm{diff}}=1-2\,p_{\mathrm{diff}}$
$\displaystyle=1-\langle|\mathrm{T}(h_{1})-\mathrm{T}(h_{2})|\rangle_{(h_{1},h_{2})\sim\mathcal{N}(0,K(c_{\mathrm{in}}))}.$

The latter expression is related to the probability
$p_{\mathrm{diff}}=\frac{1}{2}\langle|\mathrm{T}(h_{1})-\mathrm{T}(h_{2})|\rangle$
that the two neurons are in different states. This expression is of course the
same as found in (55). In the limit of $c_{\mathrm{in}}\to 1$ we thus have

$\displaystyle c_{\mathrm{out}}^{\text{disc}}(c_{\mathrm{in}})$
$\displaystyle\stackrel{{\scriptstyle
c_{\mathrm{in}}=1-\hat{\epsilon}}}{{\simeq}}1-2\,\left\langle\mathrm{T}^{\prime}\left(h\right)\right\rangle_{h\sim\mathcal{N}\left(0,q\right)}\,\sqrt{\frac{\hat{\epsilon}}{\pi}}+\mathcal{O}\left(\hat{\epsilon}^{3/2}\right)$

where $\hat{\epsilon}=\frac{\epsilon}{g^{2}}$. So the slope diverges for
$c_{\mathrm{in}}\to 1$ as

$\displaystyle\frac{dc_{\mathrm{out}}^{\mathrm{disc}}(c_{\mathrm{in}})}{dc_{\mathrm{in}}}\big{|}_{c_{\mathrm{in}}=1-\hat{\epsilon}}$
$\displaystyle\stackrel{{\scriptstyle(\ref{eq:abs_calc_Ising-1-1})}}{{=}}\frac{d}{d(-\hat{\epsilon})}\Big{[}1-2\,\left\langle\mathrm{T}^{\prime}\left(h\right)\right\rangle_{h\sim\mathcal{N}\left(0,q\right)}\,\sqrt{\frac{\hat{\epsilon}}{\pi}}+\mathcal{O}\left(\hat{\epsilon}^{3/2}\right)\Big{]}$
$\displaystyle=\left\langle\mathrm{T}^{\prime}\left(h\right)\right\rangle_{h\sim\mathcal{N}\left(0,q\right)}\,(\pi\hat{\epsilon})^{-\frac{1}{2}}+\mathcal{O}\left(\hat{\epsilon}^{1/2}\right)\,\propto(1-c_{\mathrm{in}})^{-\frac{1}{2}}.$

This divergence is present even if the gain function has a finite slope
$\mathrm{T}^{\prime}<\infty$. This is in qualitative contrast to the finite
slope found for the continuous signaling in (70).

The infinite slope for continuous signaling in the limit of a sharp activation
function, (71), can be shown to have the same form of divergence for
$c_{\mathrm{in}}\to 1$ when expanded for $c_{\text{in}}=1-\frac{1}{q}\epsilon$
in the limit of small $\epsilon\ll 1$.

###  .12 Noisy binary pattern classification task

We implement a classification task by training one linear readout

$\displaystyle S_{\alpha^{\prime}}(t)$
$\displaystyle=w_{\alpha^{\prime}}(t)^{\mathrm{T}}(x_{\alpha}(t)+\xi_{\text{pre}})+\xi_{\text{post}}$
(73)

of the network state $x_{\alpha}(t)$ at time $t$ for each of the
$\alpha^{\prime}=1,\ldots,P=50$ patterns to be detected. Here
$\xi_{\text{pre}}$ and $\xi_{\text{post}}$ are additional Gaussian readout
noises of standard deviation $\sigma_{\xi,\text{pre}}$ and
$\sigma_{\xi,\text{post}}$, respectively. $\xi_{\text{pre}}$ controls how
precisely a single neuron’s state can be read out. $\xi_{\text{post}}$
represents a noise component of the classification mechanism. Training of the
readout $w_{\alpha^{\prime}}(t)$ is performed for each time point $t$ by
linear regression (see Appendix .12.1), minimizing the quadratic error of
detecting the stimulus identity, i.e. minimizing
$(S_{\alpha^{\prime}}-\delta_{\alpha\alpha^{\prime}})^{2}$.

The patterns are presented to the network by initializing the first $L=10$ of
the $N=500$ neurons to the stimulus. All other neurons are in an initial state
corresponding to the stationary statistics. Each stimulus $\alpha$ is a random
binary pattern of length $L$ with $\\{-1,1\\}$ appearing equally likely,
superimposed with Gaussian noise of standard deviation $\sigma$. Note that
because of the noise added to the binary values, these initial states are not
strictly $\in\\{-1,1\\}$. This freedom in the initial states is just a way to
introduce the noise; after their first update the neurons’ states are strictly
$\\{-1,1\\}$ again. The resulting evolution of the network state given this
initial condition is termed $x_{\alpha}(t)$.

####  .12.1 Linear regression

Minimizing the quadratic error over all patterns amounts to linear regression;
we consider a single scalar readout target value $y_{\alpha}\in\mathbb{R}$ for
each pattern $\alpha$; in the example above $y_{\alpha}\in\\{0,1\\}$. Then
$w,x\in\mathbb{R}^{N}$ and the quadratic error is

$\displaystyle\epsilon$
$\displaystyle:=\min_{w}\sum_{\alpha=1}^{P}\big{(}w^{\mathrm{T}}x_{\alpha}-y_{\alpha}\big{)}^{2}.$
(74)

Demanding stationarity with regard to $w$ by differentiating by
$\partial_{w_{i}}$ we get $N$ equations

$\displaystyle 0$
$\displaystyle=\sum_{\alpha=1}^{P}2\,\big{(}w^{\mathrm{T}}x_{\alpha}-y_{\alpha}\big{)}\,x_{\alpha
i}\quad\forall i$ $\displaystyle
w^{\mathrm{T}}\sum_{\alpha=1}^{P}x_{\alpha}x_{\alpha}^{\mathrm{T}}$
$\displaystyle=\sum_{\alpha=1}^{P}y_{\alpha}x_{\alpha}^{\mathrm{T}}.$

The value $w^{\ast}$ to achieve stationarity is

$\displaystyle w^{\ast}$
$\displaystyle=C^{-1}\,\sum_{\alpha=1}^{P}y_{\alpha}x_{\alpha}^{\mathrm{T}},$
(75) $\displaystyle\text{with }C$
$\displaystyle:=\sum_{\alpha=1}^{P}x_{\alpha}x_{\alpha}^{\mathrm{T}},$

where we use the symmetry of $C$. Inserted into (74)

$\displaystyle\epsilon$
$\displaystyle=w^{\ast\mathrm{T}}\,C\,w^{\ast}-2\,w^{\ast\mathrm{T}}\sum_{\alpha}y_{\alpha}x_{\alpha}+\sum_{\alpha}y_{\alpha}^{2}$
$\displaystyle=\sum_{\alpha=1}^{P}y_{\alpha}^{2}-S,$ $\displaystyle S$
$\displaystyle=\big{(}\sum_{\alpha=1}^{P}y_{\alpha}x_{\alpha}^{\mathrm{T}}\big{)}\,C^{-1}\,\big{(}\sum_{\alpha=1}^{P}y_{\alpha}x_{\alpha}\big{)}.$
(76)

In the case of classification, the latter expression simplifies even further:
The first term is a constant $\sum_{\alpha=1}^{P}y_{\alpha}^{2}=1$ for labels
$y_{\alpha}\in\\{0,1\\}$, where $y_{\alpha}=1$ if the presented pattern
$\alpha$ is the pattern $\alpha^{\prime}$ to be detected and $y_{\alpha}=0$
else, $y_{\alpha}=\delta_{\alpha\alpha^{\prime}}$. The second term is then
identical to the definition (73) for $\xi=0$, obtained by inserting $w^{\ast}$
from (75). The expression shows that the signal amplitude
$S_{\alpha^{\prime}}$ actually depends on the signal-to-noise ratio, the
length of the selected vector $x_{\alpha^{\prime}}$ measured with regard to
the variability $C$ across all patterns

$\displaystyle S_{\alpha^{\prime}}$
$\displaystyle=x_{\alpha^{\prime}}^{\mathrm{T}}\,C^{-1}\,x_{\alpha^{\prime}}.$
(77)

The generalization to stochastic realizations of $x_{\alpha}$, for example due
to the presentation of noisy patterns is straightforward. We need to replace
$\sum_{\alpha=1}^{P}\ldots$ by $\sum_{\alpha=1}^{P}\langle\ldots\rangle$ in
the measure for the error (74) and thus throughout this calculation, where
$\langle\ldots\rangle$ is the expectation over the noise realizations.

####  .12.2 Approximation of orthogonal patterns and uniform noise

If the patterns $x$ are sufficiently orthogonal in the signal subspace, we can
think of the entries $k$ of any state vector $x_{\alpha}$ to be drawn
independently. So the $x_{\alpha,k}\in\\{-1,1\\}$ appear with equal
probability for those entries that lie in the subspace of dimension
$d_{\mathrm{s}}$. All remaining entries are assumed to be constant across
patterns. We may thus restrict the space to the $d_{\mathrm{s}}$ informative
components. The $kl$-th element of the covariance matrix for independently
drawn entries is

$\displaystyle C_{kl}$ $\displaystyle=\sum_{\alpha=1}^{P}x_{\alpha k}x_{\alpha
l}$ $\displaystyle\simeq P\delta_{kl}.$ (78)

The signal of the readout $\alpha^{\prime}$, following from (77), then takes
the simple form

$\displaystyle S_{\alpha^{\prime}}$ $\displaystyle\simeq
P^{-1}\,||x_{\alpha^{\prime}}||_{d_{s}}^{2}.$ (79)

If the signal is perfectly reliable, that is, if for all noise realizations
$i$ the response $x_{\alpha i}$ is equal to the stereotypical response
$x_{\alpha i}=\bar{x}_{\alpha}$, and if the dimension of the informative
subspace is $d_{s}$, so $\bar{x}\in\mathbb{R}^{d_{s}}$, we get with
$||\bar{x}||_{d_{s}}^{2}=d_{\mathrm{s}}$

$\displaystyle S_{\alpha^{\prime},\mathrm{max}}$
$\displaystyle\apprle\frac{d_{\mathrm{s}}}{P}.$

If noisy realizations of patterns cause flips in random entries of
$x_{\alpha^{\prime}}$, which is an approximation since spins are expected to
differ in their susceptibility, the responses are not perfectly reliable, so
we need to replace $x_{\alpha^{\prime}}$ by $\langle
x_{\alpha^{\prime}}\rangle$ in (77) and thus

$\displaystyle S_{\alpha^{\prime}}$ $\displaystyle\simeq P^{-1}\,||\langle
x_{\alpha^{\prime}}\rangle||_{d_{s}}^{2}.$

The above expressions thereby link the readout signal to the dimensionality of
the responses, as discussed in the main text in Section II.6.

####  .12.3 Nonzero plateau of the signal.

In the simulations, the noise distance somewhat unintuitively saturates
slightly below the signal distance. This is explainable by taking into account
that not all the initial noise realizations actually cause a crossing of the
flux-tube boundary. Instead, those realizations simply follow the unperturbed
pattern trajectory, so that $d_{n,i}=0$ in those cases. Then it is clear that
the average noise distance is smaller than the one predicted based on the
assumption of diverging trajectories:

$\displaystyle\left\langle d_{n,i}\right\rangle_{i}=$
$\displaystyle\left(1-p_{\text{no flip}}\right)d_{n}^{\text{residual}}.$

We can estimate the probability that no flip occurred due to the noise by
using the results from Appendix .9, where we calculate the average number of
flips in the network after one time constant given an additional (noise)
variance in the input of the neurons $\sigma_{\Delta h}^{2}$. In our present
case the noise is given by adding $\xi_{i}\sim\mathcal{N}(0,\sigma^{2})$ on
the output activities of the $L$ original neurons of the pattern, so that the
corresponding input variance felt by all neurons in the network is

$\displaystyle\tilde{\sigma}_{\Delta h}^{2}=$
$\displaystyle\frac{g^{2}}{N}\sigma^{2}L.$

Now we need only to consider that actually the variance is not constant for a
complete round of $N$ updates, but linearly diminishes every time one of the
$L$ source neurons is updated until none is left. Since the flip probability
depends on the square-root of the variance, there is a corrective factor
$c_{\text{dim}}=\sqrt{1-\frac{k}{N}}$ in each term of the product:

$\displaystyle p_{\text{no flip}}$
$\displaystyle=\prod_{k=0}^{N-1}\left(1-p_{\text{single flip}}(k)\right)$
$\displaystyle=\prod_{k=0}^{N-1}\left(1-\frac{1}{N}\bar{n}_{\mathrm{spawns}}\left(\tilde{\sigma}_{\Delta
h}\right)c_{\text{dim}}(k)\right)$
$\displaystyle=\prod_{k=0}^{N-1}\left(1-g\sigma\sqrt{\frac{2L}{\pi
N}}\left\langle\mathrm{T}_{p}^{\prime}(h)\right\rangle_{h\sim\mathcal{N}(R,Q_{0})}\sqrt{1-\frac{k}{N}}\right).$
(80)

This result fits well with the simulations, yielding the predicted offset of
the asymptotic average signal- and noise distances shown in Figure 7b and the
asymptotic plateau of the approximated average signal in Figure 7c.

###  .13 Description of simulations

Simulations for Figure 5 and Figure 4b,d were implemented using NEST (Linssen
_et al._ , 2018). NEST treats binary neurons in the bitlike $\\{0,1\\}$
representation. To let every neuron “see” inputs from $\\{-1,1\\}$ (Ising
spins) we add to each neuron $i$ a bias $\sum_{j}J_{ij}$ and then connect the
neurons by the connections $2J_{ij}$ instead of $J_{ij}$; thereby effectively
simulating an Ising system. To obtain the autocorrelations for the Ising case,
(41) is used, leading to the result shown in Figure 5. Furthermore, we use a
non-point-symmetric activation function $\mathrm{T}(h)=\text{tanh}(h-\Theta)$
by choosing a $\Theta$ to be nonzero for this plot. The reason is, first, that
$\Theta=0$ leads to the theoretical prediction of maximal output variance
$\llangle x^{2}\rrangle=1$ because the mean output activity $\left\langle
x\right\rangle$ is $0$. However, due to disorder, the time-averaged activity
is actually a fluctuating quantity across the population and therefore, the
population-averaged variance is always below $1$ for finite systems. This
systematic underestimation of the peak of the autocorrelation at zero time lag
can be avoided only by choosing an activation function that is not point
symmetric. Second, this choice is also natural because mean activity $0$ would
imply that neurons are active half of the time on average (Kühn and Helias,
2017, Supp. Mat. II B), which is considerably more than indicated by the low
firing rates measured in cortex (Roxin  _et al._ , 2011). For Figure 5, we
have therefore shift the working point by numerically inverting $\left\langle
x\right\rangle=\left\langle\text{tanh}(h-\Theta)\right\rangle_{h\sim\mathcal{N}(R,Q)}$
to obtain the value for $\Theta$, which in mean-field approximation
corresponds to $\left\langle x\right\rangle=-0.5$. In the simulations it turns
out that $\left\langle x\right\rangle=-0.501$, which agrees well with mean-
field theory.

For Figure 4b and d, for each point of the grid one simulation of two
identical networks was performed. After $1000\,\text{ms}$, in one replica, the
first two neurons are set to the active state and the third and fourth neurons
are set to the inactive state, after which the simulation continues for
$2500\,\text{ms}$. This method of perturbation entails the small probability
that these four neurons are already in exactly this state, so that nothing is
changed; this is the explanation for the scattered single green dots in Figure
4b,d. The advantage of the method is, however, that it guarantees the same
state of the random number generators across both replicas.

The simulations for Figure 7, Figure 9 and Figure 10 are performed using a
custom FORTRAN kernel.

The simulations of the LIF network for Figure 11 are implemented in NEST,
using the “iaf_psc_delta” neuron with its default parameter settings
$V_{\mathrm{rest}}=V_{\mathrm{reset}}=-70\,\mathrm{mV,}\,V_{th}=-55\,\mathrm{mV,}\,\tau_{m}=10\,\mathrm{ms},\,t_{\mathrm{refr}}=2\,\mathrm{ms},\,C_{m}=250\,\mathrm{pF}$.
Asynchronous-irregular firing in the inhibitory network is evoked by supplying
external excitatory Poisson input with rate of $3000\,\mathrm{Hz}$ and unit
weight, leading to an average network activity of $19\,\mathrm{Hz}$. After
simulating for $10\,\tau_{m}$ to obtain a state with stationary statistics,
the input pattern is applied by causing a spike in the corresponding neurons.
Noise is added to the pattern by additional external input spikes perturbing
the membrane potentials of the $L$ pattern neurons. Further parameters as
mentioned in figure caption:
$N=500,\,J_{ii}=-1\mathrm{mV},\,K=125,\,\tau_{m}=10\,\mathrm{ms}$.

The LSTM network for Figure 12 is simulated using the vanilla PYTORCH
implementation with a hidden layer of size $N=200$. To obtain chaotic
fluctuations over a range of timescales, the hidden weights
$W_{hi},\,W_{hf},\,W_{hg}$ are initialized as $\sim\mathcal{N}(0,5.8/N)$. This
excludes the $W_{ho}$ weights, which would cause very rapid, erratic dynamics.
Instead, these and all remaining weights and biases used the default uniformly
distributed $\sim\mathcal{U}(-1/\sqrt{N},1/\sqrt{N})$ initialization. After
simulating for $50$ time steps to obtain stationary statistics, the input
patterns are supplied using the standard input function.

Analysis of simulation data and numerical solutions are implemented in PYTHON.
The code to generate all figures is available as a Zenodo archive at
doi.org/10.5281/zenodo.4705262.

## References

  * Kempter  _et al._ (1998) R. Kempter, W. Gerstner, and J. L. Van Hemmen,  _Spike-based compared to rate-based Hebbian learning,_NIPS’98: Proceedings of the 11th International Conference on Neural Information Processing Systems , 125 (1998).
  * Pfeiffer and Pfeil (2018) M. Pfeiffer and T. Pfeil,  _Deep learning with spiking neurons: Opportunities and challenges,_Frontiers in Neuroscience 12, 774 (2018).
  * Laughlin (2001) S. Laughlin,  _Energy as a constraint on the coding and processing of sensory information,_Current Opinion in Neurobiology 11, 475 (2001).
  * Abeles (1991) M. Abeles,  _Corticonics: Neural Circuits of the Cerebral Cortex_ (Cambridge University Press, Cambridge, 1991).
  * Softky and Koch (1993) W. R. Softky and C. Koch,  _The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs,_J. Neurosci. 13, 334 (1993).
  * König  _et al._ (1996) P. König, A. K. Engel, and W. Singer,  _Integrator or coincidence detector? The role of the cortical neuron revisited,_ TINS 19, 130 (1996).
  * Shadlen and Newsome (1998) M. N. Shadlen and W. T. Newsome,  _The variable discharge of cortical neurons: Implications for connectivity, computation, and information coding,_J. Neurosci. 18, 3870 (1998).
  * Shadlen and Movshon (1999) M. N. Shadlen and A. J. Movshon,  _Synchrony unbound: A critical evaluation of the temporal binding hypothesis,_ Neuron 24, 67 (1999).
  * Rolls and Deco (2010) E. T. Rolls and G. Deco,  _The Noisy Brain: Stochastic Dynamics as a Principle_ (Oxford University Press, 2010).
  * Brette (2015) R. Brette,  _Philosophy of the spike: Rate-based vs. spike-based theories of the brain,_Frontiers in Systems Neuroscience 9, 151 (2015).
  * Denève and Machens (2016) S. Denève and C. K. Machens,  _Efficient codes and balanced networks,_Nature Neuroscience 19, 375 (2016).
  * Buonomano and Merzenich (1995) D. V. Buonomano and M. M. Merzenich,  _Temporal information transformed into a spatial code by a neural network with realistic properties,_Science 267, 1028 (1995), http://science.sciencemag.org/content/267/5200/1028.full.pdf .
  * Jaeger (2001) H. Jaeger,  _The “echo state” approach to analysing and training recurrent neural networks_ , Tech. Rep. GMD Report 148 (German National Research Center for Information Technology, St. Augustin, Germany, 2001).
  * Maass  _et al._ (2002) W. Maass, T. Natschläger, and H. Markram,  _Real-time computing without stable states: a new framework for neural computation based on perturbations,_ Neural Computation 14, 2531 (2002).
  * Vapnik (1998) V. N. Vapnik,  _Adaptive and learning systems for signal processing communications, and control,_ Statistical learning theory (1998).
  * Cover (1965) T. M. Cover,  _Geometrical and statistical properties of systems of linear inequalities with applications in pattern recognition,_IEEE Transactions on Electronic Computers EC-14, 326 (1965).
  * Toyoizumi and Abbott (2011) T. Toyoizumi and L. F. Abbott,  _Beyond the edge of chaos: Amplification and temporal integration by recurrent networks in the chaotic regime,_ Phys. Rev. E 84, 051908 (2011).
  * Bertschinger and Natschläger (2004) N. Bertschinger and T. Natschläger,  _Real-time computation at the edge of chaos in recurrent neural networks,_ Neural Computation 16, 1413 (2004).
  * Legenstein and Maass (2007a) R. Legenstein and W. Maass, What makes a dynamical system computationally powerful? in  _New Directions in Statistical Signal Processing: From System to Brains_ , edited by S. Haykin, J. C. Principe, T. J. Sejnowski, and J. G. McWhirter (MIT Press, 2007) pp. 127–154.
  * Legenstein and Maass (2007b) R. Legenstein and W. Maass,  _Edge of chaos and prediction of computational performance for neural circuit models,_Neural Networks 20, 323 (2007b).
  * Sompolinsky  _et al._ (1988) H. Sompolinsky, A. Crisanti, and H. J. Sommers,  _Chaos in random neural networks,_Phys. Rev. Lett. 61, 259 (1988).
  * Crisanti and Sompolinsky (2018) A. Crisanti and H. Sompolinsky,  _Path integral approach to random neural networks,_Phys. Rev. E 98, 062120 (2018).
  * Schuecker  _et al._ (2018) J. Schuecker, S. Goedeke, and M. Helias,  _Optimal sequence memory in driven random networks,_Phys Rev X 8, 041029 (2018).
  * van Vreeswijk and Sompolinsky (1996) C. van Vreeswijk and H. Sompolinsky,  _Chaos in neuronal networks with balanced excitatory and inhibitory activity,_Science 274, 1724 (1996).
  * van Vreeswijk and Sompolinsky (1998) C. van Vreeswijk and H. Sompolinsky,  _Chaotic balanced state in a model of cortical circuits,_Neural Comput. 10, 1321 (1998).
  * Baxendale (1992) P. H. Baxendale, in  _Diffusion Processes and Related Problems in Analysis, Volume II: Stochastic Flows_ (Birkhäuser, Boston, MA, Boston, MA, USA, 1992) pp. 3–35.
  * van Vreeswijk and Sompolinsky (1998) C. van Vreeswijk and H. Sompolinsky,  _Chaotic balanced state in a model of cortical circuits,_Neural Comput. 10, 1321 (1998).
  * Glauber (1963) R. Glauber,  _Time-dependent statistics of the Ising model,_ J. Math. Phys. 4, 294 (1963).
  * Ginzburg and Sompolinsky (1994) I. Ginzburg and H. Sompolinsky,  _Theory of correlations in stochastic neural networks,_Phys. Rev. E 50, 3171 (1994).
  * Cox (1962) D. R. Cox,  _Renewal Theory_ , Science Paperbacks (Chapman and Hall, London, 1962).
  * Papoulis and Pillai (2002) A. Papoulis and S. U. Pillai,  _Probability, Random Variables, and Stochastic Processes_ , 4th ed. (McGraw-Hill, Boston, 2002).
  * Sussillo and Abbott (2009) D. Sussillo and L. F. Abbott,  _Generating coherent patterns of activity from chaotic neural networks,_ Neuron 63, 544 (2009).
  * Kadmon and Sompolinsky (2015) J. Kadmon and H. Sompolinsky,  _Transition to chaos in random neuronal networks,_Phys. Rev. X 5, 041030 (2015).
  * Derrida and Pomeau (1986) B. Derrida and Y. Pomeau,  _Random networks of automata: a simple annealed approximation,_ EPL (Europhysics Letters) 1, 45 (1986).
  * Puelma Touzel and Wolf (2019) M. Puelma Touzel and F. Wolf,  _Statistical mechanics of spike events underlying phase space partitioning and sequence codes in large-scale models of neural circuits,_Physical Review E 99, 052402 (2019).
  * Hwang  _et al._ (2019) S. Hwang, V. Folli, E. Lanza, G. Parisi, G. Ruocco, and F. Zamponi,  _On the number of limit cycles in asymmetric neural networks,_Journal of Statistical Mechanics: Theory and Experiment 2019, 053402 (2019).
  * Turner  _et al._ (1976) M. E. Turner, E. L. Bradley, K. A. Kirk, and K. M. Pruitt,  _A theory of growth,_Mathematical Biosciences 29, 367 (1976).
  * Wainrib and Touboul (2013) G. Wainrib and J. Touboul,  _Topological and dynamical complexity of random neural networks,_ Phys. Rev. Lett. 110, 118101 (2013).
  * Engelken  _et al._ (2020) R. Engelken, F. Wolf, and L. F. Abbott, Lyapunov spectra of chaotic recurrent neural networks, (2020), arXiv:2006.02427 .
  * Hochreiter and Schmidhuber (1997) S. Hochreiter and J. Schmidhuber,  _Long short-term memory,_ Neural computation 9, 1735 (1997).
  * Monteforte and Wolf (2012) M. Monteforte and F. Wolf,  _Dynamic Flux Tubes Form Reservoirs of Stability in Neuronal Circuits,_Physical Review X 2, 041007 (2012).
  * Kauffman (1993) S. A. Kauffman,  _The Origins of Order Self-Organization and Selection in Evolution_ (Oxford University Press, 1993).
  * Snyder  _et al._ (2012) D. Snyder, A. Goudarzi, and C. Teuscher, in  _Artificial Life 13_ (MIT Press, 2012).
  * Grytskyy  _et al._ (2013) D. Grytskyy, T. Tetzlaff, M. Diesmann, and M. Helias,  _A unified view on weakly correlated recurrent networks,_Front. Comput. Neurosci. 7, 131 (2013).
  * Kühn (2020) T. Kühn,  _Path integral methods for correlated activity in neuronal networks,_RWTH Publications (2019 – 2020), 10.18154/RWTH-2020-01833.
  * Tetzlaff  _et al._ (2012) T. Tetzlaff, M. Helias, G. T. Einevoll, and M. Diesmann,  _Decorrelation of neural-network activity by inhibitory feedback,_PLOS Comput. Biol. 8, e1002596 (2012).
  * Aljadeff  _et al._ (2015) J. Aljadeff, M. Stern, and T. Sharpee,  _Transition to chaos in random networks with cell-type-specific connectivity,_Phys. Rev. Lett. 114, 088101 (2015).
  * Kuśmierz  _et al._ (2019) Ł. Kuśmierz, S. Ogawa, and T. Toyoizumi,  _Edge of chaos and scale-free avalanches in neural networks with heavy-tailed synaptic disorder,_ arXiv (2019), arXiv:1910.05780 .
  * Martin  _et al._ (1973) P. Martin, E. Siggia, and H. Rose,  _Statistical dynamics of classical systems,_ Phys. Rev. A 8, 423 (1973).
  * De Dominicis (1976) C. De Dominicis,  _Techniques de renormalisation de la théorie des champs et dynamique des phénomènes critiques,_J. Phys. Colloques 37, C1 (1976).
  * Chow and Buice (2015) C. Chow and M. Buice,  _Path integral methods for stochastic differential equations,_ J Math. Neurosci 5, 8 (2015).
  * Hertz  _et al._ (2017) J. A. Hertz, Y. Roudi, and P. Sollich,  _Path integral methods for the dynamics of stochastic and disordered systems,_Journal of Physics A: Mathematical and Theoretical 50, 033001 (2017).
  * Helias and Dahmen (2020) M. Helias and D. Dahmen,  _Statistical Field Theory for Neural Networks_, Vol. 970 (Springer International Publishing, 2020) p. 203.
  * Doi (1976) M. Doi,  _Second quantization representation for classical many-particle system,_Journal of Physics A: Mathematical and General 9, 1465 (1976).
  * Peliti (1985) L. Peliti,  _Path integral approach to birth-death processes on a lattice,_J. Phys. France 46, 1469 (1985).
  * Buice and Cowan (2007) M. A. Buice and J. D. Cowan,  _Field-theoretic approach to fluctuation effects in neural networks,_Phys. Rev. E 75, 051919 (2007).
  * Sommers (1987) H. Sommers,  _Path-integral approach to ising spin-glass dynamics,_ Phys. Rev. Lett. 58, 1268 (1987).
  * Andreanov  _et al._ (2006) A. Andreanov, G. Biroli, J.-P. Bouchaud, and A. Lefevre,  _Field theories and exact stochastic equations for interacting particle systems,_Phys. Rev. E 74, 030101(R) (2006).
  * Lefevre and Biroli (2007) A. Lefevre and G. Biroli,  _Dynamics of interacting particle systems: stochastic process and field theory,_Journal of Statistical Mechanics: Theory and Experiment 2007, P07024 (2007).
  * Büsing  _et al._ (2010) L. Büsing, B. Schrauwen, and R. Legenstein,  _Connectivity, dynamics, and memory in reservoir computing with binary and analog neurons,_ Neural Comput. 22, 1272 (2010).
  * Tetzlaff  _et al._ (2004) T. Tetzlaff, A. Morrison, T. Geisel, and M. Diesmann,  _Consequences of realistic network size on the stability of embedded synfire chains,_ Neurocomputing 58–60, 117 (2004).
  * Shea-Brown  _et al._ (2008) E. Shea-Brown, K. Josic, J. de la Rocha, and B. Doiron,  _Correlation and synchrony transfer in integrate-and-fire neurons: basic properties and consequences for coding,_ Phys. Rev. Lett. 100, 108102 (2008).
  * Tchumatchenko  _et al._ (2010) T. Tchumatchenko, A. Malyshev, T. Geisel, M. Volgushev, and F. Wolf,  _Correlations and synchrony in threshold neuron models,_ Phys. Rev. Lett. 104, 058102 (2010).
  * Tchumatchenko  _et al._ (2011) T. Tchumatchenko, T. Geisel, M. Volgushev, and F. Wolf,  _Spike correlations – what can they tell about synchrony?_Frontiers in Neuroscience 5, 68 (2011).
  * Tetzlaff  _et al._ (2003) T. Tetzlaff, M. Buschermöhle, T. Geisel, and M. Diesmann,  _The spread of rate and correlation in stationary cortical networks,_ Neurocomputing 52–54, 949 (2003).
  * Shea-Brown  _et al._ (2007) E. Shea-Brown, K. Josic, J. De la Rocha, and B. Doiron,  _Universal properties of correlation transfer in integrate-and-fire neurons,_ Arxiv:q-bio (2007).
  * Schultze-Kraft  _et al._ (2013) M. Schultze-Kraft, M. Diesmann, S. Gruen, and M. Helias,  _Noise suppression and surplus synchrony by coincidence detection,_PLOS Comput. Biol. 9, e1002904 (2013).
  * Deniz and Rotter (2017) T. Deniz and S. Rotter,  _Solving the two-dimensional fokker-planck equation for strongly correlated neurons,_Physical Review E 95 (2017), 10.1103/physreve.95.012412.
  * Lajoie  _et al._ (2013) G. Lajoie, K. K. Lin, and E. Shea-Brown,  _Chaos and reliability in balanced spiking networks with temporal drive,_ Phys. Rev. E 87, 052901 (2013).
  * Lajoie  _et al._ (2014) G. Lajoie, J.-P. Thivierge, and E. Shea-Brown,  _Structured chaos shapes spike-response noise entropy in balanced neural networks,_Frontiers in Computational Neuroscience 8, 123 (2014).
  * Gardner and Derrida (1988) E. Gardner and B. Derrida,  _Optimal storage properties of neural network models,_Journal of Physics A: Mathematical and General 21, 271 (1988).
  * Litwin-Kumar  _et al._ (2012) A. Litwin-Kumar, M. J. Chacron, and B. Doiron,  _The spatial structure of stimuli shapes the timescale of correlations in population spiking activity,_ PLOS Comput. Biol. 8, e1002667 (2012).
  * Poole  _et al._ (2016) B. Poole, S. Lahiri, M. Raghu, J. Sohl-Dickstein, and S. Ganguli, in  _Advances in Neural Information Processing Systems 29_, edited by D. D. Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett (Curran Associates, Inc., 57 Morehouse Lane; Red Hook, NY 12571, 2016) pp. 3360–3368.
  * Recanatesi  _et al._ (2019) S. Recanatesi, M. Farrell, M. Advani, T. Moore, G. Lajoie, and E. Shea-Brown,  _Dimensionality compression and expansion in Deep Neural Networks,_ArXiv e-prints (2019), 1906.00443 .
  * Farrell  _et al._ (2019) M. Farrell, S. Recanatesi, T. Moore, G. Lajoie, and E. Shea-Brown,  _Recurrent neural networks learn robust representations by dynamically balancing compression and expansion,_bioRxiv , 564476 (2019), 564476 .
  * Hennequin  _et al._ (2012) G. Hennequin, T. Vogels, and W. Gerstner,  _Non-normal amplification in random balanced neuronal networks,_ Phys. Rev. E 86, 011909 (2012).
  * Kerg  _et al._ (2019) G. Kerg, K. Goyette, M. Puelma Touzel, G. Gidel, E. Vorontsov, Y. Bengio, and G. Lajoie, in  _Advances in Neural Information Processing Systems 32_, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 57 Morehouse Lane; Red Hook, NY 12571, 2019) pp. 13613–13623.
  * Tarnowski (2020) W. Tarnowski,  _Transient amplification in balanced neural networks,_ArXiv e-prints (2020), 2011.08215 .
  * Bahri  _et al._ (2020) Y. Bahri, J. Kadmon, J. Pennington, S. S. Schoenholz, J. Sohl-Dickstein, and S. Ganguli,  _Statistical mechanics of deep learning,_Annual Review of Condensed Matter Physics 11, 501 (2020).
  * Thorpe  _et al._ (1996) S. Thorpe, D. Fize, and C. Marlot,  _Speed of processing in the human visual system,_ Nature 381, 520 (1996).
  * Hung  _et al._ (2005) C. P. Hung, G. Kreiman, T. Poggio, and J. J. DiCarlo,  _Fast readout of object identity from macaque inferior temporal cortex,_Science 310, 863 (2005).
  * Lajoie  _et al._ (2016) G. Lajoie, K. K. Lin, J.-P. Thivierge, and E. Shea-Brown,  _Encoding in balanced networks: Revisiting spike patterns and chaos in stimulus-driven systems,_PLOS Computational Biology 12, e1005258 (2016).
  * Friedrich and Laurent (2001) R. Friedrich and G. Laurent,  _Dynamic optimization of odor representations by slow temporal patterning of mitral cell activity,_ Science 291, 889 (2001).
  * Mazor and Laurent (2005) O. Mazor and G. Laurent,  _Transient dynamics versus fixed points in odor representations by locust antennal lobe projection neurons,_Neuron 48, 661 (2005).
  * Cury and Uchida (2010) K. M. Cury and N. Uchida,  _Robust odor coding via inhalation-coupled transient activity in the mammalian olfactory bulb,_Neuron 68, 570 (2010).
  * Grassberger and Procaccia (1983) P. Grassberger and I. Procaccia,  _Estimation of the kolmogorov entropy from a chaotic signal,_Phys. Rev. A 28, 2591 (1983).
  * Toker  _et al._ (2020) D. Toker, F. T. Sommer, and M. D’Esposito,  _A simple method for detecting chaos in nature,_Communications Biology 3, 11 (2020).
  * Celletti and Villa (1996) A. Celletti and A. E. P. Villa,  _Low-dimensional chaotic attractors in the rat brain,_Biological Cybernetics 74, 387 (1996).
  * Churchland  _et al._ (2010) M. M. Churchland, B. M. Yu, J. P. Cunningham, L. P. Sugrue, M. R. Cohen, G. S. Corrado, W. T. Newsome, A. M. Clark, P. Hosseini, B. B. Scott, D. C. Bradley, M. A. Smith, A. Kohn, J. A. Movshon, K. M. Armstrong, T. Moore, S. W. Chang, L. H. Snyder, S. G. Lisberger, N. J. Priebe, I. M. Finn, D. Ferster, S. I. Ryu, G. Santhanam, M. Sahani, and K. V. Shenoy,  _Stimulus onset quenches neural variability: a widespread cortical phenomenon,_Nat. Neurosci. 13, 369 (2010).
  * Mazzucato  _et al._ (2016) L. Mazzucato, A. Fontanini, and G. La Camera,  _Stimuli reduce the dimensionality of cortical activity,_Frontiers in systems neuroscience 10, 11 (2016).
  * Gao  _et al._ (2017) P. Gao, E. Trautmann, B. M. Yu, G. Santhanam, S. Ryu, K. Shenoy, and S. Ganguli,  _A theory of multineuronal dimensionality, dynamics and measurement,_bioRxiv , 214262 (2017), arXiv:bioRxiv 214262 .
  * Coolen (2000) A. C. C. Coolen,  _Statistical mechanics of recurrent neural networks ii. dynamics,_ arXiv:cond-mat/0006011 (2000).
  * Helias  _et al._ (2014) M. Helias, T. Tetzlaff, and M. Diesmann,  _The correlation structure of local cortical networks intrinsically results from recurrent dynamics,_PLOS Comput. Biol. 10, e1003428 (2014).
  * Uhlenbeck and Ornstein (1930) G. E. Uhlenbeck and L. S. Ornstein,  _On the theory of the brownian motion,_ Phys. Rev. 36, 823 (1930).
  * Note (1) (note that our notational convention differs by the factor $1/2$).
  * Lindner (2009) B. Lindner, A brief introduction to some simple stochastic processes, in  _Stochastic Methods in Neuroscience_ , edited by C. Laing and G. J. Lord (OUP Oxford, 2009).
  * Frey and Weber (2017) E. Frey and M. F. Weber,  _Master equations and the theory of stochastic path integrals,_Rep. Prog. Phys. 80, 046601 (2017).
  * Rajan  _et al._ (2010) K. Rajan, L. Abbott, and H. Sompolinsky,  _Stimulus-dependent suppression of chaos in recurrent neural networks,_ Phys. Rev. E 82, 011903 (2010).
  * Martí  _et al._ (2018) D. Martí, N. Brunel, and S. Ostojic,  _Correlations between synapses in pairs of neurons slow down dynamics in randomly connected neural networks,_Phys. Rev. E 97, 062314 (2018).
  * Linssen  _et al._ (2018) C. Linssen, M. E. Lepperød, J. Mitchell, J. Pronold, J. M. Eppler, C. Keup, A. Peyser, S. Kunkel, P. Weidel, Y. Nodem, D. Terhorst, R. Deepu, M. Deger, J. Hahne, A. Sinha, A. Antonietti, M. Schmidt, L. Paz, J. Garrido, T. Ippen, L. Riquelme, A. Serenko, T. Kühn, I. Kitayama, H. Mørk, S. Spreizer, J. Jordan, J. Krishnan, M. Senden, E. Hagen, A. Shusharin, S. B. Vennemo, D. Rodarie, A. Morrison, S. Graber, J. Schuecker, S. Diaz, B. Zajzon, and H. E. Plesser, Nest 2.16.0,  (2018).
  * Kühn and Helias (2017) T. Kühn and M. Helias,  _Locking of correlated neural activity to ongoing oscillations,_PLOS Comput. Biol. 13, e1005534 (2017).
  * Roxin  _et al._ (2011) A. Roxin, N. Brunel, D. Hansel, G. Mongillo, and C. van Vreeswijk,  _On the distribution of firing rates in networks of cortical neurons,_J. Neurosci. 31, 16217 (2011).


# Effective quasiclassical models for odd-frequency superconductivity:
energy-symmetry, preserved spectral weight, and Meissner response

Hans Gløckner Giil  Center for Quantum Spintronics, Department of Physics,
Norwegian
University of Science and Technology, NO-7491 Trondheim, Norway  Jacob Linder
<EMAIL_ADDRESS>Center for Quantum Spintronics, Department of Physics,
Norwegian
University of Science and Technology, NO-7491 Trondheim, Norway

###### Abstract

The odd-frequency superconducting state appears generally in hybrid structures
consisting of conventional superconductors and other materials, and features
electrons that form temporally non-local Cooper pairs. The quasiclassical
theory of superconductivity has been extensively used to model such systems,
finding in many cases excellent agreement with experimental measurements.
Therefore, it is of interest to study effective models of odd-frequency
superconductivity to predict new phenomena associated with this form of
unconventional pairing. We establish necessary criteria that the
quasiclassical Green functions in odd-frequency superconducting systems in the
dirty limit must satisfy in order to be physically reasonable, including
conservation of spectral weight. We show that it is possible to write down
effective models which satisfy all the abovementioned criteria, but which
still predict different behavior when it comes to the density of states and
the magnetic response of the superconductor. For instance, an effective model
for the odd-frequency anomalous Green function that gives a conserved spectral
weight can yield either a peaked or gapped density of states at the Fermi
energy, and exhibit conventional, zero, or unconventional Meissner response.
This finding demonstrates the importance of carefully considering the
properties of effective models describing odd-frequency superconductivity in
order to obtain physically reasonable results.

##  I Introduction

In the theory of superconductivity, symmetry is a key component which dictates
several properties of the elementary building block of superconductors: the
Cooper pair. The function mathematically describing how the two electrons
making up the Cooper pair correlate to each other depends, even in the
simplest case of a single-band superconductor, on the position, spin, and time
coordinate of the electrons. In BCS theory [1], the time coordinate is usually
ignored after introducing a cutoff for the energy range that the electrons
attract each other within. Nevertheless, the symmetry property of a paired
Cooper pair state allows for the intriguing possibility that the two electrons
are not correlated at equal times and that their correlation is instead
established as the time separation grows. This is the case if the correlation
function is odd in the relative time coordinate of the electrons.

Berezinskii [2] predicted that a two-electron pairing correlation, with
temporal coordinates $t_{1}$ and $t_{2}$, could be odd in the relative time
coordinate $t_{1}-t_{2}$. An equivalent way of expressing this is to say that
the time-ordered correlation function is odd in frequency, where the frequency
is the Fourier-transform of the relative time coordinate. Initially considered
as an intrinsic pairing instability [5, 3, 4, 6], the field of odd-frequency
superconductivity experienced a breakthrough when it was shown [7] that such
unconventional pairing could be generated by placing a conventional BCS
superconductor in contact with a ferromagnet [8, 9]. Although no intrinsic
bulk odd-frequency superconductors are known to date, it is by now well-
established both theoretically and experimentally that such pairing exists in
hybrid structures of essentially any sort featuring at least one conventional
superconductor [38, 39, 10].

Odd-frequency superconductivity has been shown to give rise to superconducting
properties that differ fundamentally from those of conventional $s$-wave
superconductors such as Al and Nb. This includes phenomena such as spin
supercurrents [11, 12, 13, 14], paramagnetic Meissner effects [15, 16],
gapless superconductivity [17, 18], and the appearance of Majorana bound
states [19, 20]. Because of this, it is of interest to consider effective
models describing odd-frequency superconductivity in order to predict new
emergent phenomena due to this interesting type of long-range order.

Previous works [23, 24, 25, 26] have considered effective models for the
anomalous Green function of an odd-frequency superconductor. Most of these
utilized a methodology based on the quasiclassical theory of superconductivity
[27, 28, 29], which is known to compare well with experimental measurements on
mesoscopic metallic and superconducting systems. These models have satisfied
the appropriate symmetries with respect to inversion of energy [30],
fulfilling also the normalization condition $\check{g}^{2}=1$ of the
quasiclassical Green function $\check{g}$ (some authors use a different
convention with $\check{g}^{2}=-\pi^{2}$). However, a Green function
satisfying all of these criteria can still display physically problematic
behavior. This happens in at least two possible ways. One problem which
previous works have suffered from is that the spectral weight of the system is
not conserved [24] when transitioning from the normal state to the
superconducting state, i.e. a sum rule for the spectral weight is violated, as
pointed out in Ref. [25]. Typically, spectral weight is experimentally found
to be redistributed energetically when entering the superconducting state, but
the number of states should remain the same. Moreover, as we will demonstrate
below, an anomalous Green function which satisfies both the correct (i)
symmetry relation under energy-inversion for odd-frequency superconductivity,
(ii) a conserved spectral weight in the superconducting state, and (iii) a
normalized quasiclassical Green function, can still exhibit an arbitrary
Meissner response: diamagnetic, paramagnetic, or even absent. In addition, it
can give rise to either a gapped or peaked density of states at the Fermi
level.

Because of these issues, we believe it is useful to clarify the precise form
of the Green function which describes an odd-frequency superconductor which
conserves the spectral weight of the normal-state, and to also predict which
Meissner response and density of states it produces. Below, we outline four
necessary, but not sufficient, criteria which a physically correct Green
function describing an odd-frequency superconductor should satisfy. Then, we
consider different models for such Green functions and show how they produce
qualitatively different Meissner responses.

##  II Theory

###  II.1 Preliminary: properties of odd-frequency Green function in S/F
bilayers

It is instructive to start by studying a superconductor/ferromagnet (S/F)
system, in which odd-frequency arises naturally through the proximity effect.
We assume that the superconductor retains its BCS-solution, and consider the
effect on the ferromagnet. The Usadel equation [28] holds in the diffusive
limit and can, in the absence of electromagnetic fields and spin-flip and
spin-orbit scattering, be written as

$D\partial_{x}\bigl{(}\hat{g}\partial_{x}\hat{g}\bigr{)}=-\mathrm{i}[E\hat{\tau}_{3}+\hat{M},\hat{g}],$
(1)

where
$\hat{M}=\bm{h}\cdot\operatorname{diag}(\underline{\bm{\sigma}},\underline{\bm{\sigma}}^{*})$
is the ferromagnetic term with $\bm{h}$ being the spin-splitting energy,
$\hat{g}$ is the retarded impurity averaged quasiclassical Green function,
$\hat{\tau}_{3}=\operatorname{diag}(1,1,-1,-1)$, $D$ is the diffusion constant
of the system, $E$ is the energy relative to the Fermi energy, and
$\underline{\bm{\sigma}}$ is a 3-vector containing the Pauli matrices. In this
text, $\underline{A}$ denotes a $2\times 2$ matrix in spin space, while
$\hat{A}$ is a $4\times 4$ matrix in spin-particle-hole space, also called
spin-Nambu space [31]. Finally, $\check{A}$ denotes an $8\times 8$ matrix in
Keldysh-spin-Nambu space. The quasiclassical Green function has the form

$\displaystyle\check{g}=\begin{pmatrix}\hat{g}^{R}&\hat{g}^{K}\\\
\hat{0}&\hat{g}^{A}\\\ \end{pmatrix},$ (2)

where $\hat{g}^{R}$, $\hat{g}^{A}$, and $\hat{g}^{K}$ are the retarded,
advanced, and Keldysh components. In the following we will consider the
retarded Green function, dropping the suprescript. The general structure of
the retarded Green function contains the normal Green function $\underline{g}$
and the anomalous Green function $\underline{f}$,

$\hat{g}=\begin{pmatrix}\underline{g}&\underline{f}\\\
-\underline{\tilde{f}}&-\underline{\tilde{g}}\end{pmatrix}.$ (3)

Above, we have used the “tilde conjugation” operation, which for a function
$A(E)$ is defined as

$\tilde{A}(E)=[A(-E)]^{*},$ (4)

and the $2\times 2$ matrix $\underline{g}$ has components

$\underline{g}=\begin{pmatrix}g_{\uparrow\uparrow}&g_{\uparrow\downarrow}\\\
g_{\downarrow\uparrow}&g_{\downarrow\downarrow}\end{pmatrix},$ (5)

and similarly for the other matrices in spin space.

The superconducting order parameter $\Delta$ vanishes in the ferromagnet under
the assumption of zero intrinsic attractive interaction without contact with
the superconducting layer, and the anomalous Green function arises solely due
to the proximity effect that exists due to tunneling between the materials. In
the weak proximity effect regime, where $\hat{g}$ takes a form close to the
normal metal solution, we linearize the Usadel equation by writing

$\hat{g}=\hat{\tau}_{3}+\delta\hat{f},$ (6)

and neglect terms that are second order in $\delta\hat{f}$.

From the general symmetries of the quasiclassical Green functions, it can be
shown that the retarded and advanced components are related as [10]

$\displaystyle f_{\sigma\sigma^{\prime}}^{A}(E)$
$\displaystyle=-f_{\sigma\sigma^{\prime}}^{R}(-E),\text{ for triplet odd-
frequency}$ (7) $\displaystyle f_{\sigma\sigma^{\prime}}^{A}(E)$
$\displaystyle=+f_{\sigma\sigma^{\prime}}^{R}(-E),\text{ for singlet even-
frequency},$ (8)

where the superscripts indicate the retarded and advanced components of the
Green function. In the considered regime of diffusive transport, the even-
frequency correlations are spin-singlet whereas the odd-frequency correlations
are spin-triplet. Although these relations can be useful, they do nothing
towards our goal of finding criteria for the retarded Green function itself.
We will in the following derive symmetries for the retarded Green function
directly.

Introducing a singlet-triplet decomposition according to $\delta
f=(f_{s}+\bm{d}\cdot\bm{\sigma})\mathrm{i}\sigma_{2}$, where $f_{s}$ is the
singlet component and $\bm{d}$ are the three triplet components, the
linearized Usadel equation can be written as a set of four coupled
differential equations,

$\displaystyle\frac{\mathrm{i}D}{2}\partial_{x}^{2}f_{s}$
$\displaystyle=Ef_{s}+\bm{h}\cdot\bm{d}$ (9)
$\displaystyle\frac{\mathrm{i}D}{2}\partial_{x}^{2}\bm{d}$
$\displaystyle=E\bm{d}+\bm{h}f_{s}.$ (10)

In the superconductor, the solution for the anomalous Green function is given
by

$f_{BCS}=\frac{\Delta\Theta(E^{2}-\Delta^{2})\text{sgn}\\!\left(E\right)}{\sqrt{E^{2}-\Delta^{2}}}-\mathrm{i}\frac{\Delta\Theta(\Delta^{2}-E^{2})}{\sqrt{\Delta^{2}-E^{2}}},$
(11)

where $\Theta$ is the Heaviside step function and sgn is the sign function. It
is observed from these linearized equations that if $f_{s}=-\tilde{f}_{s}$,
which holds for a BCS superconductor in a real gauge as seen from Eq. (11),
then $\bm{d}=\tilde{\bm{d}}$ must hold. This relation continues to hold even
in the presence of Rashba and Dresselhaus spin-orbit coupling, in which case
[34] an additional term $2iDa^{2}\Omega(\chi){\bm{d}}$ appears on the lower
line of Eq. (9). In this term, $a$ denotes the strength of the SOC, and
$\Omega$ is a matrix with real entries where $\chi$ depends on the relative
magnitude of Rashba and Dresselhaus coupling. We proceed by using the
Kuprianov-Lukichev boundary conditions [32], which for the upper right
component of the anomalous Green function reads [8]

$\partial_{x}\underline{f}_{2}=-2\Omega_{2}\underline{f}_{1}$ (12)

when assuming the magnitudes of the components of $\underline{f}_{2}$ to be
much smaller than the magnitude of the components of $\underline{f}_{1}$. This
approximation holds in the weak proximity regime when the non-superconducting
material is in contact with a BCS superconductor at $T\ll T_{c}$. Here, 1
denotes the material to the left of the interface, and 2 the material to the
right of the interface, and we introduced the boundary and material specific
constant $\Omega_{2}=(2\zeta_{2}L_{2})^{-1}$, where $L_{2}$ is the length of
material $2$, and $\zeta_{2}$ is the ratio between the bulk resistance and the
interface resistance. Other types of boundary conditions can be used for
highly transparent interfaces [40], magnetic interfaces [41, 42], or spin-
orbit coupled interfaces [43, 44].

For concreteness, consider a magnetization in the z-direction inside the
ferromagnet. The solution for the anomalous Green function in the ferromagnet
then becomes

$\displaystyle\begin{split}d_{z}(x)&=C_{1}\mathrm{e}^{k_{-}x}+C_{2}\mathrm{e}^{-k_{-}x}+C_{3}\mathrm{e}^{k_{+}x}+C_{4}\mathrm{e}^{-k_{+}x}\\\
&\approx C_{2}\mathrm{e}^{-k_{-}x}+C_{4}\mathrm{e}^{-k_{+}x}\end{split}$ (13)
$\displaystyle\begin{split}f_{s}(x)&=-C_{1}\mathrm{e}^{k_{-}x}+-C_{2}\mathrm{e}^{-k_{-}x}+C_{3}\mathrm{e}^{k_{+}x}+C_{4}\mathrm{e}^{-k_{+}x}\\\
&\approx-C_{2}\mathrm{e}^{-k_{-}x}+C_{4}\mathrm{e}^{-k_{+}x},\end{split}$ (14)

where we defined

$\displaystyle k=\sqrt{-\frac{2\mathrm{i}E}{D}}$ $\displaystyle
k_{\pm}=\sqrt{-\frac{2\mathrm{i}(E\pm h)}{D}}.$ (15)

Above, $h$ is the magnitude of the spin-splitting of the conduction bands in
the ferromagnet, $\\{B_{j},C_{1},C_{2},C_{3},C_{4}\\}$ are constants, and we
have assumed semi-infinite ferromagnets, such that the exponentially
increasing parts of the solutions can safely be ignored. The solutions for
magnetizations in the x\- or y-direction are similar upon exchanging $d_{z}$
with $d_{x,y}$ functions. Applying Eq. (12), we find

$\displaystyle C_{2}=-\frac{\Omega_{2}f_{BCS}}{k_{-}}$ $\displaystyle
C_{4}=\frac{\Omega_{2}f_{BCS}}{k_{+}},$ (16)

from which one verifies that $\bm{d}$ is indeed symmetric under the tilde
conjugation operation, while $f_{s}$ obtains a sign change under the same
operation. These symmetries follow from the fact that $\tilde{k}=k$ and
$\tilde{k}_{\pm}=k_{\mp}$.

Importantly, we have also verified the property $\tilde{d}(E)=d(E)$ for the
odd-frequency anomalous Green function using a full numerical solution without
approximations, using the Ricatti parametrization, in S/F multilayers using a
real superconducting order parameter $\Delta=|\Delta|$ in S. In the case of a
complex $\Delta=|\Delta|\mathrm{e}^{\mathrm{i}\phi}$, the induced odd-
frequency correlator instead satisfies

$\tilde{d}(E)=d(E)\mathrm{e}^{-2\mathrm{i}\phi},$ (17)

which may alternatively be written $\tilde{d}(E,\phi)=d(E,-\phi).$

###  II.2 Required criteria for the Green function matrix

Motivated by the solution in the last section, we will now establish necessary
criteria that the retarded quasiclassical Green function should satisfy in the
odd-frequency case, including conserved spectral weight. Our goal will be to
find a model for a bulk odd-frequency superconductor that satisfy these
properties, as well as reproducing results that are expected for such
superconductors. We assume that we are in the diffusive regime of transport,
and only consider $s$-wave superconductivity, in which case triplets will be
odd in frequency and singlets will be even in frequency.

When restricting ourselves to a real order parameter, we showed in the
previous subsection that for an S/F system with a real order parameter in S,
the singlet and triplet anomalous Green functions satisfy

$\displaystyle\tilde{f}_{s}$ $\displaystyle=-f_{s}$ (18)
$\displaystyle\tilde{\bm{d}}$ $\displaystyle=\bm{d}.$ (19)

Now, the spectral weight function $\mathcal{A}(\bm{p},E)$ is defined from the
Green function $G$ as

$\mathcal{A}_{\sigma}(\bm{p},E)=-\frac{1}{\pi}{\Im\mathrm{m}}\\{G^{R}_{\sigma\sigma}\\},$
(20)

which when we apply the quasiclassical approximation,
$G=-\mathrm{i}\pi\delta(\xi_{p})g(\bm{p}_{F},E)$, becomes

$\mathcal{A}_{\sigma}(\bm{p},E)={\mathrm{Re}}\\{g_{\sigma\sigma}\delta(\xi_{p})\\},$
(21)

reflecting the fact that in the quasicassical approximation, the Green
function describes a low-energy theory for electronic states with momentum
equal to the Fermi momentum $\bm{p}_{F}$. This means that the standard sum
rule of the spectral weight, which reads

$\int_{-\infty}^{\infty}dE\mathcal{A}_{\sigma}(\bm{p},E)=1,$ (22)

must be modified accordingly. This sum rule states that the total number of
states is conserved in the system, which holds upon transitioning into the
superconducting state even though the Green function and spectral weight
changes.

Comparing a general solution to the bulk normal metal solution,
$\underline{g}=\underline{1}$, we may express this criterion using the
quasiclassical Green function as

$\int^{\omega_{c}}_{-\omega_{c}}dE(N_{\sigma}-N_{0})=0.$ (23)

Formally, this expression can be obtained by integrating Eq. (22) over momenta
to remove the delta functions, and introducing the density of states
$N_{\sigma}=\text{Re}\\{g_{\sigma\sigma}\\}$ as well as the density of states
at the Fermi surface, $N_{0}$. Above, $\omega_{c}$ is the high-energy cutoff
in quasiclassical theory, which in practice is chosen large enough to ensure
that the normal-state solution for the Green function is recovered at energies
close to $\omega_{c}$.

When deriving the Usadel equation, information about the normalization of the
Green function is lost when subtracting the left-handed and right-handed
equations of motion. We will here use the standard choice of normalization
$\check{g}^{2}=\check{1}$, which leads to for instance

$\underline{g}^{2}-\underline{f}\underline{\tilde{f}}=\underline{1}.$ (24)

This will serve as the definition for the normal component of the
quasiclassical Green function when we make ansätze for the anomalous
component. Finally, we expect the superconducting correlations to vanish for
large energies, analogous to the BCS case. This corresponds to a vanishing
anomalous Green function for large energies.

Summarizing, the four criteria we demand for the retarded Green function
describing odd-frequency superconductivity are as follows:

  * •

Energy-symmetry: ${\bm{d}}(E,\phi)=\tilde{\bm{d}}(E,-\phi)$

  * •

Conserved spectral weight:
$\int^{\omega_{c}}_{-\omega_{c}}dE(N_{\sigma}-N_{0})=0$.

  * •

Normalization: $\hat{g}^{2}=\hat{1}$.

  * •

Vanishing correlations for large energies: $\lim_{|E|\to\infty}{\bm{d}}(E)=0$.

###  II.3 Meissner response

Before proceeding to discuss the particular Meissner response exhibited by
various anomalous Green functions that satisfy the four criteria listed above,
it is useful to briefly sketch how odd-frequency pairing can fundamentally
change the (orbital) magnetic response of a superconductor.

In the London gauge, the Maxwell equation for a time-independent system reads

$\bm{\nabla}^{2}\bm{A}=-\bm{j},$ (25)

where $\bm{A}$ is the magnetic vector potential and $\bm{j}$ is the electric
current density. Inserting the London equation, $\bm{j}=-\bm{A}/\lambda_{L}$,
where $\lambda_{L}$ is the London penetration depth of the system, one finds
from the solution of the Maxwell equation an exponentially decaying magnetic
field, consistent with the Meissner effect in conventional superconductors.

The current density response to a magnetic field in a bulk superconductor can
be written as

$\bm{j}=-\frac{\mathrm{i}\bm{A}N_{0}e^{2}D}{16}\int_{-\infty}^{\infty}dE\,\text{Tr}\\!\left\\{\hat{\tau}_{3}\left(\check{g}\bigr{[}\hat{\tau}_{3},\check{g}\bigl{]}\right)^{K}\right\\},$
(26)

where the “K” superscript denotes the upper right component in Keldysh space.

This can be simplified considerably by assuming an equilibrium situation,

$\bm{j}=\bm{A}\frac{N_{0}e^{2}D}{2}I(x),$ (27)

where we introduced the temperature-dependent integral

$I(x)=\int_{-\infty}^{\infty}dE\,\tanh\\!\left(\frac{E}{2T}\right)\text{Tr}\\!\left\\{{\Im\mathrm{m}}\bigl{\\{}\underline{f}\underline{\tilde{f}}\bigr{\\}}\right\\},$
(28)

which will determine the sign of the current.

In the case of an unconventional Meissner response where
$\bm{j}=\bm{A}/\lambda_{L}$, the solution for the magnetic field becomes
oscillating [33, 15]. We now consider a one-dimensional superconducting system
in the x-direction, subject to an external magnetic field, such that the
values at the boundaries become $\bm{B}=B_{0}\bm{e}_{z}$. Gauge freedom then
allows us to choose $\bm{A}=A_{y}(x)\bm{e}_{y}$. The system length is $L$,
producing the boundary conditions

$\displaystyle\frac{dA_{y}}{dx}(x=0)=B_{0}$
$\displaystyle\frac{dA_{y}}{dx}(x=L)=B_{0}.$ (29)

It can then be shown that the conventional Meissner response produces a field

$B_{z}(x)=B_{0}\frac{\mathrm{e}^{\lambda_{L}^{-1}L}-1}{\mathrm{e}^{\lambda_{L}^{-1}L}-\mathrm{e}^{-\lambda_{L}^{-1}L}}\mathrm{e}^{-\lambda_{L}^{-1}x}+B_{0}\frac{\mathrm{e}^{-\lambda_{L}^{-1}L}-1}{\mathrm{e}^{\lambda_{L}^{-1}L}-\mathrm{e}^{-\lambda_{L}^{-1}L}}\mathrm{e}^{\lambda_{L}^{-1}x},$
(30)

while the unconventional response produces

$B_{z}(x)=B_{0}\cos(\lambda_{L}^{-1}x)-B_{0}\frac{\cos(\lambda_{L}^{-1}L)-1}{\sin(\lambda_{L}^{-1}L)}\sin(\lambda_{L}^{-1}x).$
(31)

Defining the magnetic susceptibility

$\chi(x)=\frac{B_{z}(x)}{B_{0}}-1,$ (32)

we plot in Fig. 1 the two responses. We clearly observe that the
unconventional Meissner response produces an oscillating susceptibility,
corresponding to both a local paramagnetic and diamagnetic response.

Figure 1: The magnetic susceptibility $\chi$ of a material with length
$L=14\lambda_{L}$ with a conventional (solid) and unconventional (dashed)
Meissner response. The values of $\chi<-1$ are unphysical.

The energy density stored in the magnetic field is $\bm{B}^{2}/(1+\chi)$,
meaning that for values $\chi<-1$, we find that the energy density becomes
negative, indicating an instability of the system [33]. Another way to see
that the paramagnetic Meissner effect is problematic for a bulk superconductor
is by examining the Landau expansion for the Gibbs free energy:

$\displaystyle F$
$\displaystyle=\int\Big{[}\alpha(T)|\psi|^{2}+\frac{1}{2}\beta(T)|\psi|^{4}+\gamma(T)|(\nabla+2\mathrm{i}e/\hbar
c)\psi|^{2}$ $\displaystyle+B^{2}/8\pi\Big{]}d^{3}x.$ (33)

This free energy must be minimized with respect to variations of $\psi$ and
$\bm{A}$. Minimization with respect to $\bm{A}$, $\delta F/\delta\bm{A}=0$,
gives

$\displaystyle\nabla\times\bm{B}=(4\pi/c)\bm{J},$ (34)

where $\bm{J}$ is proportional to the supercurrent:

$\displaystyle\bm{J}=-(4e\gamma/\hbar)|\psi|^{2}(\nabla\phi+2e\bm{A}/\hbar
c).$ (35)

For $\gamma>0$, we have a diamagnetic Meissner effect and the London
penetration depth is

$\displaystyle\lambda_{L}^{-2}=32\pi(e^{2}/\hbar^{2}c^{2})\gamma|\psi|^{2}.$
(36)

The superfluid density is then $n_{s}\propto\gamma|\psi|^{2}$. The above
equations demonstrate why a paramagnetic Meissner effect in a homogeneous
superconductor becomes thermodynamically unstable: to reverse the supercurrent
response due to $\bm{A}$, we must set $\gamma<0$. This makes the superfluid
density negative. The problem is now that the free energy term $\propto\gamma$
can be made increasingly negative by making $\nabla\psi$ larger and larger in
magnitude. There is no lower bound as $\psi$ becomes ever more inhomogeneous,
which is unphysical.

A bulk intrinsic odd-frequency superconductor has yet to be discovered
experimentally. In such a system, one would expect the Meissner response to be
diamagnetic in order for the superconducting phase to be thermodynamically
stable. Note that this is not required for the odd-frequency superconductivity
induced in heterostructures, such as S/F bilayers, which can exhibit
paramagnetic behavior. The coexistence of odd-frequency correlations
displaying a positive and negative superfluid density, respectively, has shown
to be unphysical [16]. Several works have shown that a diamagnetic odd-
frequency superconducting state is thermodynamically stable [21, 22], while
Ref. [16] raised concerns about the existence of a mean-field Hamiltonian
describing such a state. For a detailed discussion, see Ref. [10]. The fact
that the Meissner response can be either diamagnetic, paramagnetic, or even
absent depending on which particular ansatz one chooses for the odd-frequency
anomalous Green function $f$ emphasizes the care that must be taken when
considering effective models for odd-frequency superconductivity.

##  III Results and Discussion

There are several ways that the condition $\bm{d}=\tilde{\bm{d}}$ can be
satisfied, since the $\tilde{\ldots}$-operation involves both inversion of
energy $E\to-E$ and complex conjugation. However, choosing an anomalous Green
function which satisfies this symmetry can still lead to physically
problematic behavior, such as non-conserved spectral weight, as in Ref. [24].
Below, we give three examples which fulfill the four criteria listed
previously, yet all have different magnetic responses. For simplicity, we
consider only the $S_{z}=0$ triplet component $d_{z}$ in what follows.

###  III.1 Model I: diamagnetic Meissner response

We consider a model

$\displaystyle d_{z}=A(E)f_{\text{BCS}},$ (37)

where $A$ is a real odd function of energy that must be chosen such that Eq.
(23) holds. Using normalization, the normal Green function becomes

$\underline{g}=\sqrt{\frac{E^{2}-\Delta^{2}(1-[A(E)]^{2})}{E^{2}-|\Delta|^{2}}}\underline{1}.$
(38)

In order to calculate the Meissner response in this model, we consider the
integral from Eq. (28), starting with the trace,

$\text{Tr}\\!\left\\{{\Im\mathrm{m}}\\{\underline{f}\tilde{\underline{f}}\\}\right\\}=\frac{-2[A(E)]^{2}|\Delta|^{2}\text{sgn}\\!\left(E\right)\delta^{\prime}}{\delta^{\prime
2}+(E^{2}-|\Delta|^{2})^{2}},$ (39)

where we have used that $E\rightarrow E+\mathrm{i}\delta$ in the retarded
Green function, where $\delta\in\mathbb{R}$ is a positive infinitesimal, and
introduced $\delta^{\prime}=2|E|\delta$. Evaluating this Lorentzian as
$\delta^{\prime}\rightarrow 0$ and performing the integral, we find that the
current becomes

$\bm{j}=-[A(\Delta)]^{2}\frac{\bm{A}}{\lambda_{L}},$ (40)

where $\lambda_{L}$ is the BCS proportionality constant. This is a
conventional diamagnetic Meissner response with a renormalized London
penetration depth given in terms of the square of the antisymmetric function.
As a simple example, one could consider a function
$A(E)=\alpha\text{sgn}\\!\left(E\right)$, which for the choice $\alpha=1$
would produce the exact same Meissner response as in the BCS case. Moreover,
it can be seen that the density of states is equivalent to the BCS density of
states.

Initially, this result may seem at odds with the general result that systems
with odd-frequency pairing should have a density of states larger or equal to
one at $E=0$, which is derived in Appendix A. This apparent paradox is
understood by noting that the sign function is defined so that
$\text{sgn}\\!\left(E=0\right)=0$, ensuring $[A(E)]^{2}=0$ for $E=0$. This
causes the density of states to be nonzero in an infinitely narrow peak
centered at $E=0$, as seen from Eq. (38) when setting $A(E=0)=0$. Since this
peak is infinitely narrow, while the height of the peak is finite, the
contribution to the spectral sum rule vanishes, and we can ignore this peak in
the following. The discussion above holds also for other choices of $A(E)$
that produce a gap at the Fermi energy. This function can also be shown to
preserve the sum rule of the spectral weight for $\alpha\leq 1$. For the
special case of $\alpha=1$, the density of states have the exact same form as
for a BCS superconductor.

Figure 2: The density of states for Model I for different values of $\alpha$.

###  III.2 Model II: mixed Meissner response

Next, we study the model

$\displaystyle d_{z}=\mathrm{i}S(E)f_{\text{BCS}},$ (41)

where $S$ is a real even function of energy that must be chosen such that Eq.
(23) holds. The normal Green function of the system is given through the
normalization condition as

$\underline{g}(E)=\sqrt{\frac{E^{2}-\left(1+[S(E)]^{2}\right)\Delta^{2}}{E^{2}-\Delta^{2}}}\underline{1},$
(42)

choosing the positive root in order for the system to reduce to the normal
metal case when the order parameter vanishes.

Some general remarks can be made about $\underline{g}$ and the corresponding
density of states, even without specifying the function $S$. The density of
states is peaked at zero energy, and is larger than unity for $E<\Delta$.
Additionally, $\underline{g}$ is purely real for energies in the domains
$|E|<|\Delta|$ and $|E|>\sqrt{1+S(E)}|\Delta|$, and purely imaginary (thus
providing a gapped density of states) in between. This means that if one wants
a model without a gap in the density of states, $S$ must be chosen in such a
way that it satisfies $S(\Delta)=0$. For this model to satisfy all of our
criteria discussed above, we only need to choose a function $S(E)$ that causes
the number of states to be conserved for all choices of $\Delta$. The easiest
such choice is $S(E)=s\in\mathbb{R}$. Using elliptic integrals, it can be
shown that this choice makes the Green function satisfy the sum rule of the
spectral weight. This continues to holds for any $\Delta$, because $f_{BCS}$,
and thus $d_{z}$ in this case, depends on $\Delta$ only through the ratio
$E/\Delta$.

In Fig 3, the density of states of this model is plotted along with the BCS
density of states for comparison. We have included the effect of inelastic
scattering via a Dynes parameter [35] $E\to E+\mathrm{i}\delta$ as in the
previous section, since it is always present to some extent in actual
experiments. We have also verified numerically that adding inelastic
scattering in this way does not cause a violation of Eq. (23). As seen in Fig.
3, a very large amount of inelastic scattering $\delta\simeq\Delta$ is
required to make the density of states have a clear peak at the Fermi level
for this particular model. However, we remark that the $\delta$ appearing in
$d_{z}=\mathrm{i}S(E)f_{\text{BCS}}$ in this way does not need to correspond
to actual inelastic scattering: alternatively, $\delta$ may be considered a
fitting parameter to give the anomalous Green function the desired density of
states behavior.

By a similar calculation as in the last subsection, we find that

$\bm{j}=[S(\Delta)]^{2}\frac{\bm{A}}{\lambda_{L}},$ (43)

which is an unconventional Meissner response, producing an oscillating
susceptibility. For the special case of $S(E)=1$, this is the exact response
discussed in the last section.

Figure 3: The density of states for Model II with $S(E)=2$ when adding a Dynes
[35] inelastic scattering parameter $\epsilon\equiv\delta/\Delta$. The BCS
density of states with a Dynes parameter $\epsilon=0$ is also included for
comparison.

###  III.3 Model III: no Meissner response

For completeness, we finally include a model for the odd-frequency anomalous
Green function which does not give rise to any Meissner response. Although
unphysical, we include this for the purpose of illustrating that writing down
an ansatz for the anomalous Green function which satisfies the four criteria
outlined previously, including conservation of spectral weight, is not
sufficient to produce a Meissner response, be it diamagnetic or paramagnetic.
The ansatz reads

$\displaystyle
d_{z}=c\theta(a^{2}-E^{2})+\mathrm{i}d\text{sgn}(E)\theta(E^{2}-a^{2})\theta(b^{2}-E^{2}),$
(44)

where $a,b,c,d$ are positive real numbers. The normal Green function becomes

$\underline{g}=g\underline{1},$ (45)

with

$g=c_{2}\Theta(a^{2}-E^{2})+d_{2}\Theta(E^{2}-a^{2})\Theta(b^{2}-E^{2})+\Theta(E^{2}-b^{2}),$
(46)

where we introduced the auxiliary parameters $c_{2}\equiv\sqrt{1+c^{2}}$ and
$d_{2}\equiv\sqrt{1-d^{2}}$.

Figure 4: The density of states for Model III.

This leads to a peaked density of states, which is plotted in Fig. 4. By
fixing one of the parameters, it can be shown to satisfy the sum rule. It is,
however, readily shown that

${\Im\mathrm{m}}\\{\underline{f}\tilde{\underline{f}}\\}=0,$ (47)

which causes the Meissner response to vanish for this choice of effective odd-
frequency anomalous Green function.

###  III.4 Model IV: numerical solution in full proximity effect regime

The models considered above (I-III) give simple expressions for the density of
states and Meissner response.

(a) The real and imaginary part of the anomalous Green function for Model IV
at the vacuum edge.

(b) The normalized density of states for Model IV at the vacuum edge.

Figure 5: The proximity effect for Model IV at the vacuum edge. Parameters
used are $h=50\Delta$, $\Omega=6.7\xi^{-1}$, $L=0.75\xi$, and $\epsilon=0.01$.

It is instructive to compare these with the properties of a system where odd-
frequency superconductivity arises naturally through the proximity effect.
Here, we consider an S/F/F/N system, modelled with a BCS equilibrium
superconductor in contact with a single ferromagnet with a varying
magnetization through an interface with a boundary specific constant $\Omega$.
The origin ($x=0$) is at the S/F interface. A natural length scale in dirty
superconducting systems is the superconducting coherence length
$\xi=\sqrt{D/\Delta}$, which we use to make the Usadel equation dimensionless.
Lengths are given in units of $\xi$, and energies in units of $\Delta$. The
solutions for the quasiclassical Green functions were found using the Riccati
parametrization [36, 45]. The magnetization close to the interface is pointing
in the z-direction, and, over a short distance, rotates with a uniform rate
into pointing in the x-direction. Throughout the ferromagnet, the spin-
splitting $h$ due to the magnetization has a constant magnitude. At the F/N
interface, the spin-splitting is set to decay exponentially over a short
distance. Tuning the parameters allows us to make a system where the triplet
amplitude dominates in the normal metal part of the material. In this way, the
normal metal becomes an effective odd-frequency superconductor. We consider
the vacuum end of the normal metal, at $x=L$, and plot the real and imaginary
part of the triplet component $d_{z}$ in Fig. 5(a), as well as the density of
states in Fig. 5(b). As seen, the energy-dependence of the real and imaginary
part is not easily described by a simple trial function.

Figure 6: A qualitative fit of Model II (dashed lines) with a constant
function $S(E)$ to the results from the numerically exact solution of the
S/F/F/N junction (Model IV) (solid lines), with the same parameters as in Fig
5. The parameters $\\{S(E),\epsilon_{f},\Delta_{f}\\}$ used in Model II, to be
distinguished from the parameters $\\{\epsilon,\Delta\\}$ used in the
numerical solution of Model IV, to make the fit are $S(E)=0.88$,
$\epsilon_{f}=0.19$, and $\Delta_{f}=0.25\Delta$.

Moreover, we calculated the Meissner response in the S/F/F/N proximity system
and found that close to the superconductor, where singlet components dominate,
the Meissner response is conventional. In the N part of the system, where odd-
frequency pairing dominates, the Meissner response is unconventional. The
peaked density of states at the Fermi energy as well as the sign of the
Meissner response suggests that the function could be fitted with Model II
above. However, the non-trivial form of the solution shows that a complicated
function $S(E)$ must be chosen in order to get a very good quantitative fit.
Nevertheless, a qualitative and quantitatively comparable fit can be obtained
by simply choosing $S(E)$ to be a constant as well as letting the Dynes
parameter take the value $\epsilon_{f}$ and the order parameter take the value
$\Delta_{f}$ as curve fit parameters, which in general can take different
values than the corresponding quantities $\epsilon,\Delta$ in the actual
proximity system. Using Eq. (41), a fit of the anomalous Green function and
the density of states for the S/F/F/N proximity system are shown in Fig. 6.
Numerically, it can be verified that the density of states satisfy Eq. (23).
It is interesting to note the compact form of the effective odd-frequency
anomalous Green function $d_{z}$ in Eq. (41) used to fit the actual numerical
solution of the S/F/F/N system: it is simply proportional to the BCS anomalous
Green function, which makes it useful to work with practically.

##  IV Summary

In this work, we establish a set of necessary criteria that a physically sound
quasiclassical Green functions describing odd-frequency superconductivity in
the dirty limit must satisfy. We demonstrate that it is possible to write down
effective models for the Green function which satisfy all the abovementioned
criteria, but which nevertheless give rise to different behavior when it comes
to the density of states and the magnetic response of the superconductor. This
highlights the importance of carefully considering the properties of effective
models describing odd-frequency superconductivity and may provide a guide for
future works considering emergent phenomena associated with odd-frequency
superconductors.

###### Acknowledgements.

We thank S. Aunsmo, J. A. Ouassou, and L. G. Johnsen for useful discussions
and A. Golubov for correspondence. We acknowledge funding via the Research
Council of Norway Grant numbers 323766, as well as through its Centres of
Excellence funding scheme, project number 262633. J. L. also acknowledges
support from the Sigma2 project no. NN9577K.

##  Appendix A Proof of increased (decreased) density of states at $E=0$ for
triplets (singlets)

In this appendix, we will derive a relation describing how the symmetries of
the quasiclassical Green function determine the effect of the density of
states at the Fermi energy. Specifically, we will show that for triplet odd-
frequency pairing, the density of states must be larger or equal to one, while
for singlet even-frequency, it must be lowered. Consider a system with both
singlet and $d_{z}$ triplet pairing, described by

$\underline{f}=\begin{pmatrix}0&d_{z}+f_{s}\\\ d_{z}-f_{s}&0\end{pmatrix},$
(48)

and consider $E=0$, meaning that $\tilde{f}(0)=[f(0)]^{*}$. The proofs for the
other triplet components are completely analogous. From the normalization
condition we get, temporarily suppressing the $E=0$ argument, and introducing
the quantity $c\equiv fs^{*}d_{z}-f_{s}d_{z}^{*}$,

$\underline{g}^{2}=\begin{pmatrix}1+|d_{z}|^{2}-|f_{s}|^{2}+2\mathrm{i}{\Im\mathrm{m}}\\{c\\}&0\\\
0&1+|d_{z}|^{2}-|f_{s}|^{2}-2\mathrm{i}{\Im\mathrm{m}}\\{c\\}\end{pmatrix}.$
(49)

Taking the positive square root of this equation yields an expression for the
normal part of the quasiclassical Green function. In order to get the spin
averaged density of states, we average over the spin-up and spin-down pairing
density, and take the real part of the quasiclassical Green function,

$N(E=0)=\frac{N_{0}}{2}\text{Re}\\{g_{\uparrow\uparrow}+g_{\downarrow\downarrow}\\}.$
(50)

In the case of a weak proximity effect, we expand the square root to first
order, and get

$N(E=0)=N_{0}\left(1+\frac{|d_{z}|^{2}-|f_{s}|^{2}}{2}\right).$ (51)

This shows that in the weak regime, singlets will always lower the zero-energy
density of states, while triplets will increase it. This property was also
noted in Ref. [37] in the weak proximity effect regime.

In the case of a strong proximity effect, we cannot expand the square root,
and the expression will not reduce in the same way as in the weak proximity
regime. If we, however, consider systems with only singlets or only triplets,
we find that

$N(E=0)=\begin{cases}N_{0}\sqrt{1-|f_{s}|^{2}},&\text{ even-frequency,}\\\
N_{0}\sqrt{1+|d_{z}|^{2}},&\text{ odd-frequency}.\end{cases}$ (52)

We thus see that also in the full proximity effect regime, the triplets
increase the density of states at zero energy, while the singlets decrease it.

## References

  * [1] J. Bardeen, L. N. Cooper, and J. R. Schrieffer Phys. Rev. 108, 1175 (1957).
  * [2] V. L. Berezinskii, Pis’ma Zh. Eksp. Teor. Fiz 20, 628–631 (1974).
  * [3] D. Belitz and T. R. Kirkpatrick, Phys. Rev. B 46, 8393 (1992).
  * [4] A. Balatsky and E. Abrahams, Phys. Rev. B 45, 13125 (1992)
  * [5] T. R. Kirkpatrick and D. Belitz, Phys. Rev. Lett. 66, 1533 (1991).
  * [6] P. Coleman, E. Miranda, and A. Tsvelik, Phys. Rev. Lett. 70, 2960 (1993).
  * [7] F. S. Bergeret, A. F. Volkov, and K. B. Efetov, Phys. Rev. Lett. 86, 4096–4099 (2001).
  * [8] F. S. Bergeret, A. F. Volkov, and K. B. Efetov Rev. Mod. Phys. 77, 1321 (2005).
  * [9] A. I. Buzdin Rev. Mod. Phys. 77, 935 (2005).
  * [10] J. Linder and A. V. Balatsky, Rev. Mod. Phys. 91, 045005 (2019).
  * [11] R. S. Keizer, S. T. B. Goennenwein, T. M. Klapwijk, G. Miao, G. Xiao, and A. Gupta Nature 439, 825 (2006).
  * [12] Trupti S. Khaire, Mazin A. Khasawneh, W. P. Pratt, Jr., and Norman O. Birge Phys. Rev. Lett. 104, 137002 (2010).
  * [13] J. W. A. Robinson, J. D. S. Witt, and M. G. Blamire, Science 329, 59 (2010).
  * [14] M. Eschrig, Reports on Progress in Physics 78, 104501 (2015).
  * [15] A. Di Bernardo, Z. Salman, X. L. Wang, M. Amado, M. Egilmez, M. G. Flokstra, A. Suter, S. L. Lee, J. H. Zhao, T. Prokscha, E. Morenzoni, M. G. Blamire, J. Linder, and J. W. A. Robinson, Phys. Rev. X 5, 041021 (2015).
  * [16] Ya. V. Fominov, Y. Tanaka, Y. Asano, and M. Eschrig, Phys. Rev. B 91, 144514 (2015).
  * [17] T. Kontos, M. Aprili, J. Lesueur, and X. Grison, Phys. Rev. Lett. 86, 304 (2001).
  * [18] F. S. Bergeret, A. F. Volkov, K. B. Efetov, Phys. Rev. B 65, 134505 (2002).
  * [19] Y. Asano and Y. Tanaka, Phys. Rev. B 87, 104513 (2013).
  * [20] Z. Huang, P. Wolfle, A. V. Balatsky, Phys. Rev. B 92, 121404 (2015).
  * [21] D. Solenov, I. Martin, and D. Mozyrsky, Phys. Rev. B 79, 132502 (2009).
  * [22] H. Kusunose, Y. Fuseya, and K. Miyake, Journal of the Physical Society of Japan 80, 054702 (2011).
  * [23] Y. Fominov, JETP Letters 86, 732 (2007).
  * [24] Y. Tanaka and A. A. Golubov, Phys. Rev. Lett. 98, 037003 (2007).
  * [25] P. O. Sukhachov and A. V. Balatsky, Phys. Rev. B 100, 134516 (2019).
  * [26] L. G. Johnsen and J. Linder, Phys. Rev. B 104, 144513 (2021).
  * [27] G. Eilenberger, Z. Physik 214, 195 (1968).
  * [28] K. D. Usadel, Phys. Rev. Lett. 25, 507 (1970).
  * [29] W. Belzig, F. K. Wilhelm, C. Bruder, G. Schön, and A. D. Zaikin, Superlattices Microstruct. 25, 1251 (1999).
  * [30] Y. Asano, Y. Fominov, and Y. Tanaka, Phys. Rev. B 90, 094512 (2014).
  * [31] Y. Nambu, Phys. Rev. 117, 3 (1960).
  * [32] M. Y. Kuprianov and V. F. Lukichev, Soviet Physics-JETP, 67, 1163 (1988).
  * [33] T. Yokoyama, Y. Tanaka, and N. Nagaosa, Phys. Rev. Lett., 106, 246601 (2011).
  * [34] S. Jacobsen, J. A. Ouassou, and J. Linder, Phys. Rev. B 92, 024510 (2015).
  * [35] R. C. Dynes, V. Narayanamurti, and J. P. Garno, Phys. Rev. Lett. 41, 1509 (1978).
  * [36] N. Schopohl and K. Maki, Phys. Rev. B 52, 490 (1995).
  * [37] S. Kawabata, Y. Asano, Y. Tanaka, A. A. Golubov, J. Phys. Soc. Jpn. 82, 124702 (2013).
  * [38] M. Eschrig, T. Löfwander, T. Champel, J. C. Cuevas, J. Kopu, and G. Schön, J. Low. Temp. Phys. 147, 457 (2007).
  * [39] Y. Tanaka, A. Golubov, S. Kashiwaya, and M. Ueda, Phys. Rev. Lett. 99, 037005 (2007).
  * [40] Yu. V. Nazarov, Superlatt. and Microstruc. 25, 1221 (1999).
  * [41] A. Cottet, D. Huertas-Hernando, W. Belzig, and Yu. V. Nazarov, Phys. Rev. B 80, 184511 (2009); Erratum: Phys. Rev. B 83, 139901 (2011).
  * [42] M. Eschrig, A. Cottet, W. Belzig, and J. Linder, New J. Phys. 17, 083037 (2015).
  * [43] M. Amundsen and J. Linder, Phys. Rev. B 100, 064502 (2019); Erratum: Phys. Rev. B 103, 059902 (2021).
  * [44] J. Linder and M. Amundsen, Phys. Rev. B 105, 064506 (2022).
  * [45] A. Konstandin, J. Kopu, and M. Eschrig Phys. Rev. B 72, 140501(R) (2005).