Particle-like topologies in light

Three-dimensional (3D) topological states resemble truly localised, particle-like objects in physical space. Among the richest such structures are 3D skyrmions and hopfions, that realise integer topological numbers in their configuration via homotopic mappings from real space to the hypersphere (sphere in 4D space) or the 2D sphere. They have received tremendous attention as exotic textures in particle physics, cosmology, superfluids, and many other systems. Here we experimentally create and measure a topological 3D skyrmionic hopfion in fully structured light. By simultaneously tailoring the polarisation and phase profile, our beam establishes the skyrmionic mapping by realising every possible optical state in the propagation volume. The resulting light field’s Stokes parameters and phase are synthesised into a Hopf fibration texture. We perform volumetric full-field reconstruction of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\Pi }}}_{{{3}}}$$\end{document}Π3 mapping, measuring a quantised topological charge, or Skyrme number, of 0.945. Such topological state control opens avenues for 3D optical data encoding and metrology. The Hopf characterisation of the optical hypersphere endows a fresh perspective to topological optics, offering experimentally-accessible photonic analogues to the gamut of particle-like 3D topological textures, from condensed matter to high-energy physics.

So-called "baby skyrmions" are the two-dimensional (2D) counterpart of 3D skyrmions: fields in 2D physical space which map to, and wrap around, a 2-sphere parameter space. Their study is much more developed in theory and experiments, notably in non-singular superfluid vortices 13 including those imprinted by structured light 14 , and especially magnetic systems 15 . Here the direction of spin at each point provides the 2-sphere parameter space, and magnetic skyrmion excitations have the potential to represent topological bits for low-power computer memory and processing 15 . Recently, 2D baby skyrmion configurations were created in optical systems, as the direction of electric field vectors, or photon spin, near a material interface 16,17 , displaying dynamics similar to magnetic skyrmions 18 . In propagating laser light, optical polarization can be structured into full Poincaré beams 21 , which realise every state of elliptic polarization in the transverse plane. These beams can also be interpreted as 2D baby skyrmions 22 , since the Poincaré sphere, as the 2sphere parameter space, parametrises transverse, elliptic polarization states. However, 3D particle-like topological objects have not been considered either theoretically or experimentally in optical fields.
Optical realisations of 3D topological states can take various forms. Much interest has focused on singularity lines, such as optical vortices or polarization singularities (e.g. C lines) 22 . In structured light, with amplitude, phase, and polarization spatially varying, these can be woven into loops, links and knots 23,24 and organise Möbius strips 25 . The state of elliptic polarization is right-or left-handed circular (RH, LH) on C lines, often described as a skeleton of the complex optical polarization field 26 . Topologically structured light has a wide range of applications including enhanced free-space optical communications 27 and advanced trapping 28 , and is related to optical currents 29 and orbital angular momentum 30 . Singular lines are topologically characterised by the fundamental homotopy group Π ! . The homotopy group Π " , on the other hand, defines topological particles such as 3D hopfions and skyrmions 2 . It is natural to ask whether these 3D excitations can be created in structured light. We describe here the design, generation and measurement of a structured, propagating beam of laser light realising such a mapping, unifying particle-like 3D topologies in free space optics with those studied in high-energy physics, cosmology and various kinds of condensed matter.

The optical hypersphere of polarization and phase
Spatially extended polarized light is represented by a complex transverse electric field vector at each point " in the propagating beam. Its RH and LH components are represented by the complex-valued scalar functions # # (") and # $ ("), and the pair (# # , # $ ) which characterises the optical state at each point is assumed normalized, i.e.
Therefore, this normalized optical field defines a mapping from each point in 3D real space to a point on the 3-sphere, which we call the optical hypersphere. The optical hypersphere is conveniently parametrised using spinorial angles ., /, 0: for 0 ≤ / ≤ 8, −8 < . ≤ 8 and −28 < 0 ≤ 28. The angles ., /, 0 have a direct interpretation in terms of the polarization and phase of the electric field state: with < ! , < % , < " the normalized Stokes parameters, . = arctan(< ! , < % ) is the polarization azimuth, and cos / = < " is the polarization ellipticity; 0 = arg # # + arg # $ is the sum of the two electric field components' phases 26,31 . Further details of these parameters and their relationship with the hypersphere and the Poincaré sphere (2-sphere) parametrising polarization may be found in the Supplementary Information.
The full Poincaré sphere of polarization states can be realised in a transverse plane of a structured light field, created from the superposition of two, differently structured, LH and RH beam components, similar to a full Poincaré beam 21 . At each spatial point, the optical field has some elliptical polarization state characterised by ., /. In 3D, points of constant elliptical polarization lie on filaments, generalizing RH and LH circular polarized C lines. 3D real space is filled by the set of polarization filaments, constituting a polarization texture (Fig. 1a). Each filament corresponds to a point on the Poincaré sphere (Fig. 1b), and many filaments cross each plane (Fig. 1c). Although the polarization is fixed on the filaments, the optical phase smoothly varies along them (Fig. 1c insets). Any 3D structured light field with varying transverse polarization can be represented by such a texture.
The 3-sphere supports the Hopf fibration 19 , a fibre bundle which divides it into linked circles. In the optical hypersphere, each fixed polarization state (with ., / constant) traces out a circle as the phase 0 goes through a 48 cycle. The phase and polarization parameters therefore realise the Hopf fibration in the optical hypersphere (this is explained in detail in the Supplementary Information). The Poincaré sphere is interpreted here as the base space of the fibration 31 . We design a 3D structured beam which realises all the transverse states of light, including polarization and phase, in its focal volume (real space). It displays the 3D Hopf fibration topology in a configuration we call a skyrmionic hopfion. The skyrmionic hopfion realises, in real space, an image of the Hopf fibration in the optical hypersphere. The fixed polarization filaments can be represented as a 3D topological texture of entwined curves, in which each pair of loops are linked.

Experimentally realising the skyrmionic hopfion
We design the skyrmionic hopfion structure in light by superimposing carefully chosen combinations of vectorial Laguerre-Gauss beams LG ℓ,1 23,30 . The LH component, # $ , is chosen to be the Laguerre-Gauss beam LG *!,2 , with a negativesigned optical vortex along the beam axis 23 . The RH component, # # , is chosen as a superposition of the Laguerre-Gauss beams LG 2,2 and LG 2,! , with a circular vortex loop in the focal plane centred on the axis 23 . Therefore, the net polarization field has a RH C line along the axis, threading a LH C line loop in the focal plane. The C lines, at which / = 0, 8, organise the rest of the texture: between them are nested tori with / = constant, including the particular L surface of linear polarization at / = 8/2, analogous to vortices in other skyrmionic textures 8 . Details of the superposition optimisation are given in the Methods and Supplementary Information.
Experimentally, the RH and LH beam components are separately shaped by a spatial light modulator (SLM) (Fig. 2), before being combined in a joint beam path to shape the skyrmionic hopfion ( Supplementary Fig. 4). The total polarization state and phase of the resulting beam are measured at each point in the propagating volume via vectorial full-field reconstruction (VFFR, see Methods). For the VFFR, we combine established metrological techniques; namely, Stokes polarimetry, interferometry, and A light field with position-dependent transverse polarization and phase is created in a volume from the superposition of RH and LH circularly polarized beams, whose amplitude and phases are carefully structured. Spatial points characterised by the same state of elliptic polarization lie on filaments (a). The 3D polarization texture can be visualised by colouring the filaments according to the position of its polarization ellipse on the Poincaré sphere (b). The azimuthal angle ", representing the ellipse orientation, is coloured with the hues and the polar angle #, representing the ellipticity, is associated with the saturation levels. The sphere's poles, representing the circular polarized states, are black (LH) and white (RH). Each optical state also has a phase, represented by the position of the arrow along the polarization ellipse. In the transverse plane in (c), states of light are fully described by colours and arrowed ellipses. Along filaments of constant polarization, the phase on the ellipses varies smoothly, as shown in the insets for three representative planes. The VFFR measurements reveal the polarization Hopf fibration in the 3D light structure (Fig. 3a). The polarization ellipticity is constant on nested tori, made up of polarization filaments labelled by constant / and varying azimuth . (Fig. 3b-d). Our polarimetric resolution identifies these filaments clearly, particularly the linking between pairs of loops. This resolution compares very well with experimentally measured hopfion structures in other systems, such as cold atoms 11,12 and liquid crystals 10 . As predicted (see Supplementary Information), the two linked C lines (vortices in the superposed beams) are the topological skeleton of the hopfion structure, on which the rest of the polarization texture hangs. They are not topologically privileged-all polarization filaments are linked loops-but the C lines form the core filaments for the system of tori, including the L surface of linear polarization. We anticipate C lines to play a similar structural role in other topological 3D polarization textures.
Considering the shaped beams' phase as well as polarization allows a comparison of the measured hopfion structure in real space (Fig. 4a, with phases along the shown filaments in Fig. 4b) with the optical hypersphere (Fig. 4c), parametrised by ., /, 0. This direct comparison gives a volume-to-volume mapping (demonstrated by the grey cubes in Fig. 4a,c). The density of hypersphere volume with respect to real space volume is the topological Skyrme density F, which can be interpreted as a continuous measure of linking 33 of the polarization filaments. Characteristic of 3D skyrmions 2,5,8 , Figure 3. Visualising the topology of the focal volume. The optical texture is reconstructed from the polarization and phase measurements via the VFFR of the optical beam. The measured volume is coloured following the Poincaré sphere and reveals the topological structure of the Hopf fibration (a). Two C lines, the black loop and the threading straight white line, organise the texture into nested tori. Each toroidal surface represents points characterised by the same ellipticity. The colours wind nontrivially around each torus, and a few polarization filaments making up these tori are shown in the insets: in (b), the lighter surface (+ " = 0.398) is made of lines characterised by RH elliptic polarization; in (c), the L surface (+ " = 0) is made of lines along which the polarization state is linear 22,25 ; in (d), the darker surface (+ " = −0.775) is made of lines characterised by LH elliptic polarization. In each inset, the cyan and red filaments, corresponding to # = 0, 4 are shown to form a Hopf link. Every pair of filaments in the texture link in this way, consistent with the Hopf fibration.
the real space integral of F, concentrated around the C line loop, integrates to a value very close to unity, covering the hypersphere of hypersolid angle , i.e. a Skyrme number of 1. The Skyrme number is the degree of the mapping from 3D real space to the hypersphere, corresponding to the element of the homotopy group Π " . More details of this are provided in the Supplementary Information.
Mathematically, the Skyrme density F is the jacobian determinant of the map from real space to the hypersphere (see Supplementary Information), This is the natural 3D generalisation of the 2D topological density for 2D skyrmions 13,15 (here, full Poincare beams), # (& J K • (∇cos/ × ∇.). As the field parameters vary longitudinally as well as transversely, three parameters are needed to determine the full, continuous topological density determining the covering of the optical hypersphere, which is nonzero when the three gradient vectors are linearly independent. The topological density in equation (2) may be rewritten in terms of the normalized optical orbital current 29 2⇡ 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " c g X D z e C K 1 6 d x D n n L A J E C q I f H x V E = " > A A A B 7 X i c b V B N T w I x E J 3 F L 8 Q v 1 K O X R m L C i e w S E v V G 4 s U j J i 6 Q w E q 6 p Q u V b r t p u y Z k w 3 / w 4 k F j v P p / v P l v L L A H B V 8 y y c t 7 M 5 m Z F y a c a e O 6 3 0 5 h Y 3 N r e 6 e 4 W 9 r b P z g 8 K h + f t L V M F a E + k V y q b o g 1 5 U x Q 3 z D D a T d R F M c h p 5 1 w c j P 3 O 0 9 U a S b F v Z k m N I j x S L C I E W y s 1 K 7 3 E / Z Q H 5 Q r b s 1 d A K 0 T L y c V y N E a l L / 6 Q 0 n S m A p D O N a 6 5 7 m J C T K s D C O c z k r 9 V N M E k w k e 0 Z 6 l A s d U B 9 n i 2 h m 6 s M o Q R V L Z E g Y t 1 N 8 T G Y 6 1 n s a h 7 Y y x G e t V b y 7 + 5 / V S E 1 0 F G R N J a q g g y 0 V R y p G R a P 4 6 G j J F i e F T S z B R z N 6 K y B g r T I w N q G R D 8 F Z f X i f t e s 1 r 1 K 7 v G p V m N Y + j C G d w D l X w 4 B K a c A s t 8 I H A I z z D K 7 w 5 0 n l x 3 p 2 P Z W v B y W d O 4 Q + c z x / k Y o 6 f < / l a t e x i t >  Details are given in Supplementary Information. An analogous expression applies to 3D skyrmions in other systems 3,7 , with an appropriate current or velocity substituted. It is also the topological helicity, describing knotted fields in high energy physics 2 , superfluids 3,7 , magnetic fields and hydrodynamics 34 . Its appearance in equation (3) suggests a relation between the 3D Skyrme density of a polarization field and the Poynting vector of optical energy flow.
We determine the Skyrme density explicitly from the measured data, as shown in Fig. 4d. The sum over the measured voxels gives a Skyrme number of 0.945, which is less than unity since low intensities limit the measured volume boundary. The corresponding covering of the optical hypersphere, with the image of the real space measurement boundary, is represented in Supplementary Fig. 10. Rather than a smooth interpolation of the optical field measurements, this density is determined discretely from a simplicial cell complex of spherical tetrahedra in the optical hypersphere arising from the measured data points. Details of the technique and its implementation are in the Methods and Supplementary Information. The value of the Skyrme number of the theoretical field, with the same boundary, is 0.997, consistent with the experimental error.

