# Deploying and Evaluating LLMs to Program Service Mobile Robots

Zichao Hu1, Francesca Lucchetti2, Claire Schlesinger2, Yash Saxena1, Anders
Freeman3,
Sadanand Modak1, Arjun Guha2, Joydeep Biswas1 1Department of Computer Science,
University of Texas at Austin, {zichao, yash.saxena, sadanandm,
<EMAIL_ADDRESS>College of Computer Sciences, Northeastern
University, {lucchetti.f, schlesinger.e<EMAIL_ADDRESS>of Computer Science, Wellesley College<EMAIL_ADDRESS>
###### Abstract

Recent advancements in large language models (LLMs) have spurred interest in
using them for generating robot programs from natural language, with promising
initial results. We investigate the use of LLMs to generate programs for
service mobile robots leveraging mobility, perception, and human interaction
skills, and where _accurate sequencing and ordering_ of actions is crucial for
success. We contribute CodeBotler, an open-source robot-agnostic tool to
program service mobile robots from natural language, and RoboEval, a benchmark
for evaluating LLMs’ capabilities of generating programs to complete service
robot tasks. CodeBotler performs program generation via few-shot prompting of
LLMs with an embedded domain-specific language (eDSL) in Python, and leverages
skill abstractions to deploy generated programs on any general-purpose mobile
robot. RoboEval evaluates the correctness of generated programs by checking
execution traces starting with multiple initial states, and checking whether
the traces satisfy temporal logic properties that encode correctness for each
task. RoboEval also includes multiple prompts per task to test for the
robustness of program generation. We evaluate several popular state-of-the-art
LLMs with the RoboEval benchmark, and perform a thorough analysis of the modes
of failures, resulting in a taxonomy that highlights common pitfalls of LLMs
at generating robot programs. We release our code and benchmark at
https://amrl.cs.utexas.edu/codebotler/.

##  I Introduction

Figure 1: The system diagram of CodeBotler and RoboEval. CodeBotler receives a
user prompt and queries a large language model (LLM) to generate a robot
program (a). It can execute the program to send instructions to a robot via
ROS Actionlib (c). Separately, RoboEval evaluates programs generated from its
benchmark tasks using a symbolic simulator and a temporal trace evaluator to
determine whether each program satisfies the task constraints or not (b).

We are interested in deploying service mobile robots to perform arbitrary user
tasks from natural language descriptions. Recent advancements in large
language models (LLMs) have shown promise in related applications involving
visuomotor tasks [1, 2], planning [3, 4, 5, 6], and in this paper, we
investigate the use of LLMs to generate programs for _service mobile robots_
leveraging mobility, perception, and human interaction skills, where accurate
_sequencing and ordering of actions_ is crucial for success. We contribute
CodeBotler and RoboEval: CodeBotler is an open-source robot-agnostic tool to
generate general-purpose service robot programs from natural language, and
RoboEval is a benchmark for evaluating LLMs’ capabilities of generating
programs to complete service robot tasks.

CodeBotler leverages an embedded domain-specific language (eDSL) in Python to
abstract key robot skills, and includes robot-agnostic bindings for such tasks
using ROS Actions [7]. Given few-shot examples, CodeBotler uses LLMs to
convert natural language task descriptions into programs in the eDSL, which
are then executed on real robots using a lightweight interpreter to interface
with robot-specific skills.

While the capabilities of LLMs at producing robot programs are impressive,
they are still susceptible to a variety of failures. To understand precisely
the failure modes of LLMs in producing robot programs, we introduce the
RoboEval benchmark. Given a program generated by CodeBotler, RoboEval first
executes it in a symbolic simulator to generate multiple program traces from
different initial world states, and then checks these traces against a set of
temporal checks that define correct behavior for the task for each initial
world state.

Existing code completion benchmarks tackle simple data processing functions
[8] or low-level robot skills [1], which are amenable to simple input-output
unit tests. However, code generation for general-purpose service robot
programs cannot be evaluated just on input/output sequences. For example,
given a task _“Check how many conference rooms have no markers”_ , it is
insufficient to just test whether the LLM-generated program executes to state
the correct answer — to ensure correctness, a correct program must first visit
all conference rooms and check for markers there before arriving at the
result. We also observe that there are significant variations in the
correctness of generated programs with small variations in the phrasing of the
natural language task descriptions [9]. We thus contribute three key ideas to
test for both correctness and robustness of LLM-generated robot programs: 1)
we evaluate the _execution traces_ of programs; 2) we check whether the
execution traces satisfy _temporal logic_ properties that encode correctness
for each task; and 3) we _vary the prompts_ and to test for robustness. Fig. 1
shows the system diagram of CodeBotler and RoboEval.

We further categorize the types of failures of different LLMs in the RoboEval
benchmark and find several common categories of failures, including Python
run-time errors, errors in executing infeasible robot actions, and errors in
satisfying task requirements. We analyze the types of errors in each category
and find several common modes of failures across LLMs. We believe this
analysis will be invaluable in furthering research on LLM-guided robot program
generation. Driven by our initial findings, we include a simple rejection
sampling procedure that shows immediate improvements in reducing robot
execution errors of LLM-generated robot programs.

In summary, this paper contributes:

  1. 1.

CodeBotler, an open-source tool to generate robot programs from natural
language using LLMs, and to enable robot-agnostic deployment of such programs;

  2. 2.

RoboEval, a benchmark to evaluate LLM-generated robot programs for service
mobile robots;

  3. 3.

a comprehensive analysis and taxonomy of failures of LLM-generated robot
programs; and

  4. 4.

a rejection sampling mechanism to reduce robot execution failures of LLM-
generated robot programs.

##  II Related Work

_LLMs for Robotics Applications._ Using LLMs to perform robotics tasks [10,
11, 12, 13, 4, 6, 3, 5, 14, 15, 16, 1, 2, 17, 18, 19, 20, 21] has attracted a
lot of attention because of their impressive out-of-box and commonsense
reasoning capabilities [22]. One approach uses LLMs as task planners to break
down a free-form natural language task description into multiple sub-goals.
Language Models as Zero-Shot Planners [4] expresses these sub-goals in the
form of natural language and builds an interpreter to convert these subgoals
into robot actions. LLM+p [6] outputs these sub-goals in the form of the well-
defined planning domain definition language (PDDL). Another approach leverages
the code-writing capabilities of LLMs to generate programs for the robot to
execute. Code-as-Policies [1] defines a robot-centric formulation of language
model-generated programs (LMPs). It proposes a hierarchical method to query an
LLM and generate executable Python programs that invoke parameterized robot
primitives. Voxposer [13] builds on the Code-as-Policies’ LMP formulation and
defines a set of primitives that enables the LLMs to generate Python programs
to create voxel cost maps. Then, it plans on the voxel cost maps and carries
out the specified manipulation tasks. CodeBotler builds on the LMP
formulation, in which we define a set of 8 robot primitives specific to
service mobile robots.

_Evaluating LLMs._ The evaluation of LLMs is an ongoing effort due to the
scope and variability of model generations [23, 24, 25]. Much progress has
been made in evaluating code-writing LLMs [9, 8, 26, 27]. While these
programming tasks assess language understanding, reasoning, algorithms, and
basic mathematics, they do not address the skills of embodied agents. In the
domain of embodied agents, ProgPrompt [2] creates a dataset of 70 household
simulation tasks to evaluate the code-writing capabilities of LLMs.
Simulations are useful for capturing complex agent-environment dynamics but
not essential for verifying program logic, and creating them can be time-
consuming. To address this limitation, we create RoboEval, a lightweight
benchmark that uses a symbolic simulator to evaluate the temporal correctness
of a robot program.

##  III CodeBotler: Robot-Agnostic Deployments

CodeBotler is a lightweight open-source tool that 1) defines robot-agnostic
skills; 2) provides a web interface for users to input commands and generate
language model programs (LMPs) using LLMs; and 3) executes the LMPs on real
robots via ROS Actionlib. In the following sections, we present the design of
CodeBotler.

Figure 2: CodeBotler robot skills (a), RoboEval temporal logic (RTL) formula
(b), and the LTL specifications vs. the RTL specifications of an example task
(c). In section (c), the terms _office_ , _meet_ , and _lunch_ are used to
represent the regex patterns. The RTL specifications are simpler to express
and have improved readability.

###  III-A CodeBotler Language Design

CodeBotler leverages an embedded domain-specific language (eDSL) in Python to
abstract 8 commonly available service robot skills. Fig. 2a shows the
definition of each skill. The design of the CodeBotler eDSL encompasses skills
such as:  1. get_current_locationand get_all_rooms that inspect the robot’s
state and the world configurations, so they do not need to be hard-coded or
manually specified in the prompt, as in previous work [1, 13, 2],  2.
is_in_roomthat utilizes the zero-shot visual-language models (VLMs) for
perceptual reasoning,  3. askthat provides a structured interface for human
interaction through multiple-choice questions, facilitating both the LMP in
processing human responses and interaction with a robot via touch-screen or
audio input.  4. pick, place, and go_to that represent core robot manipulation
and navigation abilities.  These abstractions allow CodeBotler to be used for
robot-agnostic deployments. An LMP can often be reused across different maps
and robots through the user interface.

###  III-B User Interface (UI) and Robot Program Generation

The CodeBotler UI includes a task input pane, a program preview pane, and a
status monitor to track the generation of programs on the robot. When given a
user task, CodeBotler combines it with a prompt prefix containing robot skills
and a few example programs, and then queries an LLM for program generation. In
addition, CodeBotler supports many LLM interfaces, including the OpenAI API
[28], the Google PaLM API [29], and HuggingFace models (AutoModel and Text-
Generation-Inference).

###  III-C Robot Deployment With CodeBotler

CodeBotler is designed to work with the ROS system and acts as a client by
utilizing the ROS Actionlib [7]. When CodeBotler executes an LMP and
encounters a statement that invokes a robot skill (e.g., go_to("lobby")), it
publishes a goal (e.g., "lobby") to the appropriate remote topic (e.g.,
"/go_to_server") for a ROS action server on the robot to pick up. This
approach makes CodeBotler independent of any specific robot platform and
permits CodeBotler to operate both onboard and externally to the robot.
Additionally, it provides robot deployers with the flexibility to customize
the ROS action server to accommodate their specific needs for these
primitives111The code for our implementation of the robot skills can be found
at https://github.com/ut-amrl/codebotler_amrl_impl..

