Why is the Sun almost perfectly spherical?

The symmetry of the Sun has got very little to do with any symmetry in its formation.

The Sun has had plenty of time to reach an equilibrium between its self gravity and its internal pressure gradient. Any departure from symmetry would imply a difference in pressure in regions at a similar radius but different polar or azimuthal angles. The resultant pressure gradient would trigger fluid flows that would erase the asymmetry.

Possible sources of asymmetry in stars could include rapid rotation or the presence of a binary companion, both of which break the symmetry of the effective gravitational potential, even if the star were spherically symmetric. The Sun has neither of these (the centrifugal acceleration at the equator is only about 20 millionths of the surface gravity, and Jupiter is too small and far away to have an effect) and simply relaxes to an almost spherically symmetric configuration.

The relationship between oblateness/ellipticity and rotation rate is treated in some detail here for a uniform density, self-gravitating spheroid and the following analytic approximation is obtained for the ratio of equatorial to polar radius $$ \frac{r_e}{r_p} = \frac{1 + \epsilon/3}{1-2\epsilon/3}, $$ where $\epsilon$, the ellipticity is related to rotation and mass as $$\epsilon = \frac{5}{4}\frac{\Omega^2 a^3}{GM}$$ and $a$ is the mean radius, $\Omega$ the angular velocity.

Putting in numbers for the Sun (using the equatorial rotation period), I get $\epsilon=2.8\times10^{-5}$ and hence $r_e/r_p =1.000028$ or $r_e-r_p = \epsilon a = 19.5$ km. Thus this simple calculation gives the observed value to a small factor, but is obviously only an approximation because (a) the Sun does not have a uniform density and (b) rotates differentially with latitude in its outer envelope.

A final thought. The oblateness of a single star like the Sun depends on its rotation. You might ask, how typical is the (small) rotation rate of the Sun that leads to a very small oblateness? More rapidly rotating sun-like (and especially more massive) stars do exist; very young stars can rotate up to about 100 times faster than the Sun, leading to significant oblateness. However, Sun-like stars spin-down through a magnetised wind as they get older. The spin-down rate depends on the rotation rate and this means that single (or at least stars that are not in close, tidally locked binary systems) stars converge to a close-to-unique rotation-age relationship at ages beyond a billion years. Thus we expect (it remains to be proven, since stellar ages are hard to estimate) that all Sun-like stars with a similar age to the Sun should have similar rotation rates and similarly small oblateness.

I) In this answer we will only discuss the equilibrium shape. Recall that when we discussed the shape of Earth in this Phys.SE post, the gravitational quadrupole moment was important. Unlike the Earth, from a surface perspective, it is a very good approximation to assume that all the mass of the Sun sits in the center, cf. below graph.

Moreover, the Newton's shell theorem helps out here. We conclude that it is enough to consider the gravitational monopole field

$$\tag{1} g(r)~=~\frac{GM}{r^2}$$

of the Sun. From Wikipedia, we get that

$$\tag{2} G~=~ 6.674\cdot 10^{−11} {\rm Nm}^2/{\rm kg}^2\quad\text{and}\quad M~=~(1.98855 \pm 0.00025)\cdot 10^{30} {\rm kg}. $$ The equatorial radius and period are $$\tag{3} r_e~=~(696342\pm 65)~{\rm km} \quad\text{and}\quad T_e~=~25.05 ~{\rm days}, $$

respectively. The equatorial speed is

$$\tag{4} v_e~=~\omega_e r_e~=~\frac{2\pi r_e}{T_e}~\approx~2.02 ~{\rm km/s}.$$

The equatorial surface gravity is then

$$\tag{5} g_e~=~\frac{GM}{r_e^2}~\approx~274~{\rm m/s^2}. $$

Repeating Mark Eichenlaub's monopole argument for the Sun, the height difference between the equatorial and polar radius becomes

