# Localization effects due to a random magnetic field on heat transport in a
harmonic chain

Gaëtan Cane1, Junaid Majeed Bhat2, Abhishek Dhar3 and Cédric Bernardin4 1
Université Côte d’Azur, CNRS, LJAD, Parc Valrose, 06108 NICE Cedex 02, France
<EMAIL_ADDRESS>2 International Centre for Theoretical Sciences,
Bengaluru, India-560089<EMAIL_ADDRESS>3 International Centre for
Theoretical Sciences, Bengaluru, India-560089<EMAIL_ADDRESS>4
Université Côte d’Azur, CNRS, LJAD, Parc Valrose, 06108 NICE Cedex 02, France
Interdisciplinary Scientific Center Poncelet (CNRS IRL 2615), 119002 Moscow,
Russia<EMAIL_ADDRESS>
###### Abstract

We consider a harmonic chain of $N$ oscillators in the presence of a
disordered magnetic field. The ends of the chain are connected to heat baths
and we study the effects of the magnetic field randomness on heat transport.
The disorder, in general, causes localization of the normal modes due to which
a system becomes insulating. However, for this system, the localization length
diverges as the normal mode frequency approaches zero. Therefore, the low
frequency modes contribute to the transmission, $\mathcal{T}_{N}(\omega)$, and
the heat current goes down as a power law with the system size, $N$. This
power law is determined by the small frequency behaviour of some Lyapunov
exponent, $\lambda(\omega)$, and the transmission in the thermodynamic limit,
$\mathcal{T}_{\infty}(\omega)$. While it is known that in the presence of a
constant magnetic field
$\mathcal{T}_{\infty}(\omega)\sim\omega^{3/2},\leavevmode\nobreak\
\omega^{1/2}$ depending on the boundary conditions, we find that the Lyapunov
exponent for the system behaves as $\lambda(\omega)\sim\omega$ for
$\expectationvalue{B}\neq 0$ and $\lambda(\omega)\sim\omega^{2/3}$ for
$\expectationvalue{B}=0$. Therefore, we obtain different power laws for
current vs $N$ depending on $\expectationvalue{B}$ and the boundary
conditions.

Keywords: Heat conduction, Transport properties.

##  1 Introduction

In his seminal paper [1] Anderson studied the conductance of electrons and
explained how the presence of impurities in the metal could reduce drastically
the diffusive motion of the electrons up to a complete halt and thus giving
place to an insulator. This phenomenon depends strongly on the dimension and
the Metal Insulator Transition is not instantaneous with respect to the
disorder strength only in dimension greater or equal to three. Nowadays
Anderson localization is seen as a generic phenomenon present in disordered
media whereby the addition of random defects in the medium has the tendency to
localize in space the normal modes of the system. As a consequence it reduces
drastically the transport coefficient. Originally Anderson’s work took place
in a quantum context but the phenomenon he explained appears also in a
classical one. In the 70’s Lebowitz and others [2, 3, 4, 5] started to
investigate the effect of impurities (random masses) on the transport
properties of a one dimensional harmonic chain, arguing in particular that the
conductivity $\kappa(N)$ of the chain (which is proportional to the system
size $N$ for a purely harmonic chain), loses some order of magnitude because
of disorder: $\kappa(N)\sim N^{a}$, $a<1$ (see [6] for a mathematical proof).
But at the difference with respect to original Anderson localization, the
conductivity does not become exponentially small in the system size and,
depending on the physical boundary conditions and thermostats, it can vanish
($a<0$), diverge ($0<a<1$) or even converge ($a=0$) [7, 8, 9]. The reason for
this is roughly due to the fact that for disordered harmonic chains, normal
modes with frequency $\omega$ becomes localized but with a length of
localization $\ell(\omega)\sim\omega^{-2}$. The role of thermostats and
boundary conditions is more difficult to explain without going into
computational details. Thus we note that in disordered harmonic chains the low
frequency modes have still the possibility to transport energy [10]. The case
where the disorder is in the interparticle springs instead of the masses has
recently been addressed in [11, 12]. In higher dimensions the situation is
less understood [13, 14]. More recently there has been a renewed interest for
these questions with respect to the effect of nonlinearities [15, 16] or of an
energy conserving noise [17, 18, 19, 20].

In this paper we consider an ordered (constant masses) one-dimensional chain
of two-dimensional charged oscillators subject to a random transverse magnetic
field on every lattice site (or equivalently a chain of disorderly charged
oscillators subjected to a constant magnetic field on the lattice). In [21],
by using the Non-Equlibrium Greens’s Function (NEGF) formalism we obtained an
explicit expression of the heat current in the steady state and investigated
then the transport properties of such a system when all the charges are the
same. We established that transport is ballistic like for ordered harmonic
chains. The aim of this paper is therefore to describe the effect of the
charge impurities on the behaviour of the conductivity of the system. This
will require us to investigate the frequency dependence of the localization
lengths of normal modes. We show that due to charge disorder the current shows
different scaling with the system size which depends on the boundary
conditions as well as on the expectation value of the magnetic field.

This paper is structured as follows: in Sec. 2 we introduce the model and
state the results for heat current using the NEGF formalism. We also present
numerical results for the transmission function and discuss the effects of
localization due to the random magnetic field on the transmission function. In
Sec. 3, we use the Green’s function expression as a product of random matrices
to determine the Lyapunov exponents. We also present numerical results for
Lyapunov exponents which are consistent with our theoretical results. Using
the results for the Lyapunov exponents, we finally determine the size
dependence of the mean of the heat current in Sec. 4 and compare with direct
numerical calculations of the current. We conclude in Sec. 5.

##  2 The Model and heat current by NEGF

###  2.1 The model

We consider a chain of $N$ harmonic oscillators each having two transverse
degree’s of freedom so that every oscillator is free to move in a plane
perpendicular to the length of the chain. We choose the plane of motion to be
the $x-y$ plane and denote the positions and momenta of the $n^{th}$
oscillator by $(x_{n},y_{n})$ and $(p_{n}^{x},p_{n}^{y})$ respectively, with
$n=1,2,\ldots,N$. The oscillators are assumed to have unit masses and each
carry a positive unit charge. We consider a magnetic field
$\vec{B}_{n}=B_{n}{\vec{\mathbf{e}}}_{z}$ perpendicular to the plane of motion
which can be obtained from a vector potential $\vec{\bf
A}_{n}=(-B_{n}y_{n},B_{n}x_{n},0)$ at each lattice site. In this paper, we
assume that $(B_{n})_{n}$ form a sequence of independent identically
distributed random variables with average $\langle B\rangle$ and variance
$\sigma^{2}$. The Hamiltonian of the chain is given by:

$\displaystyle
H=\sum_{n=1}^{N}\frac{(p_{n}^{x}+B_{n}y_{n})^{2}+(p_{n}^{y}-B_{n}x_{n})^{2}}{2}$
$\displaystyle+\sum_{n=0}^{N}\frac{(x_{n+1}-x_{n})^{2}+(y_{n+1}-y_{n})^{2}}{2},$

where the inter particle spring constant has been fixed to $1$. We will
consider the two different boundary conditions: (i) fixed boundaries with
$x_{0}=x_{N+1}=0$ and (ii) free boundaries with
$x_{0}=x_{1},\leavevmode\nobreak\ x_{N}=x_{N+1}$. In order to study heat
current through this system, we consider the $1^{\rm st}$ and the $N^{\rm th}$
oscillators to be connected to heat reservoirs at temperatures $T_{L}$ and
$T_{R}$ respectively. The heat reservoirs are modelled using dissipative and
noise terms leading to the following Langevin equations of motion:

