Here’s the sketch of an answer coming from another perspective.
First, here’s a fact that is not discussed in Jackson’s book. When performing a Cartesian multipole expansion, the charge density $\rho$ is changed to a “generalized point source,” i.e., [1]
$$\rho\approx \tilde{\rho}=q_0\delta-\mathbf{p}\cdot\nabla\delta+\frac{1}{2}\mathbf{Q}\nabla\nabla\delta+\dots$$
where $q_0$ is the net charge, $\mathbf{p}=\iiint \mathbf{r}\rho(\mathbf{r})\text{d}^3\mathbf{r}$ is the dipole moment (a vector), $\mathbf{Q}=\iiint \mathbf{r}\mathbf{r}^\top\rho(\mathbf{r})\text{d}^3\mathbf{r}$ is the (symmetric but not trace-free) quadrupole moment tensor, and $\delta$ is the 3D Dirac $\delta$ distribution.
In other words, each Cartesian multipole moment corresponds to a particular derivative of the Dirac $\delta$ distribution, say
$$\frac{\partial^{\alpha_x+\alpha_y+\alpha_z}\delta}{\partial x^{\alpha_x}\partial y^{\alpha_y}\partial z^{\alpha_z}}$$
where $\alpha_x,\alpha_y,\alpha_z$ are integers. For example, the dipole moments correspond to all derivatives of order 1. Let’s go one step further and recognize that each derivative term can be identified with a corresponding monomial $x^{\alpha_x}y^{\alpha_y}z^{\alpha_z}$.
Also, the field $f_{(\alpha_x,\alpha_y,\alpha_z)}$ “radiated” by a Cartesian moment $\alpha_x,\alpha_y,\alpha_z$ can easily be obtained from Green's function $g$. By definition, the latter satisfies
$$\nabla^2g=\delta$$
thus
$$\nabla^2\frac{\partial^{\alpha_x+\alpha_y+\alpha_z}g}{\partial x^{\alpha_x}\partial y^{\alpha_y}\partial z^{\alpha_z}}=\frac{\partial^{\alpha_x+\alpha_y+\alpha_z}\delta}{\partial x^{\alpha_x}\partial y^{\alpha_y}\partial z^{\alpha_z}}$$
by exchanging derivatives. Thus, $f_{(\alpha_x,\alpha_y,\alpha_z)}=\frac{\partial^{\alpha_x+\alpha_y+\alpha_z}g}{\partial x^{\alpha_x}\partial y^{\alpha_y}\partial z^{\alpha_z}}$.
Now, this is not trivial (e.g., [2-3]), but it turns out that all spherical moments also correspond to a particular set of derivatives of the Dirac $\delta$ distribution. Namely, given the field $f_{(l,m)}$ radiated by a spherical moment $(l,m)$, we have that
$$\nabla^2f_{(l,m)}\propto P_{(l,m)}(\nabla)\delta\tag{1}$$
where $P_{(l,m)}$ is a solid harmonic of order $l$ (i.e., a harmonic and homogeneous polynomial of three variables) and $P_{(l,m)}(\nabla)\delta$ means that we evaluate the polynomial in $(\partial/\partial_x,\partial/\partial_y,\partial/\partial_z)$ and apply the corresponding linear combination of derivatives to the Dirac $\delta$ distribution. Now, again, we'll identify a spherical moment $(l,m)$ with its corresponding polynomial $P_{(l,m)}$. The question can be rephrased as: “How can we write an arbitrary monomial as a sum of harmonic and homogeneous polynomials?”
This is a standard question in harmonic function theory [4]. It turns out that there is always a unique decomposition ($\alpha_x+\alpha_y+\alpha_z=l$)
$$x^{\alpha_x}y^{\alpha_y}z^{\alpha_z}=\sum_{m}\beta_{(l,m)}P_{(l,m)}(\mathbf{r})+r^2\sum_{m}\beta_{(l-2,m)}P_{(l-2,m)}(\mathbf{r})+\cdots+r^{2k}\sum_{m}\beta_{(l-2k,m)}P_{(l-2k,m)}(\mathbf{r})$$
where $k$ is the largest integer such that $l-2k\geq 0$ and the $\beta_{(\cdot,\cdot)}$ are complex numbers. Going back to the multipole expansion, this means that*
$$f_{(\alpha_x,\alpha_y,\alpha_z)}=\sum_{m}\tilde{\beta}_{(l,m)}f_{(l,m)}+\nabla^2\sum_{m}\tilde{\beta}_{(l-2,m)}f_{(l-2,m)}+\cdots+(\nabla^2)^k\sum_{m}\tilde{\beta}_{(l-2k,m)}f_{(l-2k,m)}$$
Outside the origin, by Equation (1), all terms to which the Laplacian is applied vanish, and we get
$$f_{(\alpha_x,\alpha_y,\alpha_z)}=\sum_{m}\tilde{\beta}_{(l,m)}f_{(l,m)}$$
Weird! As Jackson highlights, there is no contradiction (i.e., the Cartesian expansion can be written in terms of the spherical ones -- there are no redundant terms). Still, we did not need to use spherical harmonics of lower order...
Take, for example, $l=2$ (quadrupoles). The Cartesian moments are defined by a symmetric matrix (6 degrees of freedom). This matrix can be written as the sum of a symmetric trace-free tensor (5 degrees of freedom) plus the trace times the identity tensor (1 degree of freedom). Both these matrices are “irreducible,” i.e., their components do not mingle when undergoing rotations. The reasoning above shows that the trace might be non-zero, but it does not contribute to the field (not true for the Helmholtz equation!).
Many open questions remain; I hope this helps or at least triggers the interest of someone with a cleaner answer ;)
[1] C. A. Kocher, “Point‐multipole expansions for charge and current distributions,” American Journal of Physics, vol. 46, no. 5, pp. 578–579, May 1978.
[2] Z. Idziaszek and T. Calarco, “Pseudopotential Method for Higher Partial Wave Scattering,” Phys. Rev. Lett., vol. 96, no. 1, p. 013201, Jan. 2006, doi: 10.1103/PhysRevLett.96.013201.
[3] F. Stampfer and P. Wagner, “A mathematically rigorous formulation of the pseudopotential method,” Journal of Mathematical Analysis and Applications, vol. 342, no. 1, pp. 202–212, Jun. 2008.
[4] S. Axler, P. Bourdon, and R. Wade, Harmonic Function Theory. Springer Science & Business Media, 2013.
*We have to replace $\beta$ with $\tilde{\beta}$ because of the proportionality relation in Equation (1).
