For Table of Contents Only

# Extreme Parametric Sensitivity in the Steady-State Photoisomerization of
Two-Dimensional Model Rhodopsin

Chern Chuang<EMAIL_ADDRESS>Paul Brumer<EMAIL_ADDRESS>Chemical Physics Theory Group, Department of Chemistry, and Center for Quantum
Information and Quantum Control, University of Toronto, Toronto, Ontario M5S
3H6, Canada

###### Abstract

We computationally studied the photoisomerization reaction of the retinal
chromophore in rhodopsin using a two-state two-mode model coupled to thermal
baths. Reaction quantum yields at the steady state (10 ps and beyond) were
found to be considerably different than their transient values, suggesting a
weak correlation between transient and steady-state dynamics in these systems.
Significantly, the steady-state quantum yield was highly sensitive to minute
changes in system parameters, while transient dynamics was nearly unaffected.
Correlation of such sensitivity with standard level spacing statistics of the
nonadiabatic vibronic system suggests a possible origin in quantum chaos. The
significance of this observation of quantum yield parametric sensitivity in
biological models of vision has profound conceptual and fundamental
implications.

Retinal proteins are a family of important membrane proteins responsible for
various photobiological functions in bacteria and higher life forms.1, 2, 3
Despite the variety of their functions and the host organisms, the general
structure of the retinal proteins is consistent, i.e. seven transmembrane
helices with a retinal chromophore bound to a lysine residue through a
protonated Schiff base. Two systems are frequently studied: (1) rhodopsin that
is found in animal retina, responsible for dim light vision, and (2)
bacteriorhodopsin found in certain bacteria which functions as a light-driven
proton pump. In particular, rhodopsin, with its retinal chromophore initially
in an 11-cis conformation, undergoes photoisomerization to its all-trans
counterpart within 200 fs upon pulsed laser excitation.4, 5, 6 This has been
recently revised to 60 fs. 7 The short reaction time is attributed to a
conical intersection between the ground and excited states, and has been
suggested as being responsible for the high photoisomerization quantum yield
($\approx$67%).8 These observations have prompted the important question of
the correlation between quantum coherence and biological function. Our view is
clear – such coherences are induced by the laser pulses, and are irrelevant to
nature, which operates with incoherent light.9, 10 Indeed, the natural vision
process operates in the steady state regime.

Here, we report the observation of a dramatic dependence of the steady-state
quantum yield (QY) on minute changes in system parameters. Our study is based
on a two-state two-mode (2S2M) model coupled to thermal baths initially
proposed by Stock et al.11, 12, and later modified/extended to account for
numerous experimental observations of bovine rhodopsin photoisomerization.7,
13 The dependence of the biologically important steady-state QY on the system
parameters is correlated with level spacing statistics, a well established
measure for quantum chaos in molecular systems. This parametric sensitivity
has profound implications for the biology of vision, raising issues as to
whether nature’s parameters occupy a “Goldilocks” region for optimal function.

This letter is arranged as follows. We first introduce the necessary
background and details of the theoretical treatment, which is followed by a
demonstration of the dramatic parametric sensitivity of the steady-state
photoisomerization QY. We then discuss relationships to quantum chaos, the
possibility of experimental measurement of such extreme sensitivity, and its
implications on generic photochemistry in the condensed phase and biological
light-sensing.

Models and Method— The model for the retinal rhodopsin photoisomerization
reaction is simulated numerically using the 2S2M model pioneered by Stock et
al.11, 12, 14. This has been shown to be an adequate minimal model for
studying the photoisomerization dynamics of retinal chromophores, describing
the vibronic structure of a conical intersection.15 To account for the
environmental effects associated with the protein pocket, the 2S2M system is
further coupled to two sets of harmonic bath modes, and the system dynamics is
simulated using second order quantum master equations.

A general system-bath Hamiltonian can be written as

$\displaystyle H=H_{\mathrm{s}}+H_{\mathrm{b}}+H_{\mathrm{sb}}.$ (1)

Here, the system (2S2M) Hamiltonian is comprised of two diabatic electronic
states $|0\rangle$ and $|1\rangle$, a reaction coordinate (coupling mode)
$\phi$, and a tuning mode $x$.

$\displaystyle H_{\mathrm{s}}$ $\displaystyle=$
$\displaystyle\sum_{n,n^{\prime}=0,1}\left[\left(\hat{T}+E_{n}+(-1)^{n}\frac{V_{n}}{2}(1-\cos\phi)+\frac{\omega
x^{2}}{2}+\kappa x\delta_{n,1}\right)\delta_{n,n^{\prime}}+\lambda
x(1-\delta_{n,n^{\prime}})\right]|n\rangle\langle n^{\prime}|$ (2)

where
$T=-(1/(2m))(\partial^{2}/\partial\phi^{2})+(\omega/2)(\partial^{2}/\partial
x^{2})$ is the kinetic energy operator. All parameters, updated in 2017, are
obtained from Ref.[7], and given in Table 1. We define the range
$\phi\in\left[-\pi/2,\pi/2\right)$ to be the cis-conformer, and the range
$\phi\in\left[\pi/2,3\pi/2\right)$ to be the trans-conformer. The eigenstates
of $H_{\mathrm{s}}$ can be expressed in the direct product space of the
electronic diabatic states $\\{|0\rangle,|1\rangle\\}$, the tuning mode $x$,
and the reaction coordinate $\phi$.11 The adiabatic potential energy surfaces
are plotted in Fig. 1(a). We note that the protein pocket of rhodopsin is such
that the retinal chromophore is pre-twisted and chiral, which Eq. (2) does not
account for.16

A generic property of such a Hamiltonian is that the low energy part of its
eigenstates is localized in either the cis- or the trans-wells, while higher
energy states are more delocalized over the reaction coordinate $\phi$.17 To
better demonstrate this character, of relevance later below, we define the
following quantity

$\displaystyle l_{k}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}d\phi~{}|\langle\phi|k\rangle|^{2}\cdot\frac{(1-\cos\phi)}{2}$
(3)

where $|k\rangle$ is the $k^{th}$ eigenfunction of the system Hamiltonian
$H_{\mathrm{s}}$, and $|\phi\rangle$ is the eigenfunction of the position
operator in the reaction coordinate. This quantity is a non-negative real
number with a maximal value of unity. It measures the “trans-ness” of the
eigenfunctions: $l_{k}=0$ $(1)$ corresponds to $|k\rangle$ fully localized in
the cis- (trans-) well. For intermediate values it indicates an averaged
position of the eigenstate, localized or delocalized. Fig. 1(b) shows the
distribution of $l_{k}$ and eigenstate energies $\epsilon_{k}=\langle
k|H_{\mathrm{s}}|k\rangle$. Contribution from the trans species occurs at most
energies, but becomes significant above $\sim 1.4$ eV.

Figure 1:  (a) Adiabatic potential energy surfaces of the 2S2M model. (b) The distribution of the system eigenstates with their energies as the vertical coordinates and $l_{k}$, defined in Eq. (3) as the horizontal coordinates. The state that contributes most to the Franck-Condon state is marked with an asterisk. See text for detailed description. Table 1: List of parameters, adopted from Ref. [7 ] Parameter | Value (eV)
---|---
$m^{-1}$ | 2.80$\cdot 10^{-3}$
$E_{0}$ | 0
$E_{1}$ | 2.58
$V_{0}$ | 3.56
$V_{1}$ | 1.19
$w$ | 0.19
$\kappa$ | 0.19
$\lambda$ | 0.19

System-Bath Coupling: To properly account for relaxation effects, the system
is coupled to environmental degrees of freedom, approximated by the direct
product form

$\displaystyle H_{\mathrm{sb}}$ $\displaystyle=$
$\displaystyle\sum_{\mathrm{b}}\hat{S}^{\mathrm{(b)}}\otimes\hat{B}^{\mathrm{(b)}},$
(4)

where $\hat{S}^{\mathrm{(b)}}$ and $\hat{B}^{\mathrm{(b)}}$ are operators that
depend only on the system and the bath coordinates, respectively. The
superscript $\mathrm{b}$ runs over all system-bath interactions involved. This
form of system-bath interaction Hamiltonian is then treated with the standard
Markovian-Redfield approximation that leads to a quantum master equation,
which in the system eigenbasis reads

$\displaystyle\frac{\partial}{\partial t}\rho_{ij}$ $\displaystyle=$
$\displaystyle-i\omega_{ij}\rho_{ij}+\sum_{\mathrm{b}}\sum_{kl}\mathcal{R}_{ij,kl}^{\mathrm{(b)}}\rho_{kl}$
(5) $\displaystyle\mathcal{R}_{ij,kl}^{\mathrm{(b)}}$ $\displaystyle=$
$\displaystyle-\delta_{jl}\sum_{r}\Gamma_{irrk}^{\mathrm{(b)}}+\Gamma_{ljik}^{\mathrm{(b)}}+\left(\Gamma_{ljik}^{\mathrm{(b)}}\right)^{*}-\delta_{ik}\sum_{r}\left(\Gamma_{lrrj}^{\mathrm{(b)}}\right)^{*}$
(6)

where $\omega_{ij}$ is the energy difference between the system eigenstates
$i$ and $j$, and $\hbar=1$. The damping tensors $\Gamma_{ijkl}^{\mathrm{(b)}}$
are given by

$\displaystyle\Gamma_{ijkl}^{\mathrm{(b)}}=S_{ij}^{\mathrm{(b)}}S_{kl}^{\mathrm{(b)}}J^{\mathrm{(b)}}(\omega_{kl})\bar{n}_{\mathrm{BE}}(\omega_{kl},\beta),$
(7)

where $S_{ij}^{\mathrm{(b)}}=\langle i|\hat{S}^{\mathrm{(b)}}|j\rangle$,
$\bar{n}_{\mathrm{BE}}(\omega,\beta)=1/(e^{\beta\omega}-1)$ is the Bose-
Einstein distribution at inverse temperature $\beta=1/k_{\mathrm{B}}T$ with
the Boltzmann constant $k_{\mathrm{B}}$.

The computational cost of propagating the Redfield master equation [Eq. (6)]
can be formidable given that one typically needs up to a thousand lowest
energy eigenstates, enough to converge the Franck-Condon transition from the
ground to excited state. Significantly, this is made more difficult by the
fact that we need to examine the long time steady-state properties, which
requires propagating the dynamics up to a nanosecond and beyond. To this end,
one can invoke the Bloch-secular approximation which decouples the population
dynamics from that of the coherences.14 We compare the results from the full
nonsecular treatment to those of the secular approximation in the next
section. (For issues of secular v.s. nonsecular master equations, see
Ref.[18].)

The phonon bath is described within the harmonic approximation,
$H_{\mathrm{ph}}=\sum_{k}\omega_{k}\cdot b^{\dagger}_{k}b_{k}$, where the
system-phonon interaction can be further separated into two components:

$\displaystyle H_{\mathrm{s-ph}}$ $\displaystyle=$ $\displaystyle
H_{\mathrm{s-}x}+H_{\mathrm{s}-\phi}$ (8) $\displaystyle H_{\mathrm{s-}x}$
$\displaystyle=$ $\displaystyle|1\rangle\langle
1|~{}x\cdot\sum_{k}g_{k,x}(b_{k,x}^{\dagger}+b_{k,x})$ (9) $\displaystyle
H_{\mathrm{s-}\phi}$ $\displaystyle=$ $\displaystyle|1\rangle\langle
1|(1-\cos\phi)\cdot\sum_{k}g_{k,\phi}(b_{k,\phi}^{\dagger}+b_{k,\phi})$ (10)

