\section{Introduction}
The Weyl anomaly in quantum field theory has a long history (see e.g. \cite{Duff:1993wm} for a historical overview). It states that the Weyl transformation of the metric in the curved space-time may not be a symmetry of the system even though the quantum field theory under consideration has the conformal symmetry in the flat space-time. Indeed, in most conformal field theories, the Weyl anomaly is non-vanishing, and we say that the Weyl symmetry is quantum mechanically broken in curved space-time.

The Weyl anomaly has played many important roles in theoretical physics.
The existence of the Weyl anomaly gives a constraint on the expectation values of energy-momentum tensor in curved background \cite{Page:1982fm}, and may be related to the nature of Hawking radiation \cite{Robinson:2005pd}\cite{Iso:2006wa}\cite{Kawai:2017txu}\cite{Kawai:2014afa}.  The fact that vanishing of the Weyl anomaly happens only in a limited class of theories dictates the number of space-time dimensions in critical string theory.
More recently, we find that the universal terms in the entanglement entropy of conformal field theories are given by the coefficient of the Weyl anomaly, suggesting a deep relation between geometry and information \cite{Ryu:2006ef}.

What we would like to study in this paper is to find a way to cancel the Weyl anomaly from the other source, e.g. from the position dependent coupling constant.\footnote{Sorry for the oxymoron. It is no longer constant. Probably it is Dirac who openly advocated this idea in the early days \cite{Dirac:1938mt}.} Once we put a quantum field theory on a curved manifold, it is natural to assume that coupling constants are position dependent. The position dependent coupling constants then provide an extra contribution to the Weyl anomaly so that we may attempt to cancel the entire Weyl anomaly on the curved manifold. We would like to find under which condition such a cancellation is possible.

The similar idea of cancelling more general anomalies have been implicitly assumed in many places. For example, if we try to introduce background gauge fields for chiral current operators in the curved background (e.g. in the context of supersymmetric localization), then they may be mutually inconsistent due to the 't Hooft anomaly. One way to avoid this is to cancel the anomaly of the background gauge field from the space-time curvature and vice versa. Similarly, if preserving the Weyl anomaly is the critical issue (e.g. if we try to gauge it), our new way of doing it may be another option to be considered.

The organization of the paper is as follows. In section 2, we study the cancellation of the Weyl anomaly from position dependent coupling constant in two dimensional conformal field theories. In section 3, we study it in three dimensions and in section 4, we study it in four dimensions. We supplement the holographic viewpoint in section 5 and conclude with some discussions in section 6.

\section{Two dimensions}
Let us consider a two-dimensional conformal field theory with an exactly marginal deformation denoted by $g$. For instance, we may take the Gaussian $c=1$ boson with the compactification radius as the exactly marginal deformation here.
We put the theory on a curved background with the metric $g_{\mu\nu}(x)$ and then vary the coupling constant $g(x)$ over the manifold: schematically we consider the action $S = S_0 + \int d^2x\sqrt{g} g(x) O(x)$. Even though the theory is conformal invariant in the Euclidean space $g_{\mu\nu} = \delta_{\mu\nu}$ with $g(x) = g$, it is not necessarily so after turning on the background metric and position dependent coupling constant.\footnote{We define the scale transformation  by the change of the difference of the coordinate as in \cite{Osborn:1991gm} rather than the coordinate itself \cite{Dong:2012ua}.}
This obstruction is known as the Weyl anomaly under the infinitesimal Weyl rescaling: $\delta g_{\mu\nu}(x) = 2\delta \sigma g_{\mu\nu}$.

In terms of the free energy functional $e^{-F[g_{\mu\nu}(x),g(x)]} = \int \mathcal{D} \Phi e^{-S[\Phi]}$, the Weyl anomaly for a two-dimensional conformal field theory is given by (e.g. see  \cite{Osborn:1991gm})\footnote{In this paper, we always assume that the conformal field theories under consideration preserve the CP symmetry.}
\begin{align}
\delta F_{\sigma} = \int d^2x\sqrt{g} \delta \sigma(x) (cR -\frac{1}{2} \partial^\mu g(x) \partial_\mu g(x)) \  \label{twoa}
\end{align}
in a certain renormalization scheme so that the exactly marginal deformation has a flat line metric. Otherwise, we can always redefine the coupling constant or the renormalization scheme so that it is flat. It is clear that when $g(x) = g$, the only way to cancel the Weyl anomaly is to require $c=0$, which is typically what we demand in critical string theory.

However, we see that this is not the only available option. Now, given a positive curvature $R(x)\ge 0$, we may try to cancel the curvature term in the Weyl anomaly \eqref{twoa} against the second term originating from the position dependent coupling constant by solving the equation
\begin{align}
c R = \frac{1}{2} \partial^\mu g(x) \partial_\mu g(x) \ . \label{twoc}
\end{align}
This is possible for positive curvature $R(x) \ge 0$ (assuming $c>0$ in unitary conformal field theories).

For example, if we take the Fubini-Study metric on the sphere with the complex coordinate $z$ and $\bar{z}$:
\begin{align}
ds^2 = \frac{dzd\bar{z}}{(1+|z|^2)^2} \ ,
\end{align}
the solution of \eqref{twoc} is
\begin{align}
g(x) = \sqrt{c} \cdot \mathrm{arctan}(|z|)  + \mathrm{const} \ .
\end{align}
In this way, one can cancel the Weyl anomaly on the sphere by introducing the position dependent coupling constant. Note, however, that the position dependence of the coupling constant reduces the symmetry of the sphere from $SO(3)$ down to $SO(2)$. The idea here is we gained extra ``Weyl symmetry" at the sacrifice of the rotational symmetry.\footnote{This is not necessarily a bad idea: for examle, in the supersymmetric localization, we often do not keep the full isometry of the sphere but only the $U(1)$ subgroup of it.}

The above cancellation works both for infinitesimal generic Weyl transformation $\delta g_{\mu\nu}(x) = 2\delta \sigma(x) g_{\mu\nu}(x)$ or finite but constant Weyl transformation $g_{\mu\nu}(x) \to e^{2\bar{\sigma}} g_{\mu\nu}(x)$, where $\bar{\sigma}$ is a finite constant. The latter is because the equation to be solved in \eqref{twoc} trivially scales under the constant Weyl transformation, so once it is solved then it is also solved after finite but constant Weyl transformation.
For finite generic Weyl transformation, however, the cancellation may not persist. The point is that the curvature term in the Weyl anomaly is non-trivially transforms under the Weyl transformation:
\begin{align}
R \to e^{-2\sigma(x)} (R -2D^2 \sigma) \ ,
\end{align}
where $D^2$ is the Laplacian,
while the Weyl transformation of the second term from the position dependent coupling constant $\partial^\mu g(x) \partial_\mu g(x)$ is trivial:
\begin{align}
\partial^\mu g(x) \partial_\mu g(x) \to e^{-2\sigma(x)} \partial^\mu g(x) \partial_\mu g(x)
\end{align}
Thus, even though one may solve the cancellation condition for a given $g_{\mu\nu}(x)$ with a certain position dependent coupling constant $g(x)$, the cancellation does not persist for the Weyl transformed geometry.

 Nevertheless, we realize that the cancellation is still intact if we restrict\footnote{A similar restriction on the Weyl transformation has been studied in \cite{Edery:2014nha}.} ourselves to the harmonic Weyl transformation, which satisfies $D^2 \sigma = 0$. Thus, we may construct a quantum field theory which is exactly invariant under the harmonic Weyl transformation by cancelling the Weyl anomaly from the position dependent coupling constant.

More generically, one may consider theories with several exactly marginal deformations. The Weyl anomaly has the generalized form
\begin{align}
\delta F_{\sigma} = \int d^2x\sqrt{g} \delta\sigma(x) (cR - \chi_{ij}(g) \partial^\mu g^i(x) \partial_\mu g^j(x)) + \partial_\mu \delta\sigma(x) w_i(g) \partial^\mu g^i(x)  \ , \label{twog}
\end{align}
where $\chi_{ij}(g)$ and $w_i(g)$ may depend on the exactly marginal deformations $g^i(x)$.
For a constant Weyl transformation, the condition for the cancellation is essentially the same as before since the last term in \eqref{twog} drops out.

For infinitesimal generic Weyl transformation, however, we have to think about the cancellation of the third term proportional to $\partial_\mu \delta \sigma(x)$. We did not talk about it in the single coupling case because we were able to remove it from the local counterterm $\int d^2x\sqrt{g} b(g) R$, but we have to discuss it now with several coupling constants when it has the non-trivial curvature $\partial_i w_j - \partial_j w_i$.
While the Wess-Zumino consistency condition does not say anything about the (non-)existence of this term  \cite{Osborn:1991gm}, the recent analysis in \cite{Gomis:2015yaa} tells that on the conformal manifold spanned by the exactly marginal deformations, the curvature is trivial (i.e.  $\partial_i w_j - \partial_j w_i = 0$) and can be removed by the local counterterm $\int d^2x\sqrt{g} b(g^i) R$,
 so we actually do not have to worry about its cancellation. The non-existence of the curvature $\partial_i w_j -\partial_j w_i$ is related to the gradientness of the beta functions and it may have a deep implication in renormalization group flows  \cite{Osborn:1991gm}\cite{Friedan:2009ik}\cite{Gukov:2016tnp}.

\section{Three dimensions}
There is no curvature dependent Weyl anomaly in three dimensions. The position dependent exactly marginal deformations do not introduce the additional Weyl anomaly either under the assumption of the CP symmetry \cite{Nakayama:2013wda}. Thus there is no interesting scenario we can imagine in three dimensions.

