# Position-Dependent Diffusion induced Non-monotonic decay of Certain Non-
Equilibrium Phenomena in Condensed Phase

Sagnik Ghosh  Indian Institute of Science Education and Research, Pune-411008,
India<EMAIL_ADDRESS>Alok Samanta & Swapan K Ghosh
UM-DAE Centre for Excellence in Basic Sciences, University of Mumbai, Kalina,
Santacruz (E), Mumbai-400098, India<EMAIL_ADDRESS>
(December 12, 2021)

###### Abstract

The dynamics of various optically controlled non-equilibrium phenomena in the
condensed phase are studied using the Liouville equation. We study a
projection of the same in a slow moving coordinate, identified as the Reaction
Coordinate approach, with a position dependent diffusion coefficient.
Introduction of position dependence is shown to induce non-monotonicity in
relaxations of certain Non-equilibrium correlation functions, previously
unexplored in the theoretical as well as experimental studies. This is in
contrast to the exponential relaxation of its position independent analogue,
irrespective of initial conditions. We characterize the dependence of this
non-monotonicity on the strength of spatial inhomogeneity of diffusion and on
the strength of the restoring forces and also indicate ranges of combinations
where this feature is exhibited to pave the way for its experimental
detection.

##  I Introduction

The study of non-equilibrium processes in presence of condensed phases as a
medium forms a frontier facet of condensed matter theory. The advancement of
ultra-fast spectroscopy, such as the pump-probe method, where one excites a
subsystem using a short pulse femto second laser, and observe the systems
transient dynamics and subsequent relaxation to equilibrium, usually in terms
of dynamics of its correlation functions.

These experiments have led to explorations of optical control of various non-
equilibrium processes in the condensed phase, such as solvation, symmetric and
assymetric electron-transfers, two-step kinetics etc, especially in the
context of their dependence on the excitation frequency. The early theoretical
explorations in this field is due to Onsanger Onsager (1931a, b). The
foundation of these treatments under Linear Response Theory (LRT) is based on
the assumption that the nature of the fluctuations governing the non-
equilibrium dynamics is essentially the same as of its equilibrium
counterparts. Thus in its scope the dynamics remains essentially Markovian and
fails to demarcate its dependence on the initial preparations, that are
reported in various pump-probe type experiments in the context of solvation
dynamics Sen  _et al._ (2006); Mondal  _et al._ (2006); Mandal  _et al._
(2008); Ghatak  _et al._ (2011); Adhikari  _et al._ (2007, 2009) or in
Electron Transfer Fedunov  _et al._ (2004); Nicolet  _et al._ (2005); Vauthey
(2006); Mikhailova  _et al._ (2004); Ivanov  _et al._ (2003). The Marcus
Theory of Electron Transfer Marcus (1956, 1957a, 1957b) also postulates that
the attainment of critical configuration of the electrons is due to thermal
fluctuations and hence is independent of the initial excitement. The Keldysh
Field Theory approach to the same is also identical to Marcus Theory upto the
first order of perturbation Hansen and Pullerits (2012).

Newer Theories have been developed to address the same Patra  _et al._ (2011).
Here one studies a projection of the general Liouville Dynamics to a slow
moving co-ordinate of the phase space, namely the Reaction-Coordinate (RC),
where the dynamical equation takes the form of a Heat equation in presence of
a bound potential, thus encapsulating an interplay of the opposing effects of
diffusion and restoring forces in the pseudo energy surface. Dynamics in this
RC approach provides a quantitative justification of the related concept in
theoretical chemistry, as well as is capable of resolving the issue of initial
condition dependent relaxation of various non-equilibrium correlation
functions Bhattacharjee  _et al._ (2019); Choudhury  _et al._ (2020). The
effects of introducing reversibility Gupta  _et al._ (2017); Dhole  _et al._
(2012) and anharmonicity Das  _et al._ (2020) has also been studied. However
all these studies has employed the standard heat equation, with a constant
uniform Diffusion Coefficient.

With the woke of better resolving techniques, there have been recent interest
of exploring heterogenity in free diffusion with the incorporation of a
position dependant diffusion constant Wolfson  _et al._ (2018); Olivares-Rivas
_et al._ (2013); Ljubetič  _et al._ (2014); Nagai  _et al._ (2020). Here we
investigate the same under the RC approach, under a simple Harmonic Bound
Potential. We report that this model exhibits a previously unreported phase of
solvation reactions, which is characterised by a non-monotonic decay of the
non-equilibrium solvation correlation time, as well as of the higher moments.
This phase is non-existent in homogeneous diffusion, in both the for bound and
free cases and is thus uncaptured by the aforementioned theories as well as
remained unexplored by the similar experiments. In the balance of this paper,
we will briefly revisit the relevant theory and then study the parametric
dependence of this non-monotonic phase on the strength of inhomogeneity and
the restoring force. We show that, independent of the restoring force,
inhomogeneity in diffusion is a necessary condition for the existence of this
non-monotonic phase. We construct a phase diagram of this parametric
dependence and indicate the ranges where this non-monotonic phase exists, in
reduced units, with the fondest hope that it would pave the way for its
experimental detection.

##  II Theory

The non-equilibrium solvation processes as is captured by the aforementioned
experiments are essentially a two step process. The femto-second pulse excites
the the solute in presence of the condensed media from a ground state (A) to
an solute+solvent condensed phase (A*). This process is known in the
literature as the Frank-Cordon Excitation. The LASER source is then switched
off and the system is relaxes towards equilibrium. Once the equilibrium is
achieved the system is expected to remain there modulo small fluctuations
around the same.

During the process of relaxation the excited state may undergo a de-excitation
to ground state by emitting fluorescent radiation, which is the experimentally
observed quantity in the spectroscopic experiments. The frequency of this
radiation is related to the Non-Equilibrium Solvation time correlation
function , $S(t)$ as follows,

$S(t)=\frac{\nu(t)-\nu(\infty)}{\nu(0)-\nu(\infty)}=\frac{\overline{\Delta
E}(t)-\overline{\Delta E}(\infty)}{\overline{\Delta E}(0)-\overline{\Delta
E}(\infty)}$ (1)

$\overline{\Delta E}$ is further computed from the Phase space probability
distribution $\rho(\Gamma,t|\Gamma_{0},0)$,

$\overline{\Delta E}(t)=\int
d\Gamma[V_{e}(Q)-V_{g}(Q)]\rho(\Gamma,t|\Gamma_{0},0)$ (2)

Where $V_{g},V_{e}$ respectively are the ground and excited state potential
surfaces and $Q$ represents the set of nuclear coordinate in the phase space.
The dynamics of the phase space density is governed by the Liouville Equation,

$\frac{\partial}{\partial
t}\rho(\Gamma,t|\Gamma_{0},0)=i\mathcal{L}\rho(\Gamma,t|\Gamma_{0},0)$ (3)

where $\mathcal{L}$ is the excited state Liouville operator.

The ground state here is constructed based on energy conservation principle
taking into account the excitation frequency $\lambda$,

$\rho(\Gamma,t|\Gamma_{0},0)=\frac{\exp{(-\beta
V_{e}(Q))}\delta(V_{e}(Q)-V_{g}(Q)-\frac{\hbar c}{\lambda})}{\int
dQ\exp{(-\beta V_{e}(Q))}\delta(V_{e}(Q)-V_{g}(Q)-\frac{\hbar c}{\lambda})}$
(4)

Here our main focus is on the projection of this phase space dynamics on the
one dimensional energy gap coordinate, $\epsilon=V_{e}(Q)-V_{g}(Q)$. This is
obtained by constructing the following projection operator,.

$P(x)=\int d\epsilon\frac{(x*\delta(V_{e}(Q)-V_{g}(Q)-\epsilon))}{\int
d\Gamma\exp{(-\beta H_{g}(\Gamma))}\delta(V_{e}(Q)-V_{g}(Q)-\epsilon)}$ (5)

Employing this projection on the Liuoville Equation gives rise to the Reaction
Coordinate Approach Patra  _et al._ (2011); Zwanzig (1961), and gives rise to
the following effective dynamical equation with a general diffusion kernel
$D(\epsilon,\epsilon^{\prime}|t-t^{\prime})$ and a bound potential
$V(\epsilon)$.

$\begin{split}\frac{\partial}{\partial
t}P(\epsilon,t)=\frac{\partial}{\partial\epsilon}\int_{-\infty}^{\infty}d\epsilon^{\prime}\int_{0}^{t}dt^{\prime}D(\epsilon,\epsilon^{\prime}|t-t^{\prime})\Big{\\{}\frac{\partial}{\partial\epsilon^{\prime}}+\frac{\partial}{\partial\epsilon^{\prime\prime}}V(\epsilon^{\prime\prime})\Bigr{|}_{\begin{subarray}{c}\epsilon^{\prime\prime}=\epsilon^{\prime}\end{subarray}}\Big{\\}}P(\epsilon^{\prime},t^{\prime})\end{split}$
(6)

Following the Markovian assumption of a space-time local Diffusion Kernel
ammounts to the form,

$D(\epsilon,\epsilon^{\prime}|t-t^{\prime})=D(\epsilon)\delta(\epsilon-\epsilon^{\prime})\delta(t-t^{\prime})$
(7)

This converts the general integral-differential equation of (6), to the
following PDE,

$\begin{split}\frac{\partial}{\partial
t}P(\epsilon,t)=\frac{\partial}{\partial\epsilon}D(\epsilon)\Big{\\{}\frac{\partial}{\partial\epsilon}+\frac{\partial}{\partial\epsilon^{\prime\prime}}V(\epsilon^{\prime\prime})\Bigr{|}_{\begin{subarray}{c}\epsilon^{\prime\prime}=\epsilon^{\prime}\end{subarray}}\Big{\\}}P(x,t)\end{split}$
(8)

Our first job is to solve for $P(\epsilon,t)$ numerically from this equation.
We work with a harmonic bound potential $V(\epsilon)=\frac{k}{2}\epsilon^{2}$,
characterised by the restoring force strength $k$ and the position dependent
diffusion coefficient,

$D(\epsilon)=D_{0}(1+\alpha\epsilon^{2})$ (9)

In the following we have studied the numerical solution to 8 and dynamics of
the Non-equilibrium solvation time correlation function $S(t)$ defined as
follows,

$S(t)=\frac{1}{\epsilon_{0}}\int_{-\infty}^{\infty}d\epsilon\;\epsilon
P(\epsilon,t|\epsilon_{0},t_{0})$ (10)