where $g_{k,x}$ and $g_{k,\phi}$ are the coupling strengths of the phonon mode
$k$. It is customary to represent the summation over the phonon coupling
strengths $g_{k}$ as a continuous spectra density $J(\omega)$, taken here in
accord with the work of Stock et al.,12, 14 to be an Ohmic bath with
exponential cut-off, $J(\omega)=\eta\omega e^{-\omega/\omega_{c}}$. Couplings
to these baths represent the fluctuations and dissipation of the system
dynamics along the modes $x$ and $\phi$.

In the simulation of the retinal photoisomerization, it is customary to choose
the initial condition of the dynamics to be the Franck-Condon state from the
ground state located in the cis-well. 11, 7, 19 On the other hand, the
dynamics initiated from a state similar to the Franck-Condon state but with
all coherence in the energy eigenbasis removed, referred to as the Franck-
Condon mixed state, mimics the effect of an incoherent light source such as
sunlight.17 We focus on the dynamics initiated with the Franck-Condon states
and compare the effect of incoherent Franck-Condon mixed initial state in the
Supporting Information. Similar parametric sensitivity was also observed in
the non-equilibrium steady state case, where the retinal system is
simultaneously connected to a photon bath describing the incoherent sunlight
and a phonon bath describing the protein environment. (See Supporting
Information.)

Variation of System Parameters— As noted above, the main result of this work
is our observation of the extreme sensitivity of the steady-state
photoisomerization QY to parameter variation. To demonstrate the extreme
sensitivity, consider the results of varying two representative parameters, as
manifest in the computed steady-state QY. First, we vary the inverse mass
$m^{-1}$ in the kinetic energy term of the reaction coordinate. In the 2S2M
model, $m^{-1}$ determines the characteristic frequency of the reaction
coordinate in the harmonic regime, $\omega_{\phi}=\sqrt{k/m}$ and $k=V_{0}/2$.
It is worth noting that in the original model of Hahn and Stock11
$m^{-1}=4.84\cdot 10^{-4}$ (eV), nearly six-fold smaller than the value
recently suggested by Johnson et al.7, introduced to reflect the shortening of
the rise time of the $trans-$photoproduct observed in pump-probe and 2D photon
echo measurements with better temporal resolution. Here we consider $\pm 5\%$
variation of $m^{-1}$ about the value proposed by Johnson et al. Note that
this corresponds a small change, i.e. $~{}\pm 2.5\%$, of the harmonic
frequency along the reaction coordinate.

We have also examined the effect of varying the optical gap between the two
diabatic electronic states, namely the parameter $E_{1}$. It is estimated that
the inhomogeneous contribution to the broad absorption band of retinal
rhodopsin is $\sim$1000 cm-120, $\sim 5\%$ of the optical gap, and which is
the range over which we vary the optical gap. To avoid the complication
associated with changing the energy storage ($E_{1}-V_{1}$) and the optical
gap of the $trans-$product ($V_{0}-E_{1}+V_{1}$) when varying $E_{1}$ we
simultaneously change $V_{1}$ to keep $E_{1}+V_{1}$ constant.

Variation of $m^{-1}$: Fig. 2(a) shows the time-dependent $trans$-ground state
population of the photoisomerization dynamics, which gives the “time-dependent
quantum yield”, defined as the population of diabatic electronic state
$|1\rangle$ in the $trans-$region ($\pi/2\leq\phi<3\pi/2$). The dynamics of
the original model with parameters given in Table. 1 is shown in black, and
those with modified $m^{-1}$ parameters are color-coded. The bath temperature
is kept at 0K to avoid the complication of thermally activated isomerization.
Test computations with the bath at room temperature were also carried out to
ensure that the parameter dependences reported below remained qualitatively
unchanged.

Note first that there are signs of coherent beatings at times less than 1 ps
reflecting the fact that the population is oscillating between the reactant
and the product wells, as the population in the cis-well mirrors that of the
trans- (data not shown). The photoproduct appears on the time scale of 100 fs
alongside the beatings, in general agreement with ultrafast pulsed laser
experiments,4, 5, 6, 7 signaling coherent wavepacket dynamics along the
direction of the reaction coordinate $\phi$. While the beatings disappear on
the ps time scale, the QY continues to fluctuate on a far longer time scale,
up to 100 ps. This is consistent with the experimental observations of a
closely related system, the photoisomerization of bacteriorhodopsin, which
shows a similar time-dependent photoproduct signal (up to hundreds of ps)
after the pump pulse.21 Computationally, the photoproduct signal has a non-
monotonic time dependence, a steep rise on the time scale of 1 ps, followed by
a smooth decay until 20 ps, reaching its stationary value at longer times.

Here we are varying the particular system parameter $m^{-1}$ within the range
of $\pm 5\%$ with respect to the original value. The transient dynamics, up to
$\sim$0.5 ps, show limited $m^{-1}$ dependence. However, most significantly,
the QY at longer times (and hence not observable in transient dynamics
reported in most pulsed laser experiments) depends very sensitively on
$m^{-1}$, spanning a wide range between $35\%$ and $75\%$. This is clearly
seen in the QY recorded at 10 ns as a function of $m^{-1}$ shown in Fig. 2(b).
Note also that Fig. 2 shows both secular and nonsecular simulations. In
general, these two sets of results agree qualitatively with QY sensitive to
the value of $m^{-1}$. However, it is clear that the steady-state QY depends
on the level of treatment for the system-bath interaction.

Figure 2: (a) Time dependent $trans$-photoproduct population. The full
nonsecular dynamics is shown as solid lines and those of the Bloch-secular are
shown with dashed lines. The data corresponding to the original set of
parameters given in Table 1 are shown in black. Data with varied $m^{-1}$ are
color coded as indicated by the color bar. (b) QY recorded at 10 ns. The blue
and red shaded area correspond to the connectivity graphs shown in Fig. 1,
discussed later below. The bath parameters are adopted from Ref.[14].

Variation of The Optical Gap: To corroborate the above result and ensure that
the parameter $m^{-1}$ is not particularly special with regard to the
sensitivity of steady-state QY, we also examined the system dependence on the
value of the optical gap. This quantity is known to be inhomogeneously
broadened in natural systems and easier to measure in concert with the
photoisomerization reaction itself. In our case, this amounts to changing the
parameter $E_{1}$ in the 2S2M model as described in the Models and Method
section. The results are shown in Fig. 3, and are similar to those found in
variation of $m^{-1}$, i.e. large variations in QY as a function of small
changes in the optical gap.

Figure 3: As in Fig. 2, but with $E_{1}$ is changed instead of $m^{-1}$.

We also examined other parameters such as $\omega$ (the frequency of the
harmonic, second mode in the 2S2M model) and $\lambda$ (the magnitude of the
nonadiabatic coupling, similar to that studied in Ref.[22]) and, moreover,
two-state-one-mode and two-state-two-mode models with avoided crossings
instead of a conical intersection (parameters given in Ref.[14]). In light of
recent computational study,13, 23 we simulated and found similar pararmetric
sensitivity in a two-state-three-mode model as well. In addition to the above
mentioned nonadiabatic, two-surface models, we consider a single-surface
double well (1S2M) model as well that is more suitable for describing the
adiabatic regime (strong nonadiabatic coupling) of the isomerization reaction.
Details are provided in the Supporting Information. In all cases similar
extreme sensitivities of steady-state QY are found with respect to variations
of system parameters. This suggests that the phenomenon may be an inherent
property of the multi-dimensional vibronic/vibrational systems.24, 25

System Eigenstate Statistics and Relaxation Pathways— We conjecture that the
parametric sensitivity of the steady-state QY results from the sensitivity of
chaotic system dynamics that moves between the two wells and eventually
settles into the product well by relaxation induced by the system-bath
coupling. As mentioned in the introduction, previous studies on similar
systems in the quantum chaos literature focus mainly on the statistical
properties of their eigenstates24 or short-time trajectory-based simulations
in the absence of thermal relaxation.26, 27 The system part in the current
study, Eq. (2), falls under the same category of nonadiabatic vibronic
Hamiltonians.24, 22 It is thus expected, if the open system character is not
significant, that similar quantum chaos characteristics would exist in this
system as well. Fig. 1(b) shows the distribution of eigenstates according to
their energies and $l_{k}$, a measure of the average position along the
reaction coordinate. It is clear that in the low energy and deeply
$cis$/$trans$ regions ($l_{k}\rightarrow 0$ or $l_{k}\rightarrow 1$), the
distribution of eigenstates forms a regular grid point lattice. This implies
that these states can be written as direct products of the eigenstates of the
two modes, $|k\rangle\propto|n_{\phi}\rangle\otimes|n_{x}\rangle$, where
$n_{\phi}$ and $n_{x}$ are the quantum numbers of the respective modes. On the
other hand, for states with higher energies and $l_{k}$ closer to 0.5, the
grid pattern disappears and the characters of the two diabatic electronic
states and the two modes mix.

The nearest-neighbor spacing distribution (NNSD) of energy eigenstates is a
commonly utilized measure to characterize possible chaotic nature of an
isolated system. 26, 22 Specifically, level repulsion tends to imply chaos,
most clearly manifest in a Wigner distribution of adjacent energy level
spacings. Here, the procedure described by Haller et al. is used to account
for the smoothly varying part of the density of states.26 The “unfolded”
energy levels are defined as

$\displaystyle\tilde{\epsilon}_{i+1}=\tilde{\epsilon}_{i}+(2k+1)\frac{\epsilon_{i+1}-\epsilon_{i}}{\epsilon_{j_{2}+1}-\epsilon_{j_{1}}}$
(11)

where $\epsilon_{i}$ are the energies of the eigenstates of $H_{\mathrm{s}}$
with the state index $i$ that runs from $1$ to $n$, $k$ is a small positive
integer defining the locality of the unfolding, $j_{1}=\mathrm{max}(1,i-k)$
and $j_{2}=\mathrm{min}(n-1,i+k)$. NNSD is the distribution function of
$S=\tilde{\epsilon}_{i+1}-\tilde{\epsilon}_{i}$. A non-monotonic NNSD signals
quantum chaotic behavior. Most notably, by assuming linear repulsion between
adjacent levels one arrives at the Wigner distribution for chaotic systems:

$\displaystyle P_{\mathrm{w}}(S)=\frac{\pi S}{2D^{2}}e^{-\frac{\pi
S^{2}}{4D^{2}}}$ (12)

where $D$ is the mean of the distribution.

The NNSD for our system is shown in Fig. 1(a). It is clear that the
distribution is non-monotonic and near Wignerian and hence chaotic. This is
especially true for the higher energy part of the eigenspectrum (blue), in
comparison to the lower energy counterpart (brown).

Figure 4: (a) The nearest-neighbor energy level spacing distribution of the
model system studied, using parameters in Table 1. The blue bars (higher) are
the 1240 eigenstates taken from the energy range 4 to 6.4 eV, and the brown
ones (lower) are the 686 eigenstates with energy below 4 eV. The local mean
spacing parameter is set to $k=2$. The average spacing of the unfolded levels
is $D=1.00$ meV. The curve is the Wigner distribution, given in Eq.
(12).111The adoption of the higher energy range in Fig. 4(a) is for the
purpose of emphasizing the contribution from the chaotic regime. In all other
instances in the paper the photoisomerization dynamics is well converged by
including energy eigenstates below 3 eV. (b) The connectivity tree graph of
the systems with $m^{-1}$ values 99.70% and 99.80% of that in Table. 1 [red
shaded region in Fig. 2(b)], rooted at the brightest state (asterisk). (c) The
tree graph of the systems with $m^{-1}$ 101.3% and 101.4% of that in Table. 1
[blue shaded region in Fig. 2(b)].

Quantum chaos in the context of vibronic eigenspectrum of small molecules/ions
in the gas phase has been well established. However, the situation is less
clear with systems, such as ours, coupled to a dissipative environment. As one
insightful route, we analyze the relaxation pathways from the brightest state
(the eigenstate contributing the most to the Franck-Condon state) by examining
the Bloch-secular rate matrix elements among the eigenstates.