$\displaystyle\ddot{x}_{n}$
$\displaystyle=(x_{n+1}+x_{n-1}-c_{n}x_{n})+B_{n}\dot{y}_{n}+\eta_{L}^{x}(t)\delta_{n,1}+\eta_{R}^{x}(t)\delta_{n,N}-(\gamma\delta_{n,1}+\gamma\delta_{n,N})\dot{x}_{n}\
,$ (1) $\displaystyle\ddot{y}_{n}$
$\displaystyle=(y_{n+1}+y_{n-1}-c_{n}y_{n})-B_{n}\dot{x}_{n}+\eta_{L}^{y}(t)\delta_{n,1}+\eta_{R}^{y}(t)\delta_{n,N}-(\gamma\delta_{n,1}+\gamma\delta_{n,N})\dot{y}_{n}\
.$ (2)

for $n=1,2,\ldots,N$. Here $\eta_{L}(t):=(\eta_{L}^{x}(t),\eta_{L}^{y}(t))$
and $\eta_{R}(t):=(\eta_{R}^{x}(t),\eta_{R}^{y}(t))$ are Gaussian white noise
terms acting on the $1^{\rm st}$ and $N^{\rm th}$ oscillators respectively.
These follow the regular white noise correlations,
$\expectationvalue{\eta_{L/R}(t)\eta_{L/R}(t^{\prime})}=\sqrt{2\gamma
T_{L/R}}\delta(t-t^{\prime})$ (Boltzmann’s constant is fixed to one to
simplify), where $\gamma$ is the dissipation strength at the reservoirs. The
coefficients $c_{n}$ fix the boundary conditions of the problem. For fixed
boundaries $c_{n}=2$ for all $n$, while for free boundary conditions
$c_{n}=2-\delta_{n,1}-\delta_{n,N}$.

###  2.2 Heat current

In [21], by using the non-equilibrium Green function formalism, we obtained an
exact expression for the heat current, ${J}_{N}$, in the steady state of the
chain. More exactly, let us define the processes $(f_{n}^{\pm})_{n\geq 0}$ and
$(g_{n}^{\pm})_{n\geq 0}$ as

$\begin{split}f_{n+1}^{\pm}&=(c_{n+1}-\omega^{2}\pm\omega
B_{n+1})f_{n}^{\pm}-f_{n-1}^{\pm},\quad f_{0}^{\pm}=1,\quad
f_{1}^{\pm}=c_{1}-\omega^{2}\pm\omega B_{1}\ ,\\\
g_{n+1}^{\pm}&=(c_{n+1}-\omega^{2}\pm\omega
B_{n+1})g_{n}^{\pm}-g_{n-1}^{\pm},\quad g_{0}^{\pm}=0,\quad g_{1}^{\pm}=1\
.\end{split}$ (3)

Then we introduce

$F_{N}^{\pm}:=F_{N}^{\pm}(\omega)=f_{N}^{\pm}+i{\gamma}\omega(g_{N}^{\pm}+f_{N-1}^{\pm})-{\gamma^{2}}\omega^{2}g_{N-1}^{\pm}\
,$ (4)

and the heat current is equal to

${J}_{N}=(T_{L}-T_{R})\,\int_{-\infty}^{\infty}d\omega\,{\mathcal{T}}_{N}(\omega)=2(T_{L}-T_{R})\,\int_{0}^{\infty}d\omega\,{\mathcal{T}}_{N}(\omega)$
(5)

with the net transmission function ${\mathcal{T}}_{N}$ defined for any
frequency $\omega$ by

${\mathcal{T}}_{N}(\omega):=\frac{\gamma^{2}}{\pi}\ \omega^{2}\
\left[\frac{1}{\absolutevalue{F_{N}^{+}(\omega)}^{2}}+\frac{1}{\absolutevalue{F_{N}^{-}(\omega)}^{2}}\right].$
(6)

We denote by $\langle{J}_{N}\rangle$ the expectation of the heat current with
respect to the magnetic field distribution $\langle\cdot\rangle$ and our goal
is to understand its scaling behavior in $N$.

Observe that the stochastic processes $(f_{n}^{-})_{n\geq 0}$ and
$(g_{n}^{-})_{n\geq 0}$ are defined in terms of the two dimensional discrete
time Markov chain $(U_{n})_{n\geq 0}$ given by

$U_{n+1}=\begin{pmatrix}2-\omega^{2}-\omega B_{n+1}&-1\\\ 1&0&\end{pmatrix}\
U_{n},\quad\text{where}\quad U_{n}:=\begin{pmatrix}u_{n}\\\
u_{n-1}\end{pmatrix},$ (7)

by choosing suitable initial conditions. By replacing the $B_{n}$’s by
$-B_{n}$’s in the last display, we see that $(f_{n}^{+})_{n\geq 0}$ and
$(g_{n}^{+})_{n\geq 0}$ can also be expressed in terms of $(U_{n})_{n\geq 0}$
. The state of the Markov chain is nothing but the result of a product of
$2\crossproduct 2$ product of independent and identically distributed random
matrices. Roughly, the behaviour of $F_{N}^{\pm}$ is related to the growth of
$\|U_{N}(\omega)\|$ which will be in the form $e^{2\lambda(\omega)N}$, where

$\lambda(\omega)=\lim_{n\to\infty}\cfrac{1}{2n}\left\langle\log\|U_{n}(\omega)\|\right\rangle=\lim_{n\to\infty}\cfrac{1}{n}\left\langle\log|u_{n}(\omega)|\right\rangle>0,$
(8)

with $\left\langle...\right\rangle$ denoting a disorder average, is half the
Lyapunov exponent associated to the Markov chain $(U_{n})_{n\geq 0}$, or
equivalently of the corresponding product of random matrices. The limit exists
by Furstenberg’s Theorem [22], is non-negative, independent of the initial
condition $U_{0}$ and the limit holds in fact also for any realisation and not
only by averaging over the magnetic field distribution.

For now we quickly discuss the effect of localization due to the random
magnetic field on the heat transport and the need for calculating the Lyapunov
exponent $\lambda(\omega)$ for small frequencies $\omega$.

###  2.3 Effect of localization due to random magnetic field on the net
Transmission

(a) Uniform magnetic field

(b) Random magnetic field

Figure 1: Variation of the net transmission, in units of $k_{B}=1$, with
$\omega$ for uniform magnetic field, panel (a), and random magnetic field,
panel (b). The axes are in log scale and $\gamma=0.2$. The magnetic field in
(a) is set to be $1$ on all oscillators and in (b) it was chosen uniformly
from the interval $(0,2)$. As can be seen clearly from the plots, the
localization effects cause suppression of the transmission.

Using Eq. (6), we can calculate the net transmission $\mathcal{T}_{N}(\omega)$
for any spatial configuration of the magnetic field using a computer
programme. In Fig. (1a) and Fig. (1b) we plot the net transmission function
with $\omega$ for a uniform magnetic field and for a random magnetic field for
different system sizes respectively. On comparison of the two plots, we can
see that the randomness causes suppression of the net transmission and also
the net transmission for the random magnetic field case goes down with system
size while the system size has nearly no effect on the transmission for the
uniform magnetic field. The suppression in case of random magnetic field is
due to localization of the normal modes of the system. The normal modes of
frequency $\omega$ get exponentially localized due to randomness with a
localization length given by $1/\lambda(\omega)$ where $\lambda(\omega)$ is
the Lyapunov exponent defined in Eq. (8). As a result of this they a priori do
not contribute to the transmission. However, note that the transmission for
random magnetic field is higher near $\omega=0$ and goes down as we move away
which means that the normal modes with energies closer to $\omega=0$ have a
larger localization length, i.e. $\lambda(\omega)\to 0$ as $\omega\to 0$.
Since we are eventually interested in the size dependence of the current, for
large $N$, which is the integral of the transmission over all $\omega$, we can
reduce the integration limit to values of $\omega$ for which the localization
length is greater than the system size. For the remaining $\omega$ values for
which the localization length is less than the system size, the transmission
would be negligible. Hence, we cut off the integral limit to
$\omega=\omega^{N}_{max}$ where $1/\lambda(\omega^{N}_{max})=N$ and the
current is then given by

