#  PACMAN: A pipeline to reduce and analyze Hubble Wide Field Camera 3 IR
Grism data

Sebastian Zieba  0000-0003-0562-6750 Laura Kreidberg  0000-0003-0514-1147
0000-0001-7699-3983

††margin:  DOI:
https://doi.org/10.21105/joss.04838https://doi.org10.21105/joss.04838 Software
• https://github.com/openjournals/joss-
reviews/issues/4838https://doi.orgReview • https://github.com/sebastian-
zieba/PACMANhttps://doi.orgRepository •
https://doi.org/10.5281/zenodo.7465579https://doi.orgArchive Editor:
https://juanjobazan.com/https://doi.orgJuanjo Bazán
Reviewers: • https://github.com/astrobelhttps://doi.org@astrobel •
https://github.com/aureliocarnerohttps://doi.org@aureliocarnero Submitted: 25
July 2022
Published: 21 December 2022 License
Authors of papers retain copyright and release the work under a Creative
Commons Attribution 4.0 International License
(http://creativecommons.org/licenses/by/4.0/https://doi.orgCC BY 4.0).

## Summary

The Hubble Space Telescope (HST) has become the preeminent workhorse facility
for the characterization of extrasolar planets. Launched in 1990 and never
designed for the observations of exoplanets, the STIS spectrograph on HST was
used in 2002 to detect the first atmosphere ever discovered on a planet
outside of our solar system (Charbonneau et al., 2002).

HST currently has two of the most powerful space-based tools for
characterizing exoplanets over a broad spectral range: The Space Telescope
Imaging Spectrograph (STIS; installed in 1997) in the UV and the Wide Field
Camera 3 (WFC3; installed in 2009) in the Near Infrared (NIR). With the
introduction of a spatial scan mode on WFC3 (Deming et al., 2012; McCullough &
MacKenty, 2012) where the star moves perpendicular to the dispersion direction
during an exposure, WFC3 observations have become very efficient due to the
reduction of overhead time and the possibility of longer exposures without
saturation.

For exoplanet characterization, WFC3 is used for transit and secondary eclipse
spectroscopy, and phase curve observations. The instrument has two different
grisms: G102 with a spectral range from 800 nm to up to 1150 nm and G141
encompassing 1075 nm to about 1700 nm. The spectral range of WFC3/G141 is
primarily sensitive to molecular absorption from water at approximately 1.4
microns. This led to the successful detection of water in the atmosphere of
over a dozen of exoplanets (e.g., Deming et al., 2013; Evans et al., 2016;
Fraine et al., 2014; Huitson et al., 2013; Kreidberg, Bean, Désert, Line, et
al., 2014). The bluer part of WFC3, the G102 grism, is also sensitive to water
and most notably led to the first detection of a helium exosphere (Spake et
al., 2018).

Here we present PACMAN, an end-to-end pipeline developed to reduce and analyze
HST/WFC3 data. The pipeline includes both spectral extraction and light curve
fitting. The foundation of PACMAN has been already used in numerous
publications (e.g., Kreidberg, Bean, Désert, Line, et al., 2014; Kreidberg et
al., 2018) and these papers have already accumulated hundreds of citations.

## Statement of need

Exoplanet spectroscopy with Hubble requires very precise measurements that are
beyond the scope of standard analysis tools provided by the Space Telescope
Science Institute. The data analysis is challenging, and different pipelines
have produced discrepant results in the literature (e.g., Kreidberg et al.,
2019; Teachey & Kipping, 2018). To facilitate reproducibility and
transparency, the data reduction and analysis software should be open-source.
This will enable easy comparison between different pipelines, and also lower
the barrier to entry for newcomers in the exoplanet atmosphere field.

What sets PACMAN apart from other tools provided by the community, is that it
was specifically designed to reduce and fit HST data. There are several open-
source tools that can fit time series observations of stars to model events
like transiting exoplanets, such as EXOFASTv2 (Eastman et al., 2019), juliet
(Espinoza et al., 2019), allesfitter (Günther & Daylan, 2019 , 2021),
exoplanet (Foreman-Mackey et al., 2021a, 2021b), and starry (Luger et al.,
2019). PACMAN’s source code, however, includes fitting models that can model
systematics which are characteristic to HST data, such as the orbit-long
exponential ramps due to charge trapping or the upstream-downstream effect.
This removes the need for the user to write these functions themselves. PACMAN
will also retrieve information from the header of the FITS files,
automatically detect HST orbits and visits and use this information in the
fitting models.

The only other end-to-end open source pipeline specifically developed for the
reduction and analysis of HST/WFC3 data is https://github.com/ucl-
exoplanets/Iraclishttps://doi.orgIraclis (Tsiaras et al., 2016). Another open-
source pipeline that has been for example used as an independent check of
recent results presented in Mugnai et al. (2021) and Carone et al. (2021) is
https://jbouwman.gitlab.io/CASCADe/https://doi.orgCASCADe (Calibration of
trAnsit Spectroscopy using CAusal Data). For a more detailed discussion of
CASCADe see Appendix 1 in Carone et al. (2021).

## Outline of the pipeline steps

The pipeline starts with the _ima_ data products provided by the Space
Telescope Science Institute that can be easily accessed from
https://mast.stsci.edu/search/hsthttps://doi.orgMAST. These files created by
the WFC3 calibration pipeline, calwf3, have already several calibrations
applied (dark subtraction, linearity correction, flat-fielding) to each
readout of the IR exposure.

In the following we highlight several steps in the reduction and fitting
stages of the code which are typical for HST/WFC3 observations:

  * •

Wavelength calibration: We create a reference spectrum based on the throughput
of the respective grism (G102 or G141) and a stellar model. The user can
decide if he or she wants to download a stellar spectrum from MAST or use a
black body spectrum. This template is used for the wavelength calibration of
the WFC3 spectra. We also determine the position of the star in the direct
images which are commonly taken at the start of HST orbits to create an
initial guess for the wavelength solution using the known dispersion of the
grism. Using the reference spectrum as a template, we determine a shift and
scaling in wavelength-space that minimizes the difference between the template
and the first spectrum in the visit. This first exposure in the visit is then
used as the template for the following exposures in the visit.

  * •

Optimal extraction and outlier removal: PACMAN uses an optimal extraction
algorithm as presented in Horne (1986) which iteratively masks bad pixels in
the image. We also mask bad pixels that have been flagged by calwf3 with data
quality DQ = 4 or 512111for a list of DQ flags see
https://wfc3tools.readthedocs.io/en/latest/wfc3tools/calwf3.html.

  * •

Scanning of the detector: The majority of exoplanetary HST/WFC3 observations
use the spatial scanning technique (McCullough & MacKenty, 2012) which spreads
the light perpendicular to the dispersion direction during the exposure
enabling longer integration times before saturation. The _ima_ files taken in
this observation mode consist of a number of nondestructive reads, also known
as up-the-ramp samples, each of which we treat as an independent subexposure.
Figure 1 (left panel) shows an example of the last subexposure when using
spatial scanning together with the expected position of the trace based on the
direct image.

  * •

Fitting models: PACMAN contains several functions to fit models which are
commonly used with HST data. The user can fit models like in Equation 1 to the
white light curve or to spectroscopic light curves. An example of a raw
spectroscopic light curve and fitting Equation 1 to it, can be found in Figure
2. Here are some examples of the currently implemented models for the
instrument systematics and the astrophysical signal:

    * –

systematics models:

      * *

visit-long polynomials

      * *

orbit-long exponential ramps due to charge trapping: NIR detectors like
HST/WFC3 can trap photoelectrons (Smith et al., 2008), which will cause the
number of recorded photoelectrons to increase exponentially, creating typical
hook-like features in each orbit

    * –

astrophysical models:

      * *

transit and secondary eclipse curves as implemented in batman

      * *

sinusoids for phase curve fits

      * *

a constant offset that accounts for the upstream-downstream effect (McCullough
& MacKenty, 2012) caused by forward and reverse scanning

A typical model to fit an exoplanet transit in HST data is the following
(used, for example, by Kreidberg, Bean, Désert, Line, et al., 2014):

$F(t)=T(t)\,(c\,S(t)+k\,t_{\rm{v}})\,(1-\exp(-r_{1}\,t_{\rm{orb}}-r_{2})),$
(1)

with _T(t)_ being the transit model, _c_ (_k_) a constant (slope), _S(t)_ a
scale factor equal to 1 for exposures with spatial scanning in the forward
direction, and _s_ for reverse scans, $r_{\rm{1}}$ and $r_{\rm{2}}$ are
parameters to account for the exponential ramps. $t_{\rm{v}}$ and
$t_{\rm{orb}}$ are the times from the first exposure in the visit and in the
orbit, respectively.

  * •

Parameter estimation: The user has different options to estimate best fitting
parameters and their uncertainties:

    * –

least squared: scipy.optimize

    * –

MCMC: emcee (Foreman-Mackey et al., 2013)

    * –

nested sampling: dynesty (Speagle, 2020)

  * •

Multi-visit observations

    * –

PACMAN has also an option to share parameters across visits.

  * •

Binning of the light spectrum: The user can freely specify the bin numbers or
locations. Figure 1 (right panel) shows the resulting 1D spectrum and a user-
defined binning.

Figure 1 and Figure 2 show some figures created by PACMAN during a run using
three HST visits of GJ 1214 b collected in
https://archive.stsci.edu/proposal_search.php?id=13021&mission=hsthttps://doi.orgGO
13201 (Bean, 2012). An analysis of all 15 visits was published in Kreidberg,
Bean, Désert, Line, et al., (2014). The analysis of three visits here using
PACMAN, is consistent with the published results.

Figure 1: _Left panel_ : a typical single exposure showing the raw 2D
spectrum. _Right panel_ : 1D spectrum after the use of optimal extraction
including vertical dashed lines showing the user-set binning to generate
spectroscopic light curves. Figure 2: _panel A_ : raw white light curves for
each of the three visits. One can clearly see the constant offset between two
adjacent exposures due to the spatial scanning mode. _panel B_ : white light
curve with the best astrophysical model fit using Equation 1. _panel C_ : the
transmission spectrum after fitting 11 spectroscopic light curves revealing
the flat spectrum of GJ 1214 b as published in Kreidberg, Bean, Désert, Line,
et al., (2014).

## Dependencies

PACMAN uses typical dependencies of astrophysical Python codes: numpy (Harris
et al., 2020), matplotlib (Hunter, 2007), scipy (Virtanen et al., 2020) and
astropy (Astropy Collaboration et al., 2022, 2018, 2013).