\section{Four dimensions}
Let us consider a four-dimensional conformal field theory with an exactly marginal deformation denoted by $g$. We put the theory on a curved background with the metric $g_{\mu\nu}(x)$ and then vary the coupling constant over the manifold $g(x)$.
In a certain renormalizaiton scheme, the first order Weyl transformation (i.e. the Weyl anomaly) is given by \cite{Osborn:1991gm} (See also \cite{Nakayama:2013is}\cite{Jack:2013sha}.)
\begin{align}
\delta F_{\sigma}  = \int d^4x \sqrt{g} \delta\sigma(x) &\left( c(g) \mathrm{Weyl}^2 - a \mathrm{Euler} + (D^2 g D^2g - 2G_{\mu\nu} \partial^\mu g \partial^\nu g -\frac{R}{3}\partial^\mu g\partial_\mu g)  \right. \cr
& \left. + \chi_4(g) \partial_\mu g \partial^\mu g \partial_\nu g \partial^\nu g \right) \ .  \label{foura}
\end{align}
Here $G_{\mu\nu} = R_{\mu\nu} - \frac{R g_{\mu\nu}}{2}$ is the Einstein tensor and $D^\mu$ is the covariant derivative. In addition, we have introduced $\mathrm{Weyl}^2 = R_{\mu\nu\rho\sigma}^2 -2R_{\mu\nu}^2 + \frac{1}{3}R^2$ and $\mathrm{Euler} = R_{\mu\nu\rho\sigma}^2 - 4R_{\mu\nu}^2 + R^2$.
In principle $c(g)$ can depend on $g$, but the only such theories known are constructed in somewhat artificial holographic realization \cite{Nakayama:2017oye}.