$\left\langle{J}_{N}\right\rangle\approx
2(T_{L}-T_{R})\int_{0}^{\omega^{N}_{max}}d\omega\lim\limits_{N\rightarrow\infty}\left\langle\mathcal{T}_{N}(\omega)\right\rangle=2(T_{L}-T_{R})\int_{0}^{\omega^{N}_{max}}d\omega\mathcal{T}_{\infty}(\omega)\
.$ (9)

Note that the frequency $\omega^{N}_{max}$ would be very small for large $N$,
and for such small frequencies we expect $\mathcal{T}_{\infty}(\omega)$ to
have a weak dependence on disorder [since in the recursion Eq. (7), the
randomness is multiplied by $\omega$] — hence in the above equation
$\mathcal{T}_{\infty}(\omega)$ is written without a disorder average and can
in fact be determined by considering the chain in a _constant_ magnetic field
of strength $\expectationvalue{B}$. In [21], we proved that for constant
magnetic field $\expectationvalue{B}\neq 0$,
$\mathcal{T}_{\infty}(\omega)\sim\omega^{3/2}$ and $\sim\omega^{1/2}$ for
fixed and free boundaries respectively, while for $\expectationvalue{B}=0$ it
goes as $\omega^{2}$ and $\omega^{0}$ for the two boundary conditions
respectively. To determine the size dependence of the current in addition to
the small $\omega$ behaviour of $\mathcal{T}_{\infty}(\omega)$ we also need
the small $\omega$ behaviour of $\lambda(\omega)$. We now proceed to the next
section where we discuss the Lyapunov exponents of this equation.

##  3 Analysis of the Lyapunov exponents

In this section we present theoretical and numerical results on the
asymptotics of Lyapunov exponents for small $\omega$ for the Markov processes
defined by Eq. (7). The Lyapunov exponents are independent of the boundary
conditions – so for this section we only work with fixed boundary conditions
by setting $c_{n}=2$ for all $n$ – and of the initial condition of the process
– i.e. it is the same for $f_{n}^{\pm}$ and $g_{n}^{\pm}$. We show that Eq.
(7) has three different behaviors for the Lyapunov exponent depending on the
expected value $\langle B\rangle$ of the random magnetic field. For $\langle
B\rangle>0$ the Lyapunov exponent satisfies $\lambda(\omega)\sim\omega$ and
for $\langle B\rangle<0$, $\lambda(\omega)\sim\omega^{1/2}$. However, for
$\langle B\rangle=0$, $\lambda(\omega)\sim\omega^{2/3}$. Similar Lyapunov
exponent behaviours are found for a harmonic oscillator with parametric noise,
[23] and we will see that Eq. (7) could be written exactly in this form in the
continuum limit.

###  3.1 Theoretical results for Lyapunov exponents

Let $(z_{t})_{t\geq 0}\in{\mathbb{R}}^{2}$ be the solution of the following
stochastic differential equation (with arbitrary initial condition)

$\dot{z_{t}}=A_{0}z_{t}+\varepsilon\sigma\xi_{t}\,A_{1}z_{t}\ ,$ (10)

where $\varepsilon$ is a small positive parameter, $\sigma>0$ a constant,
$\xi_{t}$ a one dimensional standard white noise and $A_{0}$ and $A_{1}$ are
$2\times 2$ matrices such that

$A_{0}=\begin{pmatrix}0&1\\\ -c&0\end{pmatrix},\quad
A_{1}=\begin{pmatrix}0&0\\\ -1&0\end{pmatrix}\quad\text{with
$c\in\mathbb{R}$}\ .$

The Lyapunov exponent $\lambda_{z}(\varepsilon)$ of the process
$(z_{t})_{t\geq 0}$ is defined by

$\lambda_{z}(\varepsilon)=\lim_{t\rightarrow\infty}\frac{1}{t}\left\langle\log\|z_{t}\|\right\rangle,$
(11)

where $\langle\cdot\rangle$ denotes the expectation with respect to the white
noise. It is proved in A that if we denote $z_{t}=(u_{t},v_{t})^{\bot}$ then
we have the Lyapounov exponent for $(u_{t})_{t\geq 0}$ is the same as for
$(z_{t})_{t\geq 0}$:

$\lambda_{z}(\varepsilon)=\lim_{t\rightarrow\infty}\frac{1}{t}\left\langle\log|u_{t}|\right\rangle.$
(12)

The following result, proved in [24], gives the behaviour of the Lyapunov
exponent $\lambda_{z}(\varepsilon)$ for small noise

  1. 1.

If $c=0$ then
$\lambda_{z}(\varepsilon)=\hat{\lambda}(\sigma)\varepsilon^{2/3}$ where
$\hat{\lambda}(\sigma)$ is defined in Eq. (27) .

  2. 2.

If $c>0$ then
$\lambda_{z}(\varepsilon)\sim\frac{\sigma^{2}}{8c}\varepsilon^{2}$ .

  3. 3.

If $c<0$ then $\lambda_{z}(\varepsilon)\sim\sqrt{|c|}$ .

A sketch of the proof of this result is given in A.

Consider now Eq. (7) defining the discrete time Markov chain
$U_{n}=(u_{n},u_{n-1})^{\top}$ and rewrite it in the following form, for small
$\omega$,

$u_{n+1}+u_{n-1}-2u_{n}=-\omega\langle B\rangle u_{n}-\omega(B_{n+1}-\langle
B\rangle)u_{n}+{\mathcal{O}}(\omega^{2})\ .$

In the continuum limit, the discrete time process $(u_{n})_{n\geq 0}$ becomes
then the continuous time process $(u_{t})_{t\geq 0}$ solution of

$\ddot{u}_{t}=-\omega\langle B\rangle u_{t}-\omega\sigma\xi_{t}u_{t}$ (13)

where $(\xi_{t})_{t\geq 0}$ is a standard white noise and $\sigma^{2}$ the
variance of the $(B_{n})_{n}$. Defining $w_{t}=(u_{t},{\dot{u}}_{t})^{\top}$
we see that the previous equation reads

$\dot{w}_{t}=\begin{pmatrix}0&1\\\ -\omega\langle
B\rangle&0\end{pmatrix}w_{t}+\sigma\omega\xi_{t}\begin{pmatrix}0&0\\\
-1&0\end{pmatrix}w_{t}\ .$ (14)

We are interested in the Lyapunov exponent of the process $(u_{t})_{t\geq 0}$
(or equivalently of the process $(w_{t})_{t\geq 0}$ as said before):

$\lambda_{w}(\omega)=\lim_{t\to\infty}\cfrac{1}{t}\left\langle\log\|w_{t}\|\right\rangle=\lim_{t\to\infty}\cfrac{1}{t}\left\langle\log|u_{t}|\right\rangle$
(15)

Eq. (14) looks similar to Eq. (10) but to fit perfectly with Eq. (10) we
perform the time scaling

${\tilde{u}}_{t}=u_{t/\sqrt{\omega}}$

in Eq. (13) wich gives by scaling invariance of white noise

$\ddot{\tilde{u}}_{t}=-\langle
B\rangle{\tilde{u}}_{t}-\omega^{1/4}\sigma\xi_{t}{\tilde{u}}_{t}$ (16)

or equivalently for
${\tilde{z}}_{t}=({\tilde{u}}_{t},\dot{\tilde{u}}_{t})^{\top}$ the equation

$\dot{\tilde{z}}_{t}=\begin{pmatrix}0&1\\\ -\langle
B\rangle&0\end{pmatrix}{\tilde{z}}_{t}+\sigma\omega^{1/4}\xi_{t}\begin{pmatrix}0&0\\\
-1&0\end{pmatrix}{\tilde{z}}_{t}\ .$ (17)

With the previous notation we have hence

