I know that galaxies are moving away from us, and so can see that it's intuitive that if space was expanding, then the astronomical observations from Earth would be the same as at all other points in the universe.
But it's that intuition often overlooks details. So, what I would like to know, is whether there is proof that smoothly expanding space in 4 dimensions, accommodates every observer point in this way.
Does anyone know?
Yes.
You can make a model where you have coordinates $t, x,y,z$ where for any $x,y,z$ the universe looks the same.
The metric ends up looking e.g. like $$ds^2=dt^2-(a(t))^2(dx^2+dy^2+dz^2)$$ and you can move your $x,y,z$ to have any value and everything looks the same (though things do look different for different values of $t$). You end up with the densities of matter and radiation being solely functions of $t$ (which is just a a coordinate it doesn't say how fast clocks tick, that will depend on the path the clocks takes and depends on the full metric above). And then you also need to solve for $a(t)$ as a function of time.
The function $a(t)$ relates coordinate distances to actual metric distances. Being a certain coordinate distance away could end up becoming much closer or much farther away based on how the scale factor $a(t)$ changes.
Now having this metric does mean everything looks the same at every point in space. But it doesn't mean your space looks the same in every direction. This is because there are still many spacetime's that have that metric.
For instance you can consider the space $$\{(t,a,A,y,z):a^2+A^2=1\}.$$ Then $x$ can go around like an angle in the $a,A$ plane. In this space we have the same metric as $\{(t,x,y,z)\}$ so locally everything looks the same but in one of them when we travel around $2\pi$ in the $x$ direction we end up back where we started and that doesn't happen if we move in the $y$ direction or the $z$ direction.
So if you want the universe to look the same at every place and in all directions you have to be more careful than just having one of the local metrics that works for that. But that does show that it is hard to check.
One way to check is to start with $\{(t,x,y,z)\}$ and $d\tau^2=dt^2-(a(t))^2(dx^2+dy^2+dz^2)$ and then do a transformation $x\mapsto X=x+\Delta x,$ $y\mapsto Y= y+\Delta y$ and $z\mapsto Z=z+\Delta z,$ for fixed $\Delta x,$ $\Delta y,$ and $\Delta z.$ Then $\{(t,x,y,z)\}=\mathbb R^4=\{(t,X,Y,Z)\}$ and $dt^2-(a(t))^2(dX^2+dY^2+dZ^2)=dt^2-(a(t))^2(dx^2+dy^2+dz^2)= d\tau^2$ so the translation changed nothing. The objects and functions are the same, just the names and labels which were arbitrary. This truly shows it is the same at every point.
Then do another transformation, this time rotating $(x,y,z)$ to $(x',y',z')$ by a fixed amount and direction and noting that $\{(t,x,y,z)\}=\mathbb R^4=\{(t,x',y',z')\}$ and $d\tau^2=dt^2-(a(t))^2(dx^2+dy^2+dz^2)=dt^2-(a(t))^2(dx'^2+dy'^2+dz'^2)\}.$ This shows it looks the same in all directions.
To be clear, this cannot be done for the non isotropic case I mentioned before: $\{(t,a,A,y,z):a^2+A^2=1\}$ with $d\tau^2=dt^2-(a(t))^2(da^2+dA^2+dy^2+dz^2).$ that is because the rotation doesn't map that set $\{(t,a,A,y,z):a^2+A^2=1\}$ to $\{(t,x,y,z)\}.$ Showing that they are different was established when going around $2\pi$ in one takes you back and that doesn't happen for the other. I was just trying to make it clear that we really have shown that the space looks the same in all directions when we have $\{(t,x,y,z)\}$ and $d\tau^2=dt^2-(a(t))^2(dx^2+dy^2+dz^2).$
You are approaching the question from the wrong end.
The expansion of the universe is described by a particular solution to Einstein's equation called the FLRW metric. To derive this metric we have to make some assumptions, and the key assumptions are that the universe is isotropic and homogeneous i.e. that it is the same everywhere.
So the universe being the same everywhere is an initial assumption that goes into describing how the universe expands. It isn't something that is derived from the way the universe expands. You can of course show that the FLRW metric implies isotropy and homogeneity (as Timaeus does) but that's because those assumptions were built into it from the beginning.