$$\tag{6} h~:=~r_e-r_p~=~\frac{v_e^2}{2g_e}~\approx~7.5 ~{\rm km}, $$

leading to a flattening

$$\tag{7} f~=~\frac{h}{r_e}~\approx~ 11 \cdot 10^{-6} .$$

This estimate overshoots by 20% the actual observed flattening, which is only $9 \cdot 10^{-6}$.

II) In the remainder of this answer, we would like to argue that the 20% difference in eq. (7) is mainly due to the fact that the Sun does not spin as a rigid body, which we implicitly assumed in Section I. The polar period

$$\tag{8} T_p~=~34.4 ~{\rm days} $$

is slower than the equatorial period (3). To proceed, let us for simplicity assume that the square $T^2$ of the period $T$ depends on the polar angle $\theta$ in the following way$^1$

$$\tag{9} T^2~=~T_p^2 + s (T_e^2-T_p^2) , \qquad s~\equiv~\sin^2\theta,\qquad \omega~\equiv~ \frac{2\pi}{T}. $$

Analogously, define for later convenience the quantity

$$\tag{10} A~:=~\frac{GM}{\omega^2}~=~ A_p + s A^{\prime}, \qquad A^{\prime}~:=~A_e-A_p~<~0, $$

which is proportional to $T^2$. The centrifugal acceleration is $$\tag{11} a_{\rm cf}~=~\omega^2 r\sin\theta.$$

Using arguments similar to my Phys.SE answer here, the total force should be perpendicular to the surface

$$\tag{12} \left(g -a_{\rm cf} \sin \theta \right)\mathrm{d}r -a_{\rm cf} \cos \theta ~r\mathrm{d}\theta~=~0.$$

The differential (12) is inexact. After multiplying with an integrating factor, we have

$$ \tag{13} \mathrm{d}U~=~\lambda(u) \left[\left(\frac{A}{r^2}-sr\right)\mathrm{d}r -\frac{r^2}{2}\mathrm{d}s\right], $$

where

$$ \tag{14} \lambda(u)~:=~\exp\left(\frac{2}{3}A^{\prime} u^3 \right) , \qquad u~\equiv~\frac{1}{r}. $$

The potential becomes

$$ \tag{15} U~=~ -A_p \int_0^u \! du^{\prime} ~ \lambda(u^{\prime}) - s\frac{\lambda(u)}{2u^2} . $$

The difference between equatorial and polar potential should be zero:

$$\tag{16} 0~=~U_e-U_p~=~ A_p \int_{u_e}^{u_p} \! du~\lambda(u) -\frac{\lambda(u_e)}{2u_e^2} , $$

or equivalently,

$$ \frac{1}{2 A_p u_e^2} ~\stackrel{(16)}{=}~ \int_{u_e}^{u_p} \! du~ \exp\left(\frac{2}{3}A^{\prime} (u^3-u_e^3) \right)~$$ $$ \tag{17}\approx~ \int_{u_e}^{u_p} \! du~ e^{2 A^{\prime} (u-u_e) u_e^2} ~=~\frac{e^{2 A^{\prime} (u_p-u_e) u_e^2}-1}{2A^{\prime} u_e^2}.$$

The height difference becomes

$$\tag{18} h~:=~r_e-r_p~\approx~\frac{u_p-u_e}{u_e^2}~\stackrel{(17)}{\approx}~\frac{r_e^4}{2A^{\prime}} \ln \left(1 + \frac{A^{\prime}}{A_p}\right)~\approx~5.3~{\rm km},$$

leading to a flattening

$$\tag{19} f~=~\frac{h}{r_e}~\approx~ 8 \cdot 10^{-6} ,$$

which is 10% below the observed flattening. Anyway, the above simple model demonstrates that it is important to take into account the non-rigid differential rotation of the Sun.

--

$^1$ Besides fulfilling the correct boundary conditions, the ansatz (9) is admittedly chosen to make the integrating factor (14) simple (rather than being based on observations or astrophysical models).