A non-separability measure for spatially disjoint vectorial fields

Complex vectorial fields are general states of structured light [1, 2], non-separable in their spatial and polarisation degrees of freedom, giving rise to exotic inhomogeneous transverse polarisation distributions [3]. These vectorial forms of structured light have shown their potential in a wide variety of applications [4–8], including optical communications, optical tweezers, optical metrology, among others. Noteworthy, their high resemblance to quantum-entangled states, by virtue of the non-separability of their component degrees of freedom, has been the subject of several studies, which have enabled simulating quantum phenomena in the classical regime. Crucially, as has been stated in different works, such analogy cannot be pushed to the quantum limit of representing nonlocal entanglement [9]. Nonetheless, it can be used as a powerful tool to understand some of the abstract concepts of quantum mechanics. The idea then of classical entanglement is not to replace or simulate quantum entangled systems with classical states, but rather to exploit the potentialities of their analogy with quantum entanglement in applications in the classical word [10]. One example is the demonstration of a 'nonlocal' teleportation with classical systems, which opens the possibility to create a robust hybrid quantum–classical communication system [11]. Another one is the demonstration that classically-entangled photons propagating in atmospheric turbulence, where only one degree of freedom is affected (the spatial), evolve in a similar way to quantum-entangled photons, when only one of the photons is sent through the atmospheric channel. Such demonstration offers the possibility of quantum error correction in free-space or fibre-based optical communication systems [12]. A more complete analysis of such analogies, as well as other examples can be found in [13–17].

Given the high interest in vector beams, in the last two decades there has been intense interest in their generation, with techniques exploiting liquid crystal wave plates [18], glass cones [19, 20], metamaterials [21], interferometric arrays [22–25], directly from lasers [26, 27], liquid crystal spatial light modulators [28–35] and digital micromirror devices [36–43]. Conversely, characterisation techniques are still somewhat limited. These include filtering [44–46] and deterministic detection [47] of the mode, as well as quantitative analysis by a quantum toolkit, for a full decomposition via tomographic projections [48] and a reduced measurement set for determining the non-separability [49–52]. Both quantitative techniques rely on the similarities between the vector modes and quantum-entangled states, for example, the latter being the classical equivalent to the well-known concurrence, C, for qubits [53, 54]. In the classical case, rather than measuring the degree of entanglement of two entangled photons, the measurement has been adapted to measure the degree of coupling between the spatial and polarisation degrees of freedom of vector modes, and has been coined the vector quality factor (VQF) [50]. Notably, such techniques have been extended to the space-time domain, thus providing a technique to measure the nonseparability of space-time electromagnetic pulses [55].

Even though these techniques perform very well at characterising some features of vector beams, such as their non-separability, they fail to account for other properties manifested by more intricate vector modes. One such case is when the spatial modes become spatially disjoint, appearing to follow independent paths although remaining as a single coherent field. This situation can arise naturally in vectorial light and has been discussed theoretically [10] and observed experimentally [42], and even engineered as a path degree of freedom in such fields [56–58]. This path-like splitting is analogous to splitting of particle mixtures, illustrated in figure 1. In the case of such splitting phenomena in vectorial light, the mutable degree of non-separability cannot be quantified using conventional measurements such as the VQF, suggesting the need for a more refined toolbox.

Zoom In Zoom Out Reset image size

Here we introduce a novel approach for the characterisation of vectorial light fields that overcomes the aforementioned shortcoming, namely, it allows to quantify the dynamic evolution of vector fields that upon propagation experience a spatial separation. This is, to the best of our knowledge, the first proposed measure that attempts to quantify the evolution of the polarisation structure in such kind of vector beams. We outline the problem theoretically and propose an amendment to the existing framework to account not only for the non-separability of the entire field, but also its component disjointness. For this, we employ the Hellinger distance [59] as a measure of the disjointness. We provide experimental evidence to verify the efficacy of this technique, using exemplar cases of Hermite- Airy-, and parabolic-Gauss vectorial beams.

In this section, we provide the theoretical background for describing and quantifying the splitting of the polarisation pattern of a non-uniformly polarised light beam into two disjoint uniformly polarised parts, which may occur during free propagation, as sketched conceptually in figure 1. After shortly reviewing in section 2.1 the basic properties of non-uniformly polarised beams of light, we give in section 2.2 a description of the above splitting phenomenon in terms of the degree of polarisation of each part of the beam. Then, in section 2.3 we introduce a propagation-dependent measure of this phenomenon, which is suitable for experimental verification. Finally, we illustrate our approach by means of a real-world example in section 2.4.

2.1. Entangled polarisation and spatial degrees of freedom in optical beams: a cursory look

Consider a monochromatic, non-uniformly polarised paraxial beam of light travelling in vacuum (or air) along the positive direction of the axis z of a Cartesian reference frame with coordinates r = (x, y, z) = ( x , z), where x = (x, y). The measurable real-valued electric field vector E( r , t) carried by such beam can be written as

$$\begin{equation}\mathbf{E}(\boldsymbol{r},t)=\mathrm{R}\mathrm{e}\left[\mathbf{\Psi }(\boldsymbol{x},z)\,{\text{e}}^{\text{i}(kz-\omega t)}\right],\end{equation} \tag{ 1 }$$

where k = 2π/λ is the wavenumber of the light of wavelength λ and frequency ω = kc, with c the speed of light in vacuum. By definition, the vector field Ψ( x , z) can always be written as

$$\begin{equation}\mathbf{\Psi }(\boldsymbol{x},z)={\hat{\boldsymbol{e}}}_{1}\,{\psi }_{1}(\boldsymbol{x},z)+{\hat{\boldsymbol{e}}}_{2}\,{\psi }_{2}(\boldsymbol{x},z),\end{equation} \tag{ 2 }$$

where ${\hat{\boldsymbol{e}}}_{1},{\hat{\boldsymbol{e}}}_{2}$ and ${\hat{\boldsymbol{e}}}_{3}$ are real-valued unit vectors parallel to the axes x, y and z, respectively, with

$$\begin{equation}{\hat{\boldsymbol{e}}}_{i}\cdot {\hat{\boldsymbol{e}}}_{j}={\delta }_{ij},\quad (i,j=1,2,3),\end{equation} \tag{ 3 }$$

and ψ₁( x , z), ψ₂( x , z) are two arbitrary solutions of the paraxial wave equation,

$$\begin{equation}\left(\frac{{\partial }^{2}}{\partial {x}^{2}}+\frac{{\partial }^{2}}{\partial {y}^{2}}+2\mathrm{i}k\frac{\partial }{\partial z}\right){\psi }_{\alpha }(\boldsymbol{x},z)=0,\quad (\alpha =1,2).\end{equation} \tag{ 4 }$$

Alternatively, as shown in appendix C, it is possible to write Ψ( x , z) in the Schmidt form [10],

$$\begin{equation}\mathbf{\Psi }(\boldsymbol{x},z)=\sqrt{{\lambda }_{1}}\,{\hat{\boldsymbol{{\epsilon}}}}_{1}\,{\phi }_{1}(\boldsymbol{x},z)+\sqrt{{\lambda }_{2}}\,{\hat{\boldsymbol{{\epsilon}}}}_{2}\,{\phi }_{2}(\boldsymbol{x},z),\end{equation} \tag{ 5 }$$

where λ₁ ⩾ λ₂ ⩾ 0, are the eigenvalues of the coherency matrix J of the beam (see equation (D4) in appendix D), defined by

$$\begin{equation}J={\int }_{{\mathbb{R}}^{2}}{\mathrm{d}}^{2}x\,\mathbf{\Psi }{\mathbf{\Psi }}^{{\dagger}}.\end{equation} \tag{ 6 }$$

The dyad ΨΨ^† is an x - and z-dependent 2 × 2 matrix which can be written with respect to the basis $\left\{{\hat{\boldsymbol{e}}}_{1},{\hat{\boldsymbol{e}}}_{2}\right\}$ as:

$$\begin{equation}\mathbf{\Psi }{\mathbf{\Psi }}^{{\dagger}}=\left[\begin{matrix}\hfill \vert {\psi }_{1}(\boldsymbol{x},z){\vert }^{2}\hfill & \hfill {\psi }_{1}(\boldsymbol{x},z){\psi }_{2}^{\ast }(\boldsymbol{x},z)\hfill \\ \hfill {\psi }_{1}^{\ast }(\boldsymbol{x},z){\psi }_{2}(\boldsymbol{x},z)\hfill & \hfill \vert {\psi }_{2}(\boldsymbol{x},z){\vert }^{2}\hfill \end{matrix}\right],\end{equation} \tag{ 7 }$$

where equation (2) has been used. We remark that J is independent of z because both ψ₁( x , z) and ψ₂( x , z) are supposedly solutions of the paraxial wave equation. The polarisation modes ${\hat{\boldsymbol{{\epsilon}}}}_{1}$ and ${\hat{\boldsymbol{{\epsilon}}}}_{2}$ are two (possibly complex-valued) unit vectors perpendicular to the z-axis and to each other,

$$\begin{equation}{\hat{\boldsymbol{{\epsilon}}}}_{\alpha }^{\ast }\cdot {\hat{\boldsymbol{e}}}_{3}=0,\qquad {\hat{\boldsymbol{{\epsilon}}}}_{\alpha }^{\ast }\cdot {\hat{\boldsymbol{{\epsilon}}}}_{\beta }={\delta }_{\alpha \beta },\end{equation} \tag{ 8 }$$

with α, β = 1, 2. Likewise, the spatial modes ϕ₁( x , z) and ϕ₂( x , z) are square-integrable solutions of the paraxial wave equation that obey the orthonormality conditions

$$\begin{equation}{\int }_{{\mathbb{R}}^{2}}{\mathrm{d}}^{2}x\,{\phi }_{\alpha }^{\ast }(\boldsymbol{x},z){\phi }_{\beta }(\boldsymbol{x},z)={\delta }_{\alpha \beta },\quad (\alpha ,\beta =1,2).\end{equation} \tag{ 9 }$$

The entanglement between the polarisation and the spatial degrees of freedom of the vector field (2) is usually quantified by the concurrence C defined by

$$\begin{equation}C=2\,\sqrt{{p}_{1}\,{p}_{2}},\quad \text{where}\quad {p}_{\alpha }\equiv \frac{{\lambda }_{\alpha }}{{\lambda }_{1}+{\lambda }_{2}},\end{equation} \tag{ 10 }$$

and α = 1, 2. By definition, the concurrence is related to the degree of polarisation P of the beam [60–62],

