Why are the results for electric and magnetic field energy densities always valid? The result for the electric field density $\frac{1}{2}εΕ^2$ is derived by writing the sum of the electric fields produced by a small element of a conductor and the rest of the conductor to be $0$ at a point just inside the conductor next to the element, and, the result for magnetic field energy density $\frac{1}{2}μB^2$ is derived by calculating $\frac{1}{2}LI^2$ for an ideal solenoid and then dividing it by its volume. Hence these results should be valid only for calculating the energy at the surface of a conductor and the magnetic field energy of an ideal solenoid respectively. Yet the net electric and magnetic field energies of a given system are often calculated by integrating the electric and magnetic field energy densities all over space. How is doing so correct?
 A: Short answer:
On the contrary. The expressions for the electric energy stored at the surface of the conductor and the magnetic energy stored inside an ideal solenoid are both derived from the expression for the energy density $\mathcal U = \frac{1}{2}\left(\epsilon_0\vec E \cdot \vec E + \frac{1}{\mu_0}\vec B \cdot \vec B\right)$ due to any arbitrary electromagnetic field,  not the other way around. The expression for $\mathcal U$ in terms of electromagnetic fields is a direct consequence of Maxwell's equations,
$$\nabla \cdot \vec E = \frac{\rho}{\epsilon_0} \tag{1}$$
$$\nabla \cdot \vec B = 0 \tag{2}$$
$$\nabla \times \vec E = -\frac{\partial \vec B}{\partial t} \tag{3}$$
$$\nabla \times \vec B = \mu_o \vec J + \mu_o\epsilon_0 \frac{\partial \vec E}{\partial t} \tag{4}$$
The main implications present in any derivation of the electromagnetic field energy density are that energy must somehow be expended in order to establish a nonzero electromagnetic field in a region of space and that the energy expended ends up stored within the field itself.
Long answer:
I've seen a few different ways to derive the electromagnetic field energy density from fundamental principles of electromagnetism, but the derivation I am most familiar and comfortable with uses a field-theoretic approach adapted from the content and practice problems of Chapter 2 of An Introduction to Quantum Field Theory by Peskin and Schroeder:
Begin with the electromagnetic Lagrangian density of free space (i.e. no charges or currents are present),
$$\mathcal L\left(A_\mu , \partial_\nu A_\mu\right) = -\frac{1}{4\mu_0}F_{\mu\nu}F^{\mu\nu} \tag{5}$$
where
$$F_{\mu\nu} \equiv \partial_\mu A_\nu - \partial_\nu A_\mu = 
\left(\begin{matrix}0 & -E_x/c & -E_y/c & -E_z/c \\
E_x/c & 0 & -B_z & B_y \\
E_y/c & B_z & 0 & -B_x \\
E_z/c & -B_y & B_x & 0 \end{matrix}\right) \tag{6}$$
$$F^{\mu\nu} \equiv \partial^\mu A^\nu - \partial^\nu A^\mu = 
\left(\begin{matrix}0 & E_x/c & E_y/c & E_z/c \\
-E_x/c & 0 & -B_z & B_y \\
-E_y/c & B_z & 0 & -B_x \\
-E_z/c & -B_y & B_x & 0 \end{matrix}\right)$$
are respectively the covariant and contravariant form of the electromagnetic field tensor, with $\partial$ the four-gradient and $A$ the electromagnetic four-potential. The Euler-Lagrange equations for this field are the Maxwell equations in free space,
$$\partial_\mu F^{\mu\nu} = 0 \tag{7}$$
$$\partial_{[\mu} F_{\nu\sigma]} = 0 \tag{8}$$
Due to symmetries inherent to the Lagrangian given in $(5)$, we can construct a tensor $T^{\mu\nu}$, defined as
$$T^{\mu\nu} \equiv \frac{\partial \mathcal L}{\partial\left(\partial_\mu A_\sigma\right)}\partial_\nu A_\sigma - \mathcal Lg^{\mu\nu} \tag{9}$$
where
$$g^{\mu\nu} = \left(\begin{matrix}1 & 0 & 0 & 0 \\
0 & -1 & 0 & 0 \\
0 & 0 & -1 & 0 \\
0 & 0 & 0 & -1 \end{matrix}\right) \tag{10}$$
is the Minkowski metric, and then use $(9)$ along with $(5)$ to construct a symmetric stress-energy tensor $\hat T^{\mu\nu}$:
$$\hat T^{\mu\nu} \equiv T^{\mu\nu} + \partial_\sigma \left(F^{\mu\sigma}A^\nu\right) = 
\frac{1}{\mu_0}\left(F^{\mu\kappa}F_{\kappa\lambda}g^{\lambda\nu} + \frac{1}{4} F_{\mu\nu}F^{\mu\nu}\right) \tag{11}$$
Thus, the electromagnetic field Hamiltonian density (i.e. energy density) $\mathcal H$ is
$$\begin{align} \mathcal H & = \hat T^{00} \\ 
& = \frac{1}{\mu_0}\left(F^{0\kappa}F_{\kappa\lambda}g^{\lambda 0} + \frac{1}{4} F_{\mu\nu}F^{\mu\nu}\right) \\
& = \frac{1}{\mu_0}\left(\frac{\vec E \cdot \vec E}{c^2} + \frac{1}{2} \left(\vec B \cdot \vec B - \frac{\vec E \cdot \vec E}{c^2}\right)\right) \\
& = \frac{1}{2}\left(\epsilon_0\vec E \cdot \vec E + \frac{1}{\mu_0}\vec B \cdot \vec B\right) \end{align} \tag{12}$$
This page provides a more in-depth explanation of the derivation above.