where $\epsilon_{0}$ is the energy gap related to the state of initial
excitation. We exhibit in the following sections that dynamics governed by the
(8) is capable of exhibiting non-monotonic relaxation of $S(t)$ for a wide
range of parameters $\alpha,k$ given a finite $D_{0}$, which arises
exclusively from the inhomogenity of the diffusion constant.

##  III Non-monotonic Relaxation

We start with an initial Gaussian Distribution, centered at $\epsilon_{0}$,
work with vanishing boundary conditions at large $|\epsilon|$ and treat the
Dirac-$\delta$ contribution of source as a limit of Gaussian, with vanishing
variance.

$P(\epsilon_{0},0)=\frac{1}{\sqrt{\pi}\sigma}\exp{\Big{(}\frac{\epsilon-\epsilon_{0}}{\sigma}\Big{)}^{2}}$
(11)

$\begin{split}P(\infty,t)=0\\\ P(-\infty,t)=0\\\ \end{split}$ (12)

For the following studies we have used a $\sigma=0.1$ for the initial
distribution, and for all practical purposes truncated spatial infinities at
$\epsilon=\pm 10$ in reduced units

$\begin{split}P(10,t)=0\\\ P(-10,t)=0\\\ \end{split}$ (13)

We numerically solve the PDE using a dynamic grid in Mathematica with maximal
grid spacing $h=0.01$ in both spatio-temporal directions. Equation (8)
describes time evolution of states in various non-equilibrium processed in
condensed phases such as solvation. In the case of solvation the physically
interesting quantity which characterises the out of equilibrium effects is the
Non-Equilibrium Solvation time correlation function,

$S(t)=\frac{1}{\epsilon_{0}}\int_{-\infty}^{\infty}d\epsilon\;\epsilon
P(\epsilon,t|\epsilon_{0},t_{0})$ (14)

At large time it is expected to relax at its equilibrium value. Figure (1(a))
reports dynamics of $S(t)$ for various initial conditions ($\epsilon^{0}$) for
three combinations of ($\alpha,k$). The first case with $\alpha=0$ corresponds
to the case of position independent diffusion. In this case the result is well
known, and $S(t)$ exhibits an exponential decay, with a relaxation coefficient
independent of the position of the initial wavelet.

Introduction of spatial in-homogeneity lifts this degeneracy (Figure 1(b)).
With $\alpha=1,k=1$, the correlation time corresponding to wavelets initiated
at $\epsilon^{0}=1,2,3,4$ reduced units respectively exhibits exponential
relaxations with distinct rates. Further away from the ground state the system
is prepared, faster is its relaxation to equilibrium value. However, at the
intermediate range, $\alpha=0.25,k=0.25$ (Figure 1(b)), for some initial
conditions closer to the ground state we see a non-monotonic relaxation of the
Non-equilibrium relaxation time.

(a) $\alpha=0,k=1$

(b) $\alpha=0,k=1$

(c) $\alpha=0.25,k=0.25$

Figure 1: Comparison plot of the Non-equilibrium solvation correlation time
$S(t)$ vs $t$ for various initial energy $\epsilon^{0}$, for the indicated
combinations of $(\alpha,k)$. Figure 1a. shows the limiting case of position
independent Diffusion Constant in a Harmonic Potential. From Linear Response
Theory the relaxation of $S(t)$ to the equilibrium value in this case is known
to be exponential and independent of the initial energy level
($\epsilon^{0}$). Figure 1b. exhibits that introduction of position dependence
in diffusion constant breaks this universality and the relaxation is faster
for higher initial energy levels. Figure 1c. reports the same for an
intermediate combination $(\alpha=0.25,k=0.25)$. This case however depicts a
non-monotonic decay, where $S(t)$ first increases, attains a maximum value
$S_{Max}$ and then relaxes to the equilibrium value, for a large range of
initial energy levels. We identify it as a new phase and delve into
characterizing the same in the balance of this paper.

This effect forms the thesis of our work presented here. As as is depicted
later in the study, this phase is prevalent for a wide range of initial
conditions, spatial inhomoeneity of the Diffusion Coefficient as well as the
restoring strengths. This effect is also robust in the higher moments of the
energy difference. In Figure (2), we compute the next six moments for
$\epsilon^{0}=0.2k=1$ and for various $\alpha$ (1,4,7,10), where the first
moment is known to exhibit the non-monotonic relaxation. It is shown that all
the odd moments also exhibit non-monotonic relaxations for the same
parameters, whereas the even moments exhibit an exponential relaxtion.

(a) $<\epsilon^{3}>$

(b) $<\epsilon^{5}>$

(c) $<\epsilon^{7}>$

(d) $<\epsilon^{2}>$

(e) $<\epsilon^{4}>$

(f) $<\epsilon^{6}>$

Figure 2: Comparison of various higher moments of the energy gap vs t for
different alpha (1,4,7,10 blue to orange) For various initial wavelet as
indicated K 1

Figure (3) summarises the behavior of S(t) for various initial preparations
$\epsilon^{0}=0.2,0.5,0.8,1.3$ for $k=1$ as the inhomogenity, $alpha$ is
varied. It is seen that the effect is most prominent if the initial condition
is prepared near ground state and diminishes with the gap, and beyond a
certain finite value of the initial condition becomes monotonic. The effect
becomes more prominent with increasing diffusion. In the following section we
delve in studying this parametric dependence in detail and construct a phase
portrait that indicates the existence of the non-monotonic phase in the
parameter space ($\alpha,k$) as well as the relative strength.

(a) $\epsilon^{0}=0.2$

(b) $\epsilon^{0}=0.5$

(c) $\epsilon^{0}=0.8$

(d) $\epsilon^{0}=1.3$

Figure 3: Solvation: Comparison of S(t) vs t for different alpha (1,4,7,10
blue to orange) For various initial wavelet as indicated K 1

##  IV The Phase Portrait

Figure (3), in previous section portrays dependence of the strength of the
non-monotonic behaviour on the strength of the inhomogenity in diffusion
($\alpha$), and the restoring force ($k$), for various initial conditions
$\epsilon^{0}$. In this section we wish to concentrate further on the same and
attempt to quantify the behaviour of $S(t)$ for a wide range of these
parameters $\alpha,k$.

To begin with we notice the following property that demarcates the non-
monotonic phase from the usual monontonic decay. For the parametric
combination of $\alpha,k$ where the former exists, it shows an initial growth
in $S(t)$ beyond the universal initial value of 1, attains a maximum (to be
denoted by $S_{Max}$) and then further relaxes to the equilibrium value. Thus
existence of a $S_{Max}$ strictly $>1$ for the initial conditions
$\epsilon^{0}$, denotes a non-monotonic decay, as opposed to the monotonic
case, where $S_{Max}$ equals 1 irrespective of all other details.

In figure (4), we plot $S_{Max}$ as a function of $\epsilon^{0}$ in the range
[0,2] for various strengths of $\alpha(1,3,5,7)$ and $k(0,0.5,1,4)$. It is
seen that the non-monotone nature of the relaxation is strengthened with
increase in inhomogeneity and decrease in restoring strength as is denoted by
larger values of $S_{Max}$ for given initial condition, as we vary these
parameters. The effect is maximised in the free diffusion case $k=0$ as long
as we have a finite inhomogeneity $\alpha$. We also notice for a given
combination of $\alpha,k$ the range of initial conditions $\epsilon^{0}$, for
which the non-monotonic phase is exhibited has a compact support between 0 and
a certain $\epsilon^{0}_{c}ritical$, beyond which the decay is again
monotonic.

This forms our basis of construction of the phase space portrait, Figure (5).
For the parameter ranges of ($\alpha,k$), that does not exhibit the non
monotonic phase for any initial condition $\epsilon^{0}$, the
$\epsilon^{0}_{c}ritical$ is evaluated to be zero (denoted by violet). This
portrait gives detailed account of the existence of the non-monotonic phase
and its relative strength as a parametric dependence of inhomogenity and
restoring strength. Firstly we observe that $\epsilon^{0}_{c}ritical$ is
identically zero in the line $\alpha=0$ which corresponds to the space
independent diffusion. Most of the theoretical as well as experimental studies
till date has been concentrated in this regime; the reason for the non
monotonic decay to be not observed before. The phase portrait also shows the
robustness of the existence of this phases across the entire spectrum of
$\alpha,k$ and its relative strength has been indicated by a temperature map,
as well as the phase boundary with the strictly monotonic phase. These
parameters have been indicated in reduced units, and can be converted to
realistic values. Figure (5) thus forms a guide for the experiments for the
range of $\alpha,k$ to look at for the detection of this predicted phase.

Figure 4:  Figure 5: Phase Portrait of $\epsilon_{critical}^{0}$ against
various values of the position dependence ($\alpha$) and the potential
strength ($k$), which exhibits the two phases and there phase boundaries.
Positive values of $\epsilon_{critical}^{0}$ (colors above violet) denotes
existence of the non-monotonic phase, while a higher value (redder color)
indicates its abundance for a given combination of ($\alpha,k$). The effect is
stronger for higher diffusion ($\alpha$) and lower potential strength
(decreasing $k$), which reasserts its origin in interplay of these two
opposing force. $\epsilon_{critical}^{0}$ is identically zero for the
$\alpha=0$ line, irrespective of potential strength $k$. This corresponds to
the space-independent case, where we expect an exponential decay from LRT .
The non-monotonic phase is an exclusive feature of introducing non-homogeneous
diffusion.

##  V Conclusion

In summary, we have studied the effects of introducing position-dependent
diffusion in non equilibrium dynamics in a bound harmonic potential. The
position dependence introduces non-monotonicity in the relaxation of non-
equilibrium correlation functions for a wide range of combinations of
inhomogenity and restoring strength. With higher inhomogenity in diffusion as
well as lower restoring forces this feature is strengthened. This feature is
robust in relaxation of all higher odd moments of the energy gap, where the
even moments still exhibit an exponential relaxation. We have studied the
criticality and constructed a phase portrait which shows that it requires a
certain minimum positional inhomogeneity in diffusion for this non-
monotonicity to be induced, which is consistent with the exponential
relaxation of correlation functions in the position-independent case. Most
theoretical as well as experimental works till date have been performed in
this limit, causing this phenomena in the context of reactions of condensed
phase remain unexplored. The fondest expectation of our reportage of the same,
as well as predictions of the parameter range indicating its abundance is to
pave way for experiments for for its unearthing.