$\lambda_{w}(\omega)=\sqrt{\omega}\ {\lambda}_{\tilde{z}}(\omega^{1/4})\ .$
(18)

Eq. (17) fits perfectly Eq. (10) with $c=\langle B\rangle$ and
$\varepsilon=\omega^{1/4}$. Then using point (i), (ii) and (iii) of Eq. (10)
and Eq. (18) we get

  1. 1.

If $\langle B\rangle=0$,
$\lambda_{w}(\omega)=\hat{\lambda}(\sigma)\omega^{2/3}$ where
$\hat{\lambda}(\sigma)$ is defined in Eq. (27) .

  2. 2.

If $\langle B\rangle>0$,
$\lambda_{w}(\omega)\sim\frac{\sigma^{2}}{8\expectationvalue{B}}\omega$ .

  3. 3.

If $\langle B\rangle<0$,
$\lambda_{w}(\omega)\sim\sqrt{|\expectationvalue{B}|}\omega^{1/2}$ .

It makes sense to believe that $\lambda(\omega)$ defined by Eq. (8) and
$\lambda_{w}(\omega)$ defined by Eq. (15) have roughly the same behaviour as
$\omega\to 0$ but a strong theoretical argument supporting this belief is
missing. However, in the case $\langle B\rangle>0$, we can obtain directly the
behaviour of $\lambda(\omega)$ by following the approach of [25] and we
observe then a good agreement at first order between $\lambda(\omega)$ and
$\lambda_{w}(\omega)$, not only at the level of the exponent in $\omega$ but
also at the level of the prefactor, see Table 3.2. Unfortunately we were not
able to carry this approach for $\langle B\rangle<0$ or $\langle B\rangle=0$
and we decided hence to not pursue this approach. However numerical results
presented in the next section support strongly the claim that
$\lambda(\omega)\sim\lambda_{w}(\omega)$ for $\omega\to 0$.

###  3.2 Numerical results for Lyapunov exponents

(a) $\expectationvalue{B}>0$

(b) $\expectationvalue{B}=0$

(c) $\expectationvalue{B}<0$

Figure 2: Variation of numerically calculated Lyapunov exponent,
$\lambda=\frac{1}{N}\langle\log|u_{N}|\rangle$, with $\omega$.
$\langle\log|u_{N}|\rangle$ denotes average of $\log|u_{N}|$ over 100
realizations of the random magnetic field. For (a), (b) and (c) the magnetic
fields were chosen randomly from the intervals $(0,1)$, $(-1,1)$ and $(-1,0)$
respectively. The solid line is the data from the simulation while the dashed
line is a power law fit, $C\omega^{s}$, to the data with $C$ and $s$ as
fitting parameters. The obtained values of the fitting parameters agree
appreciably with the theoretical values.

Case | Range of $B_{n}$ | $s:\lambda(\omega)\sim C\omega^{s}$ | $C$ | $C_{\rm{theoretical}}$
---|---|---|---|---
$\leavevmode\nobreak\ \leavevmode\nobreak\ \expectationvalue{B}>0$ $\lambda(\omega)\sim\frac{\sigma^{2}}{8\expectationvalue{B}}\omega$ | $(0,0.25)$ | 0.986 | 0.0045 | 0.0052
$(0,0.5)$ | 0.999 | 0.0102 | 0.0104
$(0,0.75)$ | 1.0005 | 0.0156 | 0.0156
$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \expectationvalue{B}<0$ $\lambda(\omega)\sim\sqrt{|\langle B\rangle|}\omega^{1/2}$ | $(-0.25,0)$ | 0.492 | 0.315 | 0.353
$(-0.5,0)$ | 0.492 | 0.444 | 0.5
$(-0.75,0)$ | 0.491 | 0.532 | 0.612
$\leavevmode\nobreak\ \expectationvalue{B}=0$ $\lambda(\omega)=\hat{\lambda}(\sigma)\omega^{2/3}$ | $(-0.25,0.25)$ | 0.658 | 0.073 | 0.079
$(-0.5,0.5)$ | 0.658 | 0.115 | 0.127
$(-0.75,0.75)$ | 0.649 | 0.136 | 0.167
Table 1: Comparison of analytical prefactor for the three cases with the
numerical prefactor. For this table, $N=10^{7}$.

We numerically calculate the Lyapunov exponents by using Eq. (3) to generate
$u_{N}$ for 100 realizations of the random magnetic field. The Lyapunov
exponent would then be given by
$\lambda=\frac{1}{N}\langle\log|u_{N}|\rangle$, where $N$ is the number of
oscillators. We plot in Fig. (2), the numerical data thus obtained for
different $\omega$ and the power law fit, $C\omega^{s}$, for the data with $C$
and $s$ as fitting parameters. We see that the values of $s$ obtained for the
three casses, $\langle B\rangle>0,\leavevmode\nobreak\ \langle B\rangle<0$ and
$\expectationvalue{B}=0$, agree reasonably well with the theoretically
expected values. The prefactor, $C$, obtained for the three cases also seems
to agree with the expected values from theory, see Table 3.2.

We now have the behaviour of the Lyapunov exponents at small $\omega$ for Eq.
(3) and we found this to be different depending on the expectation value of
the random magnetic field. The transmission is determined by $f_{N}^{+}$ as
well as $f_{N}^{-}$ and these two have different Lyapunov exponents for
$\expectationvalue{B}\neq 0$, therefore the larger of the two exponents will
dominate in the transmission. This is the Lyapunov exponent for $f_{N}^{-}$
for $\expectationvalue{B}>0$, while for $\expectationvalue{B}=0$, $f_{N}^{+}$
and $f_{N}^{-}$ have the same Lyapunov exponent. In the next section, we
determine the size dependence of the current using these results for the
Lyapunov exponents.

##  4 Size dependence of the current

We now have the small $\omega$ behaviour of $\lambda(\omega)$ for the
transmission. We found this to be different for $\expectationvalue{B}\neq 0$
and $\expectationvalue{B}=0$, so we expect different power laws for the
current for the two cases. The boundary conditions will also play a role in
the power law via the small $\omega$ behaviour of
$\mathcal{T}_{\infty}(\omega)$. We therefore take the cases
$\expectationvalue{B}\neq 0$ and $\expectationvalue{B}=0$ separately for the
two boundary conditions.

#### Fixed boundary conditions:

  * (a)

For $\expectationvalue{B}\neq 0$,
$\mathcal{T}_{\infty}(\omega)\sim\omega^{3/2}$ and
$\lambda(\omega)\sim\omega$. Therefore using these in Eq. (9) we have
$\langle{J}_{N}\rangle\sim 1/N^{5/2}$.

  * (b)

For $\expectationvalue{B}=0$, $\mathcal{T}_{\infty}(\omega)\sim\omega^{2}$ and
$\lambda(\omega)\sim\omega^{2/3}$ which gives $\langle{J}_{N}\rangle\sim
1/N^{9/2}$.

Boundary Conditions | Average magnetic field | $\mathcal{T}_{\infty}(\omega)$ | $\lambda(\omega)$ | Power law for the current $\langle{J}_{N}\rangle$
---|---|---|---|---
Fixed | $\expectationvalue{B}\neq 0$ | $\sim\omega^{3/2}$ | $\sim\omega$ | $\sim 1/N^{5/2}$
Fixed | $\expectationvalue{B}=0$ | $\sim\omega^{2}$ | $\sim\omega^{2/3}$ | $\sim 1/N^{9/2}$
Free | $\expectationvalue{B}\neq 0$ | $\sim\omega^{1/2}$ | $\sim\omega$ | $\sim 1/N^{3/2}$
Free | $\expectationvalue{B}=0$ | $\sim\omega^{0}$ | $\sim\omega^{2/3}$ | $\sim 1/N^{3/2}$

Table 2: Power law for the current for different boundary conditions and
average magnetic fields.

