Let $\Omega\subseteq\mathbb R^d$, $u\in\mathcal L^1(\Omega)$ and $$u_\varepsilon(x):=\frac 1{\varepsilon^d}\int_\Omega\rho\left(\frac{x-y}\varepsilon\right)u(y)\;{\rm d}\lambda(y)\;\;\;\text{for }\varepsilon>0\text{ and }x\in\mathbb R^d$$ where $\rho\in C^\infty(\mathbb R^d)$ has support in the unit ball and $$\int\rho\;{\rm d}\lambda=1\;.$$
We can show, that if $u\in\mathcal L^p(\Omega)$ for some $p\in[1,\infty)$, then $u_\varepsilon\in\mathcal L^p(\Omega)$ for all $\varepsilon>0$ and $$\left\|u_\varepsilon-u\right\|_{L^p(\Omega)}\stackrel{\varepsilon\to 0}\to 0\tag 1$$ (see Elliptic Partial Differential Equations of Second Order by Gilbarg and Trudinger, Lemma 7.2).
If $\Omega$ is bounded, then $u_\varepsilon\in C_c^\infty(\mathbb R^d)$ and we can conclude that $C_c^\infty(\Omega)$ is dense in $(L^p(\Omega),\left\|\;\cdot\;\right\|_{L^p(\Omega)})$.
However, I found many lecture notes (as well as questions here on the board) where they state, that this result is true for any (open) subset $\Omega$. So, which assumptions on $\Omega$ do we really need and where can I found a rigorous proof for the result in that case?
It's true for any open $\Omega$: you have to replace $u$ by the function $$ \tilde{u}(x) = \begin{cases} u(x) & x \in \Omega \\ 0 & \text{else} \end{cases} $$ Then you get bounds on the $L^p(\Omega)$-norms from the $L^p(\mathbb{R}^n)$-norms, and exactly the same proof applies. See, for example, p.64ff. of Lieb and Loss, Analysis, and particularly sections 2.16 and 2.19.