Other dependencies required for the fitting stage depending on the model and
sampler being run are: batman (Kreidberg, 2015), emcee (Foreman-Mackey et al.,
2013), dynesty (Speagle, 2020), and corner (Foreman-Mackey, 2016).

For the barycentric correction, PACMAN accesses the https://ssd-
api.jpl.nasa.gov/obsolete/horizons_batch_cgi.htmlhttps://doi.orgAPI to JPL’s
Horizons system.

If the user decides to use a stellar spectrum for the wavelength calibration,
PACMAN will download the needed fits file from the
https://archive.stsci.edu/hlsps/reference-
atlases/cdbs/grid/https://doi.org“REFERENCE-ATLASES’ ’ high level science
product hosted on the MAST archive (STScI Development Team, 2013).

## Documentation

The documentation for PACMAN can be found at
https://pacmandocs.readthedocs.io/en/latest/https://doi.orgpacmandocs.readthedocs.io
hosted on https://readthedocs.org/https://doi.orgReadTheDocs. It includes most
notably, a full explanation of every parameter in the _pacman control file_
(pcf), the API, and an example of how to download, reduce and analyze
observations of GJ 1214 b taken with HST/WFC3/G141.

## Future work

The following features are planned for future development:

  * •

The addition of fitting models like phase curves using the open-source Python
package SPIDERMAN (Louden & Kreidberg, 2018).

  * •

Orbit-long ramp fitting using the
https://recte.readthedocs.io/en/latest/https://doi.orgRECTE systematic model.

  * •

Limb darkening calculations for users wanting to fix limb darkening parameters
to theoretical models in the fitting stage.

  * •

Extension to WFC3/UVIS data reduction.

## Acknowledgements

We acknowledge B. Zawadzki for the creation of the PACMAN logo. We also
acknowledge the comments and contributions by I. Momcheva to PACMAN.

## References

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., Earl, N., Starkman,
N., Bradley, L., Shupe, D. L., Patil, A. A., Corrales, L., Brasseur, C. E.,
Nöthe, M., Donath, A., Tollerud, E., Morris, B. M., Ginsburg, A., Vaher, E.,
Weaver, B. A., Tocknell, J., Jamieson, W., … Astropy Project Contributors.
(2022). The Astropy Project: Sustaining and Growing a Community-oriented Open-
source Project and the Latest Major Release (v5.0) of the Core Package.
_Astrophysical Journal_ , _935_(2), 167.
https://doi.org/10.3847/1538-4357/ac7c74

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M.,
Lim, P. L., Crawford, S. M., Conseil, S., Shupe, D. L., Craig, M. W.,
Dencheva, N., Ginsburg, A., VanderPlas, J. T., Bradley, L. D., Pérez-Suárez,
D., de Val-Borro, M., Aldcroft, T. L., Cruz, K. L., Robitaille, T. P.,
Tollerud, E. J., … Astropy Contributors. (2018). The Astropy Project: Building
an Open-science Project and Status of the v2.0 Core Package. _Astronomical
Journal_ , _156_(3), 123. https://doi.org/10.3847/1538-3881/aabc4f

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., Greenfield, P.,
Droettboom, M., Bray, E., Aldcroft, T., Davis, M., Ginsburg, A., Price-Whelan,
A. M., Kerzendorf, W. E., Conley, A., Crighton, N., Barbary, K., Muna, D.,
Ferguson, H., Grollier, F., Parikh, M. M., Nair, P. H., … Streicher, O.
(2013). Astropy: A community Python package for astronomy. _Astronomy and
Astrophysics_ , _558_ , A33. https://doi.org/10.1051/0004-6361/201322068

Bean, J. (2012). _Revealing the Diversity of Super-Earth Atmospheres_ (p.
13021). HST Proposal ID 13021. Cycle 20.

Carone, L., Mollière, P., Zhou, Y., Bouwman, J., Yan, F., Baeyens, R., Apai,
D., Espinoza, N., Rackham, B. V., Jordán, A., Angerhausen, D., Decin, L.,
Lendl, M., Venot, O., & Henning, T. (2021). Indications for very high
metallicity and absence of methane in the eccentric exo-Saturn WASP-117b.
_Astronomy and Astrophysics_ , _646_ , A168.
https://doi.org/10.1051/0004-6361/202038620

Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L. (2002).
Detection of an Extrasolar Planet Atmosphere. _Astrophysical Journal_ ,
_568_(1), 377–384. https://doi.org/10.1086/338770

Deming, D., Wilkins, A., McCullough, P., Burrows, A., Fortney, J. J., Agol,
E., Dobbs-Dixon, I., Madhusudhan, N., Crouzet, N., Desert, J.-M., Gilliland,
R. L., Haynes, K., Knutson, H. A., Line, M., Magic, Z., Mandell, A. M.,
Ranjan, S., Charbonneau, D., Clampin, M., … Showman, A. P. (2013). Infrared
Transmission Spectroscopy of the Exoplanets HD 209458b and XO-1b Using the
Wide Field Camera-3 on the Hubble Space Telescope. _Astrophysical Journal_ ,
_774_(2), 95\. https://doi.org/10.1088/0004-637X/774/2/95

Deming, D., Wilkins, A., McCullough, P., Madhusudhan, N., Agol, E., Burrows,
A., Charbonneau, D., Clampin, M., Desert, J., Gilliland, R., Knutson, H.,
Mandell, A., Ranjan, S., Seager, S., & Showman, A. (2012). Infrared
Spectroscopy of the Transiting Exoplanets HD189733b and XO-1 Using Hubble WFC3
in Spatial Scan Mode. _American Astronomical Society Meeting Abstracts #219_ ,
_219_ , 405.05.

Eastman, J. D., Rodriguez, J. E., Agol, E., Stassun, K. G., Beatty, T. G.,
Vanderburg, A., Gaudi, B. S., Collins, K. A., & Luger, R. (2019). EXOFASTv2: A
public, generalized, publication-quality exoplanet modeling code. _arXiv
e-Prints_ , arXiv:1907.09480. http://arxiv.org/abs/1907.09480

Espinoza, N., Kossakowski, D., & Brahm, R. (2019). juliet: a versatile
modelling tool for transiting and non-transiting exoplanetary systems.
_Monthly Notices of the RAS_ , _490_(2), 2262–2283.
https://doi.org/10.1093/mnras/stz2688

Evans, T. M., Sing, D. K., Wakeford, H. R., Nikolov, N., Ballester, G. E.,
Drummond, B., Kataria, T., Gibson, N. P., Amundsen, D. S., & Spake, J. (2016).
Detection of H2O and Evidence for TiO/VO in an Ultra-hot Exoplanet Atmosphere.
_Astrophysical Journal, Letters_ , _822_(1), L4.
https://doi.org/10.3847/2041-8205/822/1/L4

Foreman-Mackey, D. (2016). Corner.py: Scatterplot matrices in python. _The
Journal of Open Source Software_ , _1_(2), 24.
https://doi.org/10.21105/joss.00024

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. (2013). emcee: The
MCMC Hammer. _Publications of the ASP_ , _125_(925), 306.
https://doi.org/10.1086/670067

Foreman-Mackey, D., Luger, R., Agol, E., Barclay, T., Bouma, L. G., Brandt, T.
D., Czekala, I., David, T. J., Dong, J., Gilbert, E. A., Gordon, T. A.,
Hedges, C., Hey, D. R., Morris, B. M., Price-Whelan, A. M., & Savel, A. B.
(2021a). _exoplanet: Gradient-based probabilistic inference for exoplanet data
& other astronomical time series_ (Version 0.5.1). Zenodo; Zenodo.
https://doi.org/10.5281/zenodo.1998447

Foreman-Mackey, D., Luger, R., Agol, E., Barclay, T., Bouma, L., Brandt, T.,
Czekala, I., David, T., Dong, J., Gilbert, E., Gordon, T., Hedges, C., Hey,
D., Morris, B., Price-Whelan, A., & Savel, A. (2021b). exoplanet: Gradient-
based probabilistic inference for exoplanet data & other astronomical time
series. _The Journal of Open Source Software_ , _6_(62), 3285.
https://doi.org/10.21105/joss.03285

Fraine, J., Deming, D., Benneke, B., Knutson, H., Jordán, A., Espinoza, N.,
Madhusudhan, N., Wilkins, A., & Todorov, K. (2014). Water vapour absorption in
the clear atmosphere of a Neptune-sized exoplanet. _Nature_ , _513_(7519),
526–529. https://doi.org/10.1038/nature13785

Günther, M. N., & Daylan, T. (2019). _allesfitter: Flexible star and exoplanet
inference from photometry and radial velocity_ (p. ascl:1903.003).
Astrophysics Source Code Library, record ascl:1903.003.

Günther, M. N., & Daylan, T. (2021). Allesfitter: Flexible Star and Exoplanet
Inference from Photometry and Radial Velocity. _Astrophysical Journal,
Supplement_ , _254_(1), 13. https://doi.org/10.3847/1538-4365/abe70e

Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P.,
Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R.,
Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J.
F., Wiebe, M., Peterson, P., … Oliphant, T. E. (2020). Array programming with
NumPy. _Nature_ , _585_(7825), 357–362.
https://doi.org/10.1038/s41586-020-2649-2

Horne, K. (1986). An optimal extraction algorithm for CCD spectroscopy.
_Publications of the ASP_ , _98_ , 609–617. https://doi.org/10.1086/131801

Huitson, C. M., Sing, D. K., Pont, F., Fortney, J. J., Burrows, A. S., Wilson,
P. A., Ballester, G. E., Nikolov, N., Gibson, N. P., Deming, D., Aigrain, S.,
Evans, T. M., Henry, G. W., Lecavelier des Etangs, A., Showman, A. P., Vidal-
Madjar, A., & Zahnle, K. (2013). An HST optical-to-near-IR transmission
spectrum of the hot Jupiter WASP-19b: detection of atmospheric water and
likely absence of TiO. _Monthly Notices of the RAS_ , _434_(4), 3252–3274.
https://doi.org/10.1093/mnras/stt1243

Hunter, J. D. (2007). Matplotlib: A 2D Graphics Environment. _Computing in
Science and Engineering_ , _9_(3), 90–95. https://doi.org/10.1109/MCSE.2007.55

Kreidberg, L. (2015). batman: BAsic Transit Model cAlculatioN in Python.
_Publications of the ASP_ , _127_(957), 1161. https://doi.org/10.1086/683602