Discussion
We have demonstrated the experimental construction of a 3D skyrmionic hopfion in the polarization and phase pattern of a propagating light beam. The Hopf fibration is realised in the natural polarization parameters from equation (1), a mapping from 3D real space to the 3D optical hypersphere, generalising the Poincaré sphere naturally by including phase.
Our experiment and analysis manifest several topological ideas not commonly emphasised in optics. Firstly, optical polarization fields in 3D can have topological textures, analogous to textures in condensed matter, high-energy physics, etc. This might lead to further insights and possibilities for topologically structured light and its applications. Secondly, as a parameter space for the full vectorial light field, the optical hypersphere goes beyond the standard Poincaré sphere. The usual approach requires a Pancharatnam-Berry phase 35,36 to be included later, ignoring the fact that the optical field parameters define a manifold as natural as the 3-sphere. It is intriguing to speculate whether the machinery of the Poincaré sphere analysis of polarization and Jones calculus may be cast in the optical hypersphere.
The 3D polarization Skyrme density F in equations (2) and (3) is a new tool to analyse optical vectorial full fields. F is the continuous topological charge density representing the abstract optical hypersphere volume covered by each real space point. F = 0 when the gradients of ellipticity, phase, and azimuth are linearly dependent, typically occurring along surfaces in 3D. The relation between F and the optical orbital current suggests a subtle interplay between the Poynting vector 29 and energy-momentum fluxes with optical hypersphere topology (explored further in Supplementary Information).
A smooth polarization texture is disrupted at point singularities in the polarization field, such as saddle points in the parameters ., /, 0. As previously observed 24 in the reconstruction of Seifert surfaces spanning knotted optical singularities, these points are experimentally hard to control and limit the effective reconstruction of textures of polarization lines. They do not affect the Skyrme number, and such points will lie on the surfaces F = 0.
Our experiments and theory demonstrate new topological invariances possible in structured light. This formulation and measurement of an optical Π " invariant will lead to new, robust topological design principles for 3D optical fields for free space optics and nanophotonics. These skyrmionic structures generalise to fields with higher degree Skyrme numbers, involving more complex superposed beams including knots and links, offering a broader gamut of topological structures and integers that can be encoded in structured optical beams. This approach to topological beam shaping will offer further analogies with cold atoms, condensed matter and high energy physics, offering the possibility of emulating, optically, exotic particle-like topologies from field theories not accessible otherwise in the laboratory.

