# Measurement-induced phase transition in a single-body tight-binding model

Tony Jin  Pritzker School of Molecular Engineering, The University of Chicago,
Chicago, IL 60637, USA  Université Côte d’Azur, CNRS, Institut de Physique de
Nice, 06200 Nice, France  David G. Martin  Kadanoff Center for Theoretical
Physics and Enrico Fermi Institute, 933 E 56th St, The University of Chicago,
Chicago, IL 60637, USA

###### Abstract

We study the statistical properties of a single free quantum particle evolving
coherently on a discrete lattice in $d$ spatial dimensions where every lattice
site is additionally subject to continuous measurement of the occupation
number. Using perturbative renormalization group (RG) analysis, we show that
the systems undergoes a Measurement-induced Phase Transition (MiPT) for $d>2$
from a delocalized to a localized phase as the measurement strength $\gamma$
is increased beyond a critical value $\gamma_{c}$. In the language of surface
growth, the delocalized phase corresponds to a smooth phase while the
localized phase corresponds to a rough phase. We support our analytical
computations with numerical analysis which are in qualitative and quantitative
agreement with the theory.

Recently, it has been discovered that quantum chaotic systems subject to
continuous or projective measurements could undergo a phase transition
characterized by a change of the scaling properties of the entanglement
entropy with time or system size, a phenomenon now referred to as Measurement-
induced Phase Transition (MiPT) [1]. MiPT constitutes a fascinating problem at
the crossroad of statistical physics and quantum information. As such, it has
attracted a tremendeous amount of interest in the recent years [2, 3, 4, 5, 6,
7, 8, 9, 10, 11, 12, 13]. MiPTs are often characterized by a transition from
an _area law_ phase, i.e a phase where the entanglement entropy (EE) of a
subsystem doesn’t scale with its size, to a _volume law_ phase, i.e a phase
where the EE scales with the system size. Such a scaling transition occurs
upon increasing the strength of the measurement and is observed in various
systems such as in 1d interacting chaotic many-body systems.

However, surprisingly, the existence of a MiPT between two non-trivially
entangled phases for free or Gaussian fermions undergoing measurements remains
an actively debated question. While the original study of entanglement in 1d
free fermions [14] showed no signs of a phase transition, more recent
numerical and theoretical investigations showed either the existence of a
phase where the EE scales as $\log L$ [15, 16, 17] or $\left(\log
L\right)^{2}$ [18] while another recent study [19] argued that the observed
transitions are in fact sharp crossovers.

In this work, we provide additional insights on this subject by studying the
simpler, yet non trivial single-body problem of a particle evolving coherently
on a discrete lattice in $d$ spatial dimensions, where every lattice site is
subject to independent, continuous measurements of its occupation number–See
Fig.1. Combining perturbative renormalization group (RG) methods and numerical
simulations, we show that, while we do not find evidence of a transition in
$d=1$, there exists a phase transition from a _smooth/delocalized_ phase to a
_rough/localized_ phase when $d>2$. Interestingly, this shows that many-body
effects are not necessary to observe a MiPT and corroborates the result
obtained in [20] for a classical random walker undergoing continuous
measurements.

xyz continuous measurement diffusionwalker

smooth/delocalizedrough/localized$\boldsymbol{\gamma}$$\boldsymbol{\gamma_{c}}$

Figure 1: Our model consists of a single quantum walker on a $d$-dimensional
square lattice of $N^{d}$ sites undergoing a unitary tight-binding evolution
in addition to independent continuous measurements of the occupation number at
every lattice site. In $d>2$, we observe a phase transition from a
smooth/delocalized phase to a rough/localized phase upon increasing the
measurement strength $\gamma$ beyond a critical value $\gamma_{c}$. The
snapshots show a typical density profile in each of these phases for $d=3$,
where one of the spatial direction has been projected out. Parameters:
$\tau=1.5$, $N=41$, $dt=0.01$. For the smooth phase $\gamma=0.9$ and for the
rough phase $\gamma=3$.

##### Model

We consider a single quantum particle on a square lattice of $N^{d}$ sites
with periodic boundary conditions. Let
$\\{\left|\boldsymbol{j}\right\rangle\\}_{\boldsymbol{j}\in[1,N]^{d}}$ denote
the position basis. The dynamics is described by a unitary tight-binding
evolution
$\hat{H}:=-\tau\sum_{\\{|\boldsymbol{e}|=1\\}}\left|\boldsymbol{j}\right\rangle\left\langle\boldsymbol{j}+\boldsymbol{e}\right|$
where $\\{|\boldsymbol{e}|=1\\}$ is the set of vectors of norm $1$. In
addition, each site undergoes continuous measurements of strength $\gamma$ of
the local occupation
$\hat{n}_{\boldsymbol{j}}:=\left|\boldsymbol{j}\right\rangle\left\langle\boldsymbol{j}\right|$
resulting in the stochastic differential equation (SDE) [21]:

$\displaystyle d\left|\psi\right\rangle=-iH\left|\psi\right\rangle dt$ (1)
$\displaystyle+\sum_{\boldsymbol{j}}\left(-\frac{\gamma}{2}(\hat{n}_{\boldsymbol{j}}-\langle\hat{n}_{\boldsymbol{j}}\rangle_{t})^{2}dt+\sqrt{\gamma}(\hat{n}_{\boldsymbol{j}}-\langle\hat{n}_{\boldsymbol{j}}\rangle_{t})dB_{t}^{\boldsymbol{j}}\right)\left|\psi\right\rangle\;,$

where $\langle\bullet\rangle_{t}:={\rm tr}(\rho_{t}\bullet)$. In (1), the
$\\{B_{t}^{\boldsymbol{j}}\\}_{\boldsymbol{j}\in[1,N]^{d}}$ are $N^{d}$
independent Brownian processes with average
$\mathbb{E}[dB_{t}^{\boldsymbol{j}}]=0$ and Itō rules
$dB_{t}^{\boldsymbol{j}}dB_{t^{\prime}}^{\boldsymbol{k}}=\boldsymbol{1}_{0}(t-t^{\prime})\delta_{\boldsymbol{j},\boldsymbol{k}}dt$
where $\boldsymbol{1}_{0}$ is the indicator function. This model was
originally introduced in [22] for the free fermionic case and has been
subsequently studied in [14, 15, 16] in the context of MiPTs – see also [23,
24] for applications to transport and thermal engines.

In terms of the basis elements $\psi_{\boldsymbol{j}}$ defined as
$\left|\psi\right\rangle=\sum_{\boldsymbol{j}}\psi_{\boldsymbol{j}}\left|\boldsymbol{j}\right\rangle$,
Eq.(1) can be written as

$\displaystyle d\psi_{\boldsymbol{j}}=$ $\displaystyle
i\tau\sum_{\\{|\boldsymbol{e}|=1\\}}\psi_{\boldsymbol{j}+\boldsymbol{e}}dt-\frac{\gamma}{2}\psi_{\boldsymbol{j}}\big{(}1-2|\psi_{\boldsymbol{j}}|^{2}+\sum_{\boldsymbol{m}}|\psi_{\boldsymbol{m}}|^{4}\big{)}dt+\sqrt{\gamma}\psi_{\boldsymbol{j}}\big{(}dB_{t}^{\boldsymbol{j}}-\sum_{\boldsymbol{m}}|\psi_{\boldsymbol{m}}|^{2}dB_{t}^{\boldsymbol{m}}\big{)}.$
(2)

Throughout the rest of the manuscript, we will fix the initial condition of
the system to be $\psi_{\boldsymbol{j}}(t=0)=N^{-d/2}$ for all
$\boldsymbol{j}$. Note that by construction,
$\sum_{\boldsymbol{j}}|\psi_{\boldsymbol{j}}|^{2}$ is preserved on every
realization of the noise.

Even though (2) describes the dynamics of a single particle, getting an exact
solution of such a SDE is in general a formidable task. One way to make
progress is to restrict it to what we will refer to as the _delocalized_
phase, i.e to assume that $|\psi_{\boldsymbol{j}}|$ is of order $N^{-d/2}$.
Under this assumption, keeping the leading order in $N^{-1}$ in (2) gives the
simpler expression

$d\psi_{\boldsymbol{j}}=\bigg{(}i\tau\sum_{\\{|\boldsymbol{e}|=1\\}}\psi_{\boldsymbol{j}+\boldsymbol{e}}-\frac{\gamma}{2}\psi_{\boldsymbol{j}}\bigg{)}dt+\sqrt{\gamma}\psi_{\boldsymbol{j}}dB_{t}^{\boldsymbol{j}}$
(3)

which is now _local_. We further take the continuous limit by introducing the
lattice spacing $b$ and the continuous quantities $\vec{r}=\boldsymbol{j}b$,
$\varphi(\vec{r}=\boldsymbol{j}b)=b^{-d/2}\psi_{\boldsymbol{j}}$,
$d\eta(\vec{r},t):=b^{-d/2}dB_{t}^{\boldsymbol{j}}$, $D:=b^{2}\tau$,
$\lambda:=\gamma b^{d}$. Up to a global phase, (3) then becomes

$d\varphi=\left(iD\nabla^{2}\varphi-\frac{\gamma}{2}\varphi\right)dt+\sqrt{\lambda}\varphi
d\eta.$ (4)

It is important to note that the noise becomes multiplicative in (4). This
allows us to draw an analogy between (4) and the Stochastic Heat Equation
(SHE), thereby relating (3) to KPZ physics [25, 26]. Such an analogy was
already fruitfully exploited in [20], where it led to an intuitive
understanding of a MiPT in a classical context. The difference with this
previous study is that we deal with an _imaginary_ diffusion term $D$ as well
as a real “mass” term $\gamma/2$. Even though these differences lead to
quantitative modifications with respect to the SHE, we will see that one of
its main feature, namely the existence of a phase transition with lower
critical dimension $2$, remains present in (4).

We expect (3) to be valid as long as $|\psi_{\boldsymbol{j}}|$ remains close
to the homogeneous profile of order $N^{-d/2}$. Such an assumption is verified
when the renormalization flow is directed towards the delocalized phase: in
this case, $|\psi_{\boldsymbol{j}}|$ is indeed driven closer to the
homogeneous profile. On the opposite, if the renormalization flow is directed
towards the localized phase, $|\psi_{\boldsymbol{j}}|$ is driven away from the
homogeneous profile and we don’t expect (3) to remain valid at long time.
However, although (3) does not describe the strongly localized regime, it
still allows us to make quantitative assertions concerning the boundary region
between the two phases.

##### Martin-Siggia-Rose (MSR) action

We now proceed to derive the MSR action [27, 28] $Z$ associated to (4). The
details of the derivation are presented in the SM [29]. Let the superscript
$a$ denote the auxiliary fields. We have that

$Z=\int{\cal
D}[\varphi,\bar{\varphi},\varphi^{a},\bar{\varphi}^{a}]e^{iS_{0}+iS_{\nu}},$
(5)

where the bar denotes complex conjugation, and $S_{0}$, $S_{\nu}$ are
respectively the quadratic and quartic part of the action:

$\displaystyle S_{0}$ $\displaystyle=\int
d^{d}\vec{r}dt(\bar{\varphi},\bar{\varphi}^{a})\boldsymbol{G}_{0}^{-1}\begin{pmatrix}\varphi\\\
\varphi^{a}\end{pmatrix},$ (6) $\displaystyle\boldsymbol{G}_{0}^{-1}$
$\displaystyle=\frac{1}{2}\begin{pmatrix}0&-\frac{\gamma}{2}+\partial_{t}-iD\nabla^{2}\\\
-\frac{\gamma}{2}-\partial_{t}+iD\nabla^{2}&0\end{pmatrix},$ (7)
$S_{\nu}=\frac{i}{8}\int d^{d}\vec{r}dt\big{(}\lambda^{{\rm
I}}\left(\bar{\varphi}^{a}\right)^{2}\varphi^{2}+\lambda^{{\rm
II}}\bar{\varphi}^{a}\bar{\varphi}\varphi^{a}\varphi+{\rm c.c}\big{)}.$ (8)

In (8), we introduced the labels I, II for the interacting terms, as they will
behave differently under renormalization. For the microscopic theory (3), we
have $\lambda^{{\rm I}}=\lambda^{{\rm II}}=\lambda=\gamma b^{d}$.

Inverting (7) yields the free propagator in momentum $\vec{q}$ and frequency
$\omega$:

$\displaystyle G_{0}^{R}(\vec{q},\omega)$
$\displaystyle=\frac{2i}{Dq^{2}-\omega-i\frac{\gamma}{2}}$ (9)

where the $R$ label refers to the retarded propagator. The advanced propagator
$A$ is obtained by complex conjugation
$G_{0}^{A}(\vec{q},\omega)=\bar{G}_{0}^{R}(\vec{q},\omega)$.

$\lambda^{\rm{I}}$$G^{R}(q,w)$$\lambda^{\rm{II}}$$G^{A}(q,w)$MassInteractions$\lambda^{\rm{I}}$$\lambda^{\rm{I}}$$G^{R}(q,w)$$G^{R}(-q,-w)$$\lambda^{\rm{II}}$$\lambda^{\rm{II}}$$G^{R}(q,w)$$G^{A}(q,w)$
Figure 2: One-loop contributions to the RG of the mass and interactions.
Dashed lines designate the auxiliary field $\varphi^{a}$ and full lines the
field $\varphi$.