In Fig. 1(b) and (c) we show the tree graphs representing the connectivity
matrix28 of the relaxation pathways among system energy eigenstates calculated
using two pairs of $m^{-1}$ values differed by only $0.1\%$, corresponding to
the red and blue shaded regions shown in Fig. 2(b), while keeping all other
parameters fixed. Each graph shows the eigenstate distributions [as in Fig.
1(b)] of the two system Hamiltonians with slightly different $m^{-1}$ values,
using different markers (black filled circles and colored stars). As shown in
Fig. 2(b), the red region corresponds to a parameter regime that shows smooth
changes of steady-state QY with changing $m^{-1}$ values, whereas the blue
region shows dramatic changes. To simplify the graphs we only connect each
higher energy state $i$ to the two lower energy states $j$ with the largest
scattering rates $\sum_{\mathrm{b}}\mathcal{R}^{\mathrm{(b)}}_{jjii}$ in Eq.
(6). The connectivity graphs are rooted at the brightest state of the
corresponding systems, defining tree graphs of degree two.

In both pairs of the tree graphs it can be seen that very small changes in
$m^{-1}$ leaves the approximated direct product structure for states close to
the cis\- ($l_{k}\approx 0$) and the trans-wells ($l_{k}\approx 1$) nearly
unaffected. That is, a 0.1% change of $m^{-1}$ corresponds to a
$\sqrt{10^{-3}}\approx 3\%$ change in the harmonic frequency along the
reaction coordinate (the vertical spacing between the lowest and the second
lowest symbols), generally not perceivable by the naked eye. By contrast, a
significant portion of the states with $l_{k}$ close to 0.5 shows a noticeable
horizontal shift, signifying the existence of level repulsion discussed above.
In addition, a change in the eigenstate structure induces a dramatic change in
the steady-state QY only if the shifted states are involved in the dominant
relaxation pathways. This is so in the case shown in Fig. 1(c) in contrast to
that in Fig. 1(b), in agreement with the respective behavior of the steady-
state QY shown in Fig. 2(b). These observations complement the parametric
sensitivity of long-time dynamics displayed in Figs. 2 and 3, supporting the
view that isolated system chaos propagates through the relaxation dynamics in
the open system.

Possible experimental verification— We note the recent work from Schnedermann
et al., where rhodopsins reconstituted with isotope-substituted retinal
chromophores show similar short-time dynamics, but very different steady-state
QY.29 This observation supports the view that a shorter rise time (the time it
takes for the first peak of photoproduct signal to appear) corresponds to
faster reaction, but not necessarily to higher reaction QY. That is, there is
a disconnect between transient and steady state behavior. However, changes
with isotope substitution are not directly comparable to the sensitivity of
steady-state QY seen in Figs. 2 and 3. First, the isotope labelings in these
experiments constitute much stronger perturbations to the system than those
studied here. Second, the experimental observations are typically, if not
exclusively, made on large ensembles of molecules. Within the ensemble, one
expects significant inhomogeneity among the parameters, which would wash out
the kind of parametric sensitivity observed in our computations.

For the real retinal system in an opsin protein pocket, the experimental
situation is expected to be much more complicated than we have modeled above.
Indeed, the established methods of measuring steady-state QY typically involve
macroscopic amounts of the sample being photoisomerized (for example loss of
$\approx 10\%$ of the optical density of the reactant over the course of a few
minutes30, 31). This being the case, and as noted above, it is likely that the
sensitivity observed here is washed out by the sheer number of retinal
molecules in a sample ($\approx 10^{20}$), each possessing a different
conformation. Hence, experimental confirmation of the sensitivity of QY will
require measurements of small ensembles or even of single molecules. Provided
that such sensitivity can be verified experimentally for retinal systems, it
may well be the case that one can find similar phenomena in other light-
induced condensed phase curve-crossing systems.

Open Theoretical Challenges— Additional investigation is needed to clarify the
relationship between quantum chaos and parametric sensitivity of meta-stable
state populations in an open system under environmental dissipation. Whereas
we continue to examine a number of possible measures of quantum chaos for open
systems proposed in the literature, such as the Loschmidt echos and the ratio
of eigenvalue differences of the Liouville operator,32, 33, 34, 35, 36 results
thus far are inconclusive. (These results are available upon reasonable
request.) For example, they suggest a need for a criteria that applies to open
systems whose isolated system component displays both integrable and chaotic
energy level regions. Further, it is worth noting that our treatment of
system-bath interactions is restricted to a Markovian and perturbative (fast
bath) limit, with other environmental effects included as static disorder
(slow bath). In doing so, we do not detect QY sensitivity to the parameters of
the system-bath coupling term or to the bath itself. In the future it would be
interesting to examine non-pertubative system-bath effects using more
sophisticated methods that also account for the intermediate regime in which
the system and the bath time scales overlap.37, 19

Furthermore, while admittedly premature, it is not difficult to suggest a
connection between this result and its possible significant biological
implications. That is, does a biological light-sensing apparatus benefit from
the high sensitivity of the underlying curve-crossing mechanism to parameter
variation? Do the system parameters have to be in a specific (“Goldilocks”)
regime of parameter space? Has natural selection tuned these parameters? These
and related questions warrant further study. What is the key here is the
highly significant possibility of quantum chaotic contributions and parametric
sensitivity in biological processes and its role in tuning biological
function.

Conclusion— The photoisomerization dynamics of retinal chromophore in
rhodopsin was studied using a two-state-two-mode model coupled to dissipative
phonon baths. Assuming the Franck-Condon state as the initial state, the
photoproduct (trans-conformer) populations first show a rapid rise with strong
oscillatory features (quantum beats) on the 100 fs time scale. The
oscillations then fade out and are replaced by a smooth time-dependent profile
leading to the stationary values on longer time scales. Time-dependent
behavior of this type agrees with available experimental observations of the
closely related system bacteriorhodopsin, suggesting that simple pictures that
correlate short time transient spectral features with steady-state quantum
yield are suspect.

Most significantly, we find that the steady-state quantum yield is highly
sensitive to slight changes in system parameters, while the transient dynamics
remains almost unaffected. Such sensitivity is also observed by changing
numerous other parameters of the model system, suggesting that it is an
inherent feature of the nonadiabatic vibronic system studied. Further,
possible experimental measurements of such sensitivity are discussed. The
implications of these results are profound in biological light-sensing systems
, and possibly in condensed phase photochemistry.

Acknowledgement— Support under AFOSR grant FA9550-19-1-0267 is gratefully
acknowledged, as are discussions with Professor Jennifer Ogilvie, University
of Michigan and Professor Arjendu Pattanayak, Carleton College.

## Supporting Information

Parametric sensitivity of steady-state photoisomerization observed with
incoherent initial conditions and incoherent excitations, parametric
sensitivity and nonsensitivity in single surface models.

## References

  * Ottolenghi 1982 Ottolenghi, M. Molecular Aspects of the Photocycles of Rhodopsin and Bacteriorhodopsin: A Comparative Overview. _Meth. Enzymol._ 1982, _88_ , 470–491
  * Wand et al. 2013 Wand, A.; Gdor, I.; Zhu, J.; Sheves, M.; Ruhman, S. Shedding New Light on Retinal Protein Photochemistry. _Ann. Rev. Phys. Chem._ 2013, _64_ , 437–458
  * Schulten 2014 Schulten, K. In _Quantum Effects in Biology_ ; Mohseni, M., Omar, Y., Engel, G. S., Plenio, M. B., Eds.; Cambridge University Press: Cambridge, UK, 2014; Chapter 11, pp 237–263
  * Schoenlein et al. 1991 Schoenlein, R.; Peteanu, L.; Mathies, R.; Shank, C. The First Step in Vision: Femtosecond Isomerization of Rhodopsin. _Science_ 1991, _254_ , 412–415
  * Wang et al. 1994 Wang, Q.; Schoenlein, R. W.; Peteanu, L. A.; Mathies, R. A.; Shank, C. V. Vibrationally Coherent Photochemistry in the Femtosecond Primary Event of Vision. _Science_ 1994, _266_ , 422–424
  * Polli et al. 2010 Polli, D.; Altoè, P.; Weingart, O.; Spillane, K. M.; Manzoni, C.; Brida, D.; Tomasello, G.; Orlandi, G.; Kukura, P.; Mathies, R. A. et al. Conical Intersection Dynamics of the Primary Photoisomerization Event in Vision. _Nature_ 2010, _467_ , 440
  * Johnson et al. 2017 Johnson, P. J.; Farag, M. H.; Halpin, A.; Morizumi, T.; Prokhorenko, V. I.; Knoester, J.; Jansen, T. L.; Ernst, O. P.; Miller, R. D. The Primary Photochemistry of Vision Occurs at the Molecular Speed Limit. _J. Phys. Chem. B_ 2017, _121_ , 4040–4047
  * Dartnall 1968 Dartnall, H. The Photosensitivities of Visual Pigments in the Presence of Hydroxylamine. _Vision Res._ 1968, _8_ , 339–358
  * Tscherbul and Brumer 2014 Tscherbul, T. V.; Brumer, P. Long-Lived Quasistationary Coherences in a V-type System Driven by Incoherent Light. _Phys. Rev. Lett._ 2014, _113_ , 113601
  * Brumer 2018 Brumer, P. Shedding (Incoherent) Light on Quantum Effects in Light-Induced Biological Processes. _J. Phys. Chem. Lett._ 2018, _9_ , 2946–2955
  * Hahn and Stock 2000 Hahn, S.; Stock, G. Quantum-Mechanical Modeling of the Femtosecond Isomerization in Rhodopsin. _J. Phys. Chem. B_ 2000, _104_ , 1146–1149
  * Hahn and Stock 2002 Hahn, S.; Stock, G. Ultrafast cis-trans Photoswitching: A Model Study. _J. Chem. Phys._ 2002, _116_ , 1085–1091
  * Marsili et al. 2019 Marsili, E.; Farag, M. H.; Yang, X.; De Vico, L.; Olivucci, M. Two-State, Three-Mode Parametrization of the Force Field of a Retinal Chromophore Model. _J. Phys. Chem. A_ 2019, _123_ , 1710–1719
  * Balzer and Stock 2005 Balzer, B.; Stock, G. Modeling of Decoherence and Dissipation in Nonadiabatic Photoreactions by an Effective-Scaling Nonsecular Redfield Algorithm. _Chem. Phys._ 2005, _310_ , 33–41
  * González-Luque et al. 2000 González-Luque, R.; Garavelli, M.; Bernardi, F.; Merchán, M.; Robb, M. A.; Olivucci, M. Computational Evidence in Favor of a Two-State, Two-Mode Model of the Retinal Chromophore Photoisomerization. _Proc. Nat. Acad. Sci._ 2000, _97_ , 9379–9384
  * Nakamichi and Okada 2006 Nakamichi, H.; Okada, T. Crystallographic Analysis of Primary Visual Photochemistry. _Angew. Chem. Int._ 2006, _45_ , 4270–4273
  * Tscherbul and Brumer 2014 Tscherbul, T. V.; Brumer, P. Excitation of Biomolecules with Incoherent Light: Quantum Yield for the Photoisomerization of Model Retinal. _J. Phys. Chem. A_ 2014, _118_ , 3100–3111
  * Dodin et al. 2018 Dodin, A.; Tscherbul, T.; Alicki, R.; Vutha, A.; Brumer, P. Secular versus Nonsecular Redfield Dynamics and Fano Coherences in Incoherent Excitation: An Experimental Proposal. _Phys. Rev. A_ 2018, _97_ , 013421
  * Sala and Egorova 2018 Sala, M.; Egorova, D. Quantum Dynamics of Multi-Dimensional Rhodopsin Photoisomerization Models: Approximate versus Accurate Treatment of the Secondary Modes. _Chem. Phys._ 2018, _515_ , 164–176
  * Hahn and Stock 2000 Hahn, S.; Stock, G. Femtosecond Secondary Emission Arising from the Nonadiabatic Photoisomerization in Rhodopsin. _Chem. Phys._ 2000, _259_ , 297–312
  * Prokhorenko et al. 2006 Prokhorenko, V. I.; Nagy, A. M.; Waschuk, S. A.; Brown, L. S.; Birge, R. R.; Miller, R. D. Coherent Control of Retinal Isomerization in Bacteriorhodopsin. _Science_ 2006, _313_ , 1257–1261
  * Fujisaki and Takatsuka 2001 Fujisaki, H.; Takatsuka, K. Chaos Induced by Quantum Effect due to Breakdown of the Born-Oppenheimer Adiabaticity. _Phys. Rev. E_ 2001, _63_ , 066221
  * Marsili et al. 2020 Marsili, E.; Olivucci, M.; Lauvergnat, D.; Agostini, F. Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model. _J. Chem. Theo. Comp._ 2020, _16_ , 6032–6048
  * Heller 1990 Heller, E. Mode Mixing and Chaos Induced by Potential Surface Crossings. _J. Chem. Phys._ 1990, _92_ , 1718
  * Leitner et al. 1996 Leitner, D. M.; Köppel, H.; Cederbaum, L. S. Statistical Properties of Molecular Spectra and Molecular Dynamics: Analysis of Their Correspondence in NO2 and C2H4. _J. Chem. Phys._ 1996, _104_ , 434
  * Haller et al. 1983 Haller, E.; Köppel, H.; Cederbaum, L. S. On the Statistical Behaviour of Molecular Vibronic Energy Levels. _Chem. Phys. Lett._ 1983, _101_ , 215
  * Santoro et al. 2007 Santoro, F.; Lami, A.; Olivucci, M. Complex Excited Dynamics around a Plateau on a Retinal-like Potential Surface: Chaos, Multi-Exponential Decays and Quantum/Classical Differences. _Theo. Chem. Acct._ 2007, _117_ , 1061–1072
  * Biggs 1993 Biggs, N. _Algebraic Graph Theory_ ; Cambridge University Press, 1993
  * Schnedermann et al. 2018 Schnedermann, C.; Yang, X.; Liebel, M.; Spillane, K.; Lugtenburg, J.; Fernández, I.; Valentini, A.; Schapiro, I.; Olivucci, M.; Kukura, P. et al. Evidence for a Vibrational Phase-Dependent Isotope Effect on the Photochemistry of Vision. _Nat. Chem._ 2018, _10_ , 449
  * Kim et al. 2001 Kim, J. E.; Tauber, M. J.; Mathies, R. A. Wavelength Dependent cis-trans Isomerization in Vision. _Biochem._ 2001, _40_ , 13774–13778
  * Sovdat et al. 2012 Sovdat, T.; Bassolino, G.; Liebel, M.; Schnedermann, C.; Fletcher, S. P.; Kukura, P. Backbone Modification of Retinal Induces Protein-like Excited State Dynamics in Solution. _J. Amer. Chem. Soc._ 2012, _134_ , 8318–8320
  * Grobe et al. 1988 Grobe, R.; Haake, F.; Sommers, H.-J. Quantum Distinction of Regular and Chaotic Dissipative Motion. _Phys. Rev. Lett._ 1988, _61_ , 1899–1902
  * Akemann et al. 2019 Akemann, G.; Kieburg, M.; Mielke, A.; Prosen, T. Universal Signature from Integrability to Chaos in Dissipative Open Quantum Systems. _Phys. Rev. Lett._ 2019, _123_ , 254101
  * Sá et al. 2020 Sá, L.; Ribeiro, P.; Prosen, T. Complex Spacing Ratios: A Signature of Dissipative Quantum Chaos. _Phys. Rev. X_ 2020, _10_ , 021019
  * Jalabert and Pastawski 2001 Jalabert, R. A.; Pastawski, H. M. Environment-Independent Decoherence Rate in Classically Chaotic Systems. _Phys. Rev. Lett._ 2001, _86_ , 2490–2493
  * Cucchietti et al. 2003 Cucchietti, F. M.; Dalvit, D. A. R.; Paz, J. P.; Zurek, W. H. Decoherence and the Loschmidt Echo. _Phys. Rev. Lett._ 2003, _91_ , 210403
  * Chen et al. 2016 Chen, L.; Gelin, M. F.; Chernyak, V. Y.; Domcke, W.; Zhao, Y. Dissipative Dynamics at Conical Intersections: Simulations with the Hierarchy Equations of Motion Method. _Faraday Discuss._ 2016, _194_ , 61–80