Topological design of optical skyrmionic hopfion
The optical skyrmionic hopfion consists of the two scalar fields # # and # $ representing the right-and left-handed field components respectively. These scalar components are appropriately structured to give the 3D topological texture described in the main text, effectively realising the topological mapping from 3D real space to the optical hypersphere. In these methods, we will refer to unnormalized field amplitudes O # and O $ rather than their normalized counterparts: the beam intensity is P = |O # | % + |O $ | % , and As discussed in the Supplementary Information, the axial optical vortex line should have a negative sign, so we choose O $ = 2\LG *!,2 (U, V, W; Y), the simplest LG mode with an axial vortex of the correct sign, with \ a constant to be found and 2 included for calculational convenience. The vortex ring can be realised by the sum of two LG modes with ℓ = 0, O # = (−] +^)LG 2,2 (U, V, W; Y) − bLG 2,! (U, V, W; Y), where ] and ^ are parameters to be found. This guarantees the vortex ring to be in the focal plane W = 0, with a radius of `]⁄ Y, provided ],^> 0. The coefficients ] and ^ determine the intensity pattern around the vortex ring as well as its radius.
For fixed value of Y the optical skyrmionic hopfion is therefore realised for a range of values of ],^ and \. The different values of the parameters give very different shapes of the structure (residing in the polarization parameters) and distribution of the overall intensity P. A preliminary exploration of these is given in the Supplementary Information. The spreading nature of the gaussian beams means that it is not possible to cover the 3-sphere completely with polarization states realised in 3-space. We therefore choose a superposition which maximises the volume of optical states in the optical hypersphere within the measured 3D volume in real space. To be effectively generated and measured in the experiment, the values of the parameters are chosen to optimise the field configuration. To aid this optimisation, we introduce an extra scale size parameter b, with ^= ^2b % and \ = \ 2 b. The 3D size of the skyrmionic hopfion scales according to b, where now the vortex ring radius is ⁄ . The remaining parameters ],^2, \ 2 , determine the particle-like field distribution's shape. The parameters ],^2, \ 2 and b were chosen to ensure the experimental skyrmionic hopfion to be localised within the measured volume, in practice a cartesian cuboid centred around the focal point. We optimised against the criteria in the following list. (i) Vortex ring radius U 2 not larger than beam waist Y (this principle is also used in the design of optical vortex knots 37,38 ). (ii) Concentrate intensity inside the measured volume, with P ≈ 0 outside the measured volume. It was especially important to localise the intensity within the transverse cross-section, so as not to lose critical polarization information. (iii) Distribute the intensity as evenly as possible within the measured volume. To maximise the quality of the measured polarizations we avoided regions of low intensity as much as possible where the polarization state changes rapidly. (iv) Concentrate the Skyrme density (continuous topological charge density) within the measured volume, i.e. F ≈ 0 outside the measured volume. The density F is given in main text equations (2) and (3), and described in detail in the Supplementary Information. This enables a measured value of the Skyrme number very close to 1, as described in the remainder of the Methods.
We proceeded by making an estimate of the parameters based on the topological 3D plots of the numerical models, and then improved these based on the quality of the experimental measurements.
The experimental setup, as described below, requires the Fourier transform of the beam superposition to be realised on the SLM, and the desired field is mathematically back-propagated through the paraxial lens system using Fourier optics 39 The coefficients do not depend on the Fourier waist Y ℱ , so the overall beam in real space scales linearly in radius U and quadratically in propagation distance W as Y ℱ is varied. This quantity is chosen so that the skyrmionic hopfion has the desired size in real space whilst fully utilising the SLM. In our optical system, e = 532 nm, the waist of the beam on the SLM is Y ℱ = 6.252 × 10 *5 m, and the imaging system given by lenses L1 and L2 ( Supplementary  Fig. 4b) halves the size of the beam. The resulting waist width is Y = 54.2 µm (giving a Rayleigh range W # = 34.7 mm). The measured volume is a cuboid, |l| ≤ l <=> , |m| ≤ m <=> , |W| ≤ W <=> , with l <=> = 3.13Y = 170 µm; m <=> = 3.91Y = 212 µm; W <=> = 0.768W # = 26.6 mm. The values for the beam parameters were optimised in this range to be ] = 3, ^2 = 1.5, \ 2 = 0.16, b = 2.5. In terms of the original parameters, this gives the values ] = 3, ^= 9.4, \ = 0.4. With these choices, the LH C line ring is at U 2 = 0.57Y = 30.6 µm. The field configuration of this model field near the C line ring is shown in Supplementary Fig. 2b, resembling the corresponding Hopf fibration configuration (e.g. Supplementary Fig. 2a) closely.

Optical system design
The experimental skyrmionic hopfion field is the superposition of two structured beams of orthogonal circular polarization, O # and O $ . Experimentally, these two scalar components are shaped by the amplitude and phase modulation of a collimated laser beam (horizontal linear polarization, expanded) performed by a reflective phase-only spatial light modulator (SLM; Holoeye Pluto phase-only, 1920 ´ 1080 px HD display), shown in Supplementary Fig. 4b. The SLM is used in split-screen mode 40,41,42 , with each half embedding the amplitude and phase information of O # and O $ respectively. To optimise the beam quality, the Fourier hologram for each polarization component is a 600 ´ 600 pixels square. This resolution was proven to produce all details of the transverse beam structure in the focal volume. The two holograms are placed so each receives approximately homogeneous illumination of the expanded input laser beam without losing too much intensity. The phase-only hologram is shown in Supplementary Fig. 4a.
To allow for amplitude modulation by a pure phase hologram, a weighted blazed grating is applied 43 . The desired scalar modes appear in the first diffraction order, which is spatially filtered by an aperture A in the conjugate plane of the SLM, generated by lens L1 (shown in Supplementary Fig. 4b). Fourier holograms are applied on the SLM, so that the desired beams are sculpted in the focus of the Fourier lens (L1), i.e. in the conjugate plane of the SLM. The hologram for each beam is normalized separately, taking advantage of the full modulation depth of the SLM for each beam individually.
The two beams are subsequently combined on-axis by an interferometric system. Before they are combined, the two beams are given orthogonal linear polarizations by a combination of a half wave plate (HWP) and a polarizing beam splitter (PBS), allowing also for the adjustment of the beams' intensity ratio. This is a critical step to realise the complex polarization structure: the HWP angle directly affects the relative strength of the two components and hence the coefficient \ in the field design described above.
After the beams are combined, a quarter wave plate (QWP) transforms the orthogonal linear polarization states into orthogonal circular polarizations. The imaging system given by lens L1 and L2 ( Supplementary Fig. 4b) halves the size of the beam and L3 performs the final Fourier transform that gives the skyrmionic hopfion in its focal volume. The focal structure is magnified by lens L4 (16´) onto a CMOS camera (Cam; uEye SE (UI-1240SE), 1280 ´ 1024 px).

Volumetric full field reconstruction
We retrieve the full field information (transverse components of the paraxial beam) by reconstructing the polarization and phase in the focal volume. Supplementary Fig.  5 shows five transverse planes at different positions in the propagation direction for the normalized Stokes parameters < ! , < % , < " and the phases p # and p $ of the RH and LH field components. The measurements in multiple transverse planes are performed via digital propagation 32 (see Supplementary Information). A detailed description of polarimetry 44 (Supplementary Fig. 4c) and transverse phase interferometry 45 ( Supplementary Fig. 4d) can be found in the Supplementary Information. The polarization measurements across different planes are unaffected by the harmonic time dependence of the optical field and are directly stored into 3D arrays. However, when stacking volumetric phase measurements, the relative phase between neighbouring planes must be retrieved. First, we describe our procedure for connecting the transverse phase measurements to their neighbouring planes, and then we present our routine to minimise the experimental error in retrieving the field components.
The measured transverse phase structure per plane is constituted of the light field's propagation term, e i@A times the superposed LG structure described in the "Topological design" Method section above. This includes a Gouy phase factor e *iBC G , where p G (W/W R ) is the W-dependent Gouy Phase, and a phase term varying radially and longitudinally (full Laguerre-Gauss modes equation is given in Supplementary Information), and a time-dependent phase offset due to the time varying phase relation between the measured and reference beams. In order to concentrate on the transverse variation, we circumvent the effect of e i@A within the measurements per Wplane, thereby avoiding the effects of undersampling the electric field oscillation, by setting the distance between two transverse planes to a multiple of the wavelength (100e), so the propagation factor e i@A is negligible. Next, we choose a transverse reference point (" :ref , W) close to the optical axis (U ≈ 0), so that the phase at this point is only affected by the W-dependent Gouy phase term of the LG beams and is unaffected by the other spatially varying phase factors. For each plane, the phase of the reference point is set to the same value, so the Gouy phase and the timedependent phase offset are subtracted. In order to finalise the missing relation between different W-planes, the theoretical Gouy phase term p G is added. Note that the Gouy phase represents an offset value per plane, only depending on the W-position but without any dependence on the transverse coordinates. Thus, the measurements themselves are not affected by this approach and, as a result, we correct for the errors in W caused by the time-dependent variations in the measurement system. Supplementary Fig. 6 shows the l = 0 plane (longitudinal cut) of the theoretically expected (left) and the reconstructed (right) 3D phase structures of p # and p $ . This figure demonstrates that the reconstructed 3D phase distributions are consistent with the theoretical predictions.
Due to experimental errors (see Supplementary information), the singularities of the differences of the phase of the two field components p $ − p # (wrapped between −8 and 8) do not coincide with those of the polarimetrically-determined arctan(< ! , < % ) as the polarization and phase measurement are independent. Observations of the 3D structure of the C lines from the polarization measurements and the phase singularities from the phase measurements allow the systematic error to be minimised by shifting the polarization measurements until the C line loop coincides with the singular loop of p # . Moreover, the overall error is reduced by redefining the Stokes parameters < ! and < % as follows: To finalise the volumetric full field reconstruction we calculate the real and imaginary parts of the beam components from # # =`(q 2 + q " )/2 e iC * , and # $ =`(q 2 − q " )/2 e iC + . The full field is used to calculate the Skyrme density of the optical field as described in the next Method.

Numerical calculation of experimental Skyrme number
We measure the Skyrme number of the optical skyrmionic hopfion directly from the discretely sampled, measured data by taking advantage of the robustness of topology. This optimises the computational speed necessary to evaluate the Skyrme number from experimental measurements. The measured polarization and phase at each point in real 3D space correspond to a point in the optical hypersphere. The 3D cubic lattice of measured voxels is mapped into a topology-preserving but distorted lattice in the optical hypersphere. An example for the ideal hopfion field (see Supplementary Information), is shown in Supplementary Fig. 9a-b. The measured Skyrme number is therefore based on this piecewise-linear mapping generated from the measured data points without interpolation. This approach can readily be used for measurements of other physical Skyrme-like maps, including lower dimensional ones (e.g. via triangular meshes).
The fully resolved experimental data in the focal volume gives real space voxels forming a cuboidal grid. We are interested in the Skyrme density of the real space volume given by a cuboid with transverse size r : = 1.84Y = 99.4 µm and longitudinal size r ∥ = 0.768W # = 26.6 mm as defined in the Supplementary Information. Since the image of the cuboidal mesh covers the volume of the hypersphere, reducing the resolution maintains this filling. The numerical routine is made more time efficient by reducing the resolution to a cubic mesh of dimension 101 ´ 101 ´ 101 in physical space, centred at the focal point. The voxels are centred at points labelled by (s, T, t) with 1 ≤ s, T, t ≤ 100. Each such point corresponds to a normalized 4D vector u v⃗ = (Re# # , Im# # , Re# $ , Im# $ ) found via the VFFR method, giving a distorted cubic 3D grid in the optical hypersphere whose vertices are the points u v⃗(s, T, t). The distortion of the experimentally measured field is significantly greater than the example in Supplementary Fig. 9a-b, as can be seen in the images of the two real space planes in the optical hypersphere in main text Fig. 4. Each elementary cube x = x J,7,@ is labelled by s, T, t, with vertices u v⃗(s, T, t), u v⃗(s + 1, T, t), u v⃗(s, T + 1, t), ... denoted by c1, …, c8 as indicated in Supplementary Fig. 9c-d. The cube x occupies a volume Vol(x) within the optical hypersphere. We numerically determine Vol(x) as follows.
An elementary topological cell in 3D is a tetrahedron (i.e. a 3-simplex 46 , in the language of simplicial topology). We convert our cubic s, T, t lattice into a 3D simplicial complex by decomposing each cube into five irregular tetrahedra. The resulting mesh of tetrahedra, where neighbours share triangular faces, edges and vertices, make up a 3D cell complex 46 . The tetrahedra can share the cube's vertices in two distinct ways, which are given by the following ordered sets of four vertices (see For any cubic lattice, cubes can be decomposed into two choices of tetrahedral mesh: cubes of type A at positions where the quantity s + T + t is even (odd) and cubes of type B where s + T + t is odd (even). We compute both types of 3D cell complex as a check of numerical accuracy. As a result, the measurement points in real space and measured values in the hypersphere define a piecewise linear map representing the physical field.
In the hypersphere, the tetrahedra are constructed so the edges joining the vertices are geodesics. Each tetrahedron's four faces are spherical triangles, and along edges, pairs of faces meet at the dihedral angles 0 < { 7 < 8, for T = 1, 2, … , 6. The formula for the 3D volume Vol(}) of an irregular spherical tetrahedron } constructed in this way can be written explicitly in terms of dihedral angles by means of Murakami's formula 47 (see Supplementary Information). The contribution to the Skyrme number comes from the signed volumes sign[det(u v⃗ K , u v⃗ L , u v⃗ M , u v⃗ N )]Vol(}), where u v⃗ ℓ with ℓ = {],^, \, Ä} are 4D unit vectors pointing to the four vertices of a spherical tetrahedron }. Only tetrahedral cells included within a 3-dimensional hemisphere, whose volume is less than 8 % , are considered. The sign of the volume comes from the ordering of the vertices with respect to the right-hand rule, where the triangular base ],^, \ follows the fingers and the vertex Ä follows the thumb. When the volume is negative, the order of the vertices in real space and that of the vertices of the tetrahedron in the optical hypersphere are inverted. This follows the standard orientation rules of a 3-simplex.
At higher resolutions of the cubic lattice, the spherical tetrahedra are smaller, and the curved edges tend to become linear and the spherical distortion can be neglected: the tetrahedron volume are better approximated by its flat-space analogue. The 101´101´101 mesh defined above is consistent with the volume of the tetrahedra being within the range allowed by numerical precision.
The hyperspherical volume of the cube x corresponds to the union of the volumes of the associated, neighbouring tetrahedra comprising x, and its volume Vol(x) is the sum of the signed spherical tetrahedra volumes Vol(}). The results for each such Vol(x J,7,@ ) are stored in two 3D arrays, one for each type of 3D cell complex. The experimental Skyrme number is found by adding together the volumes of all the hypersphere cubes with appropriate sign, normalized by the 3-sphere volume: ∑ VolÉx J,7,@ Ñ J,7,@ /(28 % ) (see Supplementary Information). The measured 3D Skyrme number corresponds to the fraction of the hypersphere volume by the image of the measured volume of real space. The sums over all the elements in the arrays give Skyrme numbers 0.94521 and 0.94528 for the two kinds of mesh. We take the experimental Skyrme number to be the mean of these two numbers, 0.94524. It is straightforward to implement the vector calculation described here in a numerical algorithm in MATLAB or Python. The volumes of the cubes in the meshes can be calculated in parallel via high performance computers. FIG. S1. Polarization ellipse. A typical polarization ellipse, bounded by a rectangle. The ellipse has a major axis given by vector a and a minor axis vector b; their ordering is consistent with the right-handed ellipse orientation (i.e. by sign(a ⇥ b) · b z). a is chosen rather than its opposite a. The major axis of the ellipse therefore has length a and the minor axis has length b, so the bounding rectangle has side lengths 2a and 2b. The normalization 1 = |E R | 2 + |E L | 2 = a 2 + b 2 implies that the bounding box for all ellipses has a fixed diagonal length of 2. The orientation of the major axis of the ellipse (and hence the longer edge of the rectangle) is ↵/2, defined over a range ⇡/2 < ↵/2  ⇡/2. Signed by ellipse handedness, half the area of the bounding rectangle is ±2ab = cos = S3. The phase angle on the ellipse is determined by /2; geometrically, this corresponds to the ellipse auxiliary angle. To determine the full range of phase, ⇡ < /2  ⇡, it is necessary to pick a reference major axis direction ±a at which = 0.
Hopfion-like field sections of constant azimuth. a, Section of constant ' of the complex field U(X, Y, Z), representing the 3-sphere stereographically projected to infinite euclidean 3-space. The cylindrical radius p X 2 + Y 2 is on the horizontal axis, and Z is vertical. The black contours denote contours of |U|, and the white lines (and colors) contours of phase arg U. The phase is singular on the nodal line U = 0, on which X 2 + Y 2 = 1 and Z = 0. This is the reference circle of toroidal coordinates ', , ⌧ , and in fact = arg U, cosh ⌧ = 1/|U|, so white and black lines represent toroidal coordinates in this plane. b, Section of constant of the model optical field in real space defined in the Methods. The LH C line crosses the plane at its radius R0 = 0.57w . White contours and colors correspond to the axisymmetric phase R (arg E R ), and black contours correspond to the axisymmetric values of S3 = cos . The values of S3 = 0.775, 0, 0.398 chosen in the main text are shown in grey, otherwise the contours of S3 are at 1, 0.9, . . . 0.9, 1, increasing from the minimum at the C point. The field strongly resembles the Hopf fibration/toroidal configuration of a. c, Model optical field with larger R and z range plotted. The contours of S3 are 1, 0.99, . . . , 0.99, 1. The investigated range is indicated by the dark cyan lines showing that the skyrmionic hopfion is mainly inside this range. The Hopf fibration structure breaks down at the saddle points of S3 which are shown as black dots. d, Amplitude of the model optical field across the measured range (see Methods). The yellow line shows the real amplitude profile at the focal plane (z = 0).      of Equation (S2.7) overlayed by a cubic lattice in x, y , z, whose vertices represent measurement points, from 2 to 2 at each 0.5. The experimentally measured lattice has much greater resolution. b, Image of cubic lattice in the optical hypersphere with The hypersphere is represented in volume-preserving projection with center corresponding to u = 1, v = 0. This is the same projection used in main text Fig. 4b. c, b, A cubic cell is divided into five tetrahedra with the same vertices: four tetrahedra at corners of the cube, and an internal tetrahedron whose edges are diagonals of the cube faces. The two (mirror image) ways this can be done are shown. For the cubic lattice, alternating between a and b guarantees neighbouring tetrahedra always share faces (i.e. giving a 3D cell complex). b, The projection is now from a LH polarization state. The real space axial C line maps to the white loop, and the white surface encloses the complement of the filled volume.