Figure 3: The RoboEval benchmark includes 16 tasks, each with 5 prompt
paraphrases. The figure displays these tasks’ names and a detailed example of
the task CountSavory.

##  IV The RoboEval Benchmark

RoboEval consists of a simulator, an evaluator, and a benchmark suite of
tasks. Given $P$, the space of natural language prompts describing service
mobile tasks, and $\Pi$, the set of possible LMPs, CodeBotler generates LMPs
$\pi\in\Pi$ given a prompt $p\in P$. The symbolic simulator accepts a world
state $w\in W$ and an LMP $\pi\in\Pi$, and produces a program trace $r\in R$.
The evaluator accepts a trace and a temporal constraint $c\in C$, and returns
whether the trace satisfies the constraint or not (SAT/UNSAT).

$\displaystyle\textsc{CodeBotler}{}:P$ $\displaystyle\rightarrow\Pi$
$\displaystyle\mathrm{Simulator}:\Pi\times W$ $\displaystyle\rightarrow R$
$\displaystyle\mathrm{Evaluator}:R\times C$
$\displaystyle\rightarrow\\{\texttt{SAT},\texttt{UNSAT}\\}$

The results derived from traces over multiple world states and multiple task
prompts yield the success rate for an LLM on a particular task. The RoboEval
benchmark thus consists of tasks $T_{i}(i\in[1,N])$, where each task consists
of $M$ prompts, and $K$ world states. Each world state has a corresponding
temporal check. Each RoboEval task thus consists of a tuple of prompts and
multiple world-states to check against a constraint (one constraint per world
state):

$\displaystyle T_{i}=\left\langle\\{p_{i}^{j}|j=1:M\\},\\{\langle
w_{i}^{k},c_{i}^{k}\rangle|k=1:K\\}\right\rangle$

We present next  1. the RoboEval simulator,  2. the RoboEval evaluator, and
3. the tasks in the RoboEval benchmark.

###  IV-A The RoboEval Simulator

For each task $T_{i}$, the RoboEval benchmark includes multiple world states
to check against. Each world state $w_{i}^{k}\in W$ consists of  1. a list of
rooms in the world that GetAllRooms() returns, and which GoTo() is valid for;
2. a list of objects in the world that IsInRoom() returns true for;  3. a list
of objects that can be manipulated using Pick() and Place(); and  4. a list of
responsive humans, their locations, and regular expressions that define their
responses to Ask().  Thus, a single LMP may produce very different traces when
simulated with different initial world states. The simulator consists of a
Python interpreter and a symbolic simulation of each robot skill, and the
result of running an LMP $\pi$ is recorded as a trace $r$ as a sequence of
robot skills that were executed, along with the parameters (e.g., the location
parameter of a GoTo call). All Python errors or robot execution errors are
logged during simulation.

###  IV-B The RoboEval Evaluator

Given a trace $r_{i}^{k}$ produced by simulating an LMP $\pi$ with an initial
world state $w_{i}^{k}$, the RoboEval evaluator checks whether $r_{i}^{k}$
satisfies the temporal check $c_{i}^{k}$ that defines correct execution of the
task for that world state. $c_{i}^{k}$ may consist of multiple conditions,
expressed in conjunctive normal form over multiple temporal constraints. We
review Linear Temporal Logic, which is well-suited to codifying such
constraints in order to check for correctness.

_Linear Temporal Logic._ An LTL formula follows the grammar shown in Fig. 2c —
it composes atomic propositions $\pi\in\Pi$ with logical operators
$\neg,\wedge,\vee$ and temporal operators
$\mathcal{F},\mathcal{G},\mathcal{U},\mathcal{N}$. Given LTL formulas
$\phi_{1},\phi_{2}$ defined over a temporal sequence, $\mathcal{F}\phi_{1}$ is
true iff $\phi_{1}$ is true eventually at some point along the sequence,
$\mathcal{G}\phi_{1}$ is true iff $\phi_{1}$ is true over the entire sequence,
and $\phi_{1}\mathcal{U}\phi_{2}$ is true iff $\phi_{1}$ for a sub-sequence
and $\phi_{2}$ is true for the remainder of the sequence after that.
$\mathcal{N}\phi$ is true for a sequence iff the next element in the sequence
satisfies $\phi$.

_RoboEval Temporal Logic_ While LTL suffices for writing robot task
specifications, these LTL formulas can become complex as task complexity
increases. For example, consider an example task $\mathcal{T}_{1}$ where a
user asks the robot, _“tell Alice in her office to meet me in the lobby if she
agrees to lunch”_. To complete this task, the robot 1) _first_ needs to go to
Alice’s office; 2) _then_ ask Alice whether she would like to have lunch; and
3) _finally_ if she agrees, tell her to meet in the lobby. Fig. 2c shows the
complete LTL specification for this task. Declaring such specifications is
quite tedious and error-prone. To address this challenge, we observe that 1)
specifying temporal logic is easier and less error-prone for specific
scenarios (e.g., one scenario for if Alice says yes, and a different scenario
for no), and 2) the temporal formulas for robot tasks necessarily depend on
the robot skills. We thus introduce the _RoboEval Temporal Language_ (RTL), a
language derived from LTL that is particularly well-suited to specifying
temporal logic formulas for robot tasks. Fig. 2b shows the grammar of RTL, and
Fig. 2c shows the corresponding RTL formula for task $\mathcal{T}_{1}$. An
additional advantage of the condition expressed in RTL vs. LTL is improved
readability.

###  IV-C The RoboEval Benchmark Tasks

The RoboEval benchmark contains a suite of 16 tasks. Fig. 3 shows the names of
the these tasks, along with a detailed example of the task CountSavory 222A
comprehensive list of the task descriptions can be found at
https://amrl.cs.utexas.edu/codebotler/ . These tasks are designed to check
whether an LMP can  1. ground language instructions to correct function calls
to robot primitives;  2. perform accurate sequencing of robot actions;  3.
handle complex control flows based on different world configurations;  4.
solve arithmetic problems;  5. comprehend open-world knowledge.  In addition,
research has shown [30] that LLMs may not be as robust as previously thought,
and trivial prompt variations could cause significant performance variations
for LLMs [31, 9]. For this reason, we provide 5 different paraphrases of the
task prompt to evaluate the robustness of an LLM in dealing with slight prompt
variations.

##  V Benchmark Results And Analysis

GPT-4GPT-3.5PaLM2CodeLLamaStarCoderSayGoodDayGetDrinkStaplerDeliveryMovieMessengerElevatorTourHalloweenListStaplerSupplyCountSavory00.250.50.751average
pass@10.970.991.000.911.001.000.910.840.750.831.000.980.780.960.581.001.000.960.400.911.001.001.000.990.040.990.910.720.800.460.810.800.620.391.000.990.900.980.360.00FindBackpackHalloweenShoppingLunchBreakGrilledCheeseSetTemperatureLunchTimeWeatherPollMailDelivery(a)
Task Performance Breakdown00.250.50.751average
<EMAIL_ADDRESS>pass@10.990.900.990.370.000.990.850.120.210.11LunchTimeSetTemperature(b)
Commonsense Reasoning00.250.50.751average
<EMAIL_ADDRESS>
Figure 4: 5 LLMs are evaluated on the RoboEval benchmark. Each benchmark task
contains 5 different prompt paraphrases and each bar represents the average
pass@1 score of an LLM for generating responses across all 5 prompts within a
given RoboEval benchmark task. Each error bar indicates the range from the
highest to the lowest pass@1 score across all prompts. On the left side (a),
the performance of LLMs on each task of the RoboEval benchmark is displayed.
On the right side (b), the chart shows the performance of LLMs on tasks that
have been adapted to exclusively evaluate the models’ proficiency in
performing commonsense reasoning.

1.00.80.60.40.20.001pass@1$\displaystyle{0}$$\displaystyle{20}$$\displaystyle{40}$$\displaystyle{60}$$\displaystyle{80}$$\displaystyle{100}$Fraction
of Prompts
(%)67.5098.7543.7577.5028.7560.0026.2563.7520.0051.25GPT-4GPT-3.5PaLM2CodeLLamaStarCoder

Figure 5: Cumulative Distribution Function (CDF) curves depict the percentage
of prompts for which each LLM can generate correct LMPs at various pass@1
score thresholds. A perfect LLM would show a horizontal line at 100%,
indicating it can generate correct LMPs for all prompts with a pass@1 score of
1. To maintain visual clarity, we limit the x-axis to $10^{-3}$ since all CDF
plots eventually reach 100%.

(a) Overall Error
BreakdownStarCoderCodeLlamaPaLM2GPT-3.5GPT-431.2%3.7%29.3%35.8%45.7%3.4%30.3%20.6%44.8%9.7%22.0%23.5%66.0%11.9%4.1%18.0%85.5%11.7%Successful
ProgramPython ErrorRobot Execution ErrorTask Completion Error(b) Python Error
BreakdownStarCoderCodeLlamaPaLM2GPT-3.5GPT-46.7%10.1%16.1%65.8%35.0%18.6%45.7%22.0%3.4%73.4%3.3%8.8%84.3%100.0%ValueErrorUnboundLocalErrorIndexErrorNameErrorZeroDivisionErrorSyntaxErrorKeyErrorTypeErrorTimeoutError(c)
Robot Execution Error
BreakdownStarCoderCodeLlamaPaLM2GPT-3.5GPT-45.7%59.7%6.6%20.8%7.2%67.3%6.7%5.9%18.1%46.0%9.3%21.7%22.6%10.2%64.5%23.5%92.9%3.6%3.6%PickInvalidObjectAskEmptyOptionsGoToInvalidLocationAskNoPersonPlaceNoObjectPickWhileHolding(d)
Task Completion Error
BreakdownStarCoderCodeLlamaPaLM2GPT-3.5GPT-49.1%33.5%9.2%42.5%7.7%28.7%5.4%4.3%51.2%3.2%12.6%23.8%6.2%5.9%43.1%4.6%22.1%73.6%31.7%68.2%Manipulation
at LocationEvent OrderingInitial/TerminalLocationExhaustive SearchAsk at
LocationSay at LocationCheck Entity at Location

Figure 6: Causes of failures for LMPs on the RoboEval benchmark.

To gain insights into the capabilities and limitations of different state-of-
the-art LLMs for generating service mobile robot LMPs, we use the RoboEval
benchmark to answer the following questions:

  1. 1.