##### Renormalization flow

We proceed to the one-loop perturbative renormalization group (RG) analysis of
(4). We employ standard momentum-shell analysis [28]. Let $\Lambda$ be the
microscopic momentum cut-off of the theory. The critical exponents associated
to $t$, $\varphi$, $\bar{\varphi}$ are named respectively $z$, $\chi$ and
$\bar{\chi}$. The flow is parametrized by $l$.

abc

Figure 3: a Phase diagram for $d=3$. Each voxel corresponds to a set of
initial conditions $\\{\gamma_{R},\lambda_{R}^{{\rm I}},\lambda^{{\rm
II}}\\}_{l=0}$ with $\lambda_{I}^{{\rm I}}(0)=\gamma_{I}(0)=0$. The blue
voxels indicate the smooth/delocalized phase and all the other points are in
the rough/localized phase. The red voxels line is the diagonal
$\gamma_{R}(0)=\lambda_{R}^{{\rm I}}(0)=\lambda^{{\rm II}}(0)$. The
intersection between the red line and the blue domain indicates the value
$\gamma_{c}$ where the phase transition occurs. b Plots of different
trajectories in the $\\{\gamma,\lambda^{{\rm II}}\\}$ plane for the simpler
case (12) in $d=1$ and in $d=3$. We see that for $d=1$ all the trajectories
eventually bend upwards, indicating the existence of a single phase. On the
opposite, in $d=3$ there exists a critical line separating a phase where
$\lambda^{{\rm II}}\to\infty$ from a phase where $\lambda^{{\rm II}}\to 0$. c
Dependence of $\gamma_{c}$ on $\tau$. The renormalization flow equations (11)
predict a simple straight line with slope $1$, which is shown in blue (passing
by the point {$\gamma_{c}=1.14$, $\tau=1$}). The red dots correspond to values
of $\gamma_{c}$ extracted from numerical simulations of (2). They can be
fitted by a slope of $\approx 1.19$ with correlation coefficient $\approx
0.999$.

At the one-loop level, there are no diagrams renormalizing the part of the
action proportional to $\partial_{t}$ and $\nabla^{2}$. Imposing the
stationarity of the corresponding prefactors under the flow gives

$\chi+\bar{\chi}+d=0,\quad z=2.$ (10)

The one-loop contribution to the renormalization flow of the "mass" term
$\gamma$ and the interaction $\lambda^{\rm I,II}$ are depicted on Fig.2.
Computing their contribution to the renormalization flow leads to [29]:

$\displaystyle\frac{d\gamma}{dl}=$ $\displaystyle 2\gamma-{\rm
sgn}\left(\gamma_{R}\right)K_{d}\left(2\lambda^{{\rm I}}+\lambda^{{\rm
II}}\right),$ $\displaystyle\frac{d\lambda^{{\rm I}}}{dl}=$
$\displaystyle(2-d)\lambda^{{\rm I}}+{\rm
sgn}\left(\gamma_{R}\right)K_{d}\frac{\left(\lambda^{{\rm
I}}\right)^{2}}{\gamma+2iD\Lambda^{2}},$ (11)
$\displaystyle\frac{d\lambda^{{\rm II}}}{dl}=$
$\displaystyle(2-d)\lambda^{{\rm II}}+{\rm
sgn}\left(\gamma_{R}\right)K_{d}\frac{\left(\lambda^{{\rm
II}}\right)^{2}}{\gamma_{R}}.$

where ${\rm sgn}$ is the sign function, $\gamma_{R}:=\Re(\gamma)$ and we
introduced $K_{d}:=\frac{\Lambda^{d}}{\Gamma(d/2)2^{d-1}\pi^{d/2}}$ with
$\Gamma$ the Euler function.

The qualitative properties of the phase diagram can be understood by
considering the simpler case of $\lambda^{{\rm I}}=0$, $\gamma=\gamma_{R}$ and
restricting the study to the domain $\gamma\geq 0$. In this case, the flow
equations can be integrated exactly to yield

$\displaystyle\gamma$
$\displaystyle=\gamma_{0}e^{(2-\frac{d}{2})l}\sqrt{\frac{1+ce^{dl}}{1+c}},\quad\lambda=\lambda_{0}e^{(2-\frac{d}{2})l}\sqrt{\frac{1+c}{1+ce^{dl}}},$
(12)

where $\gamma_{0}:=\gamma(l=0)$, $\lambda_{0}:=\lambda(l=0)$ and
$c:=\frac{d\gamma_{0}}{2K_{d}\lambda_{0}}-1$. The asymptotic behavior of these
equations depends on the sign of $c$, i.e on the dimension and the initial
ratio $\gamma_{0}/\lambda_{0}$. For $c<0$ we see that $\gamma\to 0$ and
$\lambda\to\infty$ as $l\to\frac{1}{d}\log(-c^{-1})$. For $c>0$, the
asymptotic behavior at large $l$ is given by $\gamma\propto e^{2l}$,
$\lambda\propto e^{(2-d)l}$. For $c=0$, a quick inspection of (12) shows that
the flow remains on the line of fixed $\gamma/\lambda$ with
$\gamma\propto\lambda\propto e^{\left(2-\frac{d}{2}\right)l}$. By analogy with
the physics of surface growth processes [30], we will define the system to be
_rough_ whenever $\lambda\to\infty$ and _smooth_ whenever $\lambda\to 0$.
Building on our analysis of (12), we observe that the existence of the smooth
phase is only possible when $d>2$ and depends on the initial ratio
$\gamma_{0}/\lambda_{0}$. As highlighted in Fig.1, we can interpret the rough
phase as a _localized_ phase where measurements dominate and the smooth phase
as a _delocalized_ phase where diffusive transport dominates.

We show on Fig3.b the simplified RG flows (12) in 1d (left) and 3d (right)
with the value of the microscopic cutoff $\Lambda$ fixed to $\pi/1.226595$
111We are, a priori, free to fix the microscopic cut-off $\Lambda$. To do so,
we extract the numerical value of the microscopic critical $\gamma_{c}$ for
$\tau=1$ by simulating the full microscopic evolution (2). We found
$\gamma_{c}\approx 1.14$. We then fix the value of $\Lambda$ so that the RG
flow equations (11) yield the same value of $\gamma_{c}$. Doing so, we fix
$\Lambda$ to $\pi/1.226595$. Because of the strong linearities of the flow
equations, this value has to be finely tuned.. We see that in the 1d case, all
the lines eventually lead to a divergent $\lambda^{\rm II}$ while in the 3d
case, we have a critical line separating the smooth/delocalized phase where
$\lambda^{\rm II}\to 0$ from the rough/localized phase where $\lambda^{\rm
II}\to\infty$.

To obtain the full phase diagram beyond (12), we numerically integrate the
complete flow equations (11). In Fig.3.a, we explore the 3 dimensional
parameter space $\\{\gamma_{R},\lambda_{R}^{{\rm I}},\lambda^{{\rm
II}}\\}_{l=0}$ while keeping the last two other initial conditions fixed, i.e.
$\lambda_{I}^{{\rm I}}(0)=\gamma_{I}(0)=0$. A given point in this space
belongs to the localized phase if $\lambda^{{\rm II}}\to\infty$ when starting
from this point and conversely to the delocalized phase if $\lambda^{{\rm
II}}\to 0$ (indicated by blue voxels on Fig.3.a). In the microscopic theory
(2), there is initially one parameter beside the diffusion coefficient and we
have that $\gamma_{R}=\lambda_{R}^{{\rm I}}=\lambda^{{\rm II}}=\gamma$.
Therefore, the microscopic theory (2) starts on the initial line
$\gamma_{R}(0)=\lambda_{R}^{{\rm I}}(0)=\lambda^{{\rm II}}(0)$ (red dashed
line) in the parameter space of Fig.3.a. The critical value $\gamma_{c}$ for
the microscopic parameter corresponds to the point where the line
$\gamma_{R}(0)=\lambda_{R}^{{\rm I}}(0)=\lambda^{{\rm II}}(0)$ crosses from
one domain to the other - see Fig.3.a. Note that we characterize the
transition with $\lambda^{{\rm II}}$ because, starting from the same initial
conditions, we necessarily have $|\lambda^{{\rm I}}|\leq\lambda^{{\rm II}}$
and thus, $\lambda^{{\rm II}}$ will always diverge first.

##### Numerical results

In this section, we provide numerical simulations of the complete microscopic
equations (2) to confirm the previous discussion. In order to characterize the
different phases, we introduce the height $h_{j}$ as

$h_{j}:=\frac{1}{\sqrt{\gamma}}\log\left(|\psi_{j}|^{2}\right).$ (13)

Drawing on the analogy with the classical case [20], we expect that the width
$w$ will follow a Family-Vicsek [30, 32] type scaling according to

$w:=\sqrt{\frac{1}{N}\sum_{j}(h_{j}-\langle h\rangle_{\rm s})^{2}}\propto
N^{\alpha}f\left(\frac{t}{N^{\alpha/\beta}}\right),$ (14)

where the bracket denotes the spatial average $\langle h\rangle_{\rm
s}=\frac{1}{N}\sum_{j=1}^{N}h_{j}$ and $f(x)\propto x^{\beta}$ for $x\ll 1$
while $f(x)\propto 1$ for $x\gg 1$. The universal exponents $\alpha$ and
$\beta$ characterize the dynamical phases and are respectively called the
roughness and growth exponents.

We show on Fig.4 the dependence of $\alpha$ as a function of $\gamma$ in $d=1$
and $d=3$ for different system sizes. In 1d, we see that all the curves
collapse to the value $\alpha\approx 1$ indicating a rough/localized phase
while in 3d, we see a clear crossing of the curves at a finite value of
$\gamma$ indicating a phase transition from a smooth/delocalized phase with
$\alpha\approx 0$ to a rough/localized phase with $\alpha\approx 0.32$.

Finally, we note that one simple yet non trivial prediction from (11) is that
the flow equations are invariant under the rescaling of $D,\lambda^{{\rm
I,II}}$ and $\gamma$ by the same multiplicative factor. Thus, plotting the
dependence of $\gamma_{c}$ with respect to $D$ should simply give a straight
line with slope 1. We show on Fig.3.c the values of $\gamma_{c}$ computed from
the microscopic dynamics (2) (red dots) compared to the straight line (blue)
with slope $1$ passing by {$\gamma_{c}=1.14$, $\tau=1$}. A linear regression
of the numerical data gives a slope of $\approx 1.19$ with correlation
coefficient $\approx 0.999$, slightly above the expected value. This
discrepancy could be due to higher order terms that we neglected in the
perturbative RG.

Figure 4: Critical exponent $\alpha$ as a function of $\gamma$ for $d=1$
(left) and $d=3$ (right) as a function of the measurement strength $\gamma$
obtained in numerical simulations of (2). The different curves correspond to
different system sizes. In 3d, we see the existence of a delocalized and
localized phase while there is only a localized phase for $d=1$. Parameters:
$\tau=1$, $dt=0.01$.

##### Conclusion

We have provided numerical and analytical arguments showing the existence of a
MiPT for a single free particle undergoing continuous measurement for $d=3$
and its absence for $d=1$. Our work is one of the first to demonstrate the
critical role played by dimensionality to observe the existence of a
transition or lack thereof in a quantum setting. Compared to previous studies
in the literature, it is remarkable that many-body effects play no role in the
emergence of this transition.

Our studies raise a number of interesting questions. We first note that our
characterization of the transition differs from more conventional approaches
in the literature: we do not compute the temporal or spatial scaling of EE. It
would be valuable in subsequent studies to look at the many-body case where
the EE could be defined and characterized in order to connect better with the
previously known phenomenology.

One exciting possibility is that our transition is already visible at the
level of quantities linear in the density matrix, for instance transport-
related quantities. In 1d, it is known that measurements induce a crossover
from a ballistic to diffusive transport [33, 34, 35, 36, 37] in free fermionic
system. If we associate the ballistic behavior to the delocalized phase and
the diffusive behavior to the localized one, a tempting conjecture is that
this crossover in 1d becomes a genuine phase transition in higher spatial
dimensions. Having this other route for characterizing MiPT would be
particularly interesting in the context of experiments, since the measurement
of EE is computationally very heavy, often requiring tomography of the full
quantum trajectories and/or costly post-selection procedure [38, 39].

On the theoretical side, we note that a recent interesting body of literature
has proposed non-linear sigma models as good effective descriptions of free
fermionic or spin chains under measurements [18, 17, 19]. It would be
interesting to understand if these effective field theories are compatible
with our formalism in the single-body limit.

Finally, we note that that (4) is interesting in itself, as it seems to be a
complex version of the stochastic heat equation. Performing a Cole-Hopf [40,
41] transform on (4), i.e introducing
$h:=\frac{1}{\sqrt{\lambda}}\log\varphi$, we get, up to a constant shift in
time,

$\displaystyle dh$
$\displaystyle=i\frac{D}{\sqrt{\lambda}}\left(\nabla^{2}h+\left(\nabla
h\right)^{2}\right)+d\eta,$ (15)

which can be thought of as a complex version of the celebrated KPZ equation
[25]. To the best of our knowledge the mathematical properties of such an
equation has not been explored before. As our study shows different
qualitative and quantitative behavior than the real case, one may expect that
the complex KPZ equation entails its own rich phenomenology.