$$\begin{equation}P=\sqrt{1-\frac{4\,\mathrm{det}\,J}{{\left(\mathrm{tr}\,J\right)}^{2}}}=\vert {p}_{1}-{p}_{2}\vert ,\end{equation} \tag{ 11 }$$

by the simple relation

$$\begin{equation}{C}^{2}+{P}^{2}=1.\end{equation} \tag{ 12 }$$

When either p₁ = 0 or p₂ = 0, there is not entanglement so that C = 0 and P = 1. Vice versa, when p₁ = p₂ = 1/2 the entanglement is maximal and we obtain C = 1 and P = 0. From the definition (7) of J, it follows that the concurrence C is propagation-invariant, i.e., it does not depend on the propagation distance z, even if the polarisation pattern changes upon propagation.

2.2. Polarisation-pattern splitting and degree of polarisation

When a non-uniformly polarised paraxial beam of light travels in free space, its polarisation pattern typically changes during propagation, although the concurrence and the degree of polarisation do remain constant (see, e.g., [63] and references therein). Of particular interest are those peculiar beams whose polarisation pattern splits upon propagation into two disjoint parts with uniform orthogonal polarisation, as outlined in figure 1 [42]. Let z = z₀ be the distance from the origin of the coordinates at which the split occurs, and let ${\hat{\boldsymbol{{\epsilon}}}}_{1}$ and ${\hat{\boldsymbol{{\epsilon}}}}_{2}$ denote the unit vectors characterising the orthogonal polarisations of the two parts of the beam resulting from the splitting. Note that z₀ may be either finite or infinite as, for example, in the case of the beams investigated in [42]. Then, the field of the beam at z = z₀ must have necessarily the form

$$\begin{equation}\mathbf{\Psi }(\boldsymbol{x},{z}_{0})={\hat{\boldsymbol{{\epsilon}}}}_{1}\,{\varphi }_{1}(\boldsymbol{x},{z}_{0})+{\hat{\boldsymbol{{\epsilon}}}}_{2}\,{\varphi }_{2}(\boldsymbol{x},{z}_{0}),\end{equation} \tag{ 13 }$$

with

$$\begin{equation}{\varphi }_{1}(\boldsymbol{x},{z}_{0})\,{\varphi }_{2}(\boldsymbol{x},{z}_{0})=0,\quad \forall \enspace \boldsymbol{x}\in {\mathbb{R}}^{2}.\end{equation} \tag{ 14 }$$

The last condition is implied by the requirement that at z = z₀ the two polarisation patterns do not overlap at any point in the xy-plane. In more mathematical terms, equation (14) can be written as

$$\begin{equation}{\mathcal{D}}_{1}\cap {\mathcal{D}}_{2}=\varnothing ,\end{equation} \tag{ 15 }$$

where ${\mathcal{D}}_{\alpha }$ denotes the support of the mode function φ_α at z = z₀, namely

$$\begin{equation}{\mathcal{D}}_{\alpha }=supp\left({\varphi }_{\alpha }\right)=\left\{\boldsymbol{x}\in {\mathbb{R}}^{2}:{\varphi }_{\alpha }(\boldsymbol{x},{z}_{0})\ne 0\right\}.\end{equation} \tag{ 16 }$$

Here and hereafter we assume that ${\mathcal{D}}_{1}\cup {\mathcal{D}}_{2}={\mathbb{R}}^{2}$.

From (14) it trivially follows that the spatial modes φ₁( x , z₀) and φ₂( x , z₀) are orthogonal in the sense of (9). Therefore, equation (13) is automatically a Schmidt form of the type (5), with

$$\begin{equation}\sqrt{{\lambda }_{\alpha }}\,{\phi }_{\alpha }(\boldsymbol{x},{z}_{0})={\varphi }_{\alpha }(\boldsymbol{x},{z}_{0}),\end{equation} \tag{ 17 }$$

and

$$\begin{equation}{\lambda }_{\alpha }={\int }_{{\mathbb{R}}^{2}}{\mathrm{d}}^{2}x\,{\left\vert {\varphi }_{\alpha }(\boldsymbol{x},{z}_{0})\right\vert }^{2}.\end{equation} \tag{ 18 }$$

Since paraxial propagation is a unitary process, orthogonality is preserved during propagation, so that

$$\begin{equation}{\int }_{{\mathbb{R}}^{2}}{\mathrm{d}}^{2}x\,{\varphi }_{\alpha }^{\ast }(\boldsymbol{x},z){\varphi }_{\beta }(\boldsymbol{x},z)={\lambda }_{\alpha }\,{\delta }_{\alpha \beta },\quad (\alpha ,\beta =1,2),\end{equation} \tag{ 19 }$$

where (18) has been used, and

$$\begin{equation}{\varphi }_{\alpha }(\boldsymbol{x},z)={\int }_{{\mathbb{R}}^{2}}{\mathrm{d}}^{2}{x}^{\prime }\,U\left(\boldsymbol{x}-{\boldsymbol{x}}^{\prime },z-{z}_{0}\right){\varphi }_{\alpha }({\boldsymbol{x}}^{\prime },{z}_{0}),\end{equation} \tag{ 20 }$$

with

$$\begin{equation}U\left(\boldsymbol{x}-{\boldsymbol{x}}^{\prime },z\right)=\frac{k}{2\pi \mathrm{i}}\frac{1}{z}\,\mathrm{exp}\left[\mathrm{i}\,\frac{k}{2z}{\left(\boldsymbol{x}-{\boldsymbol{x}}^{\prime }\right)}^{2}\right],\end{equation} \tag{ 21 }$$

the Fresnel paraxial propagator [64]. This implies that also the Schmidt decomposition (13) remains valid, with

$$\begin{equation}\mathbf{\Psi }(\boldsymbol{x},z)={\hat{\boldsymbol{{\epsilon}}}}_{1}\,{\varphi }_{1}(\boldsymbol{x},z)+{\hat{\boldsymbol{{\epsilon}}}}_{2}\,{\varphi }_{2}(\boldsymbol{x},z).\end{equation} \tag{ 22 }$$

Another useful decomposition of Ψ( x , z) is the following,

$$\begin{equation}\mathbf{\Psi }(\boldsymbol{x},z)={\mathbf{\Psi }}_{1}(\boldsymbol{x},z)+{\mathbf{\Psi }}_{2}(\boldsymbol{x},z),\end{equation} \tag{ 23 }$$

where

$$\begin{align}\hfill {\mathbf{\Psi }}_{\mu }(\boldsymbol{x},z)& \equiv {\left.\mathbf{\Psi }(\boldsymbol{x},z)\right\vert }_{\boldsymbol{x}\in {\mathcal{D}}_{\mu }}\hfill \\ \hfill & ={\hat{\boldsymbol{{\epsilon}}}}_{1}{\left.{\varphi }_{1}(\boldsymbol{x},z)\right\vert }_{\boldsymbol{x}\in {\mathcal{D}}_{\mu }}+{\hat{\boldsymbol{{\epsilon}}}}_{2}{\left.{\varphi }_{2}(\boldsymbol{x},z)\right\vert }_{\boldsymbol{x}\in {\mathcal{D}}_{\mu }},\hfill \end{align} \tag{ 24 }$$

with μ = 1, 2, denotes the field Ψ( x , z) restricted to the region ${\mathcal{D}}_{\mu }$ of the xy-plane. At z = z₀, equation (24) reduces to

$$\begin{equation}{\mathbf{\Psi }}_{\mu }(\boldsymbol{x},{z}_{0})={\hat{\boldsymbol{{\epsilon}}}}_{\mu }{\varphi }_{\mu }(\boldsymbol{x},{z}_{0}),\quad (\mu =1,2).\end{equation} \tag{ 25 }$$

Both fields Ψ₁( x , z₀) and Ψ₂( x , z₀) defined by (25), are fully polarised so that their degree of polarisation is equal to one. However, for z ≠ z₀, Ψ₁( x , z) and Ψ₂( x , z) are typically partially polarised because of the physical superposition between ${\left.{\varphi }_{1}(\boldsymbol{x},z)\right\vert }_{\boldsymbol{x}\in {\mathcal{D}}_{\mu }}$ and ${\left.{\varphi }_{2}(\boldsymbol{x},z)\right\vert }_{\boldsymbol{x}\in {\mathcal{D}}_{\mu }}$, which may occur because of diffraction. In this case, the degree of polarisation of each field must be less than one. This implies that the degree of polarisation of Ψ_μ ( x , z), denoted P_μ (z), will depend on the propagation distance z and can therefore be used as a measure of the polarisation-pattern splitting phenomenon. P_μ (z) can be calculated either via (11), or as $P=\sqrt{1-{C}^{2}}$, where (12) has been used. Then, using the expression (E3) for C, we can eventually write

$$\begin{equation}{P}_{\mu }(z)=\frac{\sqrt{{S}_{1\mu }^{2}(z)+{S}_{2\mu }^{2}(z)+{S}_{3\mu }^{2}(z)}}{{S}_{0\mu }(z)},\end{equation} \tag{ 26 }$$

where the measurable Stokes parameters of the field Ψ_μ ( x , z), are given by

$$\begin{equation}\begin{aligned}\hfill {S}_{0\mu }(z)& ={\int }_{{\mathcal{D}}_{\mu }}{\mathrm{d}}^{2}x\left({\left\vert {\psi }_{1}\right\vert }^{2}+{\left\vert {\psi }_{2}\right\vert }^{2}\right),\hfill \\ \hfill {S}_{1\mu }(z)& ={\int }_{{\mathcal{D}}_{\mu }}{\mathrm{d}}^{2}x\left({\left\vert {\psi }_{1}\right\vert }^{2}-{\left\vert {\psi }_{2}\right\vert }^{2}\right),\hfill \\ \hfill {S}_{2\mu }(z)& ={\int }_{{\mathcal{D}}_{\mu }}{\mathrm{d}}^{2}x\left({\psi }_{1}^{\ast }{\psi }_{2}+{\psi }_{1}{\psi }_{2}^{\ast }\right),\hfill \\ \hfill {S}_{3\mu }(z)& =\frac{1}{\mathrm{i}}{\int }_{{\mathcal{D}}_{\mu }}{\mathrm{d}}^{2}x\left({\psi }_{1}^{\ast }{\psi }_{2}-{\psi }_{1}{\psi }_{2}^{\ast }\right),\hfill \end{aligned}\end{equation} \tag{ 27 }$$

with ${\psi }_{\alpha }={\psi }_{\alpha }(\boldsymbol{x},z)={\hat{\boldsymbol{e}}}_{\alpha }\cdot \mathbf{\Psi }(\boldsymbol{x},z)$, and α = 1, 2.