First, how do different LLMs perform in generating programs for tasks in the
RoboEval benchmark?

  2. 2.

Second, when a generated service robot LMP fails, what are the causes?

To investigate these two questions, we evaluate five LLMs:  1. GPT-4 [32],  2.
GPT-3.5 [33] (text-davinci-003), and  3. PaLM2 [34] (text-bison-001) as state-
of-the-art API-only proprietary models; and  4. CodeLlama [35] (Python-34b-hf)
and  5. StarCoder [36] as open-access models.  When evaluating code generation
models, we use standard values of $T=0.2$ for temperature and $p=0.95$ for
nucleus sampling. In the following subsections, we discuss our analysis of
each question in detail.

###  V-A Performance Of LLMs On The RoboEval Benchmark

The RoboEval benchmark consists of 16 tasks, each with 5 prompt paraphrases,
totaling 80 different prompts. For each prompt, we generate 50 program
completions and calculate the pass@1 score [8], a common metric for LMP
evaluation. This score indicates the probability of an LMP being correct if an
LLM generates only one LMP for a given prompt.

_Overall Performance Analysis._ We first investigate the overall performance
of each LLM in generating LMPs. Since each LLM gets a pass@1 score for every
prompt, we compute the percentage of prompts that have a pass@1 score greater
than or equal to a threshold value, which ranges from $1$ to $0$. We present
this information in Fig. 5 as a Cumulative Distribution Function (CDF).
Although relaxing the pass@1 score threshold for each LLM increases prompt
coverage, there are still certain prompts (ranging from 48.75% for StarCoder
to 1.25% for GPT-4) where LLMs consistently fail to generate correct LMPs.

_Performance of LLMs on Individual Tasks._ We then look into the performance
of each LLM for specific tasks. We present this information in Fig. 4a. Since
each task has 5 different prompt paraphrases and a pass@1 score is calculated
for each prompt, we average the pass@1 score across 5 prompts. Additionally,
we plot an error bar to visualize the variation between the highest to the
lowest pass@1 score over all prompts for each task.

From the error bars, we note a considerable disparity between a model’s best
and worst pass@1 scores for a given task. This suggests that models are not
robust to changes in the phrasing of the prompt. Single-word changes can
result in substantial performance variations. A common example is changing the
verb “ask” to “inquire”; “Ask him about his available ingredients” in the
GrilledCheese task thus becomes “inquire about his available ingredients”.
This seems to affect some code model’s ability to call the robot ask function.

From this chart, we further note a high variation in performance across tasks.
The top row contains the high-performing tasks, where four or more models
score over 0.7 on pass@1, while the bottom row contains the low-performing
tasks where three or more models score below 0.2 pass@1. To identify why such
a high variation exists, we run an ablation experiment. We notice that tasks
involving commonsense reasoning (CountSavory, GrilledCheese, LunchTime,
SetTemperature) tend to underperform. Hence, we ablate all but commonsense RTL
checks on these tasks. For example, the full CountSavory checks require that
the robot navigates to every office as well as understanding that _beef jerky_
and _ham sandwich_ are savory options; for the ablation, we remove these
navigation checks. We plot the resulting average pass@1 scores of this
experiment on Fig. 4b. If models are failing commonsense reasoning, we expect
the performance to be unchanged after ablation. However, we notice that models
are largely improving in performance. This suggests that the problem lies
elsewhere, which we will analyze next.

###  V-B Causes of Failures of LMPs

Given that LLMs still have room for improvement on the RoboEval benchmark, we
want to understand the causes of failures for LLMs to generate robot programs.
We classify these failures into three categories:  1. Python Errors, including
syntax, runtime, and timeout errors;  2. Robot Execution Errors, that occurs
when a program attempts to execute an infeasible action, such as navigating to
a non-existent (hallucinated) location; and  3. Task Completion Errors, where
the program runs correctly in the simulator but fails RTL checks for task
completion.  We use RoboEval’s symbolic simulator to detect and classify
Python Errors and Robot Execution Errors, and we use RoboEval’s evaluator to
capture the Task Completion Errors. Fig. 6a shows the breakdown of these
failure categories for each LLM. We observe that despite having fewer
parameters, the CodeLLMs (CodeLlama and StarCoder) generally make fewer Python
errors. This suggests that LLMs trained on a larger proportion of code may be
more adept at generating successful completions in the DSL defined in the
prompt.

In the following subsections, we will analyze each error in detail.

_Python Error Analysis._ Fig. 6b shows the specific error breakdown of the
Python Error. The Python Error encompasses a multitude of errors, but their
distribution exhibits a long-tailed pattern. In GPT-3.5 and GPT-4, _NameError_
is predominant because of undefined variables. PaLM2 and CodeLLama, on the
other hand, often generate _TypeError_ due to the misuse of data types, while
Starcoder commonly encounters _TimeoutError_ when the LMPs get stuck in loops.

_Robot Execution Error._ Fig. 6c shows the error breakdown of the Robot
Execution Error. There are 6 root causes of robot execution errors:

  1. 1.

GoToInvalidLocation/PickInvalidObject: the program calls go_to or pick with a
hallucinated argument;

  2. 2.

PlaceNoObject/PickWhileHolding: the program tries to pick/place an object when
it is not/already holding one;

  3. 3.

AskNoPerson: the program calls ask at a location with no person nearby; and

  4. 4.

AskEmptyOptions: the program calls ask with an empty list of options for the
person to choose from.

We first observe that hallucination plays a substantial role in causing
errors, Specifically, the GoToInvalidLocation and PickInvalidObject errors
contribute to 33.7% of the total robot execution errors in GPT-3.5, 67.7% in
PaLM2, 73.2% in CodeLlama, and 80.5% in StarCoder. We also notice that another
important source of errors arises from the PlaceNoObject and PickWhileHolding
errors, as well as AskNoPerson errors for GPT-3.5. These errors require the
program to be aware of the internal state of the robot or keep track of the
external world state. The prevalence of these errors suggests a gap in LLM’s
ability to keep track of changing states by binding agents to current states.
Some of these errors also result from a failure to respect the behavior of our
robot functions. For example, our pick function only allows the robot to pick
an object if it is not already holding one. However, the model may not include
this information in its completion and attempt to pick up multiple objects.

$\displaystyle{65}$$\displaystyle{70}$$\displaystyle{75}$$\displaystyle{80}$$\displaystyle{85}$$\displaystyle{90}$$\displaystyle{95}$$\displaystyle{100}$1.00.80.60.40.20.001$\displaystyle{0}$GPT-4
pass@1Fraction of Prompts
(%)$\displaystyle{100}$$\displaystyle{40}$$\displaystyle{45}$$\displaystyle{50}$$\displaystyle{55}$$\displaystyle{60}$$\displaystyle{65}$$\displaystyle{70}$$\displaystyle{75}$1.00.80.60.40.20.001$\displaystyle{0}$GPT-3.5
pass@1Fraction of Prompts
(%)$\displaystyle{100}$$\displaystyle{25}$$\displaystyle{30}$$\displaystyle{35}$$\displaystyle{40}$$\displaystyle{45}$$\displaystyle{50}$$\displaystyle{55}$$\displaystyle{60}$1.00.80.60.40.20.001$\displaystyle{0}$PaLM2
pass@1Fraction of Prompts
(%)$\displaystyle{100}$$\displaystyle{25}$$\displaystyle{30}$$\displaystyle{35}$$\displaystyle{40}$$\displaystyle{45}$$\displaystyle{50}$$\displaystyle{55}$$\displaystyle{60}$$\displaystyle{65}$1.00.80.60.40.20.001$\displaystyle{0}$CodeLLama
pass@1Fraction of Prompts
(%)$\displaystyle{100}$$\displaystyle{20}$$\displaystyle{25}$$\displaystyle{30}$$\displaystyle{35}$$\displaystyle{40}$$\displaystyle{45}$$\displaystyle{50}$1.00.80.60.40.20.001$\displaystyle{0}$StarCoder
pass@1Fraction of Prompts (%)Max Retry Limit: 100Max Retry Limit: 8Max Retry
Limit: 4Max Retry Limit: 2Max Retry Limit: 0

Figure 7: Cumulative Distribution Function (CDF) of the LLMs’
performance across different max retry limits. As the max retry
limit increases, all five LLMs improve in performance.

$\displaystyle{65}$$\displaystyle{70}$$\displaystyle{75}$$\displaystyle{80}$$\displaystyle{85}$$\displaystyle{90}$$\displaystyle{95}$$\displaystyle{100}$$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$$\displaystyle{0}$$\displaystyle{95}$$\displaystyle{100}$Maximum
Retry LimitFraction of Program Completions
(%)GPT-4GPT-3.5PaLM2CodeLLamaStarCoder

Figure 8: Cumulative Distribution Function (CDF) of the fraction of program
completions that can be executed over different max retry limits.

_Task Completion Error._ Fig. 6d shows the breakdown of the Task Completion
Errors. We classify every temporal check in the RoboEval benchmark into one of
the following categories:

  1. 1.

{Say/Ask/Manipulation/CheckEntity} AtLocation: The task requires executing a
specific action (say, ask, pick, place or check) at a specific location, but
the program fails to do so;

  2. 2.

Initial/Terminal: The program does not accurately perform an initial or final
action;

  3. 3.

EventOrdering: The program does not carry out actions in the prescribed
sequence or has redundant navigation;

  4. 4.

Location: The program commands the robot to visit a location irrelevant to the
task; and,

  5. 5.

ExhaustiveSearch: The program does not visit all locations required to
complete the task.

We notice that the most common task completion errors in each LLM are the
SayAtLocation error and the Initial/Terminal error. These two types of errors
constitute the majority of task completion errors, making up $100\%$ of the
errors in GPT-4, $95.7\%$ in GPT-3.5, $66.9\%$ in PaLM2, $79.9\%$ in
CodeLlama, and $76.0\%$ in StarCoder.

Figure 9: Examples of successful (top) and failed (bottom) program completions
related to the _“come back”_ instruction.