Acknowledgments Both authors thank Xhek Turkeshi and Gabriel Artur Weiderpass
for illuminating discussions. Part of this project was developed at “ _Les
Gustins”_ summer school.

Note added.—During the completion of this manuscript, the existence of a MiPT
in free fermions in 2d undergoing random projective measurements was put
forward in [42] using field theoretical methods relying on a mapping to a non-
linear sigma model.

## References

  * Skinner  _et al._ [2019] B. Skinner, J. Ruhman, and A. Nahum, Measurement-Induced Phase Transitions in the Dynamics of Entanglement, Phys. Rev. X 9, 031009 (2019).
  * Szyniszewski  _et al._ [2019] M. Szyniszewski, A. Romito, and H. Schomerus, Entanglement transition from variable-strength weak measurements, Phys. Rev. B 100, 064204 (2019).
  * Li  _et al._ [2019] Y. Li, X. Chen, and M. P. A. Fisher, Measurement-driven entanglement transition in hybrid quantum circuits, Phys. Rev. B 100, 134306 (2019).
  * Turkeshi  _et al._ [2020] X. Turkeshi, R. Fazio, and M. Dalmonte, Measurement-induced criticality in $(2+1)$-dimensional hybrid quantum circuits, Phys. Rev. B 102, 014315 (2020).
  * Gullans and Huse [2020a] M. J. Gullans and D. A. Huse, Scalable probes of measurement-induced criticality, Phys. Rev. Lett. 125, 070606 (2020a).
  * Gullans and Huse [2020b] M. J. Gullans and D. A. Huse, Dynamical purification phase transition induced by quantum measurements, Phys. Rev. X 10, 041020 (2020b).
  * Zabalo  _et al._ [2020] A. Zabalo, M. J. Gullans, J. H. Wilson, S. Gopalakrishnan, D. A. Huse, and J. H. Pixley, Critical properties of the measurement-induced transition in random quantum circuits, Phys. Rev. B 101, 060301 (2020).
  * Bao  _et al._ [2020] Y. Bao, S. Choi, and E. Altman, Theory of the phase transition in random unitary circuits with measurements, Phys. Rev. B 101, 104301 (2020).
  * Turkeshi  _et al._ [2021] X. Turkeshi, A. Biella, R. Fazio, M. Dalmonte, and M. Schiró, Measurement-induced entanglement transitions in the quantum ising chain: From infinite to zero clicks, Phys. Rev. B 103, 224210 (2021).
  * Weinstein  _et al._ [2022] Z. Weinstein, Y. Bao, and E. Altman, Measurement-induced power law negativity in an open monitored quantum circuit (2022).
  * Minoguchi  _et al._ [2022] Y. Minoguchi, P. Rabl, and M. Buchhold, Continuous gaussian measurements of the free boson CFT: A model for exactly solvable and detectable measurement-induced dynamics, SciPost Phys. 12, 009 (2022).
  * Granet  _et al._ [2023] E. Granet, C. Zhang, and H. Dreyer, Volume-law to area-law entanglement transition in a nonunitary periodic gaussian circuit, Phys. Rev. Lett. 130, 230401 (2023).
  * Lee  _et al._ [2023] G. Lee, T. Jin, Y.-X. Wang, A. McDonald, and A. Clerk, Entanglement phase transition due to reciprocity breaking without measurement or post-selection, arXiv e-prints , arXiv:2308.14614 (2023), arXiv:2308.14614 [quant-ph] .
  * Cao  _et al._ [2019] X. Cao, A. Tilloy, and A. De Luca, Entanglement in a fermion chain under continuous monitoring, SciPost Physics 7, 024 (2019).
  * Alberton  _et al._ [2021] O. Alberton, M. Buchhold, and S. Diehl, Entanglement transition in a monitored free-fermion chain: From extended criticality to area law, Phys. Rev. Lett. 126, 170602 (2021).
  * Buchhold  _et al._ [2021] M. Buchhold, Y. Minoguchi, A. Altland, and S. Diehl, Effective theory for the measurement-induced phase transition of dirac fermions, Phys. Rev. X 11, 041004 (2021).
  * Jian  _et al._ [2023] C.-M. Jian, H. Shapourian, B. Bauer, and A. W. W. Ludwig, Measurement-induced entanglement transitions in quantum circuits of non-interacting fermions: Born-rule versus forced measurements, arXiv e-prints , arXiv:2302.09094 (2023), arXiv:2302.09094 [cond-mat.stat-mech] .
  * Fava  _et al._ [2023] M. Fava, L. Piroli, T. Swann, D. Bernard, and A. Nahum, Nonlinear sigma models for monitored dynamics of free fermions, arXiv e-prints , arXiv:2302.12820 (2023), arXiv:2302.12820 [cond-mat.stat-mech] .
  * Poboiko  _et al._ [2023a] I. Poboiko, P. Pöpperl, I. V. Gornyi, and A. D. Mirlin, Theory of free fermions under random projective measurements, arXiv e-prints , arXiv:2304.03138 (2023a), arXiv:2304.03138 [quant-ph] .
  * Jin and Martin [2022] T. Jin and D. G. Martin, Kardar-parisi-zhang physics and phase transition in a classical single random walker under continuous measurement, Phys. Rev. Lett. 129, 260603 (2022).
  * Jacobs and Steck [2006] K. Jacobs and D. A. Steck, A straightforward introduction to continuous quantum measurement, Contemporary Physics 47, 279 (2006), https://doi.org/10.1080/00107510601101934 .
  * Bernard  _et al._ [2018] D. Bernard, T. Jin, and O. Shpielberg, Transport in quantum chains under strong monitoring, EPL (Europhysics Letters) 121, 60006 (2018).
  * Jin  _et al._ [2022] T. Jin, J. a. S. Ferreira, M. Filippone, and T. Giamarchi, Exact description of quantum stochastic models as quantum resistors, Phys. Rev. Res. 4, 013109 (2022).
  * Ferreira  _et al._ [2023] J. Ferreira, T. Jin, J. Mannhart, T. Giamarchi, and M. Filippone, Exact description of transport and non-reciprocity in monitored quantum devices, arXiv e-prints , arXiv:2306.16452 (2023), arXiv:2306.16452 [quant-ph] .
  * Kardar  _et al._ [1986] M. Kardar, G. Parisi, and Y.-C. Zhang, Dynamic scaling of growing interfaces, Phys. Rev. Lett. 56, 889 (1986).
  * Bertini and Cancrini [1995] L. Bertini and N. Cancrini, The stochastic heat equation: Feynman-kac formula and intermittence, Journal of Statistical Physics 78, 1377 (1995).
  * Martin  _et al._ [1973] P. C. Martin, E. D. Siggia, and H. A. Rose, Statistical dynamics of classical systems, Phys. Rev. A 8, 423 (1973).
  * Kamenev [2011] A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ (Cambridge University Press, 2011).
  * [29] Supplementary material .
  * Barabási and Stanley [1995] A.-L. Barabási and H. E. Stanley,  _Fractal Concepts in Surface Growth_ (Cambridge University Press, 1995).
  * Note [1] We are, a priori, free to fix the microscopic cut-off $\Lambda$. To do so, we extract the numerical value of the microscopic critical $\gamma_{c}$ for $\tau=1$ by simulating the full microscopic evolution (2). We found $\gamma_{c}\approx 1.14$. We then fix the value of $\Lambda$ so that the RG flow equations (11) yield the same value of $\gamma_{c}$. Doing so, we fix $\Lambda$ to $\pi/1.226595$. Because of the strong linearities of the flow equations, this value has to be finely tuned.
  * Family and Vicsek [1985] F. Family and T. Vicsek, Scaling of the active zone in the eden process on percolation networks and the ballistic deposition model, Journal of Physics A: Mathematical and General 18, L75 (1985).
  * Žnidarič [2010a] M. Žnidarič, Dephasing-induced diffusive transport in the anisotropic heisenberg model, New Journal of Physics 12, 043001 (2010a).
  * Žnidarič [2010b] M. Žnidarič, Exact solution for a diffusive nonequilibrium steady state of an open quantum chain, Journal of Statistical Mechanics: Theory and Experiment 2010, L05002 (2010b).
  * Žnidarič and Horvat [2013] M. Žnidarič and M. Horvat, Transport in a disordered tight-binding chain with dephasing, The European Physical Journal B 86, 67 (2013).
  * Bauer  _et al._ [2017] M. Bauer, D. Bernard, and T. Jin, Stochastic dissipative quantum spin chains (I) : Quantum fluctuating discrete hydrodynamics, SciPost Phys. 3, 033 (2017).
  * Bastianello  _et al._ [2020] A. Bastianello, J. De Nardis, and A. De Luca, Generalized hydrodynamics with dephasing noise, Phys. Rev. B 102, 161110 (2020).
  * Noel  _et al._ [2022] C. Noel, P. Niroula, D. Zhu, A. Risinger, L. Egan, D. Biswas, M. Cetina, A. V. Gorshkov, M. J. Gullans, D. A. Huse, and C. Monroe, Measurement-induced quantum phases realized in a trapped-ion quantum computer, Nature Physics 18, 760 (2022).
  * Koh  _et al._ [2023] J. M. Koh, S.-N. Sun, M. Motta, and A. J. Minnich, Measurement-induced entanglement phase transition on a superconducting quantum processor with mid-circuit readout, Nature Physics 10.1038/s41567-023-02076-6 (2023).
  * Cole [1951] J. D. Cole, On a quasi-linear parabolic equation occurring in aerodynamics, Quarterly of Applied Mathematics 9, 225 (1951).
  * Hopf [1950] E. Hopf, The partial differential equation $u_{t}+uu_{x}=\mu_{xx}$, Communications on Pure and Applied Mathematics 3, 201 (1950).
  * Poboiko  _et al._ [2023b] I. Poboiko, I. V. Gornyi, and A. D. Mirlin, Measurement-induced phase transition for free fermions above one dimension, arXiv e-prints , arXiv:2309.12405 (2023b), arXiv:2309.12405 [quant-ph] .
  * Mannella [1997] R. Mannella, Numerical integration of stochastic differential equations, arXiv preprint cond-mat/9709326 (1997).
  * Moss and McClintock [1989] F. Moss and P. V. E. McClintock, eds.,  _Experiments and simulations_ , Noise in nonlinear dynamical systems No. v. 3 (Cambridge University Press, Cambridge [Cambridgeshire] ; New York, 1989).

##  Appendix A SUPPLEMENTAL MATERIAL

##  Appendix B Martin-Siggia-Rose (MSR) action for complex fields

In this section we derive the MSR action of the stochastic differential
equation (SDE) :

$d\varphi=\left(iD\nabla^{2}\varphi-\frac{\gamma}{2}\varphi\right)dt+\sqrt{\lambda}\varphi
d\eta.$ (16)

Since the field $\varphi$ is complex, we need to treat individually the real
and imaginary part. Let $X:=\Re(\varphi)$ and $Y:=\Im(\varphi)$. We have the
coupled equations

$\displaystyle dX$
$\displaystyle=\left(-D\nabla^{2}Y-\frac{\gamma}{2}X\right)dt+\sqrt{\lambda}Xd\eta,$
(17) $\displaystyle dY$
$\displaystyle=\left(D\nabla^{2}X-\frac{\gamma}{2}Y\right)dt+\sqrt{\lambda}Yd\eta.$
(18)

The partition function is defined as

$Z:=\mathbb{E}\left[\int{\cal
D}[X,Y]\delta\left(\partial_{t}X+D\nabla^{2}Y+\frac{\gamma}{2}X-\sqrt{\lambda}X\eta\right)\delta\left(\partial_{t}Y-D\nabla^{2}X+\frac{\gamma}{2}Y-\sqrt{\lambda}Y\eta\right)\right],$
(19)

where $\mathbb{E}[\cdot]$ indicates an average with respect to the noise
$\eta$. We introduce the identity resolution at every point in space-time

$\delta\left(f(X)\right):=\int{\cal
D}[X^{a},Y^{a}]e^{-i\left(X^{a}+Y^{a}\right)f(X)}\;,$ (20)

where $f(X)$ is an arbitrary test function. The partition function then reads

$Z=\mathbb{E}\left[\int{\cal D}[X,X^{a},Y,Y^{a}]e^{-i\int
drdtX^{a}\left(\partial_{t}X+D\nabla^{2}Y+\frac{\gamma}{2}X-\sqrt{\lambda}X\eta\right)}e^{-i\int
drdtY^{a}\left(\partial_{t}Y-D\nabla^{2}X+\frac{\gamma}{2}Y-\sqrt{\lambda}Y\eta\right)}\right].$
(21)

We now make use of the probability distribution of the noise

$\mathbb{E}[f(\eta)]=\int{\cal D}[\eta]e^{-\int
drdt\frac{\eta^{2}}{2}}f(\eta)$ (22)

to get

$\displaystyle Z$ $\displaystyle=\int{\cal D}[X,X^{a},Y,Y^{a}]e^{iS},$ (23)

where the action $S$ is given by

$S=\int
drdt\left[-X^{a}\left(\partial_{t}X+D\nabla^{2}Y+\frac{\gamma}{2}X\right)-Y^{a}\left(\partial_{t}Y-D\nabla^{2}X+\frac{\gamma}{2}Y\right)+i\frac{\lambda}{2}\left(XX^{a}+YY^{a}\right)^{2}\right].$
(24)

