# Optical momentum distributions in monochromatic, isotropic random vector fields

Titouan Gadeyne, Mark R Dennis

► To cite this version:

Titouan Gadeyne, Mark R Dennis. Optical momentum distributions in monochromatic, isotropic random vector fields. Journal of Optics, 2024, 26 (6), pp.065604. 10.1088/2040-8986/ad3b14. cea-04569117

HAL Id: cea-04569117

https://cea.hal.science/cea-04569117

Submitted on 6 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Optical momentum distributions in monochromatic, isotropic random vector fields

Titouan Gadeyne, Mark R Dennis

► To cite this version:

Titouan Gadeyne, Mark R Dennis. Optical momentum distributions in monochromatic, isotropic random vector fields. Journal of Optics, 2024, 26 (6), pp.065604. 10.1088/2040-8986/ad3b14. cea-04569117

HAL Id: cea-04569117

https://cea.hal.science/cea-04569117v1

Submitted on 6 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optical momentum distributions in monochromatic, isotropic random vector fields

Titouan Gadeyne<fnref n="1" />,<fnref n="2" />,<fnref n="4" /> ID and Mark R Dennis<fnref n="1" />,<fnref n="3" />,* ID

E-mail: m.r.dennis@bham.ac.uk

Received 24 November 2023, revised 17 February 2024

Accepted for publication 5 April 2024

Published 3 May 2024

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

We investigate the decomposition of the electromagnetic Poynting momentum density in three-dimensional random monochromatic fields into orbital and spin parts, using analytical and numerical methods. In sharp contrast with the paraxial case, the orbital and spin momenta in isotropic random fields are found to be identically distributed in magnitude, increasing the discrepancy between the Poynting and orbital pictures of energy flow. Spatial correlation functions reveal differences in the generic organization of the optical momenta in complex natural light fields, with the orbital current typically forming broad channels of unidirectional flow, and the spin current manifesting larger vorticity and changing direction over subwavelength distances. These results are extended to random fields with pure helicity, in relation to the inclusion of electric-magnetic democracy in the definition of optical momenta.

Supplementary material for this article is available online

Keywords: statistical optics, optical momentum, Poynting vector, spin momentum, orbital momentum

## 1. Introduction

Conservation of electromagnetic (EM) energy is determined by the well-known theorem of Poynting [1]: in the absence of charges, the rate of change of EM energy density is equal to the divergence of the Poynting vector $P = \mathcal{E} \times \mathcal{H}$, the cross product of the electric and magnetic fields. By analogy with other continuity equations, it is customary to interpret $P$ as the direction and magnitude of EM energy flow [2, 3]. However,

this choice often fails to produce an intuitive picture, even in seemingly elementary situations: for instance, the Poynting vector for two crossed plane waves [4] or in a single evanescent surface wave [5, 6] exhibits a counterintuitive component perpendicular to the direction of propagation. Similarly, the (time-averaged) radiation pressure exerted by an optical field on a subwavelength probe particle is generally not proportional to the Poynting vector [7, 8].

Divided by $c^2$, the Poynting vector also defines the linear momentum density of the EM field. It is now well understood that in monochromatic fields, the time-averaged linear momentum $P$ is the sum of orbital $P_O$ and spin $P_S$ parts, respectively generating the orbital and spin angular momenta [9–11]. In [10], these vector fields were dubbed optical currents. This Poynting vector splitting has deep foundations, as the orbital momentum in fact corresponds to the canonical

<figure><img src="image_3.png" /></figure>

Original Content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ School of Physics and Astronomy, University of Birmingham, Birmingham B15 2TT, United Kingdom</footnote>

<footnote n="2">² Département de Chimie, École Normale Supérieure, PSL University, 75005 Paris, France</footnote>

<footnote n="3">³ EPSRC Centre for Doctoral Training in Topological Design, University of Birmingham, Birmingham B15 2TT, United Kingdom</footnote>

<footnote n="4">⁴ Université Paris-Saclay, CEA, CNRS, LIDYL 91191, Gif-sur-Yvette, France</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
# Optical momentum distributions in monochromatic, isotropic random vector fields

Titouan Gadeyne<sup>1,2,4</sup> and Mark R Dennis<sup>1,3,*</sup>

E-mail: m.r.dennis@bham.ac.uk

Received 24 November 2023, revised 17 February 2024

Accepted for publication 5 April 2024

Published 3 May 2024

<figure><img src="image_3.png" /></figure>

CrossMark

## Abstract

We investigate the decomposition of the electromagnetic Poynting momentum density in three-dimensional random monochromatic fields into orbital and spin parts, using analytical and numerical methods. In sharp contrast with the paraxial case, the orbital and spin momenta in isotropic random fields are found to be identically distributed in magnitude, increasing the discrepancy between the Poynting and orbital pictures of energy flow. Spatial correlation functions reveal differences in the generic organization of the optical momenta in complex natural light fields, with the orbital current typically forming broad channels of unidirectional flow, and the spin current manifesting larger vorticity and changing direction over subwavelength distances. These results are extended to random fields with pure helicity, in relation to the inclusion of electric-magnetic democracy in the definition of optical momenta.

Supplementary material for this article is available online

Keywords: statistical optics, optical momentum, Poynting vector, spin momentum, orbital momentum

## 1. Introduction

Conservation of electromagnetic (EM) energy is determined by the well-known theorem of Poynting [1]: in the absence of charges, the rate of change of EM energy density is equal to the divergence of the Poynting vector $P = \mathcal{E} \times \mathcal{H}$, the cross product of the electric and magnetic fields. By analogy with other continuity equations, it is customary to interpret $P$ as the direction and magnitude of EM energy flow [2, 3]. However,

this choice often fails to produce an intuitive picture, even in seemingly elementary situations: for instance, the Poynting vector for two crossed plane waves [4] or in a single evanescent surface wave [5, 6] exhibits a counterintuitive component perpendicular to the direction of propagation. Similarly, the (time-averaged) radiation pressure exerted by an optical field on a subwavelength probe particle is generally not proportional to the Poynting vector [7, 8].

Divided by $c^2$, the Poynting vector also defines the linear momentum density of the EM field. It is now well understood that in monochromatic fields, the time-averaged linear momentum $P$ is the sum of orbital $P_O$ and spin $P_S$ parts, respectively generating the orbital and spin angular momenta [9–11]. In [10], these vector fields were dubbed optical currents. This Poynting vector splitting has deep foundations, as the orbital momentum in fact corresponds to the canonical

<figure><img src="image_4.png" /></figure>

Original Content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ School of Physics and Astronomy, University of Birmingham, Birmingham B15 2TT, United Kingdom</footnote>

<footnote n="2">² Département de Chimie, École Normale Supérieure, PSL University, 75005 Paris, France</footnote>

<footnote n="3">³ EPSRC Centre for Doctoral Training in Topological Design, University of Birmingham, Birmingham B15 2TT, United Kingdom</footnote>

<footnote n="4">⁴ Université Paris-Saclay, CEA, CNRS, LIDYL 91191, Gif-sur-Yvette, France</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
momentum density derived from application of Noether's theorem to translational invariance in the relativistic field theory formulation of electromagnetism [12, 13]. The orbital momentum correctly accounts for the radiation pressure on dipole particles, and can provide a more intuitive picture of energy flow than the Poynting vector in the situations mentioned above. In the field theory framework, the spin momentum corresponds to a term introduced by Belinfante [14] to guarantee symmetry and gauge-invariance to the EM stress-energy tensor [12], which does not contribute to the total linear momentum of the field when integrated over space. As such, the Belinfante spin momentum is often described as a 'virtual' quantity introduced for theoretical reasons. Nevertheless, this spin momentum has recently been evidenced experimentally, by measuring the extraordinary optical force it induced on a nano-cantilever [15]. Importantly, the couplings to the orbital and spin parts of the Poynting vector differ by orders of magnitude, highlighting their distinct physical nature.

Recent experimental and theoretical studies have thus demonstrated striking differences between the Poynting, orbital, and spin momenta, and continue to redefine our views of EM energy flow and optical forces [8, 11, 16]. Still, they have so far been limited to rather elementary, highly symmetric fields, with geometries optimized to best showcase the differences between the three optical currents. In particular, the focus on light beams ensures the orbital momentum, related to the propagation k-vector of the field, is typically large, whereas the spin momentum, related to the rate of change of polarization with position, is smaller.

In this work, we explore generic features of these optical currents, to build insight into their organization in natural, isotropic light fields: what are their properties when many independent waves interfere with no particular symmetries? To this end, we investigate their behavior in completely coherent, monochromatic, isotropic random EM vector fields, a convenient statistical model of 3D EM fields specified by only one physical parameter, the wavelength $\lambda$. Strikingly, we will see that in this model, the magnitudes of the spin and orbital currents have the same probability distribution, but that the two vector fields have different spatial correlations: the apparent weakness of the spin current in this setting is due to its failure to organize coherent correlated vector structures over large distances in space, unlike the orbital current (and Poynting vector itself). We demonstrate these facts using analytical statistics and numerical simulations for the vector random wave model.

## 2. Theoretical framework

### 2.1. Poynting, orbital and spin momenta

We work in units where $\epsilon_0 = \mu_0 = c = 1$. In a monochromatic field with frequency $\omega = ck = c(2\pi/\lambda)$, the electric

and magnetic fields are represented by complex-valued vector fields $\mathbf{E}$ and $\mathbf{H}$, where the physical fields are $\mathcal{E} = \text{Re}\{\mathbf{E}e^{-i\omega t}\}$ and $\mathcal{H} = \text{Re}\{\mathbf{H}e^{-i\omega t}\}$. The temporal cycle-averaged Poynting momentum can be written

$$ \mathbf{P} = \frac{1}{2} \text{Re}\{\mathbf{E}^* \times \mathbf{H}\}. \quad (1) $$

Using Maxwell's equations and the vector identity $\mathbf{A} \times (\nabla \times \mathbf{B}) = \mathbf{A} \cdot (\nabla)\mathbf{B} - (\mathbf{A} \cdot \nabla)\mathbf{B}$ (where we use the customary notation $[\mathbf{A} \cdot (\nabla)\mathbf{B}]_i = \sum_j A_j \partial_i B_j$), the Poynting momentum can be split into a sum of orbital and spin momenta [10],

$$ \begin{aligned} \mathbf{P} &= \frac{1}{2\omega} \text{Im}\{\mathbf{E}^* \cdot (\nabla) \mathbf{E}\} - \frac{1}{2\omega} \text{Im}\{(\mathbf{E}^* \cdot \nabla) \mathbf{E}\} \quad (2) \\ &= \mathbf{P}_O^E + \mathbf{P}_S^E. \end{aligned} $$

This splitting is not unique: expressing the Poynting momentum using the electric (resp. magnetic) field only, we obtain the electric-biased (resp. magnetic-biased) momenta $\mathbf{P}_{O,S}^E$ (resp. $\mathbf{P}_S^H$) as in (2). But another option, retaining the electric-magnetic symmetry of Maxwell's equations for free fields, is to take the mean of the two representations,

$$ \begin{aligned} \mathbf{P} &= \frac{1}{4\omega} \text{Im}\{\mathbf{E}^* \cdot (\nabla) \mathbf{E} + \mathbf{H}^* \cdot (\nabla) \mathbf{H}\} \\ &\quad - \frac{1}{4\omega} \text{Im}\{(\mathbf{E}^* \cdot \nabla) \mathbf{E} + (\mathbf{H}^* \cdot \nabla) \mathbf{H}\} \quad (3) \\ &= \mathbf{P}_O^{EH} + \mathbf{P}_S^{EH}, \end{aligned} $$

