  One of the principal topics of this paper concerns the realization of
self-adjoint operators $L_{\Theta, \Om}$ in $L^2(\Om; d^n x)^m$, $m, n \in
\bbN$, associated with divergence form elliptic partial differential
expressions $L$ with (nonlocal) Robin-type boundary conditions in bounded
Lipschitz domains $\Om \subset \bbR^n$. In particular, we develop the theory in
the vector-valued case and hence focus on matrix-valued differential
expressions $L$ which act as $$ Lu = -
\biggl(\sum_{j,k=1}^n\partial_j\bigg(\sum_{\beta = 1}^m
a^{\alpha,\beta}_{j,k}\partial_k u_\beta\bigg) \bigg)_{1\leq\alpha\leq m},
\quad u=(u_1,...,u_m). $$ The (nonlocal) Robin-type boundary conditions are
then of the form $$ \nu \cdot A D u + \Theta \big[u\big|_{\partial \Om}\big] =
0 \, \text{on $\partial \Om$}, $$ where $\Theta$ represents an appropriate
operator acting on Sobolev spaces associated with the boundary $\partial \Om$
of $\Om$, $\nu$ denotes the outward pointing normal unit vector on
$\partial\Om$, and $Du:=\bigl(\partial_j
u_\alpha\bigr)_{\substack{1\leq\alpha\leq m 1\leq j\leq n}}$.
  Assuming $\Theta \geq 0$ in the scalar case $m=1$, we prove Gaussian heat
kernel bounds for $L_{\Theta, \Om}$ by employing positivity preserving
arguments for the associated semigroups and reducing the problem to the
corresponding Gaussian heat kernel bounds for the case of Neumann boundary
conditions on $\partial \Om$. We also discuss additional zero-order potential
coefficients $V$ and hence operators corresponding to the form sum $L_{\Theta,
\Om} + V$.