Reintroducing the complex form of the fields $\varphi=X+iY$ and
$\bar{\varphi}=X-iY$, we get

$\displaystyle S$ $\displaystyle=S_{0}+S_{\nu},$ (25) $\displaystyle S_{0}$
$\displaystyle=\frac{1}{2}\int
d^{d}\vec{r}dt(\bar{\varphi},\bar{\varphi}^{a})G_{0}^{-1}\begin{pmatrix}\varphi\\\
\varphi^{a}\end{pmatrix},$ (26) $\displaystyle G_{0}^{-1}$
$\displaystyle=\begin{pmatrix}0&-\frac{\gamma}{2}+\partial_{t}^{A}-iD\nabla^{2}\\\
-\frac{\gamma}{2}-\partial_{t}^{R}+iD\nabla^{2}&0\end{pmatrix},$ (27)
$\displaystyle S_{\nu}$ $\displaystyle=\frac{i}{8}\int
d^{d}\vec{r}dt\big{(}\lambda^{{\rm
I}}\left(\bar{\varphi}^{a}\right)^{2}\varphi^{2}+2\lambda^{{\rm
II}}\bar{\varphi}^{a}\bar{\varphi}\varphi^{a}\varphi+\lambda^{{\rm
I}*}\bar{\varphi}^{2}\left(\varphi^{a}\right)^{2}\big{)}.$ (28)

In (28), we give distinct names for the different $\lambda$ coefficients as
they will follow distinct RG flows. Remark that the $4$-vertex interactions
always contain two quantum and two auxiliary fields as well as the same number
of conjugated and non-conjugated fields.

The retarded and advanced components in momentum and frequency space obtained
from inverting $G_{0}^{-1}$ are given by

$\displaystyle G_{0}^{R}(\vec{q},\omega)$ $\displaystyle:=-i\int
d^{d}\vec{r}dte^{i(\omega
t-\vec{q}.\vec{r})}\langle\varphi(r,t)\bar{\varphi}^{a}(0,0)\rangle_{S_{0}}=\frac{-2i}{\omega-
Dq^{2}+i\frac{\gamma}{2}},$ (29) $\displaystyle G_{0}^{A}(\vec{q},\omega)$
$\displaystyle:=-i\int d^{d}\vec{r}dte^{i(\omega
t-\vec{q}.\vec{r})}\langle\varphi^{a}(r,t)\bar{\varphi}(0,0)\rangle_{S_{0}}=\frac{2i}{\omega-
Dq^{2}-i\frac{\gamma}{2}},$ (30)

where $\langle\rangle_{S_{0}}$ indicates averaging over the gaussian action
$S_{0}$.

##  Appendix C One-loop perturbative renormalization

In this section we provide details on the renormalization of the mass term
$\gamma$ and the interaction term $\lambda$ of the action. We use the
momentum-shell cut-off procedure and keep the perturbation in $\lambda$ at
one-loop level. We call $\Lambda$ the momentum cut-off and $z$, $\chi$ and
$\bar{\chi}$ the dynamical exponents respectively associated to $t$,
$\varphi$, $\bar{\varphi}$. The renormalization flow is parametrized by $l$.

As explained in the main text, at the one-loop level, there are no diagrams
renormalizing the part of the action proportional to $\partial_{t}$ and
$\nabla^{2}$. Imposing the stationarity of the corresponding prefactors under
the flow gives

$\chi+\bar{\chi}+d=0,\quad z=2.$ (31)

The Feynman diagrams contributing to the renormalization of $\gamma$ and
$\lambda$ are shown graphically on Fig.5

$\lambda^{\rm{I}}$$G^{R}(q,w)$$\lambda^{\rm{II}}$$G^{A}(q,w)$MassInteractions$\lambda^{\rm{I}}$$\lambda^{\rm{I}}$$G^{R}(q,w)$$G^{R}(-q,-w)$$\lambda^{\rm{II}}$$\lambda^{\rm{II}}$$G^{R}(q,w)$$G^{A}(q,w)$
Figure 5: One-loop contributions to the RG of the mass and interactions.
Dashed lines designate the auxiliary field indexed by $a$.

We begin by providing some useful identities for the loop integrals.

###  C.1 Abacus for the one-loop integrals

As usual in the RG procedure, we separate the fields $\varphi$ and
$\bar{\varphi}$ into fast momenta $\Lambda>q>\Lambda e^{-\delta l}$ denoted
with an $f$ and slow momenta $\Lambda e^{-\delta l}>q$ denoted with an $s$. We
compute here all the fast momenta loop integrals that show up in the
renormalization of $\gamma$ and $\lambda^{{\rm I/II}}$. We begin with the sum
of the tadpoles diagrams :

$\displaystyle\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{d+1}}\left(G^{R}(\vec{q},\omega)+G^{A}(\vec{q},\omega)\right)$
$\displaystyle=-2\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{d+1}}\frac{\gamma_{R}}{\left(\omega-
Dq^{2}+i\frac{\gamma}{2}\right)\text{$\left(\omega-
Dq^{2}-i\frac{\gamma^{*}}{2}\right)$}},$ (32) $\displaystyle=-2i{\rm
sign}\left(\gamma_{R}\right)\int_{f}\frac{d^{d}\vec{q}}{(2\pi)^{d}}\left(\frac{\gamma_{R}}{i\gamma_{R}}\right),$
(33) $\displaystyle=-2{\rm sign}(\gamma_{R})K_{d}\delta l.$ (34)

Where we introduced the notations $\gamma_{R}:=\Re(\gamma)$,
$\gamma_{I}:=\Im(\gamma)$, $K_{d}:=\frac{S_{d-1}}{(2\pi)^{d}}\Lambda^{d}$ with
$S_{d-1}$ the area of a sphere of radius $1$ in $d$ dimensions, i.e.
$S_{d-1}=\frac{2\pi^{d/2}}{\Gamma(d/2)}$ with $\Gamma$ the Euler Gamma
function.

We then compute the difference between the advanced and retarded green
functions as

$\displaystyle\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{d+1}}\left(G^{R}(\vec{q},\omega)-G^{A}(\vec{q},\omega)\right)$
$\displaystyle=-4i\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{d+1}}\frac{\omega-
Dq^{2}-\frac{\gamma_{I}}{2}}{(\omega-
Dq^{2}-\frac{\gamma_{I}}{2})^{2}+\frac{\gamma_{R}{}^{2}}{4}}=0,$ (35)

where the last equality follows because the integrand is an odd function of
$w$. Combining (34) and (35) together, we will use in later computation that

$\displaystyle\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{d+1}}G^{R}(\vec{q},\omega)=\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{d+1}}G^{A}(\vec{q},\omega)$
$\displaystyle=-{\rm sign}(\gamma_{R})K_{d}\delta l.$ (36)

We now compute the two-propagators loop as

$\displaystyle\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{(d+1)}}|G^{R}(\vec{q},\omega)|^{2}$
$\displaystyle=4\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{(d+1)}}\frac{1}{\left(\omega-
Dq^{2}+i\frac{\gamma}{2}\right)\left(\omega-
Dq^{2}-i\frac{\gamma^{*}}{2}\right)},$ (37) $\displaystyle=4i{\rm
sign}(\gamma_{R})\int_{f}\frac{d^{d}\vec{q}}{(2\pi)^{d}}\frac{1}{i\gamma_{R}},$
(38) $\displaystyle=4{\rm sign}(\gamma_{R})\frac{K_{d}}{\gamma_{R}}\delta l.$
(39)

Finally, we derive the second two-propagators loop

$\displaystyle\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{(d+1)}}G^{R}(\vec{q},\omega)G^{R}(-\vec{q},-\omega)$
$\displaystyle=4\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{(d+1)}}\frac{1}{\left(\omega-
Dq^{2}+i\frac{\gamma}{2}\right)\left(\omega+Dq^{2}-i\frac{\gamma}{2}\right)},$
(40) $\displaystyle=4i{\rm
sign}(\gamma_{R})\int_{f}\frac{d^{d}\vec{q}}{(2\pi)^{d}}\frac{1}{-2Dq^{2}+i\gamma},$
(41) $\displaystyle=4{\rm sign}(\gamma_{R})K_{d}\delta
l\frac{1}{\gamma+i2D\Lambda^{2}},$ (42)

as well as its complex conjugate

$\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{(d+1)}}G^{A}(\vec{q},\omega)G^{A}(-\vec{q},-\omega)=4{\rm
sign}(\gamma_{R})K_{d}\delta l\frac{1}{\gamma-i2D\Lambda^{2}}.$ (43)

###  C.2 Flow equations

We now turn to the derivation of the flow equations. The diagrams contributing
to the renormalization of the mass and the interaction terms are shown on Fig.
5.

####  C.2.1 Mass term $\gamma$

After performing one step of renormalization, the mass term $\gamma$ is
renormalized to $\gamma^{\prime}$ according to

$\displaystyle\gamma^{\prime}$ $\displaystyle=e^{\delta
l\left(\chi+\bar{\chi}+d+z\right)}\left(\gamma+2\lambda^{{\rm
I}}\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{d+1}}G^{R}(\vec{q},\omega)+\lambda^{{\rm
II}}\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{d+1}}G^{A}(\vec{q},\omega)\right).$
(44)

For an infinitesimal step $\delta l$, (44) reads

$\frac{d\gamma}{dl}=2\gamma-{\rm sign}(\gamma_{R})K_{d}\left(2\lambda^{{\rm
I}}+\lambda^{{\rm II}}\right),$ (45)

where we made use of (31) and (36).

####  C.2.2 Interaction term

After one step of renormalization, $\lambda^{{\rm I}}$ is renormalized into
$\lambda^{{\rm I}^{\prime}}$ which is given by

$\displaystyle\lambda^{{\rm I}^{\prime}}$ $\displaystyle=e^{\delta
l\left(2\chi+2\bar{\chi}+d+z\right)}\left(\lambda^{{\rm
I}}+\frac{1}{4}\left(\lambda^{{\rm
I}}\right)^{2}\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{d+1}}G^{R}(\vec{q},\omega)G^{R}(-\vec{q},-\omega)\right),$
(46)

For an infinitesimal step $\delta l$, (46) reads

$\frac{d\lambda^{{\rm I}}}{dl}=(2-d)\lambda^{{\rm I}}+{\rm
sign}(\gamma_{R})K_{d}\frac{\left(\lambda^{{\rm
I}}\right)^{2}}{\gamma+i2D\Lambda^{2}},$ (47)

where we made use of Eqs. (31) and (42).

For $\lambda^{{\rm II}}$, one has

$\displaystyle\lambda^{{\rm II}^{\prime}}$ $\displaystyle=e^{\delta
l\left(2\chi+2\bar{\chi}+d+z\right)}\left(\lambda^{{\rm
II}}+\frac{1}{4}\left(\lambda^{{\rm
II}}\right)^{2}\int_{f}\frac{d^{d}\vec{q}d\omega}{(2\pi)^{d+1}}|G^{R}(\vec{q},\omega)|^{2}\right),$
(48)

which for an infinitesimal step becomes

$\frac{d\lambda^{{\rm II}}}{dl}=(2-d)\lambda^{{\rm II}}+{\rm
sign}(\gamma_{R})\frac{K_{d}}{\gamma_{R}}\left(\lambda^{{\rm II}}\right)^{2}.$
(49)

####  C.2.3 Summary

Regrouping all flow equations together, we end up with

$\displaystyle\frac{d\gamma}{dl}$ $\displaystyle=2\gamma-{\rm
sign}(\gamma_{R})K_{d}\left(2\lambda^{{\rm I}}+\lambda^{{\rm II}}\right),$
(50) $\displaystyle\frac{d\lambda^{{\rm I}}}{dl}$
$\displaystyle=(2-d)\lambda^{{\rm I}}+{\rm
sign}(\gamma_{R})K_{d}\frac{\left(\lambda^{{\rm
I}}\right)^{2}}{\gamma+i2D\Lambda^{2}},$ (51)
$\displaystyle\frac{d\lambda^{{\rm II}}}{dl}$
$\displaystyle=(2-d)\lambda^{{\rm II}}+{\rm
sign}(\gamma_{R})K_{d}\frac{\left(\lambda^{{\rm II}}\right)^{2}}{\gamma_{R}}.$
(52)

For convenience, we also provide the equations for the real and imaginary
parts separately :