producing the so-called democratic (or dual) momenta $\mathbf{P}_O^{EH}$, $\mathbf{P}_S^{EH}$ [10]. In general non-paraxial fields, these are all distinct quantities, and in monochromatic fields, their definition is unambiguous (otherwise the splitting is gauge-dependent). An interesting situation arises in fields with pure helicity, consisting only of circularly-polarized plane wave components of same handedness: such fields satisfy $\mathbf{E} = \pm i\mathbf{H}$, such that all biased and democratic quantities become identical.

The dual formulation of electromagnetism that treats electric and magnetic fields equally has many attractive features when working with free fields, in the absence of matter [13, 17] — for instance, democratic momenta naturally split into two independent parts associated with components of opposite helicity [18]. However, experimental measurements require material probes, which typically do not respond identically to electric and magnetic fields: a common example is the radiation pressure on a subwavelength particle responding in the electric dipole approximation, which is proportional to the electric-biased orbital momentum $\mathbf{P}_O^E$ only. We therefore choose to center our discussion on electric-biased quantities, and devote section 3.3 to observations on democratic momenta and pure helicity fields for which the distinction vanishes.

For reference, we briefly recall the typical magnitudes of the three momenta in a paraxial beam [9, 11] propagating along $z$, for which the field is approximately $\mathbf{E} \approx e^{ikz}(E_x, E_y, 0)$
momentum density derived from application of Noether's theorem to translational invariance in the relativistic field theory formulation of electromagnetism [12, 13]. The orbital momentum correctly accounts for the radiation pressure on dipole particles, and can provide a more intuitive picture of energy flow than the Poynting vector in the situations mentioned above. In the field theory framework, the spin momentum corresponds to a term introduced by Belinfante [14] to guarantee symmetry and gauge-invariance to the EM stress-energy tensor [12], which does not contribute to the total linear momentum of the field when integrated over space. As such, the Belinfante spin momentum is often described as a 'virtual' quantity introduced for theoretical reasons. Nevertheless, this spin momentum has recently been evidenced experimentally, by measuring the extraordinary optical force it induced on a nano-cantilever [15]. Importantly, the couplings to the orbital and spin parts of the Poynting vector differ by orders of magnitude, highlighting their distinct physical nature.

Recent experimental and theoretical studies have thus demonstrated striking differences between the Poynting, orbital, and spin momenta, and continue to redefine our views of EM energy flow and optical forces [8, 11, 16]. Still, they have so far been limited to rather elementary, highly symmetric fields, with geometries optimized to best showcase the differences between the three optical currents. In particular, the focus on light beams ensures the orbital momentum, related to the propagation k-vector of the field, is typically large, whereas the spin momentum, related to the rate of change of polarization with position, is smaller.

In this work, we explore generic features of these optical currents, to build insight into their organization in natural, isotropic light fields: what are their properties when many independent waves interfere with no particular symmetries? To this end, we investigate their behavior in completely coherent, monochromatic, isotropic random EM vector fields, a convenient statistical model of 3D EM fields specified by only one physical parameter, the wavelength $\lambda$. Strikingly, we will see that in this model, the magnitudes of the spin and orbital currents have the same probability distribution, but that the two vector fields have different spatial correlations: the apparent weakness of the spin current in this setting is due to its failure to organize coherent correlated vector structures over large distances in space, unlike the orbital current (and Poynting vector itself). We demonstrate these facts using analytical statistics and numerical simulations for the vector random wave model.

## 2. Theoretical framework

### 2.1. Poynting, orbital and spin momenta

We work in units where $\epsilon_0 = \mu_0 = c = 1$. In a monochromatic field with frequency $\omega = ck = c(2\pi/\lambda)$, the electric

and magnetic fields are represented by complex-valued vector fields $\mathbf{E}$ and $\mathbf{H}$, where the physical fields are $\mathcal{E} = \text{Re}\{\mathbf{E}e^{-i\omega t}\}$ and $\mathcal{H} = \text{Re}\{\mathbf{H}e^{-i\omega t}\}$. The temporal cycle-averaged Poynting momentum can be written

$$ \mathbf{P} = \frac{1}{2} \text{Re}\{\mathbf{E}^* \times \mathbf{H}\}. \quad (1) $$

Using Maxwell's equations and the vector identity $\mathbf{A} \times (\nabla \times \mathbf{B}) = \mathbf{A} \cdot (\nabla)\mathbf{B} - (\mathbf{A} \cdot \nabla)\mathbf{B}$ (where we use the customary notation $[\mathbf{A} \cdot (\nabla)\mathbf{B}]_i = \sum_j A_j \partial_i B_j$), the Poynting momentum can be split into a sum of orbital and spin momenta [10],

$$ \begin{aligned} \mathbf{P} &= \frac{1}{2\omega} \text{Im}\{\mathbf{E}^* \cdot (\nabla) \mathbf{E}\} - \frac{1}{2\omega} \text{Im}\{(\mathbf{E}^* \cdot \nabla) \mathbf{E}\} \\ &= \mathbf{P}_O^E + \mathbf{P}_S^E. \end{aligned} \quad (2) $$

This splitting is not unique: expressing the Poynting momentum using the electric (resp. magnetic) field only, we obtain the electric-biased (resp. magnetic-biased) momenta $\mathbf{P}_{O,S}^E$ (resp. $\mathbf{P}_{S,O}^H$) as in (2). But another option, retaining the electric-magnetic symmetry of Maxwell's equations for free fields, is to take the mean of the two representations,

$$ \begin{aligned} \mathbf{P} &= \frac{1}{4\omega} \text{Im}\{\mathbf{E}^* \cdot (\nabla) \mathbf{E} + \mathbf{H}^* \cdot (\nabla) \mathbf{H}\} \\ &\quad - \frac{1}{4\omega} \text{Im}\{(\mathbf{E}^* \cdot \nabla) \mathbf{E} + (\mathbf{H}^* \cdot \nabla) \mathbf{H}\} \\ &= \mathbf{P}_O^{EH} + \mathbf{P}_S^{EH}, \end{aligned} \quad (3) $$

producing the so-called democratic (or dual) momenta $\mathbf{P}_O^{EH}$, $\mathbf{P}_S^{EH}$ [10]. In general non-paraxial fields, these are all distinct quantities, and in monochromatic fields, their definition is unambiguous (otherwise the splitting is gauge-dependent). An interesting situation arises in fields with pure helicity, consisting only of circularly-polarized plane wave components of same handedness: such fields satisfy $\mathbf{E} = \pm i\mathbf{H}$, such that all biased and democratic quantities become identical.

The dual formulation of electromagnetism that treats electric and magnetic fields equally has many attractive features when working with free fields, in the absence of matter [13, 17] — for instance, democratic momenta naturally split into two independent parts associated with components of opposite helicity [18]. However, experimental measurements require material probes, which typically do not respond identically to electric and magnetic fields: a common example is the radiation pressure on a subwavelength particle responding in the electric dipole approximation, which is proportional to the electric-biased orbital momentum $\mathbf{P}_O^E$ only. We therefore choose to center our discussion on electric-biased quantities, and devote section 3.3 to observations on democratic momenta and pure helicity fields for which the distinction vanishes.

For reference, we briefly recall the typical magnitudes of the three momenta in a paraxial beam [9, 11] propagating along $z$, for which the field is approximately $\mathbf{E} \approx e^{ikz}(E_x, E_y, 0)$
with transverse amplitude and polarization profiles $E_x(x,y)$ and $E_y(x,y)$ varying over a lengthscale $W \gg \lambda$, on the order of the beam waist. We find that $\mathbf{P}$ and $\mathbf{P}_O^E$ are mostly longitudinal, that $\mathbf{P}_S^E$ is purely transverse, and the following orders of magnitude

$$ 
\begin{aligned}
|\mathbf{P}| & \sim E^2, \\
|\mathbf{P}_O^E| & \sim \frac{1}{\omega} \mathrm{Im}\{E^* \partial_z E\} \sim \frac{k}{\omega} E^2 \sim E^2, \\
|\mathbf{P}_S^E| & \sim \frac{1}{\omega} \mathrm{Im}\{E^* \partial_x E\} \sim \frac{\lambda}{W} E^2 \ll E^2.
\end{aligned}
$$

We conclude that in a regular optical beam, orbital and Poynting momenta are closely aligned, the spin momentum being small in comparison.

## 2.2. Gaussian random optical fields

We model generic, natural EM light fields as superpositions of $N \to +\infty$ plane waves, with randomly sampled propagation directions taken uniformly over the sphere of directions, random elliptical polarizations, and random overall phases. This construction aims to portray an unprepared field, akin to ambient light in a room, or thermal black-body radiation, composed of many waves emitted from independent points or having scattered off various surfaces, producing a field with no particular symmetries or preferred directions, but statistically homogeneous and isotropic. This approach builds on the long history of the study of speckle patterns [19–21] and statistical properties of spatial structures in random light fields [22–26], which revealed salient features underlying the organization of all EM fields. Physically and geometrically, these random vector fields are very different from directed, paraxial optical beams.

The complex electric and magnetic fields can be parameterized as follows,

$$ \mathbf{E} = \sqrt{\frac{2}{N}} \sum_{n=1}^{N} e^{i\mathbf{k}_n \cdot \mathbf{r} + i\psi_n} \left[ e^{i\alpha_n/2} \cos \frac{\beta_n}{2} \mathbf{e}_+ (\mathbf{k}_n) + e^{-i\alpha_n/2} \sin \frac{\beta_n}{2} \mathbf{e}_- (\mathbf{k}_n) \right], \quad (4a) $$

$$ \mathbf{H} = \sqrt{\frac{2}{N}} \sum_{n=1}^{N} e^{i\mathbf{k}_n \cdot \mathbf{r} + i\psi_n} \frac{1}{i} \left[ e^{i\alpha_n/2} \cos \frac{\beta_n}{2} \mathbf{e}_+ (\mathbf{k}_n) - e^{-i\alpha_n/2} \sin \frac{\beta_n}{2} \mathbf{e}_- (\mathbf{k}_n) \right], \quad (4b) $$

where the sum runs over the $N \gg 1$ plane waves, with wavevectors $\mathbf{k}_n$ sampled uniformly on the sphere of directions with spherical angles $(\theta_n, \phi_n)$ and identical magnitudes $k$, and polarizations sampled uniformly on the Poincaré sphere with angles $(\beta_n, \alpha_n)$, and uniformly sampled global phases $\psi_n$. $\mathbf{e}_\pm(\mathbf{k}) = [\mathbf{e}_1 \pm i\mathbf{e}_2]/\sqrt{2}$ are helicity basis vectors, with $\{\mathbf{e}_1, \mathbf{e}_2\}$

a basis of two real orthogonal unit vectors transverse to $\mathbf{k}$ (see the SI for explicit expressions and alternative parameterizations). We introduce the following notation for the real and imaginary parts of the fields [2, 26]

$$ \mathbf{E} = \mathbf{p}^E + i\mathbf{q}^E, \quad \mathbf{H} = \mathbf{p}^H + i\mathbf{q}^H, $$

since statistics are convenient with real quantities only. Ensemble-averaging over many random fields is denoted by brackets and amounts to integrating over the five random angles

$$ \langle \bullet \rangle = \prod_{n=1}^{N} \left[ \frac{1}{32\pi^3} \int_0^\pi \sin\theta_n d\theta_n \int_0^{2\pi} d\phi_n \times \int_0^\pi \sin\beta_n d\beta_n \int_0^{2\pi} d\alpha_n \int_0^{2\pi} d\psi_n \right] \bullet. \quad (5) $$