By construction, at z = z₀ we must have P₁(z₀) = 1 = P₂(z₀). However, for z ≠ z₀ we expect to find P_μ (z) < 1. An explicit example of application of equation (26) will be presented in section 2.4.

2.3. The Hellinger distance

Let us look again at the field (22), that we rewrite here as

$$\begin{equation}\mathbf{\Psi }(\boldsymbol{x},z)=\sqrt{{\lambda }_{1}}\,{\hat{\boldsymbol{{\epsilon}}}}_{1}\,{\phi }_{1}(\boldsymbol{x},z)+\sqrt{{\lambda }_{2}}\,{\hat{\boldsymbol{{\epsilon}}}}_{2}\,{\phi }_{2}(\boldsymbol{x},z),\end{equation} \tag{ 28 }$$

where, according to equations (17) and (20), we have defined

$$\begin{equation}{\phi }_{\alpha }(\boldsymbol{x},z)=\frac{{\varphi }_{\alpha }(\boldsymbol{x},z)}{\sqrt{{\lambda }_{\alpha }}},\quad (\alpha =1,2),\end{equation} \tag{ 29 }$$

with

$$\begin{equation}{\phi }_{\alpha }(\boldsymbol{x},{z}_{0}){\phi }_{\beta }^{\ast }(\boldsymbol{x},{z}_{0})={\delta }_{\alpha \beta }\vert {\phi }_{\alpha }(\boldsymbol{x},{z}_{0}){\vert }^{2}.\end{equation} \tag{ 30 }$$

Now, consider the non-negative quantity

$$\begin{equation}h[{\phi }_{1},{\phi }_{2}](z)={\int }_{{\mathbb{R}}^{2}}{\mathrm{d}}^{2}x\,\left\vert {\phi }_{1}(\boldsymbol{x},z)\right\vert \left\vert {\phi }_{2}(\boldsymbol{x},z)\right\vert \geqslant 0,\end{equation} \tag{ 31 }$$

which is a functional of ϕ₁ and ϕ₂, and a function of z. The minimum value

$$\begin{equation}\mathrm{min}\left\{h[{\phi }_{1},{\phi }_{2}](z)\right\}=h[{\phi }_{1},{\phi }_{2}]({z}_{0})=0,\end{equation} \tag{ 32 }$$

is achieved at z = z₀ because of equation (30). The maximum value

$$\begin{equation}\mathrm{max}\left\{h[{\phi }_{1},{\phi }_{2}](z)\right\}=h[{\phi }_{1},{\phi }_{1}](z)=1,\end{equation} \tag{ 33 }$$

is obtained for $\left\vert {\phi }_{1}(\boldsymbol{x},z)\right\vert =\left\vert {\phi }_{2}(\boldsymbol{x},z)\right\vert$, because of the normalization condition (19). Next, we note that the non-negative function f_α ( x , z), defined by

$$\begin{equation}{f}_{\alpha }(\boldsymbol{x},z)\equiv {\left\vert {\phi }_{\alpha }(\boldsymbol{x},z)\right\vert }^{2}\geqslant 0,\quad (\alpha =1,2),\end{equation} \tag{ 34 }$$

has all the properties of the probability density function of a continuous real-valued two-dimensional random vector X = (X₁, X₂), with (X₁, X₂) ∼ f_α (x₁, x₂, z). In probability theory and statistics, there are several ways to quantify the distance between two probability distributions, say f(x) and g(x). Here, considering the properties (32) and (33), it is convenient to use the Hellinger distance H[f, g] [59], defined by

$$\begin{align}\hfill {H}^{2}[f,g]& =\frac{1}{2}\int \mathrm{d}x\,{\left[\sqrt{f(x)}-\sqrt{g(x)}\right]}^{2}\hfill \\ \hfill & =1-\int \mathrm{d}x\,\sqrt{f(x)g(x)}.\hfill \end{align} \tag{ 35 }$$

It is not difficult to show that such distance is a proper metric, that is for all non-negative, normalized smooth functions g₁, g₂, g₃, the Hellinger distance H[g₁, g₂] satisfies the following conditions:

• (a)  
  Positivity condition: H[g₁, g₂] ⩾ 0;
• (b)  
  Symmetry property: H[g₁, g₂] = H[g₂, g₁];
• (c)  
  Identity property: H[g₁, g₂] = 0 if and only if g₁ = g₂;
• (d)  
  Triangular inequality: H[g₁, g₃] ⩽ H[g₁, g₂] + H[g₂, g₃].

In our two-dimensional case the Hellinger distance takes the form

$$\begin{align}\hfill H[\vert {\phi }_{1}\vert ,\vert {\phi }_{2}\vert ](z)& =\sqrt{1-{\int }_{{\mathbb{R}}^{2}}\left\vert {\phi }_{1}(\boldsymbol{x},z)\right\vert \left\vert {\phi }_{2}(\boldsymbol{x},z)\right\vert \,{\mathrm{d}}^{2}x}\hfill \\ \hfill & =\sqrt{1-h[{\phi }_{1},{\phi }_{2}](z)}.\hfill \end{align} \tag{ 36 }$$

We remark that this distance does not depend directly on the amount of the entanglement of the beam, that it is independent of λ₁ and λ₂. However, it does crucially depend on the propagation distance z. For illustration, H[|ϕ₁|, |ϕ₂|](z) is calculated explicitly for an exemplary beam in section 2.4. We would also like to remark that our choice of the Hellinger distance to quantify the separation between ϕ₁ and ϕ₂, is not the only possible one. There are several statistical distances available in the literature like, e.g., the Kullback–Leibler divergence and the Jensen–Shannon distance [65]. However, we have chosen the Hellinger distance because it is a metric (conversely, the Kullback–Leibler divergence is not a metric), and it is simple to calculate. Moreover, it is automatically z-dependent.

2.4. Polarisation-pattern splitting phenomenon: a case study

To illustrate the use of equations (26) and (36), we consider here a simple non-uniformly polarised paraxial beam of light, whose polarisation pattern is split into two disjoint uniformly polarised parts at z = z₀ = 0. Imagine to have a beam of light prepared at z = 0 in the Hermite–Gauss mode u₀₁( x , 0) of Rayleigh length z_R and waist w₀ (see appendix A). Then, suppose to place at z = 0 a thin spherical lens of focal length f = z_R, characterised by the transmission function [64],

$$\begin{equation}t(\boldsymbol{x})={\mathrm{e}}^{-\mathrm{i}k\,\frac{{\left\vert \boldsymbol{x}\right\vert }^{2}}{2f}}={\mathrm{e}}^{-\mathrm{i}\,\frac{{\left\vert \boldsymbol{x}\right\vert }^{2}}{{w}_{0}^{2}}}.\end{equation} \tag{ 37 }$$

By construction, the amplitude ϕ( x ) of the beam at z = 0 immediately after the lens is

$$\begin{equation}\phi (\boldsymbol{x})=t(\boldsymbol{x})\,{u}_{01}(\boldsymbol{x},0)=\sqrt{\frac{8}{\pi }}\,\frac{y}{{w}_{0}^{2}}\,{\mathrm{e}}^{-(1+\mathrm{i})\,\frac{{\left\vert \boldsymbol{x}\right\vert }^{2}}{{w}_{0}^{2}}}.\end{equation} \tag{ 38 }$$

Using ϕ( x ), we can define the two orthonormal mode functions ϕ₁ and ϕ₂ at z = 0, as

$$\begin{equation}{\phi }_{1}(\boldsymbol{x},0)=\,\sqrt{2}\,\phi (\boldsymbol{x})\,\theta (y),\end{equation} \tag{ 39a }$$

$$\begin{equation}{\phi }_{2}(\boldsymbol{x},0)=\,\sqrt{2}\,\phi (\boldsymbol{x})\,\theta (-y),\end{equation} \tag{ 39b }$$

where θ(y) denotes the Heaviside step function [66]. Finally, imagine to send the upper (lower) mode ϕ₁ (ϕ₂) through a linear polarizer and a quarter wave plate suitably oriented, so that it is circularly polarised to the right (left). In this way, we obtain a non-uniformly polarised beam with vector amplitude

$$\begin{equation}\mathbf{\Psi }(\boldsymbol{x},0)=\phi (\boldsymbol{x})\left[{\hat{\boldsymbol{{\epsilon}}}}_{1}\,\theta (y)+{\hat{\boldsymbol{{\epsilon}}}}_{2}\,\theta (-y)\right],\end{equation} \tag{ 40 }$$

where the unit vectors

$$\begin{equation}{\hat{\boldsymbol{{\epsilon}}}}_{\alpha }=\frac{{\hat{\boldsymbol{e}}}_{1}+\mathrm{i}{(-1)}^{\alpha }{\hat{\boldsymbol{e}}}_{2}}{\sqrt{2}},\quad (\alpha =1,2),\end{equation} \tag{ 41 }$$

describe right- (α = 1) and left-circular (α = 2) polarisation. By definition, this field is maximally entangled and C = 1. The vector field Ψ( x , z) at distance z from the origin can be calculated as in equation (22).

The polarisation and intensity patterns of Ψ( x , z) at several distances z ⩾ 0 are shown in figures 2(a)–(d). At z = 0 (see figure 2(a)) the beam presents two disjoint polarisation patterns. Upon propagation by z, the upper and lower lobes of the field (41) spread and overlap, and at z = z_R the polarisation pattern becomes completely scrambled, as shown in figure 2(c).

Zoom In Zoom Out Reset image size

In figure 2(e) we show the degree of polarisation P₁(z) of the beam ${\mathbf{\Psi }}_{1}(\boldsymbol{x},z)={\left.\mathbf{\Psi }(\boldsymbol{x},z)\right\vert }_{\boldsymbol{x}\in {\mathcal{D}}_{1}}$, where ${\mathcal{D}}_{1}$ coincides with the upper half plane y ⩾ 0. For this beam P₁(z) = P₂(z). As expected, at z = z₀ = 0 the fields Ψ₁( x , 0) and Ψ₂( x , 0) are fully polarized. Then, the degree of polarisation decrease with z to reach the minimum value at z = z_R, where ${\left.{\varphi }_{1}(\boldsymbol{x},z)\right\vert }_{\boldsymbol{x}\in {\mathcal{D}}_{\mu }}$ and ${\left.{\varphi }_{2}(\boldsymbol{x},z)\right\vert }_{\boldsymbol{x}\in {\mathcal{D}}_{\mu }}$ are maximally overlapped. Finally, for z > z_R the fields tend to separate again and P₁(z) raises monotonically.