$\displaystyle\frac{d\gamma_{R}}{dl}=$ $\displaystyle
2\gamma_{R}-K_{d}\left(\lambda^{{\rm II}}+2\lambda_{R}^{{\rm I}}\right),$ (53)
$\displaystyle\frac{d\gamma_{I}}{dl}=$ $\displaystyle
2\gamma_{I}-2K_{d}\lambda_{I}^{{\rm I}},$ (54)
$\displaystyle\frac{d\lambda_{R}^{{\rm I}}}{dl}=$
$\displaystyle(2-d)\lambda_{R}^{{\rm
I}}+K_{d}\left(\frac{\left(\left(\lambda_{R}^{{\rm
I}}\right)^{2}-\left(\lambda_{I}^{{\rm
I}}\right)^{2}\right)\gamma_{R}+2(D_{e}+\gamma_{I})\lambda_{R}^{{\rm
I}}\lambda_{I}^{{\rm I}}}{\gamma_{R}^{2}+(\gamma_{I}+D_{e})^{2}}\right),$ (55)
$\displaystyle\frac{d\lambda_{I}^{{\rm I}}}{dl}=$
$\displaystyle(2-d)\lambda_{I}^{{\rm
I}}+K_{d}\left(\frac{\left(\left(\lambda_{I}^{{\rm
I}}\right)^{2}-\left(\lambda_{R}^{{\rm
I}}\right)^{2}\right)(\gamma_{I}+D_{e})+2\lambda_{R}^{{\rm
I}}\lambda_{I}^{{\rm
I}}\gamma_{R}}{\gamma_{R}^{2}+(\gamma_{I}+D_{e})^{2}}\right),$ (56)
$\displaystyle\frac{d\lambda^{{\rm II}}}{dl}=$
$\displaystyle(2-d)\lambda^{{\rm II}}+K_{d}\frac{\left(\lambda^{{\rm
II}}\right)^{2}}{\gamma_{R}}.$ (57)

##  Appendix D Details on the numerics

We hereafter describe the methods employed to numerically integrate (2). We
first note the peculiar structure of the noise term $dB_{t}^{\boldsymbol{j}}$
in (2): it is a multiplicative multi-dimensional noise. Thus, (2) falls into
the class of stochastic differential equations taking the form

$\frac{dx_{j}}{dt}=f_{j}(\\{\textbf{x}\\})+\sum_{i=1}^{N}g_{ji}(\\{\textbf{x}\\})\xi_{i}$
(58)

where $\xi_{i}$’s are gaussian white noises such that
$\langle\xi_{i}(s)\xi_{k}(s^{\prime})\rangle=\delta_{ik}\delta(s-s^{\prime})$,
and $f_{j}$ and $g_{ji}$ are functions of the set of position
$\\{\textbf{x}\\}=\\{x_{k},\ k\in[1,..,N]\\}$. Numerical integration schemes
for SDEs of type (58) have been discussed in [43] and chapter 7 of [44]. The
combination of the multiplicative and multi-dimensional nature of the noise in
(58) renders usual higher order Runge-Kutta-based SDE algorithms inoperative.
As described in [44], the two numerical schemes available for integrating (58)
are both of order $dt$ at maximum. The first one is an Euler-Maruyama scheme
(i.e. simple forward Euler), which allows for a straightforward integration of
(58) in Itō prescription. The second one is a first order Runge-Kutta scheme
with an approximate closure valid up to $dt$: it allows for numerical
integration of (58) directly in Stratonovitch prescription. As we studied (2)
within Itō formalism in the main text, we naturally choose the former Euler-
Maruyama algorithm to perform our numerical integrations. In addition, to
ensure the correct normalization of the field $\psi_{\boldsymbol{j}}$ at each
time (i.e. that $\sum_{\boldsymbol{j}}|\psi_{\boldsymbol{j}}|^{2}=1$), we
projected as
$\psi_{\boldsymbol{j}}\to\psi_{\boldsymbol{j}}\big{/}\sqrt{\sum_{\boldsymbol{j}}|\psi_{\boldsymbol{j}}|^{2}}$
after each step forward in time. To check the convergence of the algorithm, we
divided the time step $dt$ by two and verified the stability of our results.
The width $w$ at fixed $\gamma$ and fixed system size $N$ was obtained by
averaging over at least 500 realizations of (2). We made sure that simulations
ran long enough for $w$ to effectively reach its plateau value at large time.
To compute the roughening exponent $\alpha$ at fixed $\gamma$, we performed a
linear fit of the width’s plateau value $w(t=\infty)$ as a function of the
system size $N$ in Log-Log. From (14), we indeed have that
$\log(w(t=\infty))\sim\alpha\log(N)$: the coefficient of the later linear fit
gives $\alpha$.


kroeger-vink

###### Abstract

Li-containing argyrodites represent a promising family of Li-ion conductors
with several derived compounds exhibiting room-temperature ionic conductivity
$>$ 1 mS/cm and making them attractive as potential candidates as electrolytes
in solid-state Li-ion batteries. Starting from the parent phase Li7PS6,
several cation and anion substitution strategies have been attempted to
increase the conductivity of Li ions. Nonetheless, a detailed understanding of
the thermodynamics of native defects and doping of Li argyrodite and their
effect on the ionic conductivity of Li is missing. Here, we report a
comprehensive computational study of defect chemistry of the parent phase
Li7PS6 in both intrinsic and extrinsic regimes, using a newly developed
workflow to automate the computations of several defect formation energies in
a thermodynamically consistent framework. Our findings agree with known
experimental findings, rule out several unfavorable aliovalent dopants,
narrowing down the potential promising candidates that can be tested
experimentally. We also find that cation-anion co-doping can provide a
powerful strategy to further optimize the composition of argyrodite. In
particular, Si-F co-doping is predicted to be thermodynamically favorable;
this could lead to the synthesis of the first F-doped Li-containing
argyrodite. Finally, using DeePMD neural networks, we have mapped the ionic
conductivity landscape as function of the concentration of the most promising
cation and anion dopants identified from the defect calculations, and
identified the most promising region in the compositional space with high Li
conductivity that can be explored experimentally.

##  I Introduction

A global decarbonization of the world economy requires a massive transition to
renewable energy sources, such as solar and wind as well as an electrification
of the world transportation sectors, which in turn requires reliable and
affordable technologies to store electricity. Conventional lithium-ion
batteries using liquid electrolytes are a huge technological and commercial
success; nevertheless, their performance rapidly approaching the theoretical
limit and the use of highly flammable organic solvents as liquid electrolytes
poses significant safety risks Janek . One promising technology that has been
actively researched and developed for the next-generation electricity-storage
devices is that of all-solid-state batteries, with inorganic solid ion
conductors which are intrinsically non flammable and have superior thermal
stability. These would provide improved safety and better device integration
by eliminating bulky components and the packaging required to cool the battery
and to avoid leakage of liquid electrolytes. Moreover, the use of solid-state
Li-ion conductors is believed to enable the replacement of graphite anodes by
Li-metal anodes, which will further increase the energy density of the
batteries SSB_review ; ChenC6MH00218H ; ZHENG2018198 .

Figure 1: Crystal structure of the cubic argyrodite of Li6PS5X, with the
halide X sites forming a closed packed host lattice in which PS${}^{3-}_{4}$
tetrahedra occupy the octahedral voids and half of the tetrahedral voids are
filled with the free S2-. When X = Cl or Br, there is a partial occupancy of
the free S2- and X on the 4d and 4a sites which is believed to have
significant influence on the long-range diffusion of Li ions Ohno_review

.

Several family of Li-containing compounds have been studied in the past
decades as Li-ion conductors, with several compounds exhibiting room-
temperature (RT) ionic conductivities comparable to that of liquid
electrolytes muy_review ; Ohno_review ; YU2021105858 . One prominent family of
Li-ion conductors in which several highly conducting compounds were discovered
is the Li Argyrodite family, with the general composition
Li${}^{+}_{12-2m-y-x}$M${}_{m}^{2+}$(Yy+Ch${}^{2-}_{4}$)Ch${}^{2-}_{2-x}$Xx-
(M=Ca, Mg; Y=P, As, Ge, Si, Sn, Sb; Ch=O, S, Se; X=Cl, Br, I, BH4)
YU2021105858 ; Deiseroth ; Ohno ; Kraft ; Kong_Li6PO5Cl ; Kong_As ;
Kong_Li7PS6 ; Bernges ; Sakuda ; Zhou_Sb . These can be derived from the
parent compound Li7PS6 by substitution on the P-sites or S-sites or both
Deiseroth_Li7PS6 ; Kong_Li7PS6 . The low-temperature phase of the parent
compound (space group $Pna2_{1}$) has only a moderate RT ionic conductivity,
on the order of $1.610^{-3}$ mS/cm Deiseroth_Li7PS6 . The high-temperature
(HT) cubic phase (space group $F\bar{4}3m$) exhibits, on the other hand, high
Li conductivity and can be stabilized at RT by partial substitution of S by
the halide atoms (Cl, Br and I) Deiseroth ; Kraft . The (conventional) cubic
unit cell then consists of a face-centered cubic (FCC) structure formed by the
halides (Wyckoff position 4a) with MCh${}^{3}_{4-}$tetrahedra (M on position
4b) occupying the octahedral voids as seen in figure 1. Half of the
tetrahedral voids are filled by S2- (Wyckoff 4d, so called free sulfur sites)
and the partially occupied lithium-ion positions (24g and 48h) form a cage
around the free sulfur sites (see figure 1). Alternatively, the immobile host
lattice can be described as a FCC lattice of PS4 tetrahedra with the remaining
free sulfur and halide ions occupying the octahedral and tetrahedral
interstitials. The distribution of Li ions in Li6MCh5X can be well captured
having only 48h sites occupied by 50%; depending on the unit cell volume and
material composition, occupancies on the 24g site can also be found Deiseroth
; Kraft . Within this cage-like isotropic compound, different ionic conduction
pathways exist involving three possible jumps: 1) the doublet (48h–24g–48h),
2) the intra-cage (48h–48h), 3) and the inter-cage one (also 48h–48h).

To increase the ionic conductivity of the argyrodite, several doping
strategies and synthesis routes have been developed to vary the concentration
of Li ions as well as to control the distribution of sulfur and halide ions on
the 4a and 4d sites. Despite all these efforts, a detailed understanding of
the native defect chemistry as well as doping response in Li-Argyrodites is
lacking. For instance, dopants are usually assumed to affect lithium
stoichiometry by causing the formation of charge-compensating lithium
vacancies or interstitials; for example, substitution of divalent anions S2-
by monovalent anions Cl- is assumed to result in the creation of Li+ such that
the overall system remains charge-neutral. While this assumption seems to be
verified under most synthesis conditions, there is no guarantee that it is the
only doping response especially with new synthesis routes and thermal
treatments which aim to stabilize/control a particular distribution of the
sulfur and halide anions over the octahedral and tetrahedral sites in the
argyrodite.

Moreover, simple charge compensation models usually assume that defects exist
in their formal charge state for example, -1 for Li vacancies, or +2 for S
vacancies which might be a good approximation in very ionic compounds such
oxides and halides. Sulfides are much more covalent in nature, and these
formal charges might not be a good approximation of the true defect charge
states nor be consistent with the most stable charge states as a function of
the Fermi level and the chemical potentials that reflect synthesis conditions.
In fact, an accurate treatment of defect charges has to take into account a
broad range of intrinsic and extrinsic defects in all accessible charge states
whose concentration is determined self-consistently by the constraint of phase
stability at the thermodynamic equilibrium and charge neutrality
BUCKERIDGECHEMICALPOTENTIAL ; BUCKERIDGEFERMI ; FreysoldtRMP .

##  II  Methods

###  II.1  Synthesis and Experimental Characterizations

In an Ar-filled glovebox ($<$ 1 ppm H2O, $<$1 ppm O2), the precursors Li2S
(Lorad Chemical), LiCl (Sigma-Aldrich), and P2S5 (Sigma-Aldrich) were weighed
and mixed in stoichiometric proportions to obtain 4g of the desired
compositions. The mixture was then transferred into a 45 mL ZrO2 jar filled
with 5 mm zirconia (YSZ) balls, maintaining a ball-to-powder ratio of 16.5.
The sealed jar was removed from the glovebox and placed in a Fritsch Planetary
Micro Mill Pulverisette 7, where it was ball milled at 500 RPM for 2 hours,
with 15-minute breaks every 30 minutes. The powder was recovered in the Ar-
filled glovebox and the resulting powder was placed in a closed SiC crucible,
heated to 500°C at a rate of 5°C/min in a tube furnace under N2 atmosphere,
and maintained at this temperature for 12 hours. The sample was cooled to RT,
recovered in the Ar-filled glovebox, and deagglomerated in a mortar.

X-ray diffraction of the samples was collected with an Aeris Research
diffractometer from Malvern Panalytical (Cu, 600 W, 40 kV, 15 mA, Soller slit
0.02°). Diffractograms were collected in the 10° to 90° range over 2 hours.
Lattice parameters were determined by fitting the diffraction profiles using
FullProf Suite, with profile fitting performed using the F-43M space group of
argyrodite. For impedance spectroscopy, the powders were sandwiched between
two pre-dried carbon paper electrodes (Papyex soft graphite, 0.2 mm thick from
Mersen) and cold-pressed into a 6 mm die at 500 MPa (Uniaxial manual press
MS15-MDD, Eurolabo). The obtained pellets were loaded into an airtight
measurement cell (EQ-PSC from MTI). All measurements were performed under a
pressure of 40 MPa. AC impedance spectra were collected using a Biologic VMP3,
with the sample temperature controlled by a Binder thermostatic chamber. The
cell was connected to a galvanostat-potentiostat, and PEIS spectra were
recorded with a 20 mV sinusoidal perturbation around OCV from 1 MHz to 10 kHz,
with 25 points per decade, each point averaged from 50 measures. A model
circuit was used to fit the curve and extract ionic resistance. Temperature
tests were performed with 2-hour steps at each temperature, recording a
spectrum at the end of each step. The temperature ranged from -20°C to 60°C,
then back to 30°C. Electronic conductivity was measured by DC measurement at 2
volts, determined after 1 hour by the asymptotic method.

