Why is there a 1/2 in the expression for electrostatic energy U? The expression for electrostatic energy is $$U= \dfrac{1}{2} \times \int \rho\,\phi\, dV$$ where $\rho$ is the charge density and $\phi$ is the potential at that point in ($dV$) 
Let me explain what I understand about the process of calculating the $U$ intuitively. We take a small chunk of volume, multiply it by potential, to arrive at the workdone in bringing that chunk of charge to that place. You do that to all chunks in the volume and add all of it. So shouldn't $U$ just be $$U= \int \rho\,\phi \,dV$$ Why do we halve this energy ?
 A: You multiply it by $1/2$ because when performing that sum you have double-counted the potential. Each movement you have described is assuming the rest of the charge density is already there. We actually want to think about assembling the configuration from nothing. Therefore, if you assume the rest of the charge density is already present each time, you have double-counted the work you needed to do to make the charge configuration.

A simple derivation of this is found in Griffith's Introduction to Electrodynamics book.
First, you find the work done to assemble your charge configuration:
$$W=\frac{1}{4\pi\epsilon_0}\sum_{i=1}^n\sum_{j=1\\j>i}^n\frac{q_iq_j}{r_{ij}}$$
where the $j>i$ condition represents building the configuration up from nothing so that we only include each charge pair once. However, we can intentionally double count and divide by $2$:
$$W=\frac12\cdot\frac{1}{4\pi\epsilon_0}\sum_{i=1}^n\sum_{j=1\\j\neq i}^n\frac{q_iq_j}{r_{ij}}$$
The reason we do this is so that now the order of charge assembly does not matter. Therefore:
$$W=\frac12\sum_{i=1}^n q_i\left(\sum_{j=1\\j>i}^n\frac{1}{4\pi\epsilon_0}\frac{q_j}{r_{ij}}\right)$$
Since the term in parentheses is the potential due to all of the charges at the location of charge $q_i$, we have.
$$W=\frac12\sum_{i=1}^n q_iV(\mathbf r_i)$$
This is just the discrete version of your integral, and as you can see the $1/2$ came from our double-counting.
A: You can think about it by considering discrete (point-like) charges. If you have 2 charges the energy is: 
$$ U= \frac{q_1q_2}{4\pi\epsilon_0}$$
If you used your formula considering $V_i(r)=q_i/4\pi\epsilon_0$ then your total energy would be: 
$$ U= q_1\frac{q_2}{4\pi\epsilon_0} +  q_2\frac{q_1}{4\pi\epsilon_0} = 2 \frac{q_1q_2}{4\pi\epsilon_0} $$
This is two times the actual potential energy. You can see that you only need to consider the interaction one time, but using the integral:
$$U = \int \rho \phi dV$$
you count the potential energy contributions twice: one when a charge (say $q_1$) interacts with a potential of another charge (say $q_2$) and a second one when the potential of the first charge ($q_1$) interacts with the second charge ($q_2$). But the total energy between two charges is $q_1q_2/4\pi\epsilon_0$, so we only need to count it once. This reasoning can be extended to a continuous charge distribution $\rho$
However, counting it twice and then diving by 2 makes the calculation easier, otherwise, a restriction must be added to the integral such that when an interaction is counted, it is then excluded. This restriction is written in the discrete case as: 
$$ U = \sum_{i=1}^N \sum_{j>i}^N  U= \frac{q_iq_j}{4\pi\epsilon_0}$$
where $j>i$ prevents the contributions of being counted twice. 
A: Here are examples with 2 and then 3 charges.
$2$ charges:
$$
U = \frac{q_1 q_2}{4 \pi \epsilon_0 r_{12}^2} \\
= q_1 \phi_2({\bf r}_1) \\
= (\rho({\bf r}_1) dV )\phi_2({\bf r}_1)
$$
where the notation $\phi_2({\bf r}_1)$ means the potential due to 2 evaluated at location ${\bf r}_1$. But for these same two charges we can equally well write
$$
U = q_2 \phi_1({\bf r}_2)
= (\rho({\bf r}_2) dV )\phi_1({\bf r}_2) .
$$
Therefore we can if like write it as half the sum of these two equally valid expressions:
$$
U = \frac{1}{2} \left( \rho({\bf r}_1) \phi_2({\bf r}_1)
+ \rho({\bf r}_2) \phi_1({\bf r}_2) \right) dV\\
= \frac{1}{2} \sum_i \rho({\bf r}_i) \phi({\bf r}_i) dV
$$
where $\phi$ refers to the total potential at any given point produced by all charges except a point charge located exactly there. 
Next let's proceed $3$ charges. There are three pairs: (1,2), (1,3), (2,3). We adopt
the shorthand $\phi_{12} \equiv \phi_1({\bf r}_2)$. Each pairwise interaction can be written as above, as half the sum of two equally valid ways of writing it:
$$
U = \frac{1}{2} (q_1 \phi_{21} + q_2 \phi_{12}) 
+\frac{1}{2} (q_1 \phi_{31} + q_3 \phi_{13}) 
+\frac{1}{2} (q_2 \phi_{32} + q_3 \phi_{23}) 
\\
=
\frac{1}{2}\left( 
  q_1 [\phi_{31} + \phi_{21}]
+ q_2 [\phi_{32} + \phi_{12}]
+ q_3 [\phi_{13} + \phi_{23}]
\right) \\
= \frac{1}{2} \left( 
\rho({\bf r}_1)[\phi_2({\bf r}_1) + \phi_3({\bf r}_1)]
+\rho({\bf r}_2)[\phi_1({\bf r}_2) + \phi_3({\bf r}_2)]
+\rho({\bf r}_3)[\phi_1({\bf r}_3) + \phi_2({\bf r}_3)]
\right) \\
= \frac{1}{2} \sum_i \rho({\bf r}_i) \phi({\bf r}_i) dV
$$
Once you get the idea, you realise that the sum without the $1/2$ exactly double-counts all the pairwise interactions, and hence the overall result follows for any number of charges, and for the limit of a continuous charge distribution where the sum becomes an integral.
