\section{\label{Intro}Introduction}

The simulation of two-dimensional (2D) materials and prediction of their properties has become a mainstay of materials science over the past decade, with the promise and realization of valuable applications in both industrial technology and fundamental research~\cite{geim_van_2013}. Theoretical and computational methods for 2D materials have been advanced into a sophisticated machinery that enables researchers to deal with ever more realistic settings~\cite{rodriguez-vega_ground_2014}. The widely used Kohn-Sham density functional theory (KS-DFT)~\cite{KohnSham1965,Dreizler1990} presents one particularly popular ab initio approach with the capability of accurately handling hundreds of interacting particles (and up to thousands of atoms in cases where linear-scaling methods apply \cite{Goedecker1999,Soler2002,Bowler2012,Fox2014,Aarons2016,Cole2016}). The development of functionals for 2D systems has been lagging behind that of their 3D counterparts for various reasons: Some of the most heavily used 3D KS-DFT functionals are not bounded from below \cite{Chiodo2012,Kaplan2018} when the 2D limit is approached. This stems from the improper scaling behaviour of the density \cite{Chiodo2012}. Furthermore, even consistent first-order gradient corrections of the kinetic energy density (used in developing meta-generalized gradient approximations for KS-DFT) were unknown for 2D fermion systems until recently \cite{Holas1991,Brack2003,Trappe2016}. Nonetheless, considerable progress has been made and alternative derivations of some of those 2D functionals have been obtained since the turn of the millenium \cite{Pollack2000,Reimann2002,Pittalis2007,Pittalis2008,Constantin2008,Pittalis2009,Pittalis2009b,Pittalis2010,Raesaenen2010b,Vilhena2014}.

However, a systematic ab inito methodology that is universally applicable and scales favorably with particle number, thereby enabling high-throughput computations of mesoscopic systems, is not yet available --- with orbital-free density functional theory (OF-DFT) being the suspected saviour for almost a century~\cite{Thomas1927,Fermi1927,Gross2011,Burke2012,Karasiev2013,GrossBurke2015,Witt2018}. While KS-DFT scales cubically with particle number in generic settings, OF-DFT scales linearly, and sub-linear scaling can be achieved in special cases. Functional development, in particular concerning the kinetic energy functional \cite{Wang1992,DellaSala2015,Xia2015}, and implementations of OF-DFT have gained momentum in recent years \cite{Karasiev2012b,Karasiev2013,Chen2015,Das2015,Chen2016,Mi2016,Constantin2018,Witt2018}, also in conjunction with other techniques like ab initio molecular dymanics \cite{Gonzalez2008}. Chemical accuracy is approached in selected cases \cite{Xia2012,Borgoo2014,EspinosaLeal2015}. If quantum effects play a minor role or if the considered system is largely homogeneous, OF-DFT can also be used in its most basic form, the Thomas--Fermi (TF) approximation. For instance, the effect of exchange on large disordered systems with long-range interactions was studied in Ref.~[\onlinecite{Rossi2008}] using OF-DFT in TF approximation. The TF model is not only of historical significance, but presents, as an exact constraint for homogeneous systems and in the limit of infinite nuclear charges \cite{Lieb1977}, an important base line for benchmarking proposed systematic density functional improvements. To what extent then do corrections to the TF approach play a crucial role or dominate over exchange effects in 2D materials? (They do, indeed, for a number of relevant fermionic systems, ranging from atomic Fermi gases to molecules and single atoms.)

The most severe obstacle for OF-DFT in taking over as the workhorse of theoretical chemistry and materials science is the lack of accurate, reliable, systematic, and preferably universal quantum corrections to the quasi-classical TF approximation, in particular for the kinetic energy of low-dimensional systems~\cite{Holas1991,Brack2003,Salasnich2007,Trappe2016,Trappe2017}. While ad-hoc corrections to the quasi-classical limit and heuristic approximations are available for kinetic energy and particle density of low-dimensional systems~\cite{vanZyl2013,vanZyl2014}, successful derivations of systematic and consistent corrections are scarce~\cite{Burke2015,TrDeMu2015,Trappe2016,Trappe2017,Chau2018}. One promising route towards systematic orbital-free quantum corrections is provided by density-potential functional theory (DPFT)~\cite{Berge1988,Berge1992,CinalBerge1993,Trappe2016,Trappe2017,Chau2018}, a more flexible reformulation of the original Hohenberg-Kohn DFT~\cite{HohenbergKohn1964,Dreizler1990}, which circumvents the need for an explicit kinetic-energy density-functional and provides natural ways for systematic semiclassical expansions.

In this article we explore the applicability of DPFT for 2D materials by assessing quantum-corrections to the TF approximation for double-layer heterostructures of mono- and bilayer graphene. Of particular interest to us are situations that are not easily tackled with orbital-based techniques, for example 2D material sheets of mesoscopic size that are subjected to aperiodic disorder potentials. Such situations are for example of current interest in studies on Coulomb drag~\cite{narozhny_coulomb_2016} where there exists an unsettled controversy as to whether the behavior of drag measured in experiment~\cite{gorbachev_strong_2012} is due to correlation~\cite{song_energy-driven_2012} or anti-correlation~\cite{ho_theory_2018} between the density fluctuations of the layers.

Our work contributes in several ways to answering some of the questions raised above. Sections \ref{DPFT} and \ref{SCsim} provide the computational framework for obtaining quantum-corrected carrier-densities of 2D materials using DPFT. The expressions for the semiclassical particle densities developed here enable us to decide whether or not the quasi-classical TF approximation is sufficient for describing at least conglomerate properties like average inter-layer correlations of heterostructures. Section~\ref{Bilayersetup} introduces the generic double-layer system, with both layers subjected to one layer of charged impurities, while only one of the layers is strained. Charge and strain disorder potentials are expected to compete in creating correlated (from charged impurities) and anticorrelated (from strain) carrier densities in the two layers. Our model setup is designed to extract the strain strengths required for switching between correlation and anticorrelation. In Sec.~\ref{Results} we apply our new approach to double-monolayer graphene and double-bilayer graphene. We discuss whether or not the electron-hole puddles of both layers, interacting electrostatically, require a self-consistent inter-layer treatment. Finally, we analyse the effects of quantum-corrections and exchange energy on the correlations with the aid of phase diagrams that chart the correlation measures as functions of impurity density, carrier density, and ratios of strain and charge disorder. The appendix gathers background information on the units, system parameters, correlation measures, and numerical procedures employed here.

\section{\label{DPFT}Density-potential functional theory}
Instead of resorting to the computationally demanding orbital-based Kohn-Sham DFT, we make use of orbital-free density-potential functional theory (DPFT)~\cite{Berge1988,Berge1992}. It is formally equivalent to the Hohenberg-Kohn formulation, but makes systematic improvements upon the TF approximation technically feasible --- in particular for low-dimensional systems.

Specifically, by Legendre-transforming $E_{\mathrm{kin}}[n]$, the kinetic energy functional of the particle density $n(\vec r)$, w.r.t.~the new variable ${V(\vec r)=\mu-\frac{\delta E_{\mathrm{kin}}[n]}{\delta n(\vec r)}}$, we recast the total energy of an interacting quantum system with interaction energy $E_{\mathrm{int}}[n]$,
\begin{align}\label{Energynmu}
E[n,\mu]&=E_{\mathrm{kin}}[n]+E_{\mathrm{ext}}[n]+E_{\mathrm{int}}[n]\nonumber\\
&\quad+\mu\left(N-\int(\d\vec r)\,n(\vec r)\right),
\end{align}
as the density-potential functional
\begin{align}\label{EnergyVnmu}
E[V,n,\mu]&=E_1[V-\mu]-\int(\d\vec r)\,n(\vec r)\,\big(V(\vec r)-V_{\mathrm{ext}}(\vec r)\big)\nonumber\\
&\quad+E_{\mathrm{int}}[n]+\mu N.
\end{align}
Here, the external potential $V_{\mathrm{ext}}(\vec r)$ yields the external energy $E_{\mathrm{ext}}[n]$, and the particle number $N$ is enforced via the Lagrange multiplier $\mu$, viz.~the chemical potential. From Eq.~(\ref{EnergyVnmu}) we obtain the ground-state solutions of the three variables $V$, $n$, and $\mu$ by self-consistently solving
\begin{align}
n(\vec r)&=\frac{\delta E_1[V-\mu]}{\delta V(\vec r)},\label{nasdef}\\
V(\vec r)&=V_{\mathrm{ext}}(\vec r)+\frac{\delta E_{\mathrm{int}}[n]}{\delta n(\vec r)}, \label{Vdef}\\
N&=\int(\d\vec r) n(\vec r) \label{Nmu}.
\end{align}
Equation~(\ref{Nmu}) is obtained by combining $\partial E[V,n,\mu]/\partial\mu$ with Eq.~(\ref{nasdef}) and reveals the particle number constraint in Eq.~(\ref{Energynmu}).

Equations~(\ref{nasdef})--(\ref{Nmu}) are exact and reminiscent of the KS scheme, but without the need of orbitals. However, the noninteracting case aside, we have to approximate the unknown potential functional ${E_1[V-\mu]}$; here the subscript indicates that $E_1$ can be written in terms of a single-particle trace over a function of the single-particle Hamilton operator \cite{Berge1992}. We also have to provide the equally important interaction energy $E_{\mathrm{int}}[n]$ as an explicit functional of the particle density $n$. Approximate particle densities follow directly from approximations of ${E_1[V-\mu]}$ (or, rather, its functional derivative) for any given potential $V$. As is evident from Eq.~(\ref{Vdef}), $V$ constitutes an effective single-particle potential with interaction effects effectively included for any given density $n$

Following Refs.~[\onlinecite{BGE1984,Berge1992,CinalBerge1993,Trappe2016,Trappe2017,Chau2018}] we approximate $E_1$ by its noninteracting version as the single-particle trace
\begin{align}\label{E1}
E_1[V-\mu]=\tr\{(H_1-\mu)\,\eta(\mu-H_1)\},
\end{align}
where ${H_1=H_1(\vec R,\vec P)=T(\vec P)+V(\vec R)}$ is a single-particle Hamiltonian with dispersion relation $T$ and potential energy $V$, while the trace includes the degeneracy factor $g$. For example, ${g=4}$ accounts for the spin and valley multiplicity of unpolarized charge carriers in the cases of mono- and bilayer graphene. $\vec R$ and $\vec P$ are the position and momentum operators, respectively, and $\eta(\,)$ denotes the step function.

The explicit expression for $E_1$ in Eq.~(\ref{E1}) results in an explicit expression for the particle density in terms of arbitrary functions $V(\vec r)$ via Eq.~(\ref{nasdef}). The approximate nature of Eq.~(\ref{E1}) aside, the exact particle density including all quantum corrections is thereby obtained for any specified interaction energy $E_{\mathrm{int}}[n]$ and without reference to orbitals. Specifically, Eqs.~(\ref{nasdef}) and (\ref{E1}), together with the Fourier transform of the step function, yield the particle density\footnote{The contour integration circumvents the singularity at ${t=0}$ in the lower half-plane.}
\begin{align}\label{nPropagator}
n(\vec{r})=g\Int\frac{\d t}{2\pi\I t}\,\Exp{\frac{\I t}{\hbar}\mu}\bra{\vec{r}}\Exp{-\frac{\I t}{\hbar}H_1}\ket{\vec{r}},
\end{align}
see Refs.~[\onlinecite{Golden1957b,Golden1960,Light1973,Lee1975,Chau2018}].

\begin{table*}[t]
\centering
\begin{tabular}{c|c|c}
$T(\vec p)$ & $n_{\mathrm{TF}}(\vec r)$ & $n_3(\vec r)$ \\
& & \\[-2.0ex]
\hline\hline\\[-2.0ex]
$v_{\mathrm{F}}|\vec p|$ & $\frac{g}{4\pi}K(\vec r)^2$ & $\frac{g}{4\pi^2}\int(\d\vec s)\frac{K(\vec r+\vec s)^3}{|\vec s|}\,\mathrm{J}_1\big(2|\vec s|\,K(\vec r+\vec s)\big)$ \\
& & \\[-2.0ex]
\hline\\[-2.0ex]
$\frac{\vec p^2}{2m}$ & $\frac{g}{4\pi}\sigma(\vec r)$ & $\frac{g}{4\pi^2}\int(\d\vec s)\frac{\sigma(\vec r+\vec s)}{|\vec s|^2}\,\mathrm{J}_2\big(2|\vec s|\sqrt{\sigma(\vec r+\vec s)}\,\big)$ \\
\end{tabular}
\caption{\label{TableDensities}Thomas--Fermi density $n_{\mathrm{TF}}(\vec r)$ and quantum-corrected density $n_3(\vec r)$ for linear and quadratic dispersion in 2D from evaluating Eq.~(\ref{nPropagator}), with degeneracy factor $g$, Bessel functions $\mathrm{J}_\alpha(\,)$, ${\nu(\vec z)=[\mu-V(\vec z)]_+}$, ${K(\vec z)=\nu(\vec z)/(\hbarv_{\mathrm{F}})}$, ${\sigma(\vec z)=2m\,\nu(\vec z)/\hbar^2}$, and $[x]_+$ denoting $x\,\eta(x)$. We recover the TF densities from their quantum-corrected successors upon replacing ${K(\vec r+\vec s)}$ and ${\sigma(\vec r+\vec s)}$ by their local versions $K(\vec r)$ and $\sigma(\vec r)$, respectively. The case of quadratic dispersion is dealt with in Ref.~[\onlinecite{Chau2018}]; see Appendix~\ref{Appendixn3} for the derivations in the case of linear dispersion.}
\end{table*}

