We investigate the possibility of doubly heavy $qq^\prime \bar Q \bar Q^\prime$ tetraquark bound states using $n_f=2+1$ lattice QCD with pion masses $\simeq 164$, $299$ and $415$ MeV. Two types of lattice interpolating operator are chosen, reflecting first diquark-antidiquark and second meson-meson structure. Performing variational analyses using these operators and their mixings, we determine the ground and first excited states from the lattice correlators. Using non-relativistic QCD to simulate the bottom quarks and the Tsukuba formulation of relativistic heavy quarks for charm quarks, we study the $ud\bar b \bar b$, $\ell s\bar b \bar b$ as well as $ud\bar c \bar b$, channels with $\ell=u,d$. In the case of the $ud\bar b \bar b$ and $\ell s\bar b \bar b$ channels unambiguous signals for $J^P=1^+$ tetraquarks are found with binding energies $189(10)$ and $98(7)$ MeV below the corresponding free two-meson thresholds at the physical point. These tetraquarks are therefore not only strong-interaction, but also electromagnetic-interaction stable, and can decay only weakly. So far these are the first exotic hadrons predicted to have this feature. Further evidence for binding is found in the $ud\bar c \bar b$ channel, where the binding energy broadly straddles the electromagnetic stability threshold.
We also study the dependence of the tetraquark binding on heavy quark mass by considering the channels $ud\bar b' \bar b'$, $\ell s\bar b' \bar b'$ as well as $ud\bar b' \bar b$, $\ell s\bar b' \bar b$ involving a heavy $b'$ quark with mass between roughly $0.6$ and $6.3$ times the physical $b$ quark mass. The observed mass dependence of these four flavor channels is shown to follow closely a phenomenological form expected on simple physical grounds.