(a) $\expectationvalue{B}\neq 0$

(b) $\expectationvalue{B}=0$

Figure 3:  Numerically obtained power laws for the average current, averaged
over $100$ realizations of the disorder, with fixed and free boundary
conditions. For $\expectationvalue{B}>0$, $B_{n}$ is chosen from $(1,3)$ while
for $\expectationvalue{B}=0$, $B_{n}$ is chosen from $(-2,2)$.

#### Free boundary conditions:

  * (a)

For $\expectationvalue{B}\neq 0$,
$\mathcal{T}_{\infty}(\omega)\sim\omega^{1/2}$ and $\lambda(\omega)\sim\omega$
which gives $\langle{J}_{N}\rangle\sim 1/N^{3/2}$ .

  * (b)

For $\expectationvalue{B}=0$, $\mathcal{T}_{\infty}(\omega)\sim\omega^{0}$ and
$\lambda(\omega)\sim\omega^{2/3}$ which gives $\langle{J}_{N}\rangle\sim
1/N^{3/2}$ .

The results are summarized in Table 1. Fig (3) shows the numerically obtained
power laws for $\expectationvalue{B}\neq 0$ and $\expectationvalue{B}=0$.
Numerically, the power laws are obtained by calculating
$\mathcal{T}_{N}(\omega)$ for different $\omega$ and then performing the
integration numerically. We expect to see the power law behaviour at some
large enough $N$. We see a reasonable agreement with the theoretically
expected power laws except for the case with $\expectationvalue{B}=0$ and free
BC, where we get $\langle{J}_{N}\rangle\sim 1/N^{2}$ instead of the expected
$\langle{J}_{N}\rangle\sim 1/N^{3/2}$ .

The case with $\expectationvalue{B}=0$ seems to be quite subtle because of the
following reasons:

  * •

The assumption that $\mathcal{T}_{\infty}(\omega)$ may be replaced by the
transmission for the uniform case for small $\omega$ does not hold good for
$\expectationvalue{B}=0$ case. This can be clearly seen from Fig. (4), where
we show a comparison of the transmission for small $\omega$ for
$\expectationvalue{B}\neq 0$ and $\expectationvalue{B}=0$ with their
respective uniform cases. While $\expectationvalue{B}\neq 0$ shows a clear
agreement with the corresponding uniform case, $\expectationvalue{B}=0$ case
shows a clear disagreement. It is not clear how to estimate
$\mathcal{T}_{\infty}(\omega)$ for this case.

(a) $\expectationvalue{B}\neq 0$

(b) $\expectationvalue{B}=0$

Figure 4:  Comparison for the transmission for disordered and uniform cases
for the two boundary conditions. For $\expectationvalue{B}\neq 0$, $B_{n}$ is
chosen from $(1,3)$ while for $\expectationvalue{B}=0$, $B_{n}$ is chosen from
$(-1,1)$. These are compared with the transmission for the uniform cases with
$B_{n}=\expectationvalue{B}$ respectively.

  * • (a) $\expectationvalue{B}\neq 0$ Figure 5:  Scaling of lowest allowed normal mode, $\omega_{s}^{N}$ with the system size, $N$. For $\expectationvalue{B}\neq 0$, $B_{n}$ is chosen from $(1,3)$ while for $\expectationvalue{B}=0$, $B_{n}$ is chosen from $(-1,1)$. The $B=0$ plot corresponds to the ordered chain (the ordered case $B\neq 0$ is not shown and has the scaling $N^{-2}$).

Interestingly we note that the transmission coefficient has peaks at much
lower frequencies than the ordered case. These peaks correspond to the normal
modes of the isolated chain and it is then of interest to study the system
size dependence of the lowest allowed normal mode frequency, $\omega_{s}^{N}$,
for the disordered chains with $\expectationvalue{B}\neq 0$ and
$\expectationvalue{B}=0$, and the ordered case with $B=0$. In Fig. (5) we show
the scaling of $\omega_{s}^{N}$ with $N$. We see that, for
$\expectationvalue{B}\neq 0$, $\omega_{\rm max}^{N}\sim 1/N$ while
$\omega_{s}^{N}\sim 1/N^{2}$. Thus, for any finite but large $N$, we have
$\omega_{\rm max}^{N}>\omega_{s}^{N}$ and there are a suffient number of
conducting modes. On the other hand, for $\expectationvalue{B}=0$, both
$\omega_{max}^{N}$ and $\omega_{s}^{N}$ scale as $1/N^{3/2}$ and this could be
the reason why our heuristic approach for current scaling fails for this case.

##  5 Conclusion

We considered a harmonic chain of charged particles in the presence of random
magnetic fields and derived power laws for the current with respect to the
system size. The power laws were found to be sensitive to boundary conditions
and the expectation value of the magnetic field. This was understood as
arising from the different behaviour of the Lyapunov exponent
$\lambda(\omega)$ and $\mathcal{T}_{\infty}(\omega)$ for small frequency
$\omega$.

Arguing that the small $\omega$ behaviour of ${\mathcal{T}}_{\infty}(\omega)$
was the same as that for the ordered chain, we used results obtained in our
previous paper [21] using the non-equilibrium Green function approach. It was
found there that this behaviour depends strongly on the presence or not of the
magnetic field but also on the boundary conditions imposed. To estimate the
Lyapunov exponent we mapped the discrete time process which determines the
Green’s functions to the motion of a harmonic oscillator with parametric
noise. This not only revealed an interesting connection between the Lyapunov
exponents of the two systems but also showed that the Lyapunov exponent have
different behaviour for different expectation values of the magnetic field.
For $\expectationvalue{B}>0$, $\expectationvalue{B}=0$ and
$\expectationvalue{B}<0$ we find that the Lyapunov exponents were of order
$\omega,\leavevmode\nobreak\ \omega^{2/3}$ and $\omega^{1/2}$ respectively.
These behaviours of the Lyapunov exponent were also verified numerically.

Using the results for the $\mathcal{T}_{\infty}(\omega)$ and
$\lambda(\omega)$, we make analytic predictions of different system-size
dependences of the current, depending on the expectation value of the magnetic
field and the boundary conditions. For free boundary conditions the current
decreases as $1/N^{3/2}$ irrespective of the expectation value of the magnetic
field. However, for fixed boundary conditions the current decreases as
$1/N^{3/2}$ and $1/N^{9/2}$ for $\expectationvalue{B}\neq 0$ and
$\expectationvalue{B}=0$ respectively. Our direct numerical estimates show
disagreement for the case $\expectationvalue{B}=0$, and this is especially
clear for the case with free boundary conditions. We discussed possible
reasons for the disagreement, amongst which is the intriguing numerical
observation of the $1/N^{3/2}$ system-size dependence of the lowest normal
mode frequency for the $\expectationvalue{B}=0$ case. The resolution of this
issue remains an interesting outstanding problem.

##  Appendix A Lyapunov exponent for a harmonic oscillator with parametric
noise

In order to obtain an expansion of $\lambda_{z}(\varepsilon)$ we follow the
strategy developed by Pardoux et al. in [26] and by Wihstutz in [27]. The
first step of the proof is to use the ergodic theorem to obtain an explicit
formula (see Eq. (24)) for $\lambda_{z}(\varepsilon)$ instead of Eq. (11). In
the second step we perform a perturbation analysis in $\varepsilon$ with this
new expression.

First we express the solution of the $2$-dimensional SDE $(z_{t})_{t\geq 0}$
in terms of a $1$-dimensional SDE. Define $(\theta_{t})_{t\geq 0}$ to be the
solution of

$\dot{\theta}_{t}=h_{0}(\theta_{t})+\frac{1}{2}\varepsilon^{2}\partial_{\theta}h_{1}(\theta_{t})h_{1}(\theta_{t})+\varepsilon
h_{1}(\theta_{t})\xi_{t}\ ,$ (19)

with

