This is an interesting mix of misconceptions:
By definition, a photon is a one-particle state with definite momentum and helicity. In quantum mechanics, the OAM doesn't commute with momentum, and therefore a momentum eigenstate is not an OAM eigenstate. By that logic, individual photons cannot have definite OAM.
No, a correct definition of a photon doesn't need to assign it a definite momentum or helicity, or even a well-defined frequency. The most convenient basis expansions tend to have these properties, but that's not inherent to the definition of a photon.
The short of it is that when you quantize electromagnetism, you start off by finding a suitable basis of vector-valued functions $\mathbf f_n(\mathbf r)$ in which to expand the vector potential as
$$
\mathbf A(\mathbf r,t) = \sum_n\bigg[a_n(t)\mathbf f_n(\mathbf r)+a_n(t)^*\mathbf f_n(\mathbf r)^*\bigg],
$$
where $a_n(t)$ is the generalized coordinate corresponding to the mode $\mathbf f_n(\mathbf r)$, setting things up so that its Poisson bracket with its conjugate is $\{a_n,a_m^*\}=\delta_{mn}$, and then you quantize by replacing $a_n(t)$ with the annihilation operator of that mode, so that a state with a single photon on that mode is $a_n^\dagger|0⟩$.
Now, here's the important thing: there is no requirement that the mode functions $\mathbf f_n(\mathbf r)$ be plane-wave states with circular polarization. It's a convenient choice, but it's not the only possible choice. Photons are excitations of the classical mode in question. Thus, if the classical mode is a plane wave, the photon will have a well-defined linear momentum, but if it is e.g. a Laguerre-Gauss or a Bessel mode, it will have a well-defined orbital angular momentum.
And, exactly as with the basis mode functions themselves, a photon with well-defined angular momentum can be understood as a superposition of photons with well-defined linear momentum (and vice versa), in the same way that you can expand a plane wave in terms of Bessel functions and vice versa. More importantly, this extends to linear combinations of modes with different frequencies: these give single-photon wavepackets, which evolve in time and are not eigenstates of the field hamiltonian, but are still $N=1$ eigenstates of the photon-number operator and therefore equally valid as single-photon states as the single excitations of a monochromatic plane wave.
OK, thus far for the standard description of how to deal with orbital angular momentum within the larger framework of quantum electrodynamics and quantum field theory, or within the more restrictive subsets of those that are often called quantum optics. However, just because you can describe something in a quantum way doesn't mean that you need to, but unfortunately ruling out alternative possible explanations, as you ask in your second question,
Is there some physical example that cannot be explained without assuming that individual photons carry nonzero OAM in addition to spin (or helicity, to be more precise)?
is a rather hard proposition.
However, OAM in this respect is no different to any other degree of freedom of light, and for any experiment that requires photons and a quantum-mechanical description on a given coordinate, you can produce an working experiment built on OAM, from Mandel dips to Bell-inequality violations to quantum cryptography, for which a good review is
G. Molina-Terriza, J.P. Torres and L. Torner. Twisted photons. Nature Phys. 3, 305 (2007).
Now, if you want a direct mechanical detection of the angular momentum carried by a single-photon excitation of an OAM mode, then that's unlikely to be feasible - in the same way that it's likely to be infeasible for the linear momentum of that state, because both are very small and very hard to measure. In that respect, atomic absorption experiments showing modified selection rules are likely to be conceptually sufficient, but I'm unsure whether the experiment has been done as yet.
Finally, if you want a comprehensive yet readable introduction to the subject of the angular momentum of light, I would recommend
R.P. Cameron. On the angular momentum of light. PhD thesis, University of Glasgow (2014).
The situation is simpler than you think. Basically you don't need to do Clebsch-Gordan at all!
The massless irreducible representations of the (proper orthochronous) Lorentz group are $1$-dimensional, labelled by their helicity $h$. Since $h = \pm 1$ are related by parity we group them as the $\pm$ helicities of the photon.
To compute the helicity of multiparticle states you must take a tensor product of the possible helicity configurations. Then you should decompose these into irreducible representations, according to the Clebsch-Gordan procedure.
But in this case, the $+1$ and $-1$ irreps are $1$-dimensional. Therefore their tensor product is automatically an irrep! So there's no Clebsch-Gordan decomposition to do.
Thus we are left with four possibilities, namely
$$1 \otimes 1, \ 1\otimes -1,\ -1 \otimes 1,\ -1 \otimes -1$$
Considering the action of the helicity operator in each representation, you can verify that these have helicity $2,0,0,-2$ respectively.
Combining the $1\otimes -1$ and $-1\otimes 1$ representations symmetrically and antisymmetrically, you reproduce the two $0$ helicity states of Rob's answer.
Let me know if you want some more explicit detail of the calculation!
Best Answer
Here's my haltingly offered answer - I am not sure I am not overlooking a subtlety that you can see but I can't. I'm going to try to answer @dmckee 's take on your question: "How can we motivate the photon spin starting from the classical theory?"
I relate to your worries that slipping between linear and circular polarized base states might just be nonphysical co-ordinate transformations, but I believe that there are three ways that I can think of wherein Nature does show a preference for particular base states, two theoretical and one experimental; these are:
Diagonalisation of Maxwell's Equations by the Riemann-Silberstein Vectors
The Maxwell curl equations in freespace:
$$c\,\left(\begin{array}{c}\nabla \wedge \sqrt{\epsilon_0} \mathbf{E}\\ \nabla \wedge \sqrt{\mu_0} \mathbf{H}\end{array}\right) = \left(\begin{array}{cc}\mathbf{0}_{3\times3}&-\mathbf{1}_{3\times3}\\\mathbf{1}_{3\times3}&\mathbf{0}_{3\times3}\end{array}\right)\left(\begin{array}{c}\partial_t \sqrt{\epsilon_0}\mathbf{E}\\\partial_t \sqrt{\mu_0}\mathbf{H}\end{array}\right)$$
are obviously mutually coupled and can be decoupled by forming the Riemann-Silberstein vectors (arising from the diagonalisation of the block $6\times6$ matrix of ($2\times2$ matrix of four $3\times3$ scalar matrices); we get:
$$i\, \partial_t \mathbf{F}_\pm = \pm c\,\nabla \wedge \mathbf{F}_\pm$$
where:
$$\mathbf{F}_\pm = \frac{1}{\sqrt{2}}\left( \sqrt{\epsilon_0} \mathbf{E} \pm i \sqrt{\mu_0} \mathbf{H}\right)$$
(I apologize for using SI units exclusively - much of my career has been building numerical software and the only way to get such beasts debugged is by making everyone stick to the same units - now I can't think in Planck or natural units at all any more). Now, if $\mathbf{E}$ and $\mathbf{H}$ are real-valued fields, we only need one complex Riemann-Silberstein vector to encode the whole of the Maxwell equations. Equivalent information is encoded into the positive frequency parts alone of the two Riemann-Silberstein vectors. What's really neat about the second approach is that if the light is right circularly polarized, only $\mathbf{F}_+$ is nonzero; if left, only $\mathbf{F}_-$ is nonzero. So the positive frequency parts of the fields are decoupled precisely by splitting them into left and right circularly polarized components, NOT linearly polarized components. This is the first big clue that Nature indeed does show a preference for circularly polarized base states. Also note that the positive frequency (i.e. positive energy) part is the meaningful one if one is to think of Maxwell's equations as the propagation equation for the first quantized photon.
Now, in momentum (Fourier) space, the decoupled Maxwell equations become (we do spatial, not time - Fourier transform of both sides):
$$\mathrm{d}_t \tilde{\mathbf{F}}_{\pm,\mathbf{k}} = \pm c\,\mathbf{k} \wedge \tilde{\mathbf{F}}_{\pm,\mathbf{k}}$$
or, in matrix notation $\mathrm{d}_t \tilde{\mathbf{F}}_{\pm,\mathbf{k}} = \pm c\,\mathbf{K}(\mathbf{k}) \tilde{\mathbf{F}}_{\pm,\mathbf{k}}$ where $\mathbf{K}(\mathbf{k})$ is the $3\times3$ skew-Hermitian matrix corresponding to $\mathbf{k} \wedge$, i.e. the "infinitesimal" rotation in the Lie algebra $\mathfrak{so}(3)$ and the basic solutions are $\tilde{\mathbf{F}}_{\pm,\mathbf{k}} = \exp(c \mathbf{K}(\mathbf{k}) t)\tilde{\mathbf{F}}_{\pm,\mathbf{k}}(0)$, ie vectors spinning at uniform angular speed $\omega = c k$. There is a wonderful aside to the Riemann-Silberstein notation which I come back to at the end of my answer.
If you google Iwo Bialynicki-Birula and his work on the photon wave function, he has heaps more to say about such things. His personal website is http://cft.edu.pl/~birula and all his publications are downloadable therefrom. The particular scaling of the Riemann-Silberstein vectors above is Bialynicki-Birula's, and it means that $|\mathbf{F}_+|^2 + |\mathbf{F}_-|^2$ is the electromagnetic energy density. He defines the pair $(\mathbf{F}_+, \mathbf{F}_-)$, normalised so that $|\mathbf{F}_+|^2 + |\mathbf{F}_-|^2$ becomes a probability density to absorb the photon at a particular point, to be a first quantized photon wave function (without a position observable). There is special, nonlocal inner product to define the Hilbert space and in such a formalism the general Hamiltonian observable is $\hbar\, c\, \mathrm{diag}\left(\nabla\wedge, -\nabla\wedge\right)$. Please also see Arnold Neumaier's pithy summary (here) of a key result in section 7 of Bialynicki-Birula's "Photon wave function" in Progress in Optics 36 V (1996), pp. 245-294 also downloadable from arXiv:quant-ph/0508202. The Hilbert space of Riemann Silberstein vector pairs that Bialynicki-Birula defines is acted on by an irreducible unitary representation, defined by Bialynicki-Birula's observables $\hat{H}$, $\hat{\mathbf{P}}$, $\hat{\mathbf{K}}$ and $\hat{\mathbf{J}}$, of the full Poincaré group presented in the paper. So the two subspaces containing wholly right ($\mathbf{F}_-=\mathbf{0}$) and wholly left-polarized ($\mathbf{F}_+=\mathbf{0}$) states are the "particles" of the theory: you don't get the same thing with other nontrivial linear combinations of light base states (which are not eigenfunctions of the angular momentum observable).
The Classical Angular Momentum Calculation
Now we look at the classical angular momentum. The Wikipedia page on angular momentum of light gives the classical angular momentum as:
$$\frac{\epsilon_0}{2i\omega}\int \left(\mathbf{E}^\ast\wedge\mathbf{E}\right)d^{3}\mathbf{r} +\frac{\epsilon_0}{2i\omega}\sum_{i=x,y,z}\int \left({E^i}^{\ast}\left(\mathbf{r}\wedge\mathbf{\nabla}\right)E^{i}\right)d^{3}\mathbf{r}$$
when the positive frequency part alone of the fields is kept (hence the complex conjugates). The first term is the spin angular momentum, and, rewritten in positive frequency Riemann-Silberstein vectors when everything is roughly paraxial (i.e. near to a plane wave) it reads:
$$\hat{\mathbf{z}}\frac{1}{\omega}\int \left(|\mathbf{F}_+|^2-|\mathbf{F}_-|^2\right)d^{3}\mathbf{r}$$
i.e. $\frac{1}{\omega}$ times the right polarized energy density less the left polarized energy density in the direction of propagation of the light. Orbital angular momentum vanishes in the paraxial limit and so the last equation is the total angular momentum in this case. It is important to recall how this equation is derived: one imagines an electromagnetic field crossing the boundary into a conductive medium and being absorbed there, then one calculates the angular impulse exerted on the medium, exactly analogously with method 3 of the momentum calculation in my answer https://physics.stackexchange.com/a/72688/26076 . The point is that the angular momentum density $\left(|\mathbf{F}_+|^2-|\mathbf{F}_-|^2\right)/\omega$ calculated from this most basic (in the sense of fundamental) Newtonian-Maxwell physics is the difference between the intensities of the circularly polarized base states, not the linear ones. So again Nature shows her preference. This calculation says that the right and left circularly polarized components transfer angular momentum $\pm E/\omega$ in the direction of light propagation, respectively, whenever energy $E$ is absorbed. So now we see that, if the photon has energy $h\nu$, then if a high number of them are to transfer the same angular momentum as classical physics reckons, the photon's angular momentum has to be $\pm h\nu/\omega$ or $\pm \hbar$ in the direction of its propagation for right and left circularly polarized photons, respectively. I'll get to other polarization states in a moment.
General Photon Emission States and Conservation of Angular Momentum by Fluorophores
You ask about one photon emission and whether it is always circular. Absolutely NOT. General one photon states are pure quantum superpositions of circularly polarized one-photon number states. Suppose we have a linearly polarized photon and we work out the results of imparting the number observables: let $a_\pm^\dagger$ be the creation operators for right and left handed polarized states. The linear polarized states are:
$$\psi_x = \frac{1}{\sqrt{2}} \left(\psi_++\psi_-\right)$$ $$\psi_y = \frac{-i}{\sqrt{2}} \left(\psi_+-\psi_-\right)$$
where $\psi_\pm$ are the pure one photon right and left handed polarization states. Then the right and left handed number operators return:
$$n_\pm(\psi_x) = n_\pm(\psi_y) = \left<\psi|a_\pm^\dagger a_\pm | \psi\right> =\frac{1}{2}$$
and the total photon number operator returns:
$$n = \left<\psi|a_+^\dagger a_+ + a_-^\dagger a_-| \psi\right> =1$$
As you likely know, we can always change our one photon Fock state basis by any unitary transformation and the creation and annihilation operators transform likewise. The linear photon creation operators, for example, are $a_x^\dagger=\frac{1}{\sqrt{2}}\left(a_+^\dagger + a_-^\dagger\right)$ and $a_y^\dagger=\frac{-i}{\sqrt{2}}\left(a_+^\dagger - a_-^\dagger\right)$ and imparting the number operator formed from these and their respective Hermitian conjugates would return the result "1 photon" when applied to the respective polarization states and "0" when applied to the orthogonal linear polarization states. The general pure one photon state is of the form:
$$\psi(\alpha, \phi) = \alpha e^{i \frac{\phi}{2}} \psi_+ + \sqrt{1-\alpha^2} e^{-i \frac{\phi}{2}} \psi_-$$
where $\alpha \in [0, 1]$ and $\phi \in [0, 2\,\pi)$ and, from our classical calculation above, we know that its angular momentum must be:
$$\hbar \left<\psi|a_+^\dagger a_+ - a_-^\dagger a_- | \psi\right> = (2\,\alpha^2-1)\,\hbar$$
One photon emissions are always the right pure quantum superpositions that make sure angular momentum is conserved. For example, if a fluorophore absorbs linearly polarized light, and if it does not exert torque on its surroundings before spontaneous emission, the emission must be linearly polarized. Likewise for any other general polarization state absorbed by fluorophores. See the answer https://physics.stackexchange.com/a/73439/26076 for more details, and you also might like to look at the work of Gregorio Weber in the 1950s on the polarization of fluorescence. See see G. Weber, “Rotational Brownian motion and polarization of the fluorescence of solutions,” Adv. Protein Chem. 8, 415–459 (1953) and other works.
My own work in this field is summarised in J. Opt. Soc. Am. B, Vol. 24, No. 6/June 2007 p1369.
The general criterion, in my $\psi(\alpha,\phi)$ notation, where the base states are circularly polarized ones along the direction of propagation, is that if $\alpha_{in},\phi_{in}$ characterize the photon absorbed by an isolated (not transferring angular impulse) fluorophore, and $\alpha_{out},\phi_{out}$ the fluoresced photon, then $\alpha_{in} = \alpha_{out}$ and the phase angles are unrelated. Likewise, when a particle and antiparticle with opposite spins annihilate each other and only two photons are emitted, each characterized by $\alpha_{1},\phi_{1}$ and $\alpha_{2},\phi_{2}$, then $\alpha_{2} = \sqrt{1- \alpha_{1}^2}$ and the phase angles are unrelated and I'm guessing that all values of $\alpha_{1}^2$, $\phi_1$ and $\phi_2$ are equally probable
Torque Exerted on Birefringent Crystals
Now we get to birefringent crystals. It is obvious, but I believe important for this answer, to take heed that birefringent crystals are a way wherein Nature very explicitly tells the difference between linear and circular polarization. It is the linear, NOT the circularly polarized states that are the "eigenmodes" of a birefringent crystal, i.e. linearly polarized fields aligned with the fast and slow axes of the crystal are simply phase delayed and undergo no mixing. Circular polarized states are NOT eigenmodes: they mix in such crystals. The mixing induced by a quarter wave plate - in particular that which happens when the input field is linearly polarized and aligned at 45 degrees to the fast and slow axes therefore imparts a torque on the crystal: it most certainly should be possible to measure that torque and check it against classical calculations: from our calculations above, in this situation there will be a torque $P/\omega$, when the light's power is $P$. A thought experiment: a 100W, 1mm diameter linearly polarized collimated beam passes through a quarter wave plate suspended in a fluid. With infrared light at $193\mathrm{THz}$ (you can get fibre lasers at $193\mathrm{THz}$ outputting hundreds of watts) the torque will be of the order of $10^{-13}\mathrm{Nm}$, if the crystal is a millimeter in diameter and 3 millimeters or so long with a density of $3000\mathrm{kg\,m^{-3}}$, its mass moment of inertia is of the order of $4\times10^{-13}\mathrm{kg \, m^2}$, so this will be an accurately measurable effect (indeed, as the crystal spins and misaligns itself from the 45 degree position, its angular position will fulfill $\mathrm{d}_t^2 \theta = - \frac{1}{2}\Omega^2 \sin(2\theta)$ and we will have a torsional pendulum oscillating at $\Omega /(2\pi) = 0.1\mathrm{Hz}$!).
If there are experiments to observe the transfer of $\hbar$ angular momentum by one photon, then it might be to do with laser tweezer experiments: circularly polarized beams are used to spin things under the microscope in the laser trap and I believe Prof Halina Rubinsztein-Dunlop (http://physics.uq.edu.au/people/halina) was interested some years ago in what happens with such things at very low light levels - I got the impression she had an interest in directly observing one $\hbar$ of angular momentum transfer. She may know of any experiment along these lines.
Aside: More on the Riemann Silberstein Notation
I think you will like this one, Terry. The Riemann-Silberstein vectors are actually the electromagnetic (Maxwell) tensor $F^{\mu\nu}$ in disguise. We can write Maxwell's equations in a quaternion form:
$$\left(c^{-1}\partial_t + \sigma_1 \partial_x + \sigma_2 \partial_y + \sigma_3 \partial_z\right) \,\mathbf{F}_+ = {\bf 0}$$
$$\left(c^{-1}\partial_t - \sigma_1 \partial_x - \sigma_2 \partial_y - \sigma_3 \partial_z\right) \,\mathbf{F}_- = {\bf 0}$$
where $\sigma_j$ are the Pauli spin matrices and the electromagnetic field components are:
$$\begin{array}{lcl}\frac{1}{\sqrt{\epsilon}}\mathbf{F}_\pm &=& \left(\begin{array}{cc}E_z & E_x - i E_y\\E_x + i E_y & -E_z\end{array}\right) \pm i \,c\,\left(\begin{array}{cc}B_z & B_x - i B_y\\B_x + i B_y & -B_z\end{array}\right)\\ & =& E_x \sigma_1 + E_y \sigma_2+E_z\sigma_3 + i\,c\,\left(B_x \sigma_1 + B_y \sigma_2+B_z\sigma_3\right)\end{array}$$
As you know, these the Pauli spin matrices are the imaginary quaternion units reordered. When inertial reference frames are shifted by a proper Lorentz transformation:
$$L = \exp\left(\frac{1}{2}W\right)$$
where:
$$W = \left(\eta^1 + i\theta \chi^1\right) \sigma_1 + \left(\eta^2 + i\theta \chi^2\right) \sigma_2 + \left(\eta^3 + i\theta \chi^3\right) \sigma_3$$
encodes the transformation's rotation angle $\theta$, the direction cosines of $\chi^j$ of its rotation axes and its rapidities $\eta^j$, the entities $\mathbf{F}_\pm$ undergo the spinor map:
$${\bf F} \mapsto L {\bf F} L^\dagger$$
Here, we're actually dealing with the double cover $PSL(2,\mathbb{C})$ of the identity-connected component of the Lorentz group $SO(3,1)$, so we have spinor maps representing Lorentz transformations, just as we must use spinor maps to make a quaternion impart its rotation on a vector. The electromagnetic field is thus a bivector in the Clifford algebra $C\ell_3(\mathbb{R})$, and one neat thing about this notation is that polarization and spin are manifestly Lorentz covariant given the circular polarization base state interpretation of $\mathbf{F}_\pm$ afforded by the discussion above (note that the quaternion units in the de Rham derivatives also transform by the same spinor maps), something that is not so clear in other forms of Maxwell equations. Also, you can convert Bialynicki-Birula's ten observables $\hat{H}$, $\hat{\mathbf{P}}$, $\hat{\mathbf{K}}$ and $\hat{\mathbf{J}}$ also have very clean quaternionic forms.