Kreidberg, L., Bean, J. L., Désert, J.-M., Benneke, B., Deming, D., Stevenson,
K. B., Seager, S., Berta-Thompson, Z., Seifahrt, A., & Homeier, D. (2014).
Clouds in the atmosphere of the super-Earth exoplanet GJ1214b. _Nature_ ,
_505_(7481), 69–72. https://doi.org/10.1038/nature12888

Kreidberg, L., Bean, J. L., Désert, J.-M., Line, M. R., Fortney, J. J.,
Madhusudhan, N., Stevenson, K. B., Showman, A. P., Charbonneau, D.,
McCullough, P. R., Seager, S., Burrows, A., Henry, G. W., Williamson, M.,
Kataria, T., & Homeier, D. (2014). A Precise Water Abundance Measurement for
the Hot Jupiter WASP-43b. _Astrophysical Journal, Letters_ , _793_(2), L27.
https://doi.org/10.1088/2041-8205/793/2/L27

Kreidberg, L., Line, M. R., Parmentier, V., Stevenson, K. B., Louden, T.,
Bonnefoy, M., Faherty, J. K., Henry, G. W., Williamson, M. H., Stassun, K.,
Beatty, T. G., Bean, J. L., Fortney, J. J., Showman, A. P., Désert, J.-M., &
Arcangeli, J. (2018). Global Climate and Atmospheric Composition of the Ultra-
hot Jupiter WASP-103b from HST and Spitzer Phase Curve Observations.
_Astronomical Journal_ , _156_(1), 17\.
https://doi.org/10.3847/1538-3881/aac3df

Kreidberg, L., Luger, R., & Bedell, M. (2019). No Evidence for Lunar Transit
in New Analysis of Hubble Space Telescope Observations of the Kepler-1625
System. _Astrophysical Journal, Letters_ , _877_(2), L15.
https://doi.org/10.3847/2041-8213/ab20c8

Louden, T., & Kreidberg, L. (2018). SPIDERMAN: an open-source code to model
phase curves and secondary eclipses. _Monthly Notices of the RAS_ , _477_(2),
2613–2627. https://doi.org/10.1093/mnras/sty558

Luger, R., Agol, E., Foreman-Mackey, D., Fleming, D. P., Lustig-Yaeger, J., &
Deitrick, R. (2019). starry: Analytic Occultation Light Curves. _Astronomical
Journal_ , _157_(2), 64. https://doi.org/10.3847/1538-3881/aae8e5

McCullough, P., & MacKenty, J. (2012). _Considerations for using Spatial Scans
with WFC3_ (p. 8). Space Telescope WFC Instrument Science Report.

Mugnai, L. V., Modirrousta-Galian, D., Edwards, B., Changeat, Q., Bouwman, J.,
Morello, G., Al-Refaie, A., Baeyens, R., Bieger, M. F., Blain, D., Gressier,
A., Guilluy, G., Jaziri, Y., Kiefer, F., Morvan, M., Pluriel, W., Poveda, M.,
Skaf, N., Whiteford, N., … Beaulieu, J.-P. (2021). ARES. V. No Evidence For
Molecular Absorption in the HST WFC3 Spectrum of GJ 1132 b. _Astronomical
Journal_ , _161_(6), 284\. https://doi.org/10.3847/1538-3881/abf3c3

Smith, R. M., Zavodny, M., Rahmer, G., & Bonati, M. (2008). A theory for image
persistence in HgCdTe photodiodes. In D. A. Dorn & A. D. Holland (Eds.), _High
energy, optical, and infrared detectors for astronomy III_ (Vol. 7021, p.
70210J). https://doi.org/10.1117/12.789372

Spake, J. J., Sing, D. K., Evans, T. M., Oklopčić, A., Bourrier, V.,
Kreidberg, L., Rackham, B. V., Irwin, J., Ehrenreich, D., Wyttenbach, A.,
Wakeford, H. R., Zhou, Y., Chubb, K. L., Nikolov, N., Goyal, J. M., Henry, G.
W., Williamson, M. H., Blumenthal, S., … Madhusudhan, N. (2018). Helium in the
eroding atmosphere of an exoplanet. _Nature_ , _557_(7703), 68–70.
https://doi.org/10.1038/s41586-018-0067-5

Speagle, J. S. (2020). DYNESTY: a dynamic nested sampling package for
estimating Bayesian posteriors and evidences. _Monthly Notices of the RAS_ ,
_493_(3), 3132–3158. https://doi.org/10.1093/mnras/staa278

STScI Development Team. (2013). _pysynphot: Synthetic photometry software
package_ (p. ascl:1303.023).

Teachey, A., & Kipping, D. M. (2018). Evidence for a large exomoon orbiting
Kepler-1625b. _Science Advances_ , _4_(10), eaav1784.
https://doi.org/10.1126/sciadv.aav1784

Tsiaras, A., Waldmann, I. P., Rocchetto, M., Varley, R., Morello, G., Damiano,
M., & Tinetti, G. (2016). A New Approach to Analyzing HST Spatial Scans: The
Transmission Spectrum of HD 209458 b. _Astrophysical Journal_ , _832_(2), 202.
https://doi.org/10.3847/0004-637X/832/2/202

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T.,
Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der
Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R.
J., Jones, E., Kern, R., Larson, E., … SciPy 1. 0 Contributors. (2020). SciPy
1.0: fundamental algorithms for scientific computing in Python. _Nature
Methods_ , _17_ , 261–272. https://doi.org/10.1038/s41592-019-0686-2


Searching for the sources of galactic cosmic rays 2015

11institutetext: Institute of Particle and Nuclear Physics, Faculty of
Mathematics and Physics, Charles University
V Holesovickach 2, 180 00 Prague 8, Czech Republic

# Variability of VHE $\gamma$–ray emission from the binary PSR B1259–63/LS
2883

Stanislav Stefanik 11<EMAIL_ADDRESS>Dalibor Nosek 11

###### Abstract

We examine changes of the $\gamma$-ray intensity observed from the direction
of the binary system PSR B1259–63/LS 2883 during campaigns around its three
periastron passages. A simple and straightforward method is applied to the
published data obtained with the Imaging Atmospheric Cherenkov Technique.
Regardless of many issues of the detection process, the method works only with
numbers of very high energetic photons registered in the specified regions.
Within the realm of this scheme, we recognized changes attributable to the
variations of the intrinsic source activity at high levels of significance.

##  1 Introduction

PSR B1259–63/LS 2883 is one of few $\gamma$–ray binaries detected in the very
high energetic (VHE) regime. It is known to consist out of a Be star with an
equatorial disk tilted relative to the movement of the pulsar which orbits
around the star with a period $T=1237\leavevmode\nobreak\ \mathrm{days}$
Johnston . The variability of the source emission has been clearly observed in
multiple wavebands from radio to VHE with the most prominent changes occurring
around the periastron passages hess_2013 . The cause of this variations is
believed to be the interaction of the pulsar with the stellar disk.
Significant variability of the VHE $\gamma$–ray flux around periastron
hess_2005 ; hess_2009  was observed by the H.E.S.S. telescope array. Although
a considerable overlap of the pre– and post–periastron data has not been
obtained so far, a stronger evidence for the repetitive behaviour stemming
from the orbital motion was claimed based on the correspondence of the newest
data with the older lighcurve hess_2013 . Recently, detection of a GeV flare
by the _Fermi_ –LAT Abdo ; Tam  and a lack of evidence for an attributable VHE
flare hess_2013  lead to the conclusion that the high energy (HE) photons have
a different origin than the VHE ones. In this paper, we reexamine the changes
of VHE activity of PSR B1259–63/LS 2883 during campaigns targeted at periods
of periastron passages.

We deal with publicly available VHE data gathered by Cherenkov telescopes in
the on/off observation mode Berge . In this approach, the numbers of detected
events $N_{\mathrm{on}}$ and $N_{\mathrm{off}}$ are extracted from the
observed signal (on) and background (off) regions, respectively. They are
related by the on–off parameter $\alpha$ which is the ratio of exposures of
signal and background regions. Traditionally, the level of significance of the
source presence in the on–source region is determined using the Li–Ma method
LiMa . In this work, we use its modification for the investigation of changes
in the source intensity (see Nosek ; Stefanik  for details). A source
parameter $\beta>0$ is introduced to quantify the signal in the on–source
region when compared to the background ascertained from the reference regions.
The core of the modified on–off method is the statement of the null hypothesis
that a previously identified source of interest attains intensity of chosen
value $\beta$, i.e. $N_{\mathrm{on}}=\alpha\beta N_{\mathrm{off}}$ Nosek . A
level of significance for the rejection of the assumption of constant source
activity is given by Nosek ; Stefanik

$S_{\mathrm{LM}}=s\sqrt{2}\left\\{N_{\mathrm{on}}\ln{X_{1}}+N_{\mathrm{off}}\ln{X_{2}}\right\\}^{\frac{1}{2}},$
(1)

for the asymptotic Li–Ma statistics. Here, the logarithmic arguments are
$X_{1}=\frac{1+\alpha\beta}{\alpha\beta}\frac{N_{\mathrm{on}}}{N_{\mathrm{on}}+N_{\mathrm{off}}}$
and
$X_{2}=(1+\alpha\beta)\frac{N_{\mathrm{off}}}{N_{\mathrm{on}}+N_{\mathrm{off}}}$.
The $s$–term in Eq. (1) is either $+1$ or $-1$ for observations of an excess
($S_{\mathrm{LM}}>0$) or deficit ($S_{\mathrm{LM}}<0$) of events,
respectively.

Li–Ma significance given in Eq. (1) asymptotically follows the normal
distribution with zero mean and unit variance Nosek . Therefore, any
inconsistency between the sample distributions of $S_{\mathrm{LM}}$ and the
reference standard Gaussian distribution should be regarded as a sign of
change in the tested $\gamma$–ray intensity. It is also worth noting, that the
test of the null hypothesis can be advantageously inverted to derive
confidence intervals for the source parameter $\beta$ at a given level of
significance. Sequence of such confidence intervals then carries information
about changes of the observed $\gamma$–ray activity.

The correct treatment of Poisson observables is inherent in the modified
on–off method which straightforwardly exploits just the numbers of detected
photons. Its other benefit is that systematic uncertainties are fully
accounted for and are already included in the confidence intervals for the
source parameter. The application of the method is thus very simple yet the
range of possible utilizations and the interpretations of results remain
utmost general.

##  2 Data analysis

