I read in my book electrostatic energy density is $ \frac{1}{2} \epsilon_0 E^2$, where $\epsilon_0$ is permittivity of vacuum and $E$ Is electric field in that region. But I want to know whether it is equally valid in dielectric medium too or do we need to replace $\epsilon_0$ with $\epsilon_r×\epsilon_0$, where $\epsilon_r$ is dielectric constant of the medium?
[Physics] the expression of electrostatic energy density in dielectric medium
capacitancedielectricelectrostaticsenergy
Related Solutions
Dielectric constant $K$ is actually the same thing as relative permittivity, and it increases the overall permittivity $\epsilon$. So in general, whenever you see the permittivity of free space $\epsilon_0$ in an equation, if you're dealing with a dielectric, you can multiply it by the dielectric constant and see how the equation changes.
For example, since $C = \epsilon_0 \frac{A}{D}$, multiplying $\epsilon_0$ by $K$ increases capacitance by $K$. Since $\oint \mathbf E \cdot \mathrm d \mathbf A = \frac{Q}{\epsilon_0}$ by Gauss's Law, multiplying $\epsilon_0$ by $K$ decreases electric field magnitude by $K$.
There is no contradiction because there is no increase in electric flux. The dielectric decreases the electric field magnitude, which decreases the electric flux and decreases the voltage across a capacitor as well. $C = \frac{Q}{V}$ and the capacitance increases because the voltage decreases while the charge remains the same.
The Coulomb force in a medium with relative dielectric constant $\epsilon_r$ is given by your first equation. Only from this follows the electric field strength of a spherical symmetric free charge $Q$ in the dielectric with $$E=\frac{Q}{4\pi\epsilon_0\epsilon_r r^2} \tag{1}$$ which, with the electric displacement $D=\epsilon_r \epsilon_0 E$, results in the correct Gauss Law $$ \int_{sphere} \epsilon_r \epsilon_0 E da=Q \tag{2}$$ This is equivalent to the differential form of Gauss's Law, the Maxwell equation in a dielectric $$ div (\epsilon_r \epsilon_0 \vec E)=\rho$$ where $\rho$ is the free charge density.
Note added after a comment by Zhouran He: In Coulomb's Law for the electric force $F$ exerted by a free charge $q_1$ on a second (test) charge $q_2$ in a dielectric with relative permittivity $\epsilon_r$, only the charge $q_1$ as the source of the force field can be considered to be reduced by the polarization charges of the dielectric to the $q_1/\epsilon_r$ so that the vacuum Coulomb law can be used with this net charge. Even though the charge $q_2$ is also surrounded by polarization charges, the force $F$ exerted by the net charge $q_1/\epsilon_r$ works on the free charge $q_2$. One can alternatively consider $q_2/\epsilon_r$ to be the net charge exerting the force $F$ on the free (test) charge $q_1$. The free charge $q_2$ sees a net charge $q_1/\epsilon_r$ exerting a force $F$ on it according to Coulombs vacuum law. The polarization charges induced by itself around it don't exert a force on itself. The same reasoning applies with interchanged roles of the charges. Thus the second form of Coulombs Law for a dielectric is correct.
Best Answer
We have $\vec{\nabla}\phi=\vec{E}$ and $\vec{\nabla}\cdot\vec{D}=\rho$.
The energy change due to the addition of a tiny amount of charge (density) $d \rho$ is: $$d U_e=\displaystyle{\int_V \phi d \rho~dV=\int_V \phi \vec{\nabla}\cdot d\vec{D}~dV=\int_V (\vec{\nabla}\cdot(\phi d\vec{D})-d \vec{D}\cdot\vec{\nabla}\phi)~dV}$$
Then the total energy is: $\displaystyle{\int\ dU_e=\int_V\int_0^D(\vec{\nabla}\cdot(\phi d\vec{D})-d\vec{D}\cdot\vec{\nabla}\phi)~dV=\int_V\int_0^D\vec{\nabla}\cdot(\phi d\vec{D})-\int_V\int_0^Dd\vec{D}\cdot\vec{\nabla}\phi~dV}$.
The first term vanishes when we integrate over the entire space (as it does in a non-dielectric media) and the only thing left is: $\displaystyle{\int_V\int_0^Dd\vec{D}\cdot\vec{E}~dV}$. Now in a linear medium: $\vec{D}=\epsilon_0\epsilon_r\vec{E}$ so $d\vec{D}=\epsilon_0\epsilon_rd\vec{E}$. Back at the energy: $U_e=\displaystyle{\epsilon_0\epsilon_r\int_V\int_0^Ed\vec{E}\cdot\vec{E}~dV=\epsilon_0\epsilon_r\int_V\int_0^Ed\vec{E}\cdot\vec{E}~dV=\dfrac{1}{2}\epsilon_0\epsilon_r\int_V\vec{E}\cdot\vec{E}~dV}$
So indeed, you're right the density is: $\dfrac{1}{2}\epsilon_0\epsilon_rE^2$.
Note here that $\rho$ is the density of free charges.