  Let $\Omega\subset\mathbb{R}^{n+1}$, $n\ge 2$, be a 1-sided non-tangentially
accessible domain (aka uniform domain), that is, $\Omega$ satisfies the
interior Corkscrew and Harnack chain conditions, which are respectively
scale-invariant/quantitative versions of openness and path-connectedness. Let
us assume also that $\Omega$ satisfies the so-called capacity density
condition, a quantitative version of the fact that all boundary points are
Wiener regular. Consider $L_0 u=-\mathrm{div}(A_0\nabla u)$,
$Lu=-\mathrm{div}(A\nabla u)$, two real (non-necessarily symmetric) uniformly
elliptic operators in $\Omega$, and write $\omega_{L_0}$, $\omega_L$ for the
respective associated elliptic measures. The goal of this program is to find
sufficient conditions guaranteeing that $\omega_L$ satisfies an
$A_\infty$-condition or a $RH_q$-condition with respect to $\omega_{L_0}$. In
this paper we are interested in obtaining square function and non-tangential
estimates for solutions of operators as before. We establish that bounded weak
null-solutions satisfy Carleson measure estimates, with respect to the
associated elliptic measure. We also show that for every weak null-solution,
the associated square function can be controlled by the non-tangential maximal
function in any Lebesgue space with respect to the associated elliptic measure.
These results extend previous work of Dahlberg-Jerison-Kenig and are
fundamental for the proof of the perturbation results in arXiv:1901.08261.