## References

  * Onsager (1931a) L. Onsager, Phys. Rev. 37, 405 (1931a).
  * Onsager (1931b) L. Onsager, Phys. Rev. 38, 2265 (1931b).
  * Sen  _et al._ (2006) P. Sen, S. Ghosh, K. Sahu, S. K. Mondal, D. Roy, and K. Bhattacharyya, The Journal of chemical physics 124, 204905 (2006).
  * Mondal  _et al._ (2006) S. K. Mondal, S. Ghosh, K. Sahu, U. Mandal, and K. Bhattacharyya, The Journal of chemical physics 125, 224710 (2006).
  * Mandal  _et al._ (2008) U. Mandal, S. Ghosh, D. K. Das, A. Adhikari, S. Dey, and K. Bhattacharyya, Journal of Chemical Sciences 120, 15 (2008).
  * Ghatak  _et al._ (2011) C. Ghatak, V. G. Rao, S. Ghosh, S. Mandal, and N. Sarkar, The Journal of Physical Chemistry B 115, 12514 (2011).
  * Adhikari  _et al._ (2007) A. Adhikari, K. Sahu, S. Dey, S. Ghosh, U. Mandal, and K. Bhattacharyya, The Journal of Physical Chemistry B 111, 12809 (2007).
  * Adhikari  _et al._ (2009) A. Adhikari, D. K. Das, D. K. Sasmal, and K. Bhattacharyya, The Journal of Physical Chemistry A 113, 3737 (2009).
  * Fedunov  _et al._ (2004) R. G. Fedunov, S. V. Feskov, A. I. Ivanov, O. Nicolet, S. Pagès, and E. Vauthey, The Journal of chemical physics 121, 3643 (2004).
  * Nicolet  _et al._ (2005) O. Nicolet, N. Banerji, S. Pages, and E. Vauthey, The Journal of Physical Chemistry A 109, 8236 (2005).
  * Vauthey (2006) E. Vauthey, Journal of Photochemistry and Photobiology A: Chemistry 179, 1 (2006).
  * Mikhailova  _et al._ (2004) V. A. Mikhailova, A. I. Ivanov, and E. Vauthey, The Journal of chemical physics 121, 6463 (2004).
  * Ivanov  _et al._ (2003) A. Ivanov, F. Belikeev, R. G. Fedunov, and E. Vauthey, Chemical physics letters 372, 73 (2003).
  * Marcus (1956) R. A. Marcus, The Journal of chemical physics 24, 966 (1956).
  * Marcus (1957a) R. A. Marcus, The Journal of Chemical Physics 26, 867 (1957a).
  * Marcus (1957b) R. Marcus, The Journal of Chemical Physics 26, 872 (1957b).
  * Hansen and Pullerits (2012) T. Hansen and T. Pullerits, Journal of Physics B: Atomic, Molecular and Optical Physics 45, 154014 (2012).
  * Patra  _et al._ (2011) A. Patra, A. Samanta, and S. K. Ghosh, Phys. Rev. E 83, 026104 (2011).
  * Zwanzig (1961) R. Zwanzig, Phys. Rev. 124, 983 (1961).
  * Bhattacharjee  _et al._ (2019) S. Bhattacharjee, A. Patra, A. K. Pathak, and A. K. Samanta, Physica Scripta 94, 085008 (2019).
  * Choudhury  _et al._ (2020) S. Choudhury, A. Patra, A. K. Pathak, and A. K. Samanta, Molecular Physics 118 (2020), 10.1080/00268976.2020.1779365.
  * Gupta  _et al._ (2017) K. Gupta, A. Patra, K. Dhole, A. K. Samanta, and S. K. Ghosh, The Journal of Physical Chemistry Letters 8, 4545 (2017), pMID: 28872879\.
  * Dhole  _et al._ (2012) K. Dhole, N. K. Jena, A. Samanta, and S. K. Ghosh, Phys. Rev. E 85, 026105 (2012).
  * Das  _et al._ (2020) M. Das, A. Samanta, and S. K. Ghosh, Theoretical Chemistry Accounts 139, 1 (2020).
  * Wolfson  _et al._ (2018) M. Wolfson, C. Liepold, B. Lin, and S. A. Rice, The Journal of chemical physics 148, 194901 (2018).
  * Olivares-Rivas  _et al._ (2013) W. Olivares-Rivas, P. J. Colmenares, and F. López, The Journal of chemical physics 139, 074103 (2013).
  * Ljubetič  _et al._ (2014) A. Ljubetič, I. Urbančič, and J. Štrancar, The Journal of chemical physics 140, 084109 (2014).
  * Nagai  _et al._ (2020) T. Nagai, S. Tsurumaki, R. Urano, K. Fujimoto, W. Shinoda, and S. Okazaki, Journal of Chemical Theory and Computation 16, 7239 (2020).
  * Patra  _et al._ (2015) A. Patra, K. A. Acharya, and A. Samanta, The Journal of Physical Chemistry B 119, 11063 (2015).


# Metal-Insulator transition in strained Graphene: A quantum Monte carlo study

Lufeng Zhang<EMAIL_ADDRESS>School of Science, Beijing University of
Posts and Telecommunications, Beijing 100876, China  Chi Ma  Department of
Physics, Beijing Normal University, Beijing 100875, China  Tianxing Ma
<EMAIL_ADDRESS>Department of Physics, Beijing Normal University, Beijing
100875, China

###### Abstract

Motivated by the possibility of a strain tuning effect on electronic
properties of graphene, the semimetal-Mott insulator transition process on the
uniaxial honeycomb lattice is numerically studied using Determinant Quantum
Monte Carlo. As our simulations are based on the half-filled repulsive Hubbard
model, the system is sign problem free. Herein, the temperature-dependent DC
conductivity is used to characterize electronic transport properties. The data
suggest that metallic is suppressed in the presence of strain. More
interestingly, within the finite-size scaling study, a novel antiferromagnetic
phase arises at around $U\sim U_{c}$. Therefore, a phase diagram generated by
the competition between interactions and strain is established, which may help
to expand the application of strain effect on graphene.

##  I Introduction

Graphene is a 2D sheet of carbon atoms packed hexagonal structure as depicted
in Figure. 1. Since its discoveryNovoselov  _et al._ (2004), graphene has
stimulated a tremendous burst of research activities on fundamental and
applied grounds. Graphene has multiple astonishing performance in several
areas, which is the toughest 2D material ever measured Lee  _et al._ (2008),
and also has extremely high charge carriers mobility (charge carriers move at
a speed of $c/300$ ($c$ is the speed of light in vacuum)Novoselov  _et al._
(2005); Morozov  _et al._ (2008) and behave like massless particles) More than
these, graphene is impermeable to standard gasesBunch  _et al._ (2008) and it
could be a perfect thermal conductorBalandin  _et al._ (2008). Numerous other
honeycomb liked systems have attracted a lot of attention, such as
siliceneVogt  _et al._ (2012), germaneneCahangirov  _et al._ (2009) and
phosphoreneLi  _et al._ (2014), which also have excellent properties like
graphene. Thus, we expect that these honeycomb like materials will pave the
way toward a new age of material.

There is no gap in the band structure of graphene which exhibits semi-metallic
propertiesNovoselov  _et al._ (2008); Peres (2010). However, whether the
energy gap could be modified in real graphene material to realize the metal-
insulator transition (MIT)Drut and Lähde (2009); Schüler  _et al._ (2013);
Ulybyshev  _et al._ (2013) is still controversial. If it is possible to
control the opening and closing of the energy gap in real graphene, the
application of graphene in the electronic device will be greatly promoted.

Recent theoretical and experimental work shows that hydrogenationElias  _et
al._ (2009), application of stressNi  _et al._ (2008); Choi  _et al._ (2010);
Pletikosić  _et al._ (2009); Meunier  _et al._ (2016); Pereira  _et al._
(2009), forming nanoribbonsHan  _et al._ (2007), Moir$\acute{e}$-stripe
modulationHunt  _et al._ (2013); Ponomarenko  _et al._ (2011), etc. can be
used to open the band gap and achieve MIT in graphene. In this work, we
explore different approaches to achieving insulating state, considering the
interactions between electrons. Graphene is clearly a semimetal without any
electron-electron interactions, but it will turn into an insulating
antiferromagnet(AFM) under very strong interactionsKotov  _et al._ (2012).
However, it still remains unclear as to what we should expect in real graphene
materials. For instance, in the real case, it is common to find stress effect
in spite of interaction effects, and it is an amazing thing that the
electronic properties could be modified by mechanical forcesGui  _et al._
(2015); Aslani  _et al._ (2015); Fanbanrai  _et al._ (2015); Arya  _et al._
(2017); Choi  _et al._ (2009); Cosma  _et al._ (2014); Craco  _et al._ (2015).
As for graphene, defects could be introduced naturally or intentionally while
growing. The existence of these defects also attract wide interestSummerfield
_et al._ (2016); Zhang  _et al._ (2016). In the experiment, strain could be
measured and detected by atomic force microscopy. And strain of $25\%$ on the
graphene can be achievedLee  _et al._ (2008), so that the system could
completely come into the antiferromagnetic regionLee  _et al._ (2008).
Theoretically, there are also a lot of work based on the tight-binding model,
studying the effect of uniaxial strain on the properties of graphenePereira
_et al._ (2009). It seems that the gap suddenly arises when deformation is
beyond $20\%$ without a transition regionPereira  _et al._ (2009); Choi  _et
al._ (2010). To illustrate the strain effect in these heavy fermion materials,
the electron interactions should also be considered, and it is interesting to
find that the quantum critical point (QCP) strain value is suppressed by
interactionsTang  _et al._ (2015). The localization effect and mutual effect
of these two factors are the key to understanding the mechanism of MIT in
strained graphene. In this paper, we will focus on the competition between
electronic correlations and uniaxial strains, and the phase diagram is also
discussed.

Figure 1: The geometry of the $L=6$ strained honeycomb lattice with 72 sites,
where red and blue label A and B sublattice, respectively. Stress is added
along the x-direction as shown by red lines, which $t_{2}=t-\Delta{t}$. The
dark lines indicate $t_{1,3}=t$. Here $t$ is the nearest hopping integrating
term and $\Delta{t}$ represents the strength of strain.

##  II Model and numerical method

The Hubbard modelHubbard (1963) is believed to be a versatile paradigm to
study strongly correlated electrons on a lattice. It is effective to describe
electrons in partially filled energy bands. The parameters of the model, such
as hopping and electron-electron interaction, can be considered to be
generated by integrating the corresponding degrees of freedom of all energy
bands except partially filled narrow bands. Thus, this modelImada  _et al._
(1998) can be used to explain the interaction-driven transition between
conducting and insulating behaviors.

