  We introduce a nonlinear potential theory problem for the Laplacian, the
solution of which characterizes the Berezin density $B(z,\cdot)$ for the
polynomial Bergman space, where the point $z\in\mathbb{C}$ is fixed. When
$z=\infty$, the Berezin density is expressed in terms of the squared modulus of
the corresponding normalized orthogonal polynomial $P$. We use an approximate
version of this characterization to study the asymptotics of the orthogonal
polynomials in the context of exponentially varying weights. This builds on
earlier works by Its-Takhtajan and by the first author on a soft
Riemann-Hilbert problem for planar orthogonal polynomials, where in place of
the Laplacian we have the $\bar\partial$-operator.
  We adapt the soft Riemann-Hilbert approach to the nonlinear potential
problem, where the nonlinearity is due to the appearance of $|P|^2$ in place of
$\overline{P}$. Moreover, we suggest how to adapt the potential theory method
to the study of the asymptotics of more general Berezin densities $B(z,w)$ in
the off-spectral regime, that is, when $z$ is fixed outside the droplet. This
is a first installment in a program to obtain an explicit global expansion
formula for the polynomial Bergman kernel, and, in particular, of the one-point
function of the associated random normal matrix ensemble.
