Spin-orbit coupling in photonic graphene

We generate experimentally a honeycomb refractive index pattern in an atomic vapor cell using electromagnetically-induced transparency. We study experimentally and theoretically the propagation of polarized light beams in such"photonic graphene". We demonstrate that an effective spin-orbit coupling appears as a correction to the paraxial beam equations because of the strong spatial gradients of the permittivity. It leads to the coupling of spin and angular momentum at the Dirac points of the graphene lattice. Our results suggest that the polarization degree plays an important role in many configurations where it has been previously neglected.


INTRODUCTION
Topological photonics [1, 2] is a rapidly growing field, combining fundamental physics and applied optics. The research in this field has brought us new understanding of the fundamental topological properties of optical systems, which are due to the photonic spin-orbit coupling present in various kinds of inhomogeneous photonic systems [3][4][5]. Photonic spin-orbit coupling (SOC) is a crucial ingredient for solving long-standing problems like optical isolation at a microscopic scale [6][7][8][9][10] required for the functioning of lasers, opening a new field of topological lasers [11][12][13][14].
Photonic graphene is a system of a particular interest. The graphene lattice, studied for more than a half-century [15], was one of the first to demonstrate the striking manifestations of the Dirac physics [16] such as the Klein tunneling [17], with enormous potential for applications, which have already found their way to the market [18]. In photonics, the Dirac points of graphene-like lattices offer extended possibilities for the manipulation of optical angular momentum [19] and for the studies of singular optical beams [20]. Recent works address various problems, such as beam conversion [21], intervalley scattering [22], and valley pseudospin dynamics [23]. Different implemen-tations of photonic graphene include coupled waveguides [24], microwave resonators [25], photorefractive nonlinear crystals [23,26], microcavities [9,[27][28][29], and atomic vapor cells [21], as in the present work. While the main feature of the graphene lattice (the presence of the Dirac cones) is present in all these implementations, other properties can be different. In particular, the SOC in 2D photonic systems such as microcavities is induced by the splitting between TE-TM polarized modes [4]. It is known to modify the dispersion at the Dirac point [7,30], leading to trigonal warping, like in bilayer graphene [31]. In coupled waveguides arrays, usually only a single polarization mode is used and the other polarization can be neglected. In nonlinear crystals and atomic vapor cells, the effects of the SOC on the photonic graphene have not been studied so far.
The evolution of a photonic beam in a spatially varying medium is a particularly important fundamental and applied problem. It is often described in the paraxial approximation of the Helmholtz equation [32], especially in the field of nonlinear optics, where it allows to determine the spatial mode profiles. The coupling of polarizations can arise either due to the anisotropy of the material or to its inhomogeneity [33]. The former usually couples circular polarizations [34] and was already shown to lead to angular momentum transfer [35,36], while the latter has not been fully studied so far. In many cases, the intrinsic coupling of polarizations is simply neglected in the paraxial approximation [33]. Taking it into account in the calculations of the beam trajectory and properties often leads to spectacular effects, such as the spin Hall effect of light [3,37,38].
The behavior of the polarization of light has been described in the limit of geometric optics in the works of Rytov [39]. The Rytov's matrix allows to predict the rotation of the linear polarization. For rays belonging to the class of planar curves (that is, lying in a plane), such rotation is absent: the transverse electric field keeps its polarization in the plane. However, if the ray trajectory becomes three-dimensional (e.g. helix trajectory), the linear polarization starts to rotate. This rotation was linked with anholonomic effects a long time ago [40] and was shown to lead to the accumulation of the Berry phase [41] shortly after its discovery [42]. However, the corresponding theory was limited to ray tracing, equivalent to considering a point-like particle (beam center of mass) instead of a wave packet (beam envelope), whereas modern research subjects, such as the photonic graphene, clearly require a complete wave theory for the transverse beam evolution. The first attempts to develop such theory have shown the necessity for SOC in the paraxial approximation, but did not lead to a self-consistent system of equations [43].
In this work, we introduce the SOC terms into the paraxial equations for the two transverse polarizations. We experimentally demonstrate that these terms play a particular role in photonic graphene, where they couple the spin and angular momentum at the Dirac points, modifying the angular momentum of the probe beam depending on its polarization.

THE MODEL
The Helmholtz equation for the electric field of an electromagnetic wave in a dielectric medium reads In a homogeneous system, this equation does not contain any SOC terms, and the polarizations are decoupled. This allows writing a paraxial equation for a scalar amplitude, corresponding to a single chosen polarization of light (which is conserved). It is very well known that this paraxial equation is equivalent to the time-dependent Schrodinger equation for a wave function of a scalar particle. However, when the spatial gradients are not negligibly small and when the polarization effects are explicitly studied, the SOC has to be taken into account. It is associated with a coupling of the two polarizations, transverse-electric and transversemagnetic [44], which become well-defined in the presence of any gradient (to which they are transverse in addition to being perpendicular to the propagation direction). In geometric optics, the SOC leads to the evolution of the transverse linear polarization along a curved beam [39]. This adiabatic evolution has been shown to lead to dramatic effects, such as the spin Hall effect of light [3,37,45]. Coming back to the analogy with quantum mechanics, the geometric optics corresponds to studying a classical particle, whereas the paraxial equation corresponds to considering an equivalent quantum wave packet. Starting from the evolution of the polarization in the spin Hall effect of light for a beam in the geometric optics limit, that is, for a classical propagating particle, we show how this term is introduced into paraxial equations for the two projections of the electric field amplitude, leading to an original type of SOC. This is especially , Ω is the rotation frequency due to the permittivity gradient ∇ ; b) Beam trajectory in the XZ plane (curved due to the gradient ∇ ). Transverse beam profiles: c) in the XY plane in the E x polarization. Shift along x is due to ∇ ; d) in the XY plane in the E y polarization entirely generated due to SOC, proportional to Ω · k ∼ k y (zero at y = 0).
important for the study of the propagation of beams in such systems as photonic graphene, whose energy bands appear from the quantum-mechanical description. As noted already by Rytov [39], the linear polarization can rotate only for a non-planar beam trajectory. Noting the beam direction as l = k/k and the vector of its rotation as Ω, the relevant quantity describing the part of its rotation which is not in a plane can be written as Ω · l. This is illustrated by a scheme in Fig. 1(a). It was shown that the rotation of the transverse polarization adds an extra term to the Helmholtz equation for the transverse electric field This term introduced in [42] can be understood as an analog of the Coriolis force, appearing in a non-inertial frame, whose noninertial nature is due precisely to the rotation Ω. The geometric optics limit of this equation has allowed to describe the spin Hall effect of light [45].
In order to include this term into the paraxial equations, we need to link Ω with the transverse field E and the dielectric permittivity or the dielectric susceptibility χ. Let us consider a beam with its main propagation direction along z, containing non-zero x and y wave vector harmonics because of its transverse profile, in presence of a gradient along x: ∂ /∂x = 0, which leads to the deviation of the beam from its initial direction: Ω y = (2n) −1 ∂ /∂x. The projection of the rotation frequency on the direction of a particular harmonic is given by In order confirm that the polarization conversion indeed takes place in such conditions, we perform a FDTD numerical simulation of the propagation of a Gaussian beam using COMSOL in the simplest system described above. The results are shown in Fig. 1(b)-(d). The trajectory of the beam, curved by the gradient ∂ /∂x, is shown in Fig. 1(b). The beam is initially excited with E x polarization only. The transverse profile in this polarization is shown in Fig. 1(c). The beam is shifted towards positive x, as in panel (b). Its shape also changes. But the most important effect is shown by Fig. 1(d), presenting the cross-polarization E y , appearing only because of the SOC. This requires a non-zero projection of the wave vector k on its rotation frequency Ω, according to Eq. (3), which is fulfilled thanks to the presence of non-zero k y in the narrow initial Gaussian beam. The intensity of converted polarization is linear in ∇ , confirming the first-order nature of the effect. We note that since the converted signal is proportional to k y , it changes sign at y = 0: the symmetry of the state is inverted (from symmetric to anti-symmetric). This will be important for the understanding of SOC effects in a periodically modulated medium (lattice). Generalizing this result to arbitrary gradient directions and substituting k (x,y) = −i∂/∂(x, y), one obtains the following SOC term: As compared with TE-TM SOC in planar cavities [4,46] and photonic crystal slabs, the double spatial derivative of the electric field is replaced by a product of the first derivatives of the permittivity and the electric field components. In atomic systems, the permittivity can be varied through the effect of electromagnetically induced transparency (EIT) [47]. Another polarization effect which obviously needs to be taken into account in presence of EIT under polarized pumping is the polarization-dependent complex permittivity. This is naturally included in the paraxial equations via x,y or susceptibility χ x,y with real and imaginary parts (marked by and correspondingly). The final paraxial equations taking into account the SOC and the polarization-dependent propagation in the linear polarization basis read: Of course, the importance of the corrective SOC term depends on the conditions of a given experiment. We will compare the predictions of these paraxial equations with the solution of the full system of Maxwell's equations and with the experimental measurements for a particular system of photonic graphene, and show that, because of the polarization-dependent decay rate, the contribution of the SOC actually becomes dominant.