###  II.2  Defect Calculations and Machine Learning Molecular Dynamics

In order to compute defect concentrations, one needs to calculate the
formation energy of those defects. The formalism of the defect formation
energy (in dilute limit) is well established and can be obtained using the
following relation Northrup_formalism ; FreysoldtRMP :

$\begin{split}E^{f}=E_{tot}[X_{q}]-E_{tot}[bulk]-\sum_{i}n_{i}\mu_{i}+\\\
q(E_{F}+E_{vbm}+\Delta V_{pot})+E_{corr}^{q}\end{split}$ (1)

where $E^{f}[X^{q}]$ is the formation energy of defect X in charge state q,
$E_{tot}[X^{q}]$ is the (DFT) energy of the defected supercell,
$E_{tot}[bulk]$ is the (DFT) energy of the pristing supercell without defect,
$n_{i}$ are the number of atoms of type _i_ added ($n_{i}>0$) or removed
($n_{i}<0$) from the supercell, $\mu_{i}$ are the chemical potentials of
atomic species _i_ , and $E_{F}$ is the Fermi level relative to the maximum of
the valence band $E_{vbm}$. $\Delta V_{pot}$ and $E_{corr}^{q}$ are post-
processing corrections which account respectively for the spurious interaction
between the charge defect and the uniform background charge density, and the
long-range electrostatic interactions between charge images in supercell
calculations with periodic-boundary condition freysoldt_PRL_2009 ;
freysoldt_2011 ; dabo_electrostatics_2008 . The chemical potential $\mu_{i}$
is the energy associated with exchanging atoms of species _i_ between the
materials and a reservoir. Their values reflect the conditions under which the
materials are synthesized and are constrained by the stability of the
materials and that of all the competing phases as given by the phase diagram
Northrup_formalism ; BUCKERIDGECHEMICALPOTENTIAL .

To determine the Fermi level, we employ the fact that at thermodynamic
equilibrium, in the presence of multiple defects with multiple charge states,
its value is determined by requiring charge neutrality for the system
BUCKERIDGEFERMI  More precisely, if we denote $n_{o}$ the concentration of
electrons, $p_{o}$ the concentration of holes and $C_{X^{q}}$ the
concentration of defect X in the charge state q, the charge neutrality
condition imposes that

$n_{o}-\sum_{X}\sum_{q}qC_{X^{q}}=p_{o}$ (2)

where:

$n_{o}=\int_{E_{F}}^{\infty}f_{e}(E)\rho(E)dE$ (3)
$p_{o}=\int_{-\infty}^{E_{F}}f_{h}(E)\rho(E)dE$ (4)
$C_{X^{q}}=N_{site}\exp{\left(-\frac{E^{f}[X^{q}]}{kT}\right)}$ (5)

and $f_{e}$ is the Fermi-Dirac distribution,
$f_{e}(E)=[1+\exp((E_{F}-E)/kT)]^{-1}$, $f_{h}=1-f_{e}$ and $\rho(E)$ the
density of states per unit volume. Each term $n_{o}$, $p_{o}$ and $C_{X^{q}}$
depends on $E_{F}$, and therefore the condition (2) constitutes a (non-linear)
equation in $E_{F}$ which can be solved numerically to determine the value of
Fermi level that corresponds to a neutral system. Implicit in this formalism
is the so-called rigid band approximation, where it is assumed that the only
response of the electronic structure to the presence of defects is the rigid
shift of the Fermi level to adjust to the total number of electrons while the
band structure of the materials remains unchanged.

For the present defect calculations, in the intrisic regime, we have
considered all the non-equivalent lithium vacancies and interstitials
V${}_{\text{Li}}$ and i${}_{\text{Li}}$, phoshorus vacancies V${}_{\text{P}}$,
sulfur vacancies on the octahedral (V${}_{\text{S}_{\text{oct}}}$) and
tetrahedral (V${}_{\text{S}_{\text{tet}}}$) sites as well as sulfurs that are
bonded to phosphorus atoms (V${}_{\text{S}_{1-3}}$). In the extrinsic regime,
we have considered mono- and divalent substitutions on Li sites
(X${}_{\text{P}}$, X = Na, Cu, Mg, Zn, and Ca), tri- and tetravalent
substitutions on P sites (X${}_{\text{P}}$, X = Al, Ga, In, Y, La, C, Si, Ge,
Sn, Y, La, Hf, Ce) and anion substitution on S sites (X${}_{\text{S}}$, X = F,
Cl, Br, I, O, N). Finally, we have also considered co-doping scenarios where 2
dopants are simultaneously incorporated into the parent phase Li7PS6.

To compute the ionic conductivity of lithium as a function of dopant
concentrations, we used DeePMD DeePMD_kit ; DeePMD_SE  a class of neural
network potentials which have been used successfully to simulate the dynamics
of Li-ions in the LGPS family DeePMD_LGPS  as well as other superionic
conductors such as Na3PS4Delaire_Na3PS4 , Cu7PSe6Delaire_Cu7PSe6  and
AgCrSe2Delaire_AgCrSe2 . To generate the training data, we have run first-
principle molecular dynamics (AIMD) of the conventional unitcell for 120 ps at
1000K using the PWscf code of the Quantum ESPRESSO distribution Giannozzi_2017
and the Standard Solid-State Pseudopotential (SSSP) PBEsol efficiency 1.2
Prandini ; Lejaeghere_SSSP  that verifies pseudopotentials from different
methods and libraries VANSETTEN201839 ; SCHLIPF201536 ; TOPSAKAL2014263 ;
DALCORSO2014337 ; NC_pseudo ; GARRITY2014446 . Based on the results from
defect formation energy calculations, we have selected Mg for substituion on
Li-sublattice and Cl on S-sublattice. We have included 9 compositions in our
training set (Li7-x-yMgxPS6-yCly; x=0, 0.5, 1.5 and y=0, 0.25, 0.5) and we
sampled a configuration every 6 fs. To investigate the sensitivity of the
potential with respect to the size of the training set, we have also trained
the potential using smaller datasets containing only 3 and 6 compositions. As
can be seen in Fig. S2 of the Supplementary Information (SI), the ionic
conductivity already ‘converges’ after training with only 6 compositions. We
also note that the potential trained with only 3 compositions shows some
instabilities for some compositions that are not in the training, especially
at 1250K which the highest temperature of our DeePMD simulations. The
diffusivities of Li ions were computed from the slope of the mean-square
displacements vs time and then converted into ionic conductivities using the
Nernst-Einstein relation he_statistical_2018  The initial configuration of the
argyrodite was generated in such a way that the occupancies of Cl and S at
zero Mg content on the octahedral and tetrahedral sites are as close as
possible to the experimental values obtained from neutron diffraction
gautam_Li7PS6  (see Fig. S3). Due to the lack of experimentally determined
occupancies for Cl and S in the presence of other dopants, we assume that they
do not change with the concentration of cation dopants.

Figure 2: An example of the provenance graph generated by AiiDA-defects for a
typical series of defect calculations showing the inputs (green box) and
outputs of the calculations along all the intermediate steps. All the nodes
associated with the DFT calculations are in the blue box while the purple box
contains all the nodes corresponding to the calculations of the correction
terms. The red box contain nodes coming from the calculation of the chemical
potentials and the orange box the self-consistent Fermi level.

As for the DeePMD training set, all DFT calculations for the defect formation
energies in this work were carried out using the Quantum ESPRESSO distribution
Giannozzi_2017  with PBESol PBESol_Perdew  as the exchange–correlation
functional and the pseudopotentials and energy cutoffs taken from the SSSP
PBEsol Efficiency 1.2 library Prandini ; Lejaeghere_SSSP . All calculations
were performed on the 112 atom 2x2x2 supercell constructed from the primitive
cell of the HT cubic phase of Li7PS6. The partial occupancy on the Li sites
was removed using the Ewald criterion as implemented in the Pymatgen library
ONG2013314 . All calculations for defect formation energies use a 2x2x2
Monkhorst-Pack kpoint mesh while for the AIMD, only the $\Gamma$ point was
included. To calculate the formation energies of various defects in the this
study, we employ the AiiDA-defects package Muy_2023  which leverages the
advanced capabilities of AiiDA AiiDA , an open-source materials informatics
infrastructure that provides workflow automation while simultaneously
preserving and storing the full data provenance in a relational database that
is queryable, thus ensuring the reproducibility of the calculations. In AiiDA-
defects, each term in Eq.1 is computed by a dedicated AiiDA workchain which
automatically generates inputs and gathers results of the calculations and
stores them in a node which can be later readily queried for further analysis.
An example of the provenance graph generated by AiiDA-defects for a small set
of defect formation energy calculations is shown in Fig. 2 and clearly
illustrates how much automation is achieved.

Figure 3: Concentrations of native defects in Li7PS6. The dominant defects in
thermodynamic equilibrium are Li vacancies and interstitials. There is also a
comparable concentration of electrons, suggesting that stoichiometric Li7PS6
exhibits a _n_ -type electronic conductivity in addition to Li-ion
conductivity. Figure 4: Concentrations of isovalent and aliovalent
substitutions on the Li, P and S sub-lattices in the Li7PS6 argyrodite. Figure
5: a) Concentration of anion substitutions on S sites as function of cation
co-dopants. Generally speaking, the presence of cations as dopants doesn’t
affect the concentration of anion dopants in any significant way, except the
case of Si co-dopant where the concentration of the anion can be increased
substantially. b) Concentration of anion substitutions on S sites as a
function of other anion co-dopants. In general, it can be seen that the
presence of other anion co-dopants, especially within the halogen family, can
significantly alter the concentration of each anion in Li7PS6 argyrodites.

##  III Results and discussions

###  III.1 Native defects

The concentrations of native defect in the HT phase of Li7PS6 at 300K are
shown in Figure 3. The dominant ionic defects are predicted to be lithium
vacancies and interstitials. There is also comparable concentration of
electrons, suggesting that stoichiometric Li7PS6 exhibits a _n_ -type
electronic conductivity in addition to Li-ion conductivity. Notice that this
is the concentration of defects corresponding to the ‘typical’ chemical
potential of each element as given by the centroid of the stability region
BUCKERIDGECHEMICALPOTENTIAL . These concentrations are expected to change when
the chemical potentials are chosen from a different point in the stability
region that corresponds for instance to more S-rich or S-poor conditions.

###  III.2 Isovalent and aliovalent substitutions

In the following section, we discuss the doping response of Li7PS6 to
isovalent and aliovalent substitutions on Li, P and S sites.

####  III.2.1 Anion substitutions on S sites

The computed concentration of halide substitutions on the S sites are shown in
Fig. 4 (more detailed breakdown of the concentrations of halide substitutions
on the each non-equivalent S sites are shown in Figure S1). It can be seen
clearly from Figure 4 that Cl, Br and I substitutions are much more favorable
than F substitution in agreement with the apparent absence of F-substituted
Li-argyrodites in the litterature Ohno_review ; muy_review ; YU2021105858 .
Moreover, the most favorable sulfur sites for halide substitutions are the
octahedral (4a) and tetrahedral (4d) sites, as opposed to the 16d sites in
which the S ions are bonded to the phosphorus ions (Fig. S1), also in
agreement with experiments. Finally, we note that our calculations suggest
that at room-temperature the concentrations of halide substitutions on the
octahedral sites are almost 5 orders of magnitude larger than that on the
tetrahedral sites which implies that at RT, the halides (Cl-, Br- and I-)
anions will only occupy the octahedral site. Experimentally, this has been
shown to be indeed the case for I substitutions but not for Cl and Br
substitutions where it is well established that they can both occupy tet and
oct sites along with S ions leading to a disorder on the sublattice which can
greatly enhance the long-range diffusion of Li-ions deKlerk ;
kraft_influence_2017 ; Kraft ; Morgan ; Stamminger ; RayavarapuMAY . This
apparent discrepancy might point to the fact that the actual partial occupancy
of Cl and Br on oct/tet site at RT is actually a metastable structure which is
kinetically trapped from the higher temperatures at which these compounds are
synthesized. This hypothesis is consistent with the fact that thermal
treatments have been used to optimize the occupancy of S and Cl/Br on the
oct/tet sites to achieve higher ionic conductivity Gautam .

####  III.2.2 Cation substitutions on Li and P sites

To study the effect of doping response to cation substitutions, we have
considered some of the most chemically plausible mono-, di-, tri- and tetra-
valent cations on the Li and P sites. The concentrations of these
substitutions as determined from the defect formation energies are displayed
in Figure. 4. For the isovalent substitutions, the calculations show that
Na${}_{\text{Li}}$ and Cu${}_{\text{Li}}$ are energetically very favorable,
which is consistent with the presence of Ag argyrodite analogue
Eulenberger1977 ; Hahn1965 . Divalent dopants such as Mg2+ and Zn2+ can also
readily replace Li+ ions, suggesting another less explored strategy to fine-
tune the concentration of Li and therefore the its ionic conductivity. A major
downside of this approach is that at high concentrations these ions might co-
diffuse with Li, which can be detrimental for battery operation. The critical
concentrations of dopants at which this co-diffusion becomes non-negligible
requires further investigations. Ca2+ is predicted to be less favorable than
Mg2+ and Zn2+ but could still be doped at practically relevant concentration.
This is indeed in agreement with the reported synthesis of
Li5.35Ca0.1PS4.5Cl1.55 with ionic conductivities of 10.2 mS/cm at RT
Adeli_Ca_doped . Al3+ substitution on Li sites is predicted to be very
unfavorable, which seems to be in contradiction with the report synthesis of
Li6.15Al0.15Si1.55S6 HUANG2019487 , suggesting again that this composition
corresponds to a metastable structure rather than a thermodynamically stable
one. Finally, among all the tri- and tetravalent substitutions on P-sites
considered in this study, only Si and Ge show favorable energetics at
thermodynamic equilibrium, consistent with the reports of successful synthesis
of Li7+xP1-xMxS6 with M = Si and Ge; $x\approx 0.3-0.35$. SCHNEIDER2017151 ;
C8TA10790D