Figure 2(f) shows the plot of the Hellinger distance H[|ϕ₁|, |ϕ₂|](z) calculated from equation (36). It is worth noticing that the plots of H[|ϕ₁|, |ϕ₂|](z) and P₁(z) display the same qualitative behaviour. The Hellinger distance is maximal when the two polarisation patterns are fully separated. Vice versa, H[|ϕ₁|, |ϕ₂|](z) is minimal at z = z_R where the patterns are maximally superimposed.

3.1. Experimental setup

To experimentally prove our proposed theory, we implemented the experimental setup previously demonstrated in [40] and schematically depicted in figure 3 in the generation stage. Here, a horizontally polarised laser beam (λ = 532 nm, 300 mW) is expanded and collimated with the telescope formed by the pair of lenses L₁ (f₁ = 20 mm) and L₂ (f₂ = 200 mm). The polarisation state of the beam is then rotated to a diagonal polarisation state through a half-wave plate (HWP) oriented at 22.5° relative to the vertical axis. Afterwards, the beam is split, using a polarising beam splitter (PBS) or a Wallastone prism, into two different beams travelling along different optical paths, one with horizontal and the other with vertical polarisation. Both beams impinge onto a polarisation-independent digital micromirror device (DMD, DLP Light Crafter 6500 from Texas Instruments), at slightly different angles $(\approx 1.5{}^{\circ})$ and at the geometric centre of the DMD, where a digital hologram is displayed. The digital hologram contains the two constituting wave fields ϕ₁(x, y) and ϕ₂(x, y) required to generate the vector mode given by equation (5). Each field is superimposed with a unique linear phase grating and both multiplexed into a single hologram. The period of each linear grating is carefully adjusted to ensure the overlap of the two first diffraction orders along a common propagation axis after the DMD, where the desired complex vector mode is generated. Afterwards, such mode is isolated from other diffraction orders by means of a spatial filter (SF), located at the focusing plane of a telescope formed by lenses L₃ and L₄, both with focal lengths f = 100 mm. Finally, a quarter-wave plate (QW) is added to change the vector mode from the linear to the circular polarisation basis.

Zoom In Zoom Out Reset image size

The generated modes where analysed in the measurement stage through the use of a circular polariser (CP) fabricated by cementing a QWP to a linear polariser with its axis at 45° with respect to the fast axis of the QWP. Here, the generated modes were focused onto the CCD camera with a lens of focal length f = 400 mm. Transverse polarisation distributions were reconstructed using Stokes polarimetry, through a series of intensity measurements in accordance with the relations [41]

$$\begin{equation}\begin{aligned}\hfill & {S}_{0}={I}_{\mathrm{H}}+{I}_{\mathrm{V}},\hspace{54.06023pt}{S}_{1}=2{I}_{\mathrm{H}}-{S}_{0},\hfill \\ \hfill & {S}_{2}=2{I}_{\mathrm{D}}-{S}_{0},\hspace{54.06023pt}{S}_{3}=2{I}_{\mathrm{R}}-{S}_{0},\hfill \end{aligned}\end{equation} \tag{ 42 }$$

where, S₀, S₁, S₂ and S₃ are the Stokes parameters and I_H, I_V, I_D and I_R represent the intensity of the horizontal, vertical, diagonal and right-handed polarisation components, respectively. Experimentally, such intensities where recorded with a charge-coupled device (CCD) camera (FLIR FL3-U3-120S3C-C from point grey) in the measurement stage of figure 3 with the help of a CP. More specifically, I_H, I_V and I_D where obtained by setting the CP to 0°, 90° and 45°, respectively, whereas I_R was obtained by flipping the CP 180° with its angle fixed to 0° (see for example [40, 41] for further details).

3.2. Transverse polarisation distribution in spatially disjoint vector fields

In this section we will show three specific cases of engineered vector beams, whose polarisation structure changes dramatically upon free space propagation, namely, Hermite-, parabolic- and Airy–Gauss vector modes. The left panels of figure 4 show the transverse polarisation distribution overlapped on the intensity profile of such modes for different propagation distances, both experimentally and theoretically. The right ones show the position of the centroid of the right- and left-handed circular polarisation components as function of the propagation distance.

Zoom In Zoom Out Reset image size

3.2.1. First case: Hermite–Gauss vector beams

As a first case, we analyse the example mentioned in the theory section, constructed from the non-separable superposition of the HG₀₁ mode, whose spatial profile features two vertical lobes of maximum intensity. Such mode was engineered in such a way that the upper lobe carries right-handed circular polarisation, while the lower one carries the orthogonal left-handed one. This mode is then focused with a lens of focal distance f = 400 mm and scanned with the CCD camera along its propagation axis. The transverse polarisation distribution of such modes was reconstructed at selected planes z = 0, z = 0.5z_R, z = 0.75z_R and z = z_R, where z_R is the Rayleigh length. Such polarisation distributions overlapped on the spatial shape of the modes are shown in figure 4(a), experiment on the top row and numerical simulations, performed using the Rayleigh–Sommerfeld diffraction theory [64], on the bottom. Figure 4(b) further illustrates how both modes approach to each other while their polarisation distribution evolves from quasi-homogeneous to non homogeneous, as seen in the approaching centroids of left- and right-handed polarisation components upon propagation.

3.2.2. Second case: parabolic-Gauss beams

As a second example, we analyse the recently introduced parabolic-Gauss vector modes, which is implemented from a non-separable superposition of orthogonal parabolic modes, natural solutions to the Helmholtz equation in parabolic cylindrical coordinates, and orthogonal left- and right-handed polarisation states [42]. Such vector modes feature an interesting behaviour as function of their propagation in free space, namely, their degree of non-separability evolves from a non-homogeneously polarised vector beam to a quasi-homogeneously polarised one. Figure 4(c) illustrates this behaviour for different propagation distances z = 0, z = 0.25z_R, z = 0.4z_R and z = z_R. Here we show the intensity distribution overlapped on the transverse polarisation distribution, experiment on top and theory on the bottom, which we reconstructed through Stokes polarimetry. Notice that for z = 0 the transverse polarisation distribution contains only linear polarisation states whereas for z = z_R = f, the mode mainly contains left- and right-handed polarisation states. In figure 4(d) we show the centroid of each of the circular polarisation components as function of their propagation distance. In this case, the separation increases from zero to a maximum separation, which is attained in the far field or the focal plane of the Fourier transforming lens (z = f).

3.2.3. Third case: Airy–Gauss beams

As a final example, we engineered a vector mode consisting on the non-separable superposition of two Airy–Gauss beams with in-plane opposite transverse accelerations and orthogonal polarisation. As it is well-known, Airy–Gauss beams propagate in free space following a parabolic trajectory, whose parameters can be easily controlled [67]. We superimposed two Airy–Gauss beams with orthogonal circular polarisation with their maximum intensity lobes facing each other. The polarisation and intensity distribution of such modes, reconstructed for specific planes, is illustrated in figure 4(e), experiment on top and numerical simulations on the bottom. As illustrated, upon propagation the upper mode with right-handed circular polarisation moves upwards to a maximum point and then downwards. Similarly, the lower mode move downwards to a minimum point and then upwards. As result, upon propagation the modes initially separate from each other up to a maximum value, after which they approach to each other again. The trajectory of the main lobe of each beam as function of the propagation distance is shown in figure 4(f). Here, we choose arbitrarily z = 0 as the plane where the two modes are the closest to each other.

3.3. Determining the non-separability of spatially disjoint vector fields

As stated earlier, the degree of concurrence or non-separability of spatially disjoint vector fields, cannot be quantified using traditional methods, such as a basis independent tomography, as exhibited in [52]. In other words, traditional methods can only quantify the global behaviour of vector fields, but not the local behaviour exhibited by spatially disjoint vector fields. Hence, we introduced the Hellinger distance (equation (36)) as a measure that quantifies the local evolution of the non-separability in spatially disjoint vector fields. To demonstrate this, we applied equation (36) to the three vector fields introduced in the previous section. Figure 5(a) shows the case of the Hermite–Gauss vector field, where the Hellinger distance is plotted as function of propagation. For the sake of clarity, the transverse polarisation distribution at three different planes are shown as insets on the top of this figure, with arrows indicating their z value. As mentioned in the theory, the Hellinger distance, which evolves from 1 to 0, clearly captures the evolution of the mode, from two modes with orthogonal circular polarisation into a single mode with linear polarisation. The second case corresponds to the parabolic-Gauss vector modes shown in figure 5(b). Here, the transverse polarisation distribution evolves upon propagation from a single mode with linear polarisation to two spatially disjoint modes with circular orthogonal polarisation. Again, this behaviour is captured by the Hellinger distance, which in this case increases monotonically from 0 to 1. The final example corresponds to the Airy–Gauss vector modes shown in figure 5(c). This case was specifically engineered to show the effect of a parabolic-like Hellinger distance directly related to the transverse positions of the independent Airy–Gauss modes. In all cases, the theoretical results are shown as solid lines and the experimental data as solid circles, featuring remarkable agreement.

Zoom In Zoom Out Reset image size

Control over light's degrees of freedom is steadily gaining traction [68], driven by the explosion of applications that structured light affords. Concomitantly, there is a need to understand and quantify such new forms of light. In the early 1990s, exotic scalar fields became possible by the emergence of diffractive optical elements, which saw the measurement toolkit expanded to include statistical approaches (intensity moments), resulting in (for example) the beam quality factor, M². This single parameter is not sufficient to fully describe such light fields, and thus the toolkit quickly matured, forming the basis of our present ISO standards for laser beams. More recently, vectorially structured light is routinely created and employed in a variety of applications, but analogous to the situation for scalar light 30 years ago, our toolkit is very much in its infancy. This time our toolkit is quantum-inspired, returning a VQF (the vectorial equivalent to the M² parameter), but with some shortcomings that require continual improvement. The first was to alter the measurement from basis dependent [49] to a basis independent [52], removing the requirement of knowledge of the field to be probed. But here we have highlighted that even though these techniques perform very well at characterising vector modes whose spatial and polarisation degrees of freedom are shape-invariant upon propagation, it fails at characterising beams that do not satisfy such condition, whose constituting spatial modes separate from each other upon propagation, becoming path-like separable. The existing approaches are insufficient to quantify such disjoint behaviour. We envisaged a way in which we could address this issue by defining a Hellinger distance, between the constituting spatial modes that form the vector beam, as a measure. We hope that this work inspires further advances in the toolkit.