After analyzing the programs that produce these errors, we find that most are
related to a recurring statement in many of our tasks: _“Come back and tell
me…”_ (Fig. 9 shows an example). This statement is often used in tasks in
which we want a service robot to go to a different location, complete a task,
and then return to us with an update on its progress. This is the most common
instruction that fails across models, and also the error that is directly
responsible for many of our low-performing tasks. Programs that fail the
_“come back”_ instruction often do not use the get_current_location primitive
to record the robot’s starting location, and as a result, cannot refer to the
starting location at the end of the program. This also causes LLMs to
hallucinate locations and produce the GoToInvalidLocation errors — rather than
referencing a non-existent start_loc variable, some LLMs generate a location
based on context clues and send the robot to that location (for example,
returning to the kitchen in GrilledCheese, although no kitchen is mentioned in
the prompt). We observe that in tasks without a _“come back”_ instruction, but
an explicit return statement like _“return to the mail room”_ do not suffer
from the same error.

##  VI Improving Robot Program Generations

Based on the previous analysis, we aim to improve the pass@1 score for LLMs in
producing service robot programs. Recognizing the breadth of potential
improvements, this study focuses on an initial step towards mitigating one
primary issue: program execution errors that arise within the simulator,
specifically Python Errors and the majority of Robot Execution Errors.

Since the RoboEval symbolic simulator is capable of detecting these program
failures, we can automatically reject such programs before deploying them on a
robot. As a result, we can adopt a simple rejection sampling strategy to keep
prompting an LLM for a new program until the program eventually can execute or
a maximum retry limit is reached. In reality, we do not want to tune the world
state inside the simulator for every specific task, so our goal is to create a
generic world state that is task-agnostic yet capable of catching the Python
Errors and most of the Robot Execution Errors. In this world state, we cannot
know where objects and humans are located and how humans respond beforehand.
Hence, we only specify a list of rooms the robot can navigate to in the world
state, and we let the world state randomly return a result when the program
invokes the is_in_room and ask skill.

We apply this strategy to all five LLMs with four different maximum retry
limits ($2,4,8,100$) and compare them with the baseline (without rejection
sampling). Fig. 8 shows the CDF curves of the percentage of prompts that can
be successfully generated given a threshold of the pass@1 score for each LLM
with respect to different maximum retry limits. We observe an improvement in
performance across all LLMs as the maximum retry limit is increased.

We then investigate how effective the rejection sampling strategy is in
eliminating the program execution errors. We compute the percentage of total
program completions that can be eventually executed over different maximum
retry limits and plot it as a CDF in Fig. 8. Interestingly, we observe that
some LLMs (PaLM2, CodeLlama, and Starcoder) never generate successful
completions for some tasks. As a result, while the rejection sampling strategy
can improve the performances of LLMs, it is not enough to resolve all program
execution errors.

This observation, coupled with prior findings of LLMs’ consistent failures in
generating correct LMPs as detailed in Section V-A, points to a systemic
challenge in LMP generation. It highlights the issue of abstraction matching
[37], which entails aligning ambiguous natural language expressions of user
intent with their precise, unambiguous code representations. In future work,
we would like to explore more sophisticated strategies, such as utilizing
grounded abstraction matching, to solve this problem.

##  VII Conclusion, Limitations, & Future Works

In this work, we present CodeBotler and RoboEval. CodeBotler is an open-source
robot-agnostic tool to generate general-purpose service robot programs from
natural language, and RoboEval is a benchmark for evaluating LLMs’
capabilities of generating programs to complete service robot tasks. We
evaluate the performance of five LLMs in generating robot programs and perform
an analysis of the causes of failures. Our analysis reveals that the errors
exhibit a long-tail distribution, with LLMs predominantly struggling with
hallucination issues and grounding the phrase _“come back”_. Finally, we
propose a rejection sampling strategy to handle program failures. This method
has led to improved performance for all five LLMs.

This work has several limitations that could be addressed in future research.
Firstly, CodeBotler does not support tasks that require fine-grained controls,
such as _“follow Alice to her office”_. Secondly, this study does not consider
the strategies for crafting prompts that could improve the performance of LLMs
in generating service robot programs. Thirdly, the creation of RTL
constraints, although intended to lessen the workload involved in writing
specifications compared to LTL constraints, still requires manual development
for particular tasks. This limitation restricts their practical application in
real-world deployments. Therefore, investigating the Tree-of-Thoughts concept
[38] to dynamically generate RTL checks with LLMs might be valuable, as this
could foresee and prevent program failures before robot deployment.

## References

  * [1] J. Liang, W. Huang, F. Xia, P. Xu, K. Hausman, B. Ichter, P. Florence, and A. Zeng, “Code as Policies: Language Model Programs for Embodied Control,” in _arXiv:2209.07753_ , 2022.
  * [2] I. Singh, V. Blukis, A. Mousavian, A. Goyal, D. Xu, J. Tremblay, D. Fox, J. Thomason, and A. Garg, “ProgPrompt: Generating Situated Robot Task Plans using Large Language Models,” in _ICRA 2023_ , 2023.
  * [3] M. Ahn, A. Brohan, N. Brown, _et al._ , “Do As I Can, Not As I Say: Grounding Language in Robotic Affordances,” 2022.
  * [4] W. Huang, P. Abbeel, D. Pathak, and I. Mordatch, “Language Models as Zero-Shot Planners: Extracting Actionable Knowledge for Embodied Agents,” 2022\.
  * [5] D. Driess, F. Xia, M. S. M. Sajjadi, _et al._ , “PaLM-E: An Embodied Multimodal Language Model,” 2023.
  * [6] B. Liu, Y. Jiang, X. Zhang, Q. Liu, S. Zhang, J. Biswas, and P. Stone, “LLM+P: Empowering Large Language Models with Optimal Planning Proficiency,” 2023.
  * [7] “ROS Actionlib.” [Online]. Available: http://wiki.ros.org/actionlib
  * [8] M. Chen, J. Tworek, H. Jun, Q. Yuan, _et al._ , “Evaluating large language models trained on code,” _arXiv:2107.03374_ , 2021.
  * [9] H. M. Babe, S. Nguyen, Y. Zi, A. Guha, M. Q. Feldman, and C. J. Anderson, “StudentEval: A benchmark of student-written prompts for large language models of code,” _arXiv:2306.04556_ , 2023.
  * [10] C. Huang, O. Mees, A. Zeng, and W. Burgard, “Visual Language Maps for Robot Navigation,” 2023.
  * [11] Y. Ding, X. Zhang, C. Paxton, and S. Zhang, “Task and Motion Planning with Large Language Models for Object Rearrangement,” 2023.
  * [12] H. Wang, G. Gonzalez-Pumariega, Y. Sharma, and S. Choudhury, “Demo2Code: From Summarizing Demonstrations to Synthesizing Code via Extended Chain-of-Thought,” 2023.
  * [13] W. Huang, C. Wang, R. Zhang, Y. Li, J. Wu, and L. Fei-Fei, “VoxPoser: Composable 3D Value Maps for Robotic Manipulation with Language Models,” 2023\.
  * [14] A. Brohan, N. Brown, J. Carbajal, _et al._ , “RT-2: Vision-Language-Action Models Transfer Web Knowledge to Robotic Control,” 2023\.
  * [15] Y. Tang, W. Yu, J. Tan, H. Zen, A. Faust, and T. Harada, “SayTap: Language to Quadrupedal Locomotion,” 2023.
  * [16] W. Yu, N. Gileadi, C. Fu, _et al._ , “Language to Rewards for Robotic Skill Synthesis,” 2023.
  * [17] D. Shah, B. Osinski, B. Ichter, and S. Levine, “LM-Nav: Robotic Navigation with Large Pre-Trained Models of Language, Vision, and Action,” 2022.
  * [18] W. Huang, F. Xia, T. Xiao, _et al._ , “Inner Monologue: Embodied Reasoning through Planning with Language Models,” 2022.
  * [19] J. Wu, R. Antonova, A. Kan, M. Lepert, A. Zeng, S. Song, J. Bohg, S. Rusinkiewicz, and T. Funkhouser, “TidyBot: Personalized Robot Assistance with Large Language Models,” _arXiv preprint arXiv:2305.05658_ , 2023.
  * [20] B. Chen, F. Xia, B. Ichter, K. Rao, K. Gopalakrishnan, M. S. Ryoo, A. Stone, and D. Kappler, “Open-vocabulary Queryable Scene Representations for Real World Planning,” 2022.
  * [21] C. H. Song, J. Wu, C. Washington, B. M. Sadler, W.-L. Chao, and Y. Su, “LLM-Planner: Few-Shot Grounded Planning for Embodied Agents with Large Language Models,” 2023.
  * [22] T. Brown, B. Mann, N. Ryder, _et al._ , “Language Models are Few-Shot Learners,” in _Advances in Neural Information Processing Systems_ , H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, Eds., vol. 33. Curran Associates, Inc., 2020, pp. 1877–1901. [Online]. Available: https://proceedings.neurips.cc/paper˙files/paper/2020/file/1457c0d6bfcb4967418bfb8ac142f64a-Paper.pdf
  * [23] Y. Chang, X. Wang, J. Wang, Y. Wu, L. Yang, K. Zhu, H. Chen, X. Yi, C. Wang, Y. Wang, W. Ye, Y. Zhang, Y. Chang, P. S. Yu, Q. Yang, and X. Xie, “A Survey on Evaluation of Large Language Models,” 2023.
  * [24] P. Liang, R. Bommasani, T. Lee, _et al._ , “Holistic Evaluation of Language Models,” 2023.
  * [25] X. Liu, H. Yu, H. Zhang, _et al._ , “AgentBench: Evaluating LLMs as Agents,” 2023.
  * [26] Y. Lai, C. Li, Y. Wang, T. Zhang, R. Zhong, L. Zettlemoyer, S. W. tau Yih, D. Fried, S. Wang, and T. Yu, “DS-1000: A Natural and Reliable Benchmark for Data Science Code Generation,” 2022.
  * [27] J. Austin, A. Odena, M. Nye, M. Bosma, H. Michalewski, D. Dohan, E. Jiang, C. Cai, M. Terry, Q. Le, and C. Sutton, “Program Synthesis with Large Language Models,” 2021.
  * [28] “OpenAI Platform.” [Online]. Available: https://platform.openai.com/docs/guides/gpt
  * [29] “Generative AI for Developers.” [Online]. Available: https://developers.generativeai.google/
  * [30] A. V. Miceli-Barone, F. Barez, I. Konstas, and S. B. Cohen, “The Larger They Are, the Harder They Fail: Language Models do not Recognize Identifier Swaps in Python,” 2023.
  * [31] K. D. Dhole, V. Gangal, S. Gehrmann, _et al._ , “NL-Augmenter: A Framework for Task-Sensitive Natural Language Augmentation,” 2022.
  * [32] OpenAI, “GPT-4 Technical Report,” 2023.
  * [33] T. B. Brown, B. Mann, N. Ryder, _et al._ , “Language Models are Few-Shot Learners,” 2020.
  * [34] R. Anil, A. M. Dai, O. Firat, _et al._ , “PaLM 2 Technical Report,” 2023\.
  * [35] B. Rozière, J. Gehring, F. Gloeckle, _et al._ , “Code Llama: Open Foundation Models for Code,” 2023.
  * [36] R. Li, L. B. Allal, Y. Zi, _et al._ , “StarCoder: may the source be with you!” 2023.
  * [37] M. X. Liu, A. Sarkar, C. Negreanu, B. Zorn, J. Williams, N. Toronto, and A. D. Gordon, ““what it wants me to say”: Bridging the abstraction gap between end-user programmers and code-generating large language models,” in _Proceedings of the 2023 CHI Conference on Human Factors in Computing Systems_. ACM, apr 2023. [Online]. Available: https://doi.org/10.1145%2F3544548.3580817
  * [38] S. Yao, D. Yu, J. Zhao, I. Shafran, T. L. Griffiths, Y. Cao, and K. Narasimhan, “Tree of Thoughts: Deliberate Problem Solving with Large Language Models,” 2023.