$h_{0}(\theta)=\sin^{2}(\theta)(c-1)-c\quad\text{and}\quad
h_{1}(\theta)=-\sigma\cos^{2}(\theta)\ .$ (20)

One can check that

$\displaystyle{z_{t}}=R_{t}\ (\cos(\theta_{t}),\sin(\theta_{t}))^{\top}$

where

$R_{t}=\|z_{0}\|\exp\left(\int_{0}^{t}\left[q_{0}\left(\theta_{\tau}\right)+\varepsilon^{2}r\left(\theta_{\tau}\right)\right]d\tau-\varepsilon\int_{0}^{t}q_{1}\left(\theta_{\tau}\right)\xi_{\tau}d\tau\right),$
(21)

with

$\displaystyle q_{0}(\theta)$ $\displaystyle=(1-c)\cos(\theta)\sin(\theta)\
,\quad q_{1}(\theta)=\sigma^{2}\cos(\theta)\sin(\theta)\ ,$ (22)
$\displaystyle r(\theta)$
$\displaystyle=\frac{\sigma^{2}\cos^{2}(\theta)}{2}\left[2\cos^{2}(\theta)-1\right]\
.$ (23)

Observe that $\|z_{t}\|=R_{t}$. Moreover, since in Eq. (19) the noise is
vanishing exactly at the points $\theta^{*}_{k}=(2k+1)\pi/2$,
$k\in\mathbb{Z}$, and that the drift in Eq. (19) at $\theta_{k}^{*}$ is equal
to $-1$, we see that starting from
$\theta_{0}\in[\theta_{k-1}^{*},\theta_{k}^{*})$ the process
$(\theta_{t})_{t\geq 0}$ will pass successively in the intervals
$\theta_{0}\in[\theta_{\ell}^{*},\theta_{\ell+1}^{*})$ for $\ell\leq k-1$
without coming back to an interval previously visited. This defines a sequence
of random times $t_{\ell}=\inf\\{t\geq
0\;;\;\theta_{t}\in[\theta_{k-\ell-1}^{*},\theta_{k-\ell}^{*})\\}$ for
$\ell\geq 0$ with $t_{0}=0$. The process is thus clearly not ergodic. A simple
way to restore this ergodicity (that will be needed later) is to consider the
process $({\tilde{\theta}}_{t})_{t\geq 0}$, living in $[-\pi/2,\pi/2)$, and
defined by $\tilde{\theta}_{t}=\theta_{t}+(k-\ell)\pi$ for
$t\in[t_{\ell},t_{\ell+1})$. The process ${\tilde{\theta}}_{t}$ satisfies the
same stochastic differential equation as $(\theta_{t})_{t\geq 0}$ but when it
reaches $-\pi/2$ it is immediately reseted to $\pi/2$. Equivalently
$({\tilde{\theta}}_{t})_{t\geq 0}$ is solution of Eq. (19) but seen as a SDE
on the torus $[-\pi/2,\pi/2)$ where the two end points of the interval have
been identified. The process $({\tilde{\theta}}_{t})_{t\geq 0}$ has now the
nice property to be ergodic. We denote by $\rho_{\varepsilon}(\theta)d\theta$
its invariant measure which is computed below. Observe moreover that Eq. (21)
still holds by replacing $\theta$ by $\tilde{\theta}$ because the functions
$q_{0},q_{1},r$ are $\pi$-periodic. In order to keep notation simple we denote
in the sequel the process $\tilde{\theta}$ by $\theta$.

By definition (11) of Lyapunov exponent and Eq. (21) we get that

$\displaystyle\lambda_{z}(\varepsilon)$ $\displaystyle=$
$\displaystyle\lim\limits_{t\rightarrow\infty}\frac{1}{t}\left\langle\int_{0}^{t}\left[q_{0}\left(\theta_{\tau}\right)+\varepsilon^{2}r\left(\theta_{\tau}\right)\right]d\tau+\varepsilon\int_{0}^{t}q_{1}\left(\theta_{\tau}\right)\xi_{\tau}d\tau\right\rangle$
$\displaystyle=$
$\displaystyle\lim\limits_{t\rightarrow\infty}\frac{1}{t}\left\langle\int_{0}^{t}\left[q_{0}\left(\theta_{\tau}\right)+\varepsilon^{2}r\left(\theta_{\tau}\right)\right]d\tau\right\rangle,$

since
$\langle\int_{0}^{t}q_{1}\left(\theta_{\tau}\right)\xi_{\tau}d\tau\rangle=0$.
Then by using the ergodic theorem we obtain

$\lambda_{z}(\varepsilon)=\int_{-\pi/2}^{\pi/2}\left[q_{0}(\theta)+\varepsilon^{2}r(\theta)\right]\rho_{\varepsilon}(\theta)d\theta\
.$ (24)

The expansion in $\varepsilon$ for $\lambda_{z}(\varepsilon)$ can then be
obtained from the expansion of $\rho_{\varepsilon}$.

Before doing this we prove Eq. (12), i.e. that the process $(z_{t})_{t\geq
0}=\left((u_{t},v_{t})^{\top}\right)_{t\geq 0}$ and the process
$(u_{t})_{t\geq 0}$ have the same Lyapunov exponent. By definition we have

$\lim\limits_{t\rightarrow\infty}\frac{1}{t}\left\langle\log|u_{t}|\right\rangle=\lim\limits_{t\rightarrow\infty}\frac{1}{t}\left\langle\log\|z_{t}\|\right\rangle+\lim\limits_{t\rightarrow\infty}\frac{1}{t}\left\langle\log|\cos(\theta_{t})|\right\rangle.$
(25)

Since $(\theta_{t})$ is an ergodic process we obtain that

$\lim\limits_{t\rightarrow\infty}\frac{1}{t}\left\langle\log|\cos(\theta_{t})|\right\rangle=\lim\limits_{t\rightarrow\infty}\frac{1}{t}\int_{-\pi/2}^{\pi/2}\rho_{\varepsilon}(\theta)\log\left(|\cos(\theta)|\right)d\theta=0\
.$

This proves the claim.

Let us now compute $\rho_{\varepsilon}$ which is the solution of the
stationary Fokker-Planck equation

$\partial_{\theta}\left[\tfrac{\varepsilon^{2}}{2}\partial_{\theta}(h_{1}^{2}\rho_{\varepsilon})-(h_{0}+\tfrac{\varepsilon^{2}}{2}h_{1}\partial_{\theta}h_{1})\rho_{\varepsilon}\right]=0\
.$ (26)

If we look for a solution such that
$\tfrac{\varepsilon^{2}}{2}\partial_{\theta}(h_{1}^{2}\rho_{\varepsilon})-(h_{0}+\tfrac{\varepsilon^{2}}{2}h_{1}\partial_{\theta}h_{1})\rho_{\varepsilon}=0$
we get
$\rho_{\varepsilon}(\theta)\propto\cos^{-2}(\theta)e^{-\frac{2\varepsilon^{-2}}{3\sigma^{2}}\tan^{3}(\theta)-\frac{2c\varepsilon^{-2}}{\sigma^{2}}\tan(\theta)}$
which is not normalisable. Hence we have to look for a normalisable solution
such that
$\tfrac{\varepsilon^{2}}{2}\partial_{\theta}(h_{1}^{2}\rho_{\varepsilon})-(h_{0}+\tfrac{\varepsilon^{2}}{2}h_{1}\partial_{\theta}h_{1})\rho_{\varepsilon}=A$
for some constant $A$. We get then that

$\rho_{\varepsilon}(\theta)=Z_{\varepsilon}^{-1}v_{\varepsilon}(\theta)\cos^{-2}(\theta)\int_{-\infty}^{\tan(\theta)}\
\exp\left(\frac{2\varepsilon^{-2}}{3\sigma^{2}}u^{3}+\frac{2c\varepsilon^{-2}}{\sigma^{2}}u\right)du$