From the definitions above, it can be seen that any component of the real or imaginary part of a field is a sum of $N$ real-valued, identically distributed random variables. The central limit theorem ensures that in the limit $N \to +\infty$, each component is a real random variable obeying Gaussian statistics [19, 20]. The same reasoning holds for all derivatives of the components. In our case these variables are all centered, hence we only require their variances and correlations to fully describe the statistics. They are obtained by direct integration using (5), and are tabulated in the SI. With these provided, an ensemble average rewrites as an integral over a set of $M$ Gaussian random variables $\mathbf{u} = (p_x^E, p_y^E, \dots)$

$$ \langle \bullet \rangle = \sqrt{\frac{\det \Sigma^{-1}}{(2\pi)^M}} \int \dots \int d^M \mathbf{u} \exp \left\{ -\frac{\mathbf{u}^\mathrm{T} \Sigma^{-1} \mathbf{u}}{2} \right\} \bullet, \quad (6) $$

where $\Sigma$ is the covariance matrix, with $\Sigma_{ij} = \langle u_i u_j \rangle$. Useful formulae and strategies for computing averages are further described in the SI, and can be found in [22, 23, 25–27].

## 2.3. Spatial correlation functions

To investigate local order in the spatial organization of the optical currents, we will average products of vector components at two different positions in space. The statistical, directional correlators in random Gaussian vector fields, here representing EM waves, are analogous to those used in the theory of isotropic turbulence in fluids [28]. For a homogeneous random vector field $\mathbf{v}$ to be isotropic requires the two-point correlation tensor to have the form

$$ \langle v_i(\mathbf{0}) v_j(\mathbf{r}) \rangle = [f(r) - g(r)] \frac{r_i r_j}{r^2} + g(r) \delta_{ij}, $$

where $f$ and $g$ are scalar functions depending only on the magnitude $r = |\mathbf{r}|$ of the separation vector. They respectively
with transverse amplitude and polarization profiles $E_x(x,y)$ and $E_y(x,y)$ varying over a lengthscale $W \gg \lambda$, on the order of the beam waist. We find that $\mathbf{P}$ and $\mathbf{P}_O^E$ are mostly longitudinal, that $\mathbf{P}_S^E$ is purely transverse, and the following orders of magnitude

$$ 
\begin{aligned}
|\mathbf{P}| & \sim E^2, \\
|\mathbf{P}_O^E| & \sim \frac{1}{\omega} \mathrm{Im}\{E^* \partial_z E\} \sim \frac{k}{\omega} E^2 \sim E^2, \\
|\mathbf{P}_S^E| & \sim \frac{1}{\omega} \mathrm{Im}\{E^* \partial_x E\} \sim \frac{\lambda}{W} E^2 \ll E^2.
\end{aligned}
$$

We conclude that in a regular optical beam, orbital and Poynting momenta are closely aligned, the spin momentum being small in comparison.

## 2.2. Gaussian random optical fields

We model generic, natural EM light fields as superpositions of $N \to +\infty$ plane waves, with randomly sampled propagation directions taken uniformly over the sphere of directions, random elliptical polarizations, and random overall phases. This construction aims to portray an unprepared field, akin to ambient light in a room, or thermal black-body radiation, composed of many waves emitted from independent points or having scattered off various surfaces, producing a field with no particular symmetries or preferred directions, but statistically homogeneous and isotropic. This approach builds on the long history of the study of speckle patterns [19–21] and statistical properties of spatial structures in random light fields [22–26], which revealed salient features underlying the organization of all EM fields. Physically and geometrically, these random vector fields are very different from directed, paraxial optical beams.

The complex electric and magnetic fields can be parameterized as follows,

$$ \mathbf{E} = \sqrt{\frac{2}{N}} \sum_{n=1}^{N} e^{i\mathbf{k}_n \cdot \mathbf{r} + i\psi_n} \left[ e^{i\alpha_n/2} \cos \frac{\beta_n}{2} \mathbf{e}_+ (\mathbf{k}_n) + e^{-i\alpha_n/2} \sin \frac{\beta_n}{2} \mathbf{e}_- (\mathbf{k}_n) \right], \quad (4a) $$

$$ \mathbf{H} = \sqrt{\frac{2}{N}} \sum_{n=1}^{N} e^{i\mathbf{k}_n \cdot \mathbf{r} + i\psi_n} \frac{1}{i} \left[ e^{i\alpha_n/2} \cos \frac{\beta_n}{2} \mathbf{e}_+ (\mathbf{k}_n) - e^{-i\alpha_n/2} \sin \frac{\beta_n}{2} \mathbf{e}_- (\mathbf{k}_n) \right], \quad (4b) $$

where the sum runs over the $N \gg 1$ plane waves, with wavevectors $\mathbf{k}_n$ sampled uniformly on the sphere of directions with spherical angles $(\theta_n, \phi_n)$ and identical magnitudes $k$, and polarizations sampled uniformly on the Poincaré sphere with angles $(\beta_n, \alpha_n)$, and uniformly sampled global phases $\psi_n$. $\mathbf{e}_\pm(\mathbf{k}) = [\mathbf{e}_1 \pm i\mathbf{e}_2]/\sqrt{2}$ are helicity basis vectors, with $\{\mathbf{e}_1, \mathbf{e}_2\}$

a basis of two real orthogonal unit vectors transverse to $\mathbf{k}$ (see the SI for explicit expressions and alternative parameterizations). We introduce the following notation for the real and imaginary parts of the fields [2, 26]

$$ \mathbf{E} = \mathbf{p}^E + i\mathbf{q}^E, \quad \mathbf{H} = \mathbf{p}^H + i\mathbf{q}^H, $$

since statistics are convenient with real quantities only. Ensemble-averaging over many random fields is denoted by brackets and amounts to integrating over the five random angles

$$ \langle \bullet \rangle = \prod_{n=1}^{N} \left[ \frac{1}{32\pi^3} \int_0^\pi \sin\theta_n d\theta_n \int_0^{2\pi} d\phi_n \times \int_0^\pi \sin\beta_n d\beta_n \int_0^{2\pi} d\alpha_n \int_0^{2\pi} d\psi_n \right] \bullet. \quad (5) $$

From the definitions above, it can be seen that any component of the real or imaginary part of a field is a sum of $N$ real-valued, identically distributed random variables. The central limit theorem ensures that in the limit $N \to +\infty$, each component is a real random variable obeying Gaussian statistics [19, 20]. The same reasoning holds for all derivatives of the components. In our case these variables are all centered, hence we only require their variances and correlations to fully describe the statistics. They are obtained by direct integration using (5), and are tabulated in the SI. With these provided, an ensemble average rewrites as an integral over a set of $M$ Gaussian random variables $\mathbf{u} = (p_x^E, p_y^E, \dots)$

$$ \langle \bullet \rangle = \sqrt{\frac{\det \Sigma^{-1}}{(2\pi)^M}} \int \dots \int d^M \mathbf{u} \exp \left\{ -\frac{\mathbf{u}^\top \Sigma^{-1} \mathbf{u}}{2} \right\} \bullet, \quad (6) $$

where $\Sigma$ is the covariance matrix, with $\Sigma_{ij} = \langle u_i u_j \rangle$. Useful formulae and strategies for computing averages are further described in the SI, and can be found in [22, 23, 25–27].

## 2.3. Spatial correlation functions

To investigate local order in the spatial organization of the optical currents, we will average products of vector components at two different positions in space. The statistical, directional correlators in random Gaussian vector fields, here representing EM waves, are analogous to those used in the theory of isotropic turbulence in fluids [28]. For a homogeneous random vector field $\mathbf{v}$ to be isotropic requires the two-point correlation tensor to have the form

$$ \langle v_i(\mathbf{0}) v_j(\mathbf{r}) \rangle = [f(r) - g(r)] \frac{r_i r_j}{r^2} + g(r) \delta_{ij}, $$

where $f$ and $g$ are scalar functions depending only on the magnitude $r = |\mathbf{r}|$ of the separation vector. They respectively
describe longitudinal and lateral autocorrelations of a given vector component

$$f(r) = \langle v_i(\mathbf{0}) v_i(\mathbf{r}e_i) \rangle, \quad g(r) = \langle v_i(\mathbf{0}) v_i(\mathbf{r}e_j) \rangle \quad (i \neq j).$$

where the separation vector $re_i$ is taken along some chosen direction $i = x, y, z$. If, in addition, the field is solenoidal ($\nabla \cdot \mathbf{v} = 0$), $f$ and $g$ are related, such that the full correlation tensor can be determined from, for example, the longitudinal correlation function $f$ only,

$$\langle v_i(\mathbf{0}) v_j(\mathbf{r}) \rangle = -\frac{r f'(r)}{2} \frac{r_i r_j}{r^2} + \delta_{ij} \left[ f(r) + \frac{r f'(r)}{2} \right]. \quad (7)$$

Since there are no charges in the model field, cycle-averaging Poynting's theorem yields $\nabla \cdot \mathbf{P} = 0$. As the spin momentum itself is the curl of a vector field [10], it is divergenceless $\nabla \cdot \mathbf{P}_S = 0$, and consequently we also have $\nabla \cdot \mathbf{P}_O = 0$. Hence all momenta are isotropic homogeneous solenoidal random fields, to which the above results apply. They also apply to the complex electric and magnetic fields themselves. In our calculations, we will be able to express all correlation functions using the longitudinal and lateral autocorrelation functions of the electric field [29], that we respectively denote $L$ and $T$

$$L(r) = \langle p_x^E(\mathbf{0}) p_x^E(\mathbf{r}e_x) \rangle = \frac{\sin(R) - R \cos(R)}{R^3} \\ T(r) = \langle p_x^E(\mathbf{0}) p_x^E(\mathbf{r}e_y) \rangle = \frac{R \cos(R) - (1 - R^2) \sin(R)}{2R^3},$$

where $R = kr$. Further useful strategies and elementary correlation functions are provided in the SI.

### 3. Results and discussion

All analytical derivations can be found in great detail in the SI. We mostly state final results here, except when intermediate steps are useful for understanding how a result comes about.

#### 3.1. Magnitudes and relative directions of the optical momenta

We begin by deriving the fundamental statistical distributions for the magnitudes of the Poynting, orbital and spin momenta. In terms of real and imaginary field components, the Poynting momentum writes

$$\mathbf{P} = \frac{1}{2} [\mathbf{p}^E \times \mathbf{p}^H + \mathbf{q}^E \times \mathbf{q}^H].$$

Each component of $\mathbf{P}$ is a sum of products of two Gaussian random variables. As detailed in the SI, isotropy allows us to retrieve the magnitude distribution $D(P)$ from that of the

$x$-component $D_x(P_x)$ only [26]. We briefly outline this first derivation, to see the main steps involved:

$$D_x(P_x) = \left\langle \delta \left( P_x - \sum_{j,k} \frac{\epsilon_{xjk}}{2} [p_j^E p_k^H + q_j^E q_k^H] \right) \right\rangle \\ = \int \frac{ds}{2\pi} e^{-isP_x} \left\langle \exp \left\{ is \frac{1}{2} p_y^E p_z^H \right\} \right\rangle^4 \\ = \int \frac{ds}{2\pi} e^{-isP_x} \left[ \frac{1}{1 + s^2 \sigma_x^4 / 4} \right]^2 \\ = \frac{1 + 2|P_x|}{2\sigma_x^2} \exp \left\{ -\frac{2|P_x|}{\sigma_x^2} \right\},$$

where $\sigma_x^2 = \langle (p_x^E)^2 \rangle = 1/3$ (see tabulated variances in the SI). The second step involves factorization of the average using the statistical independence of field components, the third step uses (6), and the last step is an integration in the complex plane. The distribution for the magnitude of the Poynting momentum is then