# Violation of Bell’s inequality for Mathieu-Gauss vector modes

Edgar Medina-Segura1, Francisco I. Mecillas-Hernández1, Thomas Konrad2,
Carmelo Rosales-Guzmán1,† and Benjamin Perez-Garcia3,∗ 1 Centro de
Investigaciones en Óptica, A.C., Loma del Bosque 115, Colonia Lomas del
campestre, C.P. 37150 León, Guanajuato, Mexico.  2 School of Chemistry and
Physics, University of Kwazulu-Natal, Private Bag X54001, Durban 4000, South
Africa.  3 Photonics and Mathematical Optics Group, Tecnologico de Monterrey,
Monterrey 64849, Mexico<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Vector beams display varying polarisation over planes transversal to their
direction of propagation. The variation of polarisation implies that the
electric field cannot be expressed as a product of a spatial mode and its
polarisation. This non-separability has been analysed for particular vector
beams in terms of non–quantum entanglement between the spatial and the
polarisation-degrees of freedom, and equivalently, with respect to the degree
of polarisation of light. Here we demonstrate theoretically and experimentally
that Mathieu-Gauss vector modes violate a Bell-like inequality known as the
Clauser-Horn-Shimony-Holt-Bell (CHSH-Bell) inequality. This demonstration
provides new insights into that fact that a more general class of vector modes
with elliptical symmetry also violate Bell inequalities.

Keywords: Non-diffraction beams, Mathieu modes, complex vector beams, Bell
inequalities

##  1 Introduction

Entanglement constitutes one of the cornerstones of quantum mechanics that
gives rise to exceptional properties. For example, entangled particles seem to
be able to communicate instanteneously with each other, which would violate
the principles of relativity. This paradoxical behaviour, dismissively
referred as “spooky actions at a distance” by Einstein [1], is generated by
the superposition of different combinations of states of two or more systems,
which according to classical physics could only occur if the choice of
observable measured on one system would influence the outcome of measurements
on the other system(s), regardless of their separation – a property known as
nonlocality. The concepts of entanglement and nonlocality are closely related,
challenging our understanding of physical reality. In 1964, Bell proposed a
test for all theories that are local (i.e. no instantaneous communication
between remote systems) and realistic (systems have detected properties before
measurement) by deriving an inequality for correlations based on these
assumptions [2]. It has then become a standard procedure to test the
“quantumness” of a system through the violation of Bell’s inequality [3],
which for optical systems usually takes the form of Clauser-Horne-Shimony-
Holt-Bell (CHSH-Bell) inequality [4]. However, entanglement occurs also
between different degrees of freedom of the same quantum system if its state
consists of a superposition of different combinations of states in the two
degrees of freedom, which cannot be factorised, i.e. written as a product of
the state of the first degree of freedom and the state of the second degree of
freedom. For example an electromagnetic wave with elementary excitation (a
single photon) can be written as a vector mode of the form ${\bf u}({\bf
r})=\cos(\alpha/2)u^{+}({\bf r})\hat{e}_{R}+\sin(\alpha/2)u^{-}({\bf
r})\hat{e}_{L}$, which is a superposition of a certain mode $u^{+}({\bf r})$
with right-circular ${\bf\hat{e}}_{R}$ polarisation and a different mode,
$u^{+}({\bf r})$, with left-circular polarisation ${\bf\hat{e}}_{L}$. Such a
vector mode differs from a scalar light field, which can be represented by a
single product of a spatial mode and the polarisation state $u({\bf
r}){\bf\hat{e}}$ and corresponds to an entangled state that violates a Bell
inequality, as opposed to a single photon of a scalar light field with
separable spatial mode and polarisation.

A major contribution to the understanding of entanglement was the realisation
that it occurs in similar form in classical field theories, such as
electromagnetism, and there leads to features formerly thought to be reserved
to quantum systems. For example, vector beams in classical optics
(corresponding in quantum optics to coherent states with classical properties)
can also violate Bell inequalities [5, 6, 7, 8, 9, 10, 11], or they can be
used to show that not all optical operations (Mueller matrices) can be
realised in the lab, but only those that are completely positive [12].

The entanglement between different degrees of freedom of the same quantum
system or electromagnetic wave can of course not be used to refute local
realistic theories via the violation of a Bell inequality, - optics is a local
realistic theory after all. The reason for this is, that the measurements on
different degrees of freedom of the same system cannot be spatially separated
and, therefore, do not allow to exclude the influence of the measurement of
one variable on the outcomes the measurement of another. On the other hand,
strong violation of a Bell inequality by classical light can indicate strong
entanglement of classical degrees of freedom, and vice versa, i.e.
entanglement measures such as Concurrence can be used to quantify the
“vectorness” of a light beam [13, 14].

In this manuscript we study the CHSH-Bell inequalities for Helical-Mathieu
Gauss Vector modes (HMGVM), which are characterised by their elliptical
spatial profile controlled by the ellipticity parameter $\epsilon\in[0,1]$,
which changes the shape of the mode, from circular ($\epsilon=0$) to
elliptical ($0<\epsilon<1$). Here, we derive a general expression for the
CHSH-Bell parameter S that depends on $\epsilon$, which for $\epsilon=0$
reduces to the well-know expression derived for Laguerre-Gaussian modes [13].
The manuscript is organised as follows, first in section 2.1 we briefly
describe the HMGVM, as first introduced in [15], but in this case we use
Dirac’s notation to highlight their similarity with entangled states.
Afterwards, in section 2.2 we derive a theoretical expression for the CHSH-
Bell parameter. In section 3 we describe the experimental techniques behind
the generation of HMGVM as well as the technique to measure CHSH-Bell
inequalities. Finally, we conclude this manuscript with a discussion section.

##  2 Mathematical framework

###  2.1 Helical Mathieu-Gauss vector beams

Helical Mathieu-Gauss Vector (HMGV) beams are generated as a non-separable
superposition of the polarisation degree of freedom and the Helical Mathieu-
Gauss (HMG) beams, which are expressed as [15]

${\bf
u}_{m}(\mathbf{r};q)=\cos\left(\frac{\alpha}{2}\right)u_{m}^{+}(\mathbf{r};q)\hat{\bf
e}_{R}+\sin\left(\frac{\alpha}{2}\right)\exp(i\beta)u_{m}^{-}(\mathbf{r};q)\>\hat{\bf
e}_{L},$ (1)

where $\mathbf{r}=(r,\varphi,z)=(\xi,\eta,z)$ defines the polar– and
elliptical–cylindrical coordinates, respectively; $\alpha\in[0,\pi/2]$ is a
weighting coefficient that determines the contribution of each polarisation
component and the exponential term $\exp(i\beta)$ ($\beta\in[0,\pi]$) is an
inter-modal phase between both polarisation components. The terms $\hat{\bf
e}_{R}$ and $\hat{\bf e}_{L}$ represent the right and left handed unitary
vectors of the circular polarisation basis, respectively. The basis in the
spatial degree of freedom is defined through the modes
$u_{m}^{\pm}(\mathbf{r};q)=\text{MG}^{e}_{m}(\mathbf{r};q)\pm
i\text{MG}^{o}_{m}(\mathbf{r};q)$ [16] with

$\displaystyle\text{MG}_{m}^{e}(\mathbf{r};q)$
$\displaystyle=\exp\left({-\frac{ik_{t}^{2}}{2k}\frac{z}{\mu}}\right)\text{GB}({\bf
r})\text{Je}_{m}(\tilde{\xi},q)\text{ce}_{m}(\tilde{\eta},q),$ (2)
$\displaystyle\text{MG}_{m}^{o}(\mathbf{r};q)$
$\displaystyle=\exp\left({-\frac{ik_{t}^{2}}{2k}\frac{z}{\mu}}\right)\text{GB}({\bf
r})\text{Jo}_{m}(\tilde{\xi},q)\text{se}_{m}(\tilde{\eta},q),$ (3)

where $z$ is the propagation distance, and the complex elliptic coordinates
$(\tilde{\xi},\tilde{\eta})$ are determined by

$\displaystyle x$
$\displaystyle=f_{0}(1+iz/z_{R})\cosh{\tilde{\xi}}\cos{\tilde{\eta}},$ (4)
$\displaystyle y$
$\displaystyle=f_{0}(1+iz/z_{R})\sinh{\tilde{\xi}}\sin{\tilde{\eta}}.$ (5)

The parameter $f_{0}$ is the semifocal separation of the elliptical
coordinates at $z=0$, which relates to the major and minor axis $a$ and $b$,
respectively, as $f^{2}=a^{2}-b^{2}=a\varepsilon$ and $q=(fk_{t}/2)^{2}$ is a
dimensionless parameter. In addition, the function $\text{GB}({\bf r})$ is the
fundamental Gaussian beam defined as

$\text{GB}({\bf
r})=\exp\left({-\frac{r^{2}}{\mu\omega_{0}^{2}}}\right)\frac{\exp(ikz)}{\mu},\\\
$ (6)