EXPERIMENTAL IMPLEMENTATION OF PHOTONIC GRAPHENE
Whilst the SOC plays a significant role in a large variety of optical systems, we investigate its effect in a particularly in-teresting process: the vortex generation at the Dirac points in photonic graphene. Our experiments are performed in atomic vapors in the regime of EIT. Figure 2 shows the scheme of the experiment for exciting a given valley in the photonic graphene lattice. Three Gaussian coupling beams E 2 , E 2 and E 2 (wavelength λ 2 ≈ 794.97 nm, frequency ω 2 , vertical polarization, Rabi frequencies Ω 2 , Ω 2 and Ω 2 , respectively) from the same external cavity diode laser (EDCL2) symmetrically propagate along the z direction and intersect at the center of the vapor cell with the same angle of 2θ ≈ 1.2 • between any two of them to form a hexagonal optical lattice [see Fig. 2(d)], acting as the coupling field with an intensity of |Ω c | 2 . This lattice appears as a result of interference, and thus does not suffer from any broadening due to diffraction. The 5 cm long atomic vapor cell wrapped with µmetal sheets is heated to 80 • by a heat tape. The co-propagating probe field and coupling field interact with an Λ-type threelevel 85 Rb atomic system [see Fig. 2(b)], which consists of two hyperfine states F = 2 (level |1 ) and F = 3 (|2 ) of the ground state 5S 1/2 , and an excited state 5P 1/2 (|3 ). Here, the probe field With the two-photon resonant condition ∆ 1 − ∆ 2 = 0 satisfied, an EIT window [47] can occur on the transmitted probe spectrum of the probe field. For an EIT configuration, the susceptibility χ = χ + iχ experienced by the probe field is χ ∼ |Ω c | −2 [48], and the resulted susceptibility (both real part χ and imaginary part χ ) exhibits a honeycomb profile with a lattice constant a ≈ 25 µm, which corresponds to an inverted hexagonal |Ω c | 2 [21]. For a three-level EIT atomic system, the real and imaginary parts of the refractive index are described by n R = χ /2 and n I = χ /2, respectively [49,50]. As a result, a photonic graphene lattice governed by n R (x, y) is effectively "written" inside the medium. The refractive index variation and absorption spectrum are polarization-dependent. As presented in Fig. 3 of Ref. [51], for the co-polarized configuration of the coupling and probing beams, the absorption of the probe beam around zero detuning is greatly reduced by the depletion of the ground state with optical pumping, whilst hardly any EIT effect occurs. For the cross-polarized configuration, there is clearly the EIT effect around zero detuning, while the absorption of the probe beam is also relatively strong. As a result, around zero detuning, χ and χ are approximately 5 and 20 times smaller for the co-polarized component (depleted, y) than for the cross-polarized one (EIT, x). The characteristic scales are χ x ∼ 10 −3 , χ x ∼ 10 −4 and χ y ∼ 2 × 10 −4 , χ y ∼ 5 × 10 −6 (for a y-polarized pump, i.e. the coupling beam).
The A-and B-type sublattices of the optically induced twodimensional photonic graphene are marked in Fig. 2(e). By selectively covering only the A or the B sublattice with the periodic probe field as in Fig. 2(c), the K or K valley in the momentum space can be effectively excited, which, due to the beam conversion between the sublattices, introduces an orbital angular momentum (OAM) to one of the output probe beams [21,23], confirmed by a fork-like feature in the interference pattern with a Gaussian reference beam (from the same laser as probe beams). Both the transmitted probe beam and phase (interference pattern) are monitored by a charge coupled device camera (CCD). To investigate the influence of the polarization . The EIT spectrum is generated by injecting beams E 1 and E 2 (from the same laser as E 1 and E 2 , respectively) into the Rb cell 2 and received by an APD.
of the probe beams on the OAM creation from the valley pseudospin, a quarter-wave plate is applied on the path of the probe beams before entering the atomic vapor cell, allowing varying the probe polarization. It should be noted that regardless of the probe beam polarization, the CCD always detects a single linear polarization component (x) defined by the PBS before it.

