Why is the Force of Gravitational Attraction between two “Extended” bodies proportional to the product of their masses?

The statement

...the force between Any two rigid masses is only proportional to the product of their masses

is not true in general, or at least it is misleading. The shapes of the mass distributions and their relative positions matter when computing the gravitational force.

It is true that once you hold constant the shapes of the mass distributions and their relative positions, then the force will be proportional to the product of the total masses of the bodies.

There are certain situations where treating two extended massive bodies as point sources can be exactly correct (in the context of Newtonian gravity). For a spherically symmetric mass distribution, the gravitational potential outside of it is the same as that arising from a point source of the same mass. This is an application of Gauss' law.

In general, one can build up an increasingly good approximation of the gravitational potential arising from a given mass distribution via a multi-pole expansion. . The leading-order term, which drops off least rapidly with distance (force $\propto r^{-2}$), is that of a monopole like what arises for a point mass or outside a spherically symmetric system. But a general mass distribution will have contributions from higher-order terms (dipole, quadropole, octopole...), all of which drop off increasingly rapidly with distance. As one considers two bodies at increasing separation, reducing them both to their monopole terms becomes increasingly more accurate.

Finally, the fact that the gravitational force of attraction on an extended body due to another body can vary with position is essential when considering phenomena such as tidal forces.

I edited the answer to make it more readable.

It is true in general that the force between two separate and rigid bodies is equal to $ G M_1 M_2 f$ with $f$ depending on the details of the mass distributions. The main reason because this is true is that the mass is an extensive quantity. However, it is also true that the force is not always proportional to the inverse of the square of the distance. This means that in general, changing the mass distribution will change the mutual force, especially in the case where the two bodies are very close. If we also assume that the distance between the two bodies is much larger than their individual sizes, then one can also recover the inverse square law, because in this case one has that $F\approx G M_1 M_2/R^2$. In the following I demonstrate that

1) $F= G M_1 M_2 f$ with $f$ not depending on the masses but only on the details of the distributions

2) if the distance between the bodies is much larger than the individual sizes, then $F\approx G M_1 M_2/R^2$.

3) I show a counterexample where the mutual force is not $\propto 1/R^2$ but still proportional to the masses.

1) $F= G M_1 M_2 f$

For point masses, one has $$ \mathbf{F}= G M_1 M_2 \frac{\mathbf r_1- \mathbf r_2}{|\mathbf r_1- \mathbf r_2|^3} $$ So, assume that the mass distributions of the two bodies are rigid and separated in space. The total force between the two bodies can be written as an integral over the two mass distributions $$ \mathbf F_{12}= \int d r_1^3 d r_2^3 G \rho_1(\mathbf r_1) \rho_2(\mathbf r_2) \frac{\mathbf r_1- \mathbf r_2}{|\mathbf r_1- \mathbf r_2|^3} $$ The mass distributions satisfy $$ \int d r^3 \rho_{1,2}(\mathbf r) =M_{1,2} $$ In order to get rid of the mass dependences, one can factorize the masses in the density distributions, by defining $$ \rho_{1,2}(\mathbf r) = M_{1,2} \psi_{1,2}(\mathbf r) $$ The crucial point here is that the functions $\psi_1(\mathbf r)$ and $\psi_2(\mathbf r)$ do not depend on the masses $M_1$ and $M_2$, but only on the geometrical distribution of their densities. They satisfy the property $$ \int d r^3 \psi_{1,2}(\mathbf r) =1 $$ One can imagine them as ''normalized'' and adimensional density distributions. So, why these functions do not depend on the mass? Imagine a rigid body or a mass distribution in space with mass $M$ and volume V. Imagine that this body is made of, for example, wood. Then imagine one could change every infinitesimal portion of this body with, for example, iron. The final result is a body with the same ''normalized'' mass distribution $\psi$ but with a mass much greater than the initial one $M'\gg M$. Then, the function $\psi$ do not depend on the mass but only on the geometrical distribution. The same argument can be applied for rigid bodies but also for more complicated mass distributions, for example nonuniform or deformable objects, fluids, etc.

Using the functions $\psi$ one has $$ \mathbf F_{12}= G M_1 M_2 \int d r_1^3 d r_2^3 \psi_1(\mathbf r_1) \psi_2(\mathbf r_2) \frac{\mathbf r_1- \mathbf r_2}{|\mathbf r_1- \mathbf r_2|^3}=$$$$= G M_1 M_2 \cdot \mathbf f_{12}(\psi_1,\psi_2) $$ where the function $$ \mathbf f_{12}(\psi_1,\psi_2)= \int d r_1^3 d r_2^3 \psi_1(\mathbf r_1) \psi_2(\mathbf r_2) \frac{\mathbf r_1- \mathbf r_2}{|\mathbf r_1- \mathbf r_2|^3} $$ depends only in the geometrical details of the distributions and not on the masses $M_!$ and $M_2$.