where $\mu=\mu(z)=1+iz/z_{R}$, $z_{R}=k\omega_{0}^{2}/2$ is the Rayleigh range
of a Gaussian beam with waist radius $\omega_{0}$. The parameters $k_{z}$ and
$k_{t}$ are the longitudinal and transverse components of the wave vector $\bf
k$, whose magnitude $k=2\pi/\lambda$ satisfies the relation
$k^{2}=k_{t}^{2}+k_{z}^{2}$. Further, the functions $\text{Je}_{m}$ and
$\text{Jo}_{m}$, are the $m$th-order even and odd radial Mathieu functions,
respectively, and $\text{ce}_{m}$ and $\text{se}_{m}$ are the even and odd
$m$th-order angular Mathieu functions [17]. By way of example, Fig. 1
schematically illustrates the intensity and phase profiles of the two
orthogonal scalar Mathieu-Gauss beams, which are recombined in a nonseparable
fashion to produce its vector counterpart. Here, the HMG modes
$u_{m}^{\pm}(\mathbf{r};q)$, whose transverse intensity and phase profile are
shown in Figs. 1(a), 1(d) and 1(b) and 1(e), respectively, are coaxialy
superposed in a non-separable fashion to generate the vector mode shown in
Fig. 1 (c). The phase of the complex Stokes field $S=S_{1}+iS_{2}$ is shown in
Fig. 1 (f). The specific parameters for this example are $k_{t}=4$ mm-1 ,
$e=0.9$, $a=1.5$ and $m=3$.

Figure 1: Theoretical intensity (first row) and phase (second row)
distribution of the circular polarisation components and Stokes field
($S=S_{1}+iS_{2}$), for $\mathbf{u}_{3}(\mathbf{r};q)$. (a) Intensity and (d)
phase of the right–handed circular polarisation component. (b) Intensity and
(e) phase of the left–handed circular polarisation component. (c) Depicts the
intensity of the Stokes field overlapped with the non–homogeneous distribution
of polarisation states (white lines). (d) Shows the phase of the Stokes field.

Importantly, a careful use of Dirac’s notation, allows us to express Eq. 1 as

$\displaystyle\left|\psi_{m}\right\rangle_{q}=\cos\left(\frac{\alpha}{2}\right)\left|m_{q}^{+}\right\rangle\left|R\right\rangle+\sin\left(\frac{\alpha}{2}\right)\exp(i\beta)\left|m_{q}^{-}\right\rangle\left|L\right\rangle,$
(7)

where the kets
$\left|m_{q}^{\pm}\right\rangle=(\left|m_{q}^{e}\right\rangle\pm
i\left|m_{q}^{o}\right\rangle)/\sqrt{2}$ represent the unit vectors in an
infinite dimensional Hilbert space $\mathcal{H}_{\infty}$, describing the
spatial mode on a transverse plane at a fixed $z=z_{0}$. On the other hand,
$\left|R\right\rangle$ and $\left|L\right\rangle$ represent the right– and
left– handed polarization basis (two–dimensional Hilbert space). In addition

$\displaystyle\left\langle\tilde{\xi},\tilde{\eta}|m_{q}^{+}\right\rangle$
$\displaystyle=u_{m}^{+}(\tilde{\xi},\tilde{\eta},z_{0};q),$ (8)
$\displaystyle\left\langle\tilde{\xi},\tilde{\eta}|m_{q}^{-}\right\rangle$
$\displaystyle=u_{m}^{-}(\tilde{\xi},\tilde{\eta},z_{0};q),$ (9)
$\displaystyle\left\langle\tilde{\xi},\tilde{\eta}|m_{q}^{e}\right\rangle$
$\displaystyle=\text{MG}^{e}(\tilde{\xi},\tilde{\eta},z_{0};q),$ (10)
$\displaystyle\left\langle\tilde{\xi},\tilde{\eta}|m_{q}^{o}\right\rangle$
$\displaystyle=\text{MG}^{o}(\tilde{\xi},\tilde{\eta},z_{0};q),$ (11)

where the following orthogonality rules are in place:
$\left\langle{m^{\prime}}_{q}^{\sigma^{\prime}}|m_{q}^{\sigma}\right\rangle=\delta_{m,m^{\prime}}\delta_{\sigma,\sigma^{\prime}}$
with $\sigma,\sigma^{\prime}=\\{e,o\\}$; and the inner products

$\displaystyle\left\langle m_{q^{\prime}}^{e}|m_{q}^{o}\right\rangle$
$\displaystyle=0,$ (12) $\displaystyle\left\langle
m_{q^{\prime}}^{e}|m_{q}^{e}\right\rangle$
$\displaystyle=\begin{cases}\sum_{j=0}^{\infty}A_{2j}^{2n}(q^{\prime})A_{2j}^{2n}(q),\text{for
$m=2n$},\\\
\sum_{j=0}^{\infty}A_{2j+1}^{2n+1}(q^{\prime})A_{2j+1}^{2n+1}(q),\text{for
$m=2n+1$},\\\ \end{cases}$ (13) $\displaystyle\left\langle
m_{q^{\prime}}^{o}|m_{q}^{o}\right\rangle$
$\displaystyle=\begin{cases}\sum_{j=0}^{\infty}B_{2j+1}^{2n+1}(q^{\prime})B_{2j+1}^{2n+1}(q),\text{for
$m=2n+1$},\\\
\sum_{j=0}^{\infty}B_{2j+2}^{2n+2}(q^{\prime})B_{2j+2}^{2n+2}(q),\text{for
$m=2n+2$},\\\ \end{cases}$ (14)

where $A_{j}^{m}(q)$ and $B_{j}^{m}(q)$ are the coefficients of the Fourier
expansion, corresponding to the angular even and odd Mathieu functions,
respectively (see A).

###  2.2 Bell inequalities

As explained before, the helical Mathieu-Gauss vector modes are constructed as
a non–separable superposition of the spatial and polarisation degrees of
freedom, giving rise to classically-entangled states. Even though a formal
demonstration that such beams are entangled in the classical sense has not
been provided yet, it can be demonstrated through the Clauser-Horne-Shimony-
Holt (CHSH) inequality, which is one of the most commonly used Bell-like
inequality for optical systems [18]. According to [13], the CHSH-Bell
parameter $S$ is defined as

$S=E(\theta_{1},\theta_{2})-E(\theta_{1},\theta_{2}^{\prime})+E(\theta_{1}^{\prime},\theta_{2})+E(\theta_{1}^{\prime},\theta_{2}^{\prime}).$
(15)

which, for any classical theory based on hidden variables the Bell parameter
has an upper bound $|S|\leq 2$. In the context of non-separable states of
light, $\theta_{1}$ and $\theta_{2}$ are the orientation angles of the half-
wave plate and the encoded hologram, respectively, as schematically shown in
Fig. 2 for the specific cases $\theta_{1}=0,\pi/4,\pi/2,3\pi/4,\pi$ and
$\theta_{2}=0,\pi/4,\pi/2,3\pi/4,\pi$.

Figure 2: Schematic representation of the projection onto the spatial (top
panels) and polarisation (bottom panels) DoFs.

Further, and also in the context of non-separable classical states
$E(\theta_{1},\theta_{2})$ is calculated from far-field intensity measurements
along the optical axis as

$\begin{split}&E(\theta_{1},\theta_{2})=\\\
&\frac{I(\theta_{1},\theta_{2})+I(\theta_{1}^{\perp},\theta_{2}^{\perp})-I(\theta_{1}^{\perp},\theta_{2})-I(\theta_{1},\theta_{2}^{\perp})}{I(\theta_{1},\theta_{2})+I(\theta_{1}^{\perp},\theta_{2}^{\perp})+I(\theta_{1}^{\perp},\theta_{2})+I(\theta_{1},\theta_{2}^{\perp})},\end{split}$
(16)

where $\theta_{i}^{\perp}=\theta_{i}+\pi/2$. Finally, the projection angles of
Eq. 15 are not arbitrary but restricted to the condition

$\theta\equiv\theta_{2}-\theta_{1}=\theta_{2}^{\prime}+\theta_{1}^{\prime}=-\theta_{2}-\theta_{1}^{\prime}.$
(17)

To measure the intensity $I$, we consider the “classically entangled” state as
described by Eq. 7. For clarity, let us remind the following inner products

$\displaystyle\left\langle m_{q^{\prime}}^{e}|m_{q}^{e}\right\rangle$
$\displaystyle\equiv\frac{1}{\pi}\int_{0}^{2\pi}\text{ce}_{m}(\eta,q^{\prime})\text{ce}_{m}(\eta,q)d\eta,$
$\displaystyle=\begin{cases}\sum_{j=0}^{\infty}A_{2j}^{2n}(q^{\prime})A_{2j}^{2n}(q),\text{
for $m=2n$},\\\
\sum_{j=0}^{\infty}A_{2j+1}^{2n+1}(q^{\prime})A_{2j+1}^{2n+1}(q),\text{ for
$m=2n+1$},\end{cases}$ (18) $\displaystyle\left\langle
m_{q^{\prime}}^{o}|m_{q}^{o}\right\rangle$
$\displaystyle\equiv\frac{1}{\pi}\int_{0}^{2\pi}\text{se}_{m}(\eta,q^{\prime})\text{se}_{m}(\eta,q)d\eta,$
$\displaystyle=\begin{cases}\sum_{j=0}^{\infty}B_{2j+1}^{2n+1}(q^{\prime})B_{2j+1}^{2n+1}(q),\text{
for $m=2n+1$},\\\
\sum_{j=0}^{\infty}B_{2j+2}^{2n+2}(q^{\prime})B_{2j+2}^{2n+2}(q),\text{ for
$m=2n+2$},\end{cases}$ (19)

where we have used Eqs. 32 and 33.

Then, we calculate the intensity associated to Eq. 7 by projecting the
polarization degree of freedom onto
$\left|R\right\rangle+\exp(2i\theta_{1})\left|L\right\rangle$, and the spatial
degree of freedom on
$\left|m^{+}_{q^{\prime}}\right\rangle+\exp(2i\theta_{2})\left|m^{-}_{q^{\prime}}\right\rangle$,
which after some algebra leads to

$\displaystyle I(\theta_{1},\theta_{2})=$
$\displaystyle\frac{1}{4}Q(m,q,q^{\prime})[1+\sin
2\alpha\cos(2\theta_{1}+2\theta_{2}+\beta)],$ (20)

where