To shed light on the physical problems involved in the strained-correlation
fermions, we use determinant quantum Monte Carlo (DQMC) to study the half-
filled repulsive Hubbard model on the honeycomb lattice. In the Dirac
electronic correlation system, it is already known that there exist a QCP for
MIT at around $Uc\sim 3.89t$Paiva  _et al._ (2005); Sorella  _et al._ (2012);
Assaad and Herbut (2013); Arya  _et al._ (2015). In this strength of
interaction region, our numerical DQMC method can conduct an accurate
simulation. We consider the strained Hubbard Hamiltonian on the honeycomb
lattice as

$\displaystyle H=$
$\displaystyle-\sum_{\mathbf{i}\eta\sigma}t^{\phantom{\dagger}}_{\eta}a^{\dagger}_{\mathbf{i}\sigma}b^{\phantom{\dagger}}_{\mathbf{i}+\eta\sigma}+h.c.-\mu\sum_{\mathbf{i}\sigma}\left(n_{\mathbf{a}\mathbf{i}\sigma}+n_{\mathbf{b}\mathbf{i}\sigma}\right)$
(1)
$\displaystyle+U\sum_{\mathbf{i}}\left(n_{\mathbf{a}\mathbf{i}\uparrow}n_{\bf{ai}\downarrow}+n_{\bf{bi}\uparrow}n_{\bf{bi}\downarrow}\right).$

Here $a^{\dagger}_{\mathbf{i}\sigma}$ ($a_{\mathbf{i}\sigma}$) are the
spin-$\sigma$ electron creation (annihilation) operators at site $\mathbf{i}$
on sublattice A. And $b^{\dagger}_{\mathbf{i}\sigma}$ ($b_{\mathbf{i}\sigma}$)
are operators acting on sublattice B. $U>0$ is the on-site Coulomb repulsion.
$t_{\eta}$ denotes the hopping integral between two nearest-neighbor sites.
The chemical potential $\mu$ determines the average density of the system.
$n_{\bf{ai}\sigma}=a^{\dagger}_{\mathbf{i}\sigma}a^{\phantom{\dagger}}_{\mathbf{i}\sigma}$
and
$n_{\bf{bi}\sigma}=b^{\dagger}_{\mathbf{i}\sigma}b^{\phantom{\dagger}}_{\mathbf{i}\sigma}$
are the number operators. Strain is introduced through the hopping matrix
elements $t_{\eta}$ along the x-direction which is shown by red lines in Fig.
1. $t_{1,3}=t$ and $t_{2}=t-\Delta{t}$. Here $\Delta{t}$ represents a measure
of stress strength. We take $t=1$ as the scale of energy in our system. In
this work, we set $\mu=0$, thus, the system is half-filled by electrons. And
in this case, the Hamiltonian will be particle-hole symmetric even with strain
on it, so that the Hamiltonian is sign-problem free to solve in our DQMC
simulations.

In the finite-temperature DQMC method, it follows the strategy that take the
partition function as an integral over all possible configurations which is
carried out by random Monte Carlo sampling. To be more specifically, we can
get static and dynamic observables at an given temperature $T$ through this
approach. As the system is sign-problem free due to the particle-hole
symmetry, our calculations can still be in high numerical precision at large
enough $\beta=1/T$, thus, we can converge the data to get the ground state. In
this work, 4000 sweeps were used to get to the equilibrium state, and 12000
additional steps were made to measure the system. The periodic boundary
conditions were used in the simulation, and all the results are obtained on
the $2\times 12^{2}$ honeycomb lattice with is large enough to consider the
size effect. Fig. 1 shows the $L=6$ geometry with strain on it.

Figure 2: Phase diagram of the strained Hubbard model on the honeycomb lattice
at half-filling. $\Delta t$ labels the strain strength and $U$ represents the
local Coulomb repulsion. The metallic phase boundary is determined by the
temperature dependence of the conductivity $\sigma_{dc}$. The metallic region
is filled with green color and the (grey) hexagonal mark point is obtained
from our DQMC simulation results. The blue region on the left indicates the
existence of long-range anti-ferromagnetic order.

We mainly discussed the MIT in this article. The temperature-dependent DC
conductivity $\sigma_{dc}(T)$ is an intuitive physical observation to help
characterize the MIT. So that we focus on discussing $\sigma_{dc}(T)$ behavior
in our case. According to the fluctuation-dissipation theorem, in the zero
frequency limit, $\sigma_{dc}(T)$ is related to the current-current
correlation function. As we could do imaginary-time simulations using DQMC,
the real-frequency quantities can be easily obtained through analytic
continuation methods. In our case, we adopted an approximation Trivedi and
Randeria (1995) to define $\sigma_{dc}(T)$, which has been widely used as a
benchmark in previous work Trivedi  _et al._ (1996); Trivedi and Randeria
(1995); Denteneer  _et al._ (1999a).

The definition of DC conductivity is:

$\sigma_{dc}(T)=\frac{\beta^{2}}{\pi}\Lambda_{xx}({\bf q}=0,\tau=\beta/2).$
(2)

In Equation (2), $\Lambda_{xx}$ is the current-current correlation function,
and it can be expressed as $\Lambda_{xx}({\bf q},\tau)={\langle j_{x}({\bf
q},\tau)\,j_{x}(-{\bf q},0)\rangle}$, where $j_{x}({\bf q},\tau)$ is the
Fourier transform of $j_{x}(\mathbf{r},\tau)$ (time-dependent current
operator) along the $x$-direction. When we focus on studying the system in the
temperature lower than the energy scale, in which the density of states has a
significant structure, the DC conductivity definition, in Equation(2), can be
very practicable considering such approximation. The applicability has been
checked in our recent work of the Hubbard model on the honeycomb latticeMa
_et al._ (2018). We can get the MIT critical strength value $U_{c}$Denteneer
_et al._ (1999b) within the formula from the change of low-$T$ behavior of
$\sigma_{dc}(T)$. We also investigate the magnetic properties by carrying out
measurements of spin-spin correlation functions and the spin structure factor,

$S_{AF}=\frac{1}{N_{c}}\langle(\sum_{\bf r\in A}\hat{S_{\bf r}^{z}}-\sum_{\bf
r\in B}\hat{S_{\bf r}^{z}})^{2}\rangle.$ (3)

where $N_{c}$ is the number of sites in the lattice, and $A$ and $B$
represents the sublattices of the honeycomb lattice. Here $\hat{S_{\bf
r}^{z}}$ is the $z$ component spin operator. The results are normalized by
$\langle...\rangle$ including all the lattice sites. The antiferromagnetic
results present in this work are concluded from the constrained-path quantum
Monte Carlo (CPQMC) methodZhang  _et al._ (1997, 1995); Ma  _et al._ (2011),
which uses the constrained-path approximation in the simulation procedure
without a sign problem.

Figure 3: Temperature dependence of the DC conductivity $\sigma_{dc}$ measured
on the $L=12$ lattice with strain along the x-direction. $N=2\times L^{2}$ is
the total site numbers. Panels correspond to different couplings: (a) $U=1.0$,
(b) $U=2.0$, (c) $U=3.0$ and (d) $U=4.0$. In each figure, lines are guides to
the eyes. $\Delta{t}$ represents for the strain strength. Metallic and
insulating behaviors are indicated by solid and dashed lines respectively. The
MIT behavior can be easily observed from the low-$T$ behavior of
$\sigma_{dc}$, whereas in panel (d) the insulating behavior is preserved by
strain effect.

##  III Results and discussion

Our simulations are mainly computed on the $L=12$ lattice considering the size
effect. $N=2\times L^{2}$ is the total site numbers. The quantity of immediate
interest in this work is to study the possible MIT of the strained honeycomb
lattice, which can be measured by the $conductivity$ and its $T$ dependence
behavior. So we show the $\sigma_{dc}(T)$ behavior at low temperature here in
Figure. 3. For some curves, $d\sigma_{dc}(T)/dT<0$ and $\sigma_{dc}(T)$
diverges as $T$ decreases to the $T\rightarrow 0$ limit, which suggesting
metallic behavior of the system. In contrast, the system is insulating when
$d\sigma_{dc}(T)/dT>0$ at $T\rightarrow 0$ limit. Thus, we can easily tell the
MIT point of the strained system. In panels (a)-(b) of Fig. 3, we described
$\sigma_{dc}(T)$ behaviors at low coupling strengths under several strain
strengths. To understand the role of strain on the conductivity, we can take a
careful analyzation of panel (a) in Fig. 3. In spite of $\Delta{t}$ strength,
the conductivity increases as temperature increases until $T\succeq 0.25$.
However when $U=1.0$ and with a strong enough strain strength
$\Delta{t}=0.20$, the curve bent down when temperature drops and
$\sigma_{dc}(T)$ approaches zero as $T\rightarrow 0$, suggesting insulating
behavior of the system. In contrast, when we take a look at $\Delta{t}=0.10$,
the red line diverges with feature of $d\sigma_{dc}(T)/dT<0$ for $T\leq 0.2$,
characteristic the metallic behavior of the system. Thus, we can clearly get
the strain-driven MIT critical point from the panels. At $U=1.0$, the critical
$\Delta{t}$ value can be read as $\Delta{t_{c}}\sim 0.15\pm 0.01$. As for
panel (b) in Fig. 3 with $U=2.0$. At $\Delta{t}=0.12$, $\sigma_{dc}(T)$ curve
is concave and $d\sigma_{dc}(T)/dT<0$ for $T\leq 0.2$. When $T$ drops to
around $0.1$, $\sigma_{dc}(T)$ increases rapidly, in other words, the system
is metallic. At $\Delta{t}=0.14$, in contrast, the conductivity decreases as
the temperature decreases, which indicating insulating behaviors. Thus, there
exist a MIT critical strain strength at around $\Delta{t}=0.14\pm 0.01$ for
$U=2.0$ on the graphene lattice. Then it moves to the larger interaction case,
as shoen in Fig. 3(c) and (d). When $U=3.0$, the system is still metallic
without considering strain effect. So that under small strain strength, it is
still possible that the graphene can stay in metallic region. At low $T$ for
$T<0.1$, the sign of $d\sigma_{dc}(T)/dT$ is negative when $\Delta{t}\leq
0.12$, but it switches to positive when $\Delta{t}$ is larger than $0.13$.
Therefore, the MIT critical $\Delta{t}$ strengths can be read as
$\Delta{t_{c}}\sim 0.13\pm 0.01$ for $U=3.0$. As for $U=4.0$ which is larger
than the critical $U_{c}=3.89$ for MIT in the graphene lattice without strain,
the system is an insulator. From panel (d), we can see all the
$\sigma_{dc}(T)$ curves are bent down when temperature decreased for
$T\rightarrow 0$. The value of $\sigma_{dc}(T)$ gets suppressed when strength
of strain increases, that is to say the insulating behavior is preserved by
the strain effect in this case. To describe the metal-insulator phase boundary
more precisely, we closely studied the conductivity behaviour around $U\sim
3.8$. Thus, according to Fig. 6 in Appendix, we make the metal-insulator phase
boundary more smooth, as shown in the phase diagram Fig. 2.

