The angular momentum of light is a bit tricky to get right, particularly because the calculation you have attempted, using a plane wave, is somewhat misleading.
This is because, as you point out, the angular momentum density given by $\mathbf r \times (\mathbf E\times\mathbf B)$ integrates to zero over any reasonable integration volume when the fields are plane waves. On the other hand, it is perfectly possible for a plane-wave field to carry angular momentum if it is circularly polarized, so this looks like a strong contradiction.
The resolution to this contradiction is the fact that plane waves are unphysical, because they are unbounded and carry infinite energy. For a realistic beam you will have a limited spatial extent, and a circularly polarized field will give you a nonzero integral when you include the contributions from the edges of the beam.
However, this is not a particularly satisfactory resolution, and it leaves open the question of whether there is some better formulation for this quantity.
As it turns out, this is indeed the case, and you can do it by some nifty vector calculus manipulations. You know that the total angular momentum is due to the contributions from the linear momentum density, which is proportional to the Poynting vector, so you have
$$
\mathbf J = \frac{1}{\mu_0c^2} \int\mathrm d\mathbf r \: \mathbf r \times (\mathbf E\times\mathbf B).
\tag 1
$$
The trick here is essentially to fix a gauge, i.e., to express $\mathbf B=\nabla \times \mathbf A$ and then work primarily in terms of the vector potential. If you do this, then using some simple vector calculus you can expand the vector-product-of-a-curl in the form
$$
\mathbf E\times\mathbf B
= \mathbf E\times (\nabla\times\mathbf A)
= E_i \nabla \! A_i - (\mathbf E\cdot\nabla)\mathbf A
,
$$
where repeated indices are summed over. With that decomposition in hand, the next level out in the cross products gives you
$$
\mathbf J
= \frac{1}{\mu_0c^2} \int\mathrm d\mathbf r \: \mathbf r \times (E_i \nabla \! A_i )
- \frac{1}{\mu_0c^2} \int\mathrm d\mathbf r \: \mathbf r \times ((\mathbf E\cdot\nabla)\mathbf A)
\tag 2,
$$
which looks more complicated than what we started with, but we're getting there. To simplify this further, you then perform an integration by parts on the second term: using component notation, it reads
\begin{align}
\int\mathrm d\mathbf r \: \mathbf r \times ((\mathbf E\cdot\nabla)\mathbf A)
& =
\epsilon_{ijk}\mathbf e_i\int\mathrm d\mathbf r \: x_j E_n\frac{\partial A_k}{\partial x_n}
\\ & =
\epsilon_{ijk}\mathbf e_i\int\mathrm d\mathbf r \left[
\frac{\partial }{\partial x_n}\left(x_j E_nA_k\right)
- \frac{\partial x_j}{\partial x_n}E_nA_k
- x_j \frac{\partial E_n}{\partial x_n}A_k
\right]
.
\end{align}
Here the third term vanishes in vacuum, since the derivative is simply $\nabla\cdot\mathbf{E}$, and you also want the first term, i.e. the total derivative, to vanish: for a finite volume $V$, you have
\begin{align}
\int_V\mathrm d\mathbf r
\frac{\partial }{\partial x_n}\left(x_j E_nA_k\right)
& =
\int_V\mathrm d\mathbf r \:
\nabla \cdot\left(x_j A_k \mathbf E\right)
=
\int_{\partial V}\mathrm d\mathbf S \cdot\left(x_j A_k \mathbf E\right)
,
\tag 3
\end{align}
by the divergence theorem, and for a physical field with a bounded support for $\mathbf E$ and $\mathbf A$ the boundary terms should vanish when $\partial V$ is far away from the support of the fields.
What you're left with, then, is a drastically simplified expression:
\begin{align}
\int\mathrm d\mathbf r \: \mathbf r \times ((\mathbf E\cdot\nabla)\mathbf A)
& =
-\epsilon_{ijk}\mathbf e_i\int\mathrm d\mathbf r \:
\delta_{jn}E_nA_k
=
-\int\mathrm d\mathbf r \: \mathbf E\times \mathbf A
,
\end{align}
which you can then put back into $(2)$ to get
$$
\mathbf J
= \frac{1}{\mu_0c^2} \int\mathrm d\mathbf r \: (E_i (\mathbf r \times\nabla) A_i )
+ \frac{1}{\mu_0c^2} \int\mathrm d\mathbf r \: \mathbf E\times \mathbf A
\tag 4.
$$
Here the first term is known as the orbital angular momentum of the fields, and the second term gives you the spin angular momentum; this is backed up by the fact that a circularly polarized field will contribute to the second term and 'optical vortex' beams, like Laguerre-Gauss and similar beamshapes, will contribute to the first term.
Moreover, as if by magic, if you now apply these yardsticks to your original plane-wave state, you'll get zero for the orbital angular momentum, but now you can get a nonzero spin contribution from the second term if the light is circularly polarized. This should seem contradictory, given that our $\mathbf J$ was originally zero, but if you look closely you'll see that our step of neglecting the total derivative in $(3)$ does not really apply to a plane wave, since the fields do not vanish at infinity.
So... which is right? Well, we want our plane waves to model physical fields, but if the unphysical characteristics get us into trouble, then the model is wrong, and we should keep the derivations that work for physical fields. Plane waves should always be handled with care, but in this case $(4)$ is the correct expression.
For more a more classic reference on this, see §II of
What is spin? Hans C. Ohanian, Am. J. Phys. 54 500 (1986), eprint
though the Cameron et al. paper you cite has a much more sophisticated approach.
