The idea is that because the electrostatic field has zero curl, it can be written as the gradient of a potential $\phi$ via
$$\mathbf{E} = - \nabla \phi.$$
Now, using this starting point, the energy of the electrostatic field of a system of charged conductors can then be expressed in terms of the potential $\phi$ as
$$U = \frac{1}{8 \pi} \int \mathbf{E}^2 dV = \frac{1}{2} \sum_i e_i \phi_i$$
where $\phi_i$ is the potential on the $i$'th conductor, and the $e_i$ are defined explicitly by just working the integral out directly
\begin{align}
U &= \frac{1}{8 \pi} \int \mathbf{E} \cdot \mathbf{E} dV \\
&= - \frac{1}{8 \pi} \int \mathbf{E} \cdot \nabla \phi dV \\
&= - \frac{1}{8 \pi} \int \nabla \cdot (\mathbf{E} \phi ) dV + \frac{1}{8 \pi} \int (\nabla \cdot \mathbf{E}) \phi dV \\
&= - \frac{1}{8 \pi} \int \phi \mathbf{E} \cdot \mathbf{n} dS + 0 \\
&= + \frac{1}{8 \pi} \int \phi E_n dS \\
&= + \sum_i \frac{1}{8 \pi} \int \phi_i E_n dS \\
&= + \sum_i \phi_i \frac{1}{8 \pi} \int E_n dS \\
&= + \frac{1}{2} \sum_i \phi_i e_i
\end{align}
where these steps should be obvious (e.g. write $\mathbf{E} = - \nabla \phi$ or use Maxwell's equation $\nabla \cdot \mathbf{E} = 0$ or treat $\phi_i$ constant on the surface conductor $i$ so it can be pulled out of the integral, or note $\mathbf{n}$ is in the direction opposite to $\mathbf{E}$ so the sign changes when we form $E_n$, with the definition of $e_i$ as an integral clear) see the reference below for more details if necessary). Note how this relies on the zero curl property stated in your post.
Next, since the charges and potentials on the conductors are not independent, they must be related. Since the electrostatic Maxwell equations in a vacuum are linear, one assumes the relation between the $e_i$'s and $\phi_j$'s in the form
$$e_i = C_{ij} \phi_j,$$
which is the same idea as that which leads to equation (1) of your link, but now we can express the energy of the electrostatic field in terms of $C_{ij}$ as
$$U = \frac{1}{2} \sum_i \phi_i e_i = \frac{1}{2} \sum_i \phi_i \sum_j C_{ij} \phi_j = \frac{1}{2} \sum_{i,j} \phi_i C_{ij} \phi_j.$$
This doesn't yet show that $C_{ij}$ should be symmetric, however (see the reference below for a proof) one can also use Maxwell's equations to separately derive
$$
\delta U = \sum_i \phi_i \delta e_i$$
which implies
$$\frac{\partial U}{\partial e_i} = \phi_i$$
and derive
$$
\delta U = \sum_i e_i \delta \phi_i$$
which implies
$$\frac{\partial U}{\partial \phi_i} = e_i$$
(note these results are not immediately obvious by blindly differentiating any of the above forms of $U$) and using these latter results we can directly show
$$
\frac{\partial^2 U}{\partial \phi_i \phi_j} = \frac{\partial }{\partial \phi_i} [ \frac{\partial U}{\partial \phi_j} ] = \frac{\partial }{\partial \phi_i} e_j = \frac{\partial }{\partial \phi_i} \sum_k C_{jk} \phi_k = C_{ji}$$
then from
$$
\frac{\partial^2 U}{\partial \phi_i \phi_j} = \frac{\partial^2 U}{\partial \phi_j \phi_i} = C_{ij} $$
you have
$$C_{ij} = C_{ji}.$$
Reference: Landau, Electrodynamics of Continuous Media, section 2.