We seek to approximate the time evolution operator ${U=\mathrm{e}^{-\frac{\I t}{\hbar}H_1}}$ systematically via split-operator methods, for example of the Suzuki-Trotter type~\cite{Hatano2005,Chau2018}. The quasi-classical approximation ${U\approx U_2=\mathrm{e}^{-\frac{\I t}{\hbar}T(\vec P)} \mathrm{e}^{-\frac{\I t}{\hbar}V(\vec R)}}$ recovers the quasi-classical TF density $n_{\mathrm{TF}}$, while ${U_3=\mathrm{e}^{-\frac{\I t}{2\hbar}T(\vec P)} \mathrm{e}^{-\frac{\I t}{\hbar}V(\vec R)}\mathrm{e}^{-\frac{\I t}{2\hbar}T(\vec P)}}$ produces the first quantum-corrected density $n_3$ in a series of expressions that utilize higher-order factorizations~\cite{Chau2018} of $U$. We give the corresponding 2D densities for linear and quadratic dispersion in Table~\ref{TableDensities} and outline the derivation of $n_3$ for the case of linear dispersion in Appendix~\ref{Appendixn3}. In contrast to $n_{\mathrm{TF}}(\vec r)$, which is restricted to classically allowed regions of the potential and only depends on the local value $V(\vec r)$, $n_3(\vec r)$ samples $V$ in an extended region and exhibits evanescent tails beyond the quantum-classical border.

\section{\label{SCsim}Self-consistent simulation of disordered 2D materials}

We target 2D systems with chemical potential in the vicinity of the Dirac point (the point where valence and conduction bands touch) for graphene (effective bilayer graphene). The usual tight binding approach absorbs the lattice structure in an effective Hamiltonian and yields noninteracting quasiparticles in a homogeneous (${V_{\mathrm{ext}}(\vec r)=0}$) environment. The electronic structures, viz.~atoms, of the materials are thus not modeled explicitly here. The single-particle energies, viz.~band structure, associated with such quasiparticles is the dispersion relation $T(\vec p)$ whose operator version appears in Eq.~(\ref{E1}). $T(\vec p)$ can be an arbitrary function, but for the purpose of this work we shall restrict ourselves to the analytically more tractable cases of linear and quadratic dispersion\footnote{Hence, the scope of our investigation also extends to energies further away from the Dirac point (band gap) as long as $T(\vec p)$ can be considered linear (quadratic).}.

In the following we outline the procedures involved for arriving at the ground state solutions of Eqs.~(\ref{nasdef})--(\ref{Nmu}); further details are provided in Appendix~\ref{AppendixSCloop}. Upon adding external potentials $V_\mathrm{C}(\vec r)$ and $V_\mathrm{S}(\vec r)$ that model charged impurites and strain, respectively, we initiate the self-consistent loop of Eqs.~(\ref{nasdef})--(\ref{Nmu}) by evaluating the density (denoted $n_-(\vec r)$ for the quasiparticles that follow the dispersion of the conduction band) with the external potential
\begin{align}\label{V-}
V_{\mathrm{ext}}(\vec r)=V_\mathrm{C}(\vec r)+V_\mathrm{S}(\vec r)
\end{align}
for these conduction quasiparticles, see Appendix~\ref{AppendixDisorder} for details. Since no interactions are included at this stage, the effective potential is ${V_-(\vec r)=V_{\mathrm{ext}}(\vec r)}$. The density $n_+(\vec r)$ of valence quasiparticles, which follow the inverted dispersion, e.g. $T(\vec p)=-v_{\mathrm{F}}|\vec p|$ in the case of graphene, is built from the same density expression as $n_-(\vec r)$ but takes as an input the inverted potential $2\,(\mu+\Delta)-V_{\mathrm{ext}}(\vec r)$, with an optional bandgap $\Delta$ (in case of mono- and bilayer graphene, we have ${\Delta=0}$). The effective potential for the valence quasiparticles reads
\begin{align}\label{V+}
V_+(\vec r)=2\,(\mu+\Delta)-V_-(\vec r)
\end{align}
and equals the external potential for the valence quasiparticles if interactions are omitted. The such obtained carrier density ${n(\vec r)=n_-(\vec r)-n_+(\vec r)}$ updates the effective potential via the interaction contribution in Eq.~(\ref{Vdef}).

