Optical momentum distributions in monochromatic, isotropic random vector fields

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.


I. 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 = E × 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 timeaveraged linear momentum P is the sum of orbital P O and spin P S parts, respectively generating the orbital and spin angular momenta [9][10][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 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 rea-sons.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 polarisation 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 behaviour 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 λ.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 organise 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.We work in units where ε 0 = µ 0 = c = 1.In a monochromatic field with frequency ω = ck = c(2π/λ), the electric and magnetic fields are represented by complex-valued vector fields E and H, where the physical fields are E = Re{Ee −iωt } and H = Re{He −iωt }.The temporal cycle-averaged Poynting momentum can be written Using Maxwell's equations and the vector identity (where we use the customary notation [A • (∇)B] i = j A j ∂ i B j ), the Poynting momentum can be split into a sum of orbital and spin momenta [10], ( This splitting is not unique : expressing the Poynting momentum using the electric (resp.magnetic) field only, we obtain the electric-biased (resp.magneticbiased ) momenta P E O,S (resp.P H O,S ) 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, producing the so-called democratic (or dual ) momenta ].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 circularlypolarized plane wave components of same handedness : such fields satisfy E = ±iH, 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 P E O only.We therefore choose to center our discussion on electric-biased quantities, and devote section III C 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 E ≈ e ikz (Ex, Ey, 0) with transverse amplitude and polarization profiles Ex(x, y) and Ey(x, y) varying over a lengthscale W ≫ λ, on the order of the beam waist.We find that P and P E O are mostly longitudinal, that P E S is purely transverse, and the following orders of magnitude We conclude that in a regular optical beam, orbital and Poynting momenta are closely aligned, the spin momentum being small in comparison.

B. Gaussian random optical fields
We model generic, natural EM light fields as superpositions of N → +∞ 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][20][21] and statistical properties of spatial structures in random light fields [22][23][24][25][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, where the sum runs over the N ≫ 1 plane waves, with wavevectors kn sampled uniformly on the sphere of directions with spherical angles (θn, ϕn) and identical magnitudes k, and polarizations sampled uniformly on the Poincaré sphere with angles (βn, αn), and uniformly sampled global phases ψn. e ± (k) = [e 1 ± ie 2 ]/ √ 2 are helicity basis vectors, with {e 1 , e 2 } a basis of two real orthogonal unit vectors transverse to 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]  (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 → +∞, 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 u = (p E x , p E y . ..) where Σ is the covariance matrix, with Σ ij = u i u j .Useful formulae and strategies for computing averages are further described in the SI, and can be found in references [22,23,[25][26][27].

C. 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 v to be isotropic requires the two-point correlation tensor to have the form where f and g are scalar functions depending only on the magnitude r = |r| of the separation vector.They respectively describe longitudinal and lateral autocorrelations of a given vector component where the separation vector re i is taken along some chosen direction i = x, y, z.If, in addition, the field is solenoidal (∇ • 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, Since there are no charges in the model field, cycleaveraging Poynting's theorem yields ∇ • P = 0.As the spin momentum itself is the curl of a vector field [10], it is divergenceless ∇ • P S = 0, and consequently we also have ∇ • 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 where R = kr.Further useful strategies and elementary correlation functions are provided in the SI.

III. 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.

A. 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 Each component of 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 Dx(Px) only [26].We briefly outline this first derivation, to see the main steps involved: where σ 2 x = (p E x ) 2 = 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 The electric-biased orbital and spin momenta read 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 and the distribution for the magnitude is 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 and since ∂xq E y and ∂yq E x are both uncorrelated to p E y 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 wavelengthindependent, 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 ∼ ik, whereas transverse gradients are only of order ∼ 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 E O and P E S 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 realisation 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.

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 To evaluate these averages we make use of Isserlis' theorem for moments of Gaussian variables [19].f P (r) is obtained as and the correlation tensor is The strategy is similar for the orbital and spin momenta.We find giving the correlation tensor And for the spin momentum, with the correlation tensor For each momentum, the (normalized) autocorrelation function of the x-component for separation vectors 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 realisation of the random field, and colour 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 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 behaviours 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 P O homogeneously points along the bisector, whereas 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 behaviour, 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 Ψ = p Ψ +iq Ψ , defined by dropping the polarization term in brackets in (4).Ψ obeys the Helmholtz equation with wavevector k, and there is a single, divergenceless current J = 1 2ω Im{Ψ * ∇Ψ}.Its longitudinal autocorrelation function is given by 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 The correlation behaviour 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 behaviour 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, neighbouring 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 λ (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.

Vorticity of the currents
The tendency of the spin current to "turn" more can be further quantified by deriving statistical distribu- tions for the vorticities of the optical currents, that were discussed in previous studies [10,11] 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  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.

C. Democratic momenta and fields with pure helicity
Throughout this work, we have focused on electricbiased 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 Figure 4.b), we find that the angle between P EH O and P EH S 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 σ = ±1.This amounts to fixing βn to 0 or π in (4), and enforces H = −iσE such that biased and democratic quantities become equal (we denote them by a σ 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 P = [P + + 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 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.

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 behaviour of the orbital and spin currents interrogates what our approach would reveal in other types of waves.Indeed, the fieldtheoretic 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 electromagnetic 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 nongeneric 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 understanding of the spin-orbit decomposition of optical momentum, which as of today remains largely confined to the monochromatic case.
In this supplementary document, we detail the following aspects of the main text: • Parameterization of the random field.Explicit formulae for computing elementary averages are given and an alternative parameterization is discussed.
• Gaussian statistics.Elementary variances and correlations are tabulated for fields with random polarization or pure helicity, and computation strategies are outlined.
• Derivations for all analytical results of the main text.We hope the amount of detail provided can be helpful to students interested in performing similar calculations.

I. PARAMETERIZATION OF THE RANDOM FIELD
In this section we give explicit expressions for the random field, required for computing elementary ensemble averages.We briefly compare the intuitive parameterization given in the main text to the one used in the recent paper by Berry and Shukla [1], which has practical advantages.
In the main text we give the following parameterization of the random electromagnetic field where the sum runs over the N plane wave components.The chosen normalization fixes the average electric intensity Re |E| 2 /2 = 1.k n is the wavevector, parameterized by polar and azimuthal angles (θ n , ϕ n ) on the sphere of directions.In the laboratory Cartesian basis {e x , e y , e z } it writes ψ n is the global phase of the plane wave.α n and β n are respectively the azimuthal and polar angles on the usual Poincaré sphere of polarization, with poles corresponding to circularly polarized states represented by the complex helicity basis vectors e σ (k n ) transverse to k n (with σ = ±1).There are different ways to parameterize this basis, one simple possibility being Circularly-polarized plane waves are eigenfunctions of the curl operator, ∇ × (e ik•r e σ ) = σk(e ik•r e σ ) and we have for a circularly polarized plane wave component H σ (k) = −iσE σ (k), which can be used to check that the form (1) is consistent with Maxwell's equations.
In their recent study [1], Berry and Shukla used an alternative parameterization of the same field, using orthogonal real basis vectors e 1 and e 2 e σ = −iσe iσς e 1 + iσe 2 √ 2 , e 1 (ς) = sin ςe θ − cos ςe ϕ , e 2 (ς) = cos ςe θ + sin ςe ϕ Varying the parameter ς simply amounts to rotating the two real vectors around the axis of propagation.They express the random field in the form which also encompasses all polarizations states and phases.An advantage of using this parameterization is that the real and imaginary parts (p and q) of the complex fields have simpler expressions than when using complex basis vectors.This can be helpful for checking integrals by hand, or when using a symbolic calculation software.It should however be made clear that the angles γ n and µ n are no longer the usual angles on the Poincaré sphere.
In fact, they correspond to angles on a rotated sphere : it can be checked that µ n is a polar angle defined with respect to the diameter having the orthogonal linear polarization states e 1 and e 2 at its ends, and 2γ n is an azimuthal angle around this axis.The remaining degree of freedom is spanned by ς n , which affects both the phase and the orientation of the polarization ellipse.The physical meaning of these angles is arguably less intuitive than with the other parameterization, for which α n rotates the polarization ellipse (thus participating in implementing isotropy), ψ n is only a phase factor and β n alone controls the ellipticity.Nevertheless, ensemble averaging can equivalently be defined as follows The corresponding definition in [1] (Eq. 6.8) appears to be missing a factor 2, leading to the peculiar normalization ⟨1⟩ = 1/2.This likely is a typographical error, since all averages given in the cited paper are actually consistent with the formula above and with numerical simulations.

II. GAUSSIAN STATISTICS A. Table of local variances and correlations
As explained in the main text, all components of the real and imaginary fields (as well as their derivatives) are centered Gaussian random variables.Ensemble averages become integrals of the form over the vector u of all Gaussian variables involved in the quantity •.All we require is Σ, the symmetric covariance matrix, with Σ ij = ⟨u i u j ⟩.For local averages, all u i are evaluated at the same point in space (which can be set to r = 0 owing to stationarity).For the isotropic, randomly polarized field, all variances and non-zero correlations involving field components and first spatial derivatives are listed in Table I.They are consistent with those given by Berry and Shukla [1], however we do include a correlation between derivatives of the electric field that was not mentioned by them.
Variances and correlations involving field components and their first derivatives, for an isotropic, randomly polarized monochromatic field with frequency ω.Averages involving spatial derivatives are rescaled by 1/ω to produce frequencyindependent values.Other averages of pairs of variables are zero (in the randomly polarized case).
These values are easily obtained by direct integration over random angles in the plane wave representation of (1), and it is enlightening to actually proceed through some of them by hand, to understand what these values become in isotropic fields with pure helicity.All variables can be put in the form Re e iψn U n and averaging over the global phases ψ n has the effect of removing cross-terms between different plane waves where we dropped the index since we only have to average angles for a single plane wave from now on.U is always proportional to [e iα/2 cos(β/2)e u ± e −iα/2 sin(β/2)e * u ] with an appropriate complex scalar e u , such that The first situation arises in all averages in Table I.It is clear that the value of these averages is actually independent of β due to the cancellation, and will therefore remains unchanged in the particular case of a field with pure helicity.
The second situation clearly depends on whether β is integrated over (giving zero) or fixed to a single value (giving ±1 for pure helicity fields).In the latter case, new non-zero correlations therefore appear, which can be constructed by considering the replacements p E = ±q H and q E = ∓p H .These should then be considered in statistical calculations for pure helicity fields.

B. Computing statistical distributions
We now review strategies for computing statistical distributions.The probability distribution of a quantity A(u) is expressed as where we used a Fourier representation of the delta.The average ⟨exp{isA(u)}⟩ (known as the characteristic function of the variable A) is easily evaluated if A is quadratic in the variables u, in which cases the following Gaussian integration formula is used it should be noted that, since many of our random variables are uncorrelated and obey the same statistics, the average ⟨exp{isA(u)}⟩ will oftentimes be factorizable as a product of averages over smaller sets of random variables.
The remaining integral over s can be evaluated by standard contour integration techniques (Jordan's lemma, residue theorem).a. Exploiting isotropy As explained in [1], the probability distribution D(|V|) of the norm of a three-dimensional, isotropically-distributed vector can be obtained by the sole knowledge of the distribution of one of its components D i (V i ), using the tomography formula Since components of the momenta are already quadratic in the field variables, their magnitudes involve products of four variables, leading to non-Gaussian integrals.Therefore, it is very advantageous to exploit isotropy by computing D i first, which will involve Gaussian integrals only.Isotropy also enforces that D i (V i ) be even, which allows to perform the contour integrations only for the case V i ≥ 0.

C. Computing spatial correlation functions
We now consider averages involving variables evaluated at two different points in space.We tabulate elementary spatial correlation functions.In the limit r → 0, these functions can be used to recover the averages of Table I.
a. Stationarity An important property of our statistics is stationarity: averages of two quantities are independent of their absolute positions, and only depend on the separation vector r between the two Correlators involving derivatives Making use of the fact that the order of derivation and averaging can be interchanged, correlators involving space derivatives can expressed as derivatives of known correlation functions [2,3].Writing G(r = r 2 − r 1 ) = ⟨u(r 1 )v(r 2 )⟩, we have An important observation is the introduction of a minus sign for derivatives of quantities at the reference position r 1 .
c. Two-point correlation functions of EM field components The real and imaginary parts of both the electric and magnetic fields are isotropic and divergenceless.The full correlation tensor for the three components of a given field is entirely given by the longitudinal or lateral correlation functions, that we respectively denote L and T .They are and can be obtained by direct integration of the plane-wave representation of the variables, or from known results of isotropic turbulence theory (e.g.3.4.16 in [2]).These results evidently hold for q E , p H and q H . Due to the divergenceless nature of these fields, we have We also note the following useful relations We now turn to correlations between components of different fields.Direct integration reveals that the following correlations vanish: However, p E and q H are correlated, due to Maxwell's equations.Direct integration over angles gives This correlator changes sign when swapping i and j, p and q, or E and H.
We note that for the random complex scalar field Ψ, the elementary correlation function is Finally, we tabulate some correlators involving derivatives of field components, that will be encountered in our calculations.The first one is The case of interest will be R = (r, 0, 0), simplifying to We turn to the following average, involving two derivatives.We have, exploiting the previous result Evaluating at R = (r, 0, 0), we find III. DETAILED DERIVATIONS A. Magnitude distributions for the momenta

Poynting momentum
The time-averaged Poynting momentum writes Exploiting isotropy, we will first compute the distribution for the x-component x where σ 2 x = (p E x ) 2 = 1/3 (see tabulated variances).The second step involves factorization of the average using independence of field components (indeed, none of the non-zero correlations of Table I are involved here), the third step uses (7), and the last step involves integration in the complex plane, using the residue theorem at the second-order poles s = ±2i/σ 2 x (for many of the derivations to come, residues were obtained using a symbolic calculation software, by finding the relevant term in the Laurent series).The distribution of the norm of the Poynting vector is obtained by use of (8), yielding This distribution is independent of k, the magnitude of the momentum only scales with the overall intensity of the field.
Orbital and spin momenta (biased) The canonical and spin momenta read It is advantageous at this point to absorb the factor 1/ω in the space derivatives, because the rescaled variances We start by computing the distribution for the x-component of the orbital momentum where A = σ 2 x σ 2 xx /4 = 1/180 and B = σ 2 x σ 2 xy /4 = 1/90.The distribution for the norm is again obtained using ( 8) Perhaps surprisingly, we find that in calculating the distribution for the x-component of the spin momentum, the averages to compute are the same.Indeed, we have and since ∂ x q E y and ∂ y q E x are both uncorrelated to p E y and have the same variance (see Table I), the rest of the calculation is identical, and we find that the orbital and spin momenta have the exact same magnitude distribution.We also find that in our random field, their distribution is independent of k, like that of the Poynting momentum.

Orbital and spin momenta (democratic)
We repeat the previous calculation, using the democratic definition of the momenta The distribution of the x-component is with C = ±1/6 (see Table I) capturing a correlation between electric and magnetic variables.We find, using numerical values to reduce the length of the formula For the same reasons as before, we find very early in the derivation that the distribution of the spin momentum will again be identical to that of the orbital one.
Poynting momentum (pure helicity) Considering now a random field with pure helicity, we have E = ±iH, and the time-averaged Poynting momentum simplifies The derivation becomes and the magnitude distribution is Note how this result could have been predicted : we know P naturally splits into independent helicity components [P + +P − ]/2, that are identically distributed in the randomly polarized field.And indeed, the integrand in the Fourier transform in the second step of the derivation above is simply the square root of that for the democratic momenta in the randomly polarized case (up to the factor 2 of course), such that the latter is simply a self-convolution of the pure helicity result.Another interesting observation is that this distribution is also the one followed by the instantaneous Poynting momentum in the randomly polarized case, since we have P σ = P = E × H and none of the new correlations were involved here.
Orbital and spin momenta (pure helicity) In a pure helicity field, the momenta rewrite It can be checked that using this choice of variables, none of the new non-zero correlations will be involved.The derivation becomes and the magnitude distribution is Inspecting the steps of the derivation, it is evident that the self-convolution of this pure helicity result will produce the randomly polarized case (up to the appropriate factor 2).This can again be understood from the splitting As before, an identical derivation exists for the spin momentum.cancelled out.The distribution of the x-component is The magnitude distribution is Spin momentum vorticity (biased) The vorticity reads where in the last step, we used Maxwell's equation ∇ × E = iωH to express second derivatives of the electric variables using first derivatives of magnetic variables.This removes the need to evaluate variances and correlations involving the second derivatives.With some work, the average required for the distribution can be factorized as a product of four simpler averages with R = kr.

Orbital momentum correlation tensor
The same strategy applies.We have Considering the terms that are equal upon exchange of p and q, we have We now use (23).The first term involving local averages vanishes, and from now on we will drop the spatial dependence : it is understood that in each average, the first variable is evaluated at 0 and the second at re x .= j,k From tabulated results the third term is zero Correlators involving derivatives are non-zero only if the corresponding correlator of components is non-zero.We have The correlation tensor is Spin momentum correlation tensor The same strategy again applies.We have f S (r) = P E S,x (0)P E S,x (re x ) = Noticing again that half of the averages are the same upon exchange of p and q ω 2 f S (r) = We use (23), as before the term with local averages vanishes and we are left with Using tabulated results some averages are zero Scalar current correlation tensor The complex scalar random field Ψ can be defined as a sum of scalar waves as follows Ψ obeys Helmholtz's equation with wavevector k.The current in this case is and it is divergenceless [4], such that the previous strategy again applies.We have using the same procedure, we find arXiv:2309.07807v3 [physics.optics]13 Feb 2024 II.THEORETICAL FRAMEWORK A. Poynting, orbital and spin momenta

FIG. 1 .
FIG. 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 realisation 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.

FIG. 2 .
FIG. 2. Local structure of the optical currents.First three columns : optical currents of the vector EM field V = P E O , P E S , P.Last column : current in the complex scalar field V = J.First row : normalized analytical spatial autocorrelation functions ⟨Vx(0)Vx(r)⟩ / (Vx) 2 of the x-component of each momentum, for separation vectors r in the x − y 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 λ (dashed circle).All plots in a given row share the same colorbar.

FIG. 3 .
FIG. 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.

FIG. 4 .
FIG. 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 σ ).(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.