$$D(P) = -2P \frac{\partial D_x(P_x)}{\partial P_x} \Big|_{P_x=P} = 108P^2 \exp\{-6P\}.$$

The electric-biased orbital and spin momenta read

$$\mathbf{p}_O^E = \frac{1}{2\omega} \mathbf{p}^E \cdot (\nabla) \mathbf{q}^E - \frac{1}{2\omega} \mathbf{q}^E \cdot (\nabla) \mathbf{p}^E \\ \mathbf{p}_S^E = -\frac{1}{2\omega} (\mathbf{p}^E \cdot \nabla) \mathbf{q}^E + \frac{1}{2\omega} (\mathbf{q}^E \cdot \nabla) \mathbf{p}^E.$$

Again, each component is a sum of products of two Gaussian random variables (one field component and one space derivative), and we only have to find the distribution for the $x$-component. For the orbital momentum this is

$$D_x(P_{O,x}^E) = \left\langle \delta \left( P_{O,x}^E - \sum_j \frac{1}{2} [p_j^E \partial_x q_j^E - q_j^E \partial_x p_j^E] \right) \right\rangle \\ = \int \frac{ds}{2\pi} e^{-isP_{O,x}^E} \left\langle \exp \left\{ is \frac{1}{2} p_x^E \partial_x q_x^E \right\} \right\rangle^2 \left\langle \exp \left\{ is \frac{1}{2} p_y^E \partial_x q_y^E \right\} \right\rangle^4 \\ = \dots$$

and the distribution for the magnitude is

$$D(P_O^E) = 180 P_O^E \left[ e^{-6\sqrt{5}P_O^E} - \left[ 1 - \frac{3\sqrt{10}}{2} P_O^E \right] e^{-3\sqrt{10}P_O^E} \right].$$

Surprisingly, we find that in repeating the calculation for the spin momentum, the result is the same. Indeed, the first steps of the derivation read
describe longitudinal and lateral autocorrelations of a given vector component

$$f(r) = \langle v_i(\mathbf{0}) v_i(\mathbf{r}e_i) \rangle, \quad g(r) = \langle v_i(\mathbf{0}) v_i(\mathbf{r}e_j) \rangle \quad (i \neq j).$$

where the separation vector $re_i$ is taken along some chosen direction $i = x, y, z$. If, in addition, the field is solenoidal ($\nabla \cdot \mathbf{v} = 0$), $f$ and $g$ are related, such that the full correlation tensor can be determined from, for example, the longitudinal correlation function $f$ only,

$$\langle v_i(\mathbf{0}) v_j(\mathbf{r}) \rangle = -\frac{r f'(r)}{2} \frac{r_i r_j}{r^2} + \delta_{ij} \left[ f(r) + \frac{r f'(r)}{2} \right]. \quad (7)$$

Since there are no charges in the model field, cycle-averaging Poynting's theorem yields $\nabla \cdot \mathbf{P} = 0$. As the spin momentum itself is the curl of a vector field [10], it is divergenceless $\nabla \cdot \mathbf{P}_S = 0$, and consequently we also have $\nabla \cdot \mathbf{P}_O = 0$. Hence all momenta are isotropic homogeneous solenoidal random fields, to which the above results apply. They also apply to the complex electric and magnetic fields themselves. In our calculations, we will be able to express all correlation functions using the longitudinal and lateral autocorrelation functions of the electric field [29], that we respectively denote $L$ and $T$

$$L(r) = \langle p_x^E(\mathbf{0}) p_x^E(\mathbf{r}e_x) \rangle = \frac{\sin(R) - R \cos(R)}{R^3} \\ T(r) = \langle p_x^E(\mathbf{0}) p_x^E(\mathbf{r}e_y) \rangle = \frac{R \cos(R) - (1 - R^2) \sin(R)}{2R^3},$$

where $R = kr$. Further useful strategies and elementary correlation functions are provided in the SI.

### 3. Results and discussion

All analytical derivations can be found in great detail in the SI. We mostly state final results here, except when intermediate steps are useful for understanding how a result comes about.

#### 3.1. Magnitudes and relative directions of the optical momenta

We begin by deriving the fundamental statistical distributions for the magnitudes of the Poynting, orbital and spin momenta. In terms of real and imaginary field components, the Poynting momentum writes

$$\mathbf{P} = \frac{1}{2} [\mathbf{p}^E \times \mathbf{p}^H + \mathbf{q}^E \times \mathbf{q}^H].$$

Each component of $\mathbf{P}$ is a sum of products of two Gaussian random variables. As detailed in the SI, isotropy allows us to retrieve the magnitude distribution $D(P)$ from that of the

$x$-component $D_x(P_x)$ only [26]. We briefly outline this first derivation, to see the main steps involved:

$$D_x(P_x) = \left\langle \delta \left( P_x - \sum_{j,k} \frac{\epsilon_{xjk}}{2} [p_j^E p_k^H + q_j^E q_k^H] \right) \right\rangle \\ = \int \frac{ds}{2\pi} e^{-isP_x} \left\langle \exp \left\{ is \frac{1}{2} p_y^E p_z^H \right\} \right\rangle^4 \\ = \int \frac{ds}{2\pi} e^{-isP_x} \left[ \frac{1}{1 + s^2 \sigma_x^4 / 4} \right]^2 \\ = \frac{1 + 2|P_x| / \sigma_x^2}{2\sigma_x^2} \exp \left\{ -\frac{2|P_x|}{\sigma_x^2} \right\},$$

where $\sigma_x^2 = \langle (p_x^E)^2 \rangle = 1/3$ (see tabulated variances in the SI). The second step involves factorization of the average using the statistical independence of field components, the third step uses (6), and the last step is an integration in the complex plane. The distribution for the magnitude of the Poynting momentum is then

$$D(P) = -2P \frac{\partial D_x(P_x)}{\partial P_x} \bigg|_{P_x=P} = 108P^2 \exp \{-6P\}.$$

The electric-biased orbital and spin momenta read

$$\mathbf{p}_O^E = \frac{1}{2\omega} \mathbf{p}^E \cdot (\nabla) \mathbf{q}^E - \frac{1}{2\omega} \mathbf{q}^E \cdot (\nabla) \mathbf{p}^E \\ \mathbf{p}_S^E = -\frac{1}{2\omega} (\mathbf{p}^E \cdot \nabla) \mathbf{q}^E + \frac{1}{2\omega} (\mathbf{q}^E \cdot \nabla) \mathbf{p}^E.$$

Again, each component is a sum of products of two Gaussian random variables (one field component and one space derivative), and we only have to find the distribution for the $x$-component. For the orbital momentum this is

$$D_x(P_{O,x}^E) = \left\langle \delta \left( P_{O,x}^E - \sum_j \frac{1}{2} [p_j^E \partial_x q_j^E - q_j^E \partial_x p_j^E] \right) \right\rangle \\ = \int \frac{ds}{2\pi} e^{-isP_{O,x}^E} \left\langle \exp \left\{ is \frac{1}{2} P_x^E \partial_x q_x^E \right\} \right\rangle^2 \left\langle \exp \left\{ is \frac{1}{2} P_y^E \partial_x q_y^E \right\} \right\rangle^4 \\ = \dots$$

and the distribution for the magnitude is

$$D(P_O^E) = 180 P_O^E \left[ e^{-6\sqrt{5}P_O^E} - \left[ 1 - \frac{3\sqrt{10}}{2} P_O^E \right] e^{-3\sqrt{10}P_O^E} \right].$$

Surprisingly, we find that in repeating the calculation for the spin momentum, the result is the same. Indeed, the first steps of the derivation read
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Magnitudes and relative orientations of the Poynting, orbital and spin momenta. (a) Analytical (lines) and numerical (dots) probability distributions for the magnitudes of the Poynting, orbital and spin momenta. (b) Distribution of the angles between momenta, obtained numerically. (c) Illustration of one realization of the random field, showing the three vector fields on the faces of a cubic region of side $λ/2$, and a set of streamlines seeded near the center of the cube. Numerical data in this figure was obtained from $10^5$ realizations of the random field, each containing $N = 10^3$ plane waves.</figcaption></figure>

$$D_x(P_{S,x}^E) = \left\langle \delta \left( P_{S,x}^E - \sum_j \frac{1}{2} \left[ p_j^E \partial_j q_x^E - q_j^E \partial_j p_x^E \right] \right) \right\rangle = \int \frac{ds}{2\pi} e^{-i s P_{S,x}^E} \left\langle \exp \left\{ \frac{i s}{2} P_x^E \partial_x q_x^E \right\} \right\rangle^2 \left\langle \exp \left\{ \frac{i s}{2} P_y^E \partial_y q_x^E \right\} \right\rangle^4$$

and since $\partial_x q_y^E$ and $\partial_y q_x^E$ are both uncorrelated to $p_y^E$ and have the same variance (see tables in the SI), the rest of the calculation is strictly identical to that for the orbital momentum, and we conclude that the orbital and spin momenta obey the exact same magnitude distribution. All these distributions are wavelength-independent, and only scale with the overall intensity in the field.

They are plotted and checked against numerical estimates in figure 1(a). It is interesting to observe that the spin momentum, usually negligibly small in paraxial beams, becomes here equivalent in magnitude to the orbital momentum, responsible for the actual energy flow. The intuitive reason for the different order of magnitude is that in the fully non-paraxial case, there are waves propagating in all directions, such that all space derivatives result in a factor $\sim ik$, whereas transverse gradients are only of order $\sim 1/W$ in the paraxial case. Moving away from paraxiality, part of the linear momentum converts from an orbital to a spin nature, in

a manner strictly similar to how the angular momentum does [30].

To complete the picture of the three momenta at a given point in space, we present in figure 1(b) the distributions for the angles between each pair of momenta. They were obtained numerically, as attempting to compute analytical joint distributions of two momenta hardly leads to tractable expressions. We observe that the angle between $P_O^E$ and $P_S^E$ has a broad distribution roughly centered on $π/2$ (with a slight skew towards larger angles), indicating that they tend to point in perpendicular directions. Since they have comparable magnitudes, the resulting Poynting momentum is generically not closely aligned with either of them. This implies that the streamlines for the three optical currents tend to diverge away from one another.

In figure 1(c), we illustrate one realization of the random field, in a cubic region of side $λ/2$. We plot the three momenta on the sides of the box, and a set of streamlines seeded near the center of the cube. The three vector fields are indeed observed to generically point in different directions, and the streamlines to follow seemingly unrelated paths in space, crossing with angles in agreement with the distributions of figure 1(b). This reinforces the claim that the Poynting and orbital currents generally provide contrasting pictures of EM energy flow, both in terms of magnitude and direction. These observations could prove important for simulating optical forces in complex nanophotonics systems.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Magnitudes and relative orientations of the Poynting, orbital and spin momenta. (a) Analytical (lines) and numerical (dots) probability distributions for the magnitudes of the Poynting, orbital and spin momenta. (b) Distribution of the angles between momenta, obtained numerically. (c) Illustration of one realization of the random field, showing the three vector fields on the faces of a cubic region of side $\lambda/2$, and a set of streamlines seeded near the center of the cube. Numerical data in this figure was obtained from $10^5$ realizations of the random field, each containing $N = 10^3$ plane waves.</figcaption></figure>

$$ 
\begin{aligned}
D_x(P_{S,x}^E) &= \left\langle \delta \left( P_{S,x}^E - \sum_j \frac{1}{2} \left[ p_j^E \partial_j q_x^E - q_j^E \partial_j p_x^E \right] \right) \right\rangle \\
&= \int \frac{ds}{2\pi} e^{-i s P_{S,x}^E} \left\langle \exp \left\{ \frac{i s}{2} P_x^E \partial_x q_x^E \right\} \right\rangle^2 \left\langle \exp \left\{ \frac{i s}{2} P_y^E \partial_y q_x^E \right\} \right\rangle^4
\end{aligned}
$$