Supplementary Note 1: Parametrising the Optical Hypersphere
In this note we discuss the parameters describing the optical state. Starting with the Stokes parameters and Poincaré sphere for polarization, we generalise by including phase to the 3D parameter space we call the optical hypersphere. This is the unit 3-sphere in 4-dimensional space. Its natural parametrisation by the polarization and phase parameters defines the Hopf fibration structure. The Poincaré sphere is recovered from the optical hypersphere as the 2-sphere base space of the Hopf fibration.
The state of light at each point in a transverse, monochromatic light beam, travelling in the z-direction, is specified by in terms of the right-and left-handed (RH, LH) circular polarization complex unit vectors Throughout this work, for mathematical convenience, we assume that the field E is normalized, Equation (1) of the main text describes E R and E L in terms of three angle parameters ↵, and , We begin by describing the standard Stokes parameters and Poincaré sphere formalism, and also describe the phase parameter. This leads to a discussion of 3-sphere geometry and the Hopf fibration, which is then related to the parametrisation of equation (S1.1) using the angles ↵, , .
where 1 = ( 0 1 1 0 ) , 2 = 0 i i 0 , and 3 = 1 0 0 1 . These quadratic combinations of the four fundamental real-valued field quantities, Re E R , Im E R , Re E L , Im E L , only involve information about the polarization state -they are independent of the overall phase, typical of measurable quantum expectation values and classical constants of the motion. With these definitions, it is not difficult to show that S 2 1 + S 2 2 + S 2 3 = 1. The resulting normalized 3-vector (S 1 , S 2 , S 3 ), the Stokes vector, lies on the surface of a unit 2-sphere, the Poincaré sphere [1,2]. Every polarization state (regardless of its phase) corresponds to a point on this sphere: RH circular (E L = 0) to the north pole (0, 0, +1), LH circular (E R = 0) to the south pole (0, 0, 1), linearly polarized states (S 3 = 0) to the equator, and so on. More general elliptical polarized states with RH (LH) circulation, with S 3 > 0 (S 3 < 0) occur in the northern (southern) hemisphere. This parametrisation of elliptic states on the Poincaré sphere is demonstrated in main text Fig. 1.
The three standard Stokes parameters description may be extended by two extra, alternative parameters (mathematically defined in a similar way to S 1 and S 2 ) [3,4] which determine the phase information, (S1.4) With S 3 , these also satisfy the normalization condition S 2 3 + S 2 4 + S 2 5 = 1. This means that an alternative sphere similar to the Poincaré sphere can be defined with cartesian coordinates S 3 , S 4 , S 5 instead of S 1 , S 2 , S 3 [3,4]. We will discuss the significance of this definition below.

Polarization ellipse and phase
Apart from the overall normalization, the optical state E is determined by the normalized, complex two-component electric vector (E R , E L ), or equivalently, the normalized 4-dimensional vector (Re E R , Im E R , Re E L , Im E L ). These parameters, together with the equivalent angles ↵, , used in the main text and defined in equation (1), are related to the geometry of the polarization ellipse. The time-dependent, monochromatic electric field E Re traces out an ellipse in the xy -plane in time t with period 2⇡/!: A simple way of visualising this ellipse is in terms of its major axis a and minor axis b; the normalization requires |a| 2 + |b| 2 = 1 and the axes are perpendicular, a · b = 0. The ellipse curve traced out by Re{(a + i b) e i !t } is the same as that traced out by E Re (!t); the only difference is the starting point on the ellipse, which is related to the phase parameter .
The phase of E is related to the overall phase factor at t = 0 multiplying the complex combination of major and minor axis vectors at time t = 0, i.e. E = e i /2 (a + i b). From the discussion above, this implies that E Re (!t = /2) = a. Thus the effective phase evolves with time. It is necessary to be careful with this definition, as described in the following.
The major and minor axes correspond, with respect to t, to maxima and minima of , (S1. 6) which has been simplified using the alternative Stokes parameters (S1.4). This is a maximum at !t max = + and a minimum at !t min = , where Since arctan is a periodic function with period ⇡, there is another maximum/minimum at ± + ⇡. We therefore have In fact, these give a simple expression for the original complex vector E, E = (a + i b) e i + , (S1.9) where + determines the phase factor multiplying the complex vector of major and minor axes, with a ⇡ ambiguity. The definition of the phase parameter correctly over the 4⇡ range, i.e. as 2 + or 2 + ± 2⇡, requires some care, as described below. With this in mind, we discuss each of the angles ↵, , used in the main text and defined in main text equation (1).
This half angle cannot distinguish between ±a: the orientation angle ↵/2 of the major axis of the ellipse is in fact a director, only defined over a ⇡ range, say ⇡/2 < ↵/2  ⇡/2, consistent with our definition of the range ⇡ < ↵  ⇡. This half-angle nature of ↵ implies that the azimuth angle on the Poincaré sphere changes twice as quickly as the ellipse orientation angle in real space, and has been subject to much comment in the literature. The factor of 1/2 on ↵ in equation (1) is naturally associated with spinorial parametrisation. Its range from ⇡ to ⇡ is consistent with its definition as the argument of the product E ⇤ R E L . The angle is the polar angle on the Poincaré sphere, (S1.11) It determines the shape and the handedness of the ellipse, and it can also be shown that cos = 2b z · (a ⇥ b) = 2|a||b| sign S 3 . (S1.12) Therefore, cos is proportional to the area of the bounding rectangle of the ellipse, with the normalization condition |a| 2 + |b| 2 = 1 corresponding to fixing the main diagonal of the bounding rectangle, as shown in Supplementary Fig. 1. Since 1  S 3  1, the range of is from 0 to ⇡.

Optical phase and
The definition above of the phases ± at the maxima and minima ignored the subtlety that we did not distinguish between +a and a. In our previous discussion, we saw that this only determines the ellipse major axis as a director ±a, and the azimuth parameter ↵ (with ⇡ < ↵  ⇡) is twice the angle between the major axis and b x. In the definition of + in equation (S1.7), the choice of + = 1 2 arg(E R E L ) rather than 1 2 arg(E R E L ) + ⇡ was arbitrary, and is defined ⇡/2 < +  ⇡/2. This choice represents the phase difference between the point on E Re (0) on the ellipse and the closer of the two major axis vectors ±a. As the phase cycles around the whole ellipse, which one of ±a is closer to E Re ( /2) jumps, as E Re crosses the minor axis ±b when + = ±⇡/2. To determine the phase value unambiguously, it is necessary to fix a choice of major axis vector ±a for each polarization state, say +a. This allows the phase position of E Re (0) to be determined with respect to a, which extends the range of + to a full cycle, ⇡ < +  ⇡.
Just as the angles ↵ and are defined as the half-angles of the natural physical angle, so is the phase parameter twice the value of the "physical" phase angle, = 2 + . This implies that varies over a 4⇡ range, 2⇡ <  2⇡, with the general definition (S1. 13) In practice significant care is needed numerically to determine the correct branch cut of when E Re crosses a. A numerical value for can be determined nonlocally for instance by tracing the evolution of E Re with respect to the chosen sign of the major axis vector.  (2)). This double covering, in the ↵, , parametrisation, corresponds to the 2⇡ versus 4⇡ range of [6].