Our work also raises some interesting fundamental questions on the nature of visibility (spatial coherence) in the context of vectorially structured light. Our propagation induced disjointness is analogous to a double slit experiment, but as if the interference pattern is time reversed to return to the double slits. The question is how the visibility changes with propagation? In our analogous situation, the initial locally inhomogeneous field, with high visibility in the polarisation variation, evolves to become locally homogeneous with no visibility—we can observe (intuitively) two distinct regions made of scalar fields (our two 'slits'). Yet the global field retains its visibility during propagation, as if the original slits themselves are now the interference pattern with high visibility (1 inside the slits and 0 outside).

The situation has a quantum analogue too. While it is commonplace in quantum optics to exploit the path as a degree of freedom for entanglement, this is largely untouched in the context of classical quantum-like light [56–58]. When the paths are marked, no interference is observed. In our disjoint field, this is equivalent to measuring locally, one scalar component. When the paths are not marked, full visibility is possible. In our classical field, this is equivalent to making projections on the entire field, a global concurrence measurement. The physics of path embedded information in such vectorial light fields thus requires further probing, and will surely shed light on extending notions of spatial coherence to multi-partite, multiple degree of freedom light [69].

Here we have revisited the emerging field of quantum-like vectorial (classical) light, and highlighted the shortcomings of traditional quantum-inspired measures. We considered the general case where the local polarisation homogeneity can vary spatially across the field, from scalar to vector, by virtue of a propagation induced spatial disjointness. We proposed and demonstrated a new measure that correctly accounts for the local separability of the field, while accounting for the fact that the global vectorial field remains non-separable. Our work is relevant fundamentally, probing the notion of non-separability and spatial coherence in the context of vectorially polarised light, but also practically, where new standards to universally quantify these fields are a pressing need.

AA acknowledges support from the Deutsche Forschungsgemeinschaft Project No. 429529648-TRR 306 QuCoLiMa ('Quantum Cooperativity of Light and Matter'). CRG from the National Natural Science Foundation of China (61975047 and 11934013) and the High-Level Talents Project of Heilongjiang Province (Grant No. 2020GSP12). XBH from the Science Foundation of Zhejiang Sci-Tech University (ZSTU) under Grant No. 22062025-Y.

The data that support the findings of this study are available upon reasonable request from the authors.

The two-dimensional Hermite–Gauss modes u_nm ( x , z) = φ_n (x, z)φ_m (y, z), with n, m = 0, 1, ..., ∞, are defined by [70],

$$\begin{equation}{\varphi }_{n}(x,z)=\,\frac{1}{{\pi }^{1/4}}\frac{1}{\sqrt{{2}^{n}n!}}\,\frac{1}{\sqrt{\rho }}\,{\text{H}}_{n}\!\left(x/\rho \right){\mathrm{e}}^{-{(x/\rho )}^{2}/2}\,{\mathrm{e}}^{\frac{\mathrm{i}}{2}\frac{z}{{z}_{\mathrm{R}}}{(x/\rho )}^{2}}\,{\mathrm{e}}^{-\mathrm{i}\left(n+\frac{1}{2}\right)\mathrm{arctan}\left(z/{z}_{\mathrm{R}}\right)}.\end{equation} \tag{ A1 }$$

In equation (A1) H_n (x) denotes the nth-order Hermite polynomial, and

$$\begin{equation}\rho =\rho (z)=\frac{{w}_{0}}{\sqrt{2}}{\left(1+\frac{{z}^{2}}{{z}_{\mathrm{R}}^{2}}\right)}^{1/2},\end{equation} \tag{ A2 }$$

fixes the transverse length scale at distance z from the beam's origin, where the minimum beam radius w₀ > 0, is attained. The Rayleigh length ${z}_{\mathrm{R}}=k{w}_{0}^{2}/2$, sets the longitudinal length scale, giving the distance over which the beam can propagate without spreading significantly. From equations (A1) and (A2) it readily follows that

$$\begin{equation}{u}_{nm}(\boldsymbol{x},z)=\frac{1}{\xi (z)}{u}_{nm}\left(\frac{\boldsymbol{x}}{\xi (z)},0\right)\mathrm{exp}\left[\mathrm{i}{\chi }_{nm}(\boldsymbol{x},z)\right],\end{equation} \tag{ A3 }$$

with

$$\begin{equation}\xi (z)=\frac{\rho (z)}{\rho (0)}={\left(1+\frac{{z}^{2}}{{z}_{\mathrm{R}}^{2}}\right)}^{1/2},\end{equation} \tag{ A4 }$$

and

$$\begin{equation}{\chi }_{nm}(\boldsymbol{x},z)=\,\frac{z}{{z}_{\mathrm{R}}}{\left\vert \frac{1}{{w}_{0}}\frac{\boldsymbol{x}}{\xi (z)}\right\vert }^{2}-(n+m+1)\mathrm{arctan}\left(\frac{z}{{z}_{\mathrm{R}}}\right).\end{equation} \tag{ A5 }$$

In this appendix we quickly review the common use of a quantum-like notation for the study of vector beams. To begin with, let Ψ( x , z) be the vector amplitude of a non-uniformly polarized paraxial optical beam propagating along the z-axis of a Cartesian reference frame with coordinates (x, y, z) = ( x , z), defined by

$$\begin{equation}\mathbf{\Psi }(\boldsymbol{x},z)={\hat{\boldsymbol{e}}}_{1}\,{\psi }_{1}(\boldsymbol{x},z)+{\hat{\boldsymbol{e}}}_{2}\,{\psi }_{2}(\boldsymbol{x},z),\end{equation} \tag{ B1 }$$

where the real-valued unit vectors ${\hat{\boldsymbol{e}}}_{1}$ and ${\hat{\boldsymbol{e}}}_{2}$ are parallel to the x- and the y-axis, respectively, with ${\hat{\boldsymbol{e}}}_{1}\cdot {\hat{\boldsymbol{e}}}_{2}=0$, and ψ₁( x , z), ψ₂( x , z) are two arbitrary solutions of the paraxial wave equation, which do not need to be orthogonal and normalized. The global intensity of the beam is denoted I₀, and defined by

$$\begin{align}\hfill {I}_{0}=& {\int }_{{\mathbb{R}}^{2}}{\mathbf{\Psi }}^{\ast }(\boldsymbol{x},z)\cdot \mathbf{\Psi }(\boldsymbol{x},z)\,{\mathrm{d}}^{2}x\hfill \\ \hfill =& {\int }_{{\mathbb{R}}^{2}}\left\{{\left\vert {\psi }_{1}(\boldsymbol{x},z)\right\vert }^{2}+{\left\vert {\psi }_{2}(\boldsymbol{x},z)\right\vert }^{2}\right\}\,{\mathrm{d}}^{2}x,\hfill \end{align} \tag{ B2 }$$

which is independent of z.

Now, we introduce the discrete 'polarisation basis' {|1⟩, |2⟩}, and the continuous 'position basis' {| x ⟩ = |x, y⟩}, with $\boldsymbol{x}\in {\mathbb{R}}^{2}$. By definition, the polarisation basis {|i⟩}_i∈{1,2} spans the discrete polarisation Hilbert space ${\mathcal{H}}_{\text{pol}}$, and the position basis ${\left\{\vert \boldsymbol{x}\rangle \right\}}_{\boldsymbol{x}\in {\mathbb{R}}^{2}}$ spans the continuous position Hilbert space ${\mathcal{H}}_{\text{pos}}$. We postulate that these bases are orthonormal

$$\begin{equation}\langle \alpha \vert \beta \rangle =\,{\delta }_{\alpha \beta },\quad (\alpha ,\beta =1,2),\end{equation} \tag{ B3 }$$

$$\begin{align}\hfill \langle \boldsymbol{x}\vert {\boldsymbol{x}}^{\prime }\rangle & =\,\langle x\vert {x}^{\prime }\rangle \langle y\vert {y}^{\prime }\rangle \hfill \\ \hfill & =\,\delta (x-{x}^{\prime })\delta (y-{y}^{\prime })\hfill \\ \hfill & \equiv \,\delta (\boldsymbol{x}-{\boldsymbol{x}}^{\prime }),\quad (\boldsymbol{x},{\boldsymbol{x}}^{\prime }\in {\mathbb{R}}^{2}),\hfill \end{align} \tag{ B4 }$$

and complete:

$$\begin{equation}\sum\limits _{\alpha =1}^{2}\vert \alpha \rangle \!\langle \alpha \vert \,={\hat{I}}_{\text{pol}},\end{equation} \tag{ B5 }$$

$$\begin{equation}{\int }_{{\mathbb{R}}^{2}}\vert \boldsymbol{x}\rangle \!\langle \boldsymbol{x}\vert \,{\mathrm{d}}^{2}x\,={\hat{I}}_{\text{pos}},\end{equation} \tag{ B6 }$$

where ${\hat{I}}_{\text{pol}}$ and ${\hat{I}}_{\text{pos}}$ are the identity operators in the polarisation and position spaces, respectively. By definition, the direct product |α⟩ ⊗ | x ⟩ = |α⟩| x ⟩ = |α, x ⟩ spans the whole Hilbert space $\mathcal{H}={\mathcal{H}}_{\text{pol}}\otimes {\mathcal{H}}_{\text{pos}}$.

Using this quantum-like notation (but the physics of the problem is entirely classical), we can represent Ψ( x , z) by means of the state vector |Ψ(z)⟩ defined by

$$\begin{align}\hfill \vert {\Psi}(z)\rangle & =\vert 1\rangle \vert {\psi }_{1}(z)\rangle +\vert 2\rangle \vert {\psi }_{2}(z)\rangle \hfill \\ \hfill & =\sum\limits _{\beta =1}^{2}\vert \beta \rangle \vert {\psi }_{\beta }(z)\rangle ,\hfill \end{align} \tag{ B7 }$$

such that

$$\begin{align}\hfill \langle \alpha ,\boldsymbol{x}\vert {\Psi}(z)\rangle & =\sum\limits _{\beta =1}^{2}\langle \alpha \vert \beta \rangle \langle \boldsymbol{x}\vert {\psi }_{\beta }(z)\rangle \hfill \\ \hfill & =\langle \boldsymbol{x}\vert {\psi }_{\alpha }(z)\rangle \hfill \\ \hfill & ={\psi }_{\alpha }(\boldsymbol{x},z),\quad (\alpha =1,2),\hfill \end{align} \tag{ B8 }$$