Data sets used in our study comprise results of observations of PSR
B1259–63/LS 2883 by the H.E.S.S. telescope array during periastron passages in
2004 hess_2005 , 2007 hess_2009  and 2011 hess_2013 . Presented numbers of on–
and off–source counts together with values of the parameter $\alpha$ were used
as input. Measured counts in each studied observational period ranged between
over one hundred to several thousands. Data were taken in configurations with
different numbers of active telescopes.

We examined changes of the observed source intensity deduced from the H.E.S.S.
data by the means of confidence intervals for the source parameter $\beta$.
These intervals were constructed at a $3\sigma$ ($\approx 99.7\%$) level of
significance for the sequence of monthly triplets
$(N_{\mathrm{on}},\leavevmode\nobreak\ N_{\mathrm{off}},\leavevmode\nobreak\
\alpha)$ using Eq. (1) such that
$|S_{\mathrm{LM}}(N_{\mathrm{on}},\leavevmode\nobreak\
N_{\mathrm{off}},\leavevmode\nobreak\ \alpha;\leavevmode\nobreak\ \beta)|<3$.
Furthermore, we combined the monthly data into larger sets according to
calendar years and calculated the corresponding confidence intervals. The
results are visualized in Fig. 1 where the intervals are depicted as
MJD–ordered time series. Several non–overlapping monthly confidence intervals
(black segments) indicate variations of the intensity observed from the
source. Annual observations provide intervals (hatched bands) which are
consistent altogether. Confidence interval for 2011 campaign depicted in Fig.
1 was obtained with the modified data set (see Tab. 1 and Section 3).

Figure 1:  MJD–ordered sequence of confidence intervals at a $3\sigma$ level
for the source parameter $\beta$ obtained for the $\gamma$–ray binary PSR
B1259–63/LS 2883 hess_2005 ; hess_2009 ; hess_2013 . Segments with points
represent monthly intervals, hatched bands denote yearly intervals. The
horizontal dashed line indicates the average value of the source parameter
taken over all measurements. Interval for 2011 observations was estimated by
accounting for the energy spectrum presented by the H.E.S.S. collaboration,
for details see Tab. 1 and corresponding text.  Figure 2:  Confidence
intervals at a $1\sigma$ level as a function of observation date with respect
to the periastron occurrence which is marked by vertical dotted line. The
horizontal dashed line denotes the reference value corresponding to the
background–only dominated on–source region. Observations in 2011,
$\langle\beta_{-},\leavevmode\nobreak\ \beta_{+}\rangle=\langle
3.57,\leavevmode\nobreak\ 4.44\rangle$, are indicated by the arrow. For
further details see caption to Fig. 1.

The observed data can be further depicted in a plot for confidence intervals
at a $1\sigma$ level arranged according to the proximity of observations to
the periastron date. In Fig. 2, indications of a double–peak structure can be
recognized with a local minimum around the date of periastron.

In order to quantify the trends of the source activity emerging in the
sequence of confidence intervals we also performed test with the 2004 and 2007
data for the assumption of constant intensity. The source activity to be
tested was chosen as an average value of the ratio of observed and expected
on–source events over the observational periods in 2007, i.e.
$\beta_{07}=\langle N_{\mathrm{on}}/\alpha N_{\mathrm{off}}\rangle=1.21$. Our
results of testing changes of the source activity are presented as a
quantile–quantile (QQ) plot. For this purpose, Li–Ma significances were
calculated according to Eq. (1) and then assigned to the quantiles of the
standard normal distribution $N(0,1)$. The quantiles were chosen to be
$k/(m+1)$ where $m$ denotes the number of observed events and $k=1\dots m$.
The observational time sequence is indicated as increasing sizes of markers.

Figure 3:  QQ–plot of asymptotic Li–Ma significances for the $\gamma$–ray
events detected from the direction of PSR B1259–63/LS 2883 hess_2005 ;
hess_2009 . Average source parameter $\beta_{07}=1.21$ calculated from 2007
measurements was tested. The chronological sequence of observations is given
by increasing sizes of markers. The dashed line denotes the reference
$45^{\circ}$ line through the origin.

In Fig. 3, a QQ–plot of the distribution of Li–Ma significances is shown. The
dashed line denotes the reference $45^{\circ}$ line through the origin on
which the data points lie provided the hypothesis of constant source intensity
of given value is true. Most of the measurements conducted in 2007 (squares)
are approximately consistent with the reference line apart from the minor
downward shift which indicates that the value of intensity during these
periods was slightly lower than the chosen one. The only exceptional result in
2007 is obtained for the June data set which yields significances for the
rejection of a steady source activity assumption $S_{\mathrm{LM}}>3.6$. On
contrary, the 2004 sample significances (circles) exhibit considerable
dispersion from the reference line. Different skewness of the distributions of
sample significances compared to the reference standard normal distribution
hints that the source intensity changed. At two occasions a significant excess
of source activity is recognized with the significances $S_{\mathrm{LM}}>3.6$.
In one case a considerable deficit of events is observed with
$S_{\mathrm{LM}}\approx-3$.

##  3 Discussion

Table 1:  Numbers of detected on– and off–source events ($E>500$ GeV) and the on–off parameter $\alpha$ for the H.E.S.S. observations of PSR B1259–63/LS 2883 during 2011 hess_2013 . $\langle\beta_{-},\leavevmode\nobreak\ \beta_{+}\rangle_{3\sigma}$ are confidence intervals for the source parameter $\beta$ constructed at a $3\sigma$ level. Asterisk denotes the modified data set.  Period | $N_{\mathrm{on}}$ | $N_{\mathrm{off}}$ | $\alpha$ | $\langle\beta_{-},\leavevmode\nobreak\ \beta_{+}\rangle_{3\sigma}$
---|---|---|---|---
2011 | 112 | 365 | 0.077 | $\langle 2.85,\leavevmode\nobreak\ 5.46\rangle$
2011* | 80 | 365 | 0.077 | $\langle 1.93,\leavevmode\nobreak\ 4.07\rangle$

The H.E.S.S. collaboration stated that there is evidence for variability of
integral flux above the energy threshold 380 GeV from PSR B1259–63/LS 2883 in
2004 hess_2005  and also above 620 GeV in 2007 hess_2009 . Moreover, the
overall flux level was found to correspond to the one measured in 2004.
Observations in 2005 and 2006, i.e. during the epoch between two consecutive
periastron passages, did not provide significant detection of the source thus
indicating its variable character hess_2009 . No variability on the time scale
of individual nights during 2011 was proved hess_2013 . It was further stated
that the 2011 flux above the energy threshold 1 TeV is consistent with the
flux levels obtained in 2004 and 2007.

In the framework of the modified on–off method, we arrive at similar results
with regard to the evidence for variability of observed $\gamma$–ray intensity
in 2004 and 2007, see Figs. 1 and 3. Likewise, increase of the source emission
before and after the periastron visible in Fig. 2 agrees with the results of
the lightcurve analysis hess_2009 ; hess_2013 . Interestingly, our results
suggest that the source activity observed in 2011 was not consistent with the
overall intensity in previous campaigns. First row in Tab. 1 shows that the
$99.7\%$ confidence interval for 2011 measurements is lying considerably
higher than any other one deduced from earlier observations, compare with Fig.
1. This translates also in the significance $S_{\mathrm{LM}}=9.6$ with which
one can reject the assumption of the constant source intensity. Here, the
source parameter $\beta_{47}=1.25$ to be tested was chosen as an average value
taken from the 2004 and 2007 campaigns.

These results are not contradictory, however. Firstly, due to the smaller data
set the results of the 2011 measurements are less constraining than the
previous ones hess_2013 . Secondly, we worked also with events with energies
below 1 TeV, the threshold above which H.E.S.S. determined the integral flux
hess_2013 . Finally, it is worth noting that the H.E.S.S. experiment measured
considerably increased flux of photons at energies around 500 GeV, see Fig.1
in hess_2013 .

In an attempt to explain the inconsistencies we changed the number of photons
observed in the first interval of the energy spectrum. To this end, we assumed
the same simple power law spectral shape with the index $\Gamma=2.92$ as
reported in the H.E.S.S. analysis hess_2013 . The modified number of photons
in the lowest energy bin was set as consistent with the power law derived from
higher energy bins. The resulting 2011 confidence interval for the source
parameter $\beta$ stemming from this modification is given in the second row
in Tab. 1. This interval still indicates that the observed source activity
above the energy threshold of roughly 500 GeV was in 2011 above the
intensities deduced from previous campaigns provided the background remained
stable throughout all observation periods.

In view of the 2011 detection of a HE $\gamma$–ray flare by the _Fermi_ –LAT
Abdo ; Tam , it may be suggestive to relate this flare to the higher level of
the VHE activity as seen in our analysis. However, when testing the
‘Pre–flare’ and ‘Flare’ data sets of H.E.S.S. hess_2013  for deviations of the
source intensity from its average value, the test significances for both
periods are $|S_{\mathrm{LM}}|<1$. Thus, we cannot reject the assumption of
constant intensity during 2011 based on these data and, accordingly, the
connection with the GeV flare is questionable.

The increase of the observed VHE activity may be attributable to the truly
intrinsically variable emission or it may be some effect of the detection
process. We cannot exclude the latter possibility without further details on
the used data, particularly on the exact detection thresholds in each
observational campaign. Nevertheless, we have no firm indication of
inconsistency of the source intensity above 1 TeV in 2011 when compared to the
previous observations. Hence, we comply with the H.E.S.S. conclusions about
the overall agreement of integral fluxes above this energy.

##  4 Conclusions

We used the modified on–off method to judge whether significant changes of the
$\gamma$–ray intensity from the direction of binary PSR B1259–63/LS 2883 were
observed during its three periastron passages. Within our method, the problem
of variable source activity can be treated using solely the numbers of
detected events.

Temporal evolution of the observed source intensity is clearly recognized on
the account of the modified on–off technique. The results of our analysis of
the data gathered by the H.E.S.S. experiment clearly show that changes of the
source activity occurred during 2004 and 2007 between individual months while
the overall intensities in both seasons stayed in agreement. Our findings on
the data obtained above approximately 500 GeV in 2011 show that the observed
source intensity was enhanced in this season with respect to earlier
measurements.

This work was supported by the Czech Science Foundation grant 14-17501S.

