Are all atoms spherically symmetric? If so, why are atoms with half-filled/filled sub-shells often quoted as 'especially' spherically symmetric?

In general, atoms need not be spherically symmetric.

  • The source you've given is flat-out wrong. The wavefunction it mentions, $\varphi=\frac{1}{\sqrt3}[2p_x+2p_y+2p_z]$, is in no way spherically symmetric. This is easy to check: the wavefunction for the $2p_z$ orbital is $\psi_{2p_z}(\mathbf r)=\frac {1}{\sqrt {32\pi a_0^5}}\:z \:e^{-r/2a_{0}}$ (and similarly for $2p_x$ and $2p_y$), so the wavefunction of the combination is $$\varphi(\mathbf r)=\frac {1}{\sqrt {32\pi a_0^5}}\:\frac{x+y+z}{\sqrt 3} \:e^{-r/2a_{0}},$$ i.e., a $2p$ orbital oriented along the $(\hat{x}+\hat y+\hat z)/\sqrt3$ axis.

    This is an elementary fact and it can be verified at the level of an undergraduate text in quantum mechanics (and it was also obviously wrong in the 1960s). It is extremely alarming to see it published in an otherwise-reputable journal.

  • On the other hand, there are some states of the hydrogen atom in the $2p$ shell which are spherically symmetric, if you allow for mixed states, i.e., a classical probabilistic mixture $\rho$ of hydrogen atoms prepared in the $2p_x$, $2p_y$ and $2p_z$ states with equal probabilities. It is important to emphasize that it is essential that the mixture be incoherent (i.e. classical and probabilistic, as opposed to a quantum superposition) for the state to be spherically symmetric.

  • As a general rule, if all you know is that you have "hydrogen in the $2p$ shell", then you do not have sufficient information to know whether it is in a spherically-symmetric or an anisotropic state. If that's all the information available, the initial presumption is to take a mixed state, but the next step is to look at how the state was prepared:

    • The $2p$ shell can be prepared through isotropic processes, such as by excitation through collisions with a non-directional beam of electrons of the correct kinetic energy. In this case, the atom will be in a spherically-symmetric mixed state.
    • On the other hand, it can also be prepared via anisotropic processes, such as photo-excitation with polarized light. In that case, the atom will be in an anisotropic state, and the direction of this anisotropy will be dictated by the process that produced it.

    It is extremely tempting to think (as discussed previously e.g. here, here and here, and links therein) that the spherical symmetry of the dynamics (of the nucleus-electron interactions) must imply spherical symmetry of the solutions, but this is obviously wrong $-$ to start with, it would apply equally well to the classical problem! The spherical symmetry implies that, for any anisotropic solution, there exist other, equivalent solutions with complementary anisotropies, but that's it.

  • The hydrogen case is a bit special because the $2p$ shell is an excited state, and the ground state is symmetric. So, in that regard, it is valid to ask: what about the ground states of, say, atomic boron? If all you know is that you have atomic boron in gas phase in its ground state, then indeed you expect a spherically-symmetric mixed state, but this can still be polarized to align all the atoms into the same orientation.

    As a short quip: atoms can have nontrivial shapes, but the fact that we don't know which way those shapes are oriented does not make them spherically symmetric.

  • So, given an atom (perhaps in a fixed excited state), what determines its shape? In short: its term symbol, which tells us its angular momentum characteristics, or, in other words, how it interacts with rotations.

    • The only states with spherical symmetry are those with vanishing total angular momentum, $J=0$. If this is not the case, then there will be two or more states that are physically distinct and which can be related to each other by a rotation.
    • It's important to note that this anisotropy could be in the spin state, such as with the $1s$ ground state of hydrogen. If you want to distinguish the states with isotropic vs anisotropic charge distributions, then you need to look at the total orbital angular momentum, $L$. The charge distribution will be spherically symmetric if and only if $L=0$.

    A good comprehensive source for term symbols of excited states is the Levels section of the NIST ASD.

The reference paper is actually not really wrong, just poorly communicated for this physics audience. If you read the paper this letter is responding to, it becomes clear that the authors are referring to the probability distribution of an incoherent sum of the referenced orbitals. The letter to the editor sloppily writes (emphasis mine)

We therefore must describe the electron as (to use chemically familiar language) a “resonance hybrid” of $2p_x, 2p_y,$ and $2p_z$. In more detail, we write if $ \psi = \frac{1}{\sqrt{3}}(2p_x + 2p_y + 2p_z)$ and this is exactly a spherical distribution, as Johnson and Rettew have shown.

Where, out of context, every physicist assumes that the author specifically means a coherent sum of wavefunctions, and where both the wavefunction and the probability distribution resulting from it are not spherically symmetric. However, I believe (I'm no chemist) to a chemistry audience the common phrase "resonance hybrid" would immediately imply an incoherent superposition of the given states, as there's nothing particularly coherent about normal chemistry. The word "distribution" also hints that something is funny, as it's not typical to call the wavefunction itself a "distribution". Specifically, Johnson and Rettew showed that $\psi_{2p_x}^2 + \psi_{2p_y}^2 + \psi_{2p_z}^2$ is spherically symmetrical, which it is. Since there is basically only one equation in the referenced article, this is clearly what Cohen was referring to. The phrasing of the letter to the editor could clearly should have been clearer here, but good communication does take effort from both sides, especially when the two sides come from different fields where notation is not so well standardized or understood.

For completeness, if a shell is partially filled then there's the possibility of there being a specific angle between the orbits of electrons in the inner shell and the electrons in the outermost shell (think, e.g. of two concentric donuts rotating independently). Even in a mixed state, the outermost electron would be in a mixed state of interacting with various orientations of inner electrons, none of which are independently centric, which suggests that assuming a central field approximation will miss some important physics.

No coherent superposition of 2p orbitals is spherically symmetric. Your example $\frac{1}{\sqrt3}[2p_x+2p_y+2p_z]$ is a 2p orbital pointing in the 111 direction and is not spherical. The proper description is by a diagonal density matrix, which states that the atom is in an incoherent superposition of the three states.