where ψ₁( x , z), ψ₂( x , z) are defined by (B1). For the sake of clarity, here and hereafter we shall occasionally write |ψ_α ⟩ for |ψ_α (z)⟩ and |Ψ⟩ for |Ψ(z)⟩, the z-dependence being understood. Moreover, we shall use upper case Greek letters, as Ψ, to denote vectors in $\mathcal{H}$, while lower case Greek letters, as ψ, will denote vectors in ${\mathcal{H}}_{\text{pos}}$.

Using this notation, we can rewrite the intensity of the beam as

$$\begin{align}\hfill {I}_{0}& ={\int }_{{\mathbb{R}}^{2}}\left\{{\left\vert {\psi }_{1}(\boldsymbol{x},z)\right\vert }^{2}+{\left\vert {\psi }_{2}(\boldsymbol{x},z)\right\vert }^{2}\right\}{\mathrm{d}}^{2}x\hfill \\ \hfill & =\langle {\psi }_{1}\vert {\psi }_{1}\rangle +\langle {\psi }_{2}\vert {\psi }_{2}\rangle \hfill \\ \hfill & =\langle {\Psi}\vert {\Psi}\rangle ,\hfill \end{align} \tag{ B9 }$$

where, here and hereafter, for any $\vert \psi \rangle ,\vert \phi \rangle \in {\mathcal{H}}_{\text{pos}}$, and $\vert {\Psi}\rangle \in \mathcal{H}$, we write

$$\begin{align}\hfill \langle \psi \vert \phi \rangle & =\langle \psi \vert {\hat{I}}_{\text{pos}}\vert \phi \rangle \hfill \\ \hfill & ={\int }_{{\mathbb{R}}^{2}}\langle \psi \vert \boldsymbol{x}\rangle \langle \boldsymbol{x}\vert \phi \rangle \,{\mathrm{d}}^{2}x\hfill \\ \hfill & ={\int }_{{\mathbb{R}}^{2}}{\psi }^{\ast }(\boldsymbol{x},z)\phi (\boldsymbol{x},z)\,{\mathrm{d}}^{2}x,\hfill \end{align} \tag{ B10 }$$

and

$$\begin{align}\hfill \langle {\Psi}\vert {\Psi}\rangle & =\langle {\Psi}\vert {\hat{I}}_{\text{pol}}\otimes {\hat{I}}_{\text{pos}}\vert {\Psi}\rangle \hfill \\ \hfill & =\sum\limits _{\alpha =1}^{2}{\int }_{{\mathbb{R}}^{2}}\langle {\Psi}\vert \alpha ,\boldsymbol{x}\rangle \langle \alpha ,\boldsymbol{x}\vert {\Psi}\rangle \,{\mathrm{d}}^{2}x\hfill \\ \hfill & =\sum\limits _{\alpha =1}^{2}{\int }_{{\mathbb{R}}^{2}}\langle {\psi }_{\alpha }\vert \boldsymbol{x}\rangle \langle \boldsymbol{x}\vert {\psi }_{\alpha }\rangle \,{\mathrm{d}}^{2}x\hfill \\ \hfill & =\sum\limits _{\alpha =1}^{2}\langle {\psi }_{\alpha }\vert {\psi }_{\alpha }\rangle ,\hfill \end{align} \tag{ B11 }$$

where (B8) has been used.

Given the field

$$\begin{equation}\mathbf{\Psi }(\boldsymbol{x},z)={\hat{\boldsymbol{e}}}_{1}\,{\psi }_{1}(\boldsymbol{x},z)+{\hat{\boldsymbol{e}}}_{2}\,{\psi }_{2}(\boldsymbol{x},z),\end{equation} \tag{ C1 }$$

and the corresponding state vector

$$\begin{equation}\vert {\Psi}(z)\rangle =\vert 1\rangle \vert {\psi }_{1}(z)\rangle +\vert 2\rangle \vert {\psi }_{2}(z)\rangle ,\end{equation} \tag{ C2 }$$

the Schmidt decomposition is easy to evaluate. First, using |ψ₁⟩ and |ψ₂⟩ we build the orthonormal basis {|v₁⟩, |v₂⟩} defined by

$$\begin{equation}\vert {v}_{1}\rangle =\frac{\vert {\psi }_{1}\rangle }{\sqrt{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }},\end{equation} \tag{ C3 }$$

and

$$\begin{align}\hfill \vert {v}_{2}\rangle & =\frac{\vert {\psi }_{2}\rangle -\vert {v}_{1}\rangle \langle {v}_{1}\vert {\psi }_{2}\rangle }{\sqrt{\langle {\psi }_{2}\vert {\psi }_{2}\rangle -{\left\vert \langle {v}_{1}\vert {\psi }_{2}\rangle \right\vert }^{2}}}\hfill \\ \hfill & =\frac{\vert {\psi }_{2}\rangle -\vert {\psi }_{1}\rangle \frac{\langle {\psi }_{1}\vert {\psi }_{2}\rangle }{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }}{\sqrt{\langle {\psi }_{2}\vert {\psi }_{2}\rangle -\frac{{\left\vert \langle {\psi }_{1}\vert {\psi }_{2}\rangle \right\vert }^{2}}{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }}}.\hfill \end{align} \tag{ C4 }$$

Equations (C3) and (C4) imply

$$\begin{equation}\vert {\psi }_{1}\rangle =\vert {v}_{1}\rangle \sqrt{\langle {\psi }_{1}\vert {\psi }_{1}\rangle },\end{equation} \tag{ C5 }$$

and

$$\begin{equation}\vert {\psi }_{2}\rangle =\vert {v}_{1}\rangle \frac{\langle {\psi }_{1}\vert {\psi }_{2}\rangle }{\sqrt{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }}+\vert {v}_{2}\rangle \sqrt{\langle {\psi }_{2}\vert {\psi }_{2}\rangle -\frac{{\left\vert \langle {\psi }_{1}\vert {\psi }_{2}\rangle \right\vert }^{2}}{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }}.\end{equation} \tag{ C6 }$$

Substituting (C5) and (C6) into (C2) we find

$$\begin{align}\hfill \vert {\Psi}(z)\rangle & =\vert 1\rangle \,\vert {v}_{1}\rangle \left(\sqrt{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }\right)+\vert 2\rangle \,\vert {v}_{1}\rangle \left(\,\frac{\langle {\psi }_{1}\vert {\psi }_{2}\rangle }{\sqrt{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }}\right)\hfill \\ \hfill & \quad +\vert 2\rangle \,\vert {v}_{2}\rangle \left(\,\sqrt{\langle {\psi }_{2}\vert {\psi }_{2}\rangle -\frac{{\left\vert \langle {\psi }_{1}\vert {\psi }_{2}\rangle \right\vert }^{2}}{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }}\right)\hfill \\ \hfill & \equiv \sum\limits _{\alpha ,\beta =1}^{2}{A}_{\alpha \beta }\,\vert \alpha \rangle \vert {v}_{\beta }\rangle .\hfill \end{align} \tag{ C7 }$$

The 2 × 2 matrix A is defined by (C7) as

$$\begin{equation}A=\left[\begin{matrix}\hfill \sqrt{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }\hfill & \hfill 0\hfill \\ \hfill \frac{\langle {\psi }_{1}\vert {\psi }_{2}\rangle }{\sqrt{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }}\hfill & \hfill \sqrt{\langle {\psi }_{2}\vert {\psi }_{2}\rangle -\frac{{\left\vert \langle {\psi }_{1}\vert {\psi }_{2}\rangle \right\vert }^{2}}{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }}\hfill \end{matrix}\right],\end{equation} \tag{ C8 }$$

and it is independent of z. Substituting the singular value decomposition

$$\begin{equation}A=UD{V}^{{\dagger}},\end{equation} \tag{ C9 }$$

into (C7), we can easily calculate the Schmidt form of |Ψ(z)⟩ as:

$$\begin{align}\hfill \vert {\Psi}(z)\rangle & =\sum\limits _{\alpha =1}^{2}{D}_{\alpha }\left(\sum\limits _{\beta =1}^{2}{U}_{\beta \alpha }\vert \beta \rangle \right)\left(\sum\limits _{\gamma =1}^{2}{V}_{\alpha \gamma }^{{\dagger}}\vert {v}_{\gamma }\rangle \right)\hfill \\ \hfill & \equiv \sum\limits _{\alpha =1}^{2}\sqrt{{\lambda }_{\alpha }}\,\vert {\varepsilon }_{\alpha }\rangle \vert {\phi }_{\alpha }\rangle ,\hfill \end{align} \tag{ C10 }$$

where ${D}_{\alpha \beta }={\delta }_{\alpha \beta }{D}_{\alpha }\equiv {\delta }_{\alpha \beta }\sqrt{{\lambda }_{\alpha }}$, λ_α being the non-negative eigenvalues of A^† A, and we have defined

$$\begin{equation}\begin{aligned}\hfill \vert {\varepsilon }_{\alpha }\rangle & \equiv \sum\limits _{\beta =1}^{2}{U}_{\beta \alpha }\vert \beta \rangle ,\hfill \\ \hfill \vert {\phi }_{\alpha }\rangle & \equiv \sum\limits _{\gamma =1}^{2}{V}_{\alpha \gamma }^{{\dagger}}\vert {v}_{\gamma }\rangle .\hfill \end{aligned}\end{equation} \tag{ C11 }$$

A straightforward calculation gives

$$\begin{equation}\begin{aligned}\hfill {\lambda }_{1}& ={\lambda }_{+},\hfill \\ \hfill {\lambda }_{2}& ={\lambda }_{-},\hfill \end{aligned}\end{equation} \tag{ C12 }$$

where

$$\begin{equation}{\lambda }_{\pm }=\,\frac{1}{2}\left[\langle {\psi }_{1}\vert {\psi }_{1}\rangle +\langle {\psi }_{2}\vert {\psi }_{2}\rangle \pm \sqrt{4{\left\vert \langle {\psi }_{1}\vert {\psi }_{2}\rangle \right\vert }^{2}+{\left(\langle {\psi }_{1}\vert {\psi }_{1}\rangle -\langle {\psi }_{2}\vert {\psi }_{2}\rangle \right)}^{2}}\right].\end{equation} \tag{ C13 }$$

Note that since the state vector |Ψ(z)⟩ is not normalized, then

$$\begin{equation}{\lambda }_{1}+{\lambda }_{2}=\langle {\psi }_{1}\vert {\psi }_{1}\rangle +\langle {\psi }_{2}\vert {\psi }_{2}\rangle \ne 1.\end{equation} \tag{ C14 }$$

Finally, multiplying (C10) from left by ⟨ x | and using equations (B7) and (B8), we obtain the Schmidt form of the field (C1):