To simplify the analysis, let us focus on the regime in which the last quartic term in \eqref{foura} i.e. $\chi_4(g) \partial_\mu g \partial^\mu g \partial_\nu g \partial^\nu g$ can be neglected (e.g. in the small coupling regime).
Neglecting the quartic term, we try to solve the equation
\begin{align}
-c\mathrm{Weyl}^2 + a \mathrm{Euler} = (D^2 g D^2g - 2G_{\mu\nu} \partial^\mu g \partial^\nu g -\frac{R}{3}\partial^\mu g\partial_\mu g) \ . \label{fourc}
\end{align}
In particular, suppose that the metric $g_{\mu\nu}(x)$ is Ricci flat. Then the equation \eqref{fourc} becomes
\begin{align}
\sqrt{(a-c) R_{\mu\nu\rho\sigma}^2} = D^2 g \ ,
\end{align}
which may be solved by using Green's function for the Laplacian
\begin{align}
g(x) = \int d^4x' G(x,x') \sqrt{(a-c) R_{\mu\nu\rho\sigma}^2(x')}
\end{align}
when the manifold is non-compact (otherwise the regular solution does not exist).

As in two-dimensions, the above argument works both for infinitesimal Weyl transformation or finite but constant Weyl transformation.
For finite generic Weyl transformation, one may define the analogue of harmonic Weyl transformation. For this purpose, it is more convenient to choose a different renormalization scheme so that the Weyl anomaly takes the form (e.g. \cite{Gomis:2015yaa})
\begin{align}
\delta F_{\sigma}  = \int d^4x \sqrt{g} \delta\sigma(x) &\left( (c(g)-a) \mathrm{Weyl}^2 - 4a Q   + g \Delta_4 g \right. \cr
& \left. + \chi_4(g) \partial_\mu g \partial^\mu g \partial_\nu g \partial^\nu g \right) \ , \label{alt}
\end{align}
where $\Delta_4$ is the Fradkin-Tseytlin-Riegert-Paneitz conformal operator \cite{Fradkin:1982xc}\cite{Fradkin:1981jc}\cite{Riegert:1984kt}\cite{PA}
\begin{align}
\Delta_4 = (D^2)^2 + 2G_{\mu\nu}D^\mu D^\nu + \frac{1}{3}(D^\mu R) D_\mu + \frac{1}{3}R D^2,
\end{align}
which is Weyl covariant $\Delta_4 \to e^{-4\sigma} \Delta_4$,
and $Q$ is what is called the Q-curvature \cite{Q}:
\begin{align}
Q = \frac{-1}{6} D^2 R -\frac{1}{2}R^{\mu\nu}R_{\mu\nu} + \frac{1}{6}R^2 \
\end{align}
which has a nice mathematical property under the Weyl transformation
\begin{align}
Q \to e^{-4\sigma} (Q + \Delta_4 \sigma)
\end{align}

The advantage of this rewriting or a choice of the particular local counterterm is as follows. Suppose we cancelled the Weyl anomaly at a particular background by demanding
\begin{align}
0 = (c(g)-a) \mathrm{Weyl}^2 - 4a Q   + g \Delta_4 g + \chi_4(g) \partial_\mu g \partial^\mu g \partial_\nu g \partial^\nu g
\end{align}
Then, we are still able to cancel the Weyl anomaly on the Weyl transformed manifold whenever the Weyl rescaling is annihilated by the Fradkin-Tseytlin-Riegert-Paneitz operator:
\begin{align}
\Delta_4 \sigma = 0 \ . \label{ann}
\end{align}
This is because all the terms in \eqref{alt} except for the Q-curvature transform covariantly under the finite Weyl transformation.
If the Weyl scaling factor satisfies \eqref{ann}, the cancellation of the Weyl anomaly therefore persists even for finite Weyl transformation. This is precisely analogous to the special role of harmonic Weyl transformation in two dimensions.

Let us move on to the most generic cases with multiple coupling constants. The Weyl transformation is given by
\begin{align}
\delta F_{\sigma}  = \int d^4x \sqrt{g} \delta \sigma(x) &( c(g) \mathrm{Weyl}^2 - a \mathrm{Euler} +  \chi_{ij}(g) (D^2 g^i D^2g^j - 2G_{\mu\nu} \partial^\mu g^i \partial^\nu g^j -\frac{R}{3}\partial^\mu g^i \partial_\nu g^j)  \cr
& \left. + \chi_{ijkl}(g) \partial_\mu g^i \partial^\mu g^j \partial_\nu g^k \partial^\nu g^l \right)  \cr
&+ \partial_\mu \delta\sigma G^{\mu\nu} w_i(g) \partial_\nu g^j \ . \label{fourm}
\end{align}
For finite but constant Weyl transformation, we only have to cancel the first two lines in \eqref{fourm}, which is  essentially equivalent to what we have done in the above. On the other hand, for infinitesimal but generic Weyl transformation, we have to cancel the third line as well, which requires either $\partial_i w_j - \partial_j w_i =0$ or $G_{\mu\nu} = 0$ in the background.

To conclude the analysis, we would like to mention the other obstructions to the Weyl transformation if the dimension two operator $O(x)$ exist in the theory. If this is the case, there is a further operator Weyl anomaly such as
\begin{align}
\int d^4x \sqrt{g} \delta\sigma(x) (\eta(g) R O(x) + \epsilon(g) \partial_\mu g \partial^\mu g O + \tau(g)  D^2 O + \delta(g) D^2 g O) + \partial_\mu \delta \sigma(x) \theta(g) \partial^\mu g O \ .
\end{align}
It has been shown that such Weyl anomaly can be removed when $g(x) = g$ \cite{Jack:2013sha} (see also \cite{Nakayama:2013wda}\cite{Farnsworth:2017tbz} for similar analysis), but with the space-time dependent coupling, we need the extra cancellation to get the consistent picture. Schematically, the Wess-Zumino consistency condition demands $\eta = 0$ and  one can always remove $\theta$ and $\tau$ by local counterterms. Then we need to cancel $\epsilon$ term against the $\delta$ term. Since the existence of dimension two operator is non-generic, we will not pursue the cancellation in further details.

\section{Holographic models}
We revisit the cancellation mechanism we have studied in previous sections from  the holographic perspective. For definiteness we consider the case of four dimensional conformal field theories with the five dimensional bulk. Let us study the Einstein gravity coupled with a scalar field $\phi$ given by the minimal action
\begin{align}
S = \int d^5x \sqrt{g} \left( R + \Lambda + \frac{1}{2}\partial^M \phi \partial_M \phi  \right) .
\end{align}

In the AdS/CFT correspondence, we compute the on-shell action for a given boundary condition at $\rho= \epsilon$ (i.e. $\phi_{(0)}(x)$ and $g_{(0)\mu\nu}(x)$ below) with the expansion
\begin{align}
\phi &= \phi_{(0)}(x) + \rho \phi_{(1)}(x) + \rho^2 \phi_{(2)}(x) + \cdots \cr
g_{\mu\nu} &= g_{(0)\mu\nu}(x) + \rho g_{(1)\mu\nu}(x) + \rho^2 g_{(2) \mu\nu}(x) + \cdots
\end{align}
in the Graham-Fefferman gauge
\begin{align}
ds^2 = G_{MN} dx^M dx^N = \frac{d\rho^2}{\rho^2} +  \frac{g_{\mu\nu} dx^\mu dx^\nu}{\rho}  \ .
\end{align}

The resulting on-shell action is generically divergent in the limit $\epsilon \to 0$ from the $\rho$ integration as $\int_{\epsilon} d\rho \rho^{-1} S_{\log} = \log \epsilon S_{\log}$, leading to the holographic Weyl anomaly \cite{Henningson:1998gx}. Explicitly \cite{Nojiri:1998dh}, we have
\begin{align}
S = \log \epsilon \int d^4x &\left( \frac{1}{8} R_{\mu\nu (0)}^2 -\frac{1}{24} R_{(0)}^2  + \frac{1}{4} (D^2 \phi_{(0)})^2 \right. \cr
& \left.  - \frac{1}{2} R_{(0)}^{\mu\nu} \partial_\mu \phi \partial_\nu \phi + \frac{1}{6} R_{(0)} \partial^\mu \phi_{(0)} \partial_\nu \phi_{(0)} + \frac{1}{3}(\partial_\mu \phi_{(0)} \partial^\mu\phi_{(0)})^2  \right) \
\end{align}
which is exactly what we had in section 4 for the constant Weyl transformation.
Thus cancelling the Weyl anomaly from the position dependent coupling constant corresponds to the choice of the boundary values of $\phi(x)$ such that the on-shell gravity action is finite without the logarithmic divergence. The choice of such boundary conditions make the AdS/CFT correlation functions better behaved, so classifying such supergravity background may be of theoretical interest.

\section{Discussions}
In this paper, we have studied a novel way to cancel the Weyl anomaly from the position dependent coupling constant. Here we would like to mention further possibilities to cancel the Weyl anomaly.

First of all, if the theory under consideration possesses a conserved current $J_\mu$, one may introduce the background field strength by the coupling $\int d^4x \sqrt{g} A^\mu J_\mu$. This gives another contribution to the Weyl anomaly as $\int d^4x \delta\sigma(x) b_0 F_{\mu\nu}(x) F^{\mu\nu}(x)$, where $b_0$ is the coefficient of the one-loop beta function determined from the current two-point function (which is positive in unitary conformal field theories) and $F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu$. Then we may try to cancel the Weyl anomaly from this contribution.
Actually, simultaneous use of the gauge field and the position dependent coupling constant may not be a good idea because of the existence of the vector beta functions \cite{Nakayama:2013ssa}. Again we have to think about the cancellation of the extra operator Weyl anomaly such as $\int d^4x \delta \sigma \rho(g) \partial^\mu g J_\mu$. To avoid the appearance of the vector beta functions, we may only introduce the position dependent coupling constant which is neutral under the symmetry generated by $J_\mu$.

It could have been extremely interesting if we were able to find a novel class of Weyl gauging without demanding $c=0$ in two dimensions, and $c=a=0$ in four dimensions.\footnote{See e.g. \cite{Fradkin:1983tg} for a possibility in the context of supersymmetric Weyl gravity.}
Currently, the closest way to do this is to demand all the non-trivial Weyl anomalies vanish, say $a=0$ in four dimensions, and then try to cancel the $c \mathrm{Weyl}^2$ term against the space-time dependent coupling constants. Here, we should further employ the non-unitariness of the model (since $a=0$ from the beginning suggests it must be so) to obtain the cancellation in the Weyl anomalies. This is because unitarity demands the positivity of the both terms and the cancellation only happens by using the non-unitary property. Whether this is better than just demanding $c=a=0$ is yet to be seen in the context of quantum Weyl gravity in which we would like to gauge the Weyl symmetry exactly.\footnote{Note, however, that the introduction of the position dependent coupling constant modifies the conservation of the energy-momentum tensor, and the consistecy with the dynamical gravity must be discussed more carefully. We stress that in the main part of this paper, we have focused on the fixed gravitational background, so there is no inconsistency. The author would like to thank S.~Deser for the discussions.}

In two dimensions, we did not obtain any new possibilities to gauge the entire Weyl symmetry than demanding $c=0$ from the beginning. We are still able to gauge the harmonic Weyl symmetry, but the physical interest in such gauging (e.g. whether it defines new class of quantum gravity in two dimensions) should be discussed more in detail.

Finally, we point out that there is an alternative option. Once we know how to solve $g(x)$ to cancel the Weyl anomaly in a given metric $g_{\mu\nu}(x)$, one may introduce the extra transformation on $g(x)$ so that the Weyl anomaly is always cancelled (irrespective of the obstructions we have discussed above). This possibility requires further investigation if such transformation can be defined systematically and then whether such a generalized notion of the Weyl transformation is useful or not.

\section*{Acknowledgements}
This work is in part supported by JSPS KAKENHI Grant Number 17K14301.


\section{Introduction}
          Ground state wavefunctions of the vast majority of commonly
          encountered Hamiltonians are characterized by the so-called `area-law' of
          entanglement~\cite{hastings2007area,eisert2010,laflorencie2016quantum}. The
          entanglement entropy of a subsystem with respect to its
          complement, scales not as the volume of the subsystem in
          question, but rather as the surface area that links the
          subsystem to its environment. This is loosely justified on
          the grounds that since the couplings are local (for the most
          extensively studied Hamiltonians), quantum correlations in
          the ground state are also local in nature and therefore the
          contributions to the entanglement entropy come from
          correlations at the surface alone. Gapless models show a
          $\log L$-correction to the law~\cite{wolf2006violation,lai2013violation,korepin2004universality}
          - correlations here are stronger than area law because such
          ground states are at a critical point and quantum
          fluctuations induce long-range correlations whereby a region
          deep inside the subsystem offers a non-vanishing
          contribution to correlations with a region far outside
          it. Such mild area-law violations are also fairly
          extensively studied and accepted to be a consequence of the
          criticality of the model. Stronger violations of the
          area-law have also been
          reported~\cite{shiba2014volume,vitagliano2010volume,pouranvari,gori2015explicit}.

	  Long-range couplings are ubiquitous in real physical
          systems, quantum and classical~\cite{campa2014physics,mukamel2009notes}.  A wave
          of current interest exists in uncovering the novel physics
          that emerges when interactions are made
          long-range~\cite{latella2017longrange,mori2013phase,levin2014nonequilibrium}.
          Although the majority of such work is on classical systems,
          there is indeed plenty of interest and work on quantum
          systems. An inexhaustive list includes frustrated
          magnets~\cite{kittel,sandvik2010ground}, spin
          glasses~\cite{young,laumann2014manybody} and various
          ultra-cold atomic~\cite{walker,rajibul,hazzard} and optical
          systems~\cite{sarang}.  One of the characteristics of
          long-range couplings is that, even one-dimensional models
          can give rise to higher-dimensional physics. In quantum
          models, one of the special consequences of this would be
          that by making the couplings to die sufficiently slowly,
          there ought to be stronger violations of the area-law than
          observed in gapless systems. With this hunch in mind, we
          make a detailed study of a long-range disordered hopping
          model in one dimension, where the strength of the couplings
          fall off with distance as a power-law with exponent
          $\sigma$.

                      In the power law model, by tuning the exponent
                      $\sigma$, we are able to discern three distinct phases:
                      one in which the ground state is delocalized and displays
                      a strong area law violation, a second intermediate phase
                      in which the ground state is quasi-localized and adheres
                      to the area law for large subsystem sizes, and a third
                      short range class where the ground state is localized and
                      subscribes to the area-law. The much studied
                      Aubry-Andre-Harper (AAH) model~\cite{aubry,harper} is
                      included for comparison and contrast. The AAH model has
                      the well-known self-dual structure which gives a
                      localization-delocalization transition, with the localized
                      phase being characterized by an area-law abiding
                      entanglement entropy. The quantum phase transition point
                      has the well-known $\log{L}$ correction to the area-law
                      entanglement entropy - we find that in fact, the entire
                      delocalized phase carries the $\log{L}$
                      correction.

                   To characterize the phases, we employ several tools
                   including inverse participation ratio (IPR), level
                   spacing ratios, entanglement entropy, subsystem
                   number fluctuations, and non-equilibrium
                   wave-packet dynamics keeping track of the spatial
                   distribution of the wave-packet. For free fermionic
                   models, entanglement entropy has been argued to be
                   closely connected to subsystem number
                   fluctuations\cite{klich2009,song2010general,song2011,HFSong2012,calabrese2012exact,flindt2015}. We
                   find evidence in support of this connection, both
                   in the statics and the dynamics that we study in
                   our model. In this context, we also
                     study a recently introduced quantity called
                     `entanglement contour' which quantifies the
                     contribution from each site in the subsystem to
                     the entanglement. The advantage of this microscopic quantification
                     is that features like area-law violation and central
                     charge of the system can be obtained from a
                     single subsystem calculation, without the need for any
                     subsystem scaling as with other quantifiers
                     of entanglement~\cite{vidal}. Also its
                     relation with `fluctuation contour' that
                     originates from the number-fluctuations in the
                     subsystem, is useful as a comparative tool~\cite{frerot}.
                     Entanglement contour nicely captures the the area law and its violation
                     in the disordered long-range hopping model. Also
                     the relationship between the two contours shows
                     striking behavior across the
                     delocalization-localization transition point.

                     Non-equilibrium dynamics of a closed quantum system has
                     become a topic of great interest in current
                     research~\cite{Eisert2015,Fagotti,Aditi}. Nowadays one of
                     the key perspectives for understanding different types of
                     phases is the study of entanglement propagation in many-body
                     systems.  This can be probed by tracking quasi-particles in
                     many cases~\cite{zoller2014,cardy2006}.  Also contrasting
                     behavior of various types of transport such as the transport
                     of charge, correlation and entanglement in quantum systems
                     is being used to characterize phases in many-body
                     systems. For example, both the Anderson localized and
                     many-body localized phases show no charge
                     transport~\cite{anderson,baa}; in contrast, the former shows no
                     growth of the bipartite entanglement entropy with time but
                     the latter shows a logarithmic growth~\cite{moore,abanin}.
                     Recently charge transport and
                     entanglement transport have been contrasted in
                     bond-disordered short-range models~\cite{bera}. We study nonequilibrium dynamics
                     in our bond-disordered long-range model, finding evidence for the contrast
                     between charge and entanglement propagation.
                     Another aspect of study of long-range models is
                     the generalization of Lieb-Robinson bounds which suggest that in
                     short-range models~\cite{lieb}, the velocity with which correlation spreads is bounded and
                     hence results in a light-cone like spreading of correlation. This leads to a linear
                     growth of entanglement entropy with time following a
                     sudden global quench in short-range models as predicted by
                     related CFT\cite{cardy2006}. The light-cone picture can break down
                     in long-range models; this has been seen theoretically and
                     experimentally in ultracold ion traps for translationally
                     invariant long-range
                     models\cite{Gorshkov,zoller2014,kastner,daley,wouters}. We
                     numerically test the break-down of the light-cone picture
                     in our disordered long-range model and find different
                     results in the delocalized, quasilocalized and localized
                     regimes, which we will discuss later.

                   The paper is organized as follows. In section II, we discuss
                   the delocalization-localization transition in the disordered
                   long-range hopping model.  In section III we explore the
                   entanglement of free fermions in the model at the
                   single-particle and many-particle levels. In subsection IIIA
                   we talk about the single-particle entanglement in the
                   model. In subsection IIIB we study entanglement of fermions
                   and its connection to the number fluctuations in the
                   subsystem. In subsection IIIC we implement the idea of the
                   entanglement and fluctuation contours. In subsection IIID we
                   compare our long-range model with the short-range AAH model. In
                   section IV we investigate the non-equilibrium dynamics
                   at the single-particle and many-particle levels and finally
                   we summarize in section V.

\section{Random long-range hopping model}
We consider a Hamiltonian of the following generic type:
\begin{equation}
\mathcal{H} =  \sum\limits_{i \neq j}^{N} (t_{ij} c_{i}^{\dagger}c_{j} + h.c.) +  \sum\limits_{i}^{N} v_i c_{i}^{\dagger}c_{i},
\label{hamiltonian}
\end{equation}
where $c_{i}^{\dagger} (c_{i})$ is the single fermion creation
(annihilation) operator at the $i$th site. In the long-range random
hopping model $t_{ij}=J\frac{u_{ij}}{{r_{ij}}^\sigma}$, is the strength
of hopping and $v_i=0$. $u_{ij}$ is chosen from $[-1,1]$, a
uniform distribution of random numbers and
$r_{ij}=(N/\pi)\sin(\pi|i-j|/N)$, is the geometric chord distance
between the $i^{th}$ and $j^{th}$ sites, when the sites are arranged
in a periodic ring. Here $J$, the maximum magnitude of the hopping
term, is the unit of energy, which we put to unity $J=1$. In a very
similar model \cite{lima1,lima2}, where $r_{ij}={|i-j|}$, $\sigma=1$
has been shown to be the delocalization-localization transition point,
in close connection with the power law random banded matrix (PRBM)
\cite{mirlin1,Cuevas,mirlin2,Evers} model. For $\sigma<1 (\sigma \geq 1)$,
all the eigenstates are delocalized (localized)~\cite{lima1}.
\begin{figure}[h!]
\centering
\includegraphics[width=4.55 cm,height=4.45 cm]{fig1a.jpg}~
\includegraphics[width=4.0 cm,height=3.8 cm]{fig1b.pdf}
\includegraphics[width=6.5 cm,height=5.0 cm]{fig1c.pdf}
\caption{(a) Surface plot of IPR of the single particle eigenstates as
  a function of $\sigma$ for system size $N=4096$ and $100$
  realizations of disorder. Here $\alpha$ stands for index of single
  particle eigenstates in ascending order of energy.  The color `green' lies right
  at the bottom of the bar legend, and is barely visible due
  to its very small value ($\sim 1/N$) - it corresponds to
  the $IPR$ of delocalized eigenstates. (b) Variation of $D_2$ with
  $N$ showing change of slope around $\sigma=1$. (c) The level-spacing
  ratio $r$ as a function of $\sigma$ for increasing system sizes $N$,
  averaged over $100$ realizations of disorder. The two dashed
  horizontal lines denote $r=0.529$ and $r=0.386$ respectively. Error
  bars are of the same order as the samples shown in the three
  regions; for other data points, error-bars are suppressed to enhance
  clarity. }
\label{transition}
\end{figure}

To quantify the point of the localization transition, we compute the inverse participation ratio (IPR), which is defined as
\begin{equation}
I(\alpha) = \sum\limits_{i=1}^{N} {|\psi_i(\alpha)|}^4,
\end{equation}
where the coefficients are drawn from the $\alpha^{th}$ normalized single particle eigenfunction $\ket{\psi(\alpha)}=\sum\limits_{i}^{} \psi_i(\alpha) \ket{i}$ expanded in the complete set of the Wannier basis $\ket{i}$, which represents the state of a single particle localized at the site $i$ of the lattice. IPR of all the eigenstates as a function of $\sigma$
is shown in the surface plot Fig.~\ref{transition}(a). We see the presence of
localized states at the edges of the band near $\sigma=1$, which is essentially a finite size effect \cite{lima1}.
\begin{figure*}
\centering
\includegraphics[width=5.7 cm,height=4.2 cm]{fig2a.pdf}
\includegraphics[width=5.7 cm,height=4.2 cm]{fig2b.pdf}
\includegraphics[width=5.7 cm,height=4.2 cm]{fig2c.pdf}
\caption{ Probability distribution $p_i(t)$ for finding a single particle, initially at the middle of the lattice, at each site of the lattice for increasing values of time (in units of $J^{-1}$) (a) for $\sigma=0.5$ (delocalized phase); (b) for $\sigma=1.5$ (quasilocalized phase); and (c) for $\sigma=3.0$ (localized phase) respectively. For all the plots $N=2048$ and number of disorder realizations is $100$. }
\label{wfn}
\end{figure*}

We also calculate the participation moments averaged over all the
eigenstates. The $q^{th}$ participation moment is obtained by averaging over all the eigenstates and disorder configurations
\begin{equation}
P_q=\bigg \langle \frac{\sum\limits_{\alpha=1}^{N} {P_q}(\alpha)}{N} \bigg \rangle,
\end{equation}
where ${P_q}(\alpha)=1/\sum\limits_{i=1}^{N}|{\psi_i(\alpha)|}^{2q}$. However, $P_q\propto N^{D_q(q-1)}$. In a fully delocalized (localized) regime $D_q$ approaches unity (zero) as the thermodynamic limit is approached. It is evident that $\frac{\log P_q}{\log N} \propto D_q (q-1)$ and from the variation of $D_q$ with the system size, one can identify the point of transition in the thermodynamic limit. We choose $q=2$ and $D_2$ is plotted with the system size in Fig.\ref{transition}(b). The $D_2$ vs $N$ plot changes slope at $\sigma=1$, which is the point of the localization transition.

The mean of the ratio~\cite{huse2007,Atas2013} $r$ between adjacent gaps ($\delta$) in the spectrum can be used to identify a crossover from Wigner-Dyson statistics in the delocalized phase to Poisson statistics in the localized phase. Defining
\begin{equation}
r_k = \frac{\min(\delta_k,\delta_{k+1})}{\max(\delta_k,\delta_{k+1})},
\end{equation}
where $\delta_k = \epsilon_{k+1}-\epsilon_{k}$ is the $k^{th}$ energy gap, the mean
ratio is $r = \langle \overline{r} \rangle$, where the bar represents
an average over the spectrum, and the angular brackets the average
over disorder. It is known from random matrix theory that the mean
ratio $r$ is approximately $0.529$ in the delocalized phase and
$0.386$ in the localized phase~\cite{huse2007,Atas2013}.
Fig.~\ref{transition}(c), based on the finite sizes considered here,
suggests that the system is in the ergodic phase in the region $0 \leq
\sigma \leq 1$.  Then $r$ starts decreasing till it reaches the
localized phase around $\sigma=2$.  The intermediate phase showing
intermediate distributions is discussed in the following analysis.

In order to better understand the presence of different phases in the system, we have considered a wavepacket initially localized at the middle site $i_0$ of the lattice i.e. $\psi_i(t=0)=\delta_{i,i_0}$ and calculated the evolution of the spatial distribution of the wavepacket with time. The probability of finding a particle at site $i$ at a given instant $t$ is given by $p_i(t)=|\psi_i(t)|^2$. The spatial dependence of the probability distribution for increasing time is shown in Fig.\ref{wfn}. It is to be noted that in the quasilocalized phase (Fig.\ref{wfn}(b)), the central part of the wavepacket rapidly drops down to a smaller value, which then barely changes with time whereas the tails of the wavepacket keep spreading with time. In the delocalized phase, the occupancy at the initial site along with all the other sites rapidly decreases and the wavepacket takes the form of a uniform distribution(Fig.\ref{wfn}(a)) whereas in the localized phase, the dynamics of the wavepacket is almost absent and it becomes almost exponentially localized (Fig.\ref{wfn}(c)). Fig.~\ref{wfn} thus shows that the quasi-localized phase is distinct from both the delocalized and localized phases, and yet carries some character of each of these phases.
\section{Entanglement in the model}
Phase transitions in extended quantum systems are
known to be captured by different measures of entanglement~\cite{horodecki2009quantum,amico2008entanglement,eisert2010} such as
concurrence, entanglement entropy etc. In the subsequent part of this
section, we will calculate the von Neumann entanglement entropy
between a suitable subsystem and its complement, both for single particle and many particle
states. We will investigate if there is a violation of the `area law'
of the entanglement entropy and analyze our results on the basis of
the localization transition. We will discuss local
particle-number fluctuations and its relation with entanglement
entropy in the context of the transition in our model. Also  we discuss the `entanglement contour' and `fluctuation contour' in this context.
\subsection{Single-particle entanglement}
First we discuss single-particle entanglement entropy, which has been argued to be a useful resource for quantum information processing~\cite{vanEnk2005,Dashenbrook2016}. In order to calculate the entanglement entropy between two subsystems A and B for the normalized single particle states, one writes down a normalized single particle eigenstate in the following way
\begin{equation}
\ket{\psi}=\sum\limits_{i\in A} \psi_i {c_i}^\dagger \ket{0}_A\otimes\ket{0}_B + \sum\limits_{i\in B} \psi_i\ket{0}_A\otimes{c_i}^\dagger\ket{0}_B,
\end{equation}
where $\ket{0}_{A/B}$ is the vacuum state in the subsystem A/B. Then the reduced density matrix $\rho^{sp}_A=Tr_B(\ket{\psi}\bra{\psi})$ has two eigenvalues $p_A=\sum\limits_{i \in A} |\psi_i|^2$ and $p_B=1 - p_A$\cite{Chakravarty}(see Appendix \ref{appA} for more details). The single particle entanglement entropy is then given by
\begin{equation}
S^{sp}_A=-p_A\ln p_A - p_B\ln p_B.
\end{equation}
\begin{figure}
\centering
\includegraphics[width=4.275 cm,height=4.0 cm]{fig3a.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig3b.pdf}
\caption{(a) Scaling of single-particle entanglement entropy $S^{sp}_A$ with subsystem size $L$  for different values of $\sigma$. (b) Variation of $S^{sp}_A$ with $\sigma$, where $L=N/2$. The system size $N=2048$ and number of disorder realizations is $100$ for both the plots.}
\label{singleEE}
\end{figure}
This entropy is bounded between $\ln2$ and $0$. In a delocalized
eigenstate, $S^{sp}_A$ increases with $L$, the size of the subsystem
A, as $p_A=L/N$ and reaches the maximum value $\ln2$, when $L=N/2$. In
a single site localized state $S^{sp}_A$ is $0$ as $p_A=1$ or $0$ and
does not show any variation with the subsystem size. The variation of
$S^{sp}_A$ with $L$ in different phases for our model is shown in
Fig.\ref{singleEE}(a). In the quasi-localized phase, $S^{sp}_A$ varies
with $L$ but its maximum value is less than $\ln2$ and the maximum
value decreases as $\sigma$ increases towards $\sigma=2$. The curves
deviate more from the delocalized ones as $L$ increases towards $N/2$
because in the quasilocalized eigenstate the central part of the
wavefunction is more localized compared to the tails. The variation of
$S^{sp}_A$ with $\sigma$ can be seen from Fig.\ref{singleEE}(b). The
delocalized $(\sigma<1)$, quasi-localized $(1<\sigma<2)$ and localized
$(\sigma>2)$ phases are clearly seen from the plot. Also it is worth
mentioning that the quasilocalized phase shows large intrinsic
fluctuations in $S^{sp}_A$. This results in large error-bars that
cannot be significantly reduced by increasing the number of disorder
realizations. This is obvious because in the quasi-localized phase,
for an eigenstate, the probability distribution for finding a single
particle has multiple peaks and they can appear in random places in
the lattice for different realizations of disorder (not shown here)
thus making $p_A$ a highly fluctuating quantity. In the localized
phase the probability distribution is more or less singly peaked hence
$p_A$ is always close to $0$ or $1$ whereas in the delocalized phase
the probability distribution has no peak and it is a uniform one,
hence $p_A\sim L/N$ giving rise to smaller error-bars in $S^{sp}_A$.

\subsection{Fermionic entanglement and fluctuations}
In this subsection we consider noninteracting spinless fermions at
half-filling in the system and investigate signatures of the
localization transition via entanglement in many-body states. The connection between localization and entanglement is subtle. Intuitively, one would expect that the greater the delocalization, the more the entanglement and vice versa; however, this correlation is not absolute and counterexamples are available\cite{kannawadi2016persistent}. We also discuss the
relationship between subsystem number fluctuations and entanglement entropy in the model. We start with a brief discussion of the calculation of the
entanglement entropy of fermions in the
ground state\cite{peschel2003calculation,peschel2009,peschel2012special}(see
Appendix \ref{appB} for details). For the fermionic many-body ground state
$\Ket{\Psi_0}$, the density matrix can be written as
$\rho=\Ket{\Psi_0} \Bra{\Psi_0}$. The entanglement
entropy between two subsystems is then given by $S_A=-Tr(\rho_A \log
\rho_A)$, where the reduced density matrix
$\rho_{A}=Tr_{B}(\rho)$. However, for a single Slater determinant
ground state, Wick's theorem can be exploited to write the reduced
density matrix as $\rho_{A}=\frac{e^{-H_{A}}}{Z}$, where
$H_{A}=\sum\limits_{ij} H_{ij}^A c_{i}^{\dagger}c_{j}$ is called the
entanglement Hamiltonian, and $Z$ is obtained from the condition $Tr
(\rho_{A}) = 1$. The information contained in the reduced density
matrix of size $2^L\times 2^{L}$ can be captured in terms of the
correlation matrix $C$ of size $L\times
L$\cite{peschel2003calculation} within the subsystem A, where
$C_{ij}=\left\langle c_{i}^{\dagger}c_{j} \right\rangle $. The
correlation matrix and the entanglement Hamiltonian are related
by\cite{peschel2003calculation,peschel2009,peschel2012special}:
\begin{equation}
  \label{eqn:C_and_H}
  C=\frac{1}{e^{H_A}+1}.
\end{equation}
Using this relation, the entanglement entropy for free fermions is given by\cite{peschel2009,peschel2012special},
\begin{equation}
S_A=-\sum\limits_{m=1}^{L} [\lambda_m \log \lambda_m + (1-\lambda_m) \log (1-\lambda_m)],
\end{equation}
where $\lambda_m$'s are the eigenvalues of the correlation matrix $C$.
\begin{figure}
\includegraphics[width=6.5 cm,height=5.5 cm]{fig4.pdf}
\caption{Participation ratio of MEM, denoted as $P^{mem}$, as a function of $\sigma$ for increasing system sizes $N$, averaged over $100$ disorder realizations. Subsystem size $L=N/2$ for fermions at half-filling.}
\label{mem}
\end{figure}
It has been conjectured that the zero mode of the entanglement
Hamiltonian has information about topological quantum phase
transitions~\cite{LiHaldane}. The same conjecture can be extended to a
non-topological system \cite{pouranvari}. It follows from
Eqn.~\ref{eqn:C_and_H} that the zero mode of the entanglement
Hamiltonian would correspond to the eigenfunction of the correlation
matrix, whose eigenvalue is equal (closest) to $0.5$ . As this
eigenmode contributes the maximum to the entanglement entropy, it is
called the maximally entangled mode (MEM). The participation ratio of
the MEM reflects the localization transition at $\sigma=1$
[Fig.~\ref{mem}]. This is a nice example of detecting the localization
transition from the entanglement spectra without having any prior
knowledge about the original Hamiltonian.

Now we will discuss the scaling of the entanglement entropy with
subsystem size. Typically, short range models of noninteracting
fermions show logarithmic violation of the area law of entanglement
entropy i.e. $S_A \sim L^{d-1} \log L$ in $d$ dimensions~\cite{swingle}.
In our disordered long-range model we see super-logarithmic
area law violation in the delocalized phase where $0<\sigma<1$. In fact it
goes as $L^\beta$, where the exponent $\beta=1$ at $\sigma=0$ and $\beta$ decreases as $\sigma$ increases
[Fig.~\ref{EE_scaling}](a). In the quasilocalized regime $1<\sigma<2$ it shows
area law for larger subsystem sizes whereas in the localized phase
$\sigma\geq2$ it shows a strict area law.

\begin{figure}
\centering
\includegraphics[width=6.5 cm,height=5.0 cm]{fig5a.pdf}
\includegraphics[width=6.5 cm,height=5.0 cm]{fig5b.pdf}
\caption{(a) A linear-log plot showing the scaling of the entanglement entropy $S_A$ with subsystem size $L$ for increasing $\sigma$ for fermions at half-filling. (b) Similar plot for the number fluctuation $\delta^2 N_A$. For both the plots lattice size $N=2048$ and number of disorder realizations is $100$.}
\label{EE_scaling}
\end{figure}

Next we discuss entanglement and its indirect experimental measurement. It has been argued~\cite{frerot} that fluctuations of a globally conserved quantity inside a subsystem can measure entanglement entropy as the quantity shares eigenfunctions with the reduced density matrix $\rho_A$ and hence provides a good basis for Schmidt decomposition of the many-particle eigenstate (see Ref.~\onlinecite{frerot} for a rigorous proof).
In our canonical set-up, total particle number is conserved and we study fluctuations in the particle number inside the subsystem, which is also an experimentally measurable quantity\cite{pitaevski,stringari}. The particle number fluctuations inside some subsystem $A$ can be defined as
\begin{equation}
\delta^{2}N_{A} = \sum_{i\in A}\langle n_{i}^{2}\rangle - \langle n_{i}\rangle^{2}.
\end{equation}
A close connection exists between entanglement entropy and fluctuations in the
local observables in the subsystem e.g. magnetization in a spin
system or particle number in free fermionic systems\cite{klich2009,song2010general,song2011,HFSong2012,calabrese2012exact,flindt2015}. The relationship
becomes a proportionality for certain gapless models, and the
proportionality constant to leading order has also been obtained~\cite{calabrese2012exact}.

We adopt this quantity to study our long-ranged model, and look at the
scaling of the particle number fluctuations with the subsystem
size. The number fluctuations in the subsystem can be calculated using
the following relation:
\begin{equation}
\delta^2 N_A=\sum\limits_{m=1}^{L} \lambda_m(1 - \lambda_m).
\end{equation}
Fig~\ref{EE_scaling}(b) reveals that this quantity shows exactly the same scaling as $S_A$, pointing to a proportionality between them, even in this long-range off-critical model. We will see that the \emph{proportionality constant} offers a signature for
the localization-delocalization transition in the model though. Likewise, the proportionality constant shows a sudden jump at the phase transition in the AAH model as well, as will be shown at the end of this section.

\subsection{Entanglement contour and fluctuation contour}
\label{sec:contour}
In this subsection, we will define and study the `entanglement contour'\cite{vidal}
and the `fluctuation contour'\cite{frerot}. These quantities contain microscopic
details of entanglement and number fluctuations. Specifically, the
contour keeps track of the contribution from each site within the
subsystem, to the quantity under consideration. Entanglement contour
is defined as the contribution ($C_s(i) \geq 0$) from the degrees of
freedom at each site $i$ in subsystem A to the entanglement entropy
$S_A$ such that $S_A = \sum\limits_{i \in A} C_s(i)$. One can calculate $C_s(i)$ using the
following relation\cite{vidal}:
\begin{equation}
C_s(i) = \sum\limits_{m=1}^{L} g_i(m) S_m,
\end{equation}
where $S_m = -[\lambda_m \log \lambda_m + (1 - \lambda_m) \log(1 - \lambda_m)]$. Here $\lambda_m$'s are the eigenvalues of the
correlation matrix $C$ or the entanglement spectra. $g_i(m)$ describes
the spatial pattern of the $m^{th}$ normalized eigenstate $\ket{\phi(m)}$ of matrix $C$ and hence of the entanglement Hamiltonian $H_A$
i.e. $g_i(m) = |\phi_i(m)|^2$. Similarly, one can define the
contour of subsystem particle-number fluctuations (also called as
`fluctuation contour') $C_n (i)=\langle \delta n_i \delta N_A
\rangle$, which is an obvious decomposition of the particle-number
fluctuations ($\delta^2 N_A$) in $A$ such that $\delta^2
N_A=\sum\limits_{i \in A} C_n (i)$. In the canonical ensemble $\delta
N_A=-\delta N_B$. Then $C_n (i) = - \langle \delta n_i \delta N_B \rangle$.
So one can interpret $C_n (i)$ as the correlation between
number (density) fluctuations at site $i$ and those in the whole of
subsystem $B$. It can also be defined as\cite{frerot}
\begin{equation}
C_n(i) = \sum\limits_{m=1}^{L} g_i(m) \lambda_m (1 - \lambda_m),
\end{equation}
where all the terms have the same meaning as defined previously.

\begin{figure}
\centering
\includegraphics[width=7.0 cm,height=5.5 cm]{fig6a.pdf}
\includegraphics[width=7.0 cm,height=5.5 cm]{fig6b.pdf}
\caption{(a) Spatial distribution of the scaled entanglement contour $C_{s0}(i)$ in the subsystem for different $\sigma$. (b) The same for the fluctuation contour $C_{n0}(i)$ for increasing values of $\sigma$.  For all the plots $N=2048$ and number of disorder realizations is $100$. Here $i$ is the site index in the subsystem and subsystem size $L=N/2$ for half-filled fermions.}
\label{contour}
\end{figure}

It turns out that for free fermionic systems $C_n (i)$ and $C_s (i)$
show similar spatial dependence \cite{frerot}.  Spatial dependence of the $C_{s0}(i)=C_s(i)/C_s(1)$
entanglement contour and the scaled fluctuation contour $C_{n0}(i)=C_n(i)/C_n(1)$ of the random long-range
hopping model are shown in Fig.\ref{contour}(a) and in
Fig.\ref{contour}(b) respectively. Since there are two boundaries
between two subsystems in a ring and because the entanglement and the
number fluctuations decay as one moves away from the boundaries,
contours are symmetric functions of sites with respect to the
midpoint of subsystem $A$. We fit this decay with the function
$1/x^\gamma$. Since the entanglement entropy is the sum of all the
contributions of the entanglement contour, one may guess that the
entanglement entropy dependence should be given by the integral $\int
\frac{1}{x^{\gamma}}dx$, which in turn suggests that the exponent
$\beta$ should be given by $\beta \approx 1-\gamma$. Indeed, we find evidence
for this[Fig.~\ref{compare}], deep in the delocalized phase.
\begin{figure}[h!]
\includegraphics[width=7.5 cm,height=5.5 cm]{fig7.pdf}
\caption{The exponent $\beta$ calculated from the subsystem scaling of the entanglement entropy and the other exponent $\gamma$ determined from the decay of the entanglement contour in the subsystem are compared. In this log-log plot $L^\beta$ (solid lines) and $L^{1-\gamma}$ (dashed lines) are plotted to establish the relation $\beta\approx1-\gamma$ in the delocalized phase $\sigma<1$.}
\label{compare}
\end{figure}

For a finer understanding of the entanglement contour at the boundaries
and in the bulk of the subsystem, the histogram of $C_s(i)$ is plotted
in Fig.\ref{histogram}. In the delocalized regime, the entanglement contour has
a finite value at all the sites and the histogram is a sharply peaked distribution whereas the distribution gets broadened
and the peak shifts towards $0$ as one approaches the point of
quasi-localization $\sigma=1$ [Fig.\ref{histogram}(a)]. In the quasilocalized
regime the entanglement contour deep in the bulk starts vanishing [Fig.\ref{histogram}(b)], which explains the validity of
the area law for larger subsystem size. In the localized regime the
entanglement contour almost vanishes in the whole bulk region and one
gets a strict area law in this regime. This is also evident from the
histogram for $\sigma=2.0$ in Fig.\ref{histogram}(b) which shows a sharp peak at $0$ with
almost no broadening.
\begin{figure}
\centering
\includegraphics[width=4.2 cm,height=4.0 cm]{fig8a.pdf}
\includegraphics[width=4.3 cm,height=4.0 cm]{fig8b.pdf}
\caption{ Histogram of the scaled entanglement contour $C_{s0}(i)$ (a) in the delocalized phase ($\sigma<1$); (b) in (quasi)localized phase ($\sigma>1$) respectively. For both the plots spinless fermions at half-filling are considered in a system of size $N=2048$. Here i indicates sites in the subsystem, whose size $L=N/2$ and number of disorder realizations is $100$.}
\label{histogram}
\end{figure}
The fluctuation contour also shows similar behavior as the entanglement contour (not shown here).

Since the entanglement entropy and local number fluctuations are
intimately related, it is useful to study this relationship at a
microscopic level by calculating the ratio of the two contours of the
related quantities i.e. $K(i)=C_s(i)/C_n(i)$. This ratio for
increasing values of $\sigma$ in the delocalized phase is shown in
Fig.\ref{ratio1}(a). It reveals a uniform proportionality between the
two contours in the deep delocalized regime. The proportionality
becomes non-uniform as $\sigma$ approaches the transition point
$\sigma_c=1$. In the (quasi)localized regime this non-uniformity
becomes so much worse that we omit these data in the interest of
clarity. A histogram in Fig.\ref{ratio1}(b) shows a peaked
distribution of $K(i)$ for smaller $\sigma$ and the distribution gets
broadened with almost vanishing peak for larger $\sigma$.
\begin{figure}
\centering
\includegraphics[width=4.275cm,height=4.0cm]{fig9a.pdf}
\includegraphics[width=4.275cm,height=4.0cm]{fig9b.pdf}
\caption{(a) Spatial distribution of the ratio of two kinds of contour $K(i)$ of half-filled fermions for increasing $\sigma$ in the delocalized phase. (b) The corresponding histogram of $K(i)$. Here lattice size $N=2048$ and $i$ is the site index within subsystem $L=N/2$ and number of disorder realizations is $100$.}
\label{ratio1}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=4.275 cm,height=4.0 cm]{fig10a.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig10b.pdf}
\caption{(a) Variation of the entanglement entropy $S_A$ and number fluctuation in the subsystem $\delta^2 N_A$ of half-filled fermions with $\sigma$. Here the system size $N=2048$. (b) The ratio of the two quantities $K$ as a function of $\sigma$ for different system sizes $N$. In both the plots subsystem size $L=N/2$ and number of disorder realizations is $100$.}
\label{ratio2}
\end{figure}

Next we study the proportionality constant $K$ of the relationship
$S_A = K \delta^2 N_A$ for free fermionic models. In a gapless system,
$S_A \propto \log L$ and $K=\pi^2/3$ for a 1D Fermi gas as shown in a
recent article\cite{calabrese2012exact}. However, in a gapped system
$K$ is not known in general; furthermore, $K$ is not believed to be a
universal quantity. This motivates us to investigate $K$ in our
long-range model. Though entanglement entropy and number fluctuations
in the subsystem vary in a similar fashion with $\sigma$
[Fig.\ref{ratio2}(a)], near the transition they conspire in such a way
that the ratio of them leaves a signature for the transition in the
model [Fig.\ref{ratio2}(b)]. The proportionality constant $K$ shows a
maximum at $\sigma_c=1$ and becomes almost constant in the localized
phase $(\sigma>2)$. Large error bars in the (quasi)localized regime in
Fig.\ref{ratio2}(b) are a reflection of the largely broadened
distribution of $K(i)$ in the same regime.
\subsection{Comparison with AAH model}
In the following, we have done a similar study as above in the AAH model which is a
short-range model that shows a sharp localization-delocalization
transition at finite disorder.
\begin{figure}
\centering
\includegraphics[width=4.275 cm,height=4.0 cm]{fig11a.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig11b.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig11c.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig11d.pdf}
\caption{Results from the 1D AAH model. (a) Scaling of the entanglement entropy $S_A$ with the subsystem size $L$ for increasing $\lambda$. The x axis is plotted in log-scale. (b) The ratio of entanglement contour to fluctuation contour $K(i)$ in the subsystem for different $\lambda$. (c) The corresponding histogram of $K(i)$. For the plots((a)-(c) $N=610$ (d) Proportionality constant $K$ as a function of $\lambda$ for increasing $N$. The inset is a fit to the $\Delta K_c$ vs $1/N$ data points(black square), where $\Delta K_c=K(\lambda=2.025)-K(\lambda=1.975)$. The red curve represents the fitting curve given by $\Delta K_c=1.12-5.47/N^{0.51}$. For all the plots, the subsystem size $L=N/2$ (except for figure (a)) for fermions at half-filling.}
\label{harper}
\end{figure}
The AAH model can be described by a
Hamiltonian of the same form as Eqn.\ref{hamiltonian} where
$t_{ij}=\delta_{i,j+1}$ and $v_i=\lambda \cos(2 \pi \eta
i)$. Here $\eta$ is a `Diophantine number' (e.g. $\frac{\sqrt{5} -
  1}{2}$, inverse of the `golden mean') and $\lambda$ is the strength
of the quasi-periodic disorder\cite{aubry,harper}. All the single particle eigenstates get
localized at $\lambda_c=2$\cite{mathieu}.

Our results for the Harper model are summarized in Fig.\ref{harper}.
In the tight binding model without any disorder $S_A \sim \log L$. As
the quasi-periodic disorder is turned on, in the delocalized phase
($\lambda<2$) $S_A$ retains the factor of $\log L$ which is a
modulated area-law behavior and in the localized phase ($\lambda \geq
2$) $S_A$ shows a strict area law, as shown in Fig.\ref{harper}(a).
In the delocalized regime $K(i)$ is close to $\pi^2/3$ in the bulk
whereas as one enters the localized regime it is no more a constant
and starts fluctuating [Fig.\ref{harper}(b)].  This is also evident
from the histogram of the same quantity. The distribution gets
broadened and the peak almost disappears in the localized phase
Fig.\ref{harper}(c). Also $K$ shows a jump at the transition point
$\lambda_c=2$ [Fig.\ref{harper}(d)]. We define $\Delta
  K_c=K_{\lambda=2+\delta \lambda}-K_{\lambda=2-\delta \lambda}$ near
  the quantum critical point $\lambda_c=2$, the scaling of which with
  the system size $N$ is well fitted by the functional form
  $\Delta K_c=1.12-5.47/N^{0.51}$ for $\delta \lambda=0.025$ (the inset of
  Fig.\ref{harper}(d)). As $N\rightarrow \infty$, $\Delta
  K_c=1.12$. So when $\delta \lambda\rightarrow 0$, $dK/d\lambda$ will
  diverge to $\infty$ at $\lambda_c=2$ and hence the $K$ vs $\lambda$
  plot will become vertical at $\lambda_c=2$ in the thermodynamic
  limit. The proportionality constant $K$ indeed captures transitions
in the system although it changes differently in the two models
studied here.
\section{Non-equilibrium dynamics}
Having studied the static quantities to analyze different phases, in
this section we investigate the dynamical properties of the model. A
non-equilibrium situation can be created by changing a parameter of
the Hamiltonian, locally or globally, through adiabatic or sudden
processes. Here we study the dynamics of entanglement entropy post a
sudden global quench in the bond-disordered long-range model and
compare the results with those of charge transport in the
system. We also briefly discuss correlation transport
  in the system in the context of velocity bounds on transport
  and the related light-cone picture~\cite{lieb}. We calculate the growth of
bipartite entanglement entropy $S_A(t)= -Tr(\rho_A(t)
\ln(\rho_A(t)))$, between two halves of the system A and B for our
model at half-filling. The data we present are with an initial state
of the density-wave(DW) type $\ket{\Psi}=\prod\limits_{i}
{c_{2i}}^\dagger \ket{0}$, which is evolved under the Hamiltonian at a
particular $\sigma$\cite{Fazio2006}. The DW state can be
  achieved by turning on an additional strong repulsive nearest neighbor interaction and
  then suddenly turning it off. We have checked that qualitatively similar results are
obtained when the initial state is the many-body ground state of
half-filled fermions corresponding to the Hamiltonian at $\sigma=2.5$,
with a quench carried out to various other values of $\sigma$.  To
calculate entanglement entropy, we use standard free fermion
techniques \cite{peschel2003calculation,Fazio2006}(see Appendix
\ref{appC} for details). Variation of $S_A(t)$ with time for the DW
type of initial state is shown in Fig.\ref{dw}(a). The entanglement
entropy varies with time in faster-than-linear fashion for $\sigma<1$
before it saturates, indicating the existence of a non-equilibrium
steady state. In the (quasi)localized regime ($\sigma>1$), after a
super-ballistic transient $S_A(t)$ goes in a sub-linear fashion with
time before reaching a saturating steady state. In the delocalized
phase, the saturation value $S_{A}^{\infty}$ barely changes with
$\sigma$; however in the quasi-localized phase, $S_{A}^{\infty}$
decreases with increasing $\sigma$. In the localized phase
($\sigma>2$) the entanglement growth becomes substantially suppressed as compared to the corresponding translationally
  invariant nearest neighbor model~\cite{cardy2005}, where the entanglement entropy
  reaches the saturation at a time $t_{sat}\sim L/2$. Also the
saturation values of $S_A$ in the localized phase are negligibly
small. The number fluctuations $\delta^2 N_A$, which are essentially
density-density correlations, reveal similar dynamics as $S_A(t)$
[Fig.\ref{dw}(b)].

\begin{figure}
\centering
\includegraphics[width=4.275 cm,height=4.0 cm]{fig12a.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig12b.pdf}
\caption{(a) Quench dynamics of the entanglement entropy $S_A$ of half-filled fermions with time (in units of $J^{-1}$) for increasing $\sigma$ from an initial DW type state. (b) Similar plot for subsystem number fluctuations $\delta^2 N_A$. For both the log-log plots $N=2048$, $L=N/2$ and number of disorder realizations is $100$. The thick line segment shows linear dependence on time for comparison.}
\label{dw}
\end{figure}
\begin{figure*}
\stackunder{\includegraphics[width=5.5 cm,height=4.5 cm]{fig13a.jpg}}{(a)}
\stackunder{\includegraphics[width=5.5 cm,height=4.5 cm]{fig13b.jpg}}{(b)}
\stackunder{\includegraphics[width=5.5 cm,height=4.5 cm]{fig13c.jpg}}{(c)}
\caption{(a)-(c) Surface plot showing the spreading of the amount of
  correlation $\log_{10}|C_d(t)|$ in time (in units of $J^{-1}$) and
  the lattice distance $d$ for $\sigma = 0.8, 1.3, 2.1$
  respectively. Here $d=|i-j|$, is the lattice distance between sites
  $i$ and $j$. The colors represent different ranges of values of the
  correlation, indicated in the bar legend attached to each
  figure - the numbers mentioned in the bar legend are the powers to be which ten is raised.
  For all the plots half-filled fermions are studied with $N=2048$, $L=N/2$ and number of disorder realizations is $100$. The initial state is the DW state,
  described in the text.}
\label{corr}
\end{figure*}

  In short-range models with translational invariance,
  following a global quench correlation transport happens with a
  constant velocity, defined as the Lieb-Robinson bound\cite{lieb},
  giving rise to a sharp causal light-cone like view of the
  correlation transport in space-time, outside of which correlations
  are exponentially suppressed\cite{cardy2006}. This leads to
  linear growth of entanglement entropy in such
  models\cite{cardy2005}. Breaking of translation invariance in
  short-range models can give rise to a much slower light cone, e.g. a
  logarithmic light-cone in the Anderson-localized phase\cite{burrell}
  and hence the entanglement entropy also shows a slow growth.  More
  than a linear growth of the entanglement entropy with time indicates
  the violation of the picture based on Lieb-Robinson bounds,
  which also bound the rate of growth of the entanglement. This kind
  of violation has actually been seen very recently in ultracold
  ionic experiments with translationally invariant long-range
  interacting spin models\cite{Gorshkov,zoller2014}. Also theoretical
  investigations have been carried out for translationally invariant
  long-range free fermionic models in this
  context\cite{kastner,daley,wouters}. To test the validity of the
  light-cone picture for correlation transport in our long-range
  free fermionic model with disordered hopping, we calculate the
  two-point correlation function $C_{d}=\langle{c_i}^\dagger
  c_j\rangle$ as function of time ($tJ$) and distance ($d=|i-j|$)
  between the sites $i$ and $j$ inside the subsystem as depicted in
  the surface plot in Fig.~\ref{corr}. At time $tJ=0$ the correlation
  matrix is diagonal with zero off-diagonal elements due to the
  product state structure of the initial DW state and the entanglement
  entropy is zero. At later times, different sites at distance
  $d=|i-j|$ start getting correlated. The correlation transport is
  more than linear or super-ballistic in nature within very short
  time-scales $tJ\sim 1$, which shows up as a transient in the
  quasi-localized $(1<\sigma<2)$ and localized $(\sigma>2)$ phases
  whereas in the delocalized phase, super-ballistic part is
  predominant as the time-scale for entanglement growth till it
  reaches the saturation is shorter ($t_{sat}J\sim 1$) in this
  case. This explains the super-ballistic entanglement growth in the
  system and violation of the picture based on Lieb-Robinson bounds.
  However later time dynamics of the correlation reveals
  different behaviors of the light-cone picture in three different
  phases as we detail it in the following. As we can see from
  Fig.~\ref{corr}(a) for $\sigma=0.8$, one can still perceive
  sub-linear light-cones in the delocalized regime. Sub-linearity
  indicates a decreasing velocity of the correlation transport with
  time as opposed to a constant velocity in the linear light-cone
  picture. The light-cone becomes more prominent and more sub-linear in
  the quasi-localized regime as can be seen in Fig.\ref{corr}(b). Very
  sharp sub-linear light-cones are visible in the localized regime
  (Fig.\ref{corr}(c)), where velocities of the correlation transport
  depend on the threshold values of the correlation. Sub-linearity of
  light-cones is more in this regime and hence the growth of entanglement
  entropy is very slow in the same regime. Such a change in light-cone picture from
  less prominent to more prominent can be seen in the three regimes
  $(\sigma<1, 1<\sigma<2$ and $\sigma>2)$ also for the corresponding
  translationally invariant long-range hopping model with initial DW
  state\cite{kastner}. In contrast to our model though, in the non-disordered model,
  all the light cones look linear and the related velocity bounds on the correlation
  transport decrease as $\sigma$ decreases.

Next we will compare entanglement transport with charge transport in
the system at the single-particle and many-particle
levels. Single-particle entanglement entropy $S^{sp}_A$ is calculated
by choosing the subsystem $A$ such that it continues to be half the
size of the total system, but it is now taken to be centered around the
initial localized wavepacket in the middle of the lattice. The dynamics of $S^{sp}_A$ reveals
super-ballistic nature in the delocalized phase but in the
(quasi)localized phase the initial super-ballistic behavior is
followed by a ballistic part before saturation
[Fig.~\ref{single}(a)]. Also the average width of the initially
localized wavepacket is calculated, which is defined as
\begin{equation}
w^{sp}(t) = \sqrt{ \sum\limits_{i} (i-i_0)^2 p_i(t)},
\end{equation}
where $p_i(t)=|\psi_i(t)|^2$ as also mentioned earlier and $i_0$ is the center of the lattice. The dynamics of the width in
different phases is shown in Fig.\ref{single}(b).  In the (quasi)localized phase, after a ballistic
transient, $w^{sp}(t)$ goes sub-linearly before it reaches a saturation
value and the exponent of the sub-linear variation decreases as
$\sigma$ increases. This signals a sharp
contrast between charge transport and entanglement dynamics
even within the single-particle picture. Although both the quantities reach saturation at the same time, the saturation values decrease abruptly with $\sigma$ in the quasi-localized phase and becomes vanishingly small in the localized phase.
\begin{figure}
\centering
\includegraphics[width=4.275 cm,height=4.0 cm]{fig14a.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig14b.pdf}
\caption{(a) The single-particle entropy $S^{sp}_A$ as a function of time (in units of $J^{-1}$) for increasing $\sigma$. (b) Similar plot for the width $w^{sp}$ of the single particle wavepacket. The thick solid line shows the linear dependence on time for comparison. For all the log-log plots $N=2048$, $L=N/2$ and number of disorder realizations is $100$.}
\label{single}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=4.275 cm,height=4.0 cm]{fig15a.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig15b.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig15c.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig15d.pdf}
\caption{(a)-(b) Variation of width $w$ of the many-particle wavepacket and the occupation density of the initially occupied sites $n_{loc}$ respectively with time (in units of $J^{-1}$) for increasing $\sigma$. In the delocalized phase $n_{loc}$ is just the filling fraction in the long-time limit. (c)-(d) Variation of the entanglement entropy $S_A$ and number fluctuations in the subsystem $\delta^2 N_A$ respectively with time(in units of $J^{-1}$) for increasing $\sigma$. The thick solid line shows linear dependence on time for comparison. For all the log-log plots $N=2048$, $L=N/2$ and number of disorder realizations is $100$ for fermions with filling fraction $0.1$ .}
\label{many}
\end{figure}

We also study the expansion dynamics of a cloud of fermions of a given filling and initial state in which fermions sit around the center of the lattice. This type of initial state can be prepared by switching on a trap potential and suddenly switching it off to study the evolution of the system under the quenched Hamiltonian. We calculate the expansion of the width of the many-particle cloud, which can be quantified by\cite{schneider2012}
\begin{equation}
w(t) = \sqrt{\frac{1}{N_p} \sum\limits_{i} (i-i_0)^2 \langle n_i(t)\rangle - \frac{1}{N_p} \sum\limits_{i} (i-i_0)^2 \langle n_i(0)\rangle},
\end{equation}
where $N_p$ is the total number of particles and $\langle n_i\rangle$ is the average occupation at site $i$ whereas $i_0$ is the center of the lattice. Simultaneously another quantity $n_{loc}$, which is the sum of the occupation densities at the initially occupied sites, is also investigated as a function of time. This quantity is defined as\cite{masud}:
\begin{equation}
n_{loc}(t)= \frac{1}{N_p}\sum\limits_{i}^{in. occ.} \langle n_i(t)\rangle.
\end{equation}
The width of the many-particle wavepacket $w$ in different phases is
shown in Fig.\ref{many}(a) and it shows the same qualitative feature
as $w^{sp}$. The variation of $n_{loc}$ with time nicely matches with
the dynamics of $w$ [Fig.\ref{many}(b)]. It decreases rapidly to the
saturation value, which is filling fraction in the delocalized phase
and barely changes in the localized phase. In the
  quasi-localized phase, it saturates to an intermediate value, which
  increases abruptly as $\sigma$ increases in the same phase. Also
we calculate the entanglement entropy for the same initially localized
many-particle state by choosing a subsystem of $L=\frac{N}{2}$
consecutive sites, whose center coincides with the center of the
lattice. It shows the same qualitative feature as
$S^{sp}_A$[compare Fig.~\ref{single}(a) and Fig.\ref{many}(c)]. Therefore,
similar to the single-particle picture, there is a contrast between charge transport and entanglement
propagation in the many-particle picture. The number fluctuations also
show similar dependence on time but it is smoother than $S_A$
[Fig.\ref{many}(d)]. The roughness of $S_A$ and $S^{sp}_A$ may be an
artifact to the special choice of the subsystem. This whole analysis
has been carried out at a filling of $0.1$; however, we have verified
that there is no qualitative dependence of these results on the
filling fraction since there is no mobility edge in the energy
spectra.
\begin{figure}[h!]
\centering
\includegraphics[width=4.275 cm,height=4.0 cm]{fig16a.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig16b.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig16c.pdf}
\includegraphics[width=4.275 cm,height=4.0 cm]{fig16d.pdf}
\caption{(a)-(b) Scaling of the saturation values of the entanglement entropy $S_A^{\infty}$ and width of the many-particle wave-packet $w^{\infty}$ respectively with the system size $N$. (c)-(d) Variation of $S_A^{\infty}$ and $w^{\infty}$ respectively with $\sigma$. For Fig.(c),(d) $N=2048$. For all the plots $L=N/2$ and number of disorder realizations is $100$ for fermions with filling fraction $0.1$ .}
\label{inf}
\end{figure}
The saturation values of the many-particle entanglement entropy and
the width of the wavepacket show similar variation with the system
sizes $N$ [Fig.\ref{inf}(a-b)]. In the delocalized phase both the
quantities go linearly with $N$ whereas in the quasi-localized phase
the dependence is sub-linear and they become almost independent of $N$
in the localized phase. This is quite expected as it
  reflects the sensitivity of the three phases to the boundaries of
  the system. The variation of these two quantities with $\sigma$ is
  shown in [Fig.\ref{inf}(c-d)] and they show similar dependences. In
  the delocalized phase both the quantities have almost constant and
  very high values whereas in the quasi-localized phase their values
  decrease abruptly with $\sigma$ and for large $\sigma$ in the localized phase, become tiny  and
  almost $\sigma$-independent.

\section{Conclusion}
To summarize, in this paper we study many static and dynamical
quantities to investigate the link between the
delocalization-localization transition and entanglement of spinless
fermions in a random long-range hopping model. Within the system sizes
used for numerical analysis, the system shows a delocalized phase for
$\sigma<1$ and a localized phase for $\sigma>2$. One also obtains a
quasi-localized phase for $1<\sigma<2$, as reflected by the
level-spacing ratio and wave-packet dynamics, but this phase may
vanish in the thermodynamic limit as hinted in the plots of
level-spacing ratio for different system sizes. Scaling of the
entanglement entropy with subsystem size reveals strong area-law
violation in the delocalized phase whereas the (quasi)localized phase
seems to adhere (for larger subsystems) strictly to the area law. In addition
to the eigenvalues of the entanglement Hamiltonian, the maximally
entangled mode or the zero mode of the entanglement Hamiltonian,
also captures the localization transition, despite it being a non-topological
system. The entanglement contour, which is constructed out of both the
eigenvalues and the eigenfunctions of the entanglement Hamiltonian,
gives a picture of the spatial distribution of entanglement inside the
subsystem and nicely explains the violation of the area-law in the
system. Particle-number fluctuations in the subsystem have similar
dependence on space and time as the entanglement entropy. The ratio of
these two quantities shows a sharp signature at the point of the
localization transition. However, the nature of this signature is dependent
on the model in question as it is different in the AAH model from our
long-range model. The distribution of the ratio of the entanglement
contour to the fluctuation contour is sharply peaked in the
delocalized phase but the peak starts vanishing as one goes into the
(quasi)localized phase.

Also we study quench dynamics and wave-packet dynamics of fermions at
the single-particle and many-particle levels. At both the levels the
entanglement propagation and the charge transport show a sharp
contrast. Entanglement entropy shows super-ballistic behavior both in the
delocalized phase and the (quasi)localized phase, although this appears
only as a transient in the latter. This super-ballistic behavior is
attributed to the picture based on the Lieb-Robinson bounds for the spreading of correlation post a global quench. Contrastingly, the width of the wave-packet varies ballistically with time in the delocalized phase
while in the (quasi)localized phase after ballistic transient it
shows a sub-ballistic behavior with time before it saturates. In a
short-range model with disorder, the light cone picture is valid, and therefore
the time dependence of entanglement entropy is always sub-ballistic in general. However, in our model
long-range couplings give rise to super-ballistic behavior.
The saturation values of the width and entanglement entropy
show similar dependence as a function of the system size and $\sigma$
reflecting the presence of three phases in finite systems.

In our study, we have been able to explain the strong area law
violation in our long-range model by implementing the idea of
entanglement contour and connect them to the
delocalization-localization transition in the system by studying
quench and wave-packet dynamics. We hope that our results regarding
the relationship between entanglement entropy and number fluctuations
will help boost the possibility of indirect measurement of
entanglement in experiments. Also we have shown explicitly the
contrast between charge and entanglement transport, which is one of
the current topics of interest. As a future possibility, one can also
look for many-body localized phases in an interacting version of this
model. We hope that our work can trigger experimental studies of the disordered long-range model
in ongoing ionic trap experiments.

\section*{Acknowledgements}
We are grateful to the High Performance Computing (HPC) facility at IISER Bhopal, where large-scale calculations in this project were run.
A.S is grateful to Simone Paganelli and Andrea Trombettoni for helpful discussions, and to SERB for the startup grant (File Number: YSS/2015/001696).
N.R acknowledges Sourin Das for bringing to his attention useful references, and CSIR-UGC, India for his Ph.D fellowship.