I have recently read the proof of the following theorem, i.e.,
Theorem Let $I$ be an open subset of $\mathbb{R}$. Let $f \in L^1_{\mathrm{loc}} (I)$ such that $\int_I f \varphi' =0$ for all $\varphi \in C^\infty_c (I)$. Then $f$ is constant a.e. in $I$.
Proof We fix a function $\psi \in C^\infty_c (I)$ such that $\int_I \psi = 1$. Let $w \in C^\infty_c (I)$ and $h := w - ( \int_I w) \psi$. Clearly, $h \in C^\infty_c (I)$. We define $\varphi \in C^\infty (I)$ by $\varphi (x) := \int_{-\infty}^x h(s) \, \mathrm d s$. By FTC, $\varphi' = h$. Because $\int_I h =0$, we get $\varphi \in C^\infty_c (I)$. Then $$ \begin{align*} 0 &= \int_I f \varphi' \\ &= \int_I f h \\ &= \int_I f \left ( w - \left ( \int_I w \right ) \psi \right) \\ &= \int_I \left ( f - \int_I f \psi \right) w. \end{align*} $$ Because $w$ is arbitrary in $C^\infty_c (I)$, we get $f = \int_I f \psi$ a.e. on $I$.
I would like to ask if above theorem holds in higher dimension, i.e.,
Let $\Omega$ be an open subset of $\mathbb{R}^d$. Let $f \in L^1_{\mathrm{loc}} (\Omega)$ such that $\int_\Omega f \nabla \varphi =0$ for all $\varphi \in C^\infty_c (\Omega)$. Then $f$ is constant a.e. in $\Omega$.
Thank you so much for your elaboration!
I would try it with smoothing by convolutions. Here is a sketch:
Let $\rho\in C_{c}^{\infty}\left( B\right) $ be a mollifier with $\rho\geq0$, $\rho \left( -x\right) =\rho \left( x\right) $ and $\int_{B}\rho\,\mathrm{d}x=1$, where $B$ is the unit ball around the origin. Set $\rho _{k}\left( x\right) =k^{n}\rho\left( kx\right) $. Then the convolution $\rho_{k}\ast f$ is smooth and $$ \rho_{k}\ast f\rightarrow f\quad\text{in }L_{\mathrm{loc}}^{1}\left( \Omega\right) , $$ in particular, $$ \left( \rho_{k}\ast f\right) \left( x\right) \rightarrow f\left( x\right) \quad\text{for a.e. }x\in\Omega. $$
So $$ 0=\int_{\Omega}f\nabla\left( \rho_{k}\ast\varphi\right) \,\mathrm{d}% x=\int_{\Omega}f\left( \rho_{k}\ast\nabla\varphi\right) \,\mathrm{d}% x=\int_{\Omega}\left( \rho_{k}\ast f\right) \nabla\varphi\,\mathrm{d}% x=-\int_{\Omega}\nabla\left( \rho_{k}\ast f\right) \varphi\,\mathrm{d}x $$ for all $\varphi\in C_{c}^{\infty}\left( \Omega\right) $, implying that $$ \nabla\left( \rho_{k}\ast f\right) =0\quad\Rightarrow\quad\rho_{k}\ast f=\mathrm{const}\quad\text{if }\Omega\text{ is connected.}% $$ By the above, also $$ f\left( x\right) =\mathrm{const}\quad\text{for a.e. }x\in\Omega. $$
Another proof uses the Sobolev space $W^{1,1}$ and the Poincaré inequality $$ \int_{B_{R}}\left\vert f-f_{B_{R}}\right\vert \,\mathrm{d}x\leq c\left( d\right) R\int_{B_{R}}\left\vert \nabla f\right\vert \,\mathrm{d}x $$ for $f\in W^{1,1}\left( B_{R}\right) $, where $B_{R}$ is a ball of radius $R$. Under your hypothesis, $$ f\in W_{\mathrm{loc}}^{1,1}\left( \Omega\right) \quad\text{with}\quad\nabla f=0, $$ so that $f=\mathrm{const}=f_{B_{R}}$ on any ball $B_{R}\Subset\Omega$.
As a reference on smoothing by convolutions and on Sobolev spaces, I recommend
Gilbarg & Trudinger, Elliptic Partial Differential Equations of 2nd Order, Chapter 7.