# Photonic Spin-Orbit Coupling Induced by Deep-Subwavelength Structured Light

Xin Zhang1,2, Guohua Liu1,2, Yanwen Hu1,2,3∗, Haolin Lin1,2, Zepei Zeng1,2,
Xiliang Zhang1,2, Zhen Li1,2,3, Zhenqiang Chen1,2,3, and Shenhe Fu1,2,3
huyanwen<EMAIL_ADDRESS>1Department of Optoelectronic
Engineering, Jinan University, Guangzhou 510632, China
2Guangdong Provincial Key Laboratory of Optical Fiber Sensing and
Communications, Guangzhou 510632, China
3Guangdong Provincial Engineering Research Center of Crystal and Laser
Technology, Guangzhou 510632, China

###### Abstract

We demonstrate both theoretically and experimentally beam-dependent photonic
spin-orbit coupling in a two-wave mixing process described by an equivalent of
the Pauli equation in quantum mechanics. The considered structured light in
the system is comprising a superposition of two orthogonal spin-orbit-coupled
states defined as spin up and spin down equivalents. The spin-orbit coupling
is manifested by prominent pseudo spin precession as well as spin-transport-
induced orbital angular momentum generation in a photonic crystal film of
wavelength thickness. The coupling effect is significantly enhanced by using a
deep-subwavelength carrier envelope, different from previous studies which
depend on materials. The beam-dependent coupling effect can find intriguing
applications; for instance, it is used in precisely measuring variation of
light with spatial resolution up to 15 nm.

##  I Introduction

Spin-orbit coupling (SOC), which refers to interaction of a quantum particle’s
spin with its momentum, is a fundamentally important concept. It has been
extensively investigated in condensed matter physics Winkler2003 ; Goldman2018
, atomic and molecular physics Spielman2011 ; Spielman2013  and contributes to
exciting phenomena such as the spin Hall effect Kato2004  and topological
insulators Hasan2010 ; Taylor2011 . Analogous photonic SOC is also
demonstrated in a variety of settings Zayats2015 . The photonic SOC refers to
an interaction between the momentum of light, which also includes spin angular
momentum and orbital angular momentum (SAM and OAM). Whereas the SAM is
associated with photon circular polarization Padgett2002 , the OAM is relevant
to a helical wavefront of light characterized by a topological number $\ell$
Allen1992 . The photonic SOC is crucial for the optical Hall effects
Nagaosa2004 ; Glazov2005 ; Kwiat2008 ; Fu2019 , spin-to-orbital angular
momentum conversions Marrucci2006 ; Chiu2007 , spin-orbit photonic devices
Gorodetski2008 ; Shitrit2013 ; Wang2021 , etc.
The SOC can be engineered in appropriately designed materials. For examples,
engineering a tensional strain in graphene shifts the electronic dispersions
and induces a controllable vector potential for the electronic SOC Crommie2010
; Low2010 ; Juan2011 ; Guinea2010 ; Bolotin2009 . Analogous strategy can be
applied to engineer the photonic SOC, by using strained evanescently coupled
waveguide arrays Rechtsman2013 ; Lumer2019 . Other approaches for manipulating
the photonic SOCs are demonstrated by appropriately designing microcavities
RO2011 ; Longhi2013 ; Fang2012 ; Yu2012 ; SF2013 ; Szczytko2021 ,
metamaterials Liu2015 ; Sala2015 ; Zhang2023 ; Lu2023 , photonic crystals
Yesharim2022 ; Karnieli2018 ; Liu2023 , twisted optical fibers Chiao1986 ;
Liberman1992 , dual-core waveguides Torner2015 ; Malomed2016 , etc. The
resultant SOCs are material-dependent, determined by geometric configurations
of the materials which are often difficult to be tuned once fixed by designs.
As a consequence, a tunable photonic SOC process remains elusive. Recently,
several engineered photonic SOC schemes have been reported, by either
embedding a strained honeycomb metasurface inside a cavity waveguide Mann2020
or using an optical cavity filled with controllable liquid crystals
Szczytko2019 . However, the resultant photonic SOCs remain material-dependent.
In this work, we report theoretically and experimentally a new mechanism for
engineering the photonic SOC. We demonstrate this by exploiting analogy
between quantum description of a spin-1/2 system and a spin-orbit Hamiltonian
derived for structured light in a photonic crystal. The obtained Hamiltonian
is closely relevant to structured light, which means that the SOC can be
engineered by controlling carrier envelope rather than the structures of
materials. Strong SOC is achieved by using deep-subwavelength structured
light, as manifested by clear pseudo spin precessions. Although the structured
light has been extensively investigated in recent year Forbes2021 ; Shen2019 ;
Fu2023 ; Femius2015 ; Hu2021 ; Fu2020 , the dependence of the photonic SOC on
its spatial structure remains unnoticed.

##  II Theoretical model

Figure 1: (a) Geometrical representation of spin precession in the presence of
SOC. $\hat{R}$ and $\hat{L}$ define the spin-up and spin-down equivalents in
the $B_{2}$ direction, respectively; whereas $\Phi_{+}$ and $\Phi_{-}$ denotes
two spin eigenstates in the direction of synthetic field B. Spin precession is
initiated by a mixing spin $\Phi=1/\sqrt{2}(\Phi_{+}+i\Phi_{-})$ located at
$z_{0}$. (b) Bloch-sphere representation of spin-1/2 system, in the presence
of external field B. $\Phi_{\uparrow}$ and $\Phi_{\downarrow}$ are spin up and
spin down in the $z$ direction, while $\Phi_{+}^{1/2}$ and $\Phi_{-}^{1/2}$
denotes eigenstates of the system, corresponding to direction of B. (c)
Polarization states mapped on a longitude line of the first-order ($\ell=1$)
sphere in (a). (d) Corresponding spin vectors to (c). (e) SOC strength as a
function of beam width $r_{0}$. (f)-(i) Theoretical results for the spin
vectors under actions of LG beam with different widths.

We consider a two-wave mixing process involving two interacting photonic
states. The SOC takes place in a crystal, represented by its principal
refractive index: $n_{x}$, $n_{y}$, and $n_{z}$. With an approximation of the
slowly varying envelope along optical axis $z$, a coupled-wave equation for
the process is given by Guo2018

$\displaystyle\begin{aligned} 2i\beta_{x}\frac{\partial E_{x}}{\partial
z}+\frac{n_{x}^{2}}{n_{z}^{2}}\frac{\partial^{2}E_{x}}{\partial
x^{2}}+\frac{\partial^{2}E_{x}}{\partial
y^{2}}=\gamma_{y}\frac{\partial^{2}E_{y}}{\partial y\partial
x}\exp\left(+i\Delta\beta\cdot z\right)\\\ 2i\beta_{y}\frac{\partial
E_{y}}{\partial
z}+\frac{n_{y}^{2}}{n_{z}^{2}}\frac{\partial^{2}E_{y}}{\partial
y^{2}}+\frac{\partial^{2}E_{y}}{\partial
x^{2}}=\gamma_{x}\frac{\partial^{2}E_{x}}{\partial x\partial
y}\exp\left(-i\Delta\beta\cdot z\right)\end{aligned}$ (1)