with

$v_{\varepsilon}(\theta)=\exp\left\\{-\frac{2\varepsilon^{-2}}{3\sigma^{4}}\tan^{3}(\theta)-\frac{2c\varepsilon^{-2}}{\sigma^{4}}\tan(\theta)\right\\}$

and $Z_{\varepsilon}$ the partition function making $\rho_{\varepsilon}$ a
probability. Injecting this in Eq. (24) we may derive the results claimed by a
careful saddle point analysis. We prefer instead to rely on a more heuristic
analysis to bypass boring computations.

It is natural to expect that as $\varepsilon\to 0$ the stationary measure
$\rho_{\varepsilon}(\theta)d\theta$ will converge to the one of
${\dot{\theta}}_{t}=h_{0}(\theta_{t})$ (i.e. Eq. (19) with $\varepsilon=0$).
However as we will see this deterministic dynamical system has different
behaviours depending on the value of $c$ and that in some cases we have also
to compute the next order corrections.

If $c>0$, the deterministic dynamical system has a unique invariant state
$\rho_{0}(\theta)d\theta$ with
$\rho_{0}(\theta)=-\frac{\sqrt{c}}{\pi}h_{0}^{-1}(\theta)$ because $h_{0}$
never vanishes on $[-\pi/2,\pi/2)$. Hence $\rho_{\varepsilon}\to\rho_{0}$ as
$\varepsilon\to 0$. However, since
$\int_{0}^{\pi}q_{0}(\theta)\rho_{0}(\theta)d\theta=0$, we have to expand
$\rho_{\varepsilon}$ at order $\varepsilon^{2}$ to obtain the behavior of
$\lambda_{\varepsilon}$ in Eq. (24). Let us assume that
$\rho_{\varepsilon}=\rho_{0}+\varepsilon^{2}\delta\rho_{0}+o(\varepsilon^{2})$,
inject this in Eq. (26) and identify the powers in $\varepsilon$. We obtain
that

$\partial_{\theta}[h_{0}\
(\delta\rho_{0})]=\tfrac{1}{2}\partial_{\theta}\left[\partial_{\theta}(h_{1}^{2}\rho_{0})-(h_{1}\partial_{\theta}h_{1})\rho_{0}\right]$

which implies, since $\int_{-\pi/2}^{\pi/2}(\delta\rho_{0})(\theta)d\theta=0$
that

$\delta\rho_{0}=\tfrac{A}{h_{0}}+\tfrac{1}{2h_{0}}\left[\partial_{\theta}(h_{1}^{2}\rho_{0})-(h_{1}\partial_{\theta}h_{1})\rho_{0}\right].$

We deduce that

$\delta\rho_{0}=\frac{A}{h_{0}}+\frac{\sigma^{2}\sqrt{c}}{\pi}\left(\frac{\sin(\theta)\cos^{3}(\theta)}{h_{0}^{2}}+(c-1)\left(\frac{\cos^{5}(\theta)\sin(\theta)}{h_{0}^{3}(\theta)}\right)\right).$

Since $\int_{-\pi/2}^{\pi/2}(\delta\rho_{0})(\theta)d\theta=0$ we obtain $A=0$
and

$\delta\rho_{0}=\frac{\sigma^{2}\sqrt{c}}{\pi}\left(\frac{\sin(\theta)\cos^{3}(\theta)}{h_{0}^{2}}+(c-1)\left(\frac{\cos^{5}(\theta)\sin(\theta)}{h_{0}^{3}(\theta)}\right)\right).$

Hence we get that

$\lambda_{z}(\varepsilon)=\varepsilon^{2}\int_{-\pi/2}^{\pi/2}\left(r(\theta)\rho_{0}(\theta)+q_{0}(\theta)\delta\rho_{0}(\theta)\right)d\theta+o(\varepsilon^{2})\
.$

By the change of variable $x=\tan(\theta)$ we get

$\displaystyle\int_{-\pi/2}^{\pi/2}r(\theta)\rho_{0}(\theta)d\theta$
$\displaystyle=\frac{\sigma^{2}\sqrt{c}}{2\pi}\int_{-\infty}^{\infty}\frac{x^{2}-1}{(1+x^{2})(x^{2}+c)}dx=\frac{\sigma^{2}}{2(\sqrt{c}+1)^{2}}.$
$\displaystyle\int_{-\pi/2}^{\pi/2}q_{0}(\theta)\delta\rho_{0}(\theta)d\theta$
$\displaystyle=\sigma^{2}\left(\frac{(4\sqrt{c}+1)(c-1)^{2}}{8(\sqrt{c}+1)^{4}c}+\frac{1-c}{2(\sqrt{c}+1)^{3}}\right).$

Hence we finally get

$\lambda_{z}(\varepsilon)=\varepsilon^{2}\frac{\sigma^{2}}{8c}+o(\varepsilon^{2}).$

This proves case (ii).

If $c<0$ then $\tfrac{c}{c-1}\in(0,1)$ and the function $h_{0}$ vanishes on
$[-\pi/2,\pi/2)$ if and only if $\theta\in[-\pi/2,\pi/2)$ is solution of
sin^2(θ) = cc-1 . There are two solutions $\theta^{*}>0$ and $-\theta^{*}<0$.
The deterministic dynamical has two extremal invariant probability measures
given by $\delta_{\pm\theta^{\star}}$. Since
$h^{\prime}_{0}(\theta^{*})<0<h^{\prime}_{0}(-\theta^{*})$,
$\delta_{-\theta^{*}}$ is unstable while $\delta_{\theta_{*}}$ is stable. By
introducing noise in this dynamical system the stable stationary state is
selected when the intensity of the noise is sent to zero afterwards, i.e.
$\rho_{\varepsilon}(\theta)d\theta\to\delta_{\theta^{*}}$. We conclude that

$\lim_{\varepsilon\to 0}\lambda_{z}(\varepsilon)=q_{0}(\theta^{*})=\sqrt{|c|}\
.$

This proves case (iii).

The case $c=0$ is more delicate. Since $h_{0}(\cdot)=-\sin^{2}(\cdot)$, the
unique invariant measure for the deterministic dynamical system is
$\delta_{0}$ (stable) and we expect that
$\rho_{\varepsilon}(\theta)d\theta\to\delta_{0}$ as $\varepsilon\to 0$.
Observe however that $q_{0}(0)=0$ so that we have to find the first correction
to the approximation of $\rho_{\varepsilon}$ to $\delta_{0}$. Due to the
singularity of the Dirac mass we cannot perform an expansion analysis in
$\varepsilon$. Hence we will use another argument to get item (i). Consider
the following linear transformation
$T_{\varepsilon}=\left(\begin{smallmatrix}\varepsilon^{2/3}&0\\\
0&1\end{smallmatrix}\right)$ which is such that
$\varepsilon^{2/3}\|z\|\leq\|T_{\varepsilon}z\|\leq\|z\|$ for any
$z\in\mathbb{R}^{2}$ and $\varepsilon\leq 1$. This implies that
$(z_{t})_{t\geq 0}$ and $(T_{\varepsilon}z_{t})_{t\geq 0}$ have the same
Lyapunov exponent. Expressing as we did before

${\hat{z}}_{t}:=T_{\varepsilon}z_{t}=\|{\hat{z}}_{t}\|(\cos{\hat{\theta}}_{t},\sin{\hat{\theta}}_{t})^{\top}$

we notice that

$\dot{\hat{\theta}}_{t}=\varepsilon^{2/3}\left(-\sin^{2}({\hat{\theta}}_{t})-\sigma^{2}\sin({\hat{\theta}}_{t})\cos^{3}({\hat{\theta}}_{t})\right)-\varepsilon^{1/3}\sigma\cos^{2}({\hat{\theta}}_{t})\xi_{t}\
,$