$\displaystyle
Q(m,q,q^{\prime})=\bigg{|}\frac{1}{2}\sum_{j}[A_{j}^{m}(q^{\prime})A_{j}^{m}(q)+B_{j}^{m}(q^{\prime})B_{j}^{m}(q)]\bigg{|}^{2}.$

Here, it is worth highlighting that the case $q=q^{\prime}$ leads to
$Q(m,q,q)=1$ and therefore the intensity $I(\theta_{1},\theta_{2})$ reduces to
the well–known case reported in [13]. The Bell parameter $S$ reaches a maximum
value $2\sqrt{2}$ when $\theta=\pi/8$, where we obtain a maximally non-
separable state. Crucially, the value of $S$ does not depend neither on $q$
nor $q^{\prime}$, this can be seen by plugging Eq. 20 in Eq. 16, and realise
that $Q(m,q,q^{\prime})$ appears in both, the numerator and the denominator,
and cancels out from this equation. This is a crucial result that will be
demonstrated experimentally in the next section. As a final comment, it is
common to use the angles $\theta_{1}=0$, $\theta_{1}^{\prime}=-\pi/4$,
$\theta_{2}=\pi/8$ and $\theta_{2}^{\prime}=3\pi/8$ to violate the CHSH
inequality and these are precisely the angles we will use in the experiments
described in the next section.

##  3 Experimental results

###  3.1 Experimental setup

In order to corroborate experimentally our theoretical results, we implemented
a highly stable optical setup to generate arbitrary vector beams, which is
schematically illustrated in Fig. 3 and detailed in [19]. The setup starts
with a horizontally polarised HeNe laser beam ($\lambda=$ 632.8 nm),
collimated and expanded to cover from edge to edge the liquid crystal screen
of a Spatial Light Modulator (SLM1), which in our case is the Pluto LCOS phase
only from Holoeye with a spatial resolution of 1920$\times$1080 pixels and a
pixel size of 8 $\mu$m. The screen of the SLM is divided in two halves, each
of which is addressed with a digital hologram for the independent generation
and manipulations of an optical field. Each holograms is generated using
complex amplitude modulation, as detailed in [20], and superimposed with a
linear blazed grating to separate all diffraction orders and spatially filter
the first diffraction order using a telescope system formed by lenses L1 and
L2 ($f=500$ mm) and the spatial filter SF1. The beams emerge from the SLM with
a horizontal polarisation state which is then rotated to diagonal using a
Half-Wave Plate (HWP1). Afterwards, they are redirected to a common–path
triangular interferometer of the Sagnac type, which starts with a Polarising
Beam Splitter (PBS) and incorporates two mirrors. In this way, each beam is
separated by the PBS into its horizontal and vertical components travelling
inside the interferometer along opposite directions and after a round trip all
four beams exit from the port adjacent to the input one. The mirrors of the
interferometer are then adjusted to co-axially recombine the beam with
horizontal polarisation that emerges from one half of the SLM with the beam
with vertical polarisation that emerges from the other half, ensuring in this
way orthogonality in both, the polarisation and spatial Degrees of Freedom
(DOFs). The generated Helical Mathieu vector beam is then sent to our analysis
stage which consist on a Half-Wave plate (HWP2), a second spatial light
modulator (SLM2) and a Charge-Coupled Device (CCD) camera.

Figure 3: Schematic representation of the optical setup implemented to
generate helical Mathieu vector beams and perform projection measurements. BE:
Beam Expander; L1–L3: Lenses; PBS: Polarizing Beam Splitter, SLM: Spatial
Light Modulator; HWP1 and HWP2: Half-Wave plates; QWP: Quarter-Wave Plate;
M1–M3: Mirrors; SF1 and SF2: Spatial Filters; CCD: Charge Coupled Device
Camera

As a first step, we perform a simultaneous local correlation measurement on
both DOFs, polarisation and spatial mode [13]. The polarisation DOF is
measured by sending the vector beam through a polarisation filter, implemented
with a linear polariser oriented at different polarisation angles. The spatial
DOF is measured by modal decomposition using a series of modal filters
(digital holograms) encoded on a second spatial light modulator (SLM2) [20].
Here, it is important to note that since SLMs are polarisation-dependent the
diffracted light in the first diffraction order is essentially polarised in
the horizontal direction, acting as a horizontal linear polariser. Hence, to
measure the polarisation DOF, we can in principle rotate the SLM, an almost
impossible task, or alternatively insert a HWP before the SLM and rotate this
instead, a more practical approach. In regards to the spatial DOF, the input
vector field is projected into a series of phase holograms given by
$u_{m}^{+}(\vec{r};q^{\prime})+\exp(i2\theta_{2})u_{m}^{-}(\vec{r};q^{\prime})$,
for $\theta_{2}\in[0,\pi]$ rad. For the sake of clarity, Fig. 4 (a) shows
exemplary images of the projecting phases for
$\theta_{2}=0,\pi/4,\pi/2,3\pi/4$ and $\pi$. Afterwards, we measure the on-
axis intensity in the far field, achieved in the focal plane of a third lens
(L3). Figure 4 (b) shows example images of the far field intensity for the
different values of $\theta_{2}$ given above and for the specific case of
$\theta_{1}=0$, with the on-axis values marked with a red cross. The peak on-
axis intensity values of all images are normalised to the highest, and stored
for processing. For the sake of clarity, Fig. 4 (c) shows a table with such
values whereby, a maximum value of 1.00 was obtained for $\theta_{2}=0$ and
$\theta_{2}=\pi$, intermediate values of 0.48 and 0.53 were obtained for
$\theta_{2}=\pi/4$ and $\theta_{2}=3\pi/4$, respectively, and a minimum value
of 0.01 was obtained for $\theta_{2}=\pi/2$.

Figure 4: Far field intensity and the corresponding phase distribution on SLM2
for different $\theta_{2}$ values.

Figure 5 (a), shows the normalised intensity as a function of the projection
angle $\theta_{2}$ (the hologram displayed on the SLM) for three
configurations of the polarisation angle $\theta_{1}$, namely,
$\theta_{1}=0,\pi/8,\pi/4,3\pi/8$ and for the specific case
$q=q^{\prime}=0.09$ in Eq. 20. From this data and using Eqs. 15 and 16, the
calculation of the Bell parameter is straight forward and for our particular
example yields the value $S=2.88\pm 0.12$. Importantly, as predicted by Eq.
20, the maximum value of the intensity $I(\theta_{1},\theta_{2})$ is a
function of $q$ and $q^{\prime}$. More precisely, it decreases directly
proportional to the difference between $q$ and $q^{\prime}$ by the factor
$Q(m,q,q^{\prime})$. By way of example, we plotted the intensity given by Eq.
20 for the specific case $\theta_{1}=0$ and $q=0.09$, for
$q^{\prime}\in[0,6]$, which is shown in Fig. 5 (b). We corroborated this
theoretical prediction experimentally for the specific cases
$q^{\prime}=0.09,1.44,3.24$ and 4.41, which are also shown in 5 (b) as curves
with different colours, overlapping the theoretical curve.

Figure 5: (a) Bell–type curve from an input helical Mathieu–Gauss vector mode
with $m=3$ and $q=q^{\prime}=0.09$. The solid lines represent the theoretical
predictions, and the scattered points are the experimental intensities as a
function of $\theta_{2}$. (b) Bell–type surface as a function of $q^{\prime}$
and the orientation angle $\theta_{2}$ displayed on the hologram on SLM2. The
grey surface was computed through Eq. 20, whereas the coloured dots were
experimentally measured. The input HMGV beam was generated with $m=3$ and
$q=0.09$.

Crucially, such decrease in intensity does not affect the measurement of the
Bell parameter, as already mentioned above and observed by direct inspection
of Eq. 16, which evinces that the factor $Q(m,q,q^{\prime})$ is common to the
numerator and the denominator and therefore cancels out. To corroborate this
experimentally, we determined the Bell parameter $S$ for the specific case
$q=0.09$ and $q^{\prime}=5.76$. The curves of intensity as a function of
$\theta_{2}$, from which $S$ can be determined are shown in Fig. 6. Notice
that the maximum intensity value decreased significantly, but the shape of the
curves remains identical. Computation of $S$ from Eq. 15 yields the value
$S=2.51\pm 0.23$, which again indicates a violation of the CHSH-Bell
inequality. One of the main messages we can extract from this result is that
the CHSH-Bell inequality is also violated by Helical Mathieu-Gauss vector
beams, but more importantly, such violation can be determined regardless of
the ellipticity of the projection modes, given by the parameter $q$. A clear
limitation will be the sensitivity of the specific photodetector, since the
maximum and minimum intensity difference dramatically reduces as the
difference between $q$ and $q^{\prime}$ increases.

Figure 6: Bell–type curve from an input helical Mathieu–Gauss vector mode with
$m=3$ and $q=0.09$ $q^{\prime}=5.76$. The solid lines represent the
theoretical predictions, and the scattered points are the experimental
intensities as a function of $\theta_{2}$.

##  4 Conclusions

In summary, in this paper we demonstrated that Helical Mathieu-Gauss Vector
(HMGV) modes, characterised by a elliptical symmetry, are locally correlated
(entangled) in the classical sense, even though the distinguishing property of
entangled systems, non-locality, is absent. We demonstrated this by deriving a
general expression, in terms of the ellipticity of the HMGV modes that allows
to measure violations to the Bell inequality for optical systems, the Clauser-
Horn-Shimony-Holt-Bell (CHSH-Bell) inequality. In this context, the Bell $S$
parameter has an upper bound $|S|\leq 2$ but in our experiments we measure
values of $S\approx 2.88$, which shows a clear violation to the Bell
parameter. In the process we realised that $S$ is independent of the
ellipticity of the projection mode, which implies that the violation of the
same can be observed even when the ellipticity of the measured mode is not the
same as the projecting mode. Finally, our demonstration that HMGV modes also
violate Bell inequalities, generalises the concept of classical entanglement
to vector modes with arbitrary symmetry.

## Data availability statement

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

## Disclosures

The authors declare that there are no conflicts of interest related to this
article.

##  Appendix A Mathieu functions

The theory of Mathieu functions has been extensively studied in the past [21,
22, 23]. For instance, McLachlan [24] treats the theory and applications of
such functions in great detail. Even so, we will briefly outline some basic
concepts of Mathieu functions. Let us start by writing Helmholtz equation

$\displaystyle\nabla_{\perp}^{2}U+k^{2}U=0,$ (21)