$$\begin{equation}\mathbf{\Psi }(\boldsymbol{x},z)=\sqrt{{\lambda }_{1}}\,{\hat{\boldsymbol{\varepsilon }}}_{1}\,{\phi }_{1}(\boldsymbol{x},z)+\sqrt{{\lambda }_{2}}\,{\hat{\boldsymbol{\varepsilon }}}_{2}\,{\phi }_{2}(\boldsymbol{x},z).\end{equation} \tag{ C15 }$$

The amount of the entanglement in the vector state |Ψ(z)⟩ can be quantified by the so-called Schmidt number, denoted K and defined by

$$\begin{align}\hfill K& =\frac{{\left({\lambda }_{1}+{\lambda }_{2}\right)}^{2}}{{\lambda }_{1}^{2}+{\lambda }_{2}^{2}}\hfill \\ \hfill & =\frac{{\left(\langle {\psi }_{1}\vert {\psi }_{1}\rangle +\langle {\psi }_{2}\vert {\psi }_{2}\rangle \right)}^{2}}{{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }^{2}+2{\left\vert \langle {\psi }_{1}\vert {\psi }_{2}\rangle \right\vert }^{2}+{\langle {\psi }_{2}\vert {\psi }_{2}\rangle }^{2}}.\hfill \end{align} \tag{ C16 }$$

For a bipartite two-dimensional system the Schmidt number K and the concurrence C are simply related:

$$\begin{equation}K=\,\frac{2}{2-{C}^{2}},\end{equation} \tag{ C17a }$$

$$\begin{equation}C=\,\sqrt{2\left(1-\frac{1}{K}\right)}=2\,\frac{\sqrt{{\lambda }_{1}{\lambda }_{2}}}{{\lambda }_{1}+{\lambda }_{2}}.\end{equation} \tag{ C17b }$$

We present an explicit calculation of C in appendix E.

The density matrix operator $\hat{\rho }$ for the pure state |Ψ⟩ is defined as the projector

$$\begin{equation}\hat{\rho }=\frac{\vert {\Psi}\rangle \!\langle {\Psi}\vert }{\langle {\Psi}\vert {\Psi}\rangle }.\end{equation} \tag{ D1 }$$

The reduced polarisation density matrix operator ${\hat{\rho }}_{\text{pol}}$, is calculated by tracing $\hat{\rho }$ with respect to the position degrees of freedom, as follows:

$$\begin{align}\hfill {\hat{\rho }}_{\text{pol}}& ={\int }_{{\mathbb{R}}^{2}}\langle \boldsymbol{x}\vert \hat{\rho }\vert \boldsymbol{x}\rangle \,{\mathrm{d}}^{2}x\hfill \\ \hfill & =\frac{1}{\langle {\Psi}\vert {\Psi}\rangle }{\int }_{{\mathbb{R}}^{2}}\langle \boldsymbol{x}\vert {\Psi}\rangle \langle {\Psi}\vert \boldsymbol{x}\rangle \,{\mathrm{d}}^{2}x\hfill \\ \hfill & =\frac{1}{\langle {\Psi}\vert {\Psi}\rangle }\sum\limits _{\alpha ,\beta =1}^{2}\vert \alpha \rangle \!\langle \beta \vert {\int }_{{\mathbb{R}}^{2}}\langle {\psi }_{\beta }\vert \boldsymbol{x}\rangle \langle \boldsymbol{x}\vert {\psi }_{\alpha }\rangle \,{\mathrm{d}}^{2}x\hfill \\ \hfill & =\frac{1}{\langle {\Psi}\vert {\Psi}\rangle }\sum\limits _{\alpha ,\beta =1}^{2}\vert \alpha \rangle \!\langle \beta \vert \langle {\psi }_{\beta }\vert {\psi }_{\alpha }\rangle .\hfill \end{align} \tag{ D2 }$$

It is easy to check that $\mathrm{t}\mathrm{r}\,{\hat{\rho }}_{\text{pol}}=1$. This equations shows that in the polarisation basis {|1⟩, |2⟩}, ${\hat{\rho }}_{\text{pol}}$ is represented by the 2 × 2 normalized coherency matrix of the beam, denoted ρ_pol = J/tr J, and defined by

$$\begin{align}\hfill {\hat{\rho }}_{\text{pol}}& \doteq {\rho }_{\text{pol}}\hfill \\ \hfill & =\left[\begin{matrix}\hfill \langle 1\vert {\hat{\rho }}_{\text{pol}}\vert 1\rangle \hfill & \hfill \langle 1\vert {\hat{\rho }}_{\text{pol}}\vert 2\rangle \hfill \\ \hfill \langle 2\vert {\hat{\rho }}_{\text{pol}}\vert 1\rangle \hfill & \hfill \langle 2\vert {\hat{\rho }}_{\text{pol}}\vert 2\rangle \hfill \end{matrix}\right]\hfill \\ \hfill & =\frac{1}{\langle {\Psi}\vert {\Psi}\rangle }\left[\begin{matrix}\hfill \langle {\psi }_{1}\vert {\psi }_{1}\rangle \hfill & \hfill \langle {\psi }_{2}\vert {\psi }_{1}\rangle \hfill \\ \hfill \langle {\psi }_{1}\vert {\psi }_{2}\rangle \hfill & \hfill \langle {\psi }_{2}\vert {\psi }_{2}\rangle \hfill \end{matrix}\right]\hfill \\ \hfill & =\frac{1}{\mathrm{t}\mathrm{r}\,J}\,J,\hfill \end{align} \tag{ D3 }$$

where tr J = ⟨Ψ|Ψ⟩ = ⟨ψ₁|ψ₁⟩ + ⟨ψ₂|ψ₂⟩. A direct calculation shows that

$$\begin{equation}J=A{A}^{{\dagger}},\end{equation} \tag{ D4 }$$

where A is defined by equation (C8). This implies that J and A^† A have the same eigenvalues λ₁, λ₂ given in (C12).

The coherency matrix J can be written in terms of the so-called Stokes parameters {S₀, S₁, S₂, S₃} as,

$$\begin{equation}J=\,\frac{1}{2}\left[\begin{matrix}\hfill {S}_{0}+{S}_{1}\hfill & \hfill {S}_{2}-\mathrm{i}{S}_{3}\hfill \\ \hfill {S}_{2}+\mathrm{i}{S}_{3}\hfill & \hfill {S}_{0}-{S}_{1}\hfill \end{matrix}\right].\end{equation} \tag{ D5 }$$

Comparing (D3) and (D5) we obtain

$$\begin{equation}\frac{1}{2}\left[\begin{matrix}\hfill {S}_{0}+{S}_{1}\hfill & \hfill {S}_{2}-\mathrm{i}{S}_{3}\hfill \\ \hfill {S}_{2}+\mathrm{i}{S}_{3}\hfill & \hfill {S}_{0}-{S}_{1}\hfill \end{matrix}\right]=\left[\begin{matrix}\hfill \langle {\psi }_{1}\vert {\psi }_{1}\rangle \hfill & \hfill \langle {\psi }_{2}\vert {\psi }_{1}\rangle \hfill \\ \hfill \langle {\psi }_{1}\vert {\psi }_{2}\rangle \hfill & \hfill \langle {\psi }_{2}\vert {\psi }_{2}\rangle \hfill \end{matrix}\right].\end{equation} \tag{ D6 }$$

This evidently implies

$$\begin{equation}\begin{aligned}\hfill \frac{{S}_{0}+{S}_{1}}{2}& =\langle {\psi }_{1}\vert {\psi }_{1}\rangle ,\hfill \\ \hfill \frac{{S}_{0}-{S}_{1}}{2}& =\langle {\psi }_{2}\vert {\psi }_{2}\rangle ,\hfill \\ \hfill \frac{{S}_{2}-\mathrm{i}{S}_{3}}{2}& =\langle {\psi }_{2}\vert {\psi }_{1}\rangle ,\hfill \\ \hfill \frac{{S}_{2}+\mathrm{i}{S}_{3}}{2}& =\langle {\psi }_{1}\vert {\psi }_{2}\rangle .\hfill \end{aligned}\end{equation} \tag{ D7 }$$

Solving this set of four equations with respect to {S₀, S₁, S₂, S₃}, we find

$$\begin{equation}\begin{aligned}\hfill {S}_{0}& =\langle {\psi }_{1}\vert {\psi }_{1}\rangle +\langle {\psi }_{2}\vert {\psi }_{2}\rangle ,\hfill \\ \hfill {S}_{1}& =\langle {\psi }_{1}\vert {\psi }_{1}\rangle -\langle {\psi }_{2}\vert {\psi }_{2}\rangle ,\hfill \\ \hfill {S}_{2}& =\langle {\psi }_{1}\vert {\psi }_{2}\rangle +\langle {\psi }_{2}\vert {\psi }_{1}\rangle ,\hfill \\ \hfill {S}_{3}& =\frac{\langle {\psi }_{1}\vert {\psi }_{2}\rangle -\langle {\psi }_{2}\vert {\psi }_{1}\rangle }{\mathrm{i}},\hfill \end{aligned}\end{equation} \tag{ D8 }$$

so that we can rewrite (D3) as

$$\begin{align}\hfill {\rho }_{\text{pol}}& =\frac{1}{\langle {\psi }_{1}\vert {\psi }_{1}\rangle +\langle {\psi }_{2}\vert {\psi }_{2}\rangle }\left[\begin{matrix}\hfill \langle {\psi }_{1}\vert {\psi }_{1}\rangle \hfill & \hfill \langle {\psi }_{2}\vert {\psi }_{1}\rangle \hfill \\ \hfill \langle {\psi }_{1}\vert {\psi }_{2}\rangle \hfill & \hfill \langle {\psi }_{2}\vert {\psi }_{2}\rangle \hfill \end{matrix}\right]\hfill \\ \hfill & =\frac{1}{2{S}_{0}}\left[\begin{matrix}\hfill {S}_{0}+{S}_{1}\hfill & \hfill {S}_{2}-\mathrm{i}{S}_{3}\hfill \\ \hfill {S}_{2}+\mathrm{i}{S}_{3}\hfill & \hfill {S}_{0}-{S}_{1}\hfill \end{matrix}\right].\hfill \end{align} \tag{ D9 }$$

The concurrence C of the state vector |Ψ⟩ is defined as [71],

$$\begin{equation}C=\sqrt{2\left(1-\mathrm{t}\mathrm{r}\,{\rho }_{\text{pol}}^{2}\right)},\end{equation} \tag{ E1 }$$