## References

  * (1) S. Johnston et al., ApJ 387, L37–L41 (1992)
  * (2) H.E.S.S. Collaboration, A&A 551, A94 (2013)
  * (3) F. Aharonian et al., A&A 442, 1–10 (2005)
  * (4) F. Aharonian et al., A&A 507, 389–396 (2009)
  * (5) A. A. Abdo et al., ApJ Lett. 736, L11 (2011)
  * (6) P. H. T. Tam et al., ApJ Lett. 736, L10 (2011)
  * (7) D. Berge, S. Funk, J. Hinton, A&A 446, 1219–1229 (2007)
  * (8) T.P. Li, Y.Q. Ma, ApJ 272, 317–324 (1983)
  * (9) D. Nosek, S. Stefanik, J. Noskova, ICRC 2013, arXiv:1309.6476
  * (10) S. Stefanik, D. Nosek, Nucl. Phys. B-Proc. Sup. 256, 258–263 (2014)


# Correlating space, wavelength, and polarization of light: Spatio-Spectral
Vector Beams

Lea Kopf,1 Rafael Barros,1 and Robert Fickler1
1 Tampere University, Photonics Laboratory, Physics Unit, Tampere, FI-33720,
Finland

###### Abstract

Increasing the complexity of a light field through the advanced manipulation
of its degrees of freedom (DoF) provides new opportunities for fundamental
studies and technologies. Correlating polarization with the light’s spatial or
spectral shape results in so-called spatial or spectral vector beams that are
fully polarized and have a spatially or spectrally varying polarization
structure. Here, we extend the general idea of vector beams by combining both
approaches and structuring a novel state of light in three non-separable
DoF’s, i.e. space, wavelength, and polarization. We study in detail their
complex polarization structure, show that the degree of polarization of the
field is only unveiled when the field is narrowly defined in space and
wavelength, and demonstrate the analogy to the loss of coherence in non-
separable quantum systems. Such light fields allow fundamental studies on the
non-separable nature of a classical light field and new technological
opportunities, e.g. through applications in imaging or spectroscopy.

##  I Introduction

Increasing the complexity of a light field and the control of different
degrees of freedom (DoF’s) is beneficial for advancing research and
technology. Increasingly complex structures realized by combining several
DoF’s have been studied in a myriad of experiments over the last decades and
the enhanced understanding of the interplay of the DoF’s has already enabled
novel photonic technologies [1, 2, 3, 4, 5, 6]. Initially, many experiments
have studied structuring transverse light fields as well as shaping the
temporal profile of pulses, both in its scalar forms, i.e. with a uniform
polarization structure [7, 8]. The complexity of the light field’s structure
was further increased by including the polarization domain leading to beams
with spatially non-uniform polarization distributions, i.e. spatial vector
beams, as well as a temporally varying polarization vectors across the pulse
duration. Over the last years, this approach has been extended to combine all
DoF’s, e.g. the study of advanced spatio-temporal pulses of vectorial light
fields [9], or the adaptation of established concepts such as complex
transverse spatial fields to the spatio-temporal domain [10, 11].
Interestingly, the focus in most of the research efforts has been to generate
complex polarization and spatial patterns over the temporal domain of light
with much less attention to studying structured light fields in the time’s
complementary DoF, i.e. the spectral domain. Correlating polarization with the
wavelength of a light pulse, for example, can be used for advanced sensing and
pulse characterization schemes [12, 13, 14, 15, 16].

Figure 1: Conceptual idea. A simple generation scheme consisting only of a
birefringent crystal, a quarter wave-plate (QWP), and a vortex-retarder
generates a complex light field where polarization, wavelength, and space are
correlated. The color in the insets symbolizes the different wavelength
components of the light field. a) Poincaré sphere displaying that the light
field is diagonally polarized. b) The Poincaré sphere of the light field shows
that each wavelength component has a different linear polarization. c) The
higher-order Poincaré sphere [17, 18] illustrates that each wavelength has a
different spatial polarization pattern.

Here, we extend the idea of spatial and spectral vector beams by combining all
three DoF’s, namely polarization, space, and wavelength of the light field. We
term such light pulses spatio-spectral vector beams (SSVB’s), which are light
fields with a varying polarization structure in space as well as wavelength as
shown in Fig. 1 a). Using a simple setup consisting of only three optical
elements placed along a single beam line, we are able to generate highly
complex pulses of light for which any given wavelength (or transverse angular
position) shows a different spatial (or spectral) polarization pattern. We
further show that all three DoF’s are non-separable and the complex vectorial
nature of the light field can only be observed when all DoF’s are resolved. By
integrating over the light’s transverse spatial extend or wavelength spectrum
(or both), the beam would be seemingly unpolarized. Finally, we detail the
analogy of such complex structured light fields to a tripartite
Greenberger–Horne–Zeilinger (GHZ) state, which describes three entangled
particles in quantum optics. Here, the DoF’s of the classical light field act
as the different quantum states in the GHZ description. We expect that the
complexity of SSVB’s combined with the ease of our generation technique will
trigger novel studies on structured light fields and advance their
applications in novel spectroscopy, imaging, or sensing technologies.

##  II Method

We consider a SSVB of the form

$\vec{E}(\vec{r},\lambda)=\sqrt{\frac{I_{0}}{2}}\left[S_{1}(\vec{r})F_{1}(\lambda)\hat{\epsilon}_{1}+S_{2}(\vec{r})F_{2}(\lambda)\hat{\epsilon}_{2}\right]\,,$
(1)

where $\vec{r}$ is the position vector, $\lambda$ is the wavelength, $I_{0}$
is the total field intensity, $S_{1,2}$ and $F_{1,2}$ are normalized spatial
and spectral basis functions, respectively, and $\hat{\epsilon}_{1,2}$ are
unit polarization vectors. We assume orthogonal polarization vectors with
$\hat{\epsilon}_{i}\cdot\hat{\epsilon}_{j}=\delta_{ij}$, but allow for the
spatial and spectral functions to have non-vanishing overlap. Thus, the
function can describe a light field with non-separable polarization, space,
and wavelength for zero overlap.

As a measure of the degree of tripartite non-separability of field (1), we
consider the degree of polarization (DoP) of the light field. As we detail in
appendix A, the DOP is given by

$D=\sqrt{\frac{\left(\Lambda^{S}_{11}\Lambda^{F}_{11}-\Lambda^{S}_{22}\Lambda^{F}_{22}\right)^{2}+4|\Lambda^{S}_{12}|^{2}|\Lambda^{F}_{12}|^{2}}{\left(\Lambda^{S}_{11}\Lambda^{F}_{11}+\Lambda^{S}_{22}\Lambda^{F}_{22}\right)^{2}}}\,,$
(2)

where $\Lambda^{S}_{ij}=\langle S_{i},S_{j}\rangle_{\Omega_{S}}$ and
$\Lambda^{F}_{ij}=\langle F_{i},F_{j}\rangle_{\Omega_{F}}$ ($j=1,2$) are the
inner products between the spatial and spectral basis functions over finite
regions $\Omega_{S}$ and $\Omega_{F}$ in space and wavelength, respectively.
Experimentally, the DOP can be retrieved from the Stokes parameters resulting
from intensity measurements integrated over $\Omega_{S}$ and $\Omega_{F}$.

Note that when integrating over the complete spatial and spectral DoF’s
($\Omega_{S}=\mathbb{R}^{2}$ and $\Omega_{F}=\mathbb{R}$), Eq. (2) yields
$D=|\Lambda^{S}_{12}\Lambda^{F}_{12}|$. Thus, when either the spatial or
spectral basis functions are orthogonal, the measured light field is seemingly
unpolarized. The same argument holds for space or wavelength-only
measurements, which yield reduced degrees of spatial and spectral coherence,
respectively. Only joint measurements of the three DoF’s unveil the full
coherence of the input field, which is a signature of tripartite non-
separability. The apparent decoherence of the field (1) upon partial tracing
of one DoF is a by-product of its mathematical similarity to the GHZ states
[19, 20, 21]. Classically, such a tripartite non-separable behaviour can also
be explained in an intuitive way. At any point in space where both spatial
basis functions are non-vanishing, the polarization is spectrally non-uniform,
i.e. a spectral vector beam [15]. Similarly, at any wavelength present in both
spectral basis functions, the polarization is spatially non-uniform, i.e. a
spatial vector beam [22]. Lastly, measuring the linear polarization state of
the light field leads to a strong correlation between spatial and spectral
structures. In other words, specifying one DoF can project the remaining two
in a bipartite non-separable field, which we will explore in more detail in
the results section. We note that a related behaviour has been observed using
multiple beams, each with a spatially varying polarization [23].

###  II.1 Experimental Generation

To create an SSVB, we use the process shown in Fig. 1. First, Fourier-limited
laser pulses with a duration of $\tau=220$ fs and centered at a wavelength of
780 nm are linearly polarized. The polarization is chosen to be diagonal with
respect to the fast and slow axes of a 2 mm long birefringent BaB2O4 (BBO)
crystal with the optical axis oriented at $23.4^{\circ}$ from the propagation
direction. The propagation through the birefringent crystal coherently splits
each pulse into two trailing pulses with equal intensities, linearly polarized
along the crystal’s fast ($f$) and slow ($s$) axes, producing a spectral
vector beam [15]. The field at this stage is

$\displaystyle\vec{E}(\vec{r},\lambda)=$
$\displaystyle\sqrt{\frac{I_{0}}{2}}S_{0}(\vec{r})F_{0}(\lambda)\left[e^{i\pi
c\delta t/\lambda}\hat{e}_{f}+\right.$ $\displaystyle\left.e^{-i\pi c\delta
t/\lambda}\hat{e}_{s}\right]\,,$ (3)

where $S_{0}(\vec{r})$ is a Gaussian transverse spatial profile,
$F_{0}(\lambda)$ is the spectral function, $c$ is the speed of light, and
$\delta t$ is the birefringent time delay. A subsequent quarter-wave plate
oriented at 45∘ with respect to the fast axis of the crystal transforms the
field’s polarization components $\hat{e}_{f}$ and $\hat{e}_{s}$ to left and
right circular polarizations, i.e. $\hat{e}_{R}$ and $\hat{e}_{L}$
respectively. We set the time delay to $\delta t\approx\tau$ by choosing the
appropriate crystal length and fine-tuning the crystal rotation, which ensures
that the resulting linear polarization state in the wavelength domain rotates
approximately once within one spectral bandwidth. This beam is non-separable
in wavelength and polarization, but still has only a single transverse
structure, i.e. a Gaussian mode.