As an approximate interaction energy $E_{\mathrm{int}}[n]$ for the quasiparticles, we employ the regularized Hartree term for the Coulomb energy,
\begin{align}\label{HartreeEnergy}
E_{\mathrm{H}}[n]=\frac{W}{2}\int(\d\vec r)(\d\vec r')\frac{n(\vec r)\,n(\vec r')}{\mathrm{max}(|\vec r-\vec r'|,b)},
\end{align}
where $b$ is half the lattice constant of the numerical implementation and $W=\hbarv_{\mathrm{F}} r_\mathrm{s}$, with the ratio $r_\mathrm{s}$ of Coulomb potential energy and kinetic energy for graphene. Upon functional differentiation, Eq.~(\ref{HartreeEnergy}) leads to the Hartree potential\footnote{Equation (\ref{HartreePotential}) can be regarded as a simplified version of the smooth models for the Coulomb potential discussed in Ref.~[\onlinecite{GonzalezEspinoza2016}].}
\begin{align}\label{HartreePotential}
V_{\mathrm{H}}(\vec r)=W\int(\d\vec r)\frac{n(\vec r')}{\mathrm{max}(|\vec r-\vec r'|,b)}
\end{align}
as an approximate interaction contribution in Eq.~(\ref{Vdef}). Eyeing means of comparison and higher accuracy, we may supplement $E_{\mathrm{H}}$ with an exchange energy~\cite{Dirac1930}, leading to the exchange potential $V_{\mathrm{X}}$; see Appendix~\ref{AppendixExchange} for details.

The updated effective potential $V$ then determines a new quasiparticle density $n_-(\vec r)$ via Eq.~(\ref{nasdef}), thereby closing the self-consistent loop. This process is repeated until a predefined relative precision is reached (we find $10^{-6}$ to be sufficient) when comparing the local densities of subsequent loop iterations. The chemical potential is adjusted in each iteration to enforce a given particle number, viz.~average carrier density. Figure~\ref{TFvsST3+X_monolayer} highlights the differences in the converged quasiparticle densities $n_{\mathrm{TF}}$ and $n_3$ of a single graphene layer with and without exchange. Both exchange and quantum corrections tend to decrease the peaks of the density landscape. This effect is well-known in the case of exchange~\cite{Rossi2008}, and a smoothening of the carrier density in an external potential $V_\mathrm{C}$ is to be expected when tunneling starts to play a role with the inclusion of quantum corrections. In fact, when comparing the corresponding densities in Fig.~\ref{TFvsST3+X_monolayer} we find that quantum corrections considered in this work can dominate over exchange effects.

\begin{figure}[htb!]
\begin{center}
\includegraphics[width=0.95\linewidth]{TFvsST3+X_monolayer_20190411}
\caption{\label{TFvsST3+X_monolayer} The spatial distributions of quasiparticle densities $n_{\mathrm{TF}}$ and $n_3$ of monolayer graphene visibly depend on exchange and quantum corrections. The upper left panel shows the Thomas--Fermi density $n_{\mathrm{TF}}$. The other three panels display the differences between $n_{\mathrm{TF}}$ and the quantum-corrected density $n_3$ (upper right panel), the TF density with the exchange potential of Eq.~(\ref{ExchangePotentialRossi}) included (lower left panel), and $n_3$ including exchange (lower right panel). The color bar refers to densities measured in units of ${0.001/l^2}$. The horizontal and vertical axes indicate spatial position, in the units of $l$ as labelled explicitly in the lower right panel. The same disorder realization, with average impurity density ${\bar{n}_{\mathrm{imp}} = 10^{12}/\mathrm{cm}^{2}}$, is used in all four panels. }
\end{center}
\end{figure}

The analysis of disorder averages reveals a striking instance of this observation. As observed in Fig.~\ref{TFvsST3+X_monolayer} for a single disorder realization, $n_3$ and $n_{\mathrm{TF}}$ differ in their density distribution function. For $n_3$ the integration over a finite region in the disorder landscape, see Table~\ref{TableDensities}, tends to result in smoother densities compared with $n_{\mathrm{TF}}$. This effect can be quantified by density histograms compiled from many disorder realizations. The density histograms in Fig.~\ref{DisorderAverages_Histogram}, calculated for a graphene monolayer on SiO$_2$ with ${\bar{n}_{\mathrm{imp}} = 10^{12}/\mathrm{cm}^{2}}$, corroborate the snapshot of one disorder realization in Fig.~\ref{TFvsST3+X_monolayer}. Figures~\ref{TFvsST3+X_monolayer} and \ref{DisorderAverages_Histogram} highlight the (a priori) importance of not only addressing exchange but also quantum corrections beyond the TF approximation for quantitatively viable investigations of 2D materials via orbital-free DFT: Both the inclusion of exchange and of quantum corrections indicate less pronounced peaks of the carrier-density landscape. Compared to the TF approximation, the quantum corrections exhibit an averaging effect with broader density distributions that have relatively more weight on intermediate values of the density rather than a strong distribution maximum at zero density; cf. Fig.~\ref{DisorderAverages_Histogram} (left). Although exchange also shows visibly less pronounced densities in Fig.~\ref{TFvsST3+X_monolayer}, this effect stems from a global reduction in density variance rather than a redistribution of densities from very small towards intermediate values. The distributions of the quantum-corrected densities (including exchange) are captured by their Gaussian (`G') or Lorentzian (`L') fits (with offset) more accurately than the distributions resulting from the TF approximation. This observation is in line with experimental results for graphene that point towards Gaussian distributions of density and energy in the presence of charged disorder~\cite{martin_observation_2008,xue_scanning_2011}. In Fig.~\ref{DisorderAverages_Histogram} we compare our calculations `$n_{\mathrm{TF}}$(+X)' and `$n_3$(+X)', with their Gaussian fits `G(TF)' and `G', directly with the experimental data `exp.' from scanning tunnelling spectroscopy; cf.~Ref.~[\onlinecite{xue_scanning_2011}]. We find our quantum-corrected approach to predict the experimental data much better than what can be obtained from the TF approximation --- in both qualitative and quantitative terms. In view of this stark improvement over the TF approximation, we want to stress again that our quantum-corrected density expressions are based on first principles without adjustable parameters or fits, and rely solely on controlled approximations to quantum mechanics.

\begin{widetext}

\begin{figure}[htb!]
\begin{center}
\includegraphics[width=0.405\linewidth]{DisorderAverages_Histogram_CloseUp}
\includegraphics[width=0.4\linewidth]{DisorderAverages_Histogram_Frequencies_vs_Energies}
\caption{\label{DisorderAverages_Histogram}Quantum corrections are crucial for obtaining the experimentally observed Gaussian energy distribution function for monolayer graphene on SiO$_2$ at charge neutrality. Left: Density histograms for $n_{\mathrm{TF}}$ and $n_3$, with and without exchange (X), for 500 disorder realizations with ${\mw{n_j}=0}$. We choose the bin widths for the densities such that 100 counts per bin are obtained on average. In contrast to the TF approximation, the quantum-corrected density distribution $n_3$ (including exchange) is captured reasonably well by a Lorentzian and even better by a Gaussian fit. Right: Translating local densities into energies ${E_D=\mathrm{sgn}\big(n(\vec r)\big)\,\hbar v_{\mathrm{F}}\,\sqrt{\pi|n(\vec r)|}}$, with the signum function $\mathrm{sgn}(\,)$, we find that the local quantum-corrected energies $E_D$ follow a Gaussian distribution that resembles the Gaussian fit to the experimental data `exp.' (extracted from Ref.~[\onlinecite{xue_scanning_2011}]) remarkably well. This is in stark contrast to the results of the TF approximation. We convert counts/bin into frequencies by renormalizing the histograms with the maxima of the Gaussian fits `G' and `G(TF)', respectively.}
\end{center}
\end{figure}

\end{widetext}

\section{\label{Bilayersetup}Double layer setup}

The treatment of monolayers in the previous section forms our basis for the description of more complicated heterostructures. Figure~\ref{DoubleLayerSchematics} illustrates a two-layer system, where both layers L1 and L2 are sandwiched between h-BN and subjected to a charged impurity layer from the SiO$_2$ substrate. We expose L2 to the same charge disorder that affects L1, though at a larger separation, but refrain from adding disorder on L2 in order to avoid confusing inter-layer correlation effects with effects from independent disorder on L2. For the same reason we model the strain of L1 and the charge disorder by the identical type of disorder, albeit in different realizations. The layer separation of $5\,$nm suffices to justify a merely classical electrostatic interaction between L1 and L2, i.e., inter-layer tunneling of charge carriers can be neglected --- in contrast to intra-layer tunneling through the disorder potential landscape. The latter is missed by densities in TF approximation but captured (in part) via the higher-order Suzuki-Trotter factorizations.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.9\linewidth]{DoubleLayerSchematics}
\caption{\label{DoubleLayerSchematics}Schematics of the double-layer setup investigated here. Two finite-size monolayers are embedded in h-BN at a layer separation of $5\,$nm, while the SiO$_2$ substrate provides the charged impurities at an effective distance\cite{Tan2007,Chen2008b,Samaddar2016} of $1\,$nm from the first layer L1. For the purpose of revealing the impact of strain and charge disorder on the electron-hole-puddle correlations between L1 and L2, we model L2 as an unstrained clean layer. We use the same setup for analyzing both monolayer and bilayer graphene heterostructures.}
\end{center}
\end{figure}

In what follows we address the puddle correlations between L1 and L2 as a function of the ratio $R$ between the disorder strength of the strain and that of the charged impurities. Further details are provided in Appendix~\ref{AppendixDisorder}. The sum of $V_\mathrm{C}$ and $V_\mathrm{S}$ results in electron-hole puddles within the first layer L1, whose electrostatic potential adds to the external potential $V_\mathrm{C}$ for the charge carriers in the second layer\footnote{Strictly speaking, the electrostatic effects of L1 and L2 on each other have to be treated self-consistently, but the layer separation of 5nm ensures that the backreaction of L2 on L1 is of the order of $10^{-6}$ of the external potential of L1 and therefore negligible at the level of precision we are aiming for in this work.}. We expect maximal puddle correlation if ${R=0}$, that is, when no strain can obscure the then dominating effect of the charged impurities on both layers: Owing to the intra-layer Coulomb interaction, the puddles in L1 are much reduced in weight compared with the case of noninteracting carriers. That is, the tendency of a puddle in L1 to electrostatically induce a puddle of opposite charge in L2 is overcompensated by the charge disorder, which exhibits the tendency to induce a puddle of the same charge\footnote{In the (unphysical) limit of vanishing layer separation these tendencies become certainties, but a finite layer separation in concert with the disordered potential landscape allows for local variations of correlation, thereby affecting global correlation measures quantitatively.}. Following the same line of reasoning, we expect maximal anticorrelation if ${R=\infty}$, with the transition from correlation to anticorrelation occuring at some value ${R>1}$. In the following section we substantiate these claims with quantitative predictions for graphene and bilayer graphene.

\section{\label{Results}Density correlations in double-layers of mono- and bilayer graphene}

We quantify the inter-layer correlations of electron-hole-puddles of the double-layer system described in Sec.~\ref{Bilayersetup} by solving Eqs.~(\ref{nasdef}) and (\ref{Vdef}) self-consistently\footnote{During the self-consistent loop the chemical potential $\mu$ is adjusted such that the mean carrier densities ${\mw{n_j}}$ are kept at a fixed value (zero, unless stated otherwise).} and by comparing the converged carrier densities ${n_1=n(\mathrm{L1})}$ and ${n_2=n(\mathrm{L2})}$ of layers L1 and L2 locally. To that end we calculate the two correlation measures $\xi[n_1,n_2]$ and $\xi_c[n_1,n_2]$, which yield a value of one for perfectly correlated electron-hole puddles (i.e., if the density distribution $n_2$ is proportional to $n_1$ and their values have the same sign at each position $\vec r$), minus one for perfect anticorrelation (i.e., if $n_2$ is proportional to $-n_1$), and are designed for tracking the transition between these two extremes; see Appendix~\ref{AppendixCorrelation} for details.

Figure~\ref{ST3+X_DoubleLayer} depicts potentials and densities for graphene, viz.~linear dispersion, calculated for mean carrier densities ${\mw{n_j}=0}$ and equal strengths of strain and charge disorder (ratio ${R=1}$). Due to the screening effects of the Coulomb interaction within L1, the effective potential for L1 exhibits less variability then the total external potential, i.e., the sum of the disorder potentials $V_{\mathrm{C}}$ and $V_{\mathrm{S}}$. The quasiparticle density $n_1$ of L1 , which can be viewed as resulting from this effective potential, induces an external electrostatic potential for the carrier density $n_2$ of L2. However, for ${R=1}$ the charge disorder potential dominates the total external potential for $n_2$ with magnitudes by a factor of more than 50 larger than those of the electrostatic potential caused by $n_1$. For the setting that leads to Fig.~\ref{ST3+X_DoubleLayer}, the magnitudes of the total external potential for L2 are smaller than those of L1 by a factor of 3--5.
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.85\linewidth]{ST3+X_DoubleLayer}
\caption{\label{ST3+X_DoubleLayer}Visibly correlated carrier densities $n_3$ (including exchange) for layers L1 and L2 (bottom row) from impurity and strain potentials of equal strength ${R=1}$ (top row): This figure illustrates the intermediate potentials relevant in the work flow towards carrier distributions for a double-monolayer graphene system as depicted in Fig.~\ref{DoubleLayerSchematics}. The color codes are in units of $u$ for energy and $l^{-2}$ for density, respectively, and apply to both graphics in each row. Top row: Charged impurity potential $V_{\mathrm{C}}$ (left) and strain potential $V_{\mathrm{S}}$ (right) for L1 --- Second row: Total external potential $V_{\mathrm{ext}}$ (left) and converged effective potential $V$ (right) for L1 ---  Third row:  Charged impurity potential $V_{\mathrm{C}}$ for L2 (left) and electrostatically induced potential $V_{\mathrm{ind}}$ on L2 from the charge distribution $n_1$ of L1 (right, scaled up by a factor of 50) --- Fourth row: Total external potential $V_{\mathrm{ext}}$ (left) and converged effective potential $V$ (right) for L2 --- Bottom row: Converged carrier densities $n_1$ (left, scaled by a factor of 1000) and $n_2$ (right, scaled by a factor of 5000).}
\end{center}
\end{figure}
\noindent As a result, the density fluctuations of $n_2$ are diminished compared with those of $n_1$ by a factor of 5--10. It is therefore well justified to refrain from a self-consistent treatment of the electrostatically induced potentials $V_{\mathrm{ind}}$ of both layers, and to consider only $V_{\mathrm{ind}}$ stemming from $n_1$. As is evident from the bottom row of Fig.~\ref{ST3+X_DoubleLayer}, the spatial distributions of $n_1$ and $n_2$ are correlated rather than anticorrelated for ${R=1}$.

Repeating the calculation which yields the results illustrated in Fig.~\ref{ST3+X_DoubleLayer} for different values of $R$, we find the critical value $R_0$ at the transition from inter-layer correlation to anticorrelation. The correlation measure $\xi$ signifies inter-layer puddle-correlation (anticorrelation) by taking on positive (negative) values. The corresponding diagram in Fig.~\ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_ErrorBars}, calculated with charged impurities density $\bar{n}_{\mathrm{imp}}=10^{12}\,\mathrm{cm}^{-2}$, exhibits ${R_0\approx 10\pm2}$ when $n_3$ is used, and a similar value in the case of $n_{\mathrm{TF}}$. Here, we take into account exchange effects and report a rough error estimate simply based on our numerical findings from five disorder realizations\footnote{At this point we do not calculate disorder averages, but rather showcase the predictions of a few disorder realizations (also bearing in mind the limited number of samples usually investigated in actual experiments).}. Evidently, the disparity in electron-hole puddle landscapes between $n_3$ and $n_{\mathrm{TF}}$ as seen in Fig.~\ref{TFvsST3+X_monolayer} does not translate into an appreciable difference between $\xi^{(3)+\mathrm{X}}$ (for $n_3$) and $\xi^{\mathrm{TF+X}}$ (for $n_{\mathrm{TF}}$).

Although it is not surprising per se that integrated quantities like $\xi$ are less sensitive to local differences between $n_3$ and $n_{\mathrm{TF}}$, it cannot be assumed a priori. Our quantitative analysis shows that a quasi-classical approach to inter-layer correlations in monolayer-graphene heterostructures is justified in case of rather large impurity densities like of $10^{12}\,\mathrm{cm}^{-2}$ used for Fig.~\ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_ErrorBars}. Commonly, however, the TF approximation is less reliable for smaller particle numbers, and quantum corrections can be expected to play a more dominant role as the carrier density is reduced. The disorder potential for smaller $\bar{n}_{\mathrm{imp}}$ is less pronounced and gives rise to puddles that exhibit smaller carrier densities on average. Indeed, with ${\bar{n}_{\mathrm{imp}}=10^{11}/\mathrm{cm}^2}$ employed for Fig.~\ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_LowNimp2_ErrorBars}, the quantum-corrected $\xi^{(3)+\mathrm{X}}$ can be clearly distinguished from the quasi-classical $\xi^{\mathrm{TF+X}}$. Our data shown in Fig.~\ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_LowNimp2_ErrorBars} point to the critical value ${R_0^{\mathrm{TF+X}}\approx 15\pm 2}$ with an error estimate similar to that in Fig.~\ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_ErrorBars}. An increased $R_0$ at lower $\bar{n}_{\mathrm{imp}}$ can be understood from the following simplified picture. With typical values $\left.v_\mathrm{C}\right|_d$ and $v_\mathrm{S}$ of the charge and strain potentials for L1 and neglecting interactions, we have typical values ${n_{\mathrm{TF}}(\mathrm{L1})\propto\big(\mu-v_{\mathrm{ext}}(\mathrm{L1})\big)^2\approx(\left.v_\mathrm{C}\right|_d+v_\mathrm{S})^2}$ for ${\mu\approx0}$, i.e., approximately ${n_{\mathrm{TF}}(\mathrm{L1})\propto v_\mathrm{S}^2}$ at ${R\approx10}$. Then, the typical TF density in L2 is determined by ${v_{\mathrm{ext}}(\mathrm{L2})=\left.v_\mathrm{C}\right|_{\tilde d}+v_\mathrm{ind}}$, where the typical values ${\left.v_\mathrm{C}\right|_{\tilde d}\approx\frac12\left.v_\mathrm{C}\right|_d}$ at ${\tilde{d}-d\approx5}\,$nm and ${v_\mathrm{ind}\propto v_\mathrm{S}^2}$ are roughly equal if ${\bar{n}_{\mathrm{imp}}=10^{12}/\mathrm{cm}^2}$ and ${R=10}$. For a scaled ${\bar{n}_{\mathrm{imp}}=\lambda\times10^{12}/\mathrm{cm}^2}$ and the same ${R=10}$, $v_{\mathrm{ext}}(\mathrm{L1})$ scales with $\lambda$ as well, but the typical values of ${n_{\mathrm{TF}}(\mathrm{L1})}$ then scale like $\lambda^2\,v_\mathrm{S}^2$, such that ${v_{\mathrm{ext}}(\mathrm{L2})=\lambda\left.v_\mathrm{C}\right|_{\tilde d}+\lambda^2\,v_\mathrm{ind}}$. For the case of ${\lambda=1/10}$, as represented by Fig.~\ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_LowNimp2_ErrorBars}, the strength of $v_\mathrm{ind}$, relative to $\left.v_\mathrm{C}\right|_{\tilde d}$, is diminished by a factor of ten. As strain feeds into $v_\mathrm{ind}$, not into $\left.v_\mathrm{C}\right|_{\tilde d}$, (relatively) more strain is required to counteract the effect of $\left.v_\mathrm{C}\right|_{\tilde d}$, implying a larger critical $R_0$ at lower impurity densities. Disregarding the uncertainties for the critical $R$, one could estimate ${R_0^{(3)+\mathrm{X}}\approx 30}$ for the quantum-corrected correlations. However, $\xi^{(3)+\mathrm{X}}$ rather exhibits an extended cross-over regime, where the magnitude of strain can be varied substantially with little effect on the average correlation $\xi^{(3)+\mathrm{X}}$. Our main result is thus that corrections to the TF approximation can become important even for integrated or averaged observables of 2D materials --- given the proper conditions, for example, small carrier densities.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.9\linewidth]{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_ErrorBars}
\caption{\label{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_ErrorBars} Both the TF- and quantum-corrected approach predict a transition from correlated to anticorrelated electron-hole-puddles at ${R=R_0^{(3)+\mathrm{X}}\approx R_0^{\mathrm{TF+X}}\approx 10\pm2}$, where $R$ is the ratio of the strength of strain to charged impurity disorder. The correlation measure $\xi$ is plotted as a function of $R$ for five disorder realizations of double monolayer graphene, extracted from both TF- and quantum-corrected carrier densities with exchange incorporated. The average density of charged impurities is ${\bar{n}_{\mathrm{imp}}=10^{12}/\mathrm{cm}^2}$. The error bars demarcate minimal and maximal values of $\xi$ found within the set of disorder realizations; the shaded areas and the dotted line segments guide the eye; the horizontal line at ${\xi=0}$ separates the correlation phase (above) from the phase of anticorrelation (below). We use a logarithmic horizontal axis in order to more clearly showcase the transition point and the contrast of $\xi$ over a large range of $R$.}
\end{center}
\end{figure}

\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.9\linewidth]{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_LowNimp2_ErrorBars}
\caption{\label{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_LowNimp2_ErrorBars} At lower values of ${\bar{n}}_{\mathrm{imp}}$, quantum corrections lead to stark differences in the crossover behavior. Here we repeat the plot in Fig.~\ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_ErrorBars}, but for ${\bar{n}_{\mathrm{imp}}=10^{11}/\mathrm{cm}^2}$, and an increased sheet size to ensure enough disorder statistics (640 impurities on $(800\,\mathrm{nm})^2$). The extended crossover of $\xi^{(3)+\mathrm{X}}$ from correlation to anticorrelation contrasts with the pattern observed in Fig.~\ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_ErrorBars} and heralds the emergence of quantum--effects (i.e.,  major deviations from $\xi^{\mathrm{TF+X}}$) with decreasing particle numbers.}
\end{center}
\end{figure}

We briefly turn to the impact of exchange effects in inter-layer correlations. Since quantum corrections modify the TF carrier distributions of the individual layers more profoundly than exchange does, cf.~Fig.~\ref{TFvsST3+X_monolayer}, we expect that exchange plays a minor role in determining carrier correlations of double-layer systems. This is confirmed with Table~\ref{Monolayer_graphene_Strain_vs_Charge_R_diagram_withInset_TF_xi_Xnegligible} in Appendix~\ref{AppendixExchange}, where $\xi$ is determined in TF approximation for a single disorder realization, with and without exchange. In view of the magnitude of uncertainties displayed in Figs.~\ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_ErrorBars} and \ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_LowNimp2_ErrorBars} we find negligible quantitative differences when omitting or including exchange.

\begin{figure}[htb!]
\begin{center}
\includegraphics[width=0.9\linewidth]{Bilayer_graphene_Strain_vs_Charge_MinMax_TF_vs_ST3_ErrorBars}
\caption{\label{Bilayer_graphene_Strain_vs_Charge_MinMax_TF_vs_ST3_ErrorBars} Like Fig.~\ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_ErrorBars}, but for double bilayer graphene. Using $n_3$, we extract ${R_0^{(3)}\approx 10\pm3}$ for the transition from correlated to anticorrelated carrier densities, but a somewhat smaller value ${R_0^{\mathrm{TF}}\approx 8\pm1}$ in case of the TF approximation. Similar to our observations for double-monolayer graphene, we find that quantum corrections have little effect in the case of double-bilayer graphene for $\bar{n}_{\mathrm{imp}}=10^{12}\,\mathrm{cm}^{-2}$
}
\end{center}
\end{figure}

\begin{figure}[htb!]
\begin{center}
\includegraphics[width=0.9\linewidth]{xi_n2_R_Bilayer_graphene_TF_200x200_Disorder02}
\caption{\label{xi_n2_R_Bilayer_graphene_TF_200x200_Disorder02} Determination of the critical ratio ${R_0^{\mathrm{TF}}\approx 8}$ for double bilayer graphene (found in Fig.~\ref{Bilayer_graphene_Strain_vs_Charge_MinMax_TF_vs_ST3_ErrorBars}) as a function of $\mw{n_2}$, viz.~$\mu(\mathrm{L2})$. The findings for $\xi^{\mathrm{TF}}$ in Fig.~\ref{Bilayer_graphene_Strain_vs_Charge_MinMax_TF_vs_ST3_ErrorBars} correspond to ${\mw{n_2}=0}$ here. All data points are obtained for a single fixed disorder realization
}
\end{center}
\end{figure}

Our results on double bilayer graphene in Fig.~\ref{Bilayer_graphene_Strain_vs_Charge_MinMax_TF_vs_ST3_ErrorBars} show a slight quantitative difference between the transitions of $\xi^{(3)}$ and $\xi^{\mathrm{TF}}$ for the here employed density of charged impurities ${\bar{n}_{\mathrm{imp}}=10^{12}\,\mathrm{cm}^{-2}}$. However, the overall trend for both approximations is still very similar and the spread of $\xi$, originating in multiple disorder realizations like in Figs.~\ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_ErrorBars} and \ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_LowNimp2_ErrorBars}, still implies ${R_0^{(3)}\approx R_0^{\mathrm{TF}}}$ within the depicted uncertainties. Extrapolating from the double-monolayer case the double-bilayer graphene can be expected to exhibit more pronounced differences between $\xi^{(3)}$ and $\xi^{\mathrm{TF}}$ for smaller density of charged impurities.

Figure~\ref{xi_n2_R_Bilayer_graphene_TF_200x200_Disorder02} provides another angle on the transition from correlation to anticorrelation for double bilayer graphene: For fixed ${\mw{n_1}=0}$ we chart $\xi^{\mathrm{TF}}$ as a function of $\mw{n_2}$ and find a qualitative change in the shape of the curves $\xi(\mw{n_2})$ as $R$ transits across ${R_0^{\mathrm{TF}}\approx 8}$, consistent with our findings in Fig.~\ref{Bilayer_graphene_Strain_vs_Charge_MinMax_TF_vs_ST3_ErrorBars}.

\FloatBarrier

\section{Conclusions and perspectives}

The main result of this work is thus that both exchange and quantum corrections beyond the Thomas--Fermi approximation are essential for quantitatively reliable density distributions in the context of 2D materials. While exchange plays a minor role for some integrated quantities like inter-layer correlations of double-layer systems, quantum corrections become important for smaller average carrier-densities, viz. cleaner systems. In other words, systems like graphene on boron nitride samples with ${\bar{n}_{\mathrm{imp}} \approx 10^{10} \mathrm{cm}^{-2} }$ call for a more sophisticated description than the Thomas--Fermi approximation can provide. We showed that the regime of high carrier densities may still be described quasiclassically in the computationally highly efficient Thomas--Fermi approximation. For specific systems, however, quantum corrections have to be taken into account at the quantitative level. In this work we provided one such example when analyzing the crossover from correlation to anticorrelation in bilayer heterostructures of 2D materials. We focussed on describing double-layers of mono- and bilayer graphene via recently developed techniques of calculating quantum corrections, but our approach can be easily adapted for other heterostructures.

Turning to the issue of Coulomb drag~\cite{gorbachev_strong_2012}, our analysis reveals that the strain potential has to be at least ten times stronger than the impurity potential in order for anti-correlation between the layers to occur in both monolayer and bilayer graphene. Any strain potential weaker than this will result in correlation. The DPFT framework presented here thus provides the following path towards an experimental verification of the nature of correlations in graphene heterostructures, a non-trivial task since it is impossible to directly measure the local densities of encapsulated 2D layers via local scanning tunneling microscopy. The strength of the strain potential~\cite{gibertini_electron-hole_2012} may first be extracted from local height fluctuation data obtained via atomic force microscopy~\cite{dean_boron_2010}, followed by the impurity strength from conductivity measurements~\cite{adam_self-consistent_2007}. The resulting ratio $R$ may then be compared against Figs.~\ref{Monolayer_graphene_Strain_vs_Charge_MinMax_TF+X_vs_ST3+X_ErrorBars} to \ref{Bilayer_graphene_Strain_vs_Charge_MinMax_TF_vs_ST3_ErrorBars} to determine the nature of correlation between the layers. This would be especially useful for clarifying the source of the unexpected sign changes~\cite{li_negative_2016,lee_giant_2016,simonet_anomalous_2017} in Coulomb drag experiments that have been explained by earlier works simply asserting correlation~\cite{song_energy-driven_2012,schutt_coulomb_2013} or anti-correlation~\cite{gorbachev_strong_2012,ho_theory_2018} without proof. In the case of bilayer graphene, Ref.~[\onlinecite{zarenia_multiband_2018}] theoretically demonstrated that a multiband mechanism based on the thermal smearing of Fermi energy in both layers explains the unexpected sign changes at high densities where puddles may be ignored. A complete analysis at the double neutrality where both puddles and multiband effects are important is however still lacking. In future, this may be addressed by generalizing the DPFT framework presented here to incorporate multiband effects~\cite{zarenia_multiband_2018} in the determination of puddle-induced quantities such as the interlayer correlation and density fluctuations.

Finally, the findings of this work are useful in the study of van der Waals heterostructures~\cite{geim_van_2013} consisting of a pair of two-dimensional electronic layers separated by thin dielectric spacers.
These structures serve as ideal platforms for the study of a range of interesting physical effects such as exciton condensation~\cite{li_excitonic_2017,liu_quantum_2017,burg_strongly_2018,lopez_rios_evidence_2018} and strong light-matter interaction~\cite{britnell_strong_2013,yu_highly_2013,withers_light-emitting_2015} and are thus an area of intense research activity.
It is likely that the inter-layer puddle correlations calculated within this work plays a role in the above physics.

\acknowledgments
We acknowledge the support of the National Research Foundation of Singapore under its Fellowship program (NRF-NRFF2012-01) and the National University of Singapore Young Investigator Award (R-607-000-094-133). We also acknowledge the use of the dedicated computational facilities at the Centre for Advanced 2D Materials.


\section{Introduction}

\subsection{Spin synchronization model}
We consider the problem of reconstructing $n$ independent, uniform $(\pm 1)$-valued spins $X_1,\dots,X_n \stackrel{i.i.d.}{\sim} \Rad(1/2)$ living on the vertices of an $n$-vertex graph $G$, by observing their interactions $Y_{e}$ on the edges of the graph. For each $e = (i,j) \in E(G)$, the interaction $Y_e = Y_{ij}$ depends only on $X_i \cdot X_j$. Formally, we can factor the joint distribution as $$P_{X,Y}(x,y) = \left(\prod_{(i,j) \in E(G)} P_{Y_{ij}|X_i \cdot X_j}(y_{ij}|x_i \cdot x_j)\right) \cdot P_{X}(x).$$ In other words, each $Y_{ij}$ is the output of $X_i \cdot X_j$ through a channel $Q_{ij} := P_{Y_{ij} | X_i \cdot X_j}$.

\begin{definition}
We refer to $(X,Y)$ as a spin synchronization model on a graph $G$ with edge channels $Q$.
\end{definition}

This spin synchronization model has previously appeared in \cite{our_grid}, for example. Note that in \cite{our_grid} the alphabets were not restricted to be binary, but were arbitrary groups instead.

Given a spin synchronization problem $(X,Y)$, the goal is to reconstruct the spin $X_u$ for some $u \in V(G)$, given $Y$ and $X_W$ where $W \subset V(G)$ (due to symmetry, we may as well freeze at least one reference vertex). It is thus natural to look at some measure of dependency between to generic vertices in the graph, and a natural quantity of interest is the mutual information $I(X_u; X_W, Y)$.

\paragraph*{Connection to other reconstruction models}
Depending on the choices of the graph and the edge interaction channels, the spin synchronization model captures models such as (a) broadcasting on trees \cite{evans}, (b) censored block models \cite{heimlicherlelargemassoulie,abbs}, (c) synchronization on grids \cite{our_grid}, and (d) spiked Wigner models \cite{yash_sbm}.
Thus, several information-theoretic reconstruction thresholds for these models can be obtained by studying the information-theoretic thresholds for the spin synchronization model. We refer the reader to \cite{abbe2018information,yury_percolation} for details on these connections.

\paragraph*{Connection to Ising Model} The spin synchronization model is also related to the Ising model in statistical physics; conditioned on the edge observations $Y$, the posterior distribution of the vertex spins $X$ is given by an Ising model. However, we will interested here in the average-case behavior over the edge observations in the model, so naively applying correlation decay results on Ising models (e.g., Dobrushin conditions \cite{dobrushin1968problem}) would yield weaker bounds than those we will obtain.

\subsection{Main result}

The main contribution of this paper is Theorem \ref{thm:mainseriesparallelbound}, a new impossibility result for reconstruction in the spin synchronization model.\footnote{Part of the results in this note appeared in Enric Boix's undergraduate senior thesis \cite{enric2018thesis}.} This result is incomparable to the previously-known bounds derived in \cite{abbe2018information} and \cite{yury_percolation}. We discuss in detail the relationship between the new bound and the previous bounds, and in Section \ref{sec:generalizedbot} we give examples of spin synchronization problems for which the new bound is tighter. In Section \ref{sec:botapplication}, we also use Theorem \ref{thm:mainseriesparallelbound} to rederive the impossibility result of \cite{evans} for broadcasting on trees.

For our main result, we work with the Chi-squared mutual information, $I_2$, instead of the classical KL mutual information, $I_{KL}$. See Appendix~\ref{app:chi2mutualinformation} for definitions. As we will see, the Chi-squared mutual information turns out to be a more natural choice for the proof of the theorem. In any case, for our purposes the KL and Chi-squared mutual informations are equal up to a factor of 2, because of the following well-known inequalities $\frac{1}{2}I_2(A; U) \leq I_{\mathrm{KL}}(U; A) \leq I_2(U; A)$ when $U \sim \Rad(1/2)$, reproved in Appendix \ref{app:chi2mutualinformation}.

We are ready to state our main result.

\begin{theorem}\label{thm:mainseriesparallelbound} Let $(X,Y)$ be a spin synchronization model on a series-parallel graph $G$ with terminals $u$ and $v$, and arbitrary edge channels $Q$. Then $$I_2(X_u; X_v \mid Y_{E(G)}) \leq \sum_{P \in \cP_G(u,v)} I_2(X_u; X_v \mid Y_{E(P)}).$$ Here $\cP_G(u,v)$ is the set of paths (self-avoiding walks) from $u$ to $v$ in $G$.
\end{theorem}

\subsection{Comparison to previous results}
\paragraph{Comparison to \cite{abbe2018information}} The recent work \cite{abbe2018information} derives a subadditivity bound similar to Theorem \ref{thm:mainseriesparallelbound}. While the bound of \cite{abbe2018information} applies to spin synchronization models on general graphs, it requires the edge channels of the models to have certain symmetry. In contrast, in Theorem \ref{thm:mainseriesparallelbound}, we handle general channels, while giving up generality on the base graph. For ease of comparison, we state the result of \cite{abbe2018information} below.

\begin{theorem}[ \cite{abbe2018information}]\label{thm:mainsymmetricbound}
Let $(X,Y)$ be a spin synchronization model on a graph $G$, such that the edge channels $Q$ are symmetric: for each $e \in E(G)$, there is a measurable transformation $T_e = T_e^{-1}$ such that $Q_e(\cdot | +1) = T_e \sharp Q_e(\cdot | -1)$. ($T_e\sharp$ denotes the push-forward operation.)

Then for all $u \in V(G), W \subset V(G)$, $$I_2(X_u; X_W\mid Y) \leq \mathrm{conn}_{G,\gamma}(u, W),$$ where $$\gamma((i,j)) = I_2(X_i; X_j\mid Y_{ij})$$ for all $(i,j) \in E(G)$. Here, $\mathrm{conn}_{G,\gamma}(u, W)$ denotes the probability that $u$ is in the same open component as some $w \in W$ in a bond percolation on $G$ where each edge $e$ is independently open with probability $\gamma(e)$.
\end{theorem}

This gives in turn a subadditivity property for the Chi-squared mutual information analogous to Theorem \ref{thm:mainseriesparallelbound}:

\begin{corollary}[from \cite{abbe2018information}]\label{cor:subadditivitybound}
Let $(X,Y)$ be as in Theorem \ref{thm:mainsymmetricbound}. Then for all $u \in V(G), W \subset V(G)$, $I_2$ is subadditive over paths:
\begin{align*}I_2(X_u; X_W\mid Y_{E(G)}) &\leq \sum_{P \in \cP_G(u,v)} I_2(X_u; X_v\mid Y_{E(P)}),\end{align*}
where $\cP_G(u,W)$ is the set of paths (i.e., self-avoiding walks) from $u$ to $W$ in $G$.
\end{corollary}

Notice that when the edge channels are symmetric, the subadditivity bound of Theorem \ref{thm:mainseriesparallelbound} may be weaker than the percolation bound of Theorem \ref{thm:mainsymmetricbound}. However, Theorem \ref{thm:mainseriesparallelbound} also applies to spin synchronization instances with asymmetric channels, while Theorem \ref{thm:mainseriesparallelbound} does not.

For completeness, we include a shortened proof of Theorem \ref{thm:mainsymmetricbound} and Corollary \ref{cor:subadditivitybound} in Appendix~\ref{sec:thmmainsymmetricboundproof}.

\paragraph{Comparison to \cite{yury_percolation}} An  information-percolation bound  is also obtained in \cite{yury_percolation} using an argument based on Strong Data-Processing Inequalities (SDPI). The bound of \cite{yury_percolation} applies to general graphs, to general channels, and also to a model with more general vertex labels. In order to compare the bound of \cite{yury_percolation} with Theorem \ref{thm:mainseriesparallelbound}, we state the result of \cite{yury_percolation} for binary labels:

\begin{theorem}
[Theorem 2 in \cite{yury_percolation} for spins]
\label{thm:yurytheorem}
Let $(X,Y)$ be a spin synchronization model on a graph $G$. Then, for all $u \in V(G), W \subset V(G)$, $$I_{KL}(X_u; X_W | Y) \leq \mathrm{conn}_{G,\gamma}(u,W),$$ where $$\gamma((i,j)) = \eta(Q_{ij})$$ is the SDPI constant for the edge channel $Q_{ij}$. Again, $\mathrm{conn}_{G,\gamma}(u,W)$ denotes the probability that $u$ is in the same open component as some $w \in W$ in a bond percolation on $G$ where each edge $e$ is independently open with probability $\gamma(e)$.
\end{theorem}
When the channels $Q$ are symmetric in the sense of Theorem \ref{thm:mainsymmetricbound}, the SDPI constant $\eta(Q_{ij})$ equals the Chi-squared mutual information $I_2(X_i; X_j | Y_{ij})$, and therefore Theorem \ref{thm:yurytheorem} almost generalizes Theorem \ref{thm:mainsymmetricbound}, besides for a minor multiplicative factor of 2 resulting from the relationship between $I_{KL}$ and $I_2$ (see Proposition \ref{prop:evanschi2vsklonradhalf}).

However, when the channels $Q$ are asymmetric, $\eta(Q_{ij})$ can be significantly larger than $I_2(X_i; X_j | Y_{ij})$. As a result, our subadditivity bound in Theorem \ref{thm:mainseriesparallelbound} can be used to show that the bound of Theorem \ref{thm:yurytheorem} is loose. We demonstrate this explicitly in Section \ref{sec:generalizedbot} below, by analyzing a spin synchronization problem on a tied-tree model.

\section{Synchronization on trees: a generalization of  broadcasting on trees}

In this section, we first show how the classical broadcasting on tree (BOT) problem of \cite{evans} can be expressed as a synchronization problem on a series-parallel graph consisting of a tied-tree with symmetric channels (Section \ref{sec:botapplication}). Applying our Theorem \ref{thm:mainseriesparallelbound}, we then obtain an alternative proof of the impossibility result of \cite{evans} for BOT reconstruction.

We then generalize this BOT model  (Section \ref{sec:generalizedbot}), using more general channels in the  synchronization on tree (SOT) model. For this generalized BOT, we compare the information bound given by our Theorem \ref{thm:mainseriesparallelbound} based on the Chi-squared subadditivity with the information bound given by Theorem \ref{thm:yurytheorem} from \cite{yury_percolation} based on the SDPI constant.

\subsection{Broadcasting on trees}\label{bot}
In the BOT model of \cite{evans}, a random variable is broadcast from the root down the edges of a tree, with each edge flipping the variable with some probability. The goal is to reconstruct the root variable by observing of the variables of some of its descendants. Formally, each vertex $v \in V(T)$ of a tree $T$ has a binary hidden label $\sigma_v$. The hidden labels are assigned by letting the root $\rho$ have spin $\sigma_{\rho} \sim \Rad(1/2)$, and by defining edge labels $\{\eta_e\} \stackrel{i.i.d.}{\sim} \Rad(\varepsilon)$, and letting $$\sigma_{v} = \sigma_{\rho} \prod_{e} \eta_e,$$ where the product is over the edges in the path from $\rho$ to $v$.

In \cite{evans}, the following information subadditivity inequality is proved, for any finite set of vertices $W \subset V(T)$: \begin{equation}\label{eq:evanssubadditivity}I_{KL}(\sigma_{\rho}; \sigma_W) \leq \sum_{w \in W} I_{KL}(\sigma_{\rho}; \sigma_w).\end{equation} This suffices to show that if $T$ is infinite and has percolation threshold $p_c(T)$, then, when $(1-2\varepsilon)^2 < p_c(T)$, it is impossible to reconstruct $\sigma_{\rho}$ from observations of the leaf variables at infinite depth.

\subsection{Reduction from SOT to BOT}\label{sec:botapplication}

The BOT model with noise parameter $\e$ and tree $T$ that generates the random variables $(\sigma, \eta)$ is equivalent to the synchronization on tree (SOT) model with noise parameter $\e$ and tree $T$ that generates the random variables $(X,Y)$, where the vertex labels $\{X_v\}$ are i.i.d.\ $\Rad(1/2)$, and the edge channels $\{Q_{ij}\}$ are binary symmetric channels (i.e., $Y_{ij} = X_i \cdot X_j \cdot Z_{ij}$, where the $Z_{ij}$ are i.i.d $\Rad(\varepsilon)$). By equivalent, we mean the following.
\begin{proposition}\label{prop:botspinreconstruction}
In the reconstruction setting just described, for all $W \subset V(T)$, \begin{equation*}I(\sigma_{\rho}; \sigma_W) = I(X_{\rho}; X_W, Y_{E(T)}),\end{equation*} where $I$ is any $f$-mutual information (such as $I_{KL}$ or $I_2$).
\end{proposition}
\begin{proof}
First, define $\sigma'_{\rho} = X_{\rho}$ and $\sigma'_w = X_w \cdot \prod_e Y_e$ for all $w \in W$, where the product is over the edges on the path from $\rho$ to $w$. By data-processing and because $\sigma'_{W \cup \{\rho\}} \stackrel{d}{=} \sigma_{W \cup \{\rho\}}$, we have $I(X_{\rho}; X_W, Y) \geq I(\sigma'_{\rho}; \sigma'_W) = I(\sigma_{\rho}; \sigma_W)$.

Conversely, let $X'_{\rho} = \sigma_{\rho}$, $Y'_e \stackrel{i.i.d}{\sim} \Rad(1/2)$, and $X'_{w} = \sigma_{w} \cdot \prod_{e} Y'_e$. By data-processing and equality of distributions, $I(\sigma_{\rho}; \sigma_{W}) = I(\sigma_{\rho}; \sigma_{W}, Y') \geq I(X'_{\rho}; X'_{W}, Y') = I(X_{\rho}, X_W, Y)$.
\end{proof}

Since we wish to apply Theorem \ref{thm:mainseriesparallelbound} in order to bound $I(X_{\rho}; X_W, Y_{E(T)})$, and since we expressed our theorem in terms of  series-parallel graphs with two terminals, we first note that can consider an equivalent varient of SOT using a series-parallel graph by tying the leaves to a terminal with noiseless channels. Namely,  used $(X_{V(T)}, Y_{E(T)})$ to construct a spin synchronization model $(X_{V(G_W)}, Y_{E(G_W)})$ on a series-parallel graph $G_W$, where $W$ are the leaves of the tree $T$, and $G_W$ extends the tree by adding a new vertex $v$ adjacent to all $w \in W$ with noiseless edge channels between any $w$ and $v$. We call $G_W$ a ``tied-tree'' because the leaves $W$ of the tree $T$ are ``tied'' together in $G_W$ through their connections to vertex $v$.

It is easy to show via induction on the number of edges that tied-trees are indeed series-parallel graphs whose two terminals are the root of the original tree, and the new added vertex connecting the leaves. So the bound of Theorem \ref{thm:mainseriesparallelbound} applied to $G_W$ yields
\allowdisplaybreaks
\begin{align*}I_2(X_{\rho}; X_v, Y_{E(G_W)}) &= I_2(X_{\rho}; X_v | Y_{E(G_W)}) \tag{Prop. \ref{prop:botspinreconstruction}} \\ &\leq \sum_{P \in \cP_{G_W}(\rho,v)} I_2(X_{\rho}; X_v | Y_{E(P)}) \tag{Theorem \ref{thm:mainseriesparallelbound}} \\ &= \sum_{w \in W} \sum_{P \in \cP_{T}(\rho,w)} I_2(X_{\rho}; X_v | Y_{E(P) \cup (w,v)}) \\ &= \sum_{w \in W} \sum_{P \in \cP_{T}(\rho,w)} I_2(X_{\rho}; X_w | Y_{E(P) \cup (w,v)}) \tag{Since $Y_{(w,v)}$ is noiseless} \\ &= \sum_{w \in W} \sum_{P \in \cP_{T}(\rho,w)} I_2(X_{\rho}; X_w | Y_{E(P)}) \tag{Since $Y_{(w,v)} \indy X_{\rho}, X_w, Y_{E(P)}$} \\ &= \sum_{w \in W} \sum_{P \in \cP_{T}(\rho,w)} I_2(X_{\rho}; X_w | Y_{E(T)}) \tag{Since $Y_{E(T) \sm E(P)} \indy Y_{E(P)}, X_{\rho}, X_w$} \\ &= \sum_{w \in W} \sum_{P \in \cP_{T}(\rho,w)} I_2(X_{\rho}; X_w, Y_{E(T)}) \tag{Prop. \ref{prop:i2conditionequalsi2side}} \\ &= \sum_{w \in W} I_2(\sigma_{\rho}; \sigma_w) \tag{By Prop. \ref{prop:botspinreconstruction} and uniqueness of paths in tree}  \end{align*}

We conclude by noting that since $Y_{\{(w,v)\}_{w \in W}} \indy X_{\rho} | X_W, Y_{E(T)}$, $$I_2(X_{\rho}; X_v, Y_{E(G_W)}) = I_2(X_{\rho}; X_W, Y_{E(T)}, Y_{\{(w,v)\}_{w\in W}}) = I_2(X_{\rho}; X_W, Y_{E(T)}) = I_2(\sigma_{\rho}; \sigma_W).$$

Hence, we have re-derived the BOT subadditivity inequality \eqref{eq:evanssubadditivity} for the Chi-squared mutual information (and also for the KL mutual information with a factor of 2 by Proposition \ref{prop:evanschi2vsklonradhalf}). This in turn gives the impossibility result of \cite{evans}.

\subsection{Asymmetric SOT and Chi-squared v.s.\ SDPI}\label{sec:generalizedbot}
In Section \ref{sec:botapplication}, we formulated the spin broadcasting on trees problem of \cite{evans} as a spin synchronization problem on a tied-tree graph, with binary symmetric edge channels. If we relax the edge channels to allow for general (possibly asymmetric) channels, we thus obtain a generalization of the broadcasting on trees problem. Notice that the path subadditivity bound of Theorem \ref{thm:mainseriesparallelbound}  still applies to this generalized model, as the tied-tree is series-parallel.

In particular, suppose we are given a spin synchronization model $(X,Y)$ on a tied-tree $G$ with root $u$ and terminal $v$ similarly to previous section, but we now take the edge channels $Q$ on the underlying tree to be asymmetrical, i.e., $Q_{ij}(\cdot | +1) \sim \Ber(a/n)$ and $Q_{ij}(\cdot |-1) \sim \Ber(b/n)$, and the leaves of the tree are tied together at $v \in V(G)$ with noiseless edge channels. Suppose moreover that the tree is a regular tree with $d=n$ descendants at each generation and depth $t$.

For any edge $(i,j)$, we have $$I_2(X_i; X_j | Y_{ij}) = \frac{(a-b)^2}{2(a+b)n} + \frac{(a-b)^2}{n^2(1-(a+b)/(2n))} = \frac{(a-b)^2}{2(a+b)n} + o(1/n)$$ by direct calculation using Proposition \ref{prop:correqualschi2}. Hence, by Theorem \ref{thm:mainseriesparallelbound}, we have the bound \begin{align}I_2(X_u; X_v | Y_{E(G)}) &\leq \sum_{P \in \cP_{G}(u,v)} I_2(X_u; X_v | Y_{E(P)}) \tag*{} \\ &= \sum_{P \in \cP_{G}(u,v)} \left(\frac{(a-b)^2 + o(1)}{2(a+b)n}\right)^{|E(P)|}\label{eq:ourresontiedtree}\\
&=d^t \left(\frac{(a-b)^2}{2(a+b)n} + o(1/n)\right)^{t}\\
&=\left(\frac{(a-b)^2}{2(a+b)} + o_n(1)\right)^{t}\label{last1}
.\end{align}

On the other hand, as derived in \cite{polyanskiy2017strong}, the SDPI constant for the edge channels is $$\eta(Q_{ij}) = \frac{(\sqrt{a}-\sqrt{b})^2 }{n}  + o(1/n).$$ So letting $\gamma((i,j)) = \eta(Q_{i,j})$ for all $(i,j) \in E(G)$, the bound of \cite{yury_percolation} (Theorem \ref{thm:yurytheorem}) gives \begin{align}I_2(X_u; X_v | Y_{E(G)}) &\leq \mathrm{conn}_{G,\gamma}(u,v) \label{eq:yuryontiedtree} \end{align} And by inclusion-exclusion, \begin{align}\mathrm{conn}_{G,\gamma}(u,v) &\geq \sum_{P \in \cP_{G}(u,v)} \gamma^{|E(P)|} - \sum_{P_1\neq P_2 \in \cP_{G}(u,v)} \gamma^{|E(P_1) \cup E(P_2)|} \\ &= d^t \gamma^t - \frac{d^t \gamma^t}{2} \cdot \left(\sum_{k=1}^t\gamma^k (d-1) \cdot d^{k-1}\right) \\ &= d^t \gamma^t \left(1 - \frac{(d-1)(1-(\gamma d)^{t+1})}{2d(1-\gamma d)}\right) \\
&\geq \left((\sqrt{a}-\sqrt{b})^2 + o_n(1)\right)^{t} \left(1 - \frac{1}{2}\left((\sqrt{a}-\sqrt{b})^2 + o_n(1)\right)^{t+1}\right)
\end{align}
Note that for any $a,b>0$,
$$ (\sqrt{a} - \sqrt{b})^2 \ge \frac{(a-b)^2}{2(a+b)},$$ (and the inequality is strict when $a \neq b$), so when $(\sqrt{a}-\sqrt{b})^2 < 1$ and $t$ and $n$ are sufficiently large, the bound \eqref{eq:ourresontiedtree} derived from Theorem \ref{thm:mainseriesparallelbound} is tighter than the bound \eqref{eq:yuryontiedtree} implied by Theorem \ref{thm:yurytheorem} from \cite{yury_percolation}.

Further, applying a union bound on the paths in order to upper-bound the connection probability in the percolation of \cite{yury_percolation}, we have:
\begin{align}I_2(X_u; X_v | Y_{E(G)}) &\leq \mathrm{conn}_{G,\gamma}(u,v) \\ &\leq \sum_{P \in \cP_G(u,v)} \gamma^{|E(P)|} \\ &=
d^t \cdot \gamma^t \\ &= \left((\sqrt{a}-\sqrt{b})^2 + o_n(1)\right)^t \label{last22}
\end{align}
Therefore, taking now $$(\sqrt{a}-\sqrt{b})^2>1>\frac{(a-b)^2}{2(a+b)},$$
and $t$ such that $t=O(n)$ and $t = \omega_n(1)$ (e.g., $t=n$), we have that \eqref{last1} vanishes as $n$ diverges, and our bound implies that $u$ and $v$ cannot be synchronized non-trivially (a.k.a.\ reconstruction/weak recovery is not solvable). In contrast, \eqref{last22} blows up as $n$ diverges, i.e., \eqref{last22} tends to 1 (as one can always cap the upper-bound to 1).

\section{Proof of Theorem \ref{thm:mainseriesparallelbound}}
The result is obtained by using (i) a multiplicative property of the Chi-squared mutual information on paths (see Proposition \ref{prop:chi2seriesmultiplicativity}), and (ii) a subadditivity property of the Chi-squared mutual information on depth-1 trees (see Lemma \ref{lem:chi2parallelsubadditivity}). Interestingly, (i) does not hold for the classical mutual information, $I_{KL}$, making the Chi-squared mutual information, $I_2$, a natural choice for this proof:

\begin{proof}[Proof (of Theorem \ref{thm:mainseriesparallelbound})]
In the following, we implicitly use $I_2(X_kX_l; Y) = I_2(X_k; X_l | Y)$, by Proposition \ref{prop:i2equalsprod}.

The proof is by induction on $|E(G)|$. The base case, $|E(G)| = 1$, is trivial.
For the inductive step, one of two cases holds:

(Case 1) $G$ is the series composition of $H_1$ which is series-parallel with terminals $u,w$, and $H_2$, which is series-parallel with terminals $w,v$. \begin{align*}&I_2(X_u X_v ; Y_{E(G)}) = I_2((X_u X_w) \cdot (X_w X_v) ; Y_{E(H_1)}, Y_{E(H_2)}) \\ &= I_2(X_u X_w; Y_{E(H_1)}) I_2(X_w X_v; Y_{E(H_2)}) \tag*{(Prop. \ref{prop:chi2seriesmultiplicativity})} \\ &\leq \textstyle\sum_{\substack{P_1 \in \cP_{H_1}(u,w) \\ P_2 \in \cP_{H_2}(w,v)}} I_2(X_uX_w; Y_{E(P_1)}) I_2(X_wX_v; Y_{E(P_2)}) \\ &= \textstyle\sum_{\substack{P_1 \in \cP_{H_1}(u,w) \\ P_2 \in \cP_{H_2}(w,v)}} I_2(X_uX_v; Y_{E(P_1)}, Y_{E(P_2)}) \tag*{(Prop. \ref{prop:chi2seriesmultiplicativity})} \\ &= \textstyle \sum_{P \in \cP_{G}(u,v)} I_2(X_uX_v; Y_{E(G)}). \end{align*}
The inequality is by the inductive hypothesis.

(Case 2) $G$ is the parallel composition of $H_1$ and $H_2$ both series-parallel, with terminals $u,v$.
Then, \begin{align*}I_2&(X_uX_v; Y_{E(G)}) = I_2(X_uX_v; Y_{E(H_1)}, Y_{E(H_2)}) \\ &\leq I_2(X_uX_v; Y_{E(H_1)}) + I_2(X_uX_v; Y_{E(H_2)}) \tag*{(Lem. \ref{lem:chi2parallelsubadditivity})} \end{align*}
The inductive step follows by the inductive hypothesis, since $\cP_G(u,v) = \cP_{H_1}(u,v) \sqcup \cP_{H_2}(u,v)$.
\end{proof}

\section{Future directions}\label{sec:futuredirections}

\paragraph{Connections with correlation decay} As mentioned in the introduction, fixing the edge observations and applying Ising model correlation decay conditions yields bounds that are not as strong as those we proved in this paper, because the techniques in our paper allow us to deal with the average-case edge observations, while fixing the edge observations and applying the Dobrushin conditions requires us to work with the worst-case edge observations. It would nonetheless be interesting to elaborate on this connection.

\paragraph{Information subadditivity for general graphs}
A possible extension would be to prove that Theorem \ref{thm:mainseriesparallelbound} applies to spin synchronization models on general graphs $G$, that have general asymmetric edge channels.
\begin{conjecture}[Generalization of Theorem \ref{thm:mainseriesparallelbound}] \label{conj:mainseriesparallelboundgeneralization}
Let $(X,Y)$ be a spin synchronization model on {\em any} graph $G$, and arbitrary edge channels $Q$. Then, for any $u,v \in V(G)$, $$I_2(X_u; X_v \mid Y_{E(G)}) \leq \sum_{P \in \cP_G(u,v)} I_2(X_u; X_v \mid Y_{E(P)}).$$ Here $\cP_G(u,v)$ is the set of paths (self-avoiding walks) from $u$ to $v$ in $G$.
\end{conjecture}
In particular, Conjecture \ref{conj:mainseriesparallelboundgeneralization} would  directly imply the two-community stochastic block-model information-theoretic threshold for impossibility of weak recovery, without having to reduce to the broadcasting on trees model.

Another related subadditivity bound that one could hope for is:
\begin{conjecture}\label{conj:conjecturedsubadditivityvariation} Let $(X,Y)$ be a spin synchronization model on a graph $G$, and arbitrary edge channels $Q$. Then, for any $u \in V(G)$, $W \subset V(G)$,
\begin{equation*}I(X_u; X_W | Y_{E(G)}) \leq \sum_{w \in W} I(X_u; X_w | Y_{E(G)}),\end{equation*} where $I$ is either the Chi-squared mutual information, $I_2$, or the KL mutual information, $I_{KL}$.
\end{conjecture}
Applied to spin synchronization models on trees, Conjecture \ref{conj:conjecturedsubadditivityvariation} would (also) directly imply the information subadditivity inequality \eqref{eq:evanssubadditivity} from \cite{evans} for broadcasting on trees.

\paragraph{Information subadditivity for non-uniform binary spins}
The information subadditivity inequality of Theorem \ref{thm:mainseriesparallelbound} does not hold when the spin synchronization model is generalized to allow for non-uniform binary spins.

For example, consider a graph $G$ with two nodes $u,v$, and two parallel edges $e,f$ between the nodes. Suppose we have vertex labels $X_u, X_v \stackrel{i.i.d}{\sim} \Rad(\delta)$, and edge labels $Y_e = X_uX_v Z_e$ and $Y_f = X_uX_v Z_f$, where $Z_e,Z_f \stackrel{i.i.d} \sim \Rad(\epsilon)$. Then, by direct calculation
$$I_2(X_u; X_v | Y_e, Y_f) = \frac{\delta^2(1-\delta)^2(1-2\eps)^2 (\eps^2 + (1-\eps)^2)^3}{(-4\delta^2\eps^2+4\delta^2\eps-\delta^2+4\delta\eps^2-4\delta\eps+\delta+\eps^4-2\eps^3+\eps^2)^2},$$ and $$I_2(X_u; X_v | Y_e) = I_2(X_u; X_v | Y_f) = \frac{\delta^2(1-\delta)^2(1-2\eps)^2}{(-4\delta^2\eps^2+4\delta^2\eps-\delta^2+4\delta\eps^2-4\delta\eps+\delta-\eps^2+\eps)^2}.$$

When $\delta = 0.2$ and $\eps = 0.2$, one may check that $I_2(X_u; X_v | Y_e, Y_f) \geq I_2(X_u; X_v | Y_e) + I_2(X_u; X_v | Y_f)$, contradicting subadditivity.

This leaves open the question of whether there exists a natural generalization of the subadditivity theorem that holds when the vertex labels are non-uniform (possibly with additional factors).

\paragraph{Information subadditivity for general alphabets} Another way to generalize the vertex labels of the spin synchronization model is to draw them from larger alphabets than the binary alphabet. For example, as in \cite{our_grid}, one may consider vertex labels $X_v$ that are uniformly chosen from a group $\cG$, and edge labels $Y_{uv}$ that are noisy observations of the difference of the endpoints, $X_vX_u^{-1}$.

However, the subadditivity of the Chi-squared information over paths does not extend in a direct way to these larger alphabets. For example, consider a ``spoon''-shaped graph $G$ with vertices $u,v,w$, an edge $e$ between $u$ and $v$, and two parallel edges $f_1,f_2$ between $v$ and $w$. Let the vertex labels $X_u,X_v,X_w$ be i.i.d uniform in $\ZZ/4\ZZ$, and let the edge labels be as follows:
$$Y_e = \begin{cases} 0, & \mbox{if } X_u - X_v \in \{0,1\} \\ 1, & \mbox{if } X_u - X_v \in \{2,3\} \end{cases}, Y_{f_1} = \begin{cases} 0, & \mbox{if } X_v - X_w \in \{0,1\} \\ 1, & \mbox{if } X_v - X_w \in \{2,3\} \end{cases},Y_{f_2} = \begin{cases} 0, & \mbox{if } X_v - X_w \in \{0,2\} \\ 1, & \mbox{if } X_v - X_w \in \{1,3\} \end{cases}.$$

Notice that $I_2(X_u; X_w | Y_e, Y_{f_1}, Y_{f_2}) = I_2(X_u; X_v | Y_e) = 1$. Also, $I_2(X_u; X_w | Y_e, Y_{f_1}) = 1/2$ and $I_2(X_u; X_w | Y_e, Y_{f_2}) = 0$. Therefore, the subadditivity over paths $I_2(X_u; X_w | Y_e, Y_{f_1}, Y_{f_2}) \leq I_2(X_u; X_w | Y_e, Y_{f_1}) + I_2(X_u; X_w | Y_e, Y_{f_2})$ does not hold.

Nevertheless, it is an interesting problem to find a generalization of information subadditivity when the vertex labels are uniform over groups. For example, it is plausible that the scaling of the information-theoretic threshold for weak recovery in the $k$-community SBM could be recovered from such a generalization.

\paragraph{A converse to Theorem \ref{thm:mainsymmetricbound}}
When the edge channels of the spin synchronization model are symmetric, Theorem \ref{thm:mainsymmetricbound} from \cite{abbe2018information} and Theorem \ref{thm:yurytheorem} from \cite{yury_percolation} are tight on trees, so one cannot open the edges with lower probability in general. Is there a converse to Theorem \ref{thm:mainsymmetricbound}: i.e., is the mutual information lower-bounded by the connection probability on some non-trivial percolation? For example, for some bounded-degree graphs?

\appendices

\section{Properties of Chi-squared mutual information}\label{app:chi2mutualinformation}

\begin{definition} For two random variables $A,B$ with joint distribution $\nu_{A,B}$, the Chi-squared mutual information $I_2(A;B)$ is the $f$-mutual information $I_f(A; B) := D_f(\nu_{A,B} || \nu_A \times \nu_B)$ for the choice $f(t) = (t-1)^2$. The KL mutual information is the $f$-mutual information with the choice $f(t) = t \log_2 t$. Here, $D_f(\mu || \nu) := \int f\left(\frac{d\mu}{d\nu}\right) d\nu$ is the $f$-divergence introduced in \cite{csiszar1967information}, which is well-defined and has desirable properties such as nonnegativity and monotonicity when $\mu \ll \nu$, $f$ is convex, $f$ is strictly convex at 1, and $f(1) = 0$. In particular, the $f$-mutual informations have a data processing inequality.
\end{definition}

\begin{proposition}\label{prop:correqualschi2}
Let $A,U$ be jointly-distributed random variables, with $U \in \{-1,+1\}$. Then $$I_2(A; U) = \frac{\Var[\EE[U | A]]}{\Var[U]}.$$ In particular, if $U \sim \Rad(1/2)$, then $$I_2(A; U) = \EE[\EE[U | A]^2].$$
\end{proposition}
\begin{proof}
Letting $\nu_Z$ denote the distribution of $Z$, and $\Omega$ denote the sample set of $A$,
\begin{align}

        I_2(A; U) &= \int_{\Omega \times \{-1,+1\}} \left(\frac{d(\nu_{A,U})}{d(\nu_A \times \nu_U)} - 1\right)^2 d(\nu_{A} \times \nu_{U}) \tag*{} \\ &= \int_{\Omega} \sum_{u \in \{-1,+1\}} \nu_U(u) \left(\frac{1}{ \nu_U(u)} \cdot \frac{d(\nu_{A,U}(\cdot, u))}{d\nu_A(\cdot)} - 1\right)^2 d\nu_A \tag*{} \\ &= \int_{\Omega} \sum_{u \in \{-1,+1\}} \frac{1}{ \nu_U(u)} \left(\frac{d(\nu_{A,U}(\cdot, u))}{d\nu_A(\cdot)} - \nu_U(u)\right)^2 d\nu_A \tag*{}
\end{align}
So, since $$\frac{d(\nu_{A,U}(\cdot, 1))}{d\nu_A(\cdot)} - \nu_U(1) = \left(1 - \frac{d(\nu_{A,U}(\cdot, -1))}{d\nu_A(\cdot)}\right) - (1 - \nu_U(-1)) = -\left(\frac{d(\nu_{A,U}(\cdot, -1))}{d\nu_A(\cdot)} - \nu_U(-1)\right)$$ $\nu_A$-almost everywhere, we have
\begin{align*}
 I_2(A; U) &= \left( \sum_{u \in \{-1,+1\}}\frac{1}{ \nu_U(u)} \right) \cdot \int_{\Omega}  \left(\frac{d(\nu_{A,U}(\cdot, 1))}{d\nu_A(\cdot)} - \nu_U(1)\right)^2 d\nu_A \\ &= \frac{4}{\Var[U]} \cdot \int_{\Omega} \left(\PP[U = 1 | A] - \PP[U = 1]\right)^2 d\nu_A \\ &= \frac{4\Var[\PP[U = 1|A]]}{\Var[U]} = \frac{4\Var[\EE[(U/2) | A]]}{\Var[U]} = \frac{\Var[\EE[U|A]]}{\Var[U]}.
\end{align*}
When $U \sim \Rad(1/2)$, we have $\EE[U] = 0$ and $\Var[U] = 1$, so $I_2(A; U) = \EE[\EE[U | A]^2]$.
\end{proof}

\begin{proposition}\label{prop:i2equalsprod} Let $(X,Y)$ be a spin synchronization model. Then $I_2(X_uX_v; Y) = I_2(X_u; X_v, Y)$.
\end{proposition}
\begin{proof} \begin{align*}I_2(X_u; X_v, Y) &= I_2(X_uX_v; X_v, Y) \tag{Data-processing} \\ &= I_2(X_uX_v; Y) \tag{Since $X_uX_v, Y \indy X_v$}.\end{align*}
\end{proof}

\begin{proposition}\label{prop:i2conditionequalsi2side}Let $(X,Y)$ be a spin synchronization model. Then $I_2(X_u; X_W, Y) = I_2(X_u; X_W | Y).$
\end{proposition}
\begin{proof}
\begin{align*}I_2(X_u; X_W, Y) &= \EE[\EE[X_u | X_W, Y]^2] \tag{Prop. \ref{prop:correqualschi2}} \\ &= \EE[\EE[\EE[X_u | X_W, Y]^2 | Y]] \\ &= I_2(X_u; X_W | Y) \tag{Prop. \ref{prop:correqualschi2}, since $X_u \indy Y$, so $X_u | Y \sim \Rad(1/2)$} \end{align*}
\end{proof}

\begin{proposition}\label{prop:evanschi2vsklonradhalf}
Let $A,U$ be joint random variables, $U \sim \Rad(1/2)$. Then $$\frac{1}{2}I_2(A; U) \leq I_{\mathrm{KL}}(U; A) \leq I_2(U; A),$$ where $I_{KL}$ is in bits.
\end{proposition}
\begin{proof}
As shown in \cite{evans}, this follows from the inequalities $$\frac{x^2}{2} \leq \frac{1+x}{2} \log_2(1-x) + \frac{1-x}{2} \log_2(1+x) \leq x^2.$$
\end{proof}

\begin{proposition}\label{prop:chi2seriesmultiplicativity}
Let $U,V,W \stackrel{i.i.d}{\sim} \Rad(1/2)$. Let $A$ be the output of a channel on $UW$, and let $B$ be the output of a channel on $WV$. Then $$I_2( UV; A,B) = I_2(UW; A)I_2(WV; B).$$
\end{proposition}
\begin{proof}
\begin{align*}I_2(UV; A,B) &= \EE[\EE[UV|A,B]^2] \tag{Prop. \ref{prop:correqualschi2}} \\ &= \EE[\EE[UWWV| A,B]^2] \tag{Since $W^2 = 1$}\\
&= \EE[\EE[UW|A,B]^2\EE[WV|A,B]^2] \tag{Using $UW \indy WV | A,B$}\\
&= \EE[\EE[UW | A]^2\EE[WV | B]^2] \tag{Using $UW \indy B | A$ and $WV \indy A | B$} \\ &= \EE[\EE[UW | A]^2] \EE[\EE[WV | B]^2] \tag{Using $A \indy B$ because $U,V \stackrel{i.i.d}{\sim} \Rad(1/2)$.} \\ &= I_2(UW; A)I_2(WV; B) \tag{Prop. \ref{prop:correqualschi2}}.\end{align*}
\end{proof}

\begin{lemma}\label{lem:chi2parallelsubadditivity}
Let $U \sim \Rad(1/2)$. Let $A$ and $B$ be the outputs of two independent channels on $U$. Then $I_2(U; A,B) \leq I_2(U; A) + I_2(U; B).$
\end{lemma}
\begin{proof}Let $\nu_{A,B,U}$ denote the joint distribution of $A,B,U$. For simplicity, we prove the lemma when $A,B$ are discrete. For any two random variables $C,D$ with joint distribution $\nu_{C,D}$, $I_2(C; D) = D_{(1-1/t)}(\nu_C\nu_D||\nu_{C,D})$, where $D_{(1/1-t)}$ is the $(1-1/t)$-divergence, so\begin{align}
&I_2(U; A,B) - (I_2(U; A) + I_2(U; B)) \tag*{} \\ &= \EE[(\frac{\nu_{A,B,U}}{\nu_{A,B}\nu_U}-1) - (\frac{\nu_{A,U}}{\nu_A\nu_U}-1)-(\frac{\nu_{B,U}}{\nu_B\nu_U}-1)] \tag*{} \\ &= \EE[(\frac{\nu_{U|A,B}}{\nu_{U|A}}-1)(\frac{\nu_{A,U}}{\nu_A\nu_{U}} - 1)] \label{eq:term1} \\
&~~~~~~~+\EE[(\frac{\nu_{U|A,B}}{\nu_{U|A}}-1) - (\frac{\nu_{B,U}}{\nu_{B}\nu_{U}} - 1)] \label{eq:term2}
\end{align}

We claim Terms \eqref{eq:term1} and \eqref{eq:term2} are $\leq 0$, which implies the lemma statement.

We rewrite Term \eqref{eq:term1}, using $\frac{\nu_{A,U}(a,u)}{\nu_A(a)\nu_U(u)}-1 = 2 \nu_{U|A}(u|a) - 1 = \EE[U|A=a] \cdot u$:
\begin{align*}
    \mbox{\eqref{eq:term1}} &= \EE[(\frac{\nu_{U|A,B}}{\nu_{U|A}} - 1) \cdot \EE[U|A] \cdot U] \\ &= \EE[\EE[U|A] \cdot (U \cdot \frac{\nu_{U|A,B}}{\nu_{U|A}} - U)] \\ &= \EE[\EE[U|A] \cdot (U\cdot \frac{\nu_{U|A,B}}{\nu_{U|A}} - \EE[U|A])].
\end{align*}
Define $$t_a := \sum_{b} \frac{\nu_{B,U|A}(b,1|a)\nu_{B,U|A}(b,-1|a)}{\nu_{B,U|A}(b,1|a) + \nu_{B,U|A}(b,-1|a)}.$$

Note \begin{equation}\label{ineq:subadditivityofteaay}0 \leq t_a \leq \nu_{U|A}(1|a) \nu_{U|A}(-1|a)\end{equation} by the subadditivity of $f(x,y) = xy/(x+y)$ for $x,y \geq 0$. (In particular, for all $a,b,c,d \geq 0$, $f(a,b) + f(c,d) \leq f(a+c,b+d)$.) So

\begin{align*}&\EE[U \cdot \frac{\nu_{U|A,B}}{\nu_{U|A}} | A] = \sum_{u} \frac{u}{\nu_{U|A}} \sum_{b} \frac{\nu_{B,U|A} \nu_{B,U|A}}{\nu_{B|A}} \\ &= \sum_{u} \frac{u}{\nu_{U|A}} \sum_{b} \left(\frac{(\nu_{B|A} - (\nu_{B|A} -\nu_{B,U|A}))\nu_{B,U|A}}{\nu_{B|A}}\right) \\ &= \sum_{u} \frac{u}{\nu_{U|A}} (-t_A + \sum_{b} \nu_{B,U|A}) = \sum_{u} \frac{u}{\nu_{U|A}} (\nu_{U|A} - t_A) \\ &= -\sum_{u} u \frac{t_A}{\nu_{U|A}} \\ &= -\frac{t_A}{\nu_{U|A}(1|A)} + \frac{t_A}{\nu_{U|A}(-1|A)} \\ &= \frac{t_A (\nu_{U|A}(1|A) - \nu_{U|A}(-1|A))}{\nu_{U|A}(1|A)\nu_{U|A}(-1|A)} \\ &=  \frac{t_A}{\nu_{U|A}(1|A)\nu_{U|A}(-1|A)} \cdot \EE[U|A] \\
&= c_A\EE[U|A],\end{align*} for some $0 \leq c_A \leq 1$ by \eqref{ineq:subadditivityofteaay}. Thus, $\eqref{eq:term1} = \EE[U|A]^2 (c_A - 1) \leq 0$, as desired.

Now we bound Term \eqref{eq:term2}.
\begin{align*}
    \mbox{\eqref{eq:term2}} &= \EE[(\frac{\nu_{B,U}}{\nu_{B|A}\nu_U} - 1) - (\frac{\nu_{B,U}}{\nu_B\nu_{U}} - 1)] \tag{Using $B \indy A | U$}\\ &= \EE[(\frac{\nu_A\nu_{B}}{\nu_{A,B}}-1)(\frac{\nu_{B,U}}{\nu_B\nu_U} - 1)] \tag{Since $\EE[\frac{\nu_B}{\nu_{B|A}} - 1] = 0$} \\
&= \sum_{a,b,u} \nu_{A,B,U} \left(\frac{\nu_{B}}{\nu_{B \mid A}} - 1\right)\left(\frac{\nu_{B \mid U}}{\nu_{B}} - 1\right)
\end{align*}
For compactness, write $\alpha_a = \nu_{A|U}(a|1), \beta_a = \nu_{A|U}(a|-1), \gamma_b = \nu_{B|U}(b|1), \delta_b = \nu_{B|U}(b|-1)$:
\begin{align*}
\eqref{eq:term2} &= \sum_{a,b} \frac{\alpha_a\gamma_b}{2} \left(\frac{(\gamma_b + \delta_b)/2}{(\alpha_a\gamma_b + \beta_a\delta_b)/(\alpha_a + \beta_a)} - 1\right)\left(\frac{\gamma_b}{(\gamma_b + \delta_b)/2} - 1\right) \\ &~~~~~~~~~~~+ \frac{\beta_a\delta_b}{2} \left(\frac{(\gamma_b + \delta_b)/2}{(\alpha_a\gamma_b + \beta_a\delta_b)/(\alpha_a + \beta_a)} - 1\right)\left(\frac{\delta_b}{(\gamma_b + \delta_b)/2} - 1\right) \\
&= \sum_{a,b} \frac{\alpha_a\gamma_b}{2} \left(\frac{(\gamma_b + \delta_b)(\alpha_a + \beta_a)}{2(\alpha_a\gamma_b + \beta_a\delta_b)} - 1\right)\left(\frac{\gamma_b - \delta_b}{\gamma_b + \delta_b}\right) + \frac{\beta_a\delta_b}{2} \left(\frac{(\gamma_b + \delta_b)(\alpha_a + \beta_a)}{2(\alpha_a\gamma_b + \beta_a\delta_b)} - 1\right)\left(\frac{\delta_b - \gamma_b}{\gamma_b + \delta_b}\right) \\
&= \sum_{a,b}  \left(\frac{\gamma_b - \delta_b}{\gamma_b + \delta_b}\right)\left(\frac{\alpha_a\gamma_b}{2}\left(\frac{(\gamma_b + \delta_b)(\alpha_a + \beta_a)}{2(\alpha_a\gamma_b + \beta_a\delta_b)} - 1\right) - \frac{\beta_a\delta_b}{2} \left(\frac{(\gamma_b + \delta_b)(\alpha_a + \beta_a)}{2(\alpha_a\gamma_b + \beta_a\delta_b)} - 1\right)\right) \\
&= \sum_{a,b} \left(\frac{\gamma_b - \delta_b}{\gamma_b + \delta_b}\right)\left(\frac{\alpha_a\gamma_b - \beta_a\delta_b}{2}\right) \left(\frac{(\gamma_b + \delta_b)(\alpha_a + \beta_a)}{2(\alpha_a\gamma_b + \beta_a\delta_b)} - 1\right) \\
&= \sum_{a,b} \left(\frac{\gamma_b - \delta_b}{\gamma_b + \delta_b}\right)\left(\frac{\alpha_a\gamma_b - \beta_a\delta_b}{4}\right) \left(\frac{-(\alpha_a - \beta_a)(\gamma_b - \delta_b)}{\alpha_a\gamma_b + \beta_a\delta_b}\right) \\
&= \sum_{a,b} \left(\frac{(\gamma_b - \delta_b)^2}{4(\gamma_b + \delta_b)}\right) \left(-(\alpha_a - \beta_a)\frac{\alpha_a\gamma_b - \beta_a\delta_b}{\alpha_a\gamma_b + \beta_a\delta_b}\right) \\
&= \sum_{b} \left(-\frac{(\gamma_b - \delta_b)^2}{4(\gamma_b + \delta_b)}\right) \sum_{a} \left((\alpha_a - \beta_a)\frac{\alpha_a\gamma_b - \beta_a\delta_b}{\alpha_a\gamma_b + \beta_a\delta_b}\right).
\end{align*} We conclude by using \begin{equation*}\sum_a (\alpha_a - \beta_a) = \sum_a \nu_{A|U}(a|1) - \nu_{A|U}(a|-1) = 1 - 1 = 0,\end{equation*} so \begin{align*}\sum_a (\alpha_a - \beta_a) (\frac{\alpha_a \gamma_b - \beta_a \delta_b}{\alpha_a \gamma_b + \beta_a \delta_b}) &= \sum_a (\alpha_a - \beta_a) (\frac{\alpha_a \gamma_b - \beta_a \delta_b}{\alpha_a \gamma_b + \beta_a \delta_b} + 1) \\ &= \sum_a (\alpha_a - \beta_a) (\frac{2\alpha_a\gamma_b}{\alpha_a \gamma_b + \beta_a \delta_b}) \\ &\geq \sum_a (\alpha_a - \beta_a) (\frac{2\gamma_b}{\gamma_b + \delta_b}) \tag{The inequality is term-wise} \\ &= 0.\end{align*}

Since $\left(-\frac{(\gamma_b - \delta_b)^2}{4(\gamma_b + \delta_b)}\right) \leq 0$ for all $b$, this proves that Term $\eqref{eq:term2} \leq 0$.

\end{proof}

\section{Proof of Theorem \ref{thm:mainsymmetricbound} and Corollary \ref{cor:subadditivitybound}}\label{sec:thmmainsymmetricboundproof}
For the sake of completeness, we provide here a compact proof of the main result of \cite{abbe2018information}.

\begin{proof}[Proof (of Theorem \ref{thm:mainsymmetricbound})]
~
\\
\indent \textit{(Step 1) Reduce to case $|W| = 1$}: Construct the graph $G'$ by adding a new vertex $v$ adjacent to all $w \in W$, letting $X_v \sim \Rad(1/2)$ and defining $Y_{vw} = X_vX_w$ for all $w \in W$. Now $I_2(X_u; X_v | Y_{E(G')}) = I_2(X_u; X_W | Y_{E(G')}) = I_2(X_u; X_v | Y_{E(G)})$, $\mathrm{conn}_{G',\gamma'}(u,v) = \mathrm{conn}_{G,\gamma}(u,W)$, and $(X,Y)$ is a symmetric spin synchronization instance on $G'$.
Hence, it suffices to prove the bound
$I_2(X_u; X_W | Y) \leq \mathrm{conn}_{G,\gamma}(u,W)$ in the case $|W| = 1$.

\textit{(Step 2) Reduce to BSC case}:
For each edge channel $Q_e(\cdot | \pm 1)$, let $T_e$ be the symmetry transformation and define $Z_e = \{Y_e, T_e(Y_e)\}$. By the symmetry property of the edge channels, $X \indy Z$. So $\cL(X,Y|Z)$, the law of $(X,Y)$ conditioned on $Z$, is almost surely the law of a spin synchronization model $(X',Y')$ on $G$, where each of the channels $Q'_e$ is binary-valued (either $Y_e$ or $T_e(Y_e)$). Explicitly, for $z \in Z_e$, \begin{align*}Q'_e(z | +1) &= \frac{dQ_e(z | +1)}{d(Q_e(z | +1) + Q_e(T_e(z) | +1))},\end{align*} and by the symmetry property this equals \begin{align*}\frac{dQ_e(T_e(z) | -1)}{d(Q_e(T_e(z) | -1) + Q_e(z | -1))} = Q'_e(T_e(z) | -1).\end{align*} So $Q'_e$ is a binary symmetric channel.
Hence, proving Theorem \ref{thm:mainsymmetricbound} when all the edge channels are BSC yields the general bound: \begin{align*}I_2(X_u; X_v | Y) &= \EE_Z[I_2(X_u; X_v | Y,Z)] \\ &= \EE_Z[I_2(X'_u; X'_v | Y', Z)] \\ &\leq \EE_Z[\mathrm{conn}_{G,\gamma_Z}(u,v)] \tag*{(by BSC case)} \\ &= \mathrm{conn}_{G,\gamma}(u,v) \hspace{1.42cm} \mbox{(since $\gamma = \EE_Z[\gamma_Z]$)}.\end{align*} Here, $\gamma_Z((i,j)) = I_2(X_i, X_j | Y_{ij},Z_{ij}),$ and for the last equality we use the fact that $\gamma((i,j)) = I_2(X_i, X_j | Y_{ij}) = \EE_Z[\gamma_Z((i,j))]$.

\textit{(Step 3) Prove BSC case}: Now it only remains to prove $I_2(X_u; X_v | Y) \leq \mathrm{conn}_{G,\gamma}(u,v)$ when all edge channels $Q_e$ are BSC. Let the flip probability of $Q_e$ be $\varepsilon(e)$, and define $\delta(e) = (1-2\varepsilon(e))$. We can assume that $\delta(e) \in [0,1]$, because we lose no information by flipping edge labels deterministically. Also, by direct calculation $\gamma(e) = \delta(e)^2.$

The proof goes by induction on $|S_{\delta}|$, where $$S_{\delta} := \{e \in E(G) : \delta(e) \not\in \{0,1\}\}.$$

In the base case, $|S_{\delta}| = 0$, so all edge observations are completely noiseless or completely noisy. Hence, $I_2(X_u, X_v | Y) = 1$ if there is a path $P$ from $u$ to $v$ whose edges are all noiseless. If there is no such path, then $I_2(X_u, X_v | Y) = 0$. This is exactly the statement $I_2(X_u, X_v | Y) = \mathrm{conn}_{G,\gamma}(u,v)$.

For the inductive step, assume the theorem when the BSC channels are given by $\delta' : E(G) \to [0,1]$ with $|S_{\delta'}| < |S_{\delta}|$. Pick an arbitrary edge $f \in S_{\delta}$. We will now interpolate between the case in which $\delta(f) = 0$, and the case in which $\delta(f) = 1$, with the other edge channels held fixed. For any $t \in [0,1]$, let $\delta_t : E(G) \to [0,1]$ be given by $\delta_t(e) = \delta(e)$ for $e \neq f$, and $\delta_t(f) = t$. Define corresponding spin synchronization models $(X_t, Y_t)$, and also $\gamma_t = \delta_t^2$. Write $$I(t) := I_2(X_{t,u}; X_{t,v} | Y_t), \mbox{ and } C(t) := \mathrm{conn}_{G,\gamma_t}(u,v).$$ In order to prove that $I(t) \leq C(t)$ for all $t \in [0,1]$, we need the following claim:
\begin{claim}\label{claim:termwisegrowthbound} There is non-decreasing $h : [0,1] \to \RR$ such that
\begin{equation*}\label{eq:ggrowthbound}I(t) = I(0) + (I(1) - I(0)) \cdot t^2 \cdot h(t),\end{equation*}
\end{claim}
Assume the claim is true. Then, $h(1) = 1$, and since $h(t)$ is non-decreasing, $h(t) \leq 1$. Hence, \begin{align*}I(t) &\leq I(0) + (I(1) - I(0)) \cdot t^2 \\ &= I(0) \cdot (1-t^2) + I(1) \cdot t^2 \\ &\leq C(0) \cdot (1-t^2) + C(1) \cdot t^2 = C(t).\end{align*} The inequality of the last line follows because $I(0) \leq C(0)$ and $I(1) \leq C(1)$ by the inductive hypothesis. The equality of the last line follows by the linearity of the connection probability in the parameter $\gamma_t(f) = t^2$.

It only remains to prove the claim. Write $E' = E(G) \sm f$. Also write $f = (i,j)$, $A_t = X_{t,u} \cdot X_{t,v}$, and $B_t = X_{t,i} \cdot X_{t,j}$. By Proposition \ref{prop:correqualschi2}, and because $Y_{t,E'}$ is a subset of $Y_t$, $$I(t) = \EE[\EE[A_t | Y_t]^2] = \EE[\EE[\EE[A_t | Y_t]^2 | Y_{t,E'}]].$$

Since the only edge channel to change with $t$ is $Q_f$, we can couple $X_0 = X_t$, and $Y_{0,E'} = Y_{t,E'}$.
Hence, it suffices to prove that the function $$h(t; Y_{0,E'}) := \frac{1}{t^2}\left(\EE[\EE[A_t | Y_t]^2 | Y_{t,E'}] - \EE[\EE[A_t | Y_0]^2 | Y_{0,E'}]\right)$$ is non-decreasing in $t$, since then we can set $h(t) = \sum_{\sigma \in \{-1,+1\}^{E'}} h(t; \sigma) \cdot \PP[Y_{0,E'} = \sigma],$ which will also be non-decreasing in $t$.

Fix $\sigma \in \{-1,+1\}^{E'}$ such that $\PP[Y_{t,E'} = \sigma] > 0$, and let $P_{\alpha\beta} = \PP[(A_t, B_t) = (\alpha,\beta) \mid Y_{t,E'} = \sigma].$ Set $a = P_{1,1}$, $b = P_{1,-1}$, $c = P_{-1,1}$, $d = P_{-1,-1}$. Since $Y_{t,f} \indy A_t,Y_{t,E'} | B_t$, one can explicitly calculate
\begin{align*}\EE[\EE[&A_t|Y_t]^2|Y_{t,E'}] \\ = &\bigg( \frac{((a(1-t) + b(1+t)) - (c(1-t)+d(1+t)))^2}{2((a(1-t) + b(1+t)) + (c(1-t)+d(1+t)))} \\ &+\frac{((a(1+t) + b(1-t)) - (c(1+t)+d(1-t)))^2}{2((a(1+t) + b(1-t)) + (c(1+t)+d(1-t)))}\bigg),\end{align*}
Plugging this in and simplifying, if $b = d = 0$ or $a = c = 0$, then $h(t; \sigma) = 0$, which is non-decreasing because it is constant. Otherwise we get $h(t; \sigma) = \frac{16(ad - bc)^2}{1-t^2(a-b+c-d)^2},$ which is non-decreasing on $[0,1]$ because $(a-b+c-d)^2 < (a+b+c+d)^2 = 1$.
This proves the claim.

\end{proof}
\begin{lemma}\label{lem:i2onpath}
Suppose $(X,Y)$ is a spin synchronization model on a path $P$ with endpoints $u$ and $v$. Then $$\mathrm{conn}_{P,\gamma}(u,v) = I_2(X_u; X_v | Y_{E(P)}).$$
\end{lemma}
\begin{proof}
 \begin{align*}&\mathrm{conn}_{P,\gamma}(u,v) = \textstyle \prod_{(i,j) \in E(P)} I_2(X_i; X_j | Y_{ij}) \\
 &= \textstyle \prod_{(i,j) \in E(P)} I_2(X_i \cdot X_j; Y_{ij}) \tag*{(Prop. \ref{prop:i2equalsprod})} \\
 &= I_2(\textstyle \prod_{(i,j) \in E(P)} X_i \cdot X_j; Y_{E(P)}) \tag*{(Prop. \ref{prop:chi2seriesmultiplicativity})} \\
 &= I_2(X_u \cdot X_v; Y_{E(P)}) = I_2(X_u; X_v | Y_{E(P)}) \tag*{(Prop. \ref{prop:i2equalsprod})}
 \end{align*}
\end{proof}

\begin{proof}[Proof (of Corollary \ref{cor:subadditivitybound})]

The corollary follows from Theorem \ref{thm:mainsymmetricbound}, the union bound $\mathrm{conn}_{G,\gamma}(u,v) \leq \sum_{P \in \cP_G(u,v)} \mathrm{conn}_{P,\gamma}(u,v),$ and Lemma \ref{lem:i2onpath}.
\end{proof}

\section*{Acknowledgements}
This research was partly supported by NSF CAREER Award CCF-1552131.

\printbibliography

\end{document}