where

$$\begin{align}\hfill \mathrm{t}\mathrm{r}\,{\rho }_{\text{pol}}^{2}& =\frac{\mathrm{t}\mathrm{r}\left\{{\left[\begin{matrix}\hfill \langle {\psi }_{1}\vert {\psi }_{1}\rangle \hfill & \hfill \langle {\psi }_{2}\vert {\psi }_{1}\rangle \hfill \\ \hfill \langle {\psi }_{1}\vert {\psi }_{2}\rangle \hfill & \hfill \langle {\psi }_{2}\vert {\psi }_{2}\rangle \hfill \end{matrix}\right]}^{2}\right\}}{{\langle {\Psi}\vert {\Psi}\rangle }^{2}}\hfill \\ \hfill & =\frac{{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }^{2}+2\langle {\psi }_{2}\vert {\psi }_{1}\rangle \langle {\psi }_{1}\vert {\psi }_{2}\rangle +{\langle {\psi }_{2}\vert {\psi }_{2}\rangle }^{2}}{{\left(\langle {\psi }_{1}\vert {\psi }_{1}\rangle +\langle {\psi }_{2}\vert {\psi }_{2}\rangle \right)}^{2}}\hfill \\ \hfill & =\frac{{\left({S}_{0}+{S}_{1}\right)}^{2}+2{\left\vert {S}_{2}-\mathrm{i}{S}_{3}\right\vert }^{2}+{\left({S}_{0}-{S}_{1}\right)}^{2}}{4{S}_{0}^{2}}\hfill \\ \hfill & =\frac{1}{2}+\frac{{S}_{1}^{2}+{S}_{2}^{2}+{S}_{3}^{2}}{2{S}_{0}^{2}},\hfill \end{align} \tag{ E2 }$$

and equation (D7) has been used. Substituting (E2) into (E1) we obtain

$$\begin{equation}C=\sqrt{1-\frac{{S}_{1}^{2}+{S}_{2}^{2}+{S}_{3}^{2}}{{S}_{0}^{2}}}.\end{equation} \tag{ E3 }$$

Using the second line of (E2) and ⟨Ψ|Ψ⟩ = ⟨ψ₁|ψ₁⟩ + ⟨ψ₂|ψ₂⟩, we can also write the concurrence as

$$\begin{equation}C=2\,\frac{\sqrt{\langle {\psi }_{2}\vert {\psi }_{2}\rangle \langle {\psi }_{1}\vert {\psi }_{1}\rangle -\langle {\psi }_{2}\vert {\psi }_{1}\rangle \langle {\psi }_{1}\vert {\psi }_{2}\rangle }}{\langle {\psi }_{1}\vert {\psi }_{1}\rangle +\langle {\psi }_{2}\vert {\psi }_{2}\rangle }.\end{equation} \tag{ E4 }$$

Then we remember that the Schmidt decomposition of |Ψ⟩ is (see appendix C),

$$\begin{align}\hfill \vert {\Psi}\rangle & =\sqrt{{\lambda }_{1}}\,\vert {\varepsilon }_{1}\rangle \vert {\phi }_{1}\rangle +\sqrt{{\lambda }_{2}}\,\vert {\varepsilon }_{2}\rangle \vert {\phi }_{2}\rangle \hfill \\ \hfill & \equiv \vert {\varepsilon }_{1}\rangle \vert {\psi }_{1}\rangle +\vert {\varepsilon }_{2}\rangle \vert {\psi }_{2}\rangle ,\hfill \end{align} \tag{ E5 }$$

where

$$\begin{equation}\langle {\varepsilon }_{\alpha }\vert {\varepsilon }_{\beta }\rangle ={\delta }_{\alpha \beta },\quad \text{and}\quad \langle {\phi }_{\alpha }\vert {\phi }_{\beta }\rangle ={\delta }_{\alpha \beta },\end{equation} \tag{ E6 }$$

with α, β = 1, 2. Substituting (E5) into (E4) we obtain

$$\begin{align}\hfill C& =2\,\frac{\sqrt{{\lambda }_{1}{\lambda }_{2}}}{{\lambda }_{1}+{\lambda }_{2}}\sqrt{\langle {\phi }_{2}\vert {\phi }_{2}\rangle \langle {\phi }_{1}\vert {\phi }_{1}\rangle -\langle {\phi }_{2}\vert {\phi }_{1}\rangle \langle {\phi }_{1}\vert {\phi }_{2}\rangle }\hfill \\ \hfill & =2\,\frac{\sqrt{{\lambda }_{1}{\lambda }_{2}}}{{\lambda }_{1}+{\lambda }_{2}},\hfill \end{align} \tag{ E7 }$$

where (E6) has been used.

As a consistency check, we repeat now the same calculation for C, but using the equivalent definition

$$\begin{equation}C=\sqrt{2\left(1-\mathrm{t}\mathrm{r}\,{\hat{\rho }}_{\text{pos}}^{2}\right)},\end{equation} \tag{ E8 }$$

where the reduced position density matrix operator ${\hat{\rho }}_{\text{pos}}$ is calculated by tracing $\hat{\rho }$ with respect to the polarisation degrees of freedom,

$$\begin{align}\hfill {\hat{\rho }}_{\text{pos}}& =\sum\limits _{\alpha =1}^{2}\langle \alpha \vert \hat{\rho }\vert \alpha \rangle \hfill \\ \hfill & =\frac{1}{\langle {\Psi}\vert {\Psi}\rangle }\sum\limits _{\alpha =1}^{2}\langle \alpha \vert {\Psi}\rangle \langle {\Psi}\vert \alpha \rangle \hfill \\ \hfill & =\frac{1}{\langle {\Psi}\vert {\Psi}\rangle }\left(\vert {\psi }_{1}\rangle \!\langle {\psi }_{1}\vert +\vert {\psi }_{2}\rangle \!\langle {\psi }_{2}\vert \right).\hfill \end{align} \tag{ E9 }$$

Note that this quantity is not a 2 × 2 matrix, but an infinite-dimensional operator, as it can be directly seen by rewriting ${\hat{\rho }}_{\text{pos}}$ in the position basis:

$$\begin{align}\hfill {\hat{\rho }}_{\text{pos}}& ={\hat{I}}_{\text{pos}}\,{\hat{\rho }}_{\text{pos}}\,{\hat{I}}_{\text{pos}}\hfill \\ \hfill & ={\int }_{{\mathbb{R}}^{2}}{\mathrm{d}}^{2}x{\int }_{{\mathbb{R}}^{2}}{\mathrm{d}}^{2}{x}^{\prime }\,\vert \boldsymbol{x}\rangle \langle \boldsymbol{x}\vert {\hat{\rho }}_{\text{pos}}\vert {\boldsymbol{x}}^{\prime }\rangle \langle {\boldsymbol{x}}^{\prime }\vert \hfill \\ \hfill & =\frac{1}{\langle {\Psi}\vert {\Psi}\rangle }{\int }_{{\mathbb{R}}^{2}}{\mathrm{d}}^{2}x{\int }_{{\mathbb{R}}^{2}}{\mathrm{d}}^{2}{x}^{\prime }\,\vert \boldsymbol{x}\rangle \langle {\boldsymbol{x}}^{\prime }\vert \left[{\psi }_{1}(\boldsymbol{x},z){\psi }_{1}^{\ast }({\boldsymbol{x}}^{\prime },z)+{\psi }_{2}(\boldsymbol{x},z){\psi }_{2}^{\ast }({\boldsymbol{x}}^{\prime },z)\right].\hfill \end{align} \tag{ E10 }$$

From this equation it follows that

$$\begin{align}\hfill \mathrm{t}\mathrm{r}\,{\hat{\rho }}_{\text{pos}}^{2}& =\frac{1}{{\langle {\Psi}\vert {\Psi}\rangle }^{2}}\,\mathrm{t}\mathrm{r}\left[\left(\sum\limits _{\alpha =1}^{2}\vert {\psi }_{\alpha }\rangle \!\langle {\psi }_{\alpha }\vert \right)\left(\sum\limits _{\beta =1}^{2}\vert {\psi }_{\beta }\rangle \!\langle {\psi }_{\beta }\vert \right)\right]\hfill \\ \hfill & =\frac{1}{{\langle {\Psi}\vert {\Psi}\rangle }^{2}}\sum\limits _{\alpha ,\beta =1}^{2}\langle {\psi }_{\alpha }\vert {\psi }_{\beta }\rangle \mathrm{t}\mathrm{r}\left(\vert {\psi }_{\alpha }\rangle \!\langle {\psi }_{\beta }\vert \right)\hfill \\ \hfill & =\frac{1}{{\langle {\Psi}\vert {\Psi}\rangle }^{2}}\sum\limits _{\alpha ,\beta =1}^{2}{\left\vert \langle {\psi }_{\alpha }\vert {\psi }_{\beta }\rangle \right\vert }^{2}\hfill \\ \hfill & =\frac{1}{{\langle {\Psi}\vert {\Psi}\rangle }^{2}}\left[{\langle {\psi }_{1}\vert {\psi }_{1}\rangle }^{2}+2{\left\vert \langle {\psi }_{1}\vert {\psi }_{2}\rangle \right\vert }^{2}+{\langle {\psi }_{2}\vert {\psi }_{2}\rangle }^{2}\right]\hfill \\ \hfill & =\frac{1}{2}+\frac{{S}_{1}^{2}+{S}_{2}^{2}+{S}_{3}^{2}}{2{S}_{0}^{2}},\hfill \end{align} \tag{ E11 }$$

where (D7) has been used, and

$$\begin{align}\hfill \mathrm{t}\mathrm{r}\left(\vert {\psi }_{\alpha }\rangle \!\langle {\psi }_{\beta }\vert \right)& ={\int }_{{\mathbb{R}}^{2}}\langle \boldsymbol{x}\vert {\psi }_{\alpha }\rangle \langle {\psi }_{\beta }\vert \boldsymbol{x}\rangle \,{\mathrm{d}}^{2}x\hfill \\ \hfill & =\langle {\psi }_{\beta }\vert {\psi }_{\alpha }\rangle ,\quad (\alpha ,\beta =1,2).\hfill \end{align} \tag{ E12 }$$

Substituting (E11) into (E8) we obtain

$$\begin{equation}C=\sqrt{1-\frac{{S}_{1}^{2}+{S}_{2}^{2}+{S}_{3}^{2}}{{S}_{0}^{2}}}.\end{equation} \tag{ E13 }$$

As expected, we have obtained the same value for the concurrence C by using either (E1) or (E8).