RESULTS AND DISCUSSION
We are now considering polarized light described by a two component wave-function propagating in the photonic graphene lattice. In general, solving 2D paraxial equations is much more efficient than solving the full system of Maxwell's equations in 3D. We begin our analysis with the comparison of the results of the two different numerical approaches for the simple case of a Gaussian probe of a width w exciting the vicinity of the Dirac point of graphene, and neglecting the photonic SOC. In such a case the evolution of the beam is well described by the 2D Dirac equation describing the coupling between A and B sites of the honeycomb lattice. This approximation remains valid even including SOC in both excitation and detection are carried out in the same linear polarization. The simulated probe is exciting only A-sites of the graphene lattice, and its conversion to the B-sites is accompanied by the change of angular momentum according to l B = l A − 1. This conversion is due to the non-zero Berry curvature of the graphene bands in the vicinity of the Dirac point [52], as shown by previous studies [21,23]. This behavior can be obtained both with the 3D FEM and 2D paraxial equations, as shown in Fig. 7(a,c), and, naturally, with the Dirac equation (which is an extra approximation with respect to the paraxial equation). The interference patterns in panels (a,c) show the forklike dislocation which is a signature of a non-zero final angular momentum: l B = −1 (initially, l A = 0). We note that the A → B conversion period T = w/c depends on χ via the effective speed of light in the Dirac equation c.
The situation changes when the excitation contains both E x and E y (circular pumping) and the SOC starts to play a role. Figure 7(b,d) shows interference patterns between the transmitted amplitude in the x-polarization and a reference beam. It does not show any dislocations in the fringe patterns. The absence of angular momentum conversion is a joint result of 3 different effects. First, the y-polarized component exhibits a slower A → B conversion because of a smaller real susceptibility χ y ≈ 0.2χ x . So E y does not change angular momentum during the propagation time in the cell. The second effect is the smaller decay of the y-component, because of a smaller imaginary susceptibility χ y ≈ 0.05χ x . As a result, the E y component quite rapidly becomes dominant over E x and most of the x-polarized light is induced by the E y to E x conversion by the SOC. Ultimately, the E x generated by the SOC is not affected by the angular momentum conversion because the corresponding wave function is not anymore close to the Dirac point but lies in an upper band of graphene as explained below. As a result of these combined processes, the angular momentum L = 0 measured in E x in Fig. 7(b,d) is that of l A injected in E y and we can conclude that the angular momentum of the light after the cell is controlled by the incoming polarization.
Let us first make a simple estimate proving that the contribution of the polarization conversion can indeed become dominant with respect to the rapidly decaying original signal. The relevant terms in the paraxial equation for the description of the intensity contributions to the detected polarization E x become The first term provides an exponential decay with a characteristic rate k 0 χ x ≈ 10 7 × 10 −4 = 10 3 m −1 , while the second one provides a linear growth (assuming E y = const as compared with E x ). To estimate its rate, we use the experimental parameters k x /k 0 ≈ θ ≈ 10 −2 and ∂χ y /∂x ≈ 10 −4 /25 × 10 −6 = 4, which gives a rate of the order of 10 −2 . Comparing an exponential decay exp −10 3 z with a linear growth 10 −2 z for identical initial intensity (circular pumping), we see that the two contributions become equal within a propagation distance of 1 cm. We can therefore conclude that for a 5 cm vapor cell the contribution of the SOC-induced polarization conversion can indeed be important and even dominant. A more detailed discussion of the evolution of the intensities based on rate equations is given in the Supplemental materials. The other crucial effect is the absence of angular momentum conversion for the E x polarization induced by SOC. This can be understood by using arguments based on tight-binding description (while the development of a complete tight-binding model accounting for the specific SOC is beyond the scope of the present work). The SOC terms contain the first-order derivatives of the wave function. It thus inverts the mode symmetry, coupling the s (symmetric) and p (anti-symmetric) orbitals of each minimum of the effective lattice potential. Indeed, in Fig. 1 a single maximum (Gaussian, s-orbital) was converted into an anti-symmetric double-maximum structure (p-orbital). The converted polarization therefore appears at the same valley, but in the p band. The shape of the wave-function suggests that conversion occurs toward the flat p band of the honeycomb potential [27] where further evolution is completely blocked (see the Supplemental materials for more details). This is why the converted component maintains its original angular momentum. The observed signal is the superposition of the original component E x (changing the angular momentum) and the signal converted from E y (keeping the original angular momentum). The final result depends on the relative intensities of the two components, and varying their initial ratio (which is the circular polarization degree) allows to observe the vortex leaving the beam center step by step, when the circular polarization is increased (see below for the experimental results), similar to the spatial dynamics observed in Ref. [21], but in opposite direction.
The comparison of the 3D and 2D results demonstrates the correctness of the developed corrected paraxial equations. In what follows we are are going to use only the 2D model and to consider a probe with non-zero initial angular momentum. Figure 8 shows the results of numerical simulations with 2D coupled paraxial equations for the case of Gauss-Laguerre probe envelope with L = −1 (a,b) and L = +1 (c,d) with linear (a,c) and circular (b,d) probe polarization. For L = −1 and linear excitation, the beam conversion introduces an extra vortex into the pattern, giving L out = −2. For L = +1, the sign of the extra vortex is opposite to the initial angular momentum, and L out = 0. This extra vortex disappears in both cases if the initial polarization is circular, because the output beam is dominated by the E y converted to E x . A crucial advantage of using a 2D model here is that it allows a more extended treatment of the output field distributions, in particular, valley separation and interference analysis (based on Fourier transform and shifting in the reciprocal space), similar to the beam separation after the vapor cell in the experiment [21]. . Numbers correspond to circular polarization degree: from ρ c = 0 (linear probe, a/b/c1) to ρ c = 1 (circular probe, a/b/c5). In each case, the change of polarization induces the decrease of the output beam angular momentum by 1.
Finally, in Fig. 5 we present the results of experimental measurements carried out in the configuration described in the previous section, for 3 values of the angular momentum of the inbound Gauss-Laguerre probe: a) L = 0 (as in Fig. 7, b) L = −1, c) L = +1 (as in Fig. 8). The false color scale represents the intensity of the interference pattern, measured at the screen after the vapor cell in the x polarization. The interference occurs between one of the two probe beams and a reference beam. The numbered panels represent the increase of the circular polarization degree from ρ c = 0 (a1,b1,c1) to ρ c = 1 (a5,b5,c5). The experimental observations correspond to the theoretical images shown above: in each case, the extra vortex brought into the beam by the interaction with graphene lattice disappears when ρ c is increased, because signal is dominated by the converted component.
In all three cases we observe that the SOC allows controlling the angular momentum of the output beam via both the polarization and momentum of the input beam. Devices with an optically-controllable optical angular momentum are useful for practical applications, but creating them is a challenging task [53]. Our experiment provides a solution of this problem, while demonstrating the fundamental importance of SOC in the paraxial approximation, which has been neglected so far. In the future works, we will develop the tight-binding description of this original type of SOC.
To conclude, we introduce a first-order correction to paraxial equations due to the SOC linked with the conservation of the polarization plane. We demonstrate experimental evidence of the importance of this SOC in photonic graphene implemented via EIT in atomic vapors. The angular momentum of the output beam can be controlled by the circular polarization of the probe. Disclosures. The authors declare no conflicts of interest.