The Poincaré sphere and its projections
The approach to polarization and the optical state (E R , E L ) in the circular bases is very reminiscent of the mathematical spinor formalism. This is most familiar in the representation of the state of a quantum spin 1/2 particle, where the spin direction in 3D real space is represented by a point on the Bloch sphere.
Physically, a monochromatic electric field has spin 1. Since we only consider a 2-component transverse field, the complex field (E R , E L ) mathematically represents a normalized complex spinor. The optical state is not truly a spin 1/2 system, therefore, and so (E R , E L ) is properly a pseudospinor. The corresponding Poincaré sphere is defined in an abstract 3D space, unlike the Bloch sphere, whose directions correspond to those in 3D real space.
An alternative representation of polarization states is via the complex ratio (S1.14) which defines a point in the complex plane with modulus (radius) tan( /2) and azimuth ↵. This point corresponds to the stereographic projection of the Stokes vector, centred at the north pole, and including the point in the plane "at infinity" ( = ⇡, i.e. the south pole), the plane is topologically equivalent to the sphere. This is a direct illustration of one-point compactification: a euclidean plane, with a single point 1 adjoined, becomes, topologically, identical to a sphere. This procedure can be generalised directly to euclidean n-space and the n-sphere.
In principle, any positive monotonic function ⇢( ) can be used as a radial coordinate to map the 2-sphere to the plane; when ⇢( = ⇡) is finite, the map is in a disk of finite radius, and all points on the circumference ⇢(⇡) correspond to the south pole. In the example of stereographic coordinates above, ⇢ = tan( /2); this projection from the sphere to the plane is conformal, preserving angles. Another important example is ⇢( ) = 2 cos( ), Lambert equal area projection, where areas on the plane correspond to areas on the sphere. Areas are preserved here even though the south pole = ⇡ projects to a circle (radius 2).
The stereographic representation (S1.14) is based on a north pole projection, which associates the north pole = 0 with the origin of the plane, and the south pole = ⇡ with 1. With this choice, positive orientation (the local righthanded sense on the sphere) is consistent with the orientation of the increase of the sphere azimuth coordinate ↵. We could instead have chosen south pole projection where ⇢ = cot( /2), which associates the south pole with the origin and the north pole with 1. For this projection to preserve orientation, positive orientation is associated with a decrease in sphere azimuth angle ↵, so the complex coordinate associated with orientation-preserving south pole projection is cot( /2) e i ↵ . This is important for the analogous argument for the 3-sphere.
1.2. The hypersphere and the Hopf fibration 1.2.1. Parametrising the hypersphere Many of the details discussed above about 2-sphere geometry, topology and parametrisations generalise directly to the 3-sphere (hypersphere) [7,8]. For instance, we can consider 4-dimensional (4D) space with cartesian coordinates W, Z, X, Y . The reason for the non-standard ordering is as follows.
The 3-sphere is defined by all possible 4D unit vectors, represented using overarrowñ = (n W , n Z , n X , n Y ). Points on the 3-sphere can be parametrised by hyperspherical polar angles !, ✓, ' with 0  !, ✓  ⇡, ⇡ < '  ⇡ in the following way, (S1.15) Clearly these satisfy n 2 W + n 2 Z + n 2 X + n 2 Y = 1. (Note that the 4-space azimuthal angle ' is different from the 3D real space azimuthal angle used later.) Features of the topology of the 3-sphere can also be discerned from this parametrisation, by relating the XY Z coordinates to the 3D cartesian coordinates. n X , n Y , n Z resemble the coordinates of a 2-sphere, parametrised by spherical polar angles ✓ and ' as normal, and with radial coordinate sin !: all these spheres enclose a central point (at which ! = 0), and extend out to a sphere of unit radius at ! = ⇡. In this embedding, the boundary sphere corresponds to the single point ! = ⇡ in the 3-sphere (similar to the embedding of the 2-sphere in the 2-plane as discussed above). This solid ball in 3D is therefore a projection of the 3-sphere into euclidean 3D space, for which sin ! = ⇢(!) as defined in Sec. 1.1, related to the quaternion representation of rotations (where the angle of rotation is 2!) [9,10].
Just as for the 2-sphere, we can consider different "ball" projections of the 3-sphere into 3D euclidean space by different choices of ⇢(!). For instance, stereographic projection to X, Y, Z coordinates is achieved via ⇢(!) = tan 1 2 !. Explicitly, this gives (S1.16) As in 2D stereographic projection (S1.14), the 3-sphere point ! = ⇡ corresponds in euclidean 3D space to the point "1". We can also write complex coordinates U, V for the 3-sphere, (S1. 17) where in each case the second equality corresponds to stereographic projection to euclidean space (with the point with ! = ⇡, n W = 1 at infinity), and |U| 2 + |V | 2 = 1. The modulus and phase of U are shown in Supplementary Fig. 2a; as described below, with ' these correspond to toroidal coordinates in euclidean 3-space. Stereographic projection provides a very natural visualisation of the hypersphere in terms of familiar, infinite 3D euclidean space. Yet other projections of the 3-sphere into euclidean 3-space are possible [9,10]. The representations of the optical hypersphere in main text Fig. 4, are chosen to preserve 3D volume, analogous to Lambert equal area projection, for which the radius function is ⇢(!) = [ 3 4 (2! sin 2!)] 1/3 [9].

Hopf fibration
We have considered the polar angle parametrisation of the hypersphere using !, ✓, ' by analogy with the polar angles of the 2-sphere. They are not the only natural parameters of the 3-sphere, and in fact are not the natural physical parametrisation of the optical hypersphere. We identify the complex scalar U with the RH component of the optical field, i.e. U = E R , and V with the LH component, V = E L ; this is possible since both pairs of complex scalars are normalized, 1 = |U| 2 + |V | 2 = |E R | 2 + |E L | 2 . This identification gives the optical state parameters ↵, , as the parameters of the Hopf fibration in the optical hypersphere, described in the following.
Firstly, the identification of U and V with E R and E L allows us to find the relationship between ↵, , and ✓, ', !, (S1.18) The unit 4-vectorñ is identified with the 4-vector of real and imaginary parts of E R and E L considered earlier, = (cos !, sin ! cos ✓, sin ! sin ✓ cos ', sin ! sin ✓ sin ') .
We now interpret ↵, and in the 3-sphere, represented by euclidean 3D space under stereographic projection. This is possible since the Hopf fibration of the 3-sphere can in fact be constructed directly from toroidal coordinates [11][12][13] in euclidean 3-space, considered as the stereographic projection of the 3-sphere.
Consider the complex coordinates U, V as complex scalar functions of 3D euclidean space. U has a circular nodal line (phase singularity) of unit radius in the XY -plane centred on the origin, which we call the reference circle (anticipating the use of toroidal coordinates below). This is threaded by a nodal line in V along the Z-axis (which corresponds to a circle in the 3-sphere, with one point mapping to 1 in euclidean 3-space), which we call the axial line. As a stereographic representation of the optical hypersphere, the states of LH polarization ( = 0) occur on the reference circle, and the states of RH polarization ( = ⇡) on the axial line (this includes the point at 1, closing topologically into a loop). The phase varies smoothly along these lines. Representing the RH and LH states of polarization, the axial line and reference circle are special, as curves in the 3D optical hypersphere, corresponding to the points of circular polarization at the poles of the 2D Poincaré sphere.
We now construct the Hopf structure using the fact that U and V are naturally related to toroidal coordinates ', , ⌧ in X, Y, Z space [11][12][13]. The toroidal coordinate ' is the familiar azimuth angle of spherical and cylindrical polar coordinates, ' = arctan(X, Y ) = arg V = 1 2 ( + ↵), (S1.20) from equations (S1.17) and (S1.18). We use the standard convention of increase of ', but note it needs to be reversed when considering 3D orientation. The coordinate is the poloidal angle about the reference circle U = 0 (also see Supplementary Fig. 2a), and in fact, also from equations (S1.18) and (S1.17), On the punctured plane Z = 0, R > 0, = 0, and increases around the reference circle in a RH sense in a plane of fixed ', as in Supplementary Fig. 2a. The final toroidal coordinate, ⌧ , is defined over the range 0  ⌧ < 1, by which is simply a monotonic function of . Surfaces of constant ⌧ are tori centred on the axial line which enclose the reference circle, with major radius coth ⌧ and minor radius cosech ⌧ . ⌧ = 0 on the axial line, with ⌧ ! 1 on the reference circle. Now, the azimuthal and poloidal angles ' and are angles winding around the tori of constant ⌧ in flat planes; their sum and difference, ↵ = ' and = ' + , both wind around the tori as (1, 1) torus curves. Loci of constant and ↵, therefore give curves-in fact Villarceau circles-which wrap around both nodal lines in U, V as . The resulting structure of lines on which ↵, are constant define the Hopf fibration, as in Supplementary Fig. 3. Position on these Hopf circles is parametrised by on its 4⇡ cycle.
The Hopf fibration represents the 3-sphere as a fiber bundle, whose fibers are the Villarceau circles for fixed ↵, along which varies. The parameters , ↵ themselves are determined by points on the Poincaré 2-sphere (the base space of the fibration). The loci of constant (i.e. constant S 3 ) correspond to constant ⌧ , i.e. tori, around which are the nodal line loops corresponding to = 0, ⇡.
The choice of ↵, parametrising the base space and parametrising the fibers is only one possibility; it is also possible to parametrise the fibers with ↵, in which case the base space is the sphere parametrised by ↵, (with appropriate modification to the ranges). The sphere of this fibration is via the Hopf map of the alternative Stokes parameters S 3 , S 4 , S 5 introduced in equation (S1.4).

3-sphere volume
The unit 3-sphere has total volume (hypersurface area) 2⇡ 2 , as we now discuss. The normalized volume element d 3 ⌦ of the hypersphere can be derived from the hypersurface area of the 3-sphere in 4D space, in terms of either set of angle coordinates ✓, ', !, or ↵, , . Some care is needed to ensure the proper sense of orientation on this 3-manifold; we define this in terms of the usual orientation of euclidean 3D space for stereographic coordinates, which corresponds to ordering the angles !, ✓, '.
In terms of !, ✓, ', the normalized volume element (dividing by the hypersphere volume 2⇡ 2 ) is where det[ñ, @ !ñ , @ ✓ñ , @ 'ñ ] denotes the matrix of column vector partial derivatives ofñ, calculated using equation (S1. 19). This determinant plays the role, in 4D space, of the scalar triple product in 3D. With this normalization, the integral of these angles over the 3-sphere with ranges 0  !, ✓  ⇡, ⇡ < '  ⇡ indeed gives unity. We can apply this procedure to the optical state parameters ↵, , . In order to preserve the sign of the determinant (i.e. the 3-sphere orientation), it is necessary to use the ordering , , ↵ (or a cyclic permutation of them), i.e.
(S1. 24) In the last equality, was replaced by cos ; the volume element in terms of ↵, cos , has no prefactor, and is positive in lexicographical ordering (since cos increases as decreases).

Supplementary Note 2: Designing the optical skyrmionic hopfion
In this note we summarise the theoretical model of the skyrmionic hopfion field. This field configuration captures the Hopf fibration structure in real space, realising all phases and polarizations within the propagation volume. We construct solutions to the (paraxial) optical wave equation that show the desired topology, then we describe how we select the optimal superposition of Laguerre-Gauss modes representing a physically realisable optical skyrmionic hopfion.
Such a field is designed by creating RH and LH fields, as functions of real 3D space using cartesian coordinates x, y , z or cylindrical coordinates R, , z, i.e. E R (x, y , z) and E L (x, y , z). We introduce unnormalized forms of the field components R (x, y , z), L (x, y , z), such that The scalar fields = R , L are solutions of the paraxial wave equation r 2 ?
is the 2D transverse laplacian, k = 2⇡/ is the wavenumber; these are good representations of scalar components of collimated laser beams [14]. For future consideration, we define the phases of the scalar components, R ⌘ arg R , and L ⌘ arg L . As described in the main text and methods, the skyrmionic hopfion field is designed by a superposition of Laguerre-Gaussian (LG) beams, with the RH component R being zero on a circle in the focal plane, and R zero on the beam axis. The zeros of the circular components correspond to C lines, along which the state of polarization is circular [15]. In the full polarization field, these correspond to a RH C line along the axis, and a LH C line circle. The design of the optical skyrmionic hopfion field is thus based on the principle that C lines are the skeleton of the complex optical polarization field, although careful tuning of the superposition is necessary to guarantee the covering of the optical hypersphere within the focal volume.

Laguerre-Gauss modes
Laguerre-Gauss (LG) modes are a convenient set of basis functions for solutions of the paraxial wave equation. The LG modes are naturally expressed in cylindrical coordinates R, , z by [14,16] LG 3D ,p (R, , z; w ) = Here, w is the width of the gaussian envelope setting the transverse size of the beam, and z R = kw 2 , the Rayleigh range of the gaussian beam, determines the longitudinal lengthscale. L |`| p denotes an associated Laguerre polynomial [12,14]. The positive or negative (or zero) integer`is the azimuthal quantum number; when`6 = 0, there is an optical vortex along the axis of topological strength`. The nonnegative integer p is the radial order, corresponding to the order of the Laguerre polynomial, giving the number of nodal rings concentric to the axis in a typical transverse plane. The term LG`, p in the main text corresponds to LG 3D ,p (R, , z; w ) for some fixed value of w (and of course k). Although the equation (S2.2) is more useful for calculations, the form (S2.3), in which complex factors have had their modulus and phase separated, is easier to interpret physically. In particular, the LG beam has two important z-dependent phase factors: the Gouy phase e i(|`|+2p+1) arctan(z/z R ) , and the wavefront curvature e i R 2 zz R /w 2 (z 2 R +z 2 ) . Neither of these contribute in the waist plane z = 0.
For each value of z (i.e. transverse plane), the LG modes are an orthonormal basis for 2D complex amplitude distributions. The simplest such plane is the focal plane z = 0, for which As the fields are synthesised by 2D holograms in a beam's transverse plane, we approach the 3D field design primarily in terms of the field distributions in the 2D field [17]. Following standard linear space theory, any scalar solution (R, , z) of the paraxial wave equation can be written as a superposition where the generally complex coefficients are found from the 2D inner product integral

Mapping from 3D real space to the optical hypersphere
Any transverse optical field defines a mapping from 3D real space, parametrised by cartesian x, y , z or cylindrical polar R, , z, to the optical hypersphere, which is parametrised in terms of X, Y, Z, !, ✓, ', or ↵, , as discussed in Sec. 1.2, each of which are different ways of expressing the complex components E R = U, E L = V . In this section we will be careful distinguish the coordinates in the target space (i.e. the optical hypersphere), and the topological mapping to them from real space, as the notation of the complex fields E R (x, y , z) and E L (x, y , z) suggests both.
The most direct such map would be an optical field which maps each point of real space to the corresponding point in the stereographic projection of the optical hypersphere, i.e. (x, y , z) ! (X, Y, Z). Under such a map, each point of real space would correspond to a unique optical state (including the limiting point "at infinity"), wrapping around the optical hypersphere exactly once. However, such a mapping is not compatible with a physically realisable optical field, satisfying the paraxial equation.
As a step towards this, we define the following real-space counterparts to U and V of equation (S1. 17) u(x, y , z) = x 2 + y 2 + z 2 1 + 2 i z x 2 + y 2 + z 2 + 1 = R 2 + z 2 1 + 2 i z R 2 + z 2 + 1 , v(x, y , z) = 2(x + i y ) In the correspondence between x, y , z and X, Y, Z, Re u corresponds to Im U, as we explain below.
If we could realise a field with E R = u(x, y , z) and E L = v (x, y , z), we would realise a reflection of the skyrmionic hopfion (reflected in the n W = 0 plane), but as stated above, this is not possible physically with solutions of the paraxial wave equation. The vortex function 1 (x, y , z) = 2(x ± i y ), formally satisfies the paraxial equation, and can be taken as a local expansion of a physically realisable L giving the right form agreeing with the numerator of v . Furthermore, the function 2 (x, y , z) = x 2 + y 2 1 + 2 i z /k also satisfies the paraxial wave equation, so can represent a quadratic expansion of a physical field [18]. It has a circular vortex line around the axis in the z = 0 plane, and in fact agrees with the numerator of u when z is small.
The field 1 satisfies the paraxial wave equation, although its conjugate ⇤ does not; u was defined in equation (S2.7) to approximately agree with 1 , in spite of it mapping to a reflection of the hypersphere. Positive orientation can be restored to our mapping by choosing to generate ⇤ 2 in the other field, rather than 2 , i.e. a beam with a negative, rather than positive, axial vortex. This guarantees that the covering of the hypersphere by the real space beam has a positive sense, and the net mapping degree-the experimental Skyrme number-should be close to +1.
Therefore beams with L = ⇤ 1 and R = 2 , realise the correct topological mapping from real space to the hypersphere. Although not perfect realisations of the 1 1 mapping between the two 3D spaces, the physically realisable optical field we will describe is a smooth deformation of these, homotopically equivalent to the ⇧ 3 mapping around the hypersphere, except a small neighbourhood near 1 as we truncate within a finite experimental volume.