Figure 4:  Temperature dependence of the DC conductivity $\sigma_{dc}$
measured on the $L=12$ lattice with strain along different directions. Panels
correspond to different couplings and strain strengths: (a) $U=2.0$,
$\Delta{t}=0.1$ is in the metallic region, and (b) $U=4.0$, $\Delta{t}=0.1$,
(c) $U=2.0$, $\Delta{t}=0.3$, (c) $U=4.0$, $\Delta{t}=0.3$ are in the
insulating region. In each figure, lines are guides to the eyes. $\Delta{t}$
represents for the strain strength, $\alpha$ is the inclined angle between
strain and $x$-axis.

In Fig. 4, we also investigate the effect of the orientation of the strain.
Here $\alpha$ in each panel is the angle with x-axis. Panel (a) is in the
metallic regime for $U=2.0$ and $\Delta{t}=0.1$. As $\Delta{t}$ increases to
$0.3$, the system turns into the insulating regime, as panel (c) shows.
Comparing panel (b) and (d), the result at $U=4.0$ is shown that graphene
stays in the insulating state no matter how strong the strain applied.
According to the lines in Fig. 4, it seems that the behavior of $\sigma_{dc}$
is pretty similar when $\alpha$ changes. We only consider the nearest hopping
in our Hamiltonian, which is a single-shell model. The shell behavior to the
angles is very similar to the first shell of the model including up to third-
nearest-neighbor interactionsBotello-Méndez  _et al._ (2018). The angular
dependence is not very important considering atoms in the second and third
shells, thus, we get the similar behaviors of $\sigma_{dc}(T)$ here with
different orientations of strain. Therefore, we can focus our study on the
strength of strain and ignore its orientation when exploring the strain effect
on MIT.

Further more, we also studied the effect of strain on long-range magnetic
order. Data shown in Fig. 5 and Fig. 7 are AF spin structure factor on
lattices up to $L=15$. For small $U=2,3$, there is no long-range order when
$U<U_{c}\sim 3.89t$Paiva  _et al._ (2005); Sorella  _et al._ (2012); Assaad
and Herbut (2013); Arya  _et al._ (2015). Even when $\Delta{t}=0.30t$, there
is no AF order according to Fig. 7(a) and (b). As for small lattices, the
correlation strength of antiferromagnetism decreases while strain increases,
which indicates that in the weak interaction case, $S_{AF}$ are suppressed by
strain. In contrast, in Fig. 7(c) and (d), AF correlation strength is enhanced
by strain at large $U$. The AF order is protected by strain at strong
interaction ($U>U_{c}$). As Fig. 5 shows, it is interesting to note that the
AF order could exist even when $U<U_{c}$. Thus, there is a transition point
that exists for long-range order by strain effect. We obtained
(a)$U=3.2,\Delta{t}=0.30$, (b)$U=3.5,\Delta{t}=0.25$,(a)$U=3.6,\Delta{t}=0.25$
and (a)$U=3.8,\Delta{t}=0.20$ from the finite size scaling results. Based on
the analysis earlier, the paramagnetic-antiferromagnetic transition boundary
could be summarised for nonzero $\Delta{t}$, as shown in Fig. 2.

Figure 5:  Finite-size scaling results of the AF spin structure factor by
CPQMC. $S_{AF}$ is plotted as a function of $1/L$ corresponding to various
value of interaction strength. A finite $S_{AF}$ in the $L\rightarrow\infty$
limit indicates the existence of long-range magnetic order. Symbols on the
lines are simulation results on $L=3,6,9,12,15$ lattice sizes within
statistical errors. The lines are cubic polynomial fits to the data.

##  IV Summary

We have carefully studied the electronic and magnetic properties of the
strained Hubbard model on the honeycomb lattice. Using DQMC, we characterize
the $conductivity$ with the $\sigma_{dc}(T)$ behavior at low temperature to
explore the mechanism of strain-driven MIT in the honeycomb lattice. According
to our simulations, it is clear to see that the conductivity $\sigma_{dc}$ is
suppressed by strain effect. When strain is applied, the value of hopping
parameter along the strain direction decreases from $t$ to $t-\Delta{t}$,
thus, the critical electron-electron interaction strength $U$ for MIT reduces.
And the critical strain strength decreases as the local Coulomb repulsion $U$
enhanced. In other words, under the interaction of strain and Coulomb
correlation, the electrons get localized more effectively, which expands the
insulating phase and gives us a strain-interaction driven MIT phase diagram.

As for the magnetic properties, we found that strain could also enhance the AF
order at around $U_{c}$ in the clean limit. While strain is affected, we
observed the appearance of the antiferromagnetic order phase. It seems that
strain can protect the long antiferromagnetic order when $U$ is strong. In
conclusion, we show that strain-interaction-driven graphene could be a
promising rout to achieve MIT in graphene, and the phase diagram reported in
this work could be used as a guidance in the modulation of conductivity in
real graphene materials.

T. M. and L. F. Z were supported by NSFCs (No. 11774033 and 11974049) and
Beijing Natural Science Foundation (No. 1192011).