Next, we correlate the spatial DoF with polarization using an m=1 zero-order
vortex half-wave retarder. The vortex retarder has a constant half-wave
retardance across the clear aperture with the fast axis oriented along
$\theta=\phi/2$, where $\phi$ is the azimuthal angle. Upon transmission, the
left and right circular polarization components change their handedness and
acquire an azimuthally varying phase from 0 to 2$\pi$ [24] resulting in

$\displaystyle\vec{E}(\vec{r},\lambda)=$
$\displaystyle\sqrt{\frac{I_{0}}{2}}S_{0}(\vec{r})F_{0}(\lambda)\left[e^{i\phi}e^{i\pi\tau
c/\lambda}\hat{e}_{L}+\right.$ $\displaystyle\left.e^{-i\phi}e^{-i\pi\tau
c/\lambda}\hat{e}_{R}\right]\,.$ (4)

Figure 2: Spectrally resolved measurements of the generated SSVB’s. a) The
measurement setup filters space by using a spatial mask. Subsequently, it
filters polarization using a half and quarter wave plate followed by a
polarizing beam splitter (PBS). The final data acquisition is done in the
spectral domain using a spectrometer. b) Polarization-resolved spectral
measurements with a spatial mask rotation angle of 80∘. c) Reconstructed
polarization ellipses over wavelength. The colors refer to the polarization
ellipse rotation angle and are included for better visibility of the changes
in polarization. d) The SSVB at different spatial mask rotation angles.

The field (4) is exactly of the intended form of the SSVB (1), provided that
we identify $S_{1,2}(\vec{r})=S_{0}(\vec{r})e^{\pm\phi}$ and
$F_{1,2}(\lambda)=F_{0}(\lambda)e^{\pm i\pi\tau c/\lambda}$. Note that while
the spatial basis functions are mutually orthogonal ($\langle
S_{1},S_{2}\rangle_{\mathbb{R}^{2}}$=0), the overlap between the spectral
basis functions is non-vanishing ($\langle F_{1},F_{2}\rangle_{\mathbb{R}}\sim
1/e^{2}$). The spectral overlap limits the quality of our SSVB, and is the
cost for the simple method to produce it. We note that the aforementioned
limitation can be overcome by more advanced pulse shaping techniques to
generate orthogonal temporal or spectral modes [7, 25], however, they lead to
the same arguments as presented here.

##  III Results and Discussion

###  III.1 Characterization

To characterize the SSVB (4), we perform both individual and joint
measurements of its polarization, wavelength, and spatial profile. In the
first set of measurements, we use the experimental setup depicted in Fig. 2
a). We send the SSVB through a spatial mask comprising two diametrically
opposite 0.2 mm wide apertures centred at the azimuthal angles $\phi_{p}$ and
$\phi_{p}+\pi$, which roughly corresponds to two 17∘ wide slits and projects
the field onto approximately $\vec{E}(\phi_{p},\lambda)$.

We characterize the resulting field with spectrally-resolved polarization
tomography, where the Stokes parameters are measured with a spectrometer, as
shown for $\phi_{p}=80^{\circ}$ in Fig. 2 b) with the corresponding
wavelength-dependent polarization ellipses in c). In Fig. 2 d), we show
different rotation angles of the spatial mask, which display a linear shift of
the spectral-polarization structure. Note that the measured polarization
states shown in Fig. 2 c-d) are not strictly linear, as expected from (4), but
consist of elongated ellipses. This results from the non-zero overlap of the
spectral basis functions in our scheme and experimental imperfections such as
the finite area of the spatial slits.

In the second set of measurements, we project the SSVB onto wavelength instead
of space, using the setup depicted in Fig. 3 a). A monochromator consisting of
a diffraction grating, a mirror, and a tunable rectangular slit selects a
central wavelength $\lambda_{p}$, which projects the field onto approximately
$\vec{E}(\vec{r},\lambda_{p})$. The projection resolution is fixed by the
minimum bandwidth of around 2.1 nm of the monochromator. We characterize the
spectrally filtered field through spatially-resolved polarization tomography,
in which the Stokes parameters are measured with a CMOS camera. In Fig. 3 b),
we show the polarization-resolved images obtained at $\lambda_{p}=785.5\,$nm
with the corresponding spatial vector beam in c). The spatial vector beams
retrieved for different wavelengths are shown in Fig. 3 d) and, to help
visualize the wavelength-dependence of the retrieved field structures, the
related images obtained at a fixed horizontal polarization are additionally
shown in e). The polarized images show a two-lobe structure that rotates
across the spectrum, as expected from Eq. (1).

Figure 3: Spatially resolved measurements of the generated SSVB. a) The
measurement setup filters wavelength and polarization using a monochromator,
wave plates, and a polarizing beam splitter (PBS). The data is acquired in the
spatial domain with a camera. b) The polarization measurements at a wavelength
filter centred at 785.5 nm. c) The reconstructed SSVB from the polarization
measurements shown in b). The color corresponds to the rotation angle of the
polarization ellipses and acts as a guide to the eye. d) The SSVB’s at varying
wavelength filter positions. e) The projection of the vector beam on the
horizontal polarization at varying wavelength filter positions.

We further study that partially resolving the spatial and spectral domains can
be used to control the DoP of the laser beam. To this end, we perform
polarization measurements for different spatial and spectral bandwidths. We
repeat the measurements shown in Fig. 3 a) with different spectral resolutions
at $\lambda_{p}=779.75$ nm by changing the width of the rectangular slit in
the monochromator. Then, a computational spatial mask with varying angular
slit widths is added to the images in post-processing, simulating the physical
mask used in Fig. 2 a). From the measured Stokes parameters for each setting,
we retrieve the DoP shown in Fig. 4. By increasing the bandwidth of the
wavelength filter or the slit width of the spatial mask, thereby averaging
over the two DoF’s, the beam’s DoP decreases substantially. The initial DoP of
0.95, measured with a wavelength filter bandwidth of 2.1 nm and a spatial slit
mask 0.16 rad wide, drops down to 0.04 when both wavelength and spatial
filters are completely removed, showing excellent agreement with the
simulations (see Fig. A1 in the appendix). Without resolving all DoF’s, we
thus have a seemingly unpolarized light field. Moreover, note that the minimum
DoP obtained for spatially-resolved measurements is non-negligible ($\sim
0.2$), which agrees with the estimated overlap between the spectral basis
functions of $1/e^{2}$.

Figure 4: The DoP as a function of the spatial and spectral filter bandwidths.
As both filters are enlarged, the DoP decreases due to the non-separability of
space, wavelength, and polarization.

###  III.2 Analogy with quantum mechanics

As has been already recognized, classical non-separable light fields share a
mathematical similarity with entangled quantum states [26, 27, 28]. In
contrast to manifold discussions of two non-separable DoF’s of light and the
analogy to Bell states [29, 30], SSVB’s are mathematically similar to the
description of three-particle GHZ states. The GHZ state has the form

$\ket{\Psi}=\frac{1}{\sqrt{2}}\left(\ket{000}+\ket{111}\right)\,,$ (5)

where $\ket{0}$ and $\ket{1}$ denote the two possible states of each quantum
system expressed in terms of qubits in the computational or $z$-basis. The
mutually unbiased bases to $z$ are the $x$\- and $y$-bases, defined in the
Dirac notation as $\ket{\pm x}=(\ket{0}\pm\ket{1})/\sqrt{2}$ and $\ket{\pm
y}=(\ket{0}\pm i\ket{1})/\sqrt{2}$, respectively.

To demonstrate the non-separability of the SSVB, we closely follow the GHZ
argument [20], which focuses on the simultaneous measurements of all three
qubits in the $x$-basis ($xxx$ measurement). In the original GHZ argument, the
two possible measurement outcomes in each basis are ascribed with the values
$+1$ and $-1$. Together with the assumption of a local realistic theory, it is
possible to show that classical correlations require the product of the
measurement outcomes in the $x$-basis to be $+1$, while for quantum
correlations the outcome will be $-1$. By adapting these ideas to SSVB’s and
showing a similar behaviour, we can use the GHZ argumentation to demonstrate
the non-separability of the classical light field.

We assign the notation (5) to the SSVB by representing the spatial, spectral
and polarization modes as the first, second, and third entries of the state
vector, respectively. For the spatial DoF, the $z$ basis is formed by the
spatial functions $S_{1,2}(\vec{r})=S_{0}(\vec{r})e^{\pm i\phi}$, which
contains counter-rotating spiral wavefronts carrying opposite orbital angular
momenta [31], and indistinguishable intensity profiles. On the other hand, the
$x$ ($y$) basis modes have two-lobe structures oriented along $0^{\circ}$ and
$90^{\circ}$ ($45^{\circ}$ and $-45^{\circ}$) [32], and can be distinguished
by narrow slits oriented along these directions [33, 34]. We implement these
slits by using the spatial mask already presented in the measurements of Fig.
2.

In the wavelength domain, the $z$ basis is formed by the spectral functions
$F_{1,2}(\lambda)=F_{0}(\lambda)e^{\pm i\pi\tau c/\lambda}$, which cannot be
distinguished by measurements of spectral amplitude. On the other hand, the
modes composing the $x$ basis are $F_{1x}\propto F_{0}(\lambda)\cos(\pi\tau
c/\lambda)$ and $F_{2x}\propto F_{0}(\lambda)\cos(\pi\tau c/\lambda-\pi/2)$,
while those of the $y$ basis are $F_{1y}\propto F_{0}(\lambda)\cos(\pi\tau
c/\lambda-\pi/4)$ and $F_{2y}\propto F_{0}(\lambda)\cos(\pi\tau
c/\lambda-3\pi/4)$. In analogy to the spatial domain, we measure in the $x$
and $y$ bases by choosing, in post-processing, four wavelengths spaced by
$1/\tau$ at which the intensity difference between the horizontal/vertical or
diagonal/anti-diagonal polarization components are maximized.

Lastly, for the polarization DoF, we define the $z$ basis vectors as the left
and right circularly polarized components of the light field. Consequently,
the $x$ basis corresponds to horizontal/vertical linear polarizations, and the
$y$ basis to diagonal/anti-diagonal linear polarizations. A graphical overview
of the chosen bases is given in Fig. B2 in the appendix.

Figure 5: Measurements to show the non-separability of the SSVB. a)
Measurement probabilities where one DoF is measured in the $x$-basis and the
other two are measured in the $y$-basis. The blue bars correspond to the
expected states in the mathematical description of a GHZ state, the red bars
correspond to nonexistent states. b) Measured detection probabilities when all
DoF’s are measured in the $x$-basis. The measurements in the top row nicely
show, that the description of the non-separable state (middle row) is dominant
over the expected result for separable fields (lower row).