and since $\partial_x q_y^E$ and $\partial_y q_x^E$ are both uncorrelated to $p_y^E$ and have the same variance (see tables in the SI), the rest of the calculation is strictly identical to that for the orbital momentum, and we conclude that the orbital and spin momenta obey the exact same magnitude distribution. All these distributions are wavelength-independent, and only scale with the overall intensity in the field.

They are plotted and checked against numerical estimates in figure 1(a). It is interesting to observe that the spin momentum, usually negligibly small in paraxial beams, becomes here equivalent in magnitude to the orbital momentum, responsible for the actual energy flow. The intuitive reason for the different order of magnitude is that in the fully non-paraxial case, there are waves propagating in all directions, such that all space derivatives result in a factor $\sim ik$, whereas transverse gradients are only of order $\sim 1/W$ in the paraxial case. Moving away from paraxiality, part of the linear momentum converts from an orbital to a spin nature, in

a manner strictly similar to how the angular momentum does [30].

To complete the picture of the three momenta at a given point in space, we present in figure 1(b) the distributions for the angles between each pair of momenta. They were obtained numerically, as attempting to compute analytical joint distributions of two momenta hardly leads to tractable expressions. We observe that the angle between $P_O^E$ and $P_S^E$ has a broad distribution roughly centered on $\pi/2$ (with a slight skew towards larger angles), indicating that they tend to point in perpendicular directions. Since they have comparable magnitudes, the resulting Poynting momentum is generically not closely aligned with either of them. This implies that the streamlines for the three optical currents tend to diverge away from one another.

In figure 1(c), we illustrate one realization of the random field, in a cubic region of side $\lambda/2$. We plot the three momenta on the sides of the box, and a set of streamlines seeded near the center of the cube. The three vector fields are indeed observed to generically point in different directions, and the streamlines to follow seemingly unrelated paths in space, crossing with angles in agreement with the distributions of figure 1(b). This reinforces the claim that the Poynting and orbital currents generally provide contrasting pictures of EM energy flow, both in terms of magnitude and direction. These observations could prove important for simulating optical forces in complex nanophotonics systems.
## 3.2. Short-range organization of the currents

**3.2.1. Spatial correlation tensors.** Going beyond their identical magnitude distribution, we find that the orbital and spin momentum vector fields are actually arranged very differently in space. To explore this, we compute two-point spatial correlation tensors for all pairs of components of a given momentum. Each tensor will be of the form in (7), given entirely by the longitudinal autocorrelation function $f(r)$. For the Poynting momentum, this function writes

$$f_P(r) = \langle P_x(\mathbf{0}) P_x(\mathbf{r} \mathbf{e}_x) \rangle = \sum_{j,k,l,m} \left\langle \frac{1}{4} \epsilon_{jkl} \epsilon_{xlm} \left[ p_j^\mathbf{E} p_k^\mathbf{H} + q_j^\mathbf{E} q_k^\mathbf{H} \right] (\mathbf{0}) \left[ p_l^\mathbf{E} p_m^\mathbf{H} + q_l^\mathbf{E} q_m^\mathbf{H} \right] (\mathbf{r} \mathbf{e}_x) \right\rangle.$$

To evaluate these averages we make use of Isserlis' theorem for moments of Gaussian variables [19]. $f_P(r)$ is obtained as

$$f_P(r) = T^2(r) + \frac{(kr)^2}{4} L^2(r),$$

and the correlation tensor is

$$\begin{aligned} \langle P_i(\mathbf{0}) P_j(\mathbf{r}) \rangle &= \frac{r_i r_j}{r^2} \left[ 2R(R^2 - 3) \sin(2R) \right. \\ & \quad \left. + (6R^2 - 3) \cos(2R) + 2R^4 + 3 \right] / 8R^6 \\ & \quad + \delta_{ij} \left[ R(2 - R^2) \sin(2R) + (1 - 2R^2) \cos(2R) - 1 \right] / 4R^6. \end{aligned}$$

The strategy is similar for the orbital and spin momenta. We find

$$\omega^2 f_O(r) = \frac{1}{2} [L'^2(r) - L(r)L''(r)] + [T'^2(r) - T(r)T''(r)],$$

giving the correlation tensor

$$\begin{aligned} \langle P_{O,i}^{\mathbf{E}}(\mathbf{0}) P_{O,j}^{\mathbf{E}}(\mathbf{r}) \rangle &= \frac{r_i r_j}{r^2} \left[ \frac{1}{2}(R^4 - 24R^2 + 72) R \sin(2R) \right. \\ & \quad \left. + 3(R^4 - 10R^2 + 6) \cos(2R) + R^6 - 3R^4 - 6R^2 - 18 \right] / 4R^8 \\ & \quad + \delta_{ij} \left[ - (R^4 - 20R^2 + 54) R \sin(2R) \right. \\ & \quad \left. + (-5R^4 + 46R^2 - 27) \cos(2R) + 3R^4 + 8R^2 + 27 \right] / 8R^8. \end{aligned}$$

And for the spin momentum,

$$\omega^2 f_S(r) = \frac{3}{4} L'^2(r) - \frac{L(r) L''(r)}{2} - 2L'(r) \frac{T(r)}{r}$$

with the correlation tensor

$$\begin{aligned} \langle P_{S,i}^{\mathbf{E}}(\mathbf{0}) P_{S,j}^{\mathbf{E}}(\mathbf{r}) \rangle &= \frac{r_i r_j}{r^2} \left[ 3(R^4 - 14R^2 + 24) R \sin(2R) \right. \\ & \quad \left. + (16R^4 - 69R^2 + 36) \cos(2R) + 2R^4 - 3R^2 - 36 \right] / 8R^8 \\ & \quad + \delta_{ij} \left[ (-3R^4 + 32R^2 - 54) R \sin(2R) \right. \\ & \quad \left. + (-13R^4 + 52R^2 - 27) \cos(2R) - R^4 + 2R^2 + 27 \right] / 8R^8. \end{aligned}$$

For each momentum, the (normalized) autocorrelation function of the x-component for separation vectors $\mathbf{r}$ in the xy-plane is plotted in the top row of figure 2. For the orbital momentum, the degree of correlation is largely positive, and longer-ranged in the longitudinal direction. In sharp contrast, components of the spin momentum tend to change sign periodically, and more strongly so in the lateral directions. These findings hint at qualitatively distinct spatial organizations for the two currents. In the middle row of figure 2, we show 2D streamlines for each momentum in a slice through one realization of the random field, and color them according to the value of the x-component. Zero-crossings of the x-component are shown in black to better distinguish regions of 'upwards' and 'downwards' flow along $\mathbf{x}$. We observe that the orbital current keeps the same direction across relatively broad channels, with a typical size in accordance with the correlation function given above. Such structures are channels of energy flow. Conversely, the spin current changes direction more frequently, particularly along the lateral (y) direction, forming narrow pockets of oppositely directed flow. These two contrasting behaviors can seemingly be traced back to the elementary building block of the non-paraxial field, consisting of two interfering plane waves and studied in [4], in which it was found that $\mathbf{P}_O$ homogeneously points along the bisector, whereas $\mathbf{P}_S$ oscillates in the transverse direction.

Corresponding results for the Poynting current, shown in the third column of figure 2, indicate a less clear-cut behavior, close but not identical to that of the orbital current. At this point, it is enlightening to compare the currents of the vector EM field to the simpler case of a random complex scalar field $\Psi = p^\Psi + iq^\Psi$, defined by dropping the polarization term in brackets in (4). $\Psi$ obeys the Helmholtz equation with wavevector $\mathbf{k}$, and there is a single, divergenceless current $\mathbf{J} = \frac{1}{2\omega} \mathrm{Im}\{\Psi^* \nabla \Psi\}$. Its longitudinal autocorrelation function is given by