where $E_{x,y}$ are linearly polarized fields, and $\beta_{x,y}=k_{0}n_{x,y}$
denote their propagation constants. $k_{0}=2\pi/\lambda$ is free-space
wavenumber with $\lambda$ being the wavelength.
$\Delta\beta=\beta_{y}-\beta_{x}$ is a phase mismatch. We define
$\gamma_{x,y}=1-n_{x,y}^{2}/n_{z}^{2}$ as coupling parameters, related to
crystal’s polarity. The derivatives $\nabla^{2}_{xy}$ and $\nabla^{2}_{yx}$ in
Eq. (1) stem from the non-zero term $\nabla\cdot\textbf{E}\neq 0$, featuring
origin of the SOC Liberman1992 .
To address the rapid oscillation terms $\exp(\pm i\Delta\beta\cdot z)$, we
transform the wave equation to a rotating form, by defining

$\displaystyle\begin{aligned} E_{x}=\tilde{A}_{x}\exp(+i{\Delta}\beta\cdot
z/2)\\\ E_{y}=\tilde{A}_{y}\exp(-i{\Delta}\beta\cdot z/2)\end{aligned}$ (2)

respectively. Thus a Hamiltonian of the system is written as

$\textbf{H}=\frac{1}{2\bar{\beta}}\left[\begin{array}[]{cc}-\nabla_{\perp}^{2}+\bar{\gamma}\nabla^{2}_{xx},0\\\
0,-\nabla_{\perp}^{2}+\bar{\gamma}\nabla^{2}_{yy}\end{array}\right]+\left[\begin{array}[]{cc}\Delta\beta/2,\bar{\gamma}\nabla^{2}_{yx}/(2\bar{\beta})\\\
\bar{\gamma}\nabla^{2}_{xy}/(2\bar{\beta}),-\Delta\beta/2\end{array}\right]$
(3)

where $\nabla_{\perp}^{2}=\nabla^{2}_{xx}+\nabla^{2}_{yy}$ denotes the Laplace
operator. We have assumed shallow crystal birefringence, namely
$\bar{\beta}\approx(\beta_{x}+\beta_{y})/2$ and
$\bar{\gamma}\approx(\gamma_{x}+\gamma_{y})/2$. The second term in Eq. (3),
which includes the derivative operators, couples the two polarization
components. It means that the SOC is related to spatial structure of light.
We study the beam-dependent SOC in a synthetic two-level spin-orbit system. We
define right and left circularly polarized vortex states as spin up and spin
down equivalents in the $z$ direction. They are written as HuOL ; Naidoo2016 ;
Milione2011 :

$\displaystyle\begin{aligned}
\hat{R}=\exp(+i\ell\phi)(\hat{x}-i\hat{y})/\sqrt{2}\\\
\hat{L}=\exp(-i\ell\phi)(\hat{x}+i\hat{y})/\sqrt{2}\end{aligned}$ (4)

respectively, where $\hat{x}$ and $\hat{y}$ are unit vectors and
$\phi=\arctan(y/x)$. Since the pseudo spins are defined in the circular basis,
the Hamiltonian is modified by a transformation from the cartesian coordinate
to the circular basis, yielding

$\textbf{H}^{\prime}=\frac{\bar{\gamma}-2}{4\bar{\beta}}\left[\begin{array}[]{cc}\nabla_{\perp}^{2},0\\\
0,\nabla_{\perp}^{2}\end{array}\right]+\left[\begin{array}[]{cc}0,\Delta\beta/2-i\bar{\gamma}\nabla^{2}_{yx}/(2\bar{\beta})\\\
\Delta\beta/2+i\bar{\gamma}\nabla^{2}_{yx}/(2\bar{\beta}),0\end{array}\right]$
(5)

Given an overall field
$\tilde{\textbf{A}}=\tilde{A}(x,y,z)(\Phi_{R}\hat{R}+\Phi_{L}\hat{L})$, where
$\Phi_{R}$ and $\Phi_{L}$ are weights on $\hat{R}$ and $\hat{L}$,
respectively, we reduce Eq. (1) to the Schr$\ddot{\text{o}}$dinger-like
(Pauli) form

$i\frac{\partial\Phi(z)}{\partial
z}=\left(\frac{1}{2M}\textbf{P}_{\perp}^{2}{\tilde{A}}-\frac{1}{2}\sigma\cdot\textbf{B}\right)\Phi(z)$
(6)

where $\Phi=(\Phi_{R},\Phi_{L})^{\text{T}}$,
$\textbf{P}_{\perp}^{2}=[-\nabla_{\perp}^{2},0;0,-\nabla_{\perp}^{2}]$, and
$M=2\bar{\beta}\tilde{A}/(2-\bar{\gamma})$. Here $\sigma$ is the Pauli matrix
vector. The SOC is described by a term $-\sigma\cdot$B, where
$B_{1}=-\bar{\gamma}\nabla^{2}_{xy}\tilde{A}/(\bar{\beta}\tilde{A})$,
$B_{2}=0$, and $B_{3}=-\Delta\beta$. It is analogous to a coupling form which
describes interaction between a particle’s spin and its angular momentum in a
moving frame Winkler2003 ; Goldman2018 . More details refer to Appendix A.
Since $B_{2}$ is zero, the vector B lies on the purely transverse $B_{1}B_{3}$
plane. The SOC Hamiltonian admits eigenstates that point to the vector B and
comprise an equal superposition of $\hat{R}$ and $\hat{L}$, written as

$\displaystyle\begin{aligned}
\Phi_{+}=1/\sqrt{2}\left[\hat{R}+\exp(i\varphi)\hat{L}\right]\\\
\Phi_{-}=1/\sqrt{2}\left[\hat{R}-\exp(i\varphi)\hat{L}\right]\end{aligned}$
(7)

Figure 1(a) visualize the eigenstates and the pure states ($\hat{R}$ and
$\hat{L}$) in the Poincar$\acute{\text{e}}$ sphere, showing close analogies to
Bloch-sphere representation of the spin-1/2 system Scully1987 ; Berry1987
[Fig. 1(b)]. The spin states exhibit cylindrically symmetric polarization
distributions. As illustration, Fig. 1(c) displays typical polarizations of
states mapped onto a longitude line in the first-order ($\ell=1$) sphere;
while Fig. 1(d) depicts their corresponding spin vectors, represented by an
angle $\arccos(S_{2})$, where $S_{2}$ is value of polarization ellipticity.
Since the state exhibits identical polarization ellipticity in the transverse
plane, the resultant spin vectors are homogeneous.
The SOC term shows a dynamical effect, caused by the propagation-variant
envelope $\tilde{A}$. This shows sharp contrast to conventional ones which are
usually being independent terms. However, if optical diffraction is neglected,
the dynamical behavior disappears and the SOC strongly relies on the envelope.
In this scenario, a relevant beam parameter becomes an important degree of
freedom for engineering the SOC. This is demonstrated in Fig. 1(e), showing
close relationship between SOC strength and beam width in the phase-matching
condition ($\Delta\beta$=0). Here the Laguerre-Gaussian (LG) envelope is
considered as:

$\tilde{A}(r)=\frac{r}{r_{0}}\exp(-\frac{r^{2}}{r_{0}^{2}})$ (8)

where $r=(x^{2}+y^{2})^{1/2}$, and $r_{0}$ features the beam width. At the
deep-subwavelenth region ($r_{0}<\lambda/2$), the SOC strength is rapidly
increasing with a slight decrease of $r_{0}$. It becomes relatively negligible
when $r_{0}>\lambda$. This relation suggests that shrinking light to deep-
subwavelength scale significantly enhances the SOC. Although the derivations
are based on the slowly varying envelope approximation, the model can be
applied to deep-subwavelength regime at the early stage of spin evolution.
To demonstrate the deep-subwavelength-induced SOC, we set the coupling length
to be only one cycle ($z=\lambda$), such that a moderate SOC cannot cause
obvious spin transport phenomenon. On the other hand, the SOC strength can be
maintained during beam propagation, due to the short coupling length. This
results in an adiabatic spin evolution, represented as a spin precession
around B, i.e.,

$\frac{d\mathbf{S}}{dz}=\mathbf{B}\times\mathbf{S}$ (9)

where $\mathbf{S}=(S_{1},S_{2},S_{3})$ is the state vector defined as
$S_{h}=\Phi^{{\dagger}}\sigma_{h}\Phi$ ($h=1,2,3$). The spin vector is
therefore described by $S_{2}$. We initiate the spin precession from a mixing
state: $\Phi=1/\sqrt{2}[\Phi_{+}+i\Phi_{-}]$. Figure 1(f)-1(i) display
theoretical distributions of the spin vectors for different beam parameters.
Evidently, for $r_{0}$=0.05 $\mu$m, the spin rotates to an angle about
-78${}^{\text{o}}$; By comparison, increasing the parameter to $r_{0}=0.13$
$\mu$m causes less significant spin precession, manifested by a spin rotation
angle about -12${}^{\text{o}}$. This indicates that the SOC strongly depends
on the carrier envelope. Figure 1(h, i) show that a moderate SOC induced by
the relatively larger envelope cannot cause spin precession.

##  III Experimental results and discussion

Figure 2: (a) Experimental setup. BS: beam splitter; Q: q-plate; FL: flat
lens; M: mirror; OB: objective lens; TL: tube lens; QWP: quarter wave plate;
P: polarizer; CCD: charge coupled device. The laser is operating at wavelength
of $\lambda=632.8$ nm. The insert in (a) shows that an equatorial mixing spin
with equal weight on $\Phi_{R}$ and $\Phi_{L}$ is adiabatically converted to a
pure spin down in the presence of the SOC. (b) Layout of the 60-nm-thick flat
lens with NA=0.87. (c) Intensity distribution of the LG beam at the focal
plane ($z_{f}$) of the flat lens. (d) Plane-wave interference and (e)
$y$-polarization component of beam at $z_{f}$, indicating a generation of the
expected spin state $\Phi=1/\sqrt{2}(\Phi_{+}+i\Phi_{-})$. The scale bar in
(c-e) is 250 nm. In color bar, L: low; H: high.