Fig. 5 shows the measurement results for all terms of the normalized state (5)
written in the $x$ and $y$ bases. More information on the normalization
procedure can be found in the appendix. Clearly, the $xxx$ measurements show
that the product of the three outcomes mostly yields $+1$, with a small
contribution of $-1$ due to experimental limitations and imperfections. We
regard this result as an interesting way to show the non-separability between
polarization, spatial profile and wavelength by exploring the mathematical
isomorphism with the GHZ state.

##  IV Conclusion

In this article, we introduced and experimentally realized a light field with
a complex polarization structure in space as well as wavelength - a spatio-
spectral vector beam. Our measurements show that the SSVB is well defined in
its three DoF’s if, and only if, they are measured jointly, while individual
measurements result in apparent incoherence, i.e. a reduced degree of
polarization, similar to the loss of coherence in non-separable quantum
systems. In the future, it will be interesting to extend the concept further
to more complex spatio-spectral patterns, e.g. by including spatial Poincaré
beams [35] or studying SSVB’s in the context of optical skyrmions [36]. Not
only do SSVB’s contribute to a better understanding of the complex interplay
between different DoF’s, they might also enable applications benefiting from
the complex correlations and simple generation. In contrast to earlier work,
where spatial and spectral vector beams have been used for high-speed tracking
of changes in space [37] and spectrum [15], the presented light fields could
allow a fruitful combination of the two leading to novel schemes for advanced
sensing methods in all optical domains simultaneously. In addition, our
results might be applied in quantum optics experiments when complex hyper-
entangled states are utilized [38, 39]. Finally, it will be interesting to
study the possible benefits of SSVB’s in nonlinear light-matter interactions
of structured light [40].

## References

  * Rubinsztein-Dunlop  _et al._ [2016] H. Rubinsztein-Dunlop, A. Forbes, M. V. Berry, M. R. Dennis, D. L. Andrews, M. Mansuripur, C. Denz, C. Alpmann, P. Banzer, T. Bauer, E. Karimi, L. Marrucci, M. Padgett, M. Ritsch-Marte, N. M. Litchinitser, N. P. Bigelow, C. Rosales-Guzmán, A. Belmonte, J. P. Torres, T. W. Neely, M. Baker, R. Gordon, A. B. Stilgoe, J. Romero, A. G. White, R. Fickler, A. E. Willner, G. Xie, B. McMorran, and A. M. Weiner, Roadmap on structured light, Journal of Optics 19, 013001 (2016).
  * Forbes  _et al._ [2021] A. Forbes, M. de Oliveira, and M. R. Dennis, Structured light, Nature Photonics 15, 253 (2021).
  * Piccardo  _et al._ [2021] M. Piccardo, V. Ginis, A. Forbes, S. Mahler, A. A. Friesem, N. Davidson, H. Ren, A. H. Dorrah, F. Capasso, F. T. Dullo, B. S. Ahluwalia, A. Ambrosio, S. Gigan, N. Treps, M. Hiekkamäki, R. Fickler, M. Kues, D. Moss, R. Morandotti, J. Riemensberger, T. J. Kippenberg, J. Faist, G. Scalari, N. Picqué, T. W. Hänsch, G. Cerullo, C. Manzoni, L. A. Lugiato, M. Brambilla, L. Columbo, A. Gatti, F. Prati, A. Shiri, A. F. Abouraddy, A. Alù, E. Galiffi, J. B. Pendry, and P. A. Huidobro, Roadmap on multimode light shaping, Journal of Optics 24, 013001 (2021).
  * He  _et al._ [2022] C. He, Y. Shen, and A. Forbes, Towards higher-dimensional structured light, Light: Science & Applications 11, 205 (2022).
  * Shen  _et al._ [2022a] Y. Shen, Q. Zhan, L. G. Wright, D. N. Christodoulides, F. W. Wise, A. E. Willner, Z. Zhao, K.-h. Zou, C.-T. Liao, C. Hernández-García, _et al._ , Roadmap on spatiotemporal light fields, arXiv preprint arXiv:2210.11273 (2022a).
  * Wan  _et al._ [2023] Z. Wan, H. Wang, Q. Liu, X. Fu, and Y. Shen, Ultra-degree-of-freedom structured light for ultracapacity information carriers, ACS Photonics 10, 2149 (2023), https://doi.org/10.1021/acsphotonics.2c01640 .
  * Weiner [2000] A. M. Weiner, Femtosecond pulse shaping using spatial light modulators, Review of scientific instruments 71, 1929 (2000).
  * Dickey [2018] F. M. Dickey,  _Laser beam shaping: theory and techniques_ (CRC press, 2018).
  * Chen  _et al._ [2022] L. Chen, W. Zhu, P. Huo, J. Song, H. J. Lezec, T. Xu, and A. Agrawal, Synthesizing ultrafast optical pulses with arbitrary spatiotemporal control, Science Advances 8, eabq8314 (2022).
  * Jhajj  _et al._ [2016] N. Jhajj, I. Larkin, E. W. Rosenthal, S. Zahedpour, J. K. Wahlstrand, and H. M. Milchberg, Spatiotemporal optical vortices, Phys. Rev. X 6, 031037 (2016).
  * Cao  _et al._ [2022] Q. Cao, P. Zheng, and Q. Zhan, Vectorial sculpturing of spatiotemporal wavepackets, APL Photonics 7 (2022).
  * Sano  _et al._ [1980] K. Sano, K. Okada, and K. Hashimoto, Simple optical wavelength meter in 700-1200 nm wavelength region, Electronics Letters 24, 912 (1980).
  * del Toro Iniesta [2003] J. C. del Toro Iniesta,  _Introduction to spectropolarimetry_ (Cambridge university press, 2003).
  * Aspnes [2014] D. Aspnes, Spectroscopic ellipsometry—past, present, and future, Thin Solid Films 571, 334 (2014).
  * Kopf  _et al._ [2021] L. Kopf, J. R. D. Ruano, M. Hiekkamäki, T. Stolt, M. J. Huttunen, F. Bouchard, and R. Fickler, Spectral vector beams for high-speed spectroscopic measurements, Optica 8, 930 (2021).
  * Jolly  _et al._ [2021] S. W. Jolly, O. Gobert, and F. Quéré, Spatio-spectral characterization of ultrashort laser pulses with a birefringent delay line, OSA Continuum 4, 2044 (2021).
  * Holleczek  _et al._ [2010] A. Holleczek, A. Aiello, C. Gabriel, C. Marquardt, and G. Leuchs, Poincaré sphere representation for classical inseparable bell-like states of the electromagnetic field, arXiv 1007.2528, 10.48550/arXiv.1007.2528 (2010).
  * Milione  _et al._ [2011] G. Milione, H. I. Sztul, D. A. Nolan, and R. R. Alfano, Higher-order poincaré sphere, stokes parameters, and the angular momentum of light, Phys. Rev. Lett. 107, 053601 (2011).
  * Bouwmeester  _et al._ [1999] D. Bouwmeester, J.-W. Pan, M. Daniell, H. Weinfurter, and A. Zeilinger, Observation of three-photon greenberger-horne-zeilinger entanglement, Physical Review Letters 82, 1345 (1999).
  * Pan  _et al._ [2000] J. W. Pan, D. Bouwmeester, M. Daniell, H. Weinfurter, and A. Zeilinger, Experimental test of quantum nonlocality in three-photon greenberger-horne-zeilinger entanglement, Nature 403, 515– (2000).
  * Greenberger  _et al._ [1989] D. M. Greenberger, M. A. Horne, and A. Zeilinger, Going beyond bell’s theorem, in  _Bell’s Theorem, Quantum Theory and Conceptions of the Universe_, edited by M. Kafatos (Springer Netherlands, Dordrecht, 1989) pp. 69–72.
  * Zhan [2009] Q. Zhan, Cylindrical vector beams: from mathematical concepts to applications, Adv. Opt. Photon. 1, 1 (2009).
  * Balthazar  _et al._ [2016] W. F. Balthazar, C. E. R. Souza, D. P. Caetano, E. F. G. ao, J. A. O. Huguenin, and A. Z. Khoury, Tripartite nonseparability in classical optics, Opt. Lett. 41, 5797 (2016).
  * Biener  _et al._ [2002] G. Biener, A. Niv, V. Kleiner, and E. Hasman, Formation of helical beams by use of pancharatnam–berry phase optical elements, Opt. Lett. 27, 1875 (2002).
  * Monmayrant  _et al._ [2010] A. Monmayrant, S. Weber, and B. Chatel, A newcomer’s guide to ultrashort pulse shaping and characterization, Journal of Physics B: Atomic, Molecular and Optical Physics 43, 103001 (2010).
  * Spreeuw [1998] R. J. C. Spreeuw, A classical analogy of entanglement, Foundations of Physics 28, 361– (1998).
  * Aiello  _et al._ [2015] A. Aiello, F. Töppel, C. Marquardt, E. Giacobino, and G. Leuchs, Quantum-like nonseparable structures in optical beams, New J. Phys. 17, 043024 (2015).
  * Shen  _et al._ [2021] Y. Shen, I. Nape, X. Yang, X. Fu, M. Gong, D. Naidoo, and A. Forbes, Creation and control of high-dimensional multi-partite classically entangled light, Light Sci. Appl. 10, 50 (2021).
  * Karimi and Boyd [2015] E. Karimi and R. W. Boyd, Classical entanglement?, Science 350, 1172 (2015).
  * Töppel  _et al._ [2014] F. Töppel, A. Aiello, C. Marquardt, E. Giacobino, and G. Leuchs, Classical entanglement in polarization metrology, New Journal of Physics 16 (2014).
  * Allen  _et al._ [1992] L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, Orbital angular momentum of light and the transformation of laguerre-gaussian laser modes, Phys. Rev. A 45, 8185 (1992).
  * Padgett and Courtial [1999] M. J. Padgett and J. Courtial, Poincaré-sphere equivalent for light beams containing orbital angular momentum, Optics letters 24, 430 (1999).
  * Fickler  _et al._ [2012] R. Fickler, R. Lapkiewicz, W. N. Plick, M. Krenn, C. Schaeff, S. Ramelow, and A. Zeilinger, Quantum entanglement of high angular momenta, Science 338, 640 (2012).
  * Plick  _et al._ [2015] W. N. Plick, R. Fickler, R. Lapkiewicz, and S. Ramelow, Violation of an extended wigner inequality with high-angular-momentum states, Phys. Rev. A 91, 022124 (2015).
  * Beckley  _et al._ [2010] A. M. Beckley, T. G. Brown, and M. A. Alonso, Full poincaré beams, Opt. Express 18, 10777 (2010).
  * Shen  _et al._ [2022b] Y. Shen, E. C. Martínez, and C. Rosales-Guzmán, Generation of optical skyrmions with tunable topological textures, ACS Photonics 9, 296 (2022b).
  * Berg-Johansen  _et al._ [2015] S. Berg-Johansen, F. Töppel, B. Stiller, P. Banzer, M. Ornigotti, E. Giacobino, G. Leuchs, A. Aiello, and C. Marquardt, Classically entangled optical beams for high-speed kinematic sensing, Optica 2, 864 (2015).
  * Barreiro  _et al._ [2005] J. T. Barreiro, N. K. Langford, N. A. Peters, and P. G. Kwiat, Generation of hyperentangled photon pairs, Phys. Rev. Lett. 95, 260501 (2005).
  * Achatz  _et al._ [2023] L. Achatz, L. Bulla, S. Ecker, E. A. Ortega, M. Bartokos, J. C. Alvarado-Zacarias, R. Amezcua-Correa, M. Bohmann, R. Ursin, and M. Huber, Simultaneous transmission of hyper-entanglement in three degrees of freedom through a multicore fiber, npj Quantum Information 9, 45 (2023).
  * Buono and Forbes [2022] W. T. Buono and A. Forbes, Nonlinear optics with structured light, Opto-Electronic Advances 5, 210174 (2022).
  * Kagalwala  _et al._ [2015] K. H. Kagalwala, H. E. Kondakci, A. F. Abouraddy, and B. E. A. Saleh, Optical coherency matrix tomography, Scientific Reports 5, 15333 (2015).

