  We demonstrate that for gapped bilayer graphene, the nonlinear nature of the
screening of an external disorder potential and the resulting inhomogeneity of
the electron liquid are crucial for describing the electronic compressibility.
In particular, traditional diagrammatic methods of many-body theory do not
include this inhomogeneity and therefore fail to reproduce experimental data
accurately, particularly at low carrier densities. In contrast, a direct
calculation of the charge landscape via a numerical Thomas-Fermi energy
functional method along with the appropriate bulk averaging procedure captures
all the essential physics, including the interplay between the band gap and the
inhomogeneity.