####  III.2.3 Cation-anion co-dopings

Several compounds in the argyrodite family that show high ionic conductivity
at RT are actually obtained from the parent compound Li7PS6 by co-substitution
of both cation and anion as in, for example, Li6+xP1-xMxS5I (M=Si, Ge) Kraft .
It is therefore of great interest to study these co-doping scenarios in detail
in order to find novel pairs of co-dopants which are thermodynamically
favorable. The formalism used in previous sections to study a single dopant
can be readily extended to include two or even more dopants. The main
assumption here is that the defects do not interact directly with one another
and that the presence of other defects is mediated only through the charge-
neutrality condition imposed in the calculation of the Fermi level. With this
assumption, one doesn’t have to compute the energy of the supercell with two
dopants simultaneously, only the (DFT) energy of supercells with single
dopants is required and the effect of co-dopants is entirely captured in the
change of the self-consistent Fermi level as determined by Eq. 2. The
concentrations of halide, oxygen and nitrogen ions as a function of cation co-
dopants are shown in Figure 5a. One can see that in general the presence of
cation co-dopant doesn’t affect the concentration of anion dopants in any
significant way, except the case of Si co-dopant where the concentration of
the anion can be increased substantially, raising the prospect that one might
be able to synthesize the first F-doped argyrodite by co-doping the compound
with Si. The results for the anion-anion co-doping are shown in Figure 5b.
Unlike cation co-doping, the presence of other anion co-dopants, especially
within the halogen family, can significantly alter the concentration of the
other anion in Li7PS6 argyrodites.

Figure 6: The ionic conductivity and activation energy landscape of Mg-Cl and
Si-Cl co-doped Li7PS6 argyrodite as computed using the DeePMD potential
DeePMD_kit ; DeePMD_SE . The dashed lines corresponds to the compositions with
the same Li content. The region of the compositional space with the highest
ionic conductivity (lowest activation energy) are the ones with high Cl
content and low Mg/Si contents.

###  III.3 Ionic conductivity landscape

Based on the calculations of defect formation energies, Mg, Si and Cl are
predicted to be some of the most promising dopants on Li, P and S sub-lattices
respectively. In order to investigate the effects of these dopants on the
ionic conductivity of Li-ion in the (cubic) argyrodite structure, we have
trained several DeePMD potentials DeePMD_kit ; DeePMD_SE  using the protocol
previously described in Section II to run classical MD of the argyrodite over
a wide range of dopant concentrations. The ionic conductivities at 300K and
the activation energies for Mg-Cl and Si-Cl co-doping are shown in Fig. 6. The
ionic conductivities are extrapolated from higher temperatures (1250 K, 1000
K, 840 K, 715 K and 625 K) to room temperature, using Arrhenius equation. We
have verified for selected compositions that the extrapolated and the computed
values are in very good agreement (see Fig. S4). The lower bound of Cl is set
to 0.25 and 0.5 respectively for Mg and Si co-doping scenario, as it is known
that the cubic phase of pure Li7PS6 is not stable at room temperature. The
upper bound of Cl is set to 1.5 as it is was found that secondary phases start
to form beyond this limit gautam_Li7PS6 . While there are several local
maxima/minima at some particular compositions, in both cases it is clear that
the regions of the compositional space with the highest ionic conductivity
(lowest activation energy) are those with high Cl content and low Mg/Si
content. Based on this information, we have synthesized several (Mg,Cl)-doped
argyrodites Li6-2xMgxPS6-yCly with x=0, 0.125, 0.0625, 0.1 and 0.25; y = 1,
1.375 and 1.5. All compositions were found to be phase-pure except for
Li5.25Mg0.25PS4.75Cl1.25 where trace of MgS was detected. The measured ionic
conductivity of these samples between -20 and 60oC are shown in Fig. 7 and
were in good agreement with the computed results. The activation energies of
Li conduction were found to be 0.38 eV for Li6PS5Cl, 0.38 eV for
Li5.375Mg0.125PS4.625Cl1.375, 0.36 eV for Li5.375Mg0.0625PS4.5Cl1.5 and 0.38
eV for Li5.8Mg0.1PS5Cl. RT Electronic conductivity for all samples is smaller
than $10^{-8}$ S/cm patent_1 ; patent_2 . These experimental measurements
demonstrate the effectiveness of computational materials design as a
complementary tool to accelerate the discovery of new functional materials
with improved properties.

Figure 7: Experimental ionic conductivity of Li in Mg-doped Argyrodite
Li${}_{7}-2x-y$MgxPS${}_{6}-y$Cly measured between -20 and 60oC, showing good
agreement with the computed results.

##  IV Conclusion

In this study, we have provided a thorough description of the intrinsic
defects as well as dopant regimes in the Li argyrodite Li7PS6 using the AiiDA-
defects package to enable the calculation of defect formation energies in a
robust, reproducible and thermodynamically consistent way. We found that
halide dopants, with the exception of iodine readily substitute on the
octahedral and tetrahedral sites, in agreement with the experiments. We also
find that Mg, Zn, Si and Ge are some of the most promising dopants to vary the
concentration of Li-ion and its ionic conductivity, while tri-valent dopants
are predicted to be thermodynamically unfavorable. We have also trained DeePMD
potentials to perform MD simulations to investigate the effects of these
dopants on the ionic conductivities of Li, by establishing an ionic
conductivity map that indicates that the highest Li-ion conductivity can be
achieved by having high concentrations of Cl but low concentration of Mg/Si
dopants. This approach provides valuable information to guide the synthesis
and optimization of new argyrodite compositions with higher ionic
conductivity.

## Acknowledgement

The authors would like to thank the Swiss National Supercomputing Centre CSCS
(project s1073) for providing the computational ressources and Syensqo for
their financial and scientific supports for this project.