which implies by scaling invariance of the white noise that
${\hat{\theta}}_{t}=\alpha_{t\varepsilon^{2/3}}$ where

$\dot{\alpha}_{t}=\left(-\sin^{2}(\alpha_{t})-\sigma^{2}\sin(\alpha_{t})\cos^{3}(\alpha_{t})\right)-\sigma\cos^{2}({\alpha}_{t})\xi_{t}\
.$

If ${\hat{\rho}}(\alpha)d\alpha$ is the unique invariant measure for
$(\alpha_{t})_{t\geq 0}$ we have by a scaling argument that the Lyapunov
exponent satisfies

$\lambda_{z}(\varepsilon)=\varepsilon^{2/3}{\hat{\lambda}}(\sigma)$

with

$\hat{\lambda}(\sigma)=\int_{-\pi/2}^{\pi/2}\left(q_{0}(\alpha)+r(\alpha)\right)\hat{\rho}(\alpha)d\alpha,$
(27)

where $q_{0}$ and $r$ are defined respectively in Eq. (22) and Eq. (23) with
$c=0$. To obtain the value of $\hat{\lambda}(\sigma)$ it is sufficient to find
$\hat{\rho}$ which is the unique normalisable function of the Fokker-Planck
equation associated to the process$(\alpha_{t})_{t\geq 0}$, i.e.

${\hat{\rho}}(\alpha)={\hat{Z}}^{-1}\cos^{-2}(\alpha)e^{-\tfrac{2}{3\sigma^{2}}\tan^{3}(\alpha)}\
\int_{-\infty}^{\tan(\alpha)}\exp(\tfrac{2u^{3}}{3\sigma^{2}})du\ ,$

where ${\hat{Z}}$ is the normalisation constant making $\hat{\rho}$ a
probability measure.

## Acknowledgements

A.D. and J.M.B. acknowledge support of the Department of Atomic Energy,
Government of India, under Project No. RTI4001. The work of C.B. and G.C. has
been supported by the projects LSD ANR-15-CE40- 0020-01 of the French National
Research Agency (ANR), by the European Research Council (ERC) under the
European Unions Horizon 2020 research and innovative programme (grant
agreement No 715734) and the French-Indian UCA project ‘Large deviations and
scaling limits theory for non-equilibrium systems’.

## References

## References

  * [1] P. W. Anderson, “Absence of diffusion in certain random lattices,” Phys. Rev., vol. 109, pp. 1492–1505, Mar 1958.
  * [2] A. Casher and J. L. Lebowitz, “Heat flow in regular and disordered harmonic chains,” Journal of Mathematical Physics, vol. 12, no. 8, pp. 1701–1711, 1971.
  * [3] R. J. Rubin and W. L. Greer, “Abnormal lattice thermal conductivity of a one‐dimensional, harmonic, isotopically disordered crystal,” Journal of Mathematical Physics, vol. 12, no. 8, pp. 1686–1701, 1971.
  * [4] A. J. O’Connor and J. L. Lebowitz, “Heat conduction and sound transmission in isotopically disordered harmonic crystals,” J. Mathematical Phys., vol. 15, pp. 692–703, 1974.
  * [5] T. Verheggen, “Transmission coefficient and heat conduction of a harmonic chain with random masses: asymptotic estimates on products of random matrices,” Comm. Math. Phys., vol. 68, no. 1, pp. 69–82, 1979.
  * [6] O. Ajanki and F. Huveneers, “Rigorous scaling law for the heat current in disordered harmonic chain,” Comm. Math. Phys., vol. 301, no. 3, pp. 841–883, 2011.
  * [7] A. Dhar, “Heat conduction in the disordered harmonic chain revisited,” Phys. Rev. Lett., vol. 86, pp. 5882–5885, Jun 2001.
  * [8] D. Roy and A. Dhar, “Role of pinning potentials in heat transport through disordered harmonic chains,” Phys. Rev. E, vol. 78, p. 051112, Nov 2008\.
  * [9] W. De Roeck, A. Dhar, F. Huveneers, and M. Schütz, “Step density profiles in localized chains,” J. Stat. Phys., vol. 167, no. 5, pp. 1143–1163, 2017\.
  * [10] C. Bernardin, F. Huveneers, and S. Olla, “Hydrodynamic limit for a disordered harmonic chain,” Comm. Math. Phys., vol. 365, no. 1, pp. 215–237, 2019\.
  * [11] A. Amir, Y. Oreg, and Y. Imry, “Thermal conductivity in 1d: Disorder-induced transition from anomalous to normal scaling,” EPL (Europhysics Letters), vol. 124, no. 1, p. 16001, 2018.
  * [12] B. Ash, A. Amir, Y. Bar-Sinai, Y. Oreg, and Y. Imry, “Thermal conductance of one-dimensional disordered harmonic chains,” Physical Review B, vol. 101, no. 12, p. 121403, 2020.
  * [13] L. W. Lee and A. Dhar, “Heat conduction in a two-dimensional harmonic crystal with disorder,” Phys. Rev. Lett., vol. 95, p. 094302, Aug 2005.
  * [14] A. Chaudhuri, A. Kundu, D. Roy, A. Dhar, J. L. Lebowitz, and H. Spohn, “Heat transport and phonon localization in mass-disordered harmonic crystals,” Phys. Rev. B, vol. 81, p. 064301, Feb 2010.
  * [15] A. Dhar and J. L. Lebowitz, “Effect of phonon-phonon interactions on localization,” Phys. Rev. Lett., vol. 100, p. 134301, Apr 2008.
  * [16] A. Dhar and K. Saito, “Heat conduction in the disordered fermi-pasta-ulam chain,” Phys. Rev. E, vol. 78, p. 061136, Dec 2008.
  * [17] C. Bernardin, “Thermal conductivity for a noisy disordered harmonic chain,” J. Stat. Phys., vol. 133, no. 3, pp. 417–433, 2008.
  * [18] A. Dhar, K. Venkateshan, and J. L. Lebowitz, “Heat conduction in disordered harmonic lattices with energy-conserving noise,” Phys. Rev. E, vol. 83, p. 021108, Feb 2011.
  * [19] C. Bernardin and F. Huveneers, “Small perturbation of a disordered harmonic chain by a noise and an anharmonic potential,” Probab. Theory Related Fields, vol. 157, no. 1-2, pp. 301–331, 2013.
  * [20] C. Bernardin, F. Huveneers, J. L. Lebowitz, C. Liverani, and S. Olla, “Green-Kubo formula for weakly coupled systems with noise,” Comm. Math. Phys., vol. 334, no. 3, pp. 1377–1412, 2015.
  * [21] J. Bhat, G. Cane, C. Bernardin, and A. Dhar, “Heat transport in an ordered harmonic chain in presence of a uniform magnetic field,” arXiv-2106.12069, 2021.
  * [22] H. Furstenberg, “Noncommuting random products,” Trans. Amer. Math. Soc., vol. 108, pp. 377–428, 1963.
  * [23] H. Crauel, P. Imkeller, and M. Steinkamp, Bifurcations of One-Dimensional Stochastic Differential Equations.  New York, NY: Springer New York, 1999.
  * [24] V. Wihstutz, Perturbation Methods for Lyapunov Exponents, pp. 209–239.  New York, NY: Springer New York, 1999.
  * [25] H. Matsuda and K. Ishii, “Localization of normal modes and energy transport in the disordered harmonic chain,” Progr. Theoret. Phys., no. Suppl. 45, pp. 56–86, 1970.
  * [26] E. Pardoux and V. Wihstutz, “Lyapunov exponent and rotation number of two-dimensional linear stochastic systems with small diffusion,” SIAM Journal on Applied Mathematics, vol. 48, no. 2, pp. 442–457, 1988.
  * [27] V. Wihstutz, Ergodic Theory of Linear Parameter-Excited Systems, pp. 205–218.  Dordrecht: Springer Netherlands, 1981.