Topological design via C line skeleton
We thus design the structured fields E R (x, y , z) and E L (x, y , z) by matching the topology of R (x, y , z) and L (x, y , z) to the complex fields 2 (x, y , z) and ⇤ 1 (x, y , z), from Sec. 2.2 as closely as possible in the focal volume. The main design principle for the topology of the skyrmionic hopfion is therefore a vortex loop around the beam axis in the RH component (which corresponds to a LH C line loop in the full field), threaded by a straight vortex line along the z-axis in the LH component (corresponding to a RH C line).
As with other topological singular configurations designed in 3D [17,19,20], each of the RH and LH components is considered as a superposition of LG modes. The simplest way to generate a vortex ring in R is in the focal plane, by a superposition of LG 0,0 and LG 0,1 . Thus, for this component all coefficients (S2.6) are zero except c 0,0 and c 0,1 , which we redefine them in terms of new parameters a, b as: Therefore, the RH field is Here, the coefficients a and b tailor the radius of the vortex ring, the shape of the vortex core, and the related intensity. The superposition is chosen so the optical vortex ring has radius R = R 0 = p a/bw (a, b > 0) in the focal plane z = 0. In order to get a positive Skyrme number, it is necessary to use a negatively signed axial vortex in L . This is equivalent to matching the topology of v ⇤ (x, y , z), necessary to guarantee the mapping preserves the sign of 3D volume is preserved; more details are in Supplementary Note 4. For L , therefore, all c`, p are zero except for c 1,0 , which we choose to be involving the new constant c > 0. Thus, the LH component is given by This expression outside the square brackets has the same nonzero prefactor as the expression for R . From the perspective of the polarization, all that matters are the parts in the square brackets of equations (S2.9) and (S2.11).
Given fixed values of k and w , the optical skyrmionic hopfion forms in a parameter range of a, b, c. The criteria by which the parameters are optimised for experiments are reported in the Methods. When a = b = c = 1, the terms in the bracket become R 2 /w 2 1 and 2R/w e i and correspond (up to scaling), respectively, to the numerators of the complex 3-sphere coordinate fields u and v ⇤ from equation (S2.7) [21]. However, the z-dependence of the propagating light fields is different from the full u, v (v ⇤ ) fields. This is similar to the design of optical vortex knots [17], whose design is based on more general complex scalar fields (e.g. "Milnor polynomials") which have the same amplitude in the focal plane and same global 3D topology, although the functional form is different.

Supplementary Note 3: Optical generation and detection details
In this note we give an overview about the digital propagation method that is used to generate transverse planes at different propagation distances around the focus on the camera, the phase and polarization measurement technique and discuss the error rate on basis of their possible origins in the setup.

Measuring the optical skyrmionic hopfion
Digital propagation is an approach that allows for the investigation of different z-planes of a 3D structured paraxial light field without mechanically moving the observation plane. The field of interest itself is propagated by varying the Fourier hologram on a SLM [22,23] as follows. In terms of the Fourier representation, the scalar 3D light field with transverse position r ? physically propagating in the +z-direction can be represented by (r ? , z) = F 1 (F( (r ? , 0))e ikz z ) with F 1 denoting the inverse Fourier transform and, in the paraxial approximation, k z ⇡ k (k 2 x + k 2 y )/(2k). Since we display a Fourier hologram F( (r ? , 0)) per scalar mode on the SLM for generating the skyrmionic hopfion structure in the focus of L3 (see Supplementary Fig. 4b), we are able to digitally propagate the differently polarized scalar fields (r ? , z) and, thereby, the total topological texture in the z-direction, by multiplying the Fourier holograms at different transverse planes with the according phase factor e ikz z . For our analysis of the skyrmionic hopfion structure, we propagate in the range of z (see Methods) in 100 steps of 100 each, for investigating the phase as well as polarization per z-plane.
For the polarization analysis per transverse plane, a rotating QWP and fixed horizontally aligned polarizer (Pol) are mounted in front of a camera (see Supplementary Fig. 4c). Depending on the rotation angle ⇥ of the QWP, the camera records spatially resolved intensity patterns described by I(⇥) = 1 2 (s 0 + s 1 cos 2 (2⇥) + s 2 cos(2⇥) sin(2⇥) + s 3 sin(2⇥)) [24], where s 0 , s 1 , s 2 , s 3 are the unnormalized Stokes parameters, s 0 is total intensity, s j = s 0 S j , for j = 1, 2, 3. To suppress time-dependent errors, every intensity profile is the arithmetic mean of 10 intensity patterns. The Stokes parameters are measured recording the intensity patterns within the transverse plane of interest for 15 different angles ⇥. The resolution of the measurement is given by the camera (1280 x 1024 px), therefore the transverse polarization profile is determined across a single plane [25,26]. To use the full resolution of the camera, the beam structure is magnified by lens L4 with a factor of 16 onto the camera. The combination of polarimetry across transverse planes with digital propagation (see above) allows the full detection of the 3D polarization texture in the focal volume.
In order to retrieve not only the polarization but also the phase information from the realised 3D vectorial light field, in addition we implement digital holographic phase measurements [27]. This method allows for the determination of the transverse phase structure in a fixed z-plane of one orthogonal component (scalar field) of the beam relative to a reference beam. For phase metrology, we coherently interfere the scalar beam (signal beam), (r ? , z), off-axis with a reference beam, ref (r ? ) = | ref |e ik R R (plane wave). Both beams are of the same circular polarization. The reference beam illuminates the camera at an angle to the signal beam (see Supplementary Fig. 4d), generating an interference pattern whose intensity is I(r ? ) = | (r ? ) + ref (r ? )| 2 . The mixed term | || ref |e i( k R R) embedded in this interference pattern contains the desired phase information (r ? ). This term is extracted by a suited low-pass filter located at k R , the k-vector of the reference beam. By centering and performing an inverse Fourier transform of this cropped pattern we receive the pure phase profile of the signal beam. Similar to polarimetry, the beam is magnified and the intensity measurements are averaged.

Errors in measurements
In general, various systematic errors might occur within the experiment. For instance, these may include time dependent fluctuations within the system for the generation and measurement of the skyrmionic hopfion structure, and constant inaccuracies of optical components. We observe imperfections of the beam's intensity profile even before the beam is reflected on the SLM. These imperfections are from the laser source itself and the optical components that guide the beam to the SLM. To address this issue, a beam cleaning is installed immediately before the SLM to provide the best possible profile. During its way through the setup, the beam profile is affected by polarizing optics (wave plates and polarizer) and imperfections in the SLM. We minimize these systematic errors by precise adjustment, implementation of high-quality optical components, as well as phase profile corrections [28].
Within the measurement procedure of our VFFR method, minor deviation per pixel for polarization and phase measurements can be observed. The average error per pixel within the transverse plane is experimentally determined to be < 4% and < 2% for polarization and phase measurements, respectively. Since we use the arithmetic mean of ten measurements for each picture, this errors represent the standard deviation within ten measurements at the same (r ? , z)-position within the beam. Note that this average error is calculated from the errors per pixel, which vary spatially. As representative examples, Supplementary Fig. 7 shows a density plot at the focal plane (z = 0) of the deviation of the polarization component S 3 of the skyrmionic hopfion structure (a) and the phase L (b).
Errors particularly affect the polarization measurements around the nodal line of E R (i.e. the LH C line loop). The superposition of E R and E L in this area is very sensitive to small fluctuations in the intensity and phase of one of the beams. Slight fluctuations cause rotations of the polarization profile. However, observed deviations in polarization do not hinder the formation of the topology. On the other hand, the phase structure reveals a larger error close to the optical axis, around the phase singularity in E L (i.e. the axial RH C line). This effect is due to the higher phase gradient in this area, which causes small spatial vibrations in the experimental setup to have a significant effect on the phase variations between repetitive measurements. In spite of this, the phase singularity propagates close to the optical axis, as visible in Supplementary Fig. 6. Thus, it is in a good agreement with the theoretical prediction.
Any errors here are inherited in the determination of the experimental Skyrme number (see Methods). Despite visual deviations between the experimental and the theoretical structure in light and described errors, we were able to realise an optical skyrmionic hopfion with a small deviation of Skyrme number to the theoretically determined value (see Sec. 5.2).

Supplementary Note 4: Skyrme density, helicity and optical currents
Here we describe in more detail the Skyrme density, the continuous topological charge density which represents the covering of the optical hypersphere by the position-dependent phase and polarization in the optical field volume. Starting with its mathematical definition as the jacobian of the transformation between 3D configuration space and the optical hypersphere, we then put it into a more physical context related to singular optics and topological helicity. This extends to a discussion of the chosen volume over which the Skyrme number is calculated, and then further theoretical features of the model.

Skyrme density as the jacobian of the mapping from real space to optical hypersphere
We describe the basic derivation of the Skyrme density as the jacobian of the topological ⇧ 3 mapping from the experimental real space volume to the optical hypersphere. This is motivated by the analogous 2-dimensional mapping to the Poincaré sphere.
As a 2-surface, it has a standard "east-north" or "south-east" orientation. This corresponds to the sense of increase of azimuth ↵ in the northern hemisphere, and its opposite in the southern, as discussed in Sec. 1.2. This means that a mapping from the x, y -plane to the Poincaré sphere acquires a jacobian factor, which we call the 2D Skyrme density where a normalizing factor of the area of the 2-sphere has been included on the right-hand side. Therefore, including the orientation, ⌃ (2D) = @ x ↵@ y (cos ) @ y ↵@ x (cos ) = b z · r(↵) ⇥ r(cos ). From the second equality, we define the baby skyrmion density vector Thus ⌃ (2D) = ⌃ b z , and in general the density of covering the Poincaré sphere in any plane in 3D, with normal vector b u, can be found using b u · ⌃ b . This can be compared to the continuous 2D topological charge density in full Poincaré-type beams [29][30][31] and baby skyrmions in other systems [32][33][34][35][36][37][38].
Full Poincaré beams are propagating beams with position-dependent polarization such that the integral of ⌃ b z in each transverse plane integrates to an integer; this is guaranteed when, as x 2 + y 2 ! 1, the polarization state tends to a constant independent of direction. This effectively allows "point compactification" of the polarization states in the plane, and thus a topological mapping from 2D real plane to 2-sphere [39].

Skyrme density and covering the optical hypersphere
A skyrmionic hopfion is a mapping from 3D real space to the optical hypersphere, where the position-dependent optical polarization and phase parameters wrap around the hypersphere. We can find the corresponding 3D Skyrme density by analogy with the 2D case.
We therefore are considering a mapping from 3D real space with cartesian coordinates x, y , z to the optical hypersphere with polarization parameters ↵, and phase parameter (defined in equation (S1.13)). We discussed previously, in Sec. 1.2, that the 3-sphere of unit radius has total volume 2⇡ 2 , and as derived in equation (S1.24), the normalized volume element of the 3D hypersphere is (2⇡ 2 ) 1 d 3 ⌦ = (16⇡ 2 ) 1 d↵ d(cos ) d , with positive volume orientation given by (↵, cos , ). Therefore, by analogy with equation (S4.1), the 3D Skyrme density ⌃ is the jacobian where again the right hand side is normalized so the integral over the whole 3-sphere gives unity.
It would also be possible to use 3-sphere polar angles ', ✓, ! as in equation (S1.18) to define the Skyrme density, but the polarization parameters ↵, , are physically more natural to describe the optical state. These latter parameters are those of the Hopf fibration of the 3-sphere: a small interval (cos ) = S 3 (corresponding to a thin 'vertical' slice of the Poincaré sphere) corresponds to a thickened torus in the 3-sphere, on which curves of fixed and fixed ↵ wind in opposite directions as Hopf circles.
In a skyrmionic hopfion, the ↵, , parameters wind around their angles in a sense preserving overall positive orientation, guaranteeing the integral of ⌃ over 3D real space entirely covers the optical hypersphere an integer number of times, defining the Skyrme number.
This full 3D skyrmion is similar to topological textures in other systems. A Shankar monopole, for instance, is a texture in fields characterised by a 3D orthogonal frame, i.e. an orthogonal 3 ⇥ 3 matrix [40][41][42]. The space of such matrices, the Lie group SO(3), is doubly covered by the 3-sphere; in terms of our parameters, closes after a 2⇡ cycle, not 4⇡.
Other textures of interest are pure hopfions, characterised by mappings from 3D real space to the 2D sphere (such as fields of unit vectors representing magnets, or directors representing liquid crystal orientations). Their relation with our optical Skyrme density is discussed in the next section.

Polarization hopfions and optical currents
A characteristic feature of any topological texture of closed filaments filling 3D space is that they are mutually linkedinherited by the map between real space and the 3-sphere parametrised by the Hopf fibration which has this property. For polarization textures we can consider this in the context of topological features usually considered in singular optics such as optical vortices and C lines.
Optical vortices are zeros of a complex scalar field representing a polarized scalar component. Mathematically, the points r satisfying (r ) = 0, that is, the preimages of 0. These correspond to the set of vortex lines, which can be infinite lines, loops, or knots [16]. In fact, any complex constant 0 has a preimage defined by (r ) = 0 , and, just like vortex lines, these are typically filaments in 3D.
A similar argument holds for polarization filaments, where now the complex scalar is = E L /E R = L / R as defined in equation (S1.14), whose argument is ↵ and modulus is a monotonic function of and S 3 . However, on LH C lines where E R = 0, may take the value 1. (This notion may be made more mathematically rigorous with the notion of sections of complex bundles, which we avoid to remain close to standard terminology in optics.) As discussed in Supplementary Note 1, this means the parameter space of the map corresponds topologically to a sphere, allowing Hopf textures and, with phase, 3D skyrmions.
The tangent to the filament at each point is given by the vorticity field This is similar to the vorticity field ! sc ⌘ i r ⇤ ⇥ r for a normal complex scalar field [17], but here with an extra denominator regularising the divergence in the neighbourhood of = 1. Such a field is the curl of an optical current associated with the scalar , and evidently ! = r ⇥ J. Again, apart from the regularising denominator, this is the standard current associated with a complex scalar field, J sc = Im ⇤ r . In the theory of knotted fields, the linking of field lines of divergence-free vector fields are often of interest; we write such a field as a "magnetic field" B (although we are interested here in !). The continuous linking density of the field lines of B is given by the helicity density B · A = A · r ⇥ A for "vector potential" A [38,[43][44][45][46]. The integral over all 3D real space gives the total helicity of the field, for a magnetic field related to duality symmetry, and for polarization textures with tangents specified by !, related to the global linking of polarization filaments.
However, it is clear from equation (S4.5) that the corresponding polarization helicity density J · ! = 0 everywhere except where = 1, which has not been suitably regularised. An approach to this was considered in [47] (following from similar calculations in [48]), where the singularity on the LH C line in J is resolved by a singular gauge transformation This new current e J again satisfies r ⇥ e J = !, but now the helicity is nonzero. Using equation (S4.6), the new current is that is, the regularised current e J = J o , the orbital current associated with the field; the numerator is the sum of the two scalar currents associated with the two field components, appropriately normalized. Added to the spin current (which does not play a significant role in the present situation), this gives the Poynting vector of the paraxial optical field [49].
How is this related to the Skyrme density of the previous part? In terms of the angle parameters ↵, , , it is easy to see that Im E ⇤ R rE R = 1 4 (1 + cos )(r r↵) and Im E ⇤ L rE L = 1 4 (1 cos )(r + r↵), implying that that is, the vorticity of the polarization field is the baby skyrmion density vector (see equation (S4.2)). This is elegant but is not surprising, as the vorticity vector arises as the jacobian between real space and field space when counting (discrete) topological charge densities in singular optics [3,50]. As discussed in the main text, the topological role of helicity in flows has a similar expression in terms of an appropriate velocity field [32,45,47,[51][52][53]. We can summarise the preceding discussion in terms of Skyrme density as follows. The directions of the polarization filaments are given by a vector field determined by ! = ⌃ b . The natural vector potential associated with this divergencefree field is singular at the singularities of ↵ (i.e. the C lines), but when regularised by adding the phase gradient r , the polarization helicity density reduces to the 3D Skyrme density, as indeed it must: also equivalent to equation (3) in the main text.
The preceding discussion also illustrates that further gauge transformations-that is, adding nonsingular scalar fields to -does not change the overall Skyrme number. In the study of pure hopfions (i.e. maps from 3D real space to a 2-sphere), whose filament directions are given by B but there is no natural phase, it is usually necessary to construct a 'vector potential' to get the correct helicity density [44,46].
Similarly, if we choose to focus purely on optical polarization hopfions, we can choose an appropriate phase. The associated helicity equation (S4.10) integrated over all 3D space characterises the interlinking of polarization filaments, and globally the number of times the 3D field wraps around the Poincaré sphere. Any convenient choice of phase function (which is appropriately singular on C lines) suffices for this choice. Several choices are discussed in the following section, which applies the preceding analysis to our model field, as introduced in the Methods and Supplementary Note 2.
We also remark on textures for a normal (non-infinity) complex scalar field . Assuming that , as a regular complexvalued function, is nowhere infinite (not even at spatial infinity) means the total set of complex values the field takes is of finite range for a finite set, and hence the texture of its filaments (labelled by 0 ) cannot be topological in spite of some filaments being knotted and linked. If we choose to regarding the vortex lines of (i.e. preimages of zero) being special, we can choose the phase gradient r arg , analogous to the superfluid velocity, as the appropriate current vector. The curl of r arg vanishes everywhere apart from -concentrations on the vortex line; careful regularisation gives a weighting to the screw dislocation nature of the wave vortex [50]. As the global helicity of the field must be zero, the phase around a vortex loop must unwind (i.e. "self-link") whenever it is threaded (i.e. "linked") with another vortex line [54].

Integrating the Skyrme and Hopf densities of the model of skyrmionic hopfion
The topology of the optical skyrmionic hopfion resides in the map from 3D configuration space to the optical hypersphere, parametrised by ↵, , by the fields (S2.9), (S2.11) with the numerical values of the coefficients given in the Methods. As discussed in the main text in the context of the experiment, the polarization of these fields consists of a texture of filaments of constant polarization resembling the Hopf fibration (main text Fig. 3), and phase varying along them establishing a ⇧ 3 mapping between the two 3D spaces. The Skyrme density is given by ⌃ = 1 16⇡ 2 r · (r↵ ⇥ r (cos )). As described in the Methods, the coefficients of this model field have been chosen optimally to concentrate the intensity and Skyrme density. Supplementary Fig. 8a shows a cross-section of the Skyrme density for this field (equivalent to the 3D plot in the inset to Fig. 4c in the main text). The Skyrme density for the model field is strongly concentrated around the C line ring, with most of the contours shown in the plot being several orders of magnitude smaller than the maximum density.
Provided the neighbourhood of the vortex ring is included in the integration volume, we expect the integrated Skyrme density to be close to 1. Supplementary Fig. 8b demonstrates the result of this integral as a function of the transverse (L ? ) and longitudinal (L k ) dimensions of the integration volume given by |x|  L ? , |y |  L ? , |z|  L k . Most of the area shows 99% of the 3-sphere volume is filled. This justifies our claim that the optical skyrmonic hopfion is particle-likethe integration of the topological charge density is independent of the size of the sufficiently large integration box, and is very close to unity.
The evaluation of the measured Skyrme number from the experimental data is described in Methods and Supplementary Note 5. As explained there, rather than construct a smooth interpolation of ⌃ from the data then integrate numerically, we calculate the fraction of the 3-sphere filled by a piecewise linear map constructed directly from the data. Even with high-performance computers this approach requires computational effort. Justified by the particle-like nature of the skyrmionic hopfion, we limit the calculation of the experimental Skyrme number to the cuboid of transverse size L ? = 1.84w = 99.4 µm and longitudinal size L k = 0.768z R = 26.6 mm indicated in Supplementary Fig. 8, a and  b, in cyan. The numerical integral of the model field in this restricted cuboid is 0.9966, i.e. 99.7% of the 3-sphere volume. For reference, the total measured volume (see Methods) is the cuboid given by x max = 3.13w = 170 µm, y max = 3.91w = 212 µm and z max = 0.768z R = 26.6mm. The numerical integral of ⌃ for the model field in the measured volume is 0.997 and is also indicated in Supplementary Fig. 8b in cyan. Thus, to 3 significant figures, the Skyrme number in the reduced volume is no different from the integrated density in the total measured volume.
We note that there is a contour in Supplementary Fig. 8b for L ? > 2w , L k > 1.5z R where the integrated Skyrme volume equals 1. This due to the properties of gaussian beams used to realise the field, which is described, with a mathematical analysis of the asymptotic value of the field for large R and z, in the next subsection. However, for practical purposes it is challenging to measure the field's polarization and phase so far from the focal point, as the intensity is small and subject to experimental imperfections in this region.

Theoretical considerations of model skyrmionic hopfion
We now consider some of the topological details associated with our model field (S2.9), (S2.11) in more detail. For convenience, we use dimensionless cylindrical coordinates R = R/w , z = z /z R , but keeping the coefficients a, b, c. To simplify the expressions, we write 1 + z 2 = W . The Skyrme density ⌃ = 1 16⇡ 2 r · (r↵ ⇥ r (cos )) for the superposition of LG modes given by equations (S2.9) and (S2.11), in terms of these dimensionless coordinates, is . Therefore the integral over infinite radius diverges logarithmically. We will consider the meaning of this in the following.

Polarization parameters
From equation (S1.10), the gradient r↵ for the model field in cylindrical coordinates is which is singular along and circulates around the C lines: along the axis R = 0 and the circle (R, Z) = ( p a/b, 0). From (S1.11), cos = S 3 , which is given by . (S4.13) Its gradient r (cos ) is Together, these determine the polarization behavior of the skyrmionic hopfion, and the field's baby skyrmion/full Poincaré beam nature in any plane. We now consider the Skyrme density and Hopf density (polarization helicity) of the model field for R and L given in equation (S2.9) and equation (S2.11) respectively, found as superpositions of LG modes.

2D baby Skyrme density/full Poincaré beam nature
The model beam's baby Skyrme density (polarization vorticity) vector is, in cylindrical coordinates, The z-component of this vector, ⌃ b z , gives the density covering the Poincaré sphere in a transverse plane. In the focal plane z = 0, the RH axial C line corresponds to the origin, and around it is the LH C line of radius R = p a/b, equal to 0.566 with the numerical coefficients. Since the vortex sign of the LH component is negative, the C point at the origin has a 'star' form [15,16] -↵ decreases as increases in a right-handed sense around the beam axis.
With the numerical values of the coefficients, the integral of ⌃ b z in the focal plane, over the disk to radius p a/b gives 1, equivalent to covering the Poincaré sphere once in a negative sense. This is consistent with the azimuthal sense of r ? ↵ being against b while r ? (cos ) = r ? S 3 points towards the origin within this disk, ensuring b z · r ? ↵ ⇥ r ? (cos ) < 0. Outside the radius of the C line loop, the sign of r ? S 3 changes although r ? ↵ does not, so ⌃ b z > 0; since, asymptotically as R 1, the polarization state is RH, the net baby Skyrme number (degree of the map covering the Poincaré sphere) is zero.
The magnitude of the baby Skyrme density in the focal plane near the circular C line is large; within the circles where S 3 = cos = 0, at R = 0.524, 0.610 with the numerical coefficients, corresponds to covering the entire lower hemisphere, whereas the entire upper hemisphere is covered inside and outside this annulus. By topological continuity, the baby Skyrme number in every other transverse plane is also zero (the radial asymptotic state of polarization being RH circular ensures this is always well-defined). As the LH circular polarization is restricted only to the C line loop in the focal plane, the baby Skyrme number must be zero in every transverse plane except the focal plane.
The 3D Skyrme number involves the product of ⌃ b z with r . Since the C line is a phase vortex in , r changes sign on each side of the circular C line, implying the 3D Skyrme density ⌃ maintains the same sign although ⌃ b z changes sign. We briefly remark on the comparison between the 2D baby skyrmions in our transverse beam profiles and the 'Néel type' (radial) and 'Bloch type' (spiral) baby skyrmions studied in magnetic systems [34]. In a magnetic system, (baby) skyrmions are localised, nonlinear excitations of a ground state (corresponding, e.g. to LH polarization) with Skyrme number 1 in our terminology; the 'type' indicates whether azimuth angles correspond to with an offset or not. In our beams, in the focal plane, = ↵, so there is no spiral-like contribution: the baby skyrmion configuration in the focal plane is purely radial (Néel). However, in other planes, ↵ acquires an offset (which is dependent on radius) due to the phase difference between E R and E L , meaning in a typical plane it appears Bloch-type. Of course, for polarization, the distinction between radial and spiral is relative; since polarization is a pseudospinor, the relation between polarization azimuth ↵ and real space azimimuth is a matter of convention.

Polarization hopfion
We now consider our model field purely as a hopfion, i.e. how the Poincaré sphere is covered by the polarization states in the focal volume. As a map from 3D (real space) to the Poincaré 2-sphere, the preimage of each polarization state is a 1D line -these are the polarization filaments discussed in the main text, generalising C lines for circular polarization.
For an ideal hopfion, the filaments should correspond to a Hopf fibration in real space: loops wound around tori on which S 3 = cos is constant; these tori enclose the LH circular C line loop. The contours of S 3 from equation (S4.13) in a plane of fixed azimuth is shown in Supplementary Fig. 2c. Evidently, S 3 = 1 on the LH C line loop (reference circle) with R = r 0 = p a/b. The ideal hopfion structure would have the contours of constant S 3 loops around this point, topologically equivalent to the contours of |u| in Supplementary Fig.2a. With the numerical choice of a and b as discussed in Supplementary Note 2, S 3 has saddle points at (R, z) = (2 p b, ± p a + b)/ p 3b a = (1.22, ±0.702), at S 3 = 1 8c 2 /((a b) 2 + 4c 2 ) = 0.969.
At values of S 3 larger than this value, 0.969 < S 3  1, the toroidal structure beaks down, into topological cylinders extending to z ! ±1; topologically the hopfion structure breaks down here, although this counts for only 1.5% of the range of S 3 = cos .
As discussed above, the Hopf number of the hopfion can only be calculated with the choice of a suitable phase function e ; using the physical phase gives the skyrmionic hopfion, but the global integral of re · ⌃ b should not change under a gauge transformation e ! e 0 . As noted above, the physical phase fails to converge due to the LG beams' wavefront curvature, so we will first consider the physical phase gradient in more detail.

Physical phase gradient
The phase parameter is the sum of phases of the two beam components, from equation (S1.13). Its gradient therefore is a sum of phase gradients from different phase features of the full field, as superposition of LG beams, where g r is associated with the phase of the C line ring, g a with the axial C line, g G with the Gouy phase, and g c with the wavefront curvature. The separation into these terms is very natural for our model, consisting of a superposition of three LG beams, with (`, p) = (0, 0), (0, 1) and ( 1, 0). All of these have the same wavefront curvature factor, and the axial phase is due to the negative vortex in the LH component. Since the beam with (`, p) = (0, 1) dominates the superposition for large R and z, we choose the Gouy phase of this beam to represent g G ; the phase for the ring g r takes all the other phase contributions, including the contributions from the Laguerre polynomials. Only g r involves the coefficients of the superposition.
The full Skyrme density ⌃ = r · ⌃ b can therefore be considered as the contribution of separate terms, ⌃ • = g • · ⌃ b , for • = r, a, G, c respectively. The important contributions for the skyrmion topology are those phases corresponding to phase singularities in the physical field (and singularities in ↵), that is, from the C lines. These are therefore g r and g a .
The various phase gradients are readily found to be g r = 2(a b) ((a bR 2 ) 2 +2(a 2 ab(2 R 2 )+2b 2 (1 R 2 ))z 2 +(a 2b) 2 z 4 ) Evidently, from equation (S4.12), r↵ = g a g r + (1 + z 2 ) 1 b z; the final term accounts for the difference in Gouy phases of the RH and LH beams, and otherwise r↵ is the difference of the axial and ring phase gradients whilst r includes their sum, just like how, as described in Supplementary Note 1, arg U and arg V are the sum and difference of ↵ and for the Hopf fibration of the optical hypersphere.
Both ⌃ r and ⌃ a must therefore be included in any function determining the Hopf density of the full field. Numerical integration of the cut volumes of |x|, |y |  1.84, |z|  0.768 gives 0.4964 and 0.4967 for the axial and ring contributions to the Skyrme number, adding to 0.9931.
The other phase gradients g G and g c arise from the physics of propagating gaussian beams. The Gouy phase g G is purely in the z-direction, so the corresponding Hopf density ⌃ G is proportional to the transverse baby Skyrme density (as indeed would be any phase associated with the e i kz propagation); as discussed above, its integral over any transverse plane is zero, meaning the integral of ⌃ G over all space is zero. Within the investigated volume, R vol ⌃ G d 3 r = 0.0014, making an almost negligible contribution.
Unlike the other contributions to r , the wavefront curvature phase gradient g c diverges as R 1, since the phase curvature does not settle to a constant value for R 1 when z 6 = 0. In this regime, ⌃ c ⇠ c 2 /⇡ 2 b 2 R 2 , accounting for the total logarithmic divergence in the total ⌃. This contribution accounts for why the integrated Skyrme density exceeds 1 for much larger volumes, as can be seen in Supplementary Fig. 8b when the boundaries L ? and L k are quite large.
Although mathematically a divergence, the actual divergence is very slow, becoming large where the beam intensity is low, and thus hard to measure physically. In fact, the paraxial approximation of the physical beam in terms of the mathematical model breaks down for sufficiently large R and z. Within the cut volume as described above, the wavefront curvature contributes 0.0020 to the Skyrme number; slightly larger than the Gouy phase, but still very small.
The sum of these four contributions gives the value of the Skyrme number numerically integrated over the cut volume, as described above: 0.4964+0.4967+ 0.0015+0.0020 = 0.9966.

Global Hopf number
From the previous section, it is natural to calculate the global Hopf number as the integral over all space of ⌃ a + ⌃ r ; these include the important topological helicity terms and do not diverge at infinity. With the numerical values of the coefficients, this global Hopf number in fact gives zero; the nontrivial degree of the localised hopfion within the measured volume is, mathematically, unwound within the larger volumes, beyond the saddle of cos in R, z.
This does not remove the "particle-like" nature of the measured hopfion, although the precise numerical value depends on precisely which phase contributions are included; unlike the Skyrme density, which is physically defined in any volume, the more general Hopf density is only unambiguous within the infinite volume. where m, n are nonnegative integers, the first equalities define the fields in terms of cylindrical coordinates R, , z, and the second equalities in terms of toroidal coordinates [11][12][13], similar to Supplementary Note 1 but here applied to real 3D space. We perform the calculation here in toroidal coordinates, which is very natural for the ideal skyrmion field. This can be compared to the similar calculation in [55], which assumes spherical symmetry. The parameters of this definition of the field, from equations (S1.10), (S1.11) and (S1.13), are easily expressed in toroidal coordinates, = m + n , ↵ = m n , cos = sech 2n ⌧ tanh 2n ⌧ sech 2n ⌧ + tanh 2n ⌧ .

=
; (S4.20) The volume element for toroidal coordinates is sinh ⌧ d d d⌧/(cosh ⌧ cos ) 3 , so the Skyrme density ⌃ = r · r↵ ⇥ r (cos ), and dividing by the 3-sphere volume 16⇡ 2 , gives the total integrated Skyrme number where, since the net integrand is dependent only on ⌧ , the angle integrals are trivial. The net integral is identified as a pure derivative, so the overall ⌧ integral is also straightforward, leaving a net result mn from the phase windings in equation (S4.18). Choosing the conjugate of either of the fields in the definition (S4.18) gives a positive Skyrme density, as in previous discussion where m = n = 1 and a negative axial vortex was chosen. This is instructive for designing optical skyrmionic hopfions of higher Skyrme number, e.g. using superpositions of LG beams. The strategy is similar to that employed by [56], who designed scalar fields with vortex loops threaded by vortices on the axis. Increasing the axial degree m is easy, by replacing the LH component with a different azimuthal quantum number`= m. However, as described in [56], paraxial beams do not allow transverse vortex loops with |n| > 1, Realising a higher-order skyrmion with effective |n| > 1 is possible with multiple vortex loops with the same sense of vorticity, and indeed, these can be woven into various kinds of link and knot as constructed in [47,48].

Supplementary Note 5: Details of numerical methods
Given the the map from real space to the optical hypersphere described in the numerical calculation of the experimental Skyrme Number (see Methods), we can visualise the optical skyrmionic hopfion and its Skyrme density. We imagine each measured point in real space to correspond to the vertices of a cubic lattice, which is carried to an irregular lattice with the same topology in the optical hypersphere. In regions of the beam where the Skyrme density is high, such as close to the circular C line, the corresponding spherical tetrahedra (and cubes comprised of them) are large, and where the Skyrme density is small so are the spherical tetrahedra. Explicative examples are the grey cubes in Fig. 4 of the main text: in real space (Fig. 4a) the cubes have the same size, whereas in the optical hypersphere, (Fig. 4b), the cube enclosing the C line loop (in black) is much larger than the other cube. Supplementary Fig. 9, a and b, show the analogous mapping for the ideal skyrmion field. Main text Fig. 4c shows all of the sampled cubes in real space, whose colour represents Vol(C): the transparency goes to 0 for very small volumes, becomes increasingly orange for positive volumes and increasingly blue for negative volumes. Our algorithm is based on expressing the measured data as the vertices of tetrahedral meshes (i.e. 3D cell complexes) in real space and the optical hypersphere. The calculations of the spherical tetrahedra volumes in the 3-sphere are described, with a final section describing how this is found from the image of the boundary of the measured volume.

The volume of a spherical tetrahedron
The 3-sphere is embedded in four-dimensional Euclidean space. Four unit vectors,ñ`with`= {a, b, c, d}, point to the four vertices of a spherical tetrahedron T . The four faces of this spherical tetrahedron are spherical triangles, pairs of which meet along edges at the dihedral angles 0 < ' j < ⇡ for j = 1, 2, . . . , 6. The formula for Vol(T ) [57], 3D volume of the tetrahedron in the hypersphere, is in terms of the dihedral angles.
Given a j = exp(i' j ) and the dilogarithm function Li 2 (z), the volume of the spherical tetrahedron is, from Murakami's formula [57], the volume-preserving projection with the optical state at the focal point, ( 1, 0). The antipodal state ( 1, 0) therefore corresponds to the sphere boundary of the solid ball. The image of the measured boundary, being a neighbourhood of this point, therefore would appear to be a sphere close to the boundary sphere in this projection. The boundary of the measured volume, that is, the image of the six faces of the measured cubic volume, are plotted in the described projection of the optical hypersphere in Supplementary Fig. 10a. Clearly this is a surface that almost coincides with the boundary of the projection at (1, 0). The corresponding plot for the theoretical model field is shown in the inset.
In this projection it is hard to see that the volume of this neighbourhood complementing the filled volume is small. For this reason, in Supplementary Fig. 10b, we plot the optical hypersphere in a different projection, with a point on the circular C line at the center (and the opposite point on the C line loop to the boundary of the ball). The RH states on the axis in real space become a circle in this projection, and it is clear here that the volume of the optical hypersphere not covered is small (consistent with the theoretical model, inset to Supplementary Fig. 10b).