[proofplan]
We regularize the given plurisubharmonic function by convolution with a non-negative radial mollifier whose rescalings are supported in balls small enough to stay inside $\Omega$. The sub-mean value property for plurisubharmonic functions gives approximants lying above the original function, and monotonicity of spherical means gives a decreasing family as the smoothing radius decreases. Smoothness follows from differentiating the mollifier, and plurisubharmonicity is preserved because convolution is an average of translates of a plurisubharmonic function. Finally, monotone pointwise convergence together with local integrability gives convergence in $L^1_{\mathrm{loc}}$.
[/proofplan]
[step:Choose a radial mollifier supported inside the distance to the boundary]
Identify $\mathbb{C}^n$ with $\mathbb{R}^{2n}$, and let $\mathcal{L}^{2n}$ denote Lebesgue measure on $\mathbb{C}^n$. Since $\Omega_0 \subset\subset \Omega$, define
\begin{align*}
d_0 := \operatorname{dist}(\overline{\Omega_0}, \mathbb{C}^n \setminus \Omega) > 0.
\end{align*}
Choose a function $\rho \in C_c^\infty(B(0,1))$ satisfying $\rho \geq 0$, $\rho(\zeta) = \rho_0(|\zeta|)$ for a smooth radial profile $\rho_0$, $\rho_0$ non-increasing, and
\begin{align*}
\int_{\mathbb{C}^n} \rho(\zeta)\, d\mathcal{L}^{2n}(\zeta) = 1.
\end{align*}
For $0 < \varepsilon < d_0$, define
\begin{align*}
\rho_\varepsilon: \mathbb{C}^n &\to [0,\infty) \\
\zeta &\mapsto \varepsilon^{-2n}\rho(\zeta/\varepsilon).
\end{align*}
Then $\rho_\varepsilon \in C_c^\infty(B(0,\varepsilon))$ and
\begin{align*}
\int_{\mathbb{C}^n} \rho_\varepsilon(\zeta)\, d\mathcal{L}^{2n}(\zeta) = 1.
\end{align*}
For $0 < \varepsilon < d_0$, define
\begin{align*}
\varphi_\varepsilon: \Omega_0 &\to [-\infty,\infty) \\
z &\mapsto \int_{\mathbb{C}^n} \rho_\varepsilon(\zeta)\,\varphi(z-\zeta)\, d\mathcal{L}^{2n}(\zeta).
\end{align*}
The integral is well-defined because $z-\operatorname{supp}\rho_\varepsilon \subset \Omega$ for every $z \in \Omega_0$. Since $\varphi$ is plurisubharmonic and $\varphi \not\equiv -\infty$ on the domain $\Omega$, the local integrability theorem for plurisubharmonic functions gives $\varphi \in L^1_{\mathrm{loc}}(\Omega)$; in particular, $\varphi$ is finite $\mathcal{L}^{2n}$-almost everywhere on compact subsets of $\Omega$.
[/step]
[step:Show that each convolution is smooth]
Fix $0 < \varepsilon < d_0$. Let $\alpha=(\alpha_1,\dots,\alpha_{2n})$ be a multi-index in real coordinates on $\mathbb{C}^n \cong \mathbb{R}^{2n}$. Since $\rho_\varepsilon \in C_c^\infty(B(0,\varepsilon))$ and $\varphi \in L^1_{\mathrm{loc}}(\Omega)$, differentiation under the integral gives
\begin{align*}
D^\alpha \varphi_\varepsilon(z)
=
\int_{\mathbb{C}^n} D^\alpha \rho_\varepsilon(\zeta)\,\varphi(z-\zeta)\, d\mathcal{L}^{2n}(\zeta),
\end{align*}
for every $z \in \Omega_0$. The integrand is supported in $B(0,\varepsilon)$, and for $z$ varying in a compact subset of $\Omega_0$ the translated supports remain in a compact subset of $\Omega$. Hence the displayed formula depends continuously on $z$. Since $\alpha$ was arbitrary, $\varphi_\varepsilon \in C^\infty(\Omega_0)$.
[/step]
[step:Prove that convolution preserves plurisubharmonicity]
Fix $0 < \varepsilon < d_0$. For each $\zeta \in B(0,\varepsilon)$, define the translate
\begin{align*}
T_\zeta: \Omega_0 &\to \Omega \\
z &\mapsto z-\zeta.
\end{align*}
The function $\varphi \circ T_\zeta$ is plurisubharmonic on $\Omega_0$ because $T_\zeta$ is holomorphic and carries every affine complex line in $\Omega_0$ into an affine complex line in $\Omega$. Let $a \in \Omega_0$, let $v \in \mathbb{C}^n$ satisfy $|v|=1$, and let $R>0$ be such that the affine complex disc
\begin{align*}
\Delta(a,v,R):=\{a+\lambda v: \lambda \in \mathbb{C},\ |\lambda|\leq R\}
\end{align*}
is compactly contained in $\Omega_0$. For each $\zeta \in B(0,\varepsilon)$, the translated disc
\begin{align*}
\Delta(a,v,R)-\zeta=\{a-\zeta+\lambda v: \lambda \in \mathbb{C},\ |\lambda|\leq R\}
\end{align*}
is contained in the affine complex line $a-\zeta+\mathbb{C}v$ and lies in $\Omega$, so the definition of plurisubharmonicity gives the sub-[mean value inequality](/theorems/328) for the restriction of $\varphi$ to that affine line:
\begin{align*}
\varphi(a-\zeta)
\leq
\frac{1}{2\pi}\int_0^{2\pi}\varphi(a-\zeta+Re^{i\theta}v)\,d\mathcal{L}^1(\theta).
\end{align*}
The map $(\theta,\zeta)\mapsto \rho_\varepsilon(\zeta)\varphi(a-\zeta+Re^{i\theta}v)$ is absolutely integrable on $[0,2\pi]\times \operatorname{supp}\rho_\varepsilon$ because $\rho_\varepsilon$ is bounded with compact support and all points $a-\zeta+Re^{i\theta}v$ range over a compact subset of $\Omega$ on which $\varphi\in L^1$. Therefore [Fubini's theorem](/theorems/2961) is applicable after multiplying the displayed inequality by $\rho_\varepsilon(\zeta)\geq0$ and integrating with respect to $d\mathcal{L}^{2n}(\zeta)$. Using $\int_{\mathbb{C}^n}\rho_\varepsilon\,d\mathcal{L}^{2n}=1$, we obtain
\begin{align*}
\varphi_\varepsilon(a)
\leq
\frac{1}{2\pi}\int_0^{2\pi}\varphi_\varepsilon(a+Re^{i\theta}v)\,d\mathcal{L}^1(\theta).
\end{align*}
Thus the restriction of $\varphi_\varepsilon$ to every affine complex line satisfies the sub-[mean value inequality](/theorems/328) on every closed affine disc compactly contained in $\Omega_0$, so $\varphi_\varepsilon$ is plurisubharmonic.
Because $\varphi_\varepsilon$ is smooth, this is equivalently the pointwise Levi semipositivity condition
\begin{align*}
\sum_{k,\ell=1}^{n}
\frac{\partial^2 \varphi_\varepsilon}{\partial z_k\,\partial \overline{z}_\ell}(z)\,\xi_k\,\overline{\xi_\ell}
\geq 0
\end{align*}
for every $z \in \Omega_0$ and every $\xi=(\xi_1,\dots,\xi_n)\in \mathbb{C}^n$.
[/step]
[step:Use spherical means to obtain monotone approximation from above]
Let $m:=2n$, let $\sigma_{m-1}$ denote the surface measure of the Euclidean unit sphere $\partial B(0,1)\subset\mathbb{R}^{m}$, and let $\mathcal{L}^1$ denote one-dimensional Lebesgue measure on $\mathbb{R}$. Define the radial probability density
\begin{align*}
b:(0,1)&\to[0,\infty) \\
s&\mapsto \sigma_{m-1}(\partial B(0,1))\,s^{m-1}\rho_0(s).
\end{align*}
Since $\rho(s\omega)=\rho_0(s)$ for $0<s<1$ and $\omega\in\partial B(0,1)$, the polar-coordinate formula for $\mathcal{L}^{m}$ gives
\begin{align*}
\int_0^1 b(s)\,d\mathcal{L}^1(s)
&=
\int_0^1
\int_{\partial B(0,1)}s^{m-1}\rho(s\omega)\,d\sigma_{m-1}(\omega)\,d\mathcal{L}^1(s) \\
&=
\int_{B(0,1)}\rho(\eta)\,d\mathcal{L}^{m}(\eta)
=1.
\end{align*}
For $z\in\Omega_0$ and $0<r<d_0$, define the spherical mean
\begin{align*}
M_z:(0,d_0)&\to[-\infty,\infty) \\
r&\mapsto \frac{1}{\sigma_{m-1}(\partial B(0,1))}
\int_{\partial B(0,1)}\varphi(z-r\omega)\,d\sigma_{m-1}(\omega).
\end{align*}
Because $\varphi$ is plurisubharmonic on $\Omega$ and $\varphi\in L^1_{\mathrm{loc}}(\Omega)$, the distributional characterization of plurisubharmonic functions applies: the complex Hessian current $(\partial^2\varphi/\partial z_k\partial\overline z_\ell)_{k,\ell}$ is positive semidefinite. Taking its trace gives a positive distribution, and the real distributional Laplacian on $\Omega\subset\mathbb{R}^{m}$ satisfies
\begin{align*}
\Delta_{\mathbb{R}^{m}}\varphi
=
4\sum_{k=1}^{n}\frac{\partial^2\varphi}{\partial z_k\partial\overline z_k}
\geq 0
\end{align*}
in the sense of distributions. Since $\varphi$ is also upper semicontinuous, this is precisely the distributional criterion for $\varphi$ to be subharmonic as a real-valued extended function on $\Omega$. For every $z\in\Omega_0$ and every $0<r<d_0$, the closed ball $\overline{B}(z,r)$ is contained in $\Omega$ by the definition of $d_0$. Therefore the monotonicity theorem for spherical means of subharmonic functions applies to $\varphi$ on $B(z,r)$, including the case where $\varphi$ takes the value $-\infty$ at isolated or non-isolated points. It gives that $r\mapsto M_z(r)$ is non-decreasing on $(0,d_0)$ and that the sub-mean inequality gives $\varphi(z)\leq M_z(r)$ for every $0<r<d_0$.
Using the polar-coordinate substitution $\zeta=\varepsilon s\omega$, under which $d\mathcal{L}^{m}(\zeta)=\varepsilon^m s^{m-1}\,d\sigma_{m-1}(\omega)\,d\mathcal{L}^1(s)$ and $0<s<1$, $\omega\in\partial B(0,1)$, we obtain
\begin{align*}
\varphi_\varepsilon(z)
&=
\int_0^1
\int_{\partial B(0,1)}s^{m-1}\rho(s\omega)\varphi(z-\varepsilon s\omega)\,d\sigma_{m-1}(\omega)\,d\mathcal{L}^1(s) \\
&=
\int_0^1 b(s)M_z(\varepsilon s)\,d\mathcal{L}^1(s).
\end{align*}
If $0<\varepsilon_1<\varepsilon_2<d_0$, then $\varepsilon_1s\leq\varepsilon_2s$ for every $s\in(0,1)$, and the monotonicity of $M_z$ gives $M_z(\varepsilon_1s)\leq M_z(\varepsilon_2s)$. Since $b\geq0$, integration with respect to $b(s)\,d\mathcal{L}^1(s)$ yields
\begin{align*}
\varphi_{\varepsilon_1}(z) \leq \varphi_{\varepsilon_2}(z).
\end{align*}
Also, since $\varphi(z)\leq M_z(\varepsilon s)$ for every $s\in(0,1)$ and $\int_0^1b(s)\,d\mathcal{L}^1(s)=1$, we get
\begin{align*}
\varphi(z) \leq \varphi_\varepsilon(z)
\end{align*}
for every $z\in\Omega_0$ and every $0<\varepsilon<d_0$.
[/step]
[step:Identify the pointwise limit as the original function]
Fix $z \in \Omega_0$. From the preceding step, the spherical mean function
\begin{align*}
M_z:(0,d_0)&\to[-\infty,\infty) \\
r&\mapsto \frac{1}{\sigma_{m-1}(\partial B(0,1))}
\int_{\partial B(0,1)}\varphi(z-r\omega)\,d\sigma_{m-1}(\omega)
\end{align*}
is non-decreasing. The same spherical-mean theorem for subharmonic functions gives the boundary recovery property
\begin{align*}
\lim_{r\downarrow0}M_z(r)=\varphi(z),
\end{align*}
with the limit allowed to be $-\infty$.
Let first $\varphi(z)>-\infty$. Given $\delta>0$, choose $r_\delta\in(0,d_0)$ such that
\begin{align*}
M_z(r)\leq \varphi(z)+\delta
\end{align*}
whenever $0<r<r_\delta$. If $0<\varepsilon<r_\delta$, then $0<\varepsilon s<r_\delta$ for every $s\in(0,1)$, and using the formula from the preceding step gives
\begin{align*}
\varphi_\varepsilon(z)
&=
\int_0^1 b(s)M_z(\varepsilon s)\,d\mathcal{L}^1(s) \\
&\leq
(\varphi(z)+\delta)
\int_0^1 b(s)\,d\mathcal{L}^1(s)
=
\varphi(z)+\delta.
\end{align*}
Together with $\varphi_\varepsilon(z)\geq \varphi(z)$, this gives
\begin{align*}
\lim_{\varepsilon\downarrow0}\varphi_\varepsilon(z)=\varphi(z).
\end{align*}
Now suppose $\varphi(z)=-\infty$. Let $A\in\mathbb{R}$ be arbitrary. Since $M_z(r)\to-\infty$ as $r\downarrow0$, choose $r_A\in(0,d_0)$ such that $M_z(r)\leq A$ whenever $0<r<r_A$. If $0<\varepsilon<r_A$, then
\begin{align*}
\varphi_\varepsilon(z)
&=
\int_0^1 b(s)M_z(\varepsilon s)\,d\mathcal{L}^1(s) \\
&\leq
A\int_0^1 b(s)\,d\mathcal{L}^1(s)
=A.
\end{align*}
Because $A\in\mathbb{R}$ was arbitrary, $\varphi_\varepsilon(z)\to-\infty$. Thus pointwise convergence holds for every $z \in \Omega_0$.
[/step]
[step:Choose a decreasing sequence and prove local $L^1$ convergence]
Choose a sequence $(\varepsilon_j)_{j=1}^\infty$ with
\begin{align*}
d_0>\varepsilon_1>\varepsilon_2>\cdots>0,
\qquad
\lim_{j\to\infty}\varepsilon_j=0,
\end{align*}
and define
\begin{align*}
\varphi_j: \Omega_0 &\to [-\infty,\infty) \\
z &\mapsto \varphi_{\varepsilon_j}(z).
\end{align*}
By the preceding steps, each $\varphi_j$ belongs to $C^\infty(\Omega_0)$, is plurisubharmonic on $\Omega_0$, and satisfies
\begin{align*}
\varphi_1 \geq \varphi_2 \geq \cdots \geq \varphi
\end{align*}
with $\varphi_j(z)\to\varphi(z)$ for every $z\in\Omega_0$.
It remains to prove local $L^1$ convergence. Let $K\subset\Omega_0$ be compact. Choose $\eta>0$ such that
\begin{align*}
K_\eta := \{w\in\Omega : \operatorname{dist}(w,K)\leq \eta\}
\end{align*}
is compactly contained in $\Omega$, and choose $j_0$ such that $\varepsilon_j<\eta$ for all $j\geq j_0$. Since $\varphi$ is plurisubharmonic and $\varphi\not\equiv -\infty$, the local integrability theorem for plurisubharmonic functions gives $\varphi\in L^1(K_\eta)$ as a finite $\mathcal{L}^{2n}$-almost everywhere function. For $j\geq j_0$, Tonelli's theorem applies to the non-negative measurable function $(z,\zeta)\mapsto \rho_{\varepsilon_j}(\zeta)|\varphi(z-\zeta)|$. After interchanging the order of integration, make the translation change of variables $w=z-\zeta$; the Lebesgue measure is translation-invariant, so $d\mathcal{L}^{2n}(w)=d\mathcal{L}^{2n}(z)$. Since $\rho_{\varepsilon_j}$ is supported in $B(0,\varepsilon_j)$ and $\varepsilon_j<\eta$, the translated set $K-\zeta$ is contained in $K_\eta$ whenever $\rho_{\varepsilon_j}(\zeta)\neq0$. Thus
\begin{align*}
\int_K |\varphi_j(z)|\,d\mathcal{L}^{2n}(z)
&\leq
\int_K
\int_{\mathbb{C}^n}
\rho_{\varepsilon_j}(\zeta)\,|\varphi(z-\zeta)|\,d\mathcal{L}^{2n}(\zeta)
\,d\mathcal{L}^{2n}(z) \\
&=
\int_{\mathbb{C}^n}
\rho_{\varepsilon_j}(\zeta)
\left(
\int_{K-\zeta} |\varphi(w)|\,d\mathcal{L}^{2n}(w)
\right)
d\mathcal{L}^{2n}(\zeta) \\
&\leq
\left(\int_{\mathbb{C}^n}\rho_{\varepsilon_j}(\zeta)\,d\mathcal{L}^{2n}(\zeta)\right)
\int_{K_\eta}|\varphi(w)|\,d\mathcal{L}^{2n}(w) \\
&=
\int_{K_\eta}|\varphi(w)|\,d\mathcal{L}^{2n}(w)
<\infty.
\end{align*}
The sequence $\varphi_j-\varphi$ is non-negative on $K$ and decreases pointwise to $0$. The preceding estimate with $j=j_0$ gives $\varphi_{j_0}\in L^1(K)$, and $K\subset K_\eta$ gives $\varphi\in L^1(K)$. Since both functions are finite $\mathcal{L}^{2n}$-almost everywhere on $K$, their difference $\varphi_{j_0}-\varphi$ is an element of $L^1(K)$. Moreover,
\begin{align*}
0 \leq \varphi_j-\varphi \leq \varphi_{j_0}-\varphi
\end{align*}
for every $j\geq j_0$. The [Dominated Convergence Theorem](/theorems/4) applies with dominating function $\varphi_{j_0}-\varphi\in L^1(K)$ and pointwise limit $0$, so
\begin{align*}
\lim_{j\to\infty}
\int_K |\varphi_j-\varphi|\,d\mathcal{L}^{2n}(z)
=
\lim_{j\to\infty}
\int_K (\varphi_j-\varphi)\,d\mathcal{L}^{2n}(z)
=0.
\end{align*}
Because $K\subset\Omega_0$ was arbitrary, $\varphi_j\to\varphi$ in $L^1_{\mathrm{loc}}(\Omega_0)$. This completes the proof.
[/step]