## References

  * Novoselov  _et al._ (2004) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * Lee  _et al._ (2008) C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science 321, 385 (2008).
  * Novoselov  _et al._ (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, M. I. Jiang, D.and Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
  * Morozov  _et al._ (2008) S. V. Morozov, K. S. Novoselov, M. I. Katsnelson, F. Schedin, D. C. Elias, J. A. Jaszczak, and A. K. Geim, Phys. Rev. Lett. 100, 016602 (2008).
  * Bunch  _et al._ (2008) J. S. Bunch, S. S. Verbridge, J. S. Alden, A. M. Van, der Zande, J. M. Parpia, H. G. Craighead, and P. L. Mceuen, Nano Letters 8, 2458 (2008).
  * Balandin  _et al._ (2008) A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao, and C. N. Lau, Nano Letters 8, 902 (2008).
  * Vogt  _et al._ (2012) P. Vogt, P. De Padova, C. Quaresima, J. Avila, E. Frantzeskakis, M. C. Asensio, A. Resta, B. Ealet, and G. Le Lay, Phys. Rev. Lett. 108, 155501 (2012).
  * Cahangirov  _et al._ (2009) S. Cahangirov, M. Topsakal, E. Aktürk, H. Şahin, and S. Ciraci, Phys. Rev. Lett. 102, 236804 (2009).
  * Li  _et al._ (2014) L. Li, Y. Yu, G. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. Chen, and Y. Zhang, Nature Nanotechnology 9, 372 (2014).
  * Novoselov  _et al._ (2008) K. Novoselov, P. Blake, and M. Katsnelson, in  _Encyclopedia of Materials: Science and Technology_, edited by K. J. Buschow, R. W. Cahn, M. C. Flemings, B. Ilschner, E. J. Kramer, S. Mahajan, and P. Veyssière (Elsevier, Oxford, 2008) pp. 1 – 6.
  * Peres (2010) N. M. R. Peres, Physica Status Solidi 244, 4106 (2010).
  * Drut and Lähde (2009) J. E. Drut and T. A. Lähde, Phys. Rev. Lett. 102, 026802 (2009).
  * Schüler  _et al._ (2013) M. Schüler, M. Rösner, T. O. Wehling, A. I. Lichtenstein, and M. I. Katsnelson, Phys. Rev. Lett. 111, 036601 (2013).
  * Ulybyshev  _et al._ (2013) M. V. Ulybyshev, P. V. Buividovich, M. I. Katsnelson, and M. I. Polikarpov, Phys. Rev. Lett. 111, 056801 (2013).
  * Elias  _et al._ (2009) D. C. Elias, R. R. Nair, T. M. G. Mohiuddin, S. V. Morozov, P. Blake, M. P. Halsall, A. C. Ferrari, D. W. Boukhvalov, M. I. Katsnelson, A. K. Geim, and K. S. Novoselov, Science 323, 610 (2009).
  * Ni  _et al._ (2008) Z. H. Ni, T. Yu, Y. H. Lu, Y. Y. Wang, Y. P. Feng, and Z. X. Shen, ACS Nano 2, 2301 (2008), pMID: 19206396.
  * Choi  _et al._ (2010) S.-M. Choi, S.-H. Jhi, and Y.-W. Son, Phys. Rev. B 81, 081407 (2010).
  * Pletikosić  _et al._ (2009) I. Pletikosić, M. Kralj, P. Pervan, R. Brako, J. Coraux, A. T. N’Diaye, C. Busse, and T. Michely, Phys. Rev. Lett. 102, 056808 (2009).
  * Meunier  _et al._ (2016) V. Meunier, A. G. Souza Filho, E. B. Barros, and M. S. Dresselhaus, Rev. Mod. Phys. 88, 025005 (2016).
  * Pereira  _et al._ (2009) V. M. Pereira, A. H. Castro Neto, and N. M. R. Peres, Phys. Rev. B 80, 045401 (2009).
  * Han  _et al._ (2007) M. Y. Han, B. Özyilmaz, Y. Zhang, and P. Kim, Phys. Rev. Lett. 98, 206805 (2007).
  * Hunt  _et al._ (2013) B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, and R. C. Ashoori, Science 340, 1427 (2013).
  * Ponomarenko  _et al._ (2011) L. A. Ponomarenko, A. K. Geim, A. A. Zhukov, R. Jalil, S. V. Morozov, K. S. Novoselov, I. V. Grigorieva, E. H. Hill, V. V. Cheianov, and V. I. Fal’Ko, Nature Physics 7, 958 (2011).
  * Kotov  _et al._ (2012) V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and A. H. Castro Neto, Rev. Mod. Phys. 84, 1067 (2012).
  * Gui  _et al._ (2015) G. Gui, D. Morgan, J. Booske, J. Zhong, and Z. Ma, Applied Physics Letters 106, 666 (2015).
  * Aslani  _et al._ (2015) M. Aslani, C. M. Garner, S. Kumar, D. Nordlund, P. Pianetta, and Y. Nishi, Applied Physics Letters 81, 720 (2015).
  * Fanbanrai  _et al._ (2015) P. Fanbanrai, A. Hutem, and S. Boonchui, Physica B Condensed Matter 472, 84 (2015).
  * Arya  _et al._ (2017) S. Arya, M. S. Laad, and S. R. Hassan, Phys. Rev. B 96, 085126 (2017).
  * Choi  _et al._ (2009) S. M. Choi, S. H. Jhi, and Y. W. Son, Physical Review B Condensed Matter 81, 081407(R) (2009).
  * Cosma  _et al._ (2014) D. A. Cosma, M. Mucha-Kruczyński, H. Schomerus, and V. I. Fal’ko, Phys. Rev. B 90, 245409 (2014).
  * Craco  _et al._ (2015) L. Craco, D. Selli, G. Seifert, and S. Leoni, Physical Review B 91 (2015).
  * Summerfield  _et al._ (2016) A. Summerfield, A. Davies, T. S. Cheng, V. V. Korolkov, J. C. Yong, C. J. Mellor, C. T. Foxon, A. N. Khlobystov, K. Watanabe, and T. Taniguchi, Sci Rep 6, 22440 (2016).
  * Zhang  _et al._ (2016) W. Zhang, W. C. Lu, H. X. Zhang, K. M. Ho, and C. Z. Wang, Carbon 110, 330 (2016).
  * Tang  _et al._ (2015) H.-K. Tang, E. Laksono, J. N. B. Rodrigues, P. Sengupta, F. F. Assaad, and S. Adam, Phys. Rev. Lett. 115, 186602 (2015).
  * Hubbard (1963) J. Hubbard, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 276, 238 (1963).
  * Imada  _et al._ (1998) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
  * Paiva  _et al._ (2005) T. Paiva, R. T. Scalettar, W. Zheng, R. R. P. Singh, and J. Oitmaa, Phys. Rev. B 72, 085123 (2005).
  * Sorella  _et al._ (2012) S. Sorella, Y. Otsuka, and S. Yunoki, Scientific Reports 2, 992 (2012).
  * Assaad and Herbut (2013) F. F. Assaad and I. F. Herbut, Phys. Rev. X 3, 031010 (2013).
  * Arya  _et al._ (2015) S. Arya, P. V. Sriluckshmy, S. R. Hassan, and A.-M. S. Tremblay, Phys. Rev. B 92, 045111 (2015).
  * Trivedi and Randeria (1995) N. Trivedi and M. Randeria, Phys. Rev. Lett. 75, 312 (1995).
  * Trivedi  _et al._ (1996) N. Trivedi, R. T. Scalettar, and M. Randeria, Phys. Rev. B 54, R3756 (1996).
  * Denteneer  _et al._ (1999a) P. J. H. Denteneer, R. T. Scalettar, and N. Trivedi, Phys. Rev. Lett. 83, 4610 (1999a).
  * Ma  _et al._ (2018) T. Ma, L. Zhang, C.-C. Chang, H.-H. Hung, and R. T. Scalettar, Phys. Rev. Lett. 120, 116601 (2018).
  * Denteneer  _et al._ (1999b) P. J. H. Denteneer, R. T. Scalettar, and N. Trivedi, Phys. Rev. Lett. 83, 4610 (1999b).
  * Zhang  _et al._ (1997) S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev. Lett. 78, 4486 (1997).
  * Zhang  _et al._ (1995) S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev. Lett. 74, 3652 (1995).
  * Ma  _et al._ (2011) T. Ma, Z. Huang, F. Hu, and H.-Q. Lin, Phys. Rev. B 84, 121410 (2011).
  * Botello-Méndez  _et al._ (2018) A. R. Botello-Méndez, J. C. Obeso-Jureidini, and G. G. Naumis, The Journal of Physical Chemistry C 122, 15753 (2018).

Figure 6: $\sigma_{dc}$ behavior as the interaction changes at different
temperatures. Panels correspond to different strain strengths from
$\Delta{t}=0.02$ to $\Delta{t}=0.12$. The crosspoint in each panel indicates
the transition point: (a) $\Delta{t}=0.02$,$U\simeq 3.80$, (b)
$\Delta{t}=0.04$,$U\simeq 3.80$,(c) $\Delta{t}=0.06$,$U\simeq 3.75$,(d)
$\Delta{t}=0.08$,$U\simeq 3.70$,(e) $\Delta{t}=0.10$,$U\simeq 3.65$ and (f)
$\Delta{t}=0.12$,$U\simeq 3.40$.

###  .1 Computing the DC conductivity

In this work, the low-temperature behavior of DC conductivity $\sigma_{dc}$ is
used to distinguish metallic or insulating phases. We found that the
transition point changes dramatically at around $U\sim 3.80$. So in appendix
Fig. 6, we carefully checked the phase boundary by plotting the
$\sigma_{dc}(U)$ curves, which indicate a critical transition strength $U_{c}$
in each panel. While in panel (a), when $U<3.80$, $\sigma_{dc}$ gradually
increases as $\beta$ get larger, which indicates the metallic behavior. On the
other side, when $U>3.80$, $\sigma_{dc}$ decreases as $\beta$ increases, which
act as an insulator. Thus, the critical transition point for $\Delta{t}=0.02$
is $U\simeq 3.80$. And other transition points are (b)
$\Delta{t}=0.04$,$U\simeq 3.80$,(c) $\Delta{t}=0.06$,$U\simeq 3.75$,(d)
$\Delta{t}=0.08$,$U\simeq 3.70$,(e) $\Delta{t}=0.10$,$U\simeq 3.65$ and (f)
$\Delta{t}=0.12$,$U\simeq 3.40$.

Figure 7:  Study of the existence of long-range magnetic order in strained
system at various $U$ values.

###  .2 Existence of the AF order phase

To establish the phase diagram, the finite-size effect on the AF spin
structure factor $S_{AF}$ has been carefully examined in the manuscript. We
extrapolated the data to the thermodynamic limit to get the order parameters.
As shown in Fig . 7, we discussed the strain effect on the honeycomb lattice
with various interaction strengths. All the results are summarized in the
phase diagram Fig . 2.


# Concurrent Validity of Automatic Speech and Pause Measures During Passage
Reading in ALS

Saeid Alavi Naeini  Kite Research Institute, Toronto Rehabilitation Institute
– University Health Network  Institute of Biomedical Engineering, University
of Toronto, Toronto, ON, Canada  Leif Simmatis  Kite Research Institute,
Toronto Rehabilitation Institute – University Health Network  Yana Yunusova
Kite Research Institute, Toronto Rehabilitation Institute – University Health
Network  Department of Speech-Language Pathology, University of Toronto,
Toronto, Canada  Babak Taati  Kite Research Institute, Toronto Rehabilitation
Institute – University Health Network  Institute of Biomedical Engineering,
University of Toronto, Toronto, ON, Canada  Department of Computer Science,
University of Toronto, Toronto, ON, Canada  Vector Institute, Toronto, ON,
Canada

###### Abstract

The analysis of speech measures in individuals with amyotrophic lateral
sclerosis (ALS) can provide essential information for early diagnosis and
tracking disease progression. However, current methods for extracting speech
and pause features are manual or semi-automatic, which makes them
time–consuming and labour–intensive. The advent of speech-text alignment
algorithms provides an opportunity for inexpensive, automated, and accurate
analysis of speech measures in individuals with ALS. There is a need to
validate speech and pause features calculated by these algorithms against
current gold standard methods. In this study, we extracted 8 speech/pause
features from 646 audio files of individuals with ALS and healthy controls
performing passage reading. Two pretrained forced alignment models – one using
transformers and another using a Gaussian mixture / hidden Markov architecture
– were used for automatic feature extraction. The results were then validated
against semi-automatic speech/pause analysis software, with further subgroup
analyses based on audio quality and disease severity. Features extracted using
transformer-based forced alignment had the highest agreement with gold
standards, including in terms of audio quality and disease severity. This
study lays the groundwork for future intelligent diagnostic support systems
for clinicians, and for novel methods of tracking disease progression remotely
from home.

###### Index Terms:

Bulbar ALS; concurrent validity; speech; pause; forced alignment

##  I Introduction

TABLE I: Spearman correlation ($\rho$) between features calculated from the Wav2Vec2 method and from the SPA software, divided based on audio quality. The number of audio files in each category is presented in parenthesis | Good (314) | Fair (120) | Poor (212)
---|---|---|---
| $\rho$ | p | $\rho$ | p | $\rho$ | p
Pause duration | $0.82$ | $<.001$ | $0.84$ | $<.001$ | $0.62$ | $<.001$
Total duration | $0.98$ | $<.001$ | $0.96$ | $<.001$ | $0.93$ | $<.001$
Speech duration | $0.96$ | $<.001$ | $0.91$ | $<.001$ | $0.70$ | $<.001$
Pause event | $0.73$ | $<.001$ | $0.80$ | $<.001$ | $0.58$ | $<.001$
$\%$ Pause | $0.72$ | $<.001$ | $0.75$ | $<.001$ | $0.48$ | $<.001$
Mean phrase | $0.77$ | $<.001$ | $0.70$ | $<.001$ | $0.50$ | $<.001$
CV phrase duration | $-0.05$ | $0.39$ | $0.12$ | $0.18$ | $0.01$ | $0.86$
CV pause duration | $0.53$ | $<.001$ | $0.38$ | $<.001$ | $0.27$ | $<.001$
|  |  |  |  |  |
TABLE II: Spearman correlation ($\rho$) between features calculated from the MFA model and from the SPA software, divided based on audio quality. The number of audio files in each category is presented in parenthesis | Good (314) | Fair (120) | Poor (212)
---|---|---|---
| $\rho$ | p | $\rho$ | p | $\rho$ | p
Pause duration | $0.77$ | $<.001$ | $0.74$ | $<.001$ | $0.52$ | $<.001$
Total duration | $0.90$ | $<.001$ | $0.86$ | $<.001$ | $0.85$ | $<.001$
Speech duration | $0.87$ | $<.001$ | $0.86$ | $<.001$ | $0.76$ | $<.001$
Pause event | $0.74$ | $<.001$ | $0.64$ | $<.001$ | $0.56$ | $<.001$
$\%$ Pause | $0.65$ | $<.001$ | $0.62$ | $<.001$ | $0.40$ | $.04$
Mean phrase | $0.60$ | $<.001$ | $0.50$ | $<.001$ | $0.58$ | $<.001$
CV phrase duration | $-0.05$ | $.51$ | $-0.38$ | $.05$ | $0.02$ | $.86$
CV pause duration | $0.32$ | $<.001$ | $0.31$ | $.11$ | $0.27$ | $.11$
|  |  |  |  |  |

Amyotrophic lateral sclerosis (ALS) is a rapidly progressing neurodegenerative
disorder that destroys motor neurons in the cerebral cortex, brainstem, and
spinal cord [1]. The emergence of bulbar signs is an essential milestone in
the progression of ALS, as it causes the loss of speech intelligibility and,
consequently, has a significant impact on the quality of life of people with
ALS. The loss of communication, especially speech, has been reported as one of
the worst aspects of the disease [2].

There are established clinical assessment techniques for diagnosing bulbar
symptoms associated with ALS and tracking disease severity; however, these
have important limitations. Clinical impressions of range of motion, speed,
and symmetry of oral musculature are valuable, but these assessments are
typically subjective [3]. The reliability of these clinician-based judgements
have been questioned and lack of reliable assessment tools contributes to long
diagnostic delay in ALS (up to 18 months) [4]. Earlier detection of bulbar
symptoms in ALS could facilitate finding augmentative and alternative
communication referrals, and assist early bulbar symptom management [5]. There
is a substantial need to develop accessible, automated systems that can
objectively and reliably assess bulbar ALS symptoms in nonclinical settings.

Speech timing can be broadly measured using the Speech Intelligibility Test
(SIT) which provides a measure of speaking rate [6]. Speaking rate has been
recommended as the preferred objective measure of bulbar decline by speech-
language pathologists [7, 8]. Passage reading tasks are excellent alternatives
to the SIT, since they can also provide measures of speech and pause duration
[9]. Studies have shown that these features can distinguish bulbar disease at
different stages, with or without respiratory deficits [10].

Current methods of obtaining speech and pause features are either manual or
semi-automatic, which require time-intensive processing and training of
professional staff [11]. Forced alignment techniques – which temporally align
passage text and speech audio signals – could potentially overcome these
limitations. Past attempts to automatically parse and align dysarthric speech
have been error-prone [12, 13].

Transformers are a new class of deep learning models that have achieved state-
of-the-art performance in a variety of tasks, e.g. in natural language
processing or computer vision [14, 15]. The recent development of robust and
accurate transformer-based forced alignment models [16] presents a new
opportunity for fully automated speech/pause feature extraction in dysarthric
speech. However, it is an open question as to whether novel methods could
replace current gold standards for dysarthric speech alignment.

The aim of this study is to establish the concurrent validity of automatic
speech-text alignment algorithms against the Speech and Pause Analysis (SPA)
software, which is a validated semi-automated pause identification tool [11].
We investigate two open-source forced alignment models – including a new
transformer-based model [16] – and report correlations between 8 speech and
pause features calculated via SPA vs. via forced alignment. Our hypothesis is
that features calculated via automatic text-speech alignment will have strong
relationship with their SPA counterparts.

TABLE III: Spearman correlation ($\rho$) between features calculated from the Wav2Vec2 model on Good audio quality data and from the SPA software, divided based on disease severity. The number of audio files in each category is presented in parenthesis | HC (147) | Mild (70) | Moderate (59) | Severe (38)
---|---|---|---|---
| $\rho$ | p | $\rho$ | p | $\rho$ | p | $\rho$ | p
Pause duration | $0.82$ | $<.001$ | $0.84$ | $<.001$ | $0.66$ | $<.001$ | $0.94$ | $<.001$
Total duration | $0.97$ | $<.001$ | $0.95$ | $<.001$ | $0.96$ | $<.001$ | $0.99$ | $<.001$
Speech duration | $0.93$ | $<.001$ | $0.93$ | $<.001$ | $0.90$ | $<.001$ | $0.98$ | $<.001$
Pause event | $0.53$ | $<.001$ | $0.73$ | $<.001$ | $0.60$ | $<.001$ | $0.92$ | $<.001$
$\%$ Pause | $0.74$ | $<.001$ | $0.82$ | $<.001$ | $0.54$ | $<.001$ | $0.83$ | $.001$
Mean phrase | $0.76$ | $<.001$ | $0.80$ | $<.001$ | $0.68$ | $<.001$ | $0.91$ | $<.001$
CV phrase duration | $0.04$ | $.61$ | $-0.25$ | $.03$ | $0.04$ | $.60$ | $-0.10$ | $.53$
CV pause duration | $0.32$ | $<.001$ | $0.59$ | $<.001$ | $0.67$ | $<.001$ | $0.74$ | $<.001$
|  |  |  |  |  |  |  |
TABLE IV: Spearman correlation ($\rho$) between features calculated from the MFA model on Good audio quality data and from the SPA software, divided based on disease severity. The number of audio files in each category is presented in parenthesis | HC (147) | Mild (70) | Moderate (59) | Severe (38)
---|---|---|---|---
| $\rho$ | p | $\rho$ | p | $\rho$ | p | $\rho$ | p
Pause duration | $0.64$ | $<.001$ | $0.78$ | $<.001$ | $0.51$ | $<.001$ | $0.75$ | $<.001$
Total duration | $0.86$ | $<.001$ | $0.92$ | $<.001$ | $0.85$ | $<.001$ | $0.44$ | $.07$
Speech duration | $0.92$ | $<.001$ | $0.66$ | $.002$ | $0.65$ | $.002$ | $0.63$ | $.007$
Pause event | $0.52$ | $<.001$ | $0.74$ | $<.001$ | $0.77$ | $<.001$ | $0.70$ | $.002$
$\%$ Pause | $0.68$ | $<.001$ | $0.77$ | $<.001$ | $0.41$ | $.07$ | $0.41$ | $.09$
Mean phrase | $0.52$ | $<.001$ | $0.74$ | $<.001$ | $0.70$ | $<.001$ | $0.39$ | $.12$
CV phrase duration | $-0.02$ | $.83$ | $0.13$ | $.60$ | $-0.08$ | $.73$ | $-0.10$ | $.71$
CV pause duration | $0.15$ | $.14$ | $0.56$ | $.01$ | $0.53$ | $.02$ | $0.43$ | $.08$
|  |  |  |  |  |  |  |
Figure 1: The Bland–Altman plots between Wav2Vec2 and SPA features of pause
duration, speech duration, and total duration on Good audio quality data. The
dashed lines are the lower and upper limits of agreement (1.96$\times$standard
deviation) and the solid line is zero.

##  II Methods

###  II-A Participants

We used a dataset of 646 speech recordings, obtained from 3 observational
longitudinal studies of individuals with ALS and healthy controls (HC). The
recordings were of the participants reading a short (60 word) passage – the
bamboo passage [9, 10].

The total number of participants was 462 (ALS: 243, controls: 219). Additional
demographic and clinical details are presented in [11]. The time interval
between recordings varied from 3 to 6 months when more than one recording was
available from a participant. Individuals with ALS were diagnosed with
possible, probable, or definite ALS, as defined by the El Escorial Criteria
[17].

Speaking rate was obtained for each speaker and session using the SIT and was
calculated as the number of words per minute (WPM). Similar to [3], we
classified the individuals with ALS into 3 groups of mild, moderate, and
severe based on their speaking rate. The mild group comprised of individuals
with speaking rates greater than 160 WPM; the moderate group included
individuals with speaking rates between 120 and 160 WPM; and the severe group
had a speaking rate less than 120 WPM.

We also categorized the quality of audio files into 3 groups of Good, Fair,
and Poor based on signal-to-noise ratio (SNR) of each acoustic file. Following
previous research [18], audio files with SNR below 15 dB were labeled as Poor,
15–20 dB were labeled as Fair, and above 20 dB were labeled as Good. The
reason for having different quality groups were due to recordings being
obtained in a clinical (in contrast to a laboratory) context, resulting in
environmental noise and quiet/muffed recordings in some cases.

###  II-B Procedure

We applied 2 different forced alignment techniques for the automatic parsing
and speech/text alignment of the Bamboo passage - a standard passage used for
the analysis of bulbar disease severity in ALS [9]. All recordings were also
previously parsed using the semi-automated SPA software. The first forced
alignment method was the Montreal Forced Aligner (MFA) [19]. The automatic
speech recognition pipeline in MFA is based on a Gaussian mixture/hidden
Markov (GMM/HMM) architecture from the open-source Kaldi speech recognition
toolkit [19]. MFA uses Mel-frequency cepstral coefficients (MFCCs) as acoustic
input features. For the second forced alignment method, we used a pretrained
transformer-based Wav2Vec2 model [16] from the Charsiu library [20]. The model
can perform both text-dependent and text-independent phone-to-audio alignment
and maintains good performance in different settings. In all cases, we
performed speech-text alignment in a supervised manner, i.e. by providing the
text of bamboo passage as an input to the model.

###  II-C Feature Extraction

The forced alignment output comprised the time stamps associated with each
speech and pause event greater than 300 milliseconds. We then extracted the
same 8 features used in [11] based on their ability to distinguish individuals
with ALS from controls. They included:

  1. 1.

Pause duration (sec): total duration of all pause events (excluding the start
and end pauses).

  2. 2.

Total duration (sec): recording time from the onset of the first sentence to
the end of the last sentence.

  3. 3.

Speech duration (sec): total duration of all speech events (pauses excluded).

  4. 4.

Pause events: number of pauses while reading the passage.

  5. 5.

$\%$ Pause: percentage of total reading time spent pausing (only pauses
between speech events).

  6. 6.

Mean phrase (sec): average duration of a phrase. Phrases were defined as
sections of continuous speech without a pause.

  7. 7.

CV phrase duration: Coefficient of variation of phrase durations (a normalized
measure of variability).

  8. 8.

CV pause duration: Coefficient of variation of pause durations.

###  II-D Evaluation

Spearman’s correlation coefficients ($\rho$) were used to determine the
agreement between the features extracted from forced alignment algorithms and
from the SPA software. We compared these sets of 8 features across 3 different
conditions: 1) two forced alignment algorithms (MFA, and Wav2Vec2), 2) four
severity levels (HC, and ALS mild, moderate, or severe), and 3) three audio
qualities (Good, Fair, Poor). The correlation coefficients ($\rho$) and
p-values are reported (significance threshold was p$<$0.05). We also present
the agreement between a subset of the Wav2Vec2 features and their SPA
counterparts using Bland–Altman plots. The subset of 3 features was selected
as those that were previously reported to have the highest area under the
receiver operating characteristic curve (AUC-ROC) in detecting bulbar changes
in otherwise bulbar presymptomatic individuals with ALS [11].

