Here is a sketch of the "cheapest" way I know how to prove something like that. Filling in the details may still be a bit lengthy but should be essentially routine. To set things up, let's write $X$ for the space of oriented non-degenerate area $1$ triangles, so $\tau \in X$ is a pair $(v_1,v_2)$ with $v_i \in \mathbf{R}^2$ and $|v_1 \wedge v_2| = 2$. I then view your process as a Markov chain $(\tau_n,p_n)$ on $X \times \textbf{R}^2$ with the second coordinate $p_n$ being the $n$th value $p$ chosen in your description and $\tau_n$ the $n$th triangle obtained in your construction, translated and oriented in such a way that the edge $\overline{p_{n} p_{n+1}}$ of that triangle corresponds to $(v_1)_n$, say. We are then in the general setting that we have an autonomous Markov chain $\tau_n$ on $X$, a function $f\colon X \to \mathbf{R}^2$, and a random walk on $\mathbf{R}^2$ given by $p_{n+1} = p_n + f(\tau_n)$. In this case $f(v_1,v_2) = v_1$, but the argument doesn't care about $f$ and the Markov chain $\tau$ as long as they are "nice" enough as detailed below.

The space $X$ comes with a natural 'size' $D(v_1,v_2) = \max\{|v_1|,|v_2-v_1|\}$ which has compact sublevel sets. It is also not hard to check that, if we denote by $P$ our Markov transition operator on $X$, the function $D$ is such that for every $q > 1$ there exist constants $c \in (0,1)$ and $C > 0$ such that $P(D^q) \le c D^q + C$, i.e. $D^q$ is a Lyapunov function. Even better, for $\kappa >0$ small enough $\exp(\kappa D)$, one even has $P \exp(\kappa D) \le \exp(\kappa' D) + K$ for some $\kappa' < \kappa$. It is also the case that for every compact set $K \subset X$ one can find $n > 0$ and a strictly positive measure $\nu$ on $X$ such that $P^n(x,\cdot) \ge \nu$, uniformly over $x \in K$. Together with the existence of a Lyapunov function, this tells us that $P$ admits a unique invariant measure $\mu$, as well as a spectral gap in spaces of bounded functions weighted by $D^p$. (See Harris's theorem.) By uniqueness and rotation invariance of the problem, $\mu$ must be rotation invariant so $\int f\,d\mu = 0$.

We need two more ingredients inspired from stochastic homogenisation:

- a corrector $\phi\colon X \to \mathbf{R}^2$ solving $P\phi = \phi - f$ (which exists and is bounded by $D$ at infinity by the spectral gap mentioned above)
- a second-order corrector $\psi \colon X \to \mathbf{R}^2 \otimes \mathbf{R}^2$ solving $P \psi = f \otimes f + f \otimes P\phi + P\phi \otimes f - A$, where $A \in \mathbf{R}^2 \otimes \mathbf{R}^2$ is the symmetric matrix such that the average of the right hand side against $\mu$ vanishes and therefore the equation for $\psi$ is solvable.

The matrix $A$ is strictly positive definite, by some algebraic manipulation combined with Cauchy-Schwartz and the non-degeneracy of the Markov chain on $X$; then rotation invariance guarantees that it is a constant multiple of the identity. The constant $A$ turns out to be the diffusion coefficient of the limiting Brownian motion in the heuristic picture.

Let's write $Q$ for the Markov operator of the full Markov chain on $Y = X \times \mathbf{R}^2$, so that $(Qg)(\tau,p) = (Pg)(\tau,p+f(\tau))$ (with the obvious abuse of notation). To show that the process eventually covers the whole plane, note that we can find a countable sequence of open sets $U_n \subset Y$ such that if the chain visits each $U_n$ at least once, then the corresponding triangle process will have covered the plane. (Simply take a dense enough grid $y_n \in \mathbf{R^2}$ and, for each $n$ take $U_n$ to be a small enough ball around $(\Delta, y_n)$ with $\Delta$ the standard equilateral triangle.) In other words, it suffices to show that the process is recurrent. By standard criteria (see for example Section 8.4.2 in the book by Meyn and Tweedie), it is enough to find a function $W \colon Y \to \mathbf R$ with compact sublevel sets such that $Q W < W$ outside of a compact set.

For the $2D$ Brownian motion, a possible choice of such a function is given by $V(p) = (\log(1+p^2))^\alpha$ for any choice $\alpha \in (0,1)$. The idea now is to use the correctors to build a function $W$ that explicits the intuitive fact that our random walk behaves like a Brownian motion with diffusion coefficient $A$. I claim (modulo silly mistakes) that a possible choice is
$$
W(\tau,p) = V(p) + \langle \nabla V(p),\phi(\tau)\rangle + {1\over 2}\langle \nabla^2 V(p), \psi(\tau)\rangle + \exp(\kappa D(\tau))\hat V(p)\;,
$$
where $\hat V(p) = 1/(1+|p|^{5/2})$ for example. This certainly has compact sublevel sets since the middle two terms are dominated by the first and last term outside of a compact set. We also have
\begin{align}
(QW-W)(\tau,p) &= V(p+f(\tau)) + \langle \nabla V(p+f(\tau)),(P\phi)(\tau)\rangle + {1\over 2}\langle \nabla^2 V(p+f(\tau)), (P\psi)(\tau)\rangle + (PD^4)(\tau)\hat V(p+f(\tau)) - W(\tau,p)\\
&\approx {1\over 2} A\Delta V(p) + (P\exp(\kappa D(\tau)))(\tau)\hat V(p+f(\tau))-\exp(\kappa D(\tau))\hat V(p)\\
&\le {1\over 2} A\Delta V(p)
\end{align}
outside of a compact set,
which is precisely what we wanted. Here, to get from the first to the second line, we set $$V(p+f(\tau)) \approx V(p) + \langle \nabla V(p),f(\tau)\rangle + {1\over 2}\langle \nabla^2 V(p),(f\otimes f)(\tau)\rangle$$ and similarly for $\nabla V(p+f(\tau))$, and we exploit the definitions of the correctors $\phi$ and $\psi$ to get a bunch of cancellations. This approximation can be justified by absorbing the error terms into the terms $\Delta V$ and $\exp(\kappa D(\tau))\hat V(p)$.

2more comments