SUPPLEMENTAL DOCUMENTS
See Supplement 1 for supporting content.

SUPPLEMENTAL MATERIALS
In this supplemental material, we present additional results concerning the polarization conversion rate, the mixing of the s and p bands by the spin-orbit coupling, and the sensitivity of the qualitative results to the input parameters.

A. Rate equations for polarization conversion
In the main text of the manuscript, we discuss the relative contribution of the original and converted polarization components to the detected intensity. Here, we use the rate equations equivalent to the paraxial equations in order to describe the behavior of the intensity and to confirm that ultimately all measured intensity originates from converted polarization (under circular pumping).
The rate equations for the intensities of the x and y polarizations I x (t), I y (t) can be written as: This system of differential equations can be solved analytically, but in order to have less cumbersome expressions we shall apply several approximations. First, since Γ x Γ y , we can neglect Γ y . Second approximation is based on the fact that the polarization conversion rate x → y is much smaller than the decay of the x polarization: W Γ x (see the main text for the discussion of the numerical values). This allows to neglect the second-order term using W/Γ x 1. Keeping only the first order terms in W/Γ x , the solution for the detected component I x (t) can be written as: The physical meaning of this solution is clear: the first term corresponds to the original polarization, which decays even faster than just Γ x because of the additional losses due to the conversion to the other component. The second term corresponds to the intensity generated from the second component, and its decay rate is determined by W.
To get the information on the fraction of the intensity I x coming from the conversion from the other component I y , we can compare the solution Eq. (10) to the behavior of an isolated polarization component I i with the same decay rate Γ x , given by I i (t) = I 0 exp(−Γ x t). The expression for this fraction reads This expression can be simplified again by using W Γ x . At small times, the series expansion of the exponents gives and at longer times We conclude that in the regime of interest, the fraction of the converted component grows linearly with time, with a rate directly determined by the efficiency of the polarization conversion.