##  III Results

The results of correlation analysis for the transformer-based (Wav2Vec2) and
HMM-based (MFA) models are outlined in Tables I and II, respectively.
Performance is separately reported for recordings based on their audio
quality. Tables III and IV present the correlation results for Wav2Vec2 and
MFA models based on disease severity. This analysis is performed on Good
quality data to ensure disease severity is the only variable affecting the
correlation.

The Bland–Altman plots for pause duration, speech duration, and total duration
are shown in Figure 1. The plots illustrate the agreement between features
extracted using the Wav2Vec2 model vs. using the SPA software. For all three
features, there is similar agreement between Wav2Vec2 and SPA features as all
3 plots have relatively narrow limits of agreement band
($\pm$1.96$\times$standard deviation) which contains 95$\%$ of the values.

##  IV Discussion

The current gold standard for the parsing of clinical audio data and feature
extraction is the manual or semi-automatic approach, which is time–consuming
and labour–intensive. Here, we compared two novel, automated methods of audio-
based parsing to SPA, which has been previously validated for parsing clinical
audio data [9]. We demonstrated that, although both MFA and Wav2Vec2 performed
well with respect to SPA, Wav2Vec2 generalized better across clinical
severities. However, both methods were sensitive to audio data quality. The
present work represents an important step toward developing automated tools
for assessing neurodegenerative disorders, including ALS, in the home setting.