$$\omega^2 f_j(r) = \frac{1}{2} [C'^2(r) - C(r) C''(r)],$$

with $C(r) = L(r) + 2T(r) = \sin(kr)/kr$ (we remark the similarity of this expression to that for the orbital momentum), and the correlation tensor is
## 3.2. Short-range organization of the currents

**3.2.1. Spatial correlation tensors.** Going beyond their identical magnitude distribution, we find that the orbital and spin momentum vector fields are actually arranged very differently in space. To explore this, we compute two-point spatial correlation tensors for all pairs of components of a given momentum. Each tensor will be of the form in (7), given entirely by the longitudinal autocorrelation function $f(r)$. For the Poynting momentum, this function writes

$$f_P(r) = \langle P_x(\mathbf{0}) P_x(\mathbf{r}_x) \rangle \\ = \sum_{j,k,l,m} \left\langle \frac{1}{4} \epsilon_{jkl} \epsilon_{xlm} \left[ p_j^E p_k^H + q_j^E q_k^H \right] (\mathbf{0}) \left[ p_l^E p_m^H + q_l^E q_m^H \right] (\mathbf{r}_x) \right\rangle .$$

To evaluate these averages we make use of Isserlis' theorem for moments of Gaussian variables [19]. $f_P(r)$ is obtained as

$$f_P(r) = T^2(r) + \frac{(kr)^2}{4} L^2(r),$$

and the correlation tensor is

$$\langle P_i(\mathbf{0}) P_j(\mathbf{r}) \rangle = \frac{r_i r_j}{r^2} \left[ \begin{aligned} & 2R(R^2 - 3) \sin(2R) \\ & + (6R^2 - 3) \cos(2R) + 2R^4 + 3 \end{aligned} \right] / 8R^6 \\ + \delta_{ij} \left[ R(2 - R^2) \sin(2R) + (1 - 2R^2) \cos(2R) - 1 \right] / 4R^6.$$

The strategy is similar for the orbital and spin momenta. We find

$$\omega^2 f_O(r) = \frac{1}{2} [L'^2(r) - L(r)L''(r)] + [T'^2(r) - T(r)T''(r)],$$

giving the correlation tensor

$$\langle P_{O,i}^E(\mathbf{0}) P_{O,j}^E(\mathbf{r}) \rangle = \frac{r_i r_j}{r^2} \left[ \begin{aligned} & \frac{1}{2}(R^4 - 24R^2 + 72) R \sin(2R) \\ & +3(R^4 - 10R^2 + 6) \cos(2R) + R^6 - 3R^4 - 6R^2 - 18 \end{aligned} \right] / 4R^8 \\ + \delta_{ij} \left[ - (R^4 - 20R^2 + 54) R \sin(2R) \right. \\ \left. + (-5R^4 + 46R^2 - 27) \cos(2R) + 3R^4 + 8R^2 + 27 \right] / 8R^8.$$

And for the spin momentum,

$$\omega^2 f_S(r) = \frac{3}{4} L'^2(r) - \frac{L(r) L''(r)}{2} - 2L'(r) \frac{T(r)}{r}$$

with the correlation tensor

$$\langle P_{S,i}^E(\mathbf{0}) P_{S,j}^E(\mathbf{r}) \rangle = \frac{r_i r_j}{r^2} \left[ \begin{aligned} & 3(R^4 - 14R^2 + 24) R \sin(2R) \\ & + (16R^4 - 69R^2 + 36) \cos(2R) + 2R^4 - 3R^2 - 36 \end{aligned} \right] / 8R^8 \\ + \delta_{ij} \left[ (-3R^4 + 32R^2 - 54) R \sin(2R) \right. \\ \left. + (-13R^4 + 52R^2 - 27) \cos(2R) - R^4 + 2R^2 + 27 \right] / 8R^8.$$

For each momentum, the (normalized) autocorrelation function of the x-component for separation vectors $\mathbf{r}$ in the xy-plane is plotted in the top row of figure 2. For the orbital momentum, the degree of correlation is largely positive, and longer-ranged in the longitudinal direction. In sharp contrast, components of the spin momentum tend to change sign periodically, and more strongly so in the lateral directions. These findings hint at qualitatively distinct spatial organizations for the two currents. In the middle row of figure 2, we show 2D streamlines for each momentum in a slice through one realization of the random field, and color them according to the value of the x-component. Zero-crossings of the x-component are shown in black to better distinguish regions of 'upwards' and 'downwards' flow along $\mathbf{x}$. We observe that the orbital current keeps the same direction across relatively broad channels, with a typical size in accordance with the correlation function given above. Such structures are channels of energy flow. Conversely, the spin current changes direction more frequently, particularly along the lateral (y) direction, forming narrow pockets of oppositely directed flow. These two contrasting behaviors can seemingly be traced back to the elementary building block of the non-paraxial field, consisting of two interfering plane waves and studied in [4], in which it was found that $\mathbf{P}_O$ homogeneously points along the bisector, whereas $\mathbf{P}_S$ oscillates in the transverse direction.

Corresponding results for the Poynting current, shown in the third column of figure 2, indicate a less clear-cut behavior, close but not identical to that of the orbital current. At this point, it is enlightening to compare the currents of the vector EM field to the simpler case of a random complex scalar field $\Psi = p^\Psi + iq^\Psi$, defined by dropping the polarization term in brackets in (4). $\Psi$ obeys the Helmholtz equation with wavevector $\mathbf{k}$, and there is a single, divergenceless current $\mathbf{J} = \frac{1}{2\omega} \mathrm{Im}\{\Psi^* \nabla \Psi\}$. Its longitudinal autocorrelation function is given by

$$\omega^2 f_j(r) = \frac{1}{2} [C'^2(r) - C(r) C''(r)],$$

with $C(r) = L(r) + 2T(r) = \sin(kr)/kr$ (we remark the similarity of this expression to that for the orbital momentum), and the correlation tensor is
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Local structure of the optical currents. First three columns: optical currents of the vector EM field $\mathbf{V} = \mathbf{P}_O^E, \mathbf{P}_S^E, \mathbf{P}$. Last column: current in the complex scalar field $\mathbf{V} = \mathbf{J}$. First row: normalized analytical spatial autocorrelation functions $\langle V_x(\mathbf{0})V_x(\mathbf{r}) \rangle / \langle (V_x)^2 \rangle$ of the x-component of each momentum, for separation vectors $\mathbf{r}$ in the xy-plane. Second row: in-plane streamlines of each momentum in a slice through one realization of the random vector field (first three columns) and one realization of the random scalar field (last column), each containing $N = 10^3$ plane waves. Streamlines are colored according to the value of the x-component of the vector field, and zero-crossings are shown in black to better distinguish regions having a flow oppositely directed along $\mathbf{x}$. Third row: in-plane streamlines in the same slices as in the second row, after local averaging of the vector fields over a spherical volume of diameter $\lambda$ (dashed circle). All plots in a given row share the same colorbar.</figcaption></figure>

$$ \langle J_i(\mathbf{0}) J_j(\mathbf{r}) \rangle = \frac{r_i r_j}{r^2} \left[ \frac{(2R^2 + R \sin(2R) + 2 \cos(2R) - 2)}{4R^4} \right] + \delta_{ij} \left[ - \frac{(R \sin(2R) + \cos(2R) - 1)}{4R^4} \right]. $$

The correlation behavior of the scalar current, shown in the rightmost column of figure 2, appears to lie in between that of the orbital and Poynting currents, and is similar to both. This in turn emphasizes the ‘spin’ nature of $\mathbf{P}_S$, which possesses a behavior unfound in the scalar case; it raises the interesting question of how corresponding currents would behave for tensor waves describing other fundamental particles with different spin.

Finally, we discuss the experimental observability of these optical currents. As mentioned in [11], a small probe particle can hardly image subwavelength structures, as its own presence will distort the field on a comparable lengthscale. With

this in mind, it is tempting to only consider local spatial averages of the currents. Our correlation functions suggest that the orbital current will survive local averaging, as it is largely positively correlated to itself over a wavelength-sized volume. Conversely, neighboring pockets of opposite spin flow will cancel each other out. In the bottom row of figure 2, we plot the same streamlines again, but after having performed a local average of the field over a spherical volume of diameter $\lambda$ (rendered by the dashed circle). The integrated spin current indeed quickly vanishes. As a result, orbital and Poynting currents will tend to reconcile, if probed by sufficiently large particles that effectively average over the generic subwavelength inhomogeneities of the spin momentum. Consequently, we expect the difference in the orbital and Poynting streamlines highlighted in figure 1 to have its significant impact on the motion of very subwavelength objects, such as single atoms or atomic clusters.
<figure><img src="image_6.png" /><figcaption>Figure 2. Local structure of the optical currents. First three columns: optical currents of the vector EM field $V = P_O^E, P_S^E, P$. Last column: current in the complex scalar field $V = J$. First row: normalized analytical spatial autocorrelation functions $\langle V_x(0)V_x(r) \rangle / \langle (V_x)^2 \rangle$ of the x-component of each momentum, for separation vectors $r$ in the xy-plane. Second row: in-plane streamlines of each momentum in a slice through one realization of the random vector field (first three columns) and one realization of the random scalar field (last column), each containing $N = 10^3$ plane waves. Streamlines are colored according to the value of the x-component of the vector field, and zero-crossings are shown in black to better distinguish regions having a flow oppositely directed along $x$. Third row: in-plane streamlines in the same slices as in the second row, after local averaging of the vector fields over a spherical volume of diameter $\lambda$ (dashed circle). All plots in a given row share the same colorbar.</figcaption></figure>

$$ \langle J_i(\mathbf{0}) J_j(\mathbf{r}) \rangle = \frac{r_i r_j}{r^2} \left[ \frac{(2R^2 + R \sin(2R) + 2 \cos(2R) - 2)}{4R^4} \right] + \delta_{ij} \left[ - \frac{(R \sin(2R) + \cos(2R) - 1)}{4R^4} \right]. $$

The correlation behavior of the scalar current, shown in the rightmost column of figure 2, appears to lie in between that of the orbital and Poynting currents, and is similar to both. This in turn emphasizes the ‘spin’ nature of $P_S$, which possesses a behavior unfound in the scalar case; it raises the interesting question of how corresponding currents would behave for tensor waves describing other fundamental particles with different spin.

Finally, we discuss the experimental observability of these optical currents. As mentioned in [11], a small probe particle can hardly image subwavelength structures, as its own presence will distort the field on a comparable lengthscale. With

this in mind, it is tempting to only consider local spatial averages of the currents. Our correlation functions suggest that the orbital current will survive local averaging, as it is largely positively correlated to itself over a wavelength-sized volume. Conversely, neighboring pockets of opposite spin flow will cancel each other out. In the bottom row of figure 2, we plot the same streamlines again, but after having performed a local average of the field over a spherical volume of diameter $\lambda$ (rendered by the dashed circle). The integrated spin current indeed quickly vanishes. As a result, orbital and Poynting currents will tend to reconcile, if probed by sufficiently large particles that effectively average over the generic subwavelength inhomogeneities of the spin momentum. Consequently, we expect the difference in the orbital and Poynting streamlines highlighted in figure 1 to have its significant impact on the motion of very subwavelength objects, such as single atoms or atomic clusters.
<figure><img src="image_6.png" /><figcaption>Figure 3. Vorticity distributions. Analytical (lines) and numerical (dots) probability distributions for the magnitudes of the vorticities of the Poynting, orbital and spin currents. Numerical data was obtained from $10^5$ independent realizations of the random field, each containing $N = 10^3$ plane waves.</figcaption></figure>

### 3.2.2. Vorticity of the currents.
The tendency of the spin current to ‘turn’ more can be further quantified by deriving statistical distributions for the vorticities of the optical currents, that were discussed in previous studies [10, 11]

$$ \Omega_P = \nabla \times \mathbf{P} \qquad \Omega_O^E = \nabla \times \mathbf{P}_O^E \qquad \Omega_S^E = \nabla \times \mathbf{P}_S^E. \qquad (8) $$

The strategy for these calculations follows closely that for the magnitudes of the momenta themselves. We note that the additional space derivative involved now makes the distributions wavelength-dependent. The magnitude distributions for the three vorticities are

$$ D(X = \Omega_P / \omega) = \frac{9X}{78\,886\,240} \times \left[ \begin{aligned} & (237\sqrt{5}X + 172)\,819\,200e^{-6\sqrt{5}X} \\ & + (939\,752\,400X - 44\,642\,639\sqrt{10})\sinh(20X)e^{-8\sqrt{10}X} \\ & + (1\,860\,213\sqrt{10}X - 880\,640)\,160\cosh(20X)e^{-8\sqrt{10}X} \end{aligned} \right] $$

$$ D(X = \Omega_O^E / \omega) = \frac{225}{77} X \times \left[ 64e^{-(15/2)X} - 99e^{-10X} + 35e^{-6\sqrt{5}X} \right] $$

$$ D(X = \Omega_S^E / \omega) = \frac{25X}{361\,504} \times \left[ \begin{aligned} & 83\,187e^{-20X} + 9\,628\,125e^{-12X} \\ & -5\,824\,512e^{-10X} + 286\,374\sqrt{10}\sinh(20X)e^{-8\sqrt{10}X} \\ & -4\,792\,320e^{-6\sqrt{5}X} + 905\,520\cosh(20X)e^{-8\sqrt{10}X} \end{aligned} \right] $$

These distributions are shown in figure 3. Despite being identically distributed in magnitude, orbital and spin momenta have different vorticities: in agreement with the observations of the previous section, that of the spin current is statistically larger. An interesting extension of this investigation could

be to explore whether or not this relates to some difference in the density of singularities in the orbital and spin flows [11]. The geometry of these singularities, in the special case where all components of the complex electric field vanish, was recently studied in [31], where it was found that the orbital momentum always arranges in elongated ‘pseudo vortex lines’ in the vicinity of such zeros. Visual exploration of the random fields (not shown) indicates that such a coiling structure seems to occur frequently near generic zeros of both the orbital and spin momenta.

### 3.3. Democratic momenta and fields with pure helicity

Throughout this work, we have focused on electric-biased momenta. Equivalent statistics would evidently hold for the magnetic-biased quantities, but not for democratic ones. Berry and Shukla recently investigated the difference between biased and democratic quantities in similar statistical calculations [26], and concluded that as a rule of thumb, democratic quantities tend to vary more smoothly and follow narrower distributions. Indeed, including contributions from both the electric and magnetic fields (which are uncorrelated to some extent) effectively suppresses regions of strong interference, similarly to the way vector quantities built from three field components also show less interference detail than corresponding scalar quantities. We derived the magnitude distributions for the democratic momenta (see SI) and present them in figure 4(a). The distribution is still identical for the orbital and spin parts, but is indeed slightly narrower than for the biased momenta (dashed grey line). Interestingly, when computing the angle distributions numerically in 4(b), we find that the angle between $\mathbf{P}_O^{EH}$ and $\mathbf{P}_S^{EH}$ is on average narrower than that between the corresponding biased quantities. As a result, democratic momenta are (slightly) more closely aligned with the Poynting vector than their biased counterparts. Our investigations in randomly polarized fields did not reveal more striking differences between biased and democratic momenta, and we believe all qualitative descriptions given in previous sections to hold for democratic currents as well.

It is enlightening at this point to backtrack on our assumption of randomly polarized plane wave components, to consider instead random fields with pure helicity $\sigma = \pm 1$. This amounts to fixing $\beta_n$ to 0 or $\pi$ in (4), and enforces $\mathbf{H} = -i\sigma\mathbf{E}$ such that biased and democratic quantities become equal (we denote them by a $\sigma$ superscript). As detailed in the SI, this adds new non-zero correlations between variables in our statistics, though values of local averages that were already non-zero in the randomly polarized case are unaffected. Taking these new correlations into account, we can proceed through similar calculations. It is however easy to predict what the distributions will be, as democratic momenta always split into two independent terms originating from components of opposite helicity $\mathbf{P} = [\mathbf{P}^+ + \mathbf{P}^-]/2$ [18]. In a randomly polarized field, this becomes a sum of two independent identically distributed variables, whose distribution simply results from the self-convolution of the distribution for a pure helicity term. This easily appears considering the Fourier transform form of our calculations (see SI). Distributions in pure helicity fields are
<figure><img src="image_7.png" /><figcaption>Figure 3. Vorticity distributions. Analytical (lines) and numerical (dots) probability distributions for the magnitudes of the vorticities of the Poynting, orbital and spin currents. Numerical data was obtained from $10^5$ independent realizations of the random field, each containing $N = 10^3$ plane waves.</figcaption></figure>

### 3.2.2. Vorticity of the currents.
The tendency of the spin current to ‘turn’ more can be further quantified by deriving statistical distributions for the vorticities of the optical currents, that were discussed in previous studies [10, 11]

$$ \Omega_P = \nabla \times \mathbf{P} \qquad \Omega_O^E = \nabla \times \mathbf{P}_O^E \qquad \Omega_S^E = \nabla \times \mathbf{P}_S^E. \qquad (8) $$

The strategy for these calculations follows closely that for the magnitudes of the momenta themselves. We note that the additional space derivative involved now makes the distributions wavelength-dependent. The magnitude distributions for the three vorticities are

$$ \begin{aligned} D(X = \Omega_P / \omega) &= \frac{9X}{78\,886\,240} \times \left[ \begin{aligned} & (237\sqrt{5}X + 172)\,819\,200e^{-6\sqrt{5}X} \\ & + (939\,752\,400X - 44\,642\,639\sqrt{10})\sinh(20X)e^{-8\sqrt{10}X} \\ & + (1\,860\,213\sqrt{10}X - 880\,640)\,160\cosh(20X)e^{-8\sqrt{10}X} \end{aligned} \right] \end{aligned} $$

$$ D(X = \Omega_O^E / \omega) = \frac{225}{77} X \times \left[ 64e^{-(15/2)X} - 99e^{-10X} + 35e^{-6\sqrt{5}X} \right] $$

$$ D(X = \Omega_S^E / \omega) = \frac{25X}{361\,504} \times \left[ \begin{aligned} & 83\,187e^{-20X} + 9\,628\,125e^{-12X} \\ & -5\,824\,512e^{-10X} + 286\,374\sqrt{10}\sinh(20X)e^{-8\sqrt{10}X} \\ & -4\,792\,320e^{-6\sqrt{5}X} + 905\,520\cosh(20X)e^{-8\sqrt{10}X} \end{aligned} \right] $$

These distributions are shown in figure 3. Despite being identically distributed in magnitude, orbital and spin momenta have different vorticities: in agreement with the observations of the previous section, that of the spin current is statistically larger. An interesting extension of this investigation could

be to explore whether or not this relates to some difference in the density of singularities in the orbital and spin flows [11]. The geometry of these singularities, in the special case where all components of the complex electric field vanish, was recently studied in [31], where it was found that the orbital momentum always arranges in elongated ‘pseudo vortex lines’ in the vicinity of such zeros. Visual exploration of the random fields (not shown) indicates that such a coiling structure seems to occur frequently near generic zeros of both the orbital and spin momenta.

### 3.3. Democratic momenta and fields with pure helicity

Throughout this work, we have focused on electric-biased momenta. Equivalent statistics would evidently hold for the magnetic-biased quantities, but not for democratic ones. Berry and Shukla recently investigated the difference between biased and democratic quantities in similar statistical calculations [26], and concluded that as a rule of thumb, democratic quantities tend to vary more smoothly and follow narrower distributions. Indeed, including contributions from both the electric and magnetic fields (which are uncorrelated to some extent) effectively suppresses regions of strong interference, similarly to the way vector quantities built from three field components also show less interference detail than corresponding scalar quantities. We derived the magnitude distributions for the democratic momenta (see SI) and present them in figure 4(a). The distribution is still identical for the orbital and spin parts, but is indeed slightly narrower than for the biased momenta (dashed grey line). Interestingly, when computing the angle distributions numerically in 4(b), we find that the angle between $\mathbf{P}_O^{EH}$ and $\mathbf{P}_S^{EH}$ is on average narrower than that between the corresponding biased quantities. As a result, democratic momenta are (slightly) more closely aligned with the Poynting vector than their biased counterparts. Our investigations in randomly polarized fields did not reveal more striking differences between biased and democratic momenta, and we believe all qualitative descriptions given in previous sections to hold for democratic currents as well.

It is enlightening at this point to backtrack on our assumption of randomly polarized plane wave components, to consider instead random fields with pure helicity $\sigma = \pm 1$. This amounts to fixing $\beta_n$ to 0 or $\pi$ in (4), and enforces $\mathbf{H} = -i\sigma\mathbf{E}$ such that biased and democratic quantities become equal (we denote them by a $\sigma$ superscript). As detailed in the SI, this adds new non-zero correlations between variables in our statistics, though values of local averages that were already non-zero in the randomly polarized case are unaffected. Taking these new correlations into account, we can proceed through similar calculations. It is however easy to predict what the distributions will be, as democratic momenta always split into two independent terms originating from components of opposite helicity $\mathbf{P} = [\mathbf{P}^+ + \mathbf{P}^-]/2$ [18]. In a randomly polarized field, this becomes a sum of two independent identically distributed variables, whose distribution simply results from the self-convolution of the distribution for a pure helicity term. This easily appears considering the Fourier transform form of our calculations (see SI). Distributions in pure helicity fields are
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Statistics of democratic momenta in randomly polarized and pure helicity fields. (a) Analytical (lines) and numerical (dots) probability distributions for the magnitudes of the Poynting, orbital and spin (democratic) momenta in randomly polarized fields (the thin dashed curve shows the distribution for the biased momenta of figure 1), and in pure helicity fields ($P^\sigma$). (b) Distribution of the angles between democratic momenta in randomly polarized fields, obtained numerically. Numerical data in this figure was obtained from $10^5$ realizations of the random field, each containing $N = 10^3$ plane waves.</figcaption></figure>

also shown and checked against numerics in figure 4(a), and they are broader than all distributions in the randomly polarized case. This is likely explained by a weaker ‘suppression of interference’ effect, since there is now even less independence between the different components of the EM field.

Finally, it was recently shown by Aiello that for *instantaneous* (that is, not time-averaged) democratic quantities, the fast-oscillating double-frequency terms also happen to be the cross-helicity terms [32, 33]. Consequently, cycle-averaging becomes equivalent to ignoring cross-helicity terms, and has therefore no effect on democratic quantities in pure helicity fields. For this reason, the distributions derived here for pure helicity fields are also expected to be the magnitude distributions for *instantaneous* democratic momenta (for which the nature of the polarization should be irrelevant). Extending our approach to general time-dependent polychromatic fields is beyond the scope of this article, but represents an intriguing avenue, that could highlight profound relations between electric-magnetic democracy, helicity and time-averaging.

## 4. Concluding remarks

We have investigated various statistical properties of the Poynting, orbital and spin optical momenta in generic isotropic random light fields. This isotropy, equivalent to extreme non-paraxiality, was found to increase the discrepancy between Poynting and orbital flows, as the spin momentum unexpectedly becomes equivalent in magnitude to the orbital one. Deriving correlation functions, we were able to describe the distinct spatial structures of the orbital and spin currents, the former arranging in broad channels of energy flow akin to those found in a scalar random field, when the latter has higher vorticity and changes direction on a subwavelength scale. Upon local averaging over a wavelength-sized volume, the spin current rapidly averages out, leading the orbital and Poynting currents to reconcile. Still, the very different behavior of the orbital and spin currents interrogates what our approach would reveal in other types of waves. Indeed, the field-theoretic formalism decomposing the kinetic (Poynting) momentum into canonical (orbital) and Belinfante (spin) parts is of broader generality, and these investigations could be extended and compared to waves describing other particles, such as electrons described by the Dirac equation whose current decomposition into orbital and spin contributions is known as the Gordon decomposition [34, 35], but also to turbulence in acoustic [36] and gravity water waves [37], the latter extensions appearing very natural considering that results from fluid dynamics were used in the present study. The spin angular momentum density of all types of waves could also be studied, as it is arguably the more relevant quantity from a field-theory perspective, the Belinfante momentum being constructed from it.

Preliminary investigations have been extended to the random paraxial case (equivalent to the models in [24, 25]). In particular, it appears that in analogy with the isotropic case, here the *transverse* orbital and spin momenta have equivalent probability distributions ! The question about correlations appears more subtle, and depends on the chosen paraxial spectrum; upon spatial averaging we anticipate this will lead again to dominance of the orbital current. Furthermore, such correlations between the orbital and spin parts might indicate features of optical spin-orbit relations distinct from those in the ordered, structured paraxial beams previously studied, although such a discussion goes beyond the scope of this work.

Further investigations of the isotropic EM case could characterize the generic singularities of the optical currents (isolated points in 3D space) and the statistical geometry of the flows around them, something that has so far only been explored for non-generic zeros of the full complex electric field [31]. More advanced correlation functions (involving more than two positions, evaluated near extrema, etc...) could reveal finer features of the optical currents as well; random fields generally offer endless possibilities of statistical investigation [21].

Finally, there appears to be profound links to uncover in relating electric-magnetic democracy, helicity and time-averaging. This prompts the extension of our approach to general time-dependent fields, which could require introducing
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Statistics of democratic momenta in randomly polarized and pure helicity fields. (a) Analytical (lines) and numerical (dots) probability distributions for the magnitudes of the Poynting, orbital and spin (democratic) momenta in randomly polarized fields (the thin dashed curve shows the distribution for the biased momenta of figure 1), and in pure helicity fields ($P^\sigma$). (b) Distribution of the angles between democratic momenta in randomly polarized fields, obtained numerically. Numerical data in this figure was obtained from $10^5$ realizations of the random field, each containing $N = 10^3$ plane waves.</figcaption></figure>

also shown and checked against numerics in figure 4(a), and they are broader than all distributions in the randomly polarized case. This is likely explained by a weaker ‘suppression of interference’ effect, since there is now even less independence between the different components of the EM field.

Finally, it was recently shown by Aiello that for *instantaneous* (that is, not time-averaged) democratic quantities, the fast-oscillating double-frequency terms also happen to be the cross-helicity terms [32, 33]. Consequently, cycle-averaging becomes equivalent to ignoring cross-helicity terms, and has therefore no effect on democratic quantities in pure helicity fields. For this reason, the distributions derived here for pure helicity fields are also expected to be the magnitude distributions for *instantaneous* democratic momenta (for which the nature of the polarization should be irrelevant). Extending our approach to general time-dependent polychromatic fields is beyond the scope of this article, but represents an intriguing avenue, that could highlight profound relations between electric-magnetic democracy, helicity and time-averaging.

## 4. Concluding remarks

We have investigated various statistical properties of the Poynting, orbital and spin optical momenta in generic isotropic random light fields. This isotropy, equivalent to extreme non-paraxiality, was found to increase the discrepancy between Poynting and orbital flows, as the spin momentum unexpectedly becomes equivalent in magnitude to the orbital one. Deriving correlation functions, we were able to describe the distinct spatial structures of the orbital and spin currents, the former arranging in broad channels of energy flow akin to those found in a scalar random field, when the latter has higher vorticity and changes direction on a subwavelength scale. Upon local averaging over a wavelength-sized volume, the spin current rapidly averages out, leading the orbital and Poynting currents to reconcile. Still, the very different behavior of the orbital and spin currents interrogates what our approach would reveal in other types of waves. Indeed, the field-theoretic formalism decomposing the kinetic (Poynting) momentum into canonical (orbital) and Belinfante (spin) parts is of broader generality, and these investigations could be extended and compared to waves describing other particles, such as electrons described by the Dirac equation whose current decomposition into orbital and spin contributions is known as the Gordon decomposition [34, 35], but also to turbulence in acoustic [36] and gravity water waves [37], the latter extensions appearing very natural considering that results from fluid dynamics were used in the present study. The spin angular momentum density of all types of waves could also be studied, as it is arguably the more relevant quantity from a field-theory perspective, the Belinfante momentum being constructed from it.

Preliminary investigations have been extended to the random paraxial case (equivalent to the models in [24, 25]). In particular, it appears that in analogy with the isotropic case, here the *transverse* orbital and spin momenta have equivalent probability distributions ! The question about correlations appears more subtle, and depends on the chosen paraxial spectrum; upon spatial averaging we anticipate this will lead again to dominance of the orbital current. Furthermore, such correlations between the orbital and spin parts might indicate features of optical spin-orbit relations distinct from those in the ordered, structured paraxial beams previously studied, although such a discussion goes beyond the scope of this work.

Further investigations of the isotropic EM case could characterize the generic singularities of the optical currents (isolated points in 3D space) and the statistical geometry of the flows around them, something that has so far only been explored for non-generic zeros of the full complex electric field [31]. More advanced correlation functions (involving more than two positions, evaluated near extrema, etc...) could reveal finer features of the optical currents as well; random fields generally offer endless possibilities of statistical investigation [21].

Finally, there appears to be profound links to uncover in relating electric-magnetic democracy, helicity and time-averaging. This prompts the extension of our approach to general time-dependent fields, which could require introducing
the vector potentials for defining instantaneous momenta, and
the weighing of plane wave components by a power spectrum
[22, 25]. This could represent a step towards better understand-
ing of the spin-orbit decomposition of optical momentum,
which as of today remains largely confined to the monochro-
matic case.

Data availability statement

All data that support the findings of this study are included
within the article (and any supplementary files).

Acknowledgments

We are grateful to Luke Hands, Nikitas Papasimakis, Michael
Berry and Konstantin Bliokh for helpful discussions. MRD
acknowledges support from the EPSRC Centre for Doctoral
Training in Topological Design (EP/S02297X/1). The École
Normale Supérieure (Paris) is gratefully acknowledged for
TG's fellowship.

ORCID iDs

Titouan Gadeyne https://orcid.org/0000-0002-7782-5778

Mark R Dennis https://orcid.org/0000-0003-1147-1804

References

[1] Poynting J H 1884 *Phil. Trans. R. Soc.* **175** 343

[2] Born M and Wolf E 1980 *Principles of Optics* 6th edn
(Pergamon)

[3] Zangwill A 2012 *Modern Electrodynamics* (Cambridge
University Press)

[4] Bekshaev A Y, Bliokh K Y and Nori F 2015 *Phys. Rev. X*
**5** 011039

[5] Fedorov F I 1955 *J. Opt.* **15** 014002

[6] Bliokh K Y, Bekshaev A Y and Nori F 2014 *Nat. Commun.*
*5* 3300

[7] Canaguier-Durand A, Cuche A, Genet C and Ebbesen T W
2013 *Phys. Rev. A* **88** 033831

[8] Sukhov S and Dogariu A 2017 *Rep. Prog. Phys.* **80** 112001

[9] Bekshaev A Y and Soskin M S 2007 *Opt. Commun.* **271** 332

[10] Berry M V 2009 *J. Opt. A: Pure Appl. Opt.* **11** 094001
(available at: https://iopscience.iop.org/article/10.1088/
1464-4258/11/9/094001/meta)

[11] Bekshaev A, Bliokh K Y and Soskin M 2011 *J. Opt.* **13** 053001

[12] Jackson J D 1999 *Classical Electrodynamics* 3rd edn (Wiley)

[13] Bliokh K Y, Bekshaev A Y and Nori F 2013 *New J. Phys.*
*15* 033026

[14] Belinfante F 1940 *Physica* **7** 449

[15] Antognozzi M *et al* 2016 *Nat. Phys.* **12** 731

[16] Nieto-Vesperinas M and Xu X 2022 *Light Sci. Appl.* **11** 297

[17] Cameron R P and Barnett S M 2012 *New J. Phys.* **14** 123019

[18] Aiello A and Berry M V 2015 *J. Opt.* **17** 062001

[19] Goodman J W 1985 *Statistical Optics* (Wiley)

[20] Goodman J W 2007 *Speckle Phenomena in Optics: Theory and Applications* (Roberts & Company)

[21] Freund I 1998 *Waves Random Media* **8** 119

[22] Berry M V and Dennis M R 2000 *Proc. R. Soc. A* **456** 2059

[23] Berry M and Dennis M 2001 *Proc. R. Soc. A* **457** 141

[24] Dennis M R 2002 *Opt. Commun.* **213** 201

[25] Dennis M R 2007 *Eur. Phys. J. Spec. Top.* **145** 191

[26] Berry M V and Shukla P 2019 *J. Opt.* **21** 064002

[27] Dennis M R 2003 *J. Phys. A: Math. Gen.* **36** 6611

[28] Batchelor G K 1953 *The Theory of Homogeneous Turbulence*
(Cambridge University Press)

[29] Bourret R C 1960 *Il Nuovo Cimento* **18** 347

[30] Bliokh K Y, Alonso M A, Ostrovskaya E A and Aiello A 2010
*Phys. Rev. A* **82** 063825

[31] Vernon A J, Dennis M R and Rodríguez-Fortuño F J 2023 3d
zeros in electromagnetic fields (arXiv:2301.03540
[physics.optics])

[32] Aiello A 2022 *J. Phys. A: Math. Theor.* **55** 244004

[33] Aiello A 2022 *Phys. Rev. A* **106** 043519

[34] Strange P 1998 *Relativistic Quantum Mechanics: With Applications in Condensed Matter and Atomic Physics*
(Cambridge University Press)

[35] Bliokh K Y, Dennis M R and Nori F 2017 *Phys. Rev. A*
**96** 023622

[36] Burns L, Bliokh K Y, Nori F and Dressel J 2020 *New J. Phys.*
*22* 053050

[37] Bliokh K Y, Punzmann H, Xia H, Nori F and Shats M 2022
*Sci. Adv.* **8** eabm1295
the vector potentials for defining instantaneous momenta, and the weighing of plane wave components by a power spectrum [22, 25]. This could represent a step towards better understanding of the spin-orbit decomposition of optical momentum, which as of today remains largely confined to the monochromatic case.

## Data availability statement

All data that support the findings of this study are included within the article (and any supplementary files).

## Acknowledgments

We are grateful to Luke Hands, Nikitas Papasimakis, Michael Berry and Konstantin Bliokh for helpful discussions. MRD acknowledges support from the EPSRC Centre for Doctoral Training in Topological Design (EP/S02297X/1). The École Normale Supérieure (Paris) is gratefully acknowledged for TG's fellowship.

## ORCID iDs

Titouan Gadeyne https://orcid.org/0000-0002-7782-5778

Mark R Dennis https://orcid.org/0000-0003-1147-1804

## References

[1] Poynting J H 1884 *Phil. Trans. R. Soc.* **175** 343

[2] Born M and Wolf E 1980 *Principles of Optics* 6th edn (Pergamon)

[3] Zangwill A 2012 *Modern Electrodynamics* (Cambridge University Press)

[4] Bekshaev A Y, Bliokh K Y and Nori F 2015 *Phys. Rev. X* **5** 011039

[5] Fedorov F I 1955 *J. Opt.* **15** 014002

[6] Bliokh K Y, Bekshaev A Y and Nori F 2014 *Nat. Commun.* **5** 3300

[7] Canaguier-Durand A, Cuche A, Genet C and Ebbesen T W 2013 *Phys. Rev. A* **88** 033831

[8] Sukhov S and Dogariu A 2017 *Rep. Prog. Phys.* **80** 112001

[9] Bekshaev A Y and Soskin M S 2007 *Opt. Commun.* **271** 332

[10] Berry M V 2009 *J. Opt. A: Pure Appl. Opt.* **11** 094001 (available at: https://iopscience.iop.org/article/10.1088/1464-4258/11/9/094001/meta)

[11] Bekshaev A, Bliokh K Y and Soskin M 2011 *J. Opt.* **13** 053001

[12] Jackson J D 1999 *Classical Electrodynamics* 3rd edn (Wiley)

[13] Bliokh K Y, Bekshaev A Y and Nori F 2013 *New J. Phys.* **15** 033026

[14] Belinfante F 1940 *Physica* **7** 449

[15] Antognozzi M *et al* 2016 *Nat. Phys.* **12** 731

[16] Nieto-Vesperinas M and Xu X 2022 *Light Sci. Appl.* **11** 297

[17] Cameron R P and Barnett S M 2012 *New J. Phys.* **14** 123019

[18] Aiello A and Berry M V 2015 *J. Opt.* **17** 062001

[19] Goodman J W 1985 *Statistical Optics* (Wiley)

[20] Goodman J W 2007 *Speckle Phenomena in Optics: Theory and Applications* (Roberts & Company)

[21] Freund I 1998 *Waves Random Media* **8** 119

[22] Berry M V and Dennis M R 2000 *Proc. R. Soc. A* **456** 2059

[23] Berry M and Dennis M 2001 *Proc. R. Soc. A* **457** 141

[24] Dennis M R 2002 *Opt. Commun.* **213** 201

[25] Dennis M R 2007 *Eur. Phys. J. Spec. Top.* **145** 191

[26] Berry M V and Shukla P 2019 *J. Opt.* **21** 064002

[27] Dennis M R 2003 *J. Phys. A: Math. Gen.* **36** 6611

[28] Batchelor G K 1953 *The Theory of Homogeneous Turbulence* (Cambridge University Press)

[29] Bourret R C 1960 *Il Nuovo Cimento* **18** 347

[30] Bliokh K Y, Alonso M A, Ostrovskaya E A and Aiello A 2010 *Phys. Rev. A* **82** 063825

[31] Vernon A J, Dennis M R and Rodríguez-Fortuño F J 2023 3d zeros in electromagnetic fields (arXiv:2301.03540 [physics.optics])

[32] Aiello A 2022 *J. Phys. A: Math. Theor.* **55** 244004

[33] Aiello A 2022 *Phys. Rev. A* **106** 043519

[34] Strange P 1998 *Relativistic Quantum Mechanics: With Applications in Condensed Matter and Atomic Physics* (Cambridge University Press)

[35] Bliokh K Y, Dennis M R and Nori F 2017 *Phys. Rev. A* **96** 023622

[36] Burns L, Bliokh K Y, Nori F and Dressel J 2020 *New J. Phys.* **22** 053050

[37] Bliokh K Y, Punzmann H, Xia H, Nori F and Shats M 2022 *Sci. Adv.* **8** eabm1295