where $\nabla_{\perp}^{2}$ is the transverse Laplacian and $k$ is some
constant. In elliptical coordinates, a solution to Eq. 21 of the form
$U=R(\xi)\Phi(\eta)$ satisfies the differential equations [25]

$\displaystyle\frac{d^{\text{2}}\Phi}{d\eta^{\text{2}}}+(a-2q\cos
2\eta)\Phi=0,$ (22)
$\displaystyle\frac{d^{\text{2}}R}{d\xi^{\text{2}}}-(a-2q\cosh 2\xi)R=0,$ (23)

known as the angular and radial Mathieu equations, respectively. The parameter
$q$ is a dimensionless parameter and $a$ is a separation constant. The
variable $\eta\in[0,2\pi)$ is analogous to the angular coordinate, and
$\xi\in[0,\infty)$ to the radial coordinate. Solutions to Eqs. 22 and 23 are
known as the angular and radial Mathieu functions, respectively.

Angular Mathieu functions can be expanded in Fourier series as follows [17]

$\displaystyle\text{ce}_{2n}(\eta,q)$
$\displaystyle=\sum_{j=0}^{\infty}A_{2j}^{2n}(q)\cos[2j\eta],$ (24)
$\displaystyle\text{ce}_{2n+1}(\eta,q)$
$\displaystyle=\sum_{j=0}^{\infty}A_{2j+1}^{2n+1}(q)\cos[(2j+1)\eta],$ (25)
$\displaystyle\text{se}_{2n+1}(\eta,q)$
$\displaystyle=\sum_{j=0}^{\infty}B_{2j+1}^{2n+1}(q)\cos[(2j+1)\eta],$ (26)
$\displaystyle\text{se}_{2n+2}(\eta,q)$
$\displaystyle=\sum_{j=0}^{\infty}B_{2j+2}^{2n+2}(q)\cos[(2j+2)\eta],$ (27)

where $A_{j}^{m}$ and $B_{j}^{m}$ are the Fourier coefficients and can be
computed via recurrence relations [17]. It is worth mentioning that Mathieu
functions and their coefficients are readily available, for example, in the
GSL numerical library [26] and the SciPy package [27]. A change of variable
$\eta=i\xi$ allows us to obtain the Radial Mathieu functions
$\text{Je}_{m}(\xi,q)$ (even) and $\text{Jo}_{m}(\xi,q)$ (odd). The
orthogonality rules for the angular Mathieu functions are

$\displaystyle\int_{0}^{2\pi}\text{ce}_{m}(\eta,q)\text{ce}_{n}(\eta,q)d\eta$
$\displaystyle=\pi\delta_{m,n},$ (28)
$\displaystyle\int_{0}^{2\pi}\text{se}_{m}(\eta,q)\text{se}_{n}(\eta,q)d\eta$
$\displaystyle=\pi\delta_{m,n},$ (29)
$\displaystyle\int_{0}^{2\pi}\text{ce}_{m}(\eta,q)\text{se}_{n}(\eta,q)d\eta$
$\displaystyle=0.$ (30)

Importantly, we can explore the cases with different $q$ parameter, namely

$\displaystyle\int_{0}^{2\pi}\text{ce}_{m}(\eta,q^{\prime})\text{se}_{m}(\eta,q)d\eta=$
$\displaystyle\begin{cases}\int_{0}^{2\pi}\sum_{j,j^{\prime}}^{\infty}A_{2j^{\prime}}^{2n}(q^{\prime})B_{2j+2}^{2n+2}(q)\cos(2j^{\prime}\eta)\sin[(2j+2)\eta]d\eta,\\\
\int_{0}^{2\pi}\sum_{j,j^{\prime}}^{\infty}A_{2j^{\prime}+1}^{2n+1}(q^{\prime})B_{2j+1}^{2n+1}(q)\cos[(2j^{\prime}+1)\eta]\sin[(2j+1)\eta]d\eta,\end{cases}$
$\displaystyle=\begin{cases}\sum_{j,j^{\prime}}^{\infty}A_{2j^{\prime}}^{2n}(q^{\prime})B_{2j+2}^{2n+2}(q)\int_{0}^{2\pi}\cos(2j^{\prime}\eta)\sin[(2j+2)\eta]d\eta,\\\
\sum_{j,j^{\prime}}^{\infty}A_{2j^{\prime}+1}^{2n+1}(q^{\prime})B_{2j+1}^{2n+1}(q)\int_{0}^{2\pi}\cos[(2j^{\prime}+1)\eta]\sin[(2j+1)\eta]d\eta,\end{cases}$
$\displaystyle=0.$ (31)

Moreover

$\displaystyle\int_{0}^{2\pi}\text{ce}_{m}(\eta,q^{\prime})\text{ce}_{m}(\eta,q)d\eta=$
$\displaystyle\begin{cases}\int_{0}^{2\pi}\sum_{j,j^{\prime}}^{\infty}A_{2j^{\prime}}^{2n}(q^{\prime})A_{2j}^{2n}(q)\cos(2j^{\prime}\eta)\cos(2j\eta)d\eta,\\\
\int_{0}^{2\pi}\sum_{j,j^{\prime}}^{\infty}A_{2j^{\prime}+1}^{2n+1}(q^{\prime})A_{2j+1}^{2n+1}(q)\cos[(2j^{\prime}+1)\eta]\cos[(2j+1)\eta]d\eta,\end{cases}$
$\displaystyle=\begin{cases}\pi\sum_{j=0}^{\infty}A_{2j}^{2n}(q^{\prime})A_{2j}^{2n}(q),\text{
for $m=2n$},\\\
\pi\sum_{j=0}^{\infty}A_{2j+1}^{2n+1}(q^{\prime})A_{2j+1}^{2n+1}(q),\text{ for
$m=2n+1$}.\end{cases}$ (32)

Similarly, for odd functions

$\displaystyle\int_{0}^{2\pi}\text{se}_{m}(\eta,q^{\prime})\text{se}_{m}(\eta,q)d\eta=$
$\displaystyle\begin{cases}\pi\sum_{j=0}^{\infty}B_{2j+1}^{2n+1}(q^{\prime})B_{2j+1}^{2n+1}(q),\text{
for $m=2n+1$},\\\
\pi\sum_{j=0}^{\infty}B_{2j+2}^{2n+2}(q^{\prime})B_{2j+2}^{2n+2}(q),\text{ for
$m=2n+2$}.\end{cases}$ (33)

## References

## References

  * [1] Einstein A, Podolsky B and Rosen N 1935 Phys. Rev. 47(10) 777–780
  * [2] Bell J S 1964 Physics Physique Fizika 1(3) 195–200
  * [3] Bell J S 1966 Rev. Mod. Phys. 38(3) 447–452
  * [4] Clauser J F, Horne M A, Shimony A and Holt R A 1969 Phys. Rev. Lett. 23(15) 880–884
  * [5] Ndagano B, Sroor H, McLaren M, Rosales-Guzmán C and Forbes A 2016 Opt. Lett. 41 3407
  * [6] Forbes A, Aiello A and Ndagano B 2019 Classically entangled light Progress in Optics (Elsevier Ltd.) pp 99–153
  * [7] Selyem A, Rosales-Guzmán C, Croke S, Forbes A and Franke-Arnold S 2019 Phys. Rev. A 100(6) 063842
  * [8] Qian X F and Eberly J H 2011 Opt. Lett. 36 4110–4112
  * [9] Aiello A, Töppel F, Marquardt C, Giacobino E and Leuchs G 2015 New J. Phys. 17 043024
  * [10] Konrad T and Forbes A 2019 Contemporary Physics 1–22
  * [11] Yao-Li, Hu X B, Perez-Garcia B, Bo-Zhao, Gao W, Zhu Z H and Rosales-Guzmán C 2020 Applied Physics Letters 116 221105
  * [12] Simon B N, Simon S, Gori F, Santarsiero M, Borghi R, Mukunda N and Simon R 2010 Phys. Rev. Lett. 104(2) 023901 URL https://link.aps.org/doi/10.1103/PhysRevLett.104.023901
  * [13] McLaren M, Konrad T and Forbes A 2015 Phys. Rev. A 92 023833
  * [14] Otte E, Nape I, Rosales-Guzmán C, Vallés A, Denz C and Forbes A 2018 Phys. Rev. A 98(5) 053818
  * [15] Rosales-Guzmán C, Hu X, Rodríguez-Fajardo V, Hernandez-Aranda R I, Forbes A and Perez-Garcia B 2021 J. Opt. 23 034004
  * [16] Gutiérrez-Vega J C and Bandres M A 2005 J. Opt. Soc. Am. A 22 289–298
  * [17] Olver F, of Standards N I, (US) T, Lozier D, Boisvert R and Clark C 2010 NIST Handbook of Mathematical Functions Hardback and CD-ROM (Cambridge University Press) ISBN 9780521192255
  * [18] Toninelli E, Ndagano B, Vallés A, Sephton B, Nape I, Ambrosio A, Capasso F, Padgett M J and Forbes A 2019 Advances in Optics and Photonics 11 67–134
  * [19] Perez-Garcia B, López-Mariscal C, Hernandez-Aranda R I and Gutiérrez-Vega J C 2017 Appl. Opt. 56 6967–6972
  * [20] Rosales-Guzmán C and Forbes A 2017 How to shape light with spatial light modulators SPIE.SPOTLIGHT (SPIE Press)
  * [21] Kirkpatrick E 1960 Mathematics of Computation 118–129
  * [22] Gutiérrez-Vega J C, Rodríguez-Dagnino R M, Meneses-Nava M A and Chávez-Cerda S 2003 American Journal of Physics 71 233–242
  * [23] Meixner J, Schäfke F W and Wolf G 2006 Mathieu functions and spheroidal functions and their mathematical foundations: further studies vol 837 (Springer)
  * [24] Arscott F 1968 The Mathematical Gazette 52 94–95
  * [25] Mathieu E 1868 Journal de Mathemátiques Pures et Appliquées 13 137–203
  * [26] Galassi M et al, GNU Scientific Library Reference Manual (3rd Ed.)
  * [27] 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, Carey C J, Polat İ, Feng Y, Moore E W, VanderPlas J, Laxalde D, Perktold J, Cimrman R, Henriksen I, Quintero E A, Harris C R, Archibald A M, Ribeiro A H, Pedregosa F, van Mulbregt P and SciPy 10 Contributors 2020 Nature Methods 17 261–272