The general pattern than can be observed in the correlation results is that
the Wav2Vec2 model performed better with most features being very strongly (or
strongly) associated. We also observed that correlations between SPA and
automated methods appeared to change as a function of audio sample quality
(Tables I and II), with Good audio having the strongest correlations. This
agrees with the results of previous studies; for example, the performance of
forced alignment algorithms is reduced when sound to noise ratio is not
optimized [19, 20]. Our results reinforce the need for high quality audio data
in the assessment of neurological and neurodegenerative diseases.

The present results add to existing literature and support previous calls to
incorporate more flexible methods into automated and semi-automated speech
analysis. In [21], automated diadochokinetic analysis methods were explored,
and it was found that automated methods’ performance suffered as clinical
severity increased. The results of the present study suggest that this could
be alleviated by using flexible and robust methods, such as Wav2Vec2-based
forced alignment in our current passage-reading paradigm. Importantly, we
observed that correlations between SPA/Wave2Vec2 across clinical severities
(Tables 3 and 4) were almost uniformly better than they were in the SPA/MFA
comparison. This suggests that transformer-based models are more robust to
range of clinical severities and are viable replacements for current manual
and semi-automatic procedures used for clinical audio data analysis.

In conclusion, our work has demonstrated that automated forced alignment
methods perform well compared to a semi-automated ground truth, and that the
Wav2Vec2-based approach in particular is robust across a range of clinical
severity levels. Importantly, our results also highlight the need for good-
quality audio data to be used for clinical speech analyses. These results have
important implications for the development of new clinical tools for home-
based speech assessment, and also for future work that will aim to prognose
patients based on speech function. These are challenging problems [22, 23]
that will require further future work. It will also be important in future
work to extend our current findings to include other populations that
experience speech impairments, such as stroke and Parkinson’s disease.

##  V Acknowledgements

The authors would like to thank the Natural Sciences and Engineering Research
Council of Canada (NSERC), ALS Canada, Parkinson Canada, Canadian Partnership
for Stroke Recovery, Michael J. Fox Foundation and Weston Brain Institute, the
Kite Research Institute – UHN, and Mitacs for supporting our research. We
would also like to thank Dr. Madhura Kulkarni for her role in supporting data
collection and management for this project.

## References

  * [1] M. C. Kiernan, S. Vucic, B. C. Cheah, M. R. Turner, A. Eisen, O. Hardiman, J. R. Burrell, and M. C. Zoing, “Amyotrophic lateral sclerosis,” _The lancet_ , vol. 377, no. 9769, pp. 942–955, 2011.
  * [2] M. Hecht, T. Hillemacher, E. Gräsel, S. Tigges, M. Winterholler, D. Heuss, M.-J. Hilz, and B. Neundörfer, “Subjective experience and coping in als,” _Amyotrophic Lateral Sclerosis and Other Motor Neuron Disorders_ , vol. 3, no. 4, pp. 225–231, 2002.
  * [3] S. Shellikeri, J. R. Green, M. Kulkarni, P. Rong, R. Martino, L. Zinman, and Y. Yunusova, “Speech movement measures as markers of bulbar disease in amyotrophic lateral sclerosis,” _Journal of Speech, Language, and Hearing Research_ , vol. 59, no. 5, pp. 887–899, 2016.
  * [4] M. R. Turner, J. Scaber, J. A. Goodfellow, M. E. Lord, R. Marsden, and K. Talbot, “The diagnostic pathway and prognosis in bulbar-onset amyotrophic lateral sclerosis,” _Journal of the neurological sciences_ , vol. 294, no. 1-2, pp. 81–85, 2010.
  * [5] J. R. Green, Y. Yunusova, M. S. Kuruvilla, J. Wang, G. L. Pattee, L. Synhorst, L. Zinman, and J. D. Berry, “Bulbar and speech motor assessment in als: Challenges and future directions,” _Amyotrophic Lateral Sclerosis and Frontotemporal Degeneration_ , vol. 14, no. 7-8, pp. 494–500, 2013.
  * [6] M. Dorsey, K. Yorkston, D. Beukelman, and M. Hakel, “Speech intelligibility test for windows,” _Institute for Rehabilitation Science and Engineering at Madonna_ , 2007.
  * [7] K. M. Yorkston, “Speech deterioration in amyotrophic lateral sclerosis: Implications for the timing of intervention,” _Jounal of Medical Speech-Language Pathology_ , vol. 1, pp. 35–46, 1993.
  * [8] L. J. Ball, A. Willis, D. R. Beukelman, and G. L. Pattee, “A protocol for identification of early bulbar signs in amyotrophic lateral sclerosis,” _Journal of the neurological sciences_ , vol. 191, no. 1-2, pp. 43–53, 2001\.
  * [9] J. R. Green, D. R. Beukelman, and L. J. Ball, “Algorithmic estimation of pauses in extended speech samples of dysarthric and typical speech,” _Journal of medical speech-language pathology_ , vol. 12, no. 4, p. 149, 2004\.
  * [10] Y. Yunusova, N. L. Graham, S. Shellikeri, K. Phuong, M. Kulkarni, E. Rochon, D. F. Tang-Wai, T. W. Chow, S. E. Black, L. H. Zinman _et al._ , “Profiling speech and pausing in amyotrophic lateral sclerosis (als) and frontotemporal dementia (ftd),” _PloS one_ , vol. 11, no. 1, p. e0147573, 2016.
  * [11] C. Barnett, J. R. Green, R. Marzouqah, K. L. Stipancic, J. D. Berry, L. Korngut, A. Genge, C. Shoesmith, H. Briemberg, A. Abrahao _et al._ , “Reliability and validity of speech & pause measures during passage reading in als,” _Amyotrophic Lateral Sclerosis and Frontotemporal Degeneration_ , vol. 21, no. 1-2, pp. 42–50, 2020.
  * [12] M. Pleva, J. Juhár, and A. S. Thiessen, “Automatic acoustic speech segmentation in praat using cloud based asr,” in _2015 25th International Conference Radioelektronika (RADIOELEKTRONIKA)_. IEEE, 2015, pp. 172–175.
  * [13] Y. T. Yeung, K. H. Wong, and H. Meng, “Improving automatic forced alignment for dysarthric speech transcription,” in _Sixteenth Annual Conference of the International Speech Communication Association_ , 2015.
  * [14] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [15] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly _et al._ , “An image is worth 16x16 words: Transformers for image recognition at scale,” _arXiv preprint arXiv:2010.11929_ , 2020.
  * [16] A. Baevski, Y. Zhou, A. Mohamed, and M. Auli, “wav2vec 2.0: A framework for self-supervised learning of speech representations,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 12 449–12 460, 2020.
  * [17] B. R. Brooks, R. G. Miller, M. Swash, and T. L. Munsat, “El escorial revisited: revised criteria for the diagnosis of amyotrophic lateral sclerosis,” _Amyotrophic lateral sclerosis and other motor neuron disorders_ , vol. 1, no. 5, pp. 293–299, 2000.
  * [18] H. Xu and H. Zheng, “The simple SNR estimation algorithms for mpsk signals,” in _Proceedings 7th International Conference on Signal Processing, 2004. Proceedings. ICSP’04. 2004._ , vol. 2. IEEE, 2004, pp. 1781–1785.
  * [19] M. McAuliffe, M. Socolof, S. Mihuc, M. Wagner, and M. Sonderegger, “Montreal forced aligner: Trainable text-speech alignment using kaldi.” in _Interspeech_ , vol. 2017, 2017, pp. 498–502.
  * [20] J. Zhu, C. Zhang, and D. Jurgens, “Phone-to-audio alignment without text: A semi-supervised approach,” _arXiv preprint arXiv:2110.03876_ , 2021.
  * [21] C. Tanchip, D. L. Guarin, S. McKinlay, C. Barnett, S. Kalra, A. Genge, L. Korngut, J. R. Green, J. Berry, L. Zinman _et al._ , “Validating automatic diadochokinesis analysis methods across dysarthria severity and syllable task in amyotrophic lateral sclerosis,” _Journal of Speech, Language, and Hearing Research_ , vol. 65, no. 3, pp. 940–953, 2022.
  * [22] F. Kimura, C. Fujimura, S. Ishida, H. Nakajima, D. Furutama, H. Uehara, K. Shinoda, M. Sugino, and T. Hanafusa, “Progression rate of alsfrs-r at time of diagnosis predicts survival time in als,” _Neurology_ , vol. 66, no. 2, pp. 265–267, 2006.
  * [23] J. R.-J. Lee, J. F. Annegers, and S. H. Appel, “Prognosis of amyotrophic lateral sclerosis and the effect of referral selection,” _Journal of the neurological sciences_ , vol. 132, no. 2, pp. 207–215, 1995.