Why do electrostatic potentials superimpose? I've been trying to convince myself that the assertion that I've read in basic E&M books (Halliday & Resnick, Purcell), and even Griffiths, that the electrostatic potential at a point in space is equal to the sum of the potential contributions from each of the individual charges.  My hang-up has been the direction in 3D space that a path integral is brought in from infinity, given that there are multiple charges each creating their own lines with the arbitrary point, and each line extends to infinity along different directions.  (I hope my meaning is not lost here). 
Would this reasoning starting from a system of one charge, and then adding additional charges, be sufficient to explain it:
1) For a system of one charge, the electrostatic potential at an arbitrary point is the negative path integral of the field dot ds from infinity, in the direction along the line connecting the point to the charge.
2) Adding a second charge, the electrostatic potential at the same arbitrary point is the negative path integral of the sum of the fields dot ds from infinity (in the same direction as 1) to the point.
Can this sum be broken down as such:
-a) the path integral of the sum of the fields is the sum of the path integrals from the field of the original charge E1 and the path integral of the field of the second charge E2
-b) the path integral of the field of the second charge dot ds (along the original direction), because of the nature of conservative fields, is equal to the sum of 
---i) the path integral of the field of the second charge along the line 
       aligning the second charge to the arbitrary point, given in the 
       usual form: V = kq/r.  This is in a different direction than 
       the path integral for charge 1, but the change in potential is described 
       by the equation here.
---ii) the path integral of field along the connecting arc at infinite 
        distance.  This connecting arc path integral has length 
        proportional to r (going to infinity), but field strength 
        inversely proportional to r squared, so that this component 
        becomes negligible.
So the potential is then the sum of the original path integral of one charge and the path integral of the second charge along a different direction to infinity.  This logic is repeated for additional charges.
I can draw a picture of my thinking if people suggest that in the comments.   
 A: Electrostatic potentials superpose because Poisson’s equation is linear. (And it is linear because Maxwell’s equations are linear in both the fields and their sources, and the potentials and the fields have a linear relation.)
Specifically, if
$$-\nabla^2\varphi_1=4\pi\rho_1$$
and
$$-\nabla^2\varphi_2=4\pi\rho_2$$
then
$$-\nabla^2(\varphi_1+\varphi_2)=4\pi(\rho_1+\rho_2).$$
These equations are in Gaussian units but the same reasoning applies in any other units.
A: Let us start with the electrostatic field $\mathbf{E}$. One approach to this is to assume two things:


*

*The field of a point charge is given by Coulomb's law: $$\mathbf{E}(\mathbf{r})=\dfrac{1}{4\pi\epsilon_0}\dfrac{q}{r^2}\hat{\mathbf{e}}_r$$

*The field of an arbitrary collection of charges is the sum of the individual fields.
These two assumptions are experimental facts that you can use as the starting point for the theory.
In that case, the most general static charge distribution is specified by a charge density $\rho(\mathbf{r})$. Applying the two principles together you get that
$$\mathbf{E}(\mathbf{r})=\dfrac{1}{4\pi\epsilon_0}\int_{V}\rho(\mathbf{r}')\dfrac{1}{|\mathbf{r}-\mathbf{r}'|^3}(\mathbf{r}-\mathbf{r}')d^3\mathbf{r}'.$$
If you now consider this $\mathbf{E}$ and take its curl you will find that $\nabla \times \mathbf{E}=0$ because of the vector identity
$$\nabla_{\mathbf{r}}\times \dfrac{\mathbf{r}-\mathbf{r}'}{|\mathbf{r}-\mathbf{r}'|^3}=0.$$
In that case it is a mathematical fact that there is $\phi$ such that $$\mathbf{E}=-\nabla \phi$$
with the $-$ sign purely conventional since it could be just as well go inside $\phi$.
So up to this point we have seen that: (1) the two assumptions imply that $\mathbf{E}$ must have zero curl and (2) this implies mathematically the existence of $\phi$.
To answer why $\phi$ satisfies superposition, we ought understand how $\phi$ is actually constructed.
To that matter recall the fundamental theorem for line integrals:
$$\int_\mathbf{a}^\mathbf{b} \nabla \psi \cdot d\mathbf{r}=\psi(\mathbf{b})-\psi(\mathbf{a}).$$
Integrate both sides of the defining equation $\mathbf{E}=-\nabla \phi$ along a path $\mathbf{r}(t)$ starting at $\mathbf{a}$ and ending at $\mathbf{x}$. This immediately yields
$$-(\phi(\mathbf{x})-\phi(\mathbf{a}))=\int \mathbf{E}\cdot d\mathbf{r}$$
Now leaving $\mathbf{x}$ arbitrary this gives you
$$\phi(\mathbf{x})=\phi(\mathbf{a})-\int \mathbf{E}\cdot d\mathbf{r}.$$
Setting $\phi(\mathbf{a})=0$ you get $$\phi(\mathbf{x})=-\int \mathbf{E}\cdot d\mathbf{r}.$$
Now notice that (1) $\mathbf{E}$ superimposes under combinations of charges and (2) the integral is linear. In that case, $\phi$ must also satisfy superposition.
Notice we didn't need to get specific about the path. Want do it? No problem, the path is arbitrary becasue the integral on the RHS is path-independent, nonetheless the canonical choice is to use spherical coordinates and pick a radial path ($\theta,\phi$ constant along it) with $\mathbf{a}$ corresponding to $r\to \infty$. So it is one "incoming radial path". 
Again, this is just a convenient choice that exploits the spherical symmetry of $\mathbf{E}$ for point charges.
Finally let me just mention that this is just one of the approaches. You can start with Maxwell's equations for $\mathbf{E}$, namely $\nabla \cdot \mathbf{E} = \rho/\epsilon_0$ and $\nabla\cdot \mathbf{E} =0$ and then notice that they admit solutions of the form $\mathbf{E} = -\nabla \phi$ and hence that $\nabla^2\phi = -\rho/\epsilon_0$. This immediately implies $\phi$ satisfies the superposition principle as well.
A: Alternatively, potentials superpose because electric fields superpose:
$$
V=-\int_\infty^r \sum_{s}\vec E_s\cdot d\vec \ell=\sum_s \left(-\int_\infty^r \vec E_s\cdot d\vec\ell\right)=\sum_s V_s
$$
where 
$$
V_s=-\int_\infty^r \vec E_s\cdot d\vec\ell
$$
is the potential due to source $s$ creating the field $\vec E_s$ of that source.  The sum becomes an integral when the source distribution is continuous.
Note that, because the field is conservative, the integral from $\infty$ to $\vec r$ does not depend on the path taken.  Thus, from a computational perspective, it is always possible to move on a path of arbitrary large radius without changing $V$ so that the path for each individual $\vec E_s\cdot d\vec \ell$ is radial, thus simplifying each integral, irrespective of the initial starting point at $\infty$.