2) $F\approx G M_1 M_2/R^2$ at large distances

One can go further and redefine the spatial coordinates with respect to the center of mass of the two bodies 1 and 2. That means that $\mathbf r_1={\mathbf r}_1'+\mathbf R_{1}$ and $\mathbf r_2={\mathbf r}_2'+\mathbf R_{2}$ where $\mathbf R_{1}$ and $\mathbf R_{2}$ are the positions of the center of mass of the two bodies. If $ \mathbf R_{12} = \mathbf R_{1}-\mathbf R_{2}$ is their distance, then one has

$$ \mathbf f_{12}(\psi_1,\psi_2)= \int d {r}_1^{\prime 3} d {r}_2^{\prime 3} \psi_1({\mathbf r}_1') \psi_2({\mathbf r}_2') \frac{{\mathbf r}_1'- {\mathbf r}_2'+\mathbf R_{12}}{|{\mathbf r}_1'- {\mathbf r}_2' +\mathbf R_{12}|^3} $$ depends on the distance of the two masses and on the geometrical details of their distributions, but not on the value of the masses $M_1$ and $M_2$.

The dependence on the distance of the center of masses does not necessarily goes as the inverse of the square of the distance. However, if the distance between the two bodies is much larger than the sizes of the two bodies, one has that $|\mathbf R_{12}|\gg |\mathbf r_1- \mathbf r_2| $ then one has, as a first approximation $$ \mathbf f_{12}(\mathbf R_{12}) \approx \int d {r_1}^{\prime 3} d {r}_2^{\prime 3} \psi_1({\mathbf r}_1') \psi_2({\mathbf r}_2') \frac{\mathbf R_{12}}{|\mathbf R_{12}|^3} =\frac{\mathbf R_{12}}{|\mathbf R_{12}|^3} $$ The last equality follows from the fact that the functions $\psi$ are normalized to 1, which gives $$ \int d {r_1}^{\prime 3} d {r}_2^{\prime 3} \psi_1({\mathbf r}_1') \psi_2({\mathbf r}_2')= \int d {r}^{\prime 3} \psi_1({\mathbf r}') \times \int d {r}^{\prime 3} \psi_2({\mathbf r}') =1$$ Therefore in this case one has $$ |\mathbf F_{12}|\approx \frac{G M_1 M_2}{|\mathbf R_{12}|^2} $$

2) $F\neq G M_1 M_2/R^2$ but still $F= G M_1 M_2 f$

The obvious violation to the case 2) is where the distance $R$ between the center of masses of the bodies is not large. This condition is obviously violated if the bodies are one inside the other. This is realized, for example, in the famous case where one has a spherical shell containing a sphere. In this case (see wikipedia) all contributions to the mutual gravitational force cancel each other, and one has $$ \mathbf f_{12}(\psi_1,\psi_2)= \int d r_1^3 d r_2^3 \psi_1(\mathbf r_1) \psi_2(\mathbf r_2) \frac{\mathbf r_1- \mathbf r_2}{|\mathbf r_1- \mathbf r_2|^3}=0 $$ This still satisfy $F=G M_1 M_2 f$ with $f=0$.

Now consider the case of a sphere and imagine to dig a tunnel from one pole to the other and consider the gravitational force of an object inside the tunnel. In this case one has that $ \mathbf f_{12}(\psi_1,\psi_2)=1/R$ which gives $$ {F}_{12}= G M_1 M_2 \frac{1}{|\mathbf r_1- \mathbf r_2|} =\frac{G M_1 M_2}{R_{12}} $$ The dependence on the distance is now $\propto1/R$ but the force is still proportional to $M_1 M_2$ and satisfy the equation derived in 1).

Example of mass distribution

To make things more intuitive, consider the easiest example of mass distribution: A rigid sphere of radius $R$. In this case one has $$ \psi(\mathbf r)= \begin{cases} 1 & \text{for} |\mathbf r|<R\\ 0 & \text{for} |\mathbf r|>R \end{cases} $$

How does this sound?

Let's say you have two extended bodies, A and B, each made up of a number of particles. Let's consider the force on a particle in A, call it P, due to body B. Each particle in body B exerts a force on P that is proportional to that B particle's mass. The sum of such forces gives the net force on P due to B. Now let's suppose we double the mass of each B particle. That would result in doubling the force on P. But doubling the mass of each B particle is doubling the total mass of body B. So doubling the total mass of body B doubles the force on P. Hence, the force due to body B is proportional to B's total mass.

To get the total force due to body B on the body A, we have to add up the forces on all of A's particles. Again, if we doubled each particle's mass, the total force on A would double, as would its total mass. Hence that total force is proportional to A's total mass.