Do you know that the energy is associated with a charge distribution . This  associated energy is also known as Electrostatic Energy. Let us calculate this energy associated with charge distribution:

When a test charge is moved from a point of lower potential to a point of higher potential, the potential energy of the system increases. But, there is a decrease in potential energy when a test charge is moved from higher to lower potential.

The potential energy imparted to the system by external work on a charge when moved from lower to higher potential is

W=QV

Where W is the work done, Q is the test charge. To calculate the energy of a system say consisting of three points charges Q_1­,Q₂,Q₃ ,we first bring the charge Q₁ from infinity to position 1,no work is done here as the field does not exist in the system.

Thus                                              W₁=0

Similarly for q₂, work done ,W₂=Q₂V₂₁.where v_2.1 is the potential at the location of Q₂due to Q₁.Also for Q₃we have W₃=Q₃V_3.1+Q₃.V_3.2 where V_3.1 and V_3.2 are potentials at position 3 due to Q_1­ and Q₂respectively.

Total work done                                          W=W₁+W₂+W₃ (1)

=0+Q₂V₂V_2.1+Q₃V_3.1+Q₃V_3.2

If charges are brought in the reverse order then,

W=W₃+W₂+W₁=0+Q₂V_2.3+Q₁V_1.3+Q₁V_1.2 (2)

Adding equations (2) and (3),we get

2W=Q₁V_1.3+Q₁v_1.2+Q₂V₂₁+Q₂V₂₁+Q₃V_3.1+Q₃V_3.2 (3)

V_1.3+V_1.2 is the sum of potential at position Q₁due to charges Q₃ amd Q₂ respectively i.e total potential at Q₁ is V₁ .Similarly total potential at Q₂ ans Q₃ is V₂ and V₃ respectively. Thus equation (3) can be written as

2W=Q₁V₁+Q₂v₂+Q₃V₃

Or                                                          W=1/2 (Q₁V₁+Q₂V₂+Q₃V₃)

For large number of charges ,we have

W=1/2 ∑ⁿ_i=1Q_iV_i

But for a continuous charge distribution, work done in a volume V having charge density ρ_v enclosed within a sphere of radius R is

W=1/2∫ρ_v VdV                                                  (Q=ρ_vdV)

W=1/2 ∫(del.D)VdV                          [ from Gauss’s law del.D=ρ_v]         (4)

As per vector identity ,we know that

Del(del D)=(del D)V+D.(delV)

Or                                                               (Del.D)V=(del.VD)-(D.del V)

Put this equation in equation (4)

W=1/2 ∫(del.VD)dV-1/2∫(D.delV)dV                                                             (5)

Now for large R→∞,using divergence theorem

=1/2 ∫(Del.VD)dV =1/2 ∫VD.dS

At very large radius r,D=1/R²,V∞1/R so VD∞1/R³ and surface ∞R²

Thus               ½ ∫VD.dS=0                                             ( 1/2∫(1/R³)R²=0)

From equation (s),we have

W=1/2 ∫(D.delV)dV

Since                                                     E=-del V

W=1/2 ∫(D.E)dV

Also                                                       D=ε₀E

Thus ,Energy                                     W=1/2 ∫ε₀E²dV

In the above expressions,it appears that each tiny volume dV be assigned the energy content WdV where

Energy density                            =E/dV

W_E=1/2 ε₀E²

Here                                                    ε=ε₀ε_r

Energy Density,

W_E=1/2 ε₀ε_rE²