## Acknowledgements

The authors acknowledge the support of the Academy of Finland through the
Competitive Funding to Strengthen University Research Profiles (decision
301820) and the support of the Photonics Research and Innovation Flagship
(PREIN - decision 320165). LK acknowledges the support of the Vilho, Yrjö and
Kalle Väisälä Foundation of the Finnish Academy of Science and Letters. RB
acknowledges the support of the Academy of Finland through the postdoctoral
researcher funding (decision 349120). RF acknowledges the support of the
Academy of Finland through the Academy Research Fellowship (decision 332399).

##  Appendix A Simulation of the degree of polarization

The degree of polarization (DoP) of a field $\vec{E}(\vec{r},\lambda)$ is
defined as [41]

$D=\frac{1}{\langle I_{0}\rangle}\sqrt{S_{H}^{2}+S_{D}^{2}+S_{R}^{2}}$ (6)

where $\langle I_{0}\rangle$ is the total average intensity, and
$S_{j}=2\langle I_{j}\rangle-\langle I_{0}\rangle$ ($j=H,D,R$) are the so-
called Stokes parameters. The indices H, D, and R denote polarization
projections along the horizontal, diagonal, and right-circular polarizations,
without loss of generality. The average intensities $\langle I_{j}\rangle$ are
given by

$\langle
I_{j}\rangle=\int_{\Omega_{S}}d^{2}\vec{r}\int_{\Omega_{F}}d\lambda\,|\vec{E}(\vec{r},\lambda)\cdot\hat{\epsilon_{j}}|^{2},$
(7)

where $\Omega_{S}$ and $\Omega_{F}$ are the integration regions in the spatial
and spectral domains, respectively. Using the field expressed in Eq. (1) of
the main text, we obtain

$\langle
I_{j}\rangle=\frac{I_{0}}{2}\sum_{i=1,2}\Lambda^{S}_{ii}\Lambda^{F}_{ii}|\Lambda^{P}_{ij}|^{2}+I_{0}\textrm{Re}\left(\\\
\Lambda^{S}_{12}\Lambda^{F}_{12}\Lambda^{P}_{1j}\Lambda^{P*}_{2j}\right)\,,$
(8)

where $\Lambda^{P}_{ij}=\hat{\epsilon}_{i}\cdot\hat{\epsilon}_{j}^{\,*}\,$,
and where

$\displaystyle\Lambda^{S}_{ik}$ $\displaystyle=$
$\displaystyle\int_{\Omega_{S}}d^{2}\vec{r}\,S_{i}(\vec{r})S_{k}^{*}(\vec{r})\,,$
(9) $\displaystyle\Lambda^{F}_{ik}$ $\displaystyle=$
$\displaystyle\int_{\Omega_{F}}d\lambda\,F_{i}(\lambda)F_{k}^{*}(\lambda)\,,$
(10)

with $i,k=1,2$, are the overlaps of the spatial and spectral basis functions
of $\vec{E}(\vec{r},\lambda)$ over $\Omega_{S}$ and $\Omega_{F}$. Finally, by
inserting the Stokes parameters calculated from (8) into (6), we obtain

$D=\sqrt{\frac{\left(\Lambda^{S}_{11}\Lambda^{F}_{11}-\Lambda^{S}_{22}\Lambda^{F}_{22}\right)^{2}+4|\Lambda^{S}_{12}|^{2}|\Lambda^{F}_{12}|^{2}}{\left(\Lambda^{S}_{11}\Lambda^{F}_{11}+\Lambda^{S}_{22}\Lambda^{F}_{22}\right)^{2}}}\,,$
(11)

which is the expression we use in the main text.

Figure A1: Simulated DoP for different spatial and spectral measurement
bandwidths.

In Fig. A1, we show the DoP calculated with Eq. (11) as a function of the
spatial and spectral measurement bandwidths. We use the spatial and spectral
basis functions defined in Eq. (4) of the main text, while considering a
Gaussian spatial profile with approximately $w_{0}=600\,\mu$m beam waist and a
Gaussian spectral profile corresponding to a pulse duration of $\tau=220$ fs,
matching the parameters used in the experiment. The spatial bandwidth
corresponds to the opening angle of an angular mask, while the spectral
bandwidth corresponds to the width of a sharp window centred at $780$ nm. As
expected, the simulated DoP has excellent agreement with the measured DoP
shown in Fig. 4 of the main text.

##  Appendix B Analogy to the quantum description

###  B.1 Bases

Figure B2: Overview of the chosen mutually unbiased bases and assigned labels
in the three different DoF’s. The projections in the spatial and spectral
domains on infinitesimal small slits are depicted by the dotted lines.

Fig. B2 shows an overview of the mutually unbiased bases and the assigned
labels needed to show the non-separability in analogy to the quantum
description. The following list gives additional descriptions to the shown
figure.

  * •

Space domain: The computational basis in the space domain corresponds to a
light field with counter-rotating spiral phases as discussed in the main text.
Here, the basis states are two doughnut-shaped intensity structures with
opposite azimuthal phase gradients as indicated by the green arrows in Fig.
B2. The complementary $x$-bases correspond to a pair of two-lobbed structures
which are rotated by 90∘ with respect to each other. The $y$-basis consists of
two similar modes that are rotated by 45∘ with respect to the $x$-basis [32].
In the measurements, we approximate that a narrow spatial slit projects on a
superposition of orthogonal spatial modes [34]. These slits are depicted by
dotted lines in Fig. B2.

  * •

Wavelength/Time domain: The computational basis for the spectral/temporal DoF
is formed by two time-delayed pulses, which correspond to the functions
$F_{1,2}(\lambda)=F_{0}(\lambda)e^{\pm i\pi\tau c/\lambda}$ given in the main
text. The $x$ and $y$ basis are made up of four specific frequencies at which
the intensity difference between the horizontal/vertical or diagonal/anti-
diagonal polarization components are maximized (shown in Fig. B2 as dotted
lines). We follow the same argumentation as in the spatial domain: The
transmission through a very narrow slit, i.e. strongly filtering in the
frequency domain, approximates the projection on a superposition of modes.

  * •

Polarization: The computational domain is defined by the left and right
circularly polarized components of the light field. The $x$ basis projects on
horizontal/vertical linear polarization, while the $y$ basis projects on
diagonal/anti-diagonal linear polarization.

###  B.2 Normalization

In the measurements, the intensity of the states projected onto polarization
and space in the chosen bases are measured in the spectral domain. In the
Dirac notation all possible measurement outcomes in a basis should add up to a
probability amplitude of 1. Thus, the intensity measurements have to be
normalized. First, the measured intensity signal is normalized with respect to
the intensity summed up over all polarizations, i.e. with respect to the
$S_{0}$ Stokes parameter. Then, the measurements of all possible permutations
of the basis are summed up to find a normalization constant $N$. In the
$xyy$-basis the state for example reads

$\displaystyle\ket{\phi}=$
$\displaystyle\frac{1}{N}\left({\color[rgb]{0,0,1}\alpha}\ket{x,y,-y}+{\color[rgb]{1,0,0}\beta}\ket{x,y,y}+{\color[rgb]{0,0,1}\gamma}\ket{x,-y,y}\right.$
$\displaystyle\left.+{\color[rgb]{1,0,0}\delta}\ket{x,-y,-y}+{\color[rgb]{0,0,1}\epsilon}\ket{-x,y,y}+{\color[rgb]{1,0,0}\zeta}\ket{-x,y,-y}+\right.$
$\displaystyle\left.{\color[rgb]{0,0,1}\eta}\ket{-x,-y,-y}+{\color[rgb]{1,0,0}\theta}\ket{-x,-y,y}\right)$
(12)

where ideally for a non-separable state
${{\color[rgb]{1,0,0}\beta}}={{\color[rgb]{1,0,0}\delta}}={{\color[rgb]{1,0,0}\zeta}}={{\color[rgb]{1,0,0}\theta}}=0$.
The normalization factor is given by

$N=\absolutevalue{{\color[rgb]{0,0,1}\alpha}}^{2}+\absolutevalue{{\color[rgb]{1,0,0}\beta}}^{2}+\absolutevalue{{\color[rgb]{0,0,1}\gamma}}^{2}+\absolutevalue{{\color[rgb]{1,0,0}\delta}}^{2}+\absolutevalue{{\color[rgb]{0,0,1}\epsilon}}^{2}+\absolutevalue{{\color[rgb]{1,0,0}\zeta}}^{2}+\absolutevalue{{\color[rgb]{0,0,1}\eta}}^{2}+\absolutevalue{{\color[rgb]{1,0,0}\theta}}^{2}.$
(13)