## References

  * (1) Adeli, P., Bazak, J. D., Huq, A., Goward, G. R., and Nazar, L. F. Influence of aliovalent cation substitution and mechanical compression on li-ion conductivity and diffusivity in argyrodite solid electrolytes.  Chemistry of Materials 33, 1 (2021), 146–157.
  * (2) Bachman, J. C., Muy, S., Grimaud, A., Chang, H.-H., Pour, N., Lux, S. F., Paschos, O., Maglia, F., Lupart, S., Lamp, P., Giordano, L., and Shao-Horn, Y. Inorganic solid-state electrolytes for lithium batteries: Mechanisms and properties governing ion conduction.  Chemical Reviews 116, 1 (2016), 140–162.
  * (3) Bernges, T., Culver, S. P., Minafra, N., Koerver, R., and Zeier, W. G. Competing structural influences in the li superionic conducting argyrodites li6ps5–xsexbr (0 $\leq$ x $\leq$ 1) upon se substitution.  Inorganic Chemistry 57, 21 (2018), 13920–13928.
  * (4) Buckeridge, J. Equilibrium point defect and charge carrier concentrations in a material determined through calculation of the self-consistent fermi energy.  Computer Physics Communications 244 (2019), 329 – 342.
  * (5) Buckeridge, J., Scanlon, D., Walsh, A., and Catlow, C. Automated procedure to determine the thermodynamic stability of a material and the range of chemical potentials necessary for its formation relative to competing phases and compounds.  Computer Physics Communications 185, 1 (2014), 330 – 338.
  * (6) Chen, R., Qu, W., Guo, X., Li, L., and Wu, F. The pursuit of solid-state electrolytes for lithium batteries: from comprehensive insight to emerging horizons.  Mater. Horiz. 3 (2016), 487–516.
  * (7) Dabo, I., Kozinsky, B., Singh-Miller, N. E., and Marzari, N. Electrostatics in periodic boundary conditions and real-space corrections.  Physical Review B 77, 11 (Mar. 2008), 115139.
  * (8) Dal Corso, A. Pseudopotentials periodic table: From h to pu.  Computational Materials Science 95 (2014), 337–350.
  * (9) de Klerk, N. J. J., Rosłoń, I., and Wagemaker, M. Diffusion mechanism of li argyrodite solid electrolytes for li-ion batteries and prediction of optimized halogen doping: The effect of li vacancies, halogens, and halogen disorder.  Chemistry of Materials 28, 21 (2016), 7955–7963.
  * (10) Deiseroth, H.-J., Kong, S.-T., Eckert, H., Vannahme, J., Reiner, C., Zaiß, T., and Schlosser, M. Li6ps5x: A class of crystalline li-rich solids with an unusually high li+ mobility.  Angewandte Chemie International Edition 47, 4 (2008), 755–758.
  * (11) Deiseroth, H.-J., Maier, J., Weichert, K., Nickel, V., Kong, S.-T., and Reiner, C. Li7ps6 and li6ps5x (x: Cl, br, i): Possible three-dimensional diffusion pathways for lithium ions and temperature dependence of the ionic conductivity by impedance measurements.  Zeitschrift für anorganische und allgemeine Chemie 637, 10 (2011), 1287–1294.
  * (12) Ding, J., Niedziela, J. L., Bansal, D., Wang, J., He, X., May, A. F., Ehlers, G., Abernathy, D. L., Said, A., Alatas, A., Ren, Y., Arya, G., and Delaire, O. Anharmonic lattice dynamics and superionic transition in agcrse2.  Proceedings of the National Academy of Sciences 117, 8 (2020), 3930–3937.
  * (13) EMERY, A., BRAIDA, M.-D., Mercier, T. L., Marzari, N., MUY, S., and ERCOLE, L. Solid material comprising li, mg, p, s and halogen elements, jun 2023.
  * (14) EMERY, A., BRAIDA, M.-D., Mercier, T. L., Marzari, N., MUY, S., and ERCOLE, L. Solid material comprising li, mg, p, s and halogen elements, jun 2023.
  * (15) Eulenberger, G. Die kristallstruktur der tieftemperaturmodifikation von ag8ges6.  Monatshefte für Chemie / Chemical Monthly 108, 4 (Jul 1977), 901–913.
  * (16) Freysoldt, C., Grabowski, B., Hickel, T., Neugebauer, J., Kresse, G., Janotti, A., and Van de Walle, C. G. First-principles calculations for point defects in solids.  Rev. Mod. Phys. 86 (Mar 2014), 253–305.
  * (17) Freysoldt, C., Neugebauer, J., and Van de Walle, C. G. Fully Ab Initio Finite-Size Corrections for Charged-Defect Supercell Calculations.  Physical Review Letters 102, 1 (Jan. 2009), 016402.
  * (18) Freysoldt, C., Neugebauer, J., and Walle, C. G. V. d. Electrostatic interactions between charged defects in supercells.  physica status solidi (b) 248, 5 (2011), 1067–1076.
  * (19) Garrity, K. F., Bennett, J. W., Rabe, K. M., and Vanderbilt, D. Pseudopotentials for high-throughput dft calculations.  Computational Materials Science 81 (2014), 446–452.
  * (20) Gautam, A., Ghidiu, M., Suard, E., Kraft, M. A., and Zeier, W. G. On the lithium distribution in halide superionic argyrodites by halide incorporation in li7–xps6–xclx.  ACS Applied Energy Materials 4, 7 (2021), 7309–7315.  Publisher: American Chemical Society.
  * (21) Gautam, A., Sadowski, M., Prinz, N., Eickhoff, H., Minafra, N., Ghidiu, M., Culver, S. P., Albe, K., Fässler, T. F., Zobel, M., and Zeier, W. G. Rapid crystallization and kinetic freezing of site-disorder in the lithium superionic argyrodite li6ps5br.  Chemistry of Materials 31, 24 (2019), 10178–10185.
  * (22) Giannozzi, P., Andreussi, O., Brumme, T., Bunau, O., Nardelli, M. B., Calandra, M., Car, R., Cavazzoni, C., Ceresoli, D., Cococcioni, M., Colonna, N., Carnimeo, I., Corso, A. D., de Gironcoli, S., Delugas, P., DiStasio, R. A., Ferretti, A., Floris, A., Fratesi, G., Fugallo, G., Gebauer, R., Gerstmann, U., Giustino, F., Gorni, T., Jia, J., Kawamura, M., Ko, H.-Y., Kokalj, A., Küçükbenli, E., Lazzeri, M., Marsili, M., Marzari, N., Mauri, F., Nguyen, N. L., Nguyen, H.-V., de-la Roza, A. O., Paulatto, L., Poncé, S., Rocca, D., Sabatini, R., Santra, B., Schlipf, M., Seitsonen, A. P., Smogunov, A., Timrov, I., Thonhauser, T., Umari, P., Vast, N., Wu, X., and Baroni, S. Advanced capabilities for materials modelling with quantum ESPRESSO.  Journal of Physics: Condensed Matter 29, 46 (oct 2017), 465901.
  * (23) Gupta, M. K., Ding, J., Bansal, D., Abernathy, D. L., Ehlers, G., Osti, N. C., Zeier, W. G., and Delaire, O. Strongly anharmonic phonons and their role in superionic diffusion and ultralow thermal conductivity of cu7pse6.  Advanced Energy Materials 12, 23 (2022), 2200596.
  * (24) Gupta, M. K., Ding, J., Osti, N. C., Abernathy, D. L., Arnold, W., Wang, H., Hood, Z., and Delaire, O. Fast na diffusion and anharmonic phonon dynamics in superionic na3ps4.  Energy Environ. Sci. 14 (2021), 6554–6563.
  * (25) Hahn, H., Schulze, H., and Sechser, L. Über einige ternäre chalkogenide vom argyrodit-typ.  Naturwissenschaften 52, 15 (Jan 1965), 451–451.
  * (26) He, X., Zhu, Y., Epstein, A., and Mo, Y. Statistical variances of diffusional properties from ab initio molecular dynamics simulations.  npj Computational Materials 4, 1 (2018), 1–9.
  * (27) Huang, J., Zhang, L., Wang, H., Zhao, J., Cheng, J., and E, W. Deep potential generation scheme and simulation protocol for the li10gep2s12-type superionic conductors.  The Journal of Chemical Physics 154, 9 (2021), 094703.
  * (28) Huang, W., Yoshino, K., Hori, S., Suzuki, K., Yonemura, M., Hirayama, M., and Kanno, R. Superionic lithium conductor with a cubic argyrodite-type structure in the li–al–si–s system.  Journal of Solid State Chemistry 270 (2019), 487–492.
  * (29) Huber, S. P., Zoupanos, S., Uhrin, M., Talirz, L., Kahle, L., Häuselmann, R., Gresch, D., Müller, T., Yakutovich, A. V., Andersen, C. W., Ramirez, F. F., Adorf, C. S., Gargiulo, F., Kumbhar, S., Passaro, E., Johnston, C., Merkys, A., Cepellotti, A., Mounet, N., Marzari, N., Kozinsky, B., and Pizzi, G. AiiDA 1.0, a scalable computational infrastructure for automated reproducible workflows and data provenance.  Scientific Data 7, 1 (Sept. 2020), 300.
  * (30) Janek, J., and Zeier, W. A solid future for battery development.  Nature Energy 1, 9 (2016), 16141 – 16144.
  * (31) Kong, S.-T., Deiseroth, H.-J., Maier, J., Nickel, V., Weichert, K., and Reiner, C. Li6po5br and li6po5cl: The first lithium-oxide-argyrodites .  Zeitschrift für anorganische und allgemeine Chemie 636, 11 (2010), 1920–1924.
  * (32) Kong, S.-T., Deiseroth, H.-J., Reiner, C., Gün, O., Neumann, E., Ritter, C., and Zahn, D. Lithium argyrodites with phosphorus and arsenic: Order and disorder of lithium atoms, crystal chemistry, and phase transitions.  Chemistry – A European Journal 16, 7 (2010), 2198–2206.
  * (33) Kong, S.-T., Gün, O., Koch, B., Deiseroth, H.-J., Eckert, H., and Reiner, C. Structural characterisation of the li argyrodites li7ps6 and li7pse6 and their solid solutions: Quantification of site preferences by mas-nmr spectroscopy.  Chemistry – A European Journal 16, 17 (2010), 5138–5147.
  * (34) Kraft, M. A., Culver, S. P., Calderon, M., Böcher, F., Krauskopf, T., Senyshyn, A., Dietrich, C., Zevalkink, A., Janek, J., and Zeier, W. G. Influence of Lattice Polarizability on the Ionic Conductivity in the Lithium Superionic Argyrodites Li${}_{\textrm{6}}$PS${}_{\textrm{5}}$X (X = Cl, Br, I).  Journal of the American Chemical Society 139, 31 (Aug. 2017), 10909–10918.
  * (35) Kraft, M. A., Ohno, S., Zinkevich, T., Koerver, R., Culver, S. P., Fuchs, T., Senyshyn, A., Indris, S., Morgan, B. J., and Zeier, W. G. Inducing high ionic conductivity in the lithium superionic argyrodites li6+xp1–xgexs5i for all-solid-state batteries.  Journal of the American Chemical Society 140, 47 (2018), 16330–16339.
  * (36) Lejaeghere, K., Bihlmayer, G., Björkman, T., Blaha, P., Blügel, S., Blum, V., Caliste, D., Castelli, I. E., Clark, S. J., Corso, A. D., de Gironcoli, S., Deutsch, T., Dewhurst, J. K., Marco, I. D., Draxl, C., Dułak, M., Eriksson, O., Flores-Livas, J. A., Garrity, K. F., Genovese, L., Giannozzi, P., Giantomassi, M., Goedecker, S., Gonze, X., Grånäs, O., Gross, E. K. U., Gulans, A., Gygi, F., Hamann, D. R., Hasnip, P. J., Holzwarth, N. A. W., Iuşan, D., Jochym, D. B., Jollet, F., Jones, D., Kresse, G., Koepernik, K., Küçükbenli, E., Kvashnin, Y. O., Locht, I. L. M., Lubeck, S., Marsman, M., Marzari, N., Nitzsche, U., Nordström, L., Ozaki, T., Paulatto, L., Pickard, C. J., Poelmans, W., Probert, M. I. J., Refson, K., Richter, M., Rignanese, G.-M., Saha, S., Scheffler, M., Schlipf, M., Schwarz, K., Sharma, S., Tavazza, F., Thunström, P., Tkatchenko, A., Torrent, M., Vanderbilt, D., van Setten, M. J., Speybroeck, V. V., Wills, J. M., Yates, J. R., Zhang, G.-X., and Cottenier, S. Reproducibility in density functional theory calculations of solids.  Science 351, 6280 (2016), aad3000.
  * (37) Morgan, B. J. Mechanistic origin of superionic lithium diffusion in anion-disordered li6ps5x argyrodites.  Chemistry of Materials 33, 6 (2021), 2004–2018.
  * (38) Muy, S., Johnston, C., and Marzari, N. Aiida-defects: an automated and fully reproducible workflow for the complete characterization of defect chemistry in functional materials.  Electronic Structure 5, 2 (jun 2023), 024009.
  * (39) Ohno, S., Banik, A., Dewald, G. F., Kraft, M. A., Krauskopf, T., Minafra, N., Till, P., Weiss, M., and Zeier, W. G. Materials design of ionic conductors for solid state batteries.  Progress in Energy 2, 2 (mar 2020), 022001.
  * (40) Ohno, S., Helm, B., Fuchs, T., Dewald, G., Kraft, M. A., Culver, S. P., Senyshyn, A., and Zeier, W. G. Further evidence for energy landscape flattening in the superionic argyrodites li6+xp1–xmxs5i (m = si, ge, sn).  Chemistry of Materials 31, 13 (2019), 4936–4944.
  * (41) Ong, S. P., Richards, W. D., Jain, A., Hautier, G., Kocher, M., Cholia, S., Gunter, D., Chevrier, V. L., Persson, K. A., and Ceder, G. Python materials genomics (pymatgen): A robust, open-source python library for materials analysis.  Computational Materials Science 68 (2013), 314–319.
  * (42) Perdew, J. P., Ruzsinszky, A., Csonka, G. I., Vydrov, O. A., Scuseria, G. E., Constantin, L. A., Zhou, X., and Burke, K. Erratum: Restoring the density-gradient expansion for exchange in solids and surfaces [phys. rev. lett. 100, 136406 (2008)].  Phys. Rev. Lett. 102 (Jan 2009), 039902.
  * (43) Prandini, G., Marrazzo, A., Castelli, I. E., Mounet, N., and Marzari, N. Precision and efficiency in solid-state pseudopotential calculations.  NPJ COMPUTATIONAL MATERIALS 4 (DEC 6 2018).
  * (44) Rayavarapu, P. R., Sharma, N., Peterson, V. K., and Adams, S. Variation in structure and li+-ion migration in argyrodite-type li6ps5x (x = cl, br, i) solid electrolytes.  JOURNAL OF SOLID STATE ELECTROCHEMISTRY 16, 5 (MAY), 1807–1813.
  * (45) Sakuda, A., Yamauchi, A., Yubuchi, S., Kitamura, N., Idemoto, Y., Hayashi, A., and Tatsumisago, M. Mechanochemically prepared li2s–p2s5–libh4 solid electrolytes with an argyrodite structure.  ACS Omega 3, 5 (2018), 5453–5458.
  * (46) Schlipf, M., and Gygi, F. Optimization algorithm for the generation of oncv pseudopotentials.  Computer Physics Communications 196 (2015), 36–44.
  * (47) Schneider, H., Du, H., Kelley, T., Leitner, K., ter Maat, J., Scordilis-Kelley, C., Sanchez-Carrera, R., Kovalev, I., Mudalige, A., Kulisch, J., Safont-Sempere, M. M., Hartmann, P., Weiß, T., Schneider, L., and Hinrichsen, B. A novel class of halogen-free, super-conductive lithium argyrodites: Synthesis and characterization.  Journal of Power Sources 366 (2017), 151–160.
  * (48) Stamminger, A. R., Ziebarth, B., Mrovec, M., Hammerschmidt, T., and Drautz, R. Ionic conductivity and its dependence on structural disorder in halogenated argyrodites li6ps5x (x = br, cl, i).  Chemistry of Materials 31, 21 (2019), 8673–8678.
  * (49) Topsakal, M., and Wentzcovitch, R. Accurate projected augmented wave (paw) datasets for rare-earth elements (re=la–lu).  Computational Materials Science 95 (2014), 263–270.
  * (50) van Setten, M., Giantomassi, M., Bousquet, E., Verstraete, M., Hamann, D., Gonze, X., and Rignanese, G.-M. The pseudodojo: Training and grading a 85 element optimized norm-conserving pseudopotential table.  Computer Physics Communications 226 (2018), 39–54.
  * (51) Wang, H., Zhang, L., Han, J., and E, W. Deepmd-kit: A deep learning package for many-body potential energy representation and molecular dynamics.  Computer Physics Communications 228 (2018), 178–184.
  * (52) Willand, A., Kvashnin, Y. O., Genovese, L., Vázquez-Mayagoitia, A., Deb, A. K., Sadeghi, A., Deutsch, T., and Goedecker, S. Norm-conserving pseudopotentials with chemical accuracy compared to all-electron calculations.  The Journal of Chemical Physics 138, 10 (2013), 104109.
  * (53) Yu, C., Zhao, F., Luo, J., Zhang, L., and Sun, X. Recent development of lithium argyrodite solid-state electrolytes for solid-state batteries: Synthesis, structure, stability and dynamics.  Nano Energy 83 (2021), 105858.
  * (54) Zhang, L., Han, J., Wang, H., Saidi, W., Car, R., and E, W. End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems.  In Advances in Neural Information Processing Systems (2018), S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, Eds., vol. 31, Curran Associates, Inc.
  * (55) Zhang, S. B., and Northrup, J. E. Chemical potential dependence of defect formation energies in gaas: Application to ga self-diffusion.  Phys. Rev. Lett. 67 (Oct 1991), 2339–2342.
  * (56) Zhang, Z., Shao, Y., Lotsch, B., Hu, Y.-S., Li, H., Janek, J., Nazar, L. F., Nan, C.-W., Maier, J., Armand, M., and Chen, L. New horizons for inorganic solid state ion conductors.  Energy Environ. Sci. 11 (2018), 1945–1976.
  * (57) Zhang, Z., Sun, Y., Duan, X., Peng, L., Jia, H., Zhang, Y., Shan, B., and Xie, J. Design and synthesis of room temperature stable li-argyrodite superionic conductors via cation doping.  J. Mater. Chem. A 7 (2019), 2717–2722.
  * (58) Zheng, F., Kotobuki, M., Song, S., Lai, M. O., and Lu, L. Review on solid electrolytes for all-solid-state lithium-ion batteries.  Journal of Power Sources 389 (2018), 198–213.
  * (59) Zhou, L., Assoud, A., Zhang, Q., Wu, X., and Nazar, L. F. New family of argyrodite thioantimonate lithium superionic conductors.  Journal of the American Chemical Society 141, 48 (2019), 19002–19013.

Supplementary materials

Optimizing ionic conductivity of Lithium in Li7PS6 argyrodite via dopant
engineering

Figure S1: Concentration of halide substitutions on various sulfur sites in
cubic argyrodite Li7PS6. The concentration of F at equilibrium is very low
compared to other halide substitutions, consistent with the absence of reports
of F-doped Li argyrodites. On the other hand, Cl, Br and I readily substitute
S on the octahedral and tetahedral sites, in agreement with experimental
findings.

Figure S2: Comparison between the measured values of Cl occupancy on the
octahedral and tetrahedral sites (green) gautam_Li7PS6 , and those used in the
structural models in the MD simulations (red)

.

Figure S3: Ionic conductivity landscape obtained from 3 DeePMD potentials
which were trained on 3, 6 and 9 compositions respectively (shown as the blue
dots in the upper panels).

.

Figure S4: Comparison between the computed and extrapolated diffusivity of Li
in a) Li5.75Mg0.375PS4.5Cl0.5 and b) Li5.375Mg0.25PS4.875Cl1.125.