Experiments are carried out to confirm the predictions. A crucial ingredient
is to generate the required spin-orbit state at the deep-subwavelength scale.
This is challenging since the incident state cannot maintain its property
after tightly focused by the high-numerical-aperture (NA) objective lens
Dorn2003 ; Wang2008 ; Xie2014 . To overcome this problem, we fabricate a
topology-preserving high-NA flat lens (the thickness is 60 nm) according to a
technique reported in Zhang2022 . The flat lens [the layout is shown in Fig.
2(b)] has a NA up to 0.87 and a focal length of $z_{f}=8$ $\mu$m. A system
comprising an objective lens (150$\times$, NA=0.9) and a tube lens is utilized
to characterize the flat lens, see Fig. 2(a). Figure 2(c) presents recorded
intensity distribution of light at the focal plane. The focused LG beam
exhibits a parameter of $r_{0}\simeq 0.32$ $\mu$m. The recorded regular
interference [(Fig. 2(d)] and $y$-polarization component [(Fig. 2(e)] suggest
that the expected initial spin is generated. Theoretical derivation about
topology-preserving property of the flat lens (Appendix B) further confirms
the generation.
An experimental setup is built for measuring the spin procession. A linearly
polarized He-Ne laser ($\lambda=632.8$ nm) is divided by a beam splitter. A
q-plate with a charge of $q=1/2$ is applied to transform the beam into
expected spin state carrier by the LG envelope. The purity of the spin state
from the q-plate is measured as 95.2% (Appendix C). The LG beam is focused
into deep-subwavelength region by the flat lens. A c-cut lithium niobate
crystal film ($\bar{\gamma}=-0.08$) with a thickness about one wavelength is
placed at the focal plane. The emerging beam, in the presence of the SOC, is
expected to accumulate a non-trivial spin phenomenon [see the insert in Fig.
2(a)].

Figure 3: Experimental observation of spin rotation induced by the deep-
subwavelength LG beam ($r_{0}=0.32$ $\mu$m), as manifested by spin angular
momentum conversion (flipping) from right-handed one (b) to the left-handed
one (a). In comparison, a larger LG beam parameter $r_{0}$=2.2 $\mu$m is
considered, resulting in balanced left-handed (c) and right-handed (d)
components. (e, f) The measured spin vectors before and after the crystal
film, for (e) $r_{0}=0.32$ $\mu$m, and (f) $r_{0}=2.2$ $\mu$m. In the color
bar (d), L: low; H: high.

Figure 3 presents measurements confirming the spin precession. Since the spin
is relevant to the circular polarization, we measure the right- (spin
$\uparrow$) and left-handed (spin $\downarrow$) circular polarization
components. These are achieved by rotating a quarter wave plate to an angle of
$-\pi/4$ and $+\pi/4$ with respect to $x$ axis, respectively, while inserting
a linear polarizer in front of the camera. Figure 3(a) and 3(b) depicts
intensity distributions of $\Phi_{L}$ and $\Phi_{R}$, respectively. The
measured $\Phi_{L}$ component is stronger than the $\Phi_{R}$ one, indicating
a spin precession toward south pole of the sphere. Figure 3(e) shows the
measured spin rotation by an angle of -5.2${}^{\text{o}}$, compared to the
initial one SR . This approximately matches to the simulated result. However,
for a larger parameter ($r_{0}=2.2$ $\mu$m), the $\Phi_{L}$ and $\Phi_{R}$
components are approximately identical [Fig. 3(c, d)], meaning that the
induced SOC is insufficient to flip the spin [Fig. 3(f)]. Slight difference
between the experiment and theory can be mainly attributed to the imperfect LG
envelope that is closely relevant to the derivative operator $\nabla_{xy}^{2}$
SR .
We observe non-trivial spin-precession phenomenon, manifested by a generation
of the photonic OAM. Initially, both the SAM and OAM of state at the equator
are zero. Under the action of the SOC, its intrinsic OAM and SAM are separated
simultaneously. This non-trivial phenomenon is observed in Fig. 4(a), showing
a clear dislocation in the plane-wave interference fringes for the deep-
subwavelength LG beam. This is a manifestation of wavefront helicity with a
topological charge being $\ell=1$. The spin precession accompanied by the OAM
generation confirms the phenomenon of spin-orbit separation. This effect
becomes negligible for larger envelope, since the spin remains at its original
position, as indicated by the regular interference fringes [Fig. 4(c)].
Theoretical results correspondingly shown in Fig. 4(b) and 4(d) are in
accordance with the measurements.

Figure 4: Observation of the orbital-angular-momentum state induced by spin
precession. (a, c) The experimentally measured plane-wave interference
patterns, for two different LG beam parameters: (a) $r_{0}=0.32$ $\mu$m, and
(c) $r_{0}=2.2$ $\mu$m. (b, d) The simulated [based on Eq. (1)] interference
patterns corresponding to the measurements in (a, c). Experimental conditions
are kept the same as those in Fig. 3.  Figure 5:  Observation of the spin
rotation by using the deep-subwavelength BG beam ($r_{0}=0.12$ $\mu$m). (a, b)
Experimentally measured intensity distributions of the left- and right-handed
circular polarizations. (c, d) Plane-wave interference patterns obtained both
in experiment (c) and in simulation (d). (e, f) The measured output spin
rotation in comparison with the initial one: (e) experiment; (f) simulation.

We observe more prominent spin precession by considering the Bessel structured
light with deeper subwavelength feature size. The carrier envelope is replaced
by

$\tilde{A}(r)=J_{\ell}(r/r_{0})$ (10)

where $J_{\ell}$ denotes the Bessel function of order $\ell$. In practice, we
should properly truncate the ideal Bessel beam by using a Gaussian factor. The
resultant Bessel-Gaussian (BG) profile exhibits nondiffracting property over a
certain distance. We generate this BG beam using a metasurface whose geometry
exhibits cylindrical symmetry. The highly localized BG beam is a result of in-
phase interference of many high-spatial-frequency waves Fu2020 . We
demonstrate result for a beam parameter of $r_{0}=0.12$ $\mu$m, while
maintaining other parameters unchanged. Similarly, an initial balance between
the left and right-handed components is broken by the SOC [Fig. 5(a, b)]. The
output spin rotates to a larger angle of -17.1${}^{\text{o}}$, nearly in
accordance with the theoretical calculation [Fig. 5(f)]. The measured and
simulated interference patterns verify the spin-precession-induced OAM
generation, see Fig. 5(c) and 5(d), respectively.

Figure 6: (a) The simulated [based on Eq. (6)] beam-dependent spin oscillatory
modes. (b)(c) The simulated [based on Eq. (1)] phase distributions of the
output light states from the barium metaborate crystal film
($\bar{\gamma}=-0.16$), for (b) $r_{0}=95$ nm, and (c) $r_{0}=110$ nm.

Finally, we propose to using the beam-dependent SOC in precision measurement
of slight variation of structured light, with measurement accuracy up to 15
nm. This nanometric resolution is usually impossible to be reached by current
optical detectors. It requires to realize rapid oscillation between the spin
up and spin down. Specifically, we exploit the deep-subwavelength BG beam as
carrier envelope of the spin. In this scenario, the Pauli equation [Eq. (6)]
emulates a SOC process for the spin oscillation. Figure 6(a) depicts the SOC-
supported spin harmonic oscillations along with the coupling distance, for
different cases of beam widths. Obviously, the spin oscillation is very
sensitive to the change of spatial structure of light, giving rising to
ultrasensitive beam-dependent oscillatory modes. As a result, a slight change
of the beam width leads to significant spin flipping. This allows to detect
the spatial variation of light as small as 15 nm. To verify the result, we
present simulated outcomes [see Fig. 6(b) and 6(c)], clearly showing opposite
helical wavefronts of the output states (corresponding to the spin down and
spin up), for $r_{0}$=95 nm and $r_{0}$=110 nm. Note that one can further
increase the measuring sensitivity by properly reducing the beam width.

##  IV conclusion

In summary, we have demonstrated both theoretically and experimentally novel
SOC phenomena, caused by the deep-subwavelength spin-orbit structured light.
This beam-dependent SOC contrasts to those being material-dependent Mann2020 ;
Szczytko2019 . The reported SOC is closely relevant to the spatial gradient of
light field, hence it can be significantly enhanced by using the deep-
subwavelength carrier envelopes. We have qualitatively characterized this
effect, by measuring the spin precessions under different beam parameters.
Particularly, based on the deep-subwavelength Bessel beam, a significant spin
rotation about -17.1${}^{\text{o}}$, accompanied by OAM generation, was
achieved within a coupling length of only one wavelength. The influence of the
phase mismatch on the beam-dependent SOC was also discussed, see Appendix D.
These fundamental SOC phenomena may find interesting applications in different
areas Mirhosseini2013 ; JChen2021 ; Vijay2012 ; Abele2018 . As an example, we
have proposed to use such a strong SOC effect in the precise measurement of
slight spatial change of light with nanometric resolution.

##  V ACKNOWLEDGMENTS

We thank Boris Malomed from Tel Aviv University for kind discussions about the
SOC. This work was supported by the National Natural Science Fundation of
China (62175091, 12304358), and the Guangzhou science and technology project
(202201020061).

##  VI APPENDIX A: Analogy of spin-orbit coupling in spin-1/2 system and
synthetic two-level system

Table 1: Analogies between the presented synthetic spin-1/2 system in the higher-order optical regime and the spin-1/2 system in the quantum mechanics. The direct analogies between these two different settings enable us to emulate intriguing spin transport phenomena in the presence of spin-orbit coupling. Physical parameters | Spin-1/2 system | Synthetic spin-1/2 system
---|---|---
Spins |  $\Phi_{\uparrow}$ and $\Phi_{\downarrow}$ |  $\hat{R}$ and $\hat{L}$
Eigenstates |  $\Phi_{+}^{1/2}$ and $\Phi_{-}^{1/2}$ |  $\Phi_{+}$ and $\Phi_{-}$
Field vector |  B (real) |  B=$(-\bar{\gamma}/(\bar{\beta}\tilde{A})\nabla_{xy}^{2}\tilde{A},0,-\Delta\beta)$
Spin-orbit coupling term | $\mathbf{H_{1/2}}=\sigma\cdot\mathbf{B}$ | $\mathbf{H_{SOC}}=\sigma\cdot\mathbf{B}$
Space/time coordinates | $(x,y,t)$ | $(x,y,z)$
Mass | $m$ | $M=2\bar{\beta}\tilde{A}/(2-\bar{\gamma})$
Momentum operator | P${}_{\perp}^{2}=[-\nabla_{\perp}^{2},0;0,-\nabla_{\perp}^{2}]$ | P${}_{\perp}^{2}=[-\nabla_{\perp}^{2},0;0,-\nabla_{\perp}^{2}]$

The spin-1/2 dynamics in the external vector field B can be described by a
Hamiltonian term $\mathbf{H_{1/2}}\mathrm{=}\sigma*\textbf{B}$, where $\sigma$
is the Pauli matrix vector. In a normalized form, it can be expressed as

$\mathbf{H_{1/2}}=\frac{1}{2}{\begin{bmatrix}\cos\theta&\mathrm{sin}\theta\cdot\mathrm{exp}\left(-i\varphi\right)\\\
\mathrm{sin}\theta\cdot\exp\left(i\varphi\right)&-\mathrm{cos}\theta\end{bmatrix}}$
(11)

where $\theta$ and $\varphi$ are two angles that define a normalized (unit)
sphere. The vector $\mathbf{B}$ then possesses around the sphere, with
direction determined by $\theta$ and $\varphi$. This Hamiltonian
$\mathbf{H_{1/2}}$ admits two spin eigenstates that point along to
$\mathbf{B}$, written as

$\begin{gathered}\Phi_{+}^{1/2}=\cos\left(\frac{\theta}{2}\right)\Phi_{\uparrow}+\exp(i\varphi)\sin\left(\frac{\theta}{2}\right)\Phi_{\downarrow}\\\
\Phi_{-}^{1/2}=\sin\left(\frac{\theta}{2}\right)\Phi_{\uparrow}-\exp(i\varphi)\cos\left(\frac{\theta}{2}\right)\Phi_{\downarrow}\end{gathered}$
(12)

where $\Phi_{\uparrow}=[1\quad 0]^{\mathrm{T}}$ and
$\Phi_{\downarrow}=[0\quad-1]^{\mathrm{T}}$ are spin-up and spin-down states
defined in the $z$ direction. Figure 1(b) geometrically depicts this picture
onto a Bloch sphere. All possible spins of the system can now be mapped onto
the sphere, with the spin up $\Phi_{\uparrow}$ and spin down
$\Phi_{\downarrow}$ located at the north and south poles of the sphere,
respectively. In the presence of the external field $\mathbf{B}$, the initial
spin precesses around the vector $\mathbf{B}$, giving rise to many intriguing
spin transport phenomena such as the geometric phase.
In our case, we study spin-orbit coupling of structured light in a photonic
crystal. The structured light in the system is comprising a superposition of
two orthogonal spin-orbit states with non-trivial topological structures. They
can be written as $\hat{R}=\exp{(il\phi)}(\hat{x}-i\hat{y})/\sqrt{2}$ and
$\hat{L}=\exp{(-il\phi)}(\hat{x}+i\hat{y})/\sqrt{2}$, respectively. These
topological states define the spin up and spin down equivalents along the $z$
axis, respectively, but they are not eigenstates of the analogous spin-orbit
Hamiltonian $\mathbf{H_{soc}}=-\sigma*\mathbf{B}$. In the circular basis, a
similar Hamiltonian matrix can be written as

$\mathbf{{H}_{soc}}=\begin{bmatrix}B_{2}&B_{3}-iB_{1}\\\
B_{3}+iB_{1}&-B_{2}\end{bmatrix}$ (13)

In our case, since $B_{2}$ is zero (see the main text), the effective vector B
obtained here lies on the purely transverse $B_{1}$$B_{3}$ plane, as shown in
Fig. 1(a). As a result, the pseudospin eigenstates of $\mathbf{{H}_{soc}}$
that point along this transverse vector B comprise an equal superposition of
$\hat{R}$ and $\hat{L}$, express as

$\displaystyle\Phi_{+}$
$\displaystyle=\cos\left(\frac{\pi}{4}\right)\hat{R}+\exp(i\varphi)\sin\left(\frac{\pi}{4}\right)\hat{L}$
(14) $\displaystyle\Phi_{-}$
$\displaystyle=\sin\left(\frac{\pi}{4}\right)\hat{R}-\exp(i\varphi)\cos\left(\frac{\pi}{4}\right)\hat{L}$

We can now interpret these eigenstates as a mixing of $\hat{R}$ and $\hat{L}$.
Poincaré-sphere representation allows us to visualize these spin eigenstates
as well as the pure states $\hat{R}$ and $\hat{L}$. Clearly, this is analogous
to the Bloch-sphere representation for the spin-1/2 system. The spin-orbit
coupling makes this state evolves along the Poincaré sphere, which can be
described by the synthetic Pauli equation,

$i\frac{\partial\Phi}{\partial
z}=\left(\frac{1}{2M}\mathbf{P}_{\perp}^{2}\tilde{A}-\frac{1}{2}\sigma\cdot\mathbf{B}\right)\Phi$
(15)

Table I summaries analogous formulas between these two systems.

##  VII APPENDIX B: Theoretical derivation for the topology-preserving flat
lens

In this section, we theoretically prove that the flat lens used in the
experiment does not change the spin-orbit property of the LG beam after
tightly focusing. The flat lens is designed by an amplitude-only hologram
generated from an interference between an angular cosine wave and a spherical
wave (see ref. [61] in the text). When the LG beam $\tilde{A}(x,y)$ carrying a
general spin state $\Phi$ passes through the flat lens, it is modulated in
binary. As a result, the light field behind the flat lens can be expressed as

$\mathbf{E}(x,y,z=0)=\tilde{A}(x,y)*t(x,y)\big{[}\Phi_{x}(\phi)\hat{x}+\Phi_{y}(\phi)\hat{y}\big{]}$
(16)

where $t(x,y)$ denotes transmission function of the flat lens and
$\phi=\arctan\left(y/x\right)$. Within this initial condition, we solve the
diffractive problem according to the vectorial Helmholtz wave equation. The
diffractive field at the focal plane of the flat lens can be written as

$\begin{split}\mathbf{E}(x,y,z_{f})=&\frac{k}{i2\pi
z_{f}}\iint\mathbf{E}(x^{\prime},y^{\prime},z=0)\\\
&\exp\left\\{\frac{ik}{2z_{f}}[(x-x^{\prime})^{2}+(y-y^{\prime})^{2}]\right\\}dx^{\prime}dy^{\prime}\end{split}$
(17)

Note that owing to the cylindrical symmetry of the flat lens (see the layout
in the text, Fig. 2(b)), the transmission function can be also given in a
cylindrical form of $t(r)$, where $r=(x^{2}+y^{2})^{1/2}$. In this case, the
complex amplitude of the initial field is separable in the polar coordinates
$(r,\phi)$. We therefore rewrite the solution in the cylindrical coordinate
system and deal with the integrals. We finally obtain the analytical solution
for the vectorial light field at the focal plane, given by

$\mathbf{E}(x,y,z_{f})=f(r)\big{[}\Phi_{x}(\phi)\hat{x}+\Phi_{y}(\phi)\hat{y}\big{]}$
(18)

where

$\begin{split}f(r)=&-\frac{k}{z_{f}}\int_{0}^{\infty}\tilde{A}(r^{\prime})t(r^{\prime})r^{\prime}J_{1}\left(\frac{krr^{\prime}}{z_{f}}\right)\\\
&\exp\left[\frac{ik}{2z_{f}}(r^{2}+r^{\prime
2})\right]\mathrm{d}r^{\prime}\end{split}$ (19)

and $J_{1}$ indicates the first-order Bessel function. It is evident that the
diffractive field at the focal plane shares a similar analytic form to the
initial one, except for that the envelope becomes a $z$-dependent function. It
indicates that the flat lens can completely retain the initial spin state when
it is focused into the input end of the crystal. The topology-preserving flat
lens enables us to detect the pseudo spin precession caused by the deep-
subwavelength structured light, which cannot be achieved by using the
conventional high NA objective lens.

##  VIII APPENDIX C: Purity measurement of the first-order LG beam from the
q-plate

Figure 7: Modal decomposition results. (a) An experimental setup used to
measure the purity of the first-order LG beam emerging from the q-plate. The
LG beam is decomposed into LG basis modes. The linearly polarized He-Ne laser
operating at the wavelength of 632.8 nm is considered. QWP, quarter wave
plate; QP, q-plate with a topological number of $q=1/2$; SLM, spatial light
modulator; BS, beam splitter; CCD, charge-coupled device. (b) The modal
decomposition results at the basis of LG modes with topological charge ranging
from $l=$ -5 to 5.

we perform additional experiment to show that the generated first-order LG
beam from the q-plate is of high purity, which is sufficiently enough to
detect the photonic spin-orbit coupling effect. We utilize a modal
decomposition method schulze2013 ; wei2019  to measure the purities of the
output LG mode from the q-plate with a topological charge of $q=1/2$, see an
experimental setup in Fig. 7(a). Two quarter wave plates (QWPs) are used to
select a proper polarization of the generated first-order ($l$=1) LG beam that
matches to the spatial light modulator (SLM). A group of pure LG modes
generated from digital holograms by using the SLM are considered to decomposed
the LG beam. Fig. 7(b) shows the decomposing result depicted in a histogram.
It is seen that the measured purity of the first-order LG beam from the
q-plate is 95.2$\%$.

Figure 8: Controllable spin-orbit coupling by engineering the phase mismatch
in a c-cut electro-optic lithium niobate crystal. (a) Experimental scheme for
observing the electrically engineered spin-orbit coupling. (b-f)
Experimentally measured photonic spin states at different applied voltages:
(b) $U$=60 V, (c) $U$=80 V, (d) $U$=100 V, (e) $U$=120 V, and (f) $U$=140 V.
(g) The measured topological charge as a function of applied voltage. Sim:
simulation; Expt: experiment. In this experiment, the coupling length of the
crystal is set to $z$=30 $mm$.

##  IX APPENDIX D: Engineering photonic spin-orbit coupling by tuning the
phase mismatch

In addition to the beam-dependent photonic spin-orbit coupling which we have
shown in the main text, we perform additional experiments confirming that the
spin-orbit coupling can be also controlled by engineering the phase mismatch.
To this end, we consider electrically tuning the phase mismatch in a c-cut
electro-optic lithium niobate (LN) crystal, whose optical axis is in
accordance with propagation direction of the beam [see Fig. 8(a)]. In the
presence of transverse modulation, the phase mismatch can be written as

$\Delta\beta=-k_{0}n_{o}^{3}\gamma_{22}U/d$ (20)

where $k_{0}=2\pi/\lambda$ denotes wavenumber in vacuum with $\lambda$ being
wavelength, $n_{o}$ is the ambient refractive index of the crystal, $U$ is the
applied voltage, $d$ is the thickness, and $\gamma_{22}=6.8$ $\mathrm{pm/V}$
is an electro-optic coefficient of the crystal. In this case, the external
knob $U$ is utilized to finely tune the phase mismatch and the resulting spin-
orbit coupling. We use the same experimental setup and obtain a voltage-
dependent transition between different spin states in the phase mismatching
regime. Panels (b-f) in Fig. 8 show controllable spin states of light by
varying the applied voltage. Moreover, we perform detailed experiments to
measure the topological charge of the output light state as a function of
voltage, see Fig. 8(g). These results suggest another important degree of
freedom for engineering the spin-orbit coupling.

## References

  * (1) R. Winkler, “Spin-Orbit Coupling Effects in Two-Dimensional Electron and Hole Systems,” Springer Tracts Mod. Phys. 191, Springer, New York (2003).
  * (2) M. Aidelsburger, S. Nascimbene, and N. Goldman, “Artifical gauge fields in materials and engineered systems,” C. R. Phys. 19, 394 (2018).
  * (3) Y. -J. Lin, K. Jim$\acute{\text{e}}$nez-Garc$\acute{\text{i}}$a, and I. B. Spielman, “Spin-orbit-coupled Bose-Einstein condensates,” Nature 471, 83 (2011).
  * (4) V. Galitski, and I. B. Spielman, “Spin-orbit coupling in quantum gases,” Nature 494, 49 (2013).
  * (5) Y. K. Kato, R. C. Myers, A. C. Gossard, and D. D. Awschalom, “Observation of the spin Hall effect in semiconductors,” Science 306, 1910 (2004).
  * (6) M. Z. Hasan, and C. L. Kane, “Colloquium: Topological insulators,” Rev. Mod. Phys. 82, 3045 (2010).
  * (7) M. Hafezi, E. A. Demler, M. D. Lukin, and J. M. Taylor, “Robust optical delay lines with topologial protection,” Nat. Phys. 7, 907 (2011).
  * (8) K. Y. Bliokh, F. J. Rodr$\acute{\text{i}}$guez-Fortu$\tilde{\text{n}}$o, F. Nori, and A. V. Zayats, “Spin-orbit interactions of light,” Nat. Photonics 9, 796 (2015).
  * (9) A. T. O’Neil, I. MacVicar, L. Allen, and M. J. Padgett, “Intrinsic and Extrinsic Nature of the Orbital Angular Momentum of a Light Beam,” Phys. Rev. Lett. 88, 053601 (2002).
  * (10) L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, “Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes,” Phys. Rev. A 45, 8185 (1992).
  * (11) M. Onoda, S. Murakami, and N. Nagaosa, “Hall Effect of Light,” Phys. Rev. Lett. 93, 083901 (2004).
  * (12) A. Kavokin, G. Malpuech, and M. Glazov, “Optical Spin Hall Effect,” Phys. Rev. Lett. 95, 136601 (2005).
  * (13) O. Hosten, and P. Kwiat, “Observation of the Spin Hall Effect of Light via Weak Measurements,” Science 319, 787 (2008).
  * (14) S. Fu, C. Guo, G. Liu, Y. Li, H. Yin, Z. Li, and Z. Chen, “Spin-Orbit Optical Hall Effect,” Phys. Rev. Lett. 123, 243904 (2019).
  * (15) L. Marrucci, C. Manzo, and D. Paparo, “Optical Spin-to-Orbital Angular Momentum Conversion in Inhomogeneous Anisotropic Media,” Phys. Rev. Lett. 96, 163905 (2006).
  * (16) Y. Zhao, J. S. Edgar, G. D. M. Jeffries, D. McGloin, and D. T. Chiu, “Spin-to-Orbital Angular Momentum Conversion in a Strongly Focused Optical Beam,” Phys. Rev. Lett. 99, 073901 (2007).
  * (17) Y. Gorodetski, A. Niv, V. Kleiner, and E. Hasman, “Observation of the Spin-Based Plasmonic Effect in Nanoscale Structures,” Phys. Rev. Lett. 101, 043903 (2008).
  * (18) N. Shitrit, I. Yulevich, E. Maguid, D. Ozeri, D. Veksler, V. Kleiner, and E. Hasman, “Spin-Optical Metamaterial Route to Spin-Controlled Photonics,” Science 340, 724 (2013).
  * (19) L. Fang, H. Wang, Y. Liang, H. Cao, and J. Wang, “Spin-Orbit Mapping of Light,” Phys. Rev. Lett. 127, 233901 (2021).
  * (20) N. Levy, S. A. Burke, K. L. Meaker, M. Panlasigul, A. Zettl, F. Guinea, A. H. C. Neto, and M. F. Crommie, “Strain-Induced Pseudo-Magnetic Feidls Greater Than 300 Tesla in Graphene Nanobubbles,” Science 329, 544 (2010).
  * (21) T. Low, and F. Guinea, “Strain-Induced Pseudomagnetic Field for Novel Graphene Electronics,” Nano Lett. 10, 3551 (2010).
  * (22) F. de Juan, A. Cortijo, M. A. H. Vozmediano, and A. Cano, “Aharonov-Bohm interferences from local deformations in graphene,” Nat. Phys. 7, 810 (2011).
  * (23) F. Guinea, M. I. Katsnelson, and A. K. Geim, “Energy gaps and a zero-field quantum Hall effect in graphene by strain engineering,” Nat. Phys. 6, 30 (2010).
  * (24) K. I. Bolotin, F. Ghahari, M. D. Shulman, H. L. Stormer, and P. Kim, “Observation of the fractional quantum Hall effect in graphene,” Nature 462, 196 (2009).
  * (25) M. C. Rechtsman,J. M.Zeuner, A. T$\ddot{\text{u}}$nnermann ,S. Nolte, M. Segev, and A. Szameit, “Strain-induced pseudomagnetic field and photonic Landau levels in dielectric structures,” Nat. Photonics 7, 153 (2013).
  * (26) Y. Lumer, M. A. Bandres, M. Heinrich, L. J. Maczewsky, H. Herzig-Sheinfux, A. Szameit, and M. Segev, “Light guiding by artificial gauge fields,” Nat. Photonics 13, 339 (2019).
  * (27) R. O. Umucal1lar, and I. Carusotto, “Artificial gauge field for photons in coupled cavity arrays,” Phys. Rev. A 84, 043804 (2011).
  * (28) S. Longhi, “Effective magnetic field for photons in waveguide and coupled resonator lattices,” Opt. Lett. 38, 3570 (2013).
  * (29) K. Fang, Z. Yu, and S. Fan, “Realizing effective magnetic field for photons by controlling the phase of dynamic modulation,” Nat. Photonics 6, 782 (2012).
  * (30) K. Fang, Z. Yu, and S. Fan, “Photonic Aharonov-Bohm Effect Based on Dynamic Modulation,” Phys. Rev. Lett. 108, 153901 (2012).
  * (31) K. Fang, and S. Fan, “Controlling the Flow of Light Using the Inhomogeneous Effective Gauge Field that Emerges from Dynamic Modulation,” Phys. Rev. Lett. 111, 203901 (2013).
  * (32) M. Kr$\acute{\text{o}}$l, K. Rechci$\acute{\text{n}}$ska, H. Sigurdsson, P. Oliwa, R. Mazur, P. Morawiak, W. Piecek, P. Kula, P. G. Lagoudakis, M. Matuszewski, W. Bardyszewski, B. Pietka, and J. Szczytko, “Realizing Optical Persistent Spin Helix and Stern-Gerlach Defection in an Anisotropic Liquid Crystal Microcavity,” Phys. Rev. Lett. 127, 190401 (2021).
  * (33) F. Liu, and J. Li, “Gauge Field Optics with Anisotropic Media,” Phys. Rev. Lett. 114, 103902 (2015).
  * (34) V. G. Sala, D. D. Solnyshkov, I. Carusotto, T. Jacqmin, A. Lemaitre, H. Tercas, A. Nalitov, M. Abbarchi, E. Galopin, I. Sagnes, J. Bloch, G. Malpuech, and A. Amo, “Spin-Orbit Coupling for Photons and Polaritons in Microstructures,” Phys. Rev. X 5, 011034 (2015).
  * (35) S. Ma, H. Jia, Y. Bi, S. Ning, F. Guan, H. Liu, C. Wang, and S. Zhang, “Gauge Field Induced Chiral Zero Mode in Five-Dimensional Yang Monopole Metamaterials,” Phys. Rev. Lett. 130, 243801 (2023).
  * (36) Q. Yan, Z. Wang, D. Wang, R. Ma, C. Lu, G. Ma, X. Hu, and Q. Gong, “Non-Abelian gauge field in optics,” Adv. Opt. Photonics 15, 907 (2023).
  * (37) O. Yesharim, A. Karnieli, S. Jackel, G. Di Domenico, S. Trajtenberg-Mills, and A. Arie, “Observation of the all-optical Stern–Gerlach effect in nonlinear optics,” Nat. Photonics 16, 582 (2022).
  * (38) A. Karnieli, and A. Arie, “All-optical Stern-Gerlach effect,” Phys. Rev. Lett. 120, 053901 (2018).
  * (39) G. Liu, S. Fu, S. Zhu, H. Yin, Z. Li, and Z. Chen, “Higher-order optical rabi oscillations,” Fundamental Research 3, 898 (2023).
  * (40) A. Tomita, and R. Y. Chiao, “Observation of Berry’s topological phase by use of an optical fiber,” Phys. Rev. Lett. 57, 937 (1986).
  * (41) V. S. Liberman, and B. Y. Zel’dovich, “Spin-orbit interaction of a photon in an inhomogeneous medium,” Phys. Rev. A 46, 5199 (1992).
  * (42) Y. V. Kartashov, B. A. Malomed, V. V. Konotop, V. E. Lobanov, and L. Torner, “Stabilization of spatiotemporal solitons in Kerr media by dispersive coupling,” Opt. Lett. 40, 1045 (2015).
  * (43) H. Sakaguchi, and B. A. Malomed, “One- and two-dimensional solitons in PT-symmetric systems emulating spin-orbit coupling,” New. J. Phys. 18, 105005 (2016).
  * (44) C. Mann, S. A. R. Horsley, and E. Mariani, “Tunable pseudo-magnetic fields for polaritons in strained metasurfaces,” Nat. Photonics 14, 669 (2020).
  * (45) K. Rechci$\acute{\text{n}}$ska, M. Kr$\acute{\text{o}}$l, R. Mazur, P. Morawiak, R. Mirek, K. Lempicka, W. Bardyszewski, M. Matuszewski, P. Kula, W. Piecek, P. G. Lagoudakis, B. Pietka, and J. Szczytko, “Engineering spin-orbit synthetic Hamiltonians in liquid-crystal optical cavities,” Science 366, 727 (2019).
  * (46) A. Forbes, M. de Oliveira, and M. R. Dennis, “Structured light,” Nat. Photonics 15, 253 (2021).
  * (47) Y. Shen, X. Wang, Z. Xie, C. Min, X. Fu, Q. Liu, M. Gong, and X. Yuan, “Optical vortices 30 years on: OAM manipulation from topological charge to multiple singularities,” Light: Sci. Appl. 8, 90 (2019).
  * (48) J. Jia, H. Lin, Y. Liao, Z. Li, Z. Chen, and S. Fu, “Pendulum-type light beams,” Optica 10, 90 (2023).
  * (49) A. F. Koenderink, Andrea Al$\grave{\text{u}}$, and A. Polman, “Nanophotonics: Shrinking light-based technology,” Science 348, 516 (2015).
  * (50) Y. Hu, S. Wang, J. Jia, S. Fu, H. Yin, Z. Li, and Z. Chen, “Optical superoscillatory waves without side lobes along a symmetric cut,” Adv. Photonics 3, 045002 (2021).
  * (51) Y. Hu, S. Fu, H. Yin, Z. Li, Z. Li, and Z. Chen, “Subwavelength generation of nondiffracting structured light beams,” Optica 7, 1261 (2020).
  * (52) C. Guo, S. Fu, H. Lin, Z. Li, H. Yin, and Z. Chen, “Dynamic control of cylindrical vector beams via anisotropy,” Opt. Express 26, 18721 (2018).
  * (53) Y. Hu, G. Mo, Z. Ma, S. Fu, S. Zhu, H. Yin, Z. Li, and Z. Chen, “Vector vortex state preservation in Fresnel cylindrical diffraction,” Opt. Lett. 46, 1313 (2021).
  * (54) D. Naidoo, F. S. Roux, A. Dudley, I. Litvin, B. Piccirillo, L. Marrucci, and A. Forbes, “Controlled generation of higher-order Poincaré sphere beams from a laser,” Nat. Photonics 10, 327 (2016).
  * (55) G. Milione, H. I. Sztul, D. A. Nolan, and R. R. Alfano, “Higher-Order Poincar$\acute{\text{e}}$ Sphere, Stokes Parameters, and the Angular Momentum of Light,” Phys. Rev. Lett. 107, 053601 (2011).
  * (56) M. O. Scully, W. E. Lamb, and A. Barut, “Theory of the Stern-Gerlach apparatus,” Found. Phys. 17, 575 (1987).
  * (57) M. V. Berry, “The adiabatic phase and Pancharatnam’s phase for polarized light,” J. Mod. Opt. 34, 1401 (1987).
  * (58) R. Dorn, S. Quabis, and G. Leuchs, “Sharper Focus for a Radially Polarized Light Beam,” Phys. Rev. Lett. 91, 233901 (2003).
  * (59) X. Xie, Y. Chen, K. Yang, and J. Zhou, “Harnessing the Point-Spread Function for High-Resolution Far-Field Optical Microscopy,” Phys. Rev. Lett. 113, 263901 (2014).
  * (60) H. Wang, L. Shi, B. Lukyanchuk, C. Sheppard, and C. T. Chong, “Creation of a needle of longitudinally polarized light in vacuum using binary optics,” Nat. Photonics 2, 501 (2008).
  * (61) X. Zhang, Y. Hu, X. Zhang, S. Zhu, H. Yin, Z. Li, Z. Chen, and S. Fu, “A Global Phase-Modulation Mechanism for Flat-Lens Design,” Adv. Opt. Mater. 10, 2202270 (2022).
  * (62) In experiment, the spin rotation was obtained by measuring $S_{2}=|\Phi_{R}|^{2}-|\Phi_{L}|^{2}$. Ideally, the spatial distribution of $S_{2}$ is homogeneous; whereas in experiment it exhibits fluctuations due to the imperfect beam envelope. We therefore average it by using an expression: $\bar{S}_{2}=\sum_{n}^{N_{x}}\sum_{m}^{N_{y}}S_{2,n,m}/(N_{x}N_{y})$.
  * (63) M. Mirhosseini, M. Malik, Z. Shi, and R. W. Boyd, “Efficient separation of the orbital angular momentum eigenstates of light,” Nat. Commun. 4, 2781 (2013).
  * (64) J. Chen, C. Wan, and Q. Zhan, “Engineering photonic angular momentum with structured light: a reivew,” Adv. Photonics 3, 064001 (2021).
  * (65) R. Vijay, C. Macklin, D. H. Slichter, S. J. Weber, K. W. Murch, R. Naik, A. N. Korotkov, and I. Siddiqi, “Stabilizing Rabi oscillations in a superconducting qubit using quantum feedback,” Nature 490, 77 (2012).
  * (66) G. Cronenberg, P. Brax, H. Filter, P. Geltenbort, T. Jenke, G. Pignol, M. Pitschmann, M. Thalhammer, and H. Abele, “Acoustic Rabi oscillations between gravitational quantum states and impact on symmetron dark energy,” Nat. Phys. 14, 1022 (2018).
  * (67) C. Schulze, A. Dudley, D. Flamm, M. Duparr$\acute{\text{e}}$, and A. Forbes, Measurement of the orbital angular momentum density of light by modal decomposition, New J. Phys. 15, 073025 (2013).
  * (68) D. Wei, Y. Cheng, R. Ni, Y. Zhang, X. Hu, S. Zhu, and M. Xiao, Generating controllable Laguerre-Gaussian laser modes through intracavity spin-orbital angular momentum conversion of light, Phys. Rev. Appl. 11, 014038 (2019).