B. Band mixing by spin-orbit coupling
In this section, we explain the microscopic mechanism of the mixing of the s and p bands of the photonic graphene by the spin-orbit coupling appearing due to the susceptibility gradient.
As an illustration, we perform a 3D FEM simulation of the beam evolution in the photonic graphene lattice using COMSOL, as in Fig. 3(c)   the excitation of the system with orthogonal linear polarization E y , in order to elucidate clearly how it evolves in the lattice. Moreover, we simulate a 3-beam excitation instead of a 2-beam one, in order to be able to make conclusions on the band mixing by studying the intensity and phase distribution. Figure 6(a) shows the distribution of intensity E 2 y at the input of the vapor cell. The susceptibility profile responsible for the effective potential of the honeycomb lattice is shown with white contour lines. The maxima of intensity appear at the positions of the sites of a single type only (the A-sites), as expected for a state of the s-band. Figure 6(b) shows the intensity in the x-polarisation at the output of the cell. The maxima are now located at the center of each unit cell, corresponding to a state of the p-band, as expected from the symmetry of the spin-orbit coupling. To explain the phase distribution in the final state, we present a scheme showing the microscopic mechanism of its formation in Fig. 6(c). The wavefunction of the K point of the s-band is initially distributed over the A-sites with a phase forming a vortex around the center of the unit cell (phase values shown at each site). Because of the effective potential due to the susceptibility gradient, the intensity expands in the directions shown by the arrows, but the phase of the generated wave is opposite at each side of each arrow (anti-symmetric function).
The final values of the phase in each sector are shown with colored numbers. They correspond to a vortex around the empty B-site. This is confirmed by the phase arg E x extracted from the numerical simulations, shown in Fig. 6(d). This Bloch function corresponds to the same valley as the one of the s-band, but now it belongs to the p-band, because the particles are located at the barriers of the potential, and not at its minima. The shape of the orbital corresponds to the frustrated states from the flat p-band of photonic graphene studied experimentally and theoretically in Ref. [27].
We note that the vortices forming in the Bloch function at each unit cell should not be confused with the vortex forming in the envelope function! The spatial separation of the beams forming the Bloch function after the vapor cell allows to get access to the envelope function. In simulations, this is done using Fourier transform and k-space shifting. We have tested this configuration (with only E y excitation) experimentally as well, and these experiments confirm the blocking of the formation of the vortex in the envelope function, which we explain by the coupling with the flat p-band.

C. Sensitivity to parameters
The susceptibility of the atomic vapors under the effect of optical pumping is not known with a very high accuracy. This especially concerns its polarization dependence. While the qualitative ratio χ /χ ∼ 0.1, χ y /χ x ∼ 0.1 is generally accepted, the exact ratios are difficult to determine and depend on the detuning of a particular experiment. The results shown in the main text correspond to χ x = 5χ y and χ x = 20χ y . In this section, we show that qualitatively the same results are obtained if these ratios are taken equal to 10 and 10, respectively. Figure 7 shows the comparison of the 3D (FEM beam envelope, (a) and (b)) and 2D (paraxial, (c) and (d)) models for a Gaussian incident beam (L = 0) with linear ((a) and (c)) and circular ((b) and (d)) polarization. While the distribution of intensity is slightly modified with respect to the main text, the qualitative behavior is conserved: the generated vortex observed in the envelope function for the linear probe disappears for the circular probe. Figure 8 shows the results obtained by numerical simulations based on the 2D paraxial equation for the Gauss-Laguerre envelope beams L = ±1. In both cases, the angular momentum changes by 1 under linear pumping (negligible SOC), while this change does not occur under circular pumping (SOC-dominated regime). The behavior is therefore the same as in the main text.