High-performance software typically uses dynamic memory allocations and multi-threading to leverage multi-core CPUs. Model checking such software not only has to deal with state space explosion, but also with variable-length states due to dynamic allocations. Moreover, changes between states are typically small, calling for incremental updates. Many model checkers, although efficiently dealing with the latter, only support fixed-length state vectors. In this paper, we introduce , a concurrent compression tree data structure that compactly stores variable-length states while allowing partial state reconstruction and incremental updates without reconstructing states. We implemented in the multi-core model checker. We show that, for models with states of varying length, is up to 2.9 times faster and uses on average 29% less memory than state-of-the-art tools.

§ INTRODUCTION

High-performance concurrent software is complex to write and even harder to reason about.
The more threads run in parallel, the more different interleavings are possible, easily causing billions (or trillions) of reachable program states.
Programmers want to make sure that they are valid, i.e. none represent an error due to stack overflows, race conditions, buffer overruns, null pointer dereferences, or other erroneous operations.
One way to verify such properties is to model check <cit.> the program, during which the entire state space is explored.
We can remember the visited states to avoid visiting them multiple times and thus avoid doing the same work redundantly.
A complication when model checking software is dealing with variable-length states, due to dynamic memory allocations, e.g. heap memory.
However, changes between states are typically small. Model checkers can take advantage of this for efficiency.
Thus, model checking software has four major requirements:

R1) efficiently store billions or trillions of states;
R2) efficiently calculate successor states;
R3) efficiently determine whether or not a state has been visited;
R4) support states of varying length, due to dynamic memory.
These requirements have implications on how states are stored, but also on how next states are determined and communicated
between model, model checker and state storage, in a modularized implementation as shown on the right.

§.§ Related work

LTSmin <cit.>, together with its compression trees  <cit.> and  <cit.> is one example on how to approach these requirements in a modular implementation. LTSmin is a model checker with multi-core <cit.> (up to 64 threads), distributed <cit.> and symbolic implementations <cit.>.
These implementations share a common interface: PINS, the Partitioned Next-State Interface, their API between search core and model.
PINS does not lend itself well for software, because it only support fixed-length states. Moreover, it requires the complete state to be available when the model is asked for successor states.
Because these states are stored in or , these states need to be reconstructed first, which comes at a price. Communicating a next state is also done using a complete state. LTSmin can use a projection <cit.> to detect the changed parts to minimize copying.
However, this projection requires static information on which parts of the state are touched, which is not guaranteed available in software due to the dynamic nature of heap memory and spawning threads dynamically.
LTSmin does have a software front-end <cit.>, but since LTSmin only supports fixed-length states, that front-end uses a sub-optimal uncompressed chunk table to model heap memory, which also lacks projection.

<cit.> is a thread-safe state storage for fixed-length states.
It uses a binary compression tree, based on work by Blom et al. <cit.>, in order to identify common sub-vectors.
Laarman et al. show it to be efficient (R1, R2 and R3), but it does not support states of varying length (R4).
<cit.> is an evolution of , where a Cleary table <cit.> is used for its root set. We will explain and in more detail in <ref>.

SPIN version 6 <cit.> is an explicit-state model checker for models. It supports a dynamic number of processes by having a variable length root state. Moreover, it has a state compression method  <cit.>: instead of fully storing all combinations of process states in the root state, the states of the processes themselves are stored and mapped to a unique ID, and those IDs are stored in the root state instead.
Thus, a state is a two-level tree of states: a root state and the process states. This is similar to using the chunk table of LTSmin.
SPIN supports a number of search algorithms, among which is parallel BFS up to 63 threads.

version 4 <cit.> is an explicit-state model checker for LLVM IR assembly code <cit.>. It uses a graph to model the stack and heap. Each node in this graph represents an allocated section of memory, for example space allocated for a struct. Edges on this graph represent reachability, e.g. a pointer in a struct to another struct.

Other, more remote approaches to state compression, are Binary Decision Diagrams <cit.> and PTries <cit.>.
BDDs compress by sharing prefixes and suffixes of (Boolean) state vectors, and are used in symbolic model checking; in their standard form, they are restricted to static vector lengths.
PTries compress by only sharing prefixes of subvectors; they do support dynamic state lengths natively.
Laarman compares the compression and performance of BDDs and PTries to  <cit.>.

§.§ Contributions

In this paper we introduce , a state storage data structure that satisfies R1–R4: it is concurrent, compresses states, can handle variable-length states and provides partial reconstruction of a stored state and incremental updates to states without the need for reconstructing the entire state.
We implemented as part of the .

To showcase the potential of , we compare it to other state storage components using models with variable-length states. In the case of , we pad all vectors to the length of the largest vector (manually determined). Here, is able to perform up to 2.9 times faster than using 29% less memory on average, without the a priori need to manually determine the largest vector.
The advantage of increases as the difference in length between states increases.
To evaluate how much efficiency is lost to gain variable-length support, we compare using to LTSmin and SPIN using models from the BEEM database <cit.> with fixed-length states.
In this setting, is actually marginally faster than , but uses 30% more memory.
compresses fixed-length states 2.3x more than .
Compared to SPIN using , is 8.1x faster and uses 6.1x less memory.

§ TREE COMPRESSION

The layout of the indexed hash set after inserting colabcdabcd and colabcdefabcdef.

We now explain and , which inspired .
States are vectors of state slots; a variable in a model occupies one or more state slots.
For the sake of simplicity, let us assume these state slots are 32-bit.
The principle behind the compression is to insert pairs of state slots, thus 64-bit values, into an indexed hash set with a 32-bit index.
In an indexed hash set, the index is determined by hashing the value.
A pair of 32-bit indices forms another 64-bit value, which in turn can be put into the set.
This process continues, creating an ordered tree of nodes.
This is illustrated in <ref>: the vectors colabcdabcd and colabcdefabcdef share a common sub-vector whiteabcd and index white3 maps to this sub-vector.
This also illustrates the compression: common nodes are only inserted once.
The higher up the common node is in the tree, the larger the common sub-vector it maps to, the more memory is saved.
Which common sub-vectors can be identified and reused, depends on the structure of the binary tree: while whiteabcd can be identified, whitecdef cannot.
If we were to add the vector whitecdef, we would add the node white24. Similarly, we would add the node white42 for whiteefcd and the node white33 for whiteabcdabcd.

Note that the rounded boxes colabcd3 and colabcdef5 are the indices to the root nodes of the trees mapping to the states. However, an index is not enough to identify a state, as index white3 could map to white12 as well. Indeed, if we would have added the state colfour12 (with index colfour3), this would be completely correct, thus we need more information other than white3 to distinguish colabcd3 from colfour3. This gives rise to two issues: uniquely identifying a state and determining whether or not a state has already been added.

§.§ and its implementation

Laarman et al.'s compression tree schema as in their paper.
The top state adds 7 nodes, bringing the total to 16.

Laarman et al.'s compression tree schema as in their implementation.
The top state adds 10 nodes; 19 nodes in total.

An illustration of the difference between the paper and implementation of the compression tree of Laarman et al. Given a 40 byte state colabcdefghijabcdefghij (bottom tree), the next state colabcdefghijkabcdefghijk (top tree) appends colabcdefghijkk, 4 bytes.

solves the first issue by only supporting fixed-length states and having the length of the state dictate the shape of the tree. Thus, all trees are isomorphic and an index uniquely identifies a state. However, the implementation <cit.> of how the length dictates the shape of the tree differs from the how it was originally presented <cit.>.
We show both versions in <ref>, ignoring colabcdefghijkX (top) and only focusing on colabcdefghijX (bottom) for now.
The paper version recursively divides the vector in two, associating one half with the left child of a node and the other half with the right child (<ref>). For efficiency, the implementation version represents the tree as an array, where a node at index $n$ has its children at indices $2n$ and $2n+1$ (<ref>). Inserting the nodes of a vector then becomes a simple for-loop since children are always adjacent. Both versions produce a balanced tree, but their balancing is different, influencing locality. For example, the paper version has a single pair of adjacent state slots (colabcdefghijef) that requires traversing both the left and right side of the tree, whereas the implementation version has two (colabcdefghijbc and colabcdefghijfg). The importance of this will become apparent in <ref>.

To determine whether a state has been added already, the paper version of uses a single indexed hash set and a root-bit in the nodes to indicate whether a node is a root node or a tree node. This distinction is needed because otherwise a leaf node could be mistaken for a root node and it could be erroneously concluded a state was already added.
The current implementation removes the root-bit in favour of using two hash sets: a root hash set that stores the root nodes and a hash set that stores all other nodes.
Not needing the root-bit, 32-bit indices can be used instead of just 31-bit indices and this lifts the 31-bit restriction on the data in state slots as well.
In addition, it allows the root set to be larger than $2^{32}$, since its indices need not fit in 32 bits. The data set is still limited to 32-bit indices, because paired they are limited to 64 bits.
However, one specific 64-bit value in the hash map is still reserved to indicate the empty value. This value is needed to distinguish between an empty bucket and a used bucket. In the case of , this is the value $-1$, which means a model cannot pair two state slots that have a $-1$ value, because together they would form a 64-bit $-1$ value. As we will see in the experiments, this does happen.

<cit.> is the same as , with the root set replaced by a Cleary Table <cit.>. supports root set indices of fewer than 64 bits, but can only return a 64-bit root index, due to the nature of the Cleary Table. Furthermore, the paper on shows the same manner of determining the shape of the tree as the paper, but the implementation is shared with , with a flag to enable the use of a Cleary Table. This means Compact Tree suffers from the same alignment issues as .

§ DTREE

Like , we base our compression tree on the compression tree of Blom et al. <cit.> and use two indexed hash sets. We diverge on a number of crucial points from in order to meet the requirements mentioned in <ref>. Since supports only fixed-length vectors, there is no need to remember the length of each individual vector, but we do (R4).
We also introduce a different compression tree structure that benefits dynamic memory allocation (R4).
Furthermore, we extend the capabilities of the compression tree to store and reconstruct only parts of states in order to improve calculating next states (R2). Lastly, our approach supports any data, including $-1$ values in adjacent state slots.

First, we show why we need a different tree structure.
Let us consider the scenario where a model has allocated 40 bytes (10 state slots) and in the next state the model wishes to allocate 4 more bytes. The model does this by growing the state by one state slot. <Ref> shows how this scenario is supported by the two versions of .
The implementation version starts forming 64-bit sections starting at the end, so a change of 4 bytes (32 bits) to the size of the state will unalign the tree nodes compared to the previous state. This results in storing the same information twice. This effect is not limited to the last level of the tree: a change of 8 bytes in the size of the state results in the second to last level of the tree to become unaligned. A change of 16 bytes would unalign the third to last level, etc.
<Ref> shows two improvements to these structures. <Ref> describes the structure of the implementation version of , but backwards in an attempt to limit the alignment issue. This works for the leaves, but one level higher in the tree we face the same issue: instead of colabcdefghij12, colabcdefghij34 we get colabcdefghijk23, colabcdefghijk45. Moreover, this approach still has two pairs of adjacent state slots that require traversing both the left and right side of the tree (colabcdefghijde and colabcdefghijhi).

§.§ A Chain of Perfectly Balanced Binary Trees

Laarman et al.'s compression tree schema as in their implementation, but left-to-right instead of right-to-left.
The next state adds 5 nodes; 14 nodes in total.

The two states using our chain of perfectly balanced binary trees schema.
The next state adds 2 nodes, bringing the total to 11. Additionally, it shows a third state colabcdefghijklX when colabcdefghijkll is appended.

Like <ref>, given a 40 byte state colabcdefghijX, the next state colabcdefghijkX appends colabcdefghijkk, 4 bytes.

<Ref> shows our approach, named the chain of perfectly balanced binary trees. This approach does not suffer from the issues of the other tree shapes as <ref> clearly demonstrates. In a chain, the length of the sub-vector the left child of a node leads to is always $\mbox{lpst}(\mbox{len}(V))$, where $\mbox{lpst}(x)$ is the largest power-of-two smaller than $x$, $\mbox{len}(V)$ is the length of the vector $V$ and $V$ is the vector that the node itself leads to.
Thus, a left child is always perfectly balanced. The right part maps to the remainder of the vector.
For example, the state colabcdefghijkabcdefghijk of length 11 in <ref>, is shaped by perfectly balanced trees mapping to 8, 2 and 1 state slots; the state colabcdefghijklabcdefghijkl is shaped by trees mapping to 8 and 4 state slots. This can result in a less balanced tree than . For example, for the state colabcdefghijabcdefghij, (<ref>) has 2 nodes on the lowest level (colabcdefghijgh, colabcdefghijij) and 3 on the second-to-lowest level (colabcdefghijab, colabcdefghijcd, colabcdefghijef). (<ref>) on the other hand has 4 on the lowest level (colabcdefghijab, colabcdefghijcd, colabcdefghijef, colabcdefghijgh) and one on the second-to-highest level (colabcdefghijij). This means changes to colabcdefghijab, colabcdefghijcd or colabcdefghijef require an extra node and changes to colabcdefghijij require two fewer nodes, compared to .
We will investigate the impact of the difference in <ref>.
<Ref> illustrates that this approach does lend itself to appending to states.
This targets a combination of R1 and R4 from <ref>: to more efficiently handle dynamic memory allocation.
It also has a single pair of adjacent state slots requiring traversing both sides of the tree (colabcdefghijhi).

§.§ Incremental updates

The importance of minimizing such adjacent pairs becomes apparent when adding support for incremental updates. With this, a model conveys to the model checker only what changed compared to a previous state, in lieu of a complete state.
supports incremental updates by detecting unchanged parts, when inserting a complete state, but it does not support simply applying an incremental update to a state. We describe the difference with an example. Considering the state colabcdefabcdef, the model determines that in the next state colabcdefcd changes to colabghefgh. We have two options of communicating this to the model checker:

* the model can communicate the entire new state colabghefabghef; or

* the model can communicate only the change of colabghefgh at offset 2, denoted colabghefgh@2.

For small states the difference in performance is likely negligible. However, since we want to model check software, our states can be quite large and thus this could save a significant amount of copying.

When communicating only changes, the model checker also needs to know which state the change needs to be applied to. colabcdef34 (as per <ref>) alone does not identify colabcdefabcdef: we need the length of the vector as well.

We considered three ways to solve this:

* Always remember the length alongside the index. This has the downside that high-performance low-level atomic instructions such as compare-and-swap cannot be used. 128-bit atomic instructions are more expensive <cit.>.

* Add another node to the tree, containing the length and index to the remainder of the tree. Thus, the index to that node would uniquely identify a state. This adds a level to the tree, which can negatively influence performance.

* Remember index and length, but limit their combination to 64 bits. For example, 40 bits for the index and 24 for the length. This has the disadvantage of only supporting a trillion states of $2^{24}-1$ state slots each. At the current technological age, both seem upperbounds on what we need for the foreseeable future, so this should suffice.

delta() is called recursively

Applied delta colabghefgh@2 to colabcdefabcdef

The delta-vector colabghefgh is applied at offset 2 to colabcdefabcdef from <ref>.

We chose the third alternative. We call this combination the StateID and we write this for example as colabcdef56, which uniquely identifies colabcdefabcdef: 5 is the index, 6 is the length of the state.
Applying a delta (incremental update) such as colabghefgh@2 to colabcdef56 is described in <ref>.
First, in <ref>, we recursively traverse the tree to find the leaves that need change. At each node that leads to vector $V$, we check which of the children are affected. From <ref> follows that the left child is affected iff the offset is smaller than $\mbox{lpst}(\mbox{len}(V))$. The right child is affected iff the offset plus the length of the delta is greater than $\mbox{lpst}(\mbox{len}(V))$. Thus, the tree is traversed as follows:

whitegh@2 to colabcdef56 only affects the left child, since $2<\mbox{lpst}(6)$ and $2+2 \not>\mbox{lpst}(6)$

whitegh@2 to colabcdef34 only affects the left child, since $2\not<\mbox{lpst}(4)$ and $2+2 >\mbox{lpst}(4)$

whitegh@0 to colabcdef22 completely replaces the leaf node.

At this point we can traverse back, inserting new nodes along the path. First we insert the new node colabghefgh at index colabghef6 in the data set, then we use that index to create a new node colabghef16, where we copy the index of the left child from colabcdef12. Finally, we insert the root node colabghef74, providing StateID colabghef86 (<ref>).
A delta beyond the length of a state is also supported. For example, we could have applied whitegh@5 to colabcdef56, resulting in whiteabcdegh. When applying whitegh@8 to colabcdef56 it would yield whiteabcdef00gh.

§.§ Partial reconstruction of states

In a similar fashion to incremental updates, we can reconstruct (parts of) states. We recursively traverse the tree of nodes until we reach the correct leaf nodes. Then, we copy the contents of these leaf nodes into a buffer, yielding (a part of) a state.

§.§ Determining a state is new

When creating a new state, either by inserting a complete new one or by incremental update, the model checker needs to know if the state has already been visited or not. This is accomplished by use of the root set. A state has already been inserted into iff there is a corresponding root node. As an example, let us insert colabcdefabcdef into the of <ref>, which already exists. Recursively, we traverse the nodes. At the leaves, we conclude whiteab, whitecd and whiteef have already been inserted.
The tree node white12 also has already been inserted. Then we conclude that the root node white34 has been inserted as well.
Now let us insert colabcdabcd. Again, we conclude the leave nodes whiteab and whitecd are already in the tree. Even though white12 is already in the data set, it is not yet in the root set and thus we correctly conclude it is a new state. Additionally, supports inserting states/vectors that are inserted purely in the data set and not its root node in the root set. This can be useful to compress data that is not a state in the model as we will see later in <ref>. Thus, provides the following API:

StateID insert(Slot[] V, bool root)
StateID delta(StateID s, int offset, Slot[] D, bool root)
Slot[]  get(StateID s, bool root)
Slot[]  get(StateID s, int offset, int length, bool root)

Here, root indicates if the root node should indeed be placed in the root set. Note that this is also needed for get() because we need to look for the root node in the right set. The interfaces insert() and delta() insert a new state (<ref>). The two get() interfaces allow to obtain (parts of) the state (<ref>).

§ DMC MODEL CHECKER

Q = ;
IL^StateID s = Q.pop();

A basic search core that simply requests next states from a model M until all have been visited. For brevity, the on-the-fly checking of properties is left out.

void initialState()
IL^MC.insert(0,0,0,0, true);
IL># ^ implicitly pushes to Q
void nextStates(StateID s)
IL|for i in 0,1,2,3
IL|IL^Slot v = MC.get(s, i, 1);
IL|IL|v = (v + 1)
IL>IL>MC.delta(s, i, v, true);

A model with four counters.

A search core and a simple model with a state space size of $10^4$.

StateID delta(StateID s, int offset, Slot[] D, bool root)
IL^StateID s, bool isNew = DTREE.delta(s, offset, D, root);
IL|if(isNew and root) Q.push(s);
IL>return s;

An implementation of the delta() interface of the .

The modular DMC Model Checker.

We now expand on the concept of a modular model checker as described in <ref>. A model M implements the (NS) API:

void initialState()
void nextStates(StateID s)

Here, initialState() sets up the initial state and nextStates() computes the next states of the specified state s. Communicating (“uploading”) states to the model checker is done using the API, which is exposed to the model via the search core. This flow is illustrated in <ref>.

The API is called by a search core, which defines the search strategy. For example, in <ref> we define a search core that simply pops a StateID from a queue Q and requests the next states of the model M until the queue is empty. An example of such a model M is shown in <ref>, modeling four counters, going from 0 to 9, looping back to 0.
With each invocation of nextStates(), four next states are generated. For each Slot 0–3 in the state (line 5), the current value is obtained (line 6), it is incremented modulo 10 (line 7) and this delta is communicated to the model checker (line 8) using the delta() interface.
Thus, the initial state is white0000 and for example
white0974 generates the next states white1974, white0074, white0984 and white0975.
The search core wraps the interface to add new states to the queue Q.
An example of this is shown in <ref>. Line 2 inserts the potentially new state in and is returned a StateID and an indicator if the state was new. If the state s is new and a root state (we will see a use case for non-root states in <ref>), it is added to the queue Q in line 3. The insert() interface is implemented similarly, updating Q if needed. The get() interfaces are simply passed on.

§.§ An example using a tree of states

A tree of states.

Four processes counting from 0 to 9 and wrapping.

struct Process
IL^int pc;      # program counter
IL>int i;       # some variable
struct SV
IL^int n;       # number of processes
IL>StateID p[];
void nextStates(StateID rootID)
IL^SV sv = MC.get(rootID, true)                           # root state
IL|for(p in [0..sv.n])
IL|IL^int pi = MC.get(sv.p[p],    Process::i, 1, false)
IL|IL|pi = (pi + 1)
IL|IL|int newP = MC.delta(sv.p[p],    Process::i, pi, false);
IL|IL|int offset =    SV::p[p]                 # offset to p[p] within SV
IL>IL>MC.delta(rootID, offset, newP, true);              # root state
void initialState(StateID rootID)
IL^Process p = 1,0                                 # pc is 1, i is 0
IL|StateID initP = MC.insert(p, false)
IL|SV rootState = 4, initP, initP, initP, initP
IL>MC.insert(rootState, true)                             # root state

We can use the basic building blocks of uploading and downloading states for a more interesting concept: tree-structured states.
The basic idea is that a state can contain a StateID that is associated with a non-root state, a sub-state. These are called sub-states because they could well be states of processes of which the combined state is the root state.

To illustrate this, we constructed a small model using processes in <ref>, which conceptually does the same as the model in <ref>.

<Ref> illustrates what a complete state looks like. The important lines are 10–12, where we obtain the value of i of the current process (10), change it (11) and create a new sub-state with the changed value (12). Note that the get() in line 8 has root=true, indicating that the specified StateID is a root state, while the get() in line 10 does not. This matches with the corresponding insert()s and delta()s.

This is much like the method of SPIN <cit.>

as briefly touched on in <ref>. SPIN supports a varying number of processes by having exactly this kind of structure: a root state with indices to indicate the state of individual processes. Just as , this allows to leverage the fact that duplicate processes share states. Since these states are inserted separately instead of a single, large vector, their data aligns, which helps in identifying common sub-vectors. However, supports arbitrary deep hierarchies of nested states.

§ EXPERIMENTS

To be able to compare to , we have two options: 1) initialise with the length of the largest state and pad others with 0's; 2) initialise for the most common length and use a different storage component for other lengths. We implemented both options we call and $S$, where $S$ is the other storage component. Without these measures, would not support variable-length states and could not run the variable-length experiments.

While we compared to std::unordered_map of the STL library surrounded by an std::mutex, and a concurrent chaining hash map (cchm) <cit.>, we omit the results for them as in all cases they perform roughly an order of magnitude worse. This is not unexpected, since they store an order of magnitude more memory because they do not compress the inserted states. The hash map std::unordered_map specifically has a global mutex around it, causing high contention.
We could not make into a DMC storage component because it needs a 64-bit index and the DMC API reserves only 40 bits for that purpose.

In addition to comparing these storages to each other, we compare DMC to the multicore implementation of LTSmin 3.0.2 with both and , both parallel DFS and BFS, using 64 threads. We also compare to SPIN 6.5.1 using and parallel BFS with 63 threads, as that is SPIN's maximum.
This gives us a baseline of where the performance of our model checker DMC as a whole is.

We ran DMC at 64 threads for a fair comparison, but it can use any number of threads. We tried both the DFS and BFS strategy for LTSmin, but these showed only minor difference. We will show the results for BFS, as we run SPIN using BFS as well and DMC currently only has a BFS search core.

The hardware that we ran our experiments on is “caserta”, a Dell R930 with 2TiB of RAM and four E7-8890-v4 CPUs.
Each CPU has 24 cores, 60MiB of L3 cache and 512GiB of RAM, offering 96 physical cores in total. We ran our experiments on Linux 4.4.0 and all tools were compiled using GCC 9.3.0.

§.§ Variable-length state models

The design of is meant for software that uses dynamic memory. To test this we implemented three models: 1) a model that implements a concurrent hash map with a number of inserts <cit.>; 2) a model of the concurrent Michael Scott queue <cit.> with various enqueue (E) and dequeue (D) operations; 3) a model of a sorted linked list (SSL).
All these data structures are modeled by creating a model that inserts a number of elements into them. The initial state is constructed by creating a number of processes, like in the example in <ref>. These processes then run, with every interleaving explored, akin to how every interleaving of incrementing one of the four counters in the example of <ref> is explored. For the Michael-Scott queue, we made a model that runs 3 enqueue operations and 3 dequeue operations in parallel. The Sorted Linked List model models a number of processes (6) that each insert a single element of 48 bytes (12 slots), dynamically allocated. The hash map model similarly inserts 9 pointers to elements of 16 slots in parallel. Dynamic memory is modeled using a StateID at the beginning of the root state to a memory slab sub-state that can expand.

The effect these dynamic allocations and parallel insertions have on the distribution of the length of states is shown in <ref>. For example, the hash map model is most pronounced: starting with an initial global memory of 130 slots (some string data), the first process dynamically allocates 16 slots and inserts the element, yielding a number of states of length 146. If we add a second process and all the interleavings, we get many more states of length 162, etc.

MSQ 3E+3E

SSL (6)

Hashmap 8/64

The number of states (y-axis) inserted of a certain length in ints (x-axis). The length of the root state is highlighted.

Experiments that use the . Set scale X-Y(-Z) means $2^X$ root set nodes, $2^Y$ data set nodes, $2^Z$ entries for the sub-storage. cchm($L$) and ($L$) mean they are initialized for states of length $L$ (in ints).

Set Scale

black!15!white-3.6cmSLL(6 inserts)

black!15!white-3.6cmHashmap(8 inserts)

black!15!white-3.6cmHashmap(9 inserts)

The results for variable-length state models are shown in <ref>.
Overall, is clearly the faster of the three storages. It is up to 44 times faster than cchm and 1.2–2.9 times faster than . The slow times of cchm are largely caused by the slower, uncompressed cchm that is used for states of lengths other than the root state so the results have to be interpreted as such. For , we see the downside of padding with zeroes: with increased state-length variance, the overhead increases and performance drops.

Regarding compression, we notice cchm is actually better for the MSQ and Sorted Linked List models. This can be explained by looking at the variance of the distribution of state lengths: the root state length (18) dominates all other lengths. The shape of the tree of nodes of has 2 nodes on the lowest level and 7 on the second-to-lowest. on the other hand, has 8 nodes on the deepest level and 1 on the second-to-highest, because it uses a chain of balanced trees, in this case leading to 16 state slots (left) and 2 (right).
As theorised in <ref>, this difference causes that often needs one more node for even a small change. Since often the memory is changed, which is modeled by a StateID at the beginning of the root state, this is precisely what happens.
Thus, the root state length is so dominant that the overhead of uncompressed cchm entries is less than the overhead of the less balanced tree that uses.

When we look at the Hashmap model, which has a significantly higher variance, we see that cchm requires an order of magnitude more space. Here, stores only states of length 20 and all others are stored in the sub-storage cchm, which does not compress states.

In general these results show that dedicated support for variable-length outperforms padding zeroes () and offloading other-sized vectors (cchm). The results also show that an increase of variance in state-length increases the advantage of .

§.§ Fixed-length state models

DMC/dtree compared to
LTSmin/compact (*fig:experiments:beemresults:vscompacttime,
and SPIN (*fig:experiments:beemresults:vsspintime,
*fig:experiments:beemresults:vsspinbytesperstate) on runtime
*fig:experiments:beemresults:vsspintime) and bytes per state

To evaluate the cost of adding variable-length support, we also run experiments using models with fixed-length states.
We modified SpinS <cit.> to emit models that implement the API of DMC for models from the BEEM database <cit.>. For SPIN, these states could be of varying length, but SpinS emits fixed-length states models since its primary target is PINS for LTSmin. Thus, these tests do not use the delta() interface, but only insert and get complete states. Of these models, we include the results for 118 models. Other models could not be compared because either they have two adjacent state slots with -1 values (and do not support that, e.g. gear.1, in which case LTSmin aborts), the could not be translated (e.g. train-gate.1) or they are too large for all tools to finish (e.g. public_subscribe.5).

Total Time (s)
Average B/state

SPIN (collapse)

SPIN (bfspar)

DMC/    773

Totals for 118 models from the BEEM database.

The results for the comparison to LTSmin and SPIN using models from the BEEM database are shown in <ref>. When looking at <ref>, we can see that the combination of with is nearing the performance of LTSmin with time-wise. There are a number of outliers in both directions, but in the more time-consuming models one can see that LTSmin with has the edge over DMC, still. In terms of bytes per state (<ref>) we see a similar result. This shows in the total runtime and B/states as well: LTSmin/takes 20% less time and 23% less memory on average.

The comparison to follows the same trend time-wise (<ref>). clearly outperforms in terms of bytes per state (<ref>). The use of a Cleary table for the root set is a clear winner.

<Ref> compares DMC/dtree to SPIN with parallel BFS. Here, we see that DMC/dtree has the edge over SPIN, time-wise. <Ref> compares DMC/dtree to SPIN with , but without parallel BFS (these cannot be used simultaneously). Even then, SPIN with is outperformed by DMC/dtree in all benchmarks.

Comparing the storage components of DMC between themselves (<ref>), we notice all variants compress equally and also equal to LTSmin/. This is expected, as for fixed-length states, the implementations should be equivalent.
Furthermore, we see that loses in terms of bytes per state. The different structure of compression tree is again a likely cause: using the chain of perfectly balanced trees on average results in a less balanced tree than and thus on average a change (via insert() or delta()) requires more new nodes.
This is the price to pay to support variable-length states.

§ DISCUSSION AND CONCLUSION

We presented , a concurrent variable-length state storage component that stores states in a compression tree. It allows partial reconstruction of states and incremental updates to parts of states without reconstructing the entire states. To showcase , we implemented the DMC Model Checker
such that we can expose the functionality of through the .
We implemented three variable-length state vector models that make full use of the .

We compare to other state storage components, such as . In the case of , we pad all vectors to the length of the largest vector (manually determined). Here, is able to perform up to 2.9 times faster than using 29% less memory on average, without the a priori need to manually determine the largest vector. The advantage of increases as the difference in length between states increases.

To evaluate how much efficiency is lost to gain variable-length support, we compare using to LTSmin and SPIN using models from the BEEM database with fixed-length states. In this setting, is actually marginally faster than , but uses 30% more memory.
as a whole can only approach the performance of fixed-length state model checker LTSmin, sacrificing 20–23%. It is still outclassed by in terms of bytes per state, which compressed up to 2.3x more than . In the same tests, with is faster and provides a higher compression than SPIN with .
Compared to SPIN using , is 8.1x faster and uses 6.1x less memory.
Even though nor are particularly optimized for fixed-length complete-state changes, the performance overall is reasonable.

This research is part of the ongoing research towards creating a software model checker for use in a continuous integration pipeline. With the advent of [The source code of can be found at <https://github.com/bergfi/dtree>.], we are now one step closer.

§.§ Future work

We have seen that outperforms for fixed-length states. uses a Cleary Table for its root set. To improve , we can investigate if we can leverage such a table, since as is uses a 64-bit index and the DMC API currently has only 40 bits available for the state index.

DMC still needs to improve as well. For example, the search core of DMC is a simple parallel BFS, lacking a more sophisticated work-stealing algorithm. Implementing such a feature would improve the performance of DMC as a whole.

We aim to use DMC as the core for our upcoming multi-core software model checker LLMC. The purpose of LLMC is to model check LLVM IR assembly code. To model the stack and heap, lends itself perfectly. We can then compare the resulting implementation with .

§.§ Acknowledgements

The author would like to thank Arnd Hartmanns and Jaco van de Pol for their invaluable contributions and Alfons Laarman for discussions on and .
This research is sponsored by 3TU Big Software on the Run project (<http://www.3tu-bsr.nl/>).


40 1pt ARC Centre of Excellence for Transformative Meta-Optical Systems
(TMOS), Department of Electronic Materials Engineering, Research School of
Physics, The Australian National University, Canberra, ACT 2600, Australia
Fraunhofer Institute for Applied Optics and Precision Engineering, 07745 Jena,
Germany

# Robust Classical and Quantum Polarimetry with a Single Nanostructured
Metagrating

Shaun Lung<EMAIL_ADDRESS>Abbe Center of Photonics, Friedrich-
Schiller Universität, Albert-Einstein-Straße 15, 07745 Jena, Germany  Kai Wang
Department of Physics, McGill University, 3600 rue University, Montreal,
Quebec H3A 2T8, Canada  Nicolas R.H. Pedersen  Frank Setzpfandt  Abbe Center
of Photonics, Friedrich-Schiller Universität, Albert-Einstein-Straße 15, 07745
Jena, Germany  Andrey A. Sukhorukov  ARC Centre of Excellence for
Transformative Meta-Optical Systems (TMOS), Department of Electronic Materials
Engineering, Research School of Physics, The Australian National University,
Canberra, ACT 2600, Australia

###### Abstract

We formulate a new conceptual approach for one-shot complete polarization
state measurement with nanostructured metasurfaces applicable to classical
light and multi-photon quantum states, by drawing on the principles of
generalized quantum measurements based on positive operator-valued measures
(POVMs). Accurate polarization reconstruction from a combination of photon
counts or correlations from several diffraction orders is robust with respect
to even strong fabrication inaccuracies, requiring only a single classical
calibration of metasurface transmission. Furthermore, this approach operates
with a single metagrating without interleaving, allowing for the metasurface
size reduction while preserving the high transmission efficiency and output
beam quality. We theoretically obtained original metasurface designs,
fabricated the metasurface from amorphous silicon nanostructures deposited on
glass, and experimentally confirmed accurate polarization reconstruction for
laser beams. We also anticipate robust operation under changes in
environmental conditions, opening new possibilities for space-based imaging
and satellite optics.

###### keywords:

Metamaterials; Nanophotonics and photonic crystals; Polarization-selective
devices; Polarimetry

## Introduction

Single-shot optical polarimetry using ultra-thin nanostructured metasurfaces
opens up new opportunities for diverse applications 1, 2, 3, facilitating the
measurement of both classical 4, 5, 6, 7, 8 and quantum 9 polarization states.
Whereas traditional polarimetry involves multiple measurements performed via
physically varying bulk optical elements such as waveplates and polarizers 10,
single-shot approaches remove the need for reconfigurability, thereby avoiding
the associated measurement errors and facilitating real-time polarization
state monitoring 11, 12 that can be combined with spectral imaging 13, 14, 15,
16.

Similar to conventional classical 10 and quantum polarimetry 17, metasurfaces
were initially designed to perform several complimentary projection
measurements 18, such that each of the outputs corresponds to a specific
polarization 4, 1, 2, 3. This functionality can be accomplished with
metasurfaces that split particular polarization components into distinct
diffraction orders or focal spots. For the tomographic characterization of
quantum states, metasurfaces were previously developed by interleaving several
metagratings whose number scales linearly with the number of photons 9.
However, the simple interleaving (i) limits the device compactness since the
photons have to spatially overlap with multiple gratings at once and (ii)
introduces output beam distortions which reduce the detection efficiency. This
poses a question on how to overcome such limitations and thereby support the
emerging applications of metasurfaces in quantum imaging 19, 20, 21.

We reveal that efficient polarimetry with metasurfaces can be accomplished
without the commonly considered requirement of realizing close-to-perfect
polarization projection measurements. Remarkably, even if each of the
metasurface output ports represents a partial polarizer operation that by
itself provides inconclusive information on the input state, a tailored
combination of all outputs allows for very accurate polarization
reconstruction. We achieve this by adopting the framework of generalized
quantum measurements based on positive operator-valued measures (POVMs) 22,
23, 24, 25, 26, 27, 28 for the metasurface design. Such an approach
fundamentally improves the robustness with respect to nanofabrication
inaccuracies and also extends the flexibility in metasurface designs, allowing
in particular for higher efficiency and output beam shaping with small-area
metasurfaces. We present simulation results for one- and two-photon states.
Then, we experimentally demonstrate the operation with laser light that
illustrates the classical regime and also emulates a single-photon quantum
case 29.

Figure 1:  (a) Conceptual sketch of single-metagrating polarimetry.
Highlighted in yellow is a single 16 resonator unit cell. From the output
diffraction orders, any unknown input state may be reconstructed by knowing
the metasurface instrument matrix. (b) Poincaré sphere representation of two
arbitrarily chosen input states $|\Psi_{in}^{(1)}\rangle$,
$|\Psi_{in}^{(2)}\rangle$ characterized in Jones formalism by the polarization
angles and phases $(1.5,0.1)$, $(0.1,1.5)$ respectively. (c,d) Intensities of
output diffraction orders corresponding to the respective input states
$|\Psi_{in}^{(1)}\rangle$, $|\Psi_{in}^{(2)}\rangle$.

## Theory of multi-photon polarization measurements

We first formulate a general theory of polarization measurements. We consider
a metasurface, which splits an input beam into several diffraction orders as
sketched in Fig. 1(a). Then, the measurements of photon correlations between
the output ports can be used to reconstruct the input quantum state. In a
previous study 9, it was assumed that the photon polarization is specifically
selected at each diffraction order. This effectively required the metasurface
to act as a near-perfect polarization splitter, which came at a cost of larger
device area requirements due to a need to interleave several gratings and the
distortion of the output beam profile. In the following, we show that such
limitations can be removed, allowing for a compact and robust metasurface
design.

Since we consider a linear regime, the transformation of quantum states can be
expressed through the classical Jones transfer matrices $\mathbf{T}^{(m)}$ to
each of the diffraction orders number $m$,

$\psi_{out}^{(m)}=\mathbf{T}^{(m)}\psi_{in}\,,$ (1)

where $\psi_{in}$ and $\psi_{out}^{(m)}$ are the classical input and output
polarization states in the form of Jones vectors at the respective diffraction
orders. Each diffraction order functions, in general, as a partial polarizer
acting on the incoming light. It is then convenient for the following analysis
to perform a singular value decomposition (SVD) of the transfer matrices,

$\mathbf{T}^{(m)}=\textbf{U}^{(m)}\,\begin{pmatrix}\zeta_{m,1}&0\\\
0&\zeta_{m,2}\end{pmatrix}\left(\textbf{W}^{(m)}\right)^{\dagger}\,,$ (2)

where $\zeta_{m,1},\zeta_{m,2}\geq 0$ are the singular values,
$\textbf{U}^{(m)}=[\textbf{U}^{(m)}_{1},\textbf{U}^{(m)}_{2}]$ and
$\textbf{W}^{(m)}=[\textbf{W}^{(m)}_{1},\textbf{W}^{(m)}_{2}]$ are unitary
matrices, the subscripts indexing the respective columns. We choose the order
of singular values such that $\zeta_{m,1}\geq\zeta_{m,2}$. If $\zeta_{m,2}=0$,
the polarization state at the corresponding diffraction order is fixed as for
a case of a perfect polarizer. However, in a general case of $\zeta_{m,2}>0$,
the metasurface effectively acts as a partial polarizer, with a power
extinction ratio of $(\zeta_{m,2}/\zeta_{m,1})^{2}$.

We now formulate the transformation by the metasurface of the photon creation
and annihilation operators from the input, $\hat{a}_{p}^{\dagger}$ and
$\hat{a}_{p}$, where $p=\\{H,V\\}$ is the polarization state, to each of the
output diffraction orders, $\hat{b}_{m,p}^{\dagger}$ and $\hat{b}_{m,p}$. This
can be expressed through the linear transfer matrix elements as:

$\hat{b}_{m,p}=\sum_{p^{\prime}=H,V}\left(\mathbf{T}^{(m)}_{p,p^{\prime}}\right)^{\ast}\,\hat{a}_{p^{\prime}}\;,\quad\hat{b}_{m,p}^{\dagger}=\sum_{p^{\prime}=H,V}\mathbf{T}^{(m)}_{p,p^{\prime}}\,\hat{a}_{p^{\prime}}^{\dagger}\;.$
(3)

We target the tomography of quantum polarization-entangled states with a fixed
photon number $N$ at the input, which is a common practical task 17, 30, 31,
32, considering the simplest type of click detectors that cannot resolve the
number of arriving photons and cannot distinguish the photon polarization
state. Notably, the state characterization can also be generalized to the
regime when the maximum photon number is known using the approach of Ref. 33.
If no more than one photon arrives at such a detector positioned at the
diffraction order $m$, then its response is governed by the following POVM
operator 34,

$\hat{A}^{(m)}=\sum_{p=H,V}\hat{b}_{m,p}^{\dagger}\hat{b}_{m,p}=\sum_{p,p^{\prime}=H,V}\textbf{A}^{(m)}_{p,p^{\prime}}\,\hat{a}_{p}^{\dagger}\hat{a}_{p^{\prime}}\;,$
(4)

where using Eqs. (2) and (3) we calculate the matrix expression

$\mathbf{A}^{(m)}=\left(\mathbf{T}^{(m)}\right)^{\dagger}\mathbf{T}^{(m)}=\sum_{j=1,2}\zeta_{m,j}^{2}\left|\textbf{W}^{(m)}_{j}\right\rangle\left\langle\textbf{W}^{(m)}_{j}\right|\,=\left(\zeta_{m,1}^{2}-\zeta_{m,2}^{2}\right)\left|\textbf{W}^{(m)}_{1}\right\rangle\left\langle\textbf{W}^{(m)}_{1}\right|+\zeta_{m,2}^{2}\textbf{I}\,.$
(5)

We see that this is a sum of a polarization projection operator and a
polarization-insensitive detection. The presence of the latter term is a
consequence of the partial-polarizer transformation at each of the diffraction
orders. Although conventional polarimetry requires near-perfect polarizers
(i.e. $\zeta_{m,2}=0$), we apply the POVM formalism that enables unique and
accurate quantum state reconstruction in the regime of $\zeta_{m,2}>0$. Beyond
metasurfaces, we expect that the formulated approach can be also used to
enhance polarization measurements using nanowire detectors 35, overcoming the
limitations due to relatively low polarization extinction ratios 36.

After determining the detection operators, we find the probabilities of the
simultaneous detection of $N$ photons by a combination of $N$ detectors at the
diffraction orders $m_{1},m_{2},\ldots,m_{N}$, when there is exactly one
photon at each detector. Notably, if more than one photon arrives at a
particular detector, then the total number of coincidences measured across all
detectors will be less than $N$, and thus excluded from the analysis. The
$N$-detector correlations can be calculated as

$\Gamma(m_{1},m_{2},\ldots,m_{N})={\rm
Tr}\left(\rho^{(N)}\hat{A}_{m_{1}}\hat{A}_{m_{2}}\cdots\hat{A}_{m_{N}}\right).$
(6)

Here, $\rho^{(N)}$ is an input density matrix. Then, we follow an established
procedure 9, 37 to enumerate with index $q$ all the possible $N$ combinations
of $M$ detectors, $(m_{1},m_{2},\ldots,m_{N})$, and rewrite Eq. (6) in an
equivalent form,

$\Gamma_{q}=\sum_{s=1}^{S}{\bf B}_{q,s}r_{s}\;.$ (7)

Here $r_{s}$ are the independent real and imaginary parts of the input density
matrix defined according to the procedure in Ref. 32, $S=(N+3)!/(3!N!)$,
$q=1,\ldots,Q$, $Q=M!/(N!(M-N)!)$, and $M$ is the total number of detected
diffraction orders. The matrix elements ${\bf B}_{p,s}$ depend on the transfer
matrix elements, and more specifically on the vectors $\textbf{W}^{(m)}_{j}$
and singular values $\zeta_{m,p}$ according to the form of Eq. (5).

We can then reconstruct an input state from the correlation measurements by
performing a pseudo-inversion of Eq. (7), provided the number of different
correlations matches or exceeds the number of unknowns, $Q\geq S$, which is
satisfied when

$M\geq N+3\;.$ (8)

Importantly, in addition to the necessary condition in Eq. (8), it is
essential that reconstruction results are robust in presence of experimental
errors in the correlation measurements 38. This can be expressed as a
requirement to minimize the condition number $\kappa$ of matrix ${\bf B}$,
defined as a ratio of its largest and smallest singular values 39, 40, 9. The
condition number of an equation characterizes the worst-case error in the
output for a given error in the set of input parameters. In particular, the
condition number of a linear Eq. (7) characterizes the accuracy of the
reconstruction against errors in the observations.

## Experimental methods

### Optimized metagrating design

Figure 2:  Design of metagratings for (a,b) one-photon or classical and (c,d)
two-photon polarimetry. (a,c) Bottom - rendering of metagratings made of
amorphous silicon on glass substrate, illustrating different shapes and
orientations of rectangular nanopillars. Top - the phase retardances along the
two principal axes of nanopillars shown with blue and orange colors. (b,d)
Corresponding Poincaré sphere representation of partial polarizer
transformations at different diffraction orders as labeled.

To implement the multiple polarization transformations required to implement
the proposed polarimetry scheme, we leverage the flexibility of designing
metasurfaces as follows. First, we defined a metagrating consisting of $L$
rectangular nanoresonators arranged into a periodic supercell, such as shown
in Figs. 2(a,c). A transfer matrix of each nanoresonator can be represented in
the Jones formalism 41 as

$\mathbf{U}^{(l)}=\mathrm{\mathbf{R}}\left(-\theta^{(l)}\right)\cdot\begin{bmatrix}e^{i\phi_{o}^{(l)}}&0\\\
0&e^{i\phi_{e}^{(l)}}\end{bmatrix}\cdot\mathrm{\mathbf{R}}\left(\theta^{(l)}\right)\,,$
(9)

where $l$ is the resonator number, $\phi_{o}^{(l)}$ and $\phi_{e}^{(l)}$ are
the phase shifts imposed by the resonator along the ordinary and extraordinary
axes, and $\mathrm{\mathbf{R}}\left(\theta^{(l)}\right)$ is a two-by-two
rotation matrix by angle $\theta^{(l)}$. Then, we calculate the transfer
matrices of the metasurface at the output diffraction orders using the Fourier
transform,

$\mathbf{T}^{(m)}=\sum_{l=0}^{L}\mathbf{U}^{(l)}e^{-\frac{2\pi i}{N}ml}\,.$
(10)

We calculate these matrices numerically and then use Eqs. (4)-(7) to determine
the corresponding condition number $\kappa$, defined as the ratio of the
maximum and minimum singular value decomposition values of matrix ${\bf B}$
appearing in Eq. (7). By minimizing the condition number via a gradient-
descent algorithm, it is thus possible to obtain a highly robust design. This
optimization produces a set of phase parameters corresponding to each pixel of
the unit cell, consisting of phase shifts along the ordinary and extraordinary
axes as well as the orientation angles.

We present the optimized metagrating designs for one- and two-photon cases in
Figs. 2(a,c) for the minimum required number of outputs according to Eq. (8),
$M=4$ and $M=5$, respectively. The numerical values of the nanoresonator
parameters are presented in Supplementary Information Sec. S1. In these
examples, we chose to consider an array of $L=16$ nanoresonators, but it
should be observed that this choice is neither unique nor constrained, and
different numbers can be selected depending on the required angular diversion
of the diffraction orders after the metasurface.

As implied by Eq. (5), the vectors $\mathbf{W}_{j}^{(m)}$ are the basis states
of the polarization measurement. After converting them to a Stokes basis, we
plot them on a Poincaré sphere for convenient visualization of the basis
states at different diffraction orders for a given metasurface design, as
shown in Figs. 2(b,d). We define the vector lengths as

$R=1-\frac{\zeta_{m,2}^{2}}{\zeta_{m,1}^{2}},$ (11)

where $\zeta_{m,2}^{2}$ and $\zeta_{m,1}^{2}$ are by definition proportional
to the minimum and maximum powers transmitted to the relevant diffraction
order $m$ across a set of all possible polarization states under the same
input power. We note that $R=1$ corresponds to a fully polarized output, as
would be required for projective measurements. We observe that $R<1$ for
several outputs of our optimized metagratings, meaning that they act as
partial polarizers with the finite extinction ratio $(1-R)^{-1}$.
Nevertheless, this allows for efficient reconstruction of the input
polarization states.

Indeed, the corresponding inverse condition numbers, calculated for the
metagratings that are numerically optimized to maximize their values, are
around $0.37$ and $0.21$ in the spectral region of interest for single- and
two-photon states. The single-photon value is comparable to the theoretically-
best maximum 18 of $1/\sqrt{3}\simeq 0.58$. While the best possible value for
two-photon state reconstruction with click detectors, which do not resolve the
events of both photons exiting from the same port, is not known analytically,
previous numerical optimisations of integrated waveguide circuits 30, 37
reported values up to $\sim 0.25$, again close to the two-photon result above
in our current study.

### Robustness to fabrication errors

Figure 3:  Simulated effects of random phase errors, up to a maximum of
$\Delta\phi$, applied to each of the nanoresonators in the metagratings
optimized for one-photon and two-photon polarimetry as illustrated in Figs.
2(a) and (c), respectively. Shown are the probability densities, normalized to
a maximum of unity at each $\Delta\phi$, for the inverse condition numbers
$1/\kappa$ and diffraction efficiencies $\eta_{min}$. White lines indicate the
average values.

Based on the POVM formulation of the design, strong robustness against
fabrication errors is expected of the metasurface, since accurate
reconstruction is possible even if the polarization extinction ratios at
individual outputs are affected. We have carried out numerical simulations to
demonstrate this, considering as representative examples the metasurface
designs for single- and two-photon polarimetry presented above.

Under realistic fabrication scenarios, the most common deviations from the
analytical design pertain to the overall sizes of the nanoresonators. These
were modelled as variance in the phase shifts along the ordinary and
extraordinary axes of the individual nanopixels, corresponding to the
decomposition shown in Eq. (9), and thus, random errors up to $\Delta\phi$
were added to each nanopixel, and then the transmission of the altered
metasurface structure was calculated via Eq. (10). From this result, we
computed for the perturbed metasurface the overall inverse condition numbers
$1/\kappa$ and diffraction efficiency $\eta_{min}$.

The diffraction efficiency is defined as the minimum fraction of the total
input power that is diffracted to the chosen orders used in computing the
corresponding inverse condition number, specifically to the ($\pm 1$, $\pm 2$)
orders in the one- and (0, $\pm 1$, $\pm 2$) orders in the two-photon case. By
repeating such simulations over $10^{5}$ times for each value of $\Delta\phi$,
we estimated the degradation of performance vs. the degree of random errors,
as summarised in Fig. 3. The errors trialed range up to an extreme case of
$\Delta\phi=\pi/2$, which far exceeds the possible effects of the
nanofabrication errors. We also note that in the two-photon case, the zero
order diffraction was included in the computation of both $1/\kappa$ and
$\eta_{min}$, whereas the single-photon case excludes the zero order. The
additional basis state provided by the zero order diffraction was found to
slightly increase the achievable diffraction efficiency in the two-photon
case.

We find that for realistic levels of errors, below $10^{-1}$, both the inverse
condition number and the diffraction efficiencies do not drop significantly,
as shown in Fig. 3. However, as one might expect, the two-photon design is
more sensitive to errors than the single-photon case, owing to the necessary
consideration of additional diffraction orders to fully resolve multiple
photons.

## Experimental results

### Metasurface patterns: design and fabrication

We first determine the physical dimensions of nanoresonators realizing the
required phase delays, such as those summarized in Supplementary Information
Sec. S1, using an established procedure 41, 9. Specifically, we perform a
sweep of the length and width parameters of cuboidal pixels, while keeping
fixed the meta-grating period, height, and refractive indices, and calculate
the transmission using a numerical technique known as Rigorous Coupled Wave
Analysis (RCWA) 42, 43, 44. Thereby, we produce a lookup table from which
suitable designs could be simply selected. Thereby, one can design a
metasurface for operating at the desired spectral regions.

As a demonstration, we develop a metasurface for operation at the
telecommunication band, around the 1550 nm wavelength. As a platform, we
choose a dielectric metasurface made of amorphous silicon on a glass
substrate, leveraging the high transmissivity to ensure high efficiency of
operation. The height and period used were $832nm$ and $800nm$ respectively,
the former corresponding to measurements of deposited silicon that would be
used to fabricate the metasurface. The physical parameters of the metasurface
were thus designed by selecting pixels that fit desired phase parameters, up
to an arbitrary global phase. The combined metasurface structure designed was
then simulated using a commercial electrodynamics solver, CST Studio, as a
final optimization pass and to check the influence of optical couplings
between adjacent pixels that are not accounted for in the design of individual
nanoresonators. This optimization pass did not alter the design singificantly,
changing the design parameters by $<2\%$.

The metasurfaces were fabricated from a $832$ nm-thick amorphous silicon layer
prepared at the ANU node of Australian Nanofabrication Facility (ANFF) using
Plasma-Enhanced Chemical Vapor Deposition (PECVD) on a glass substrate. It was
subsequently etched at the University of Jena using Electron Beam Lithography
(EBL) and Inductively Coupled Plasma (ICP) etching. Slight variants in the
metasurfaces were prepared by intentionally varying the EBL exposure times,
allowing for a selection of best-case fabrication outcomes.

## Results and discussion

Figure 4:  (a) Schematic of the experimental setup used to classically
characterize the metasurface. Input states were prepared from a variable-
wavelength infrared laser using a fixed polarizer and motorized half- and
quarter- waveplates, before being collimated on the metasurface by lenses.
Output diffraction order intensities were collected using a CCD camera. (b)
Scanning Electron Microscope (SEM) image of the metasurface, fabricated as
$832nm$ amorphous silicon on glass via electron beam lithography. (c)
Representative reading as taken using the camera, and the processed intensity
plot obtained by slicing across the diffraction orders. Intensities were
normalized to 1.

After fabrication, we measured the transmission of classical light through the
metasurface. This enables the characterization of the metasurface
transformation from the input to distinct outputs, which can then be used for
the reconstruction of multi-photon states. We use a classical-quantum analogy
29 since both the single-photon (for $N=1$) counts at the outputs and the
classical output intensities can be described by Eqs. (6),(7). Here the
density matrix can be conveniently related to the Stokes parameters $\vec{S}$
45, 17, 46,

$\rho^{(1)}=\frac{1}{2}\sum_{j=0}^{3}S_{j}{\bf\sigma}_{j}\,=\frac{1}{2}\left(\begin{array}[]{ll}S_{0}+S_{3}&S_{1}-iS_{2}\\\
S_{1}+iS_{2}&S_{0}-S_{3}\end{array}\right),$ (12)

where ${\bf\sigma}_{j}$ are the Pauli matrices. For convenience, we choose the
independent parameters in the input density matrix as $\vec{r}=\vec{S}$. Then,
the matrix ${\bf B}$ in Eq. (7) has the meaning of an instrument matrix, since

$\Gamma=\mathbf{B}\,\vec{S}$ (13)

determines the intensities at the output diffraction orders for classical
light with the input Stokes parameters $\mathbf{S}$. Specifically, at the
output number $m$

$\Gamma_{m}=\sum_{j}\mathbf{B}_{m,j}\,S_{j}.$ (14)

On the other hand, using Eq. (6) at $N=1$ the output intensities can be
expressed as

$\mathbf{\Gamma}_{m}=\sum_{j,l}\mathbf{A}^{(m)}_{j,l}\mathbf{\rho}^{(1)}_{l,j}\,.$
(15)

We then obtain

$\mathbf{A}^{(m)}=\left(\begin{array}[]{ll}B_{m,1}+B_{m,4}&B_{m,2}-iB_{m,3}\\\
B_{m,2}+iB_{m,3}&B_{m,1}-B_{m,4}\end{array}\right).$ (16)

In this way, the classical instrument matrix $\mathbf{B}$ can be used to
determine all matrices $\mathbf{A}^{(m)}$ and thereby allow the subsequent
reconstruction of arbitrary multi-photon states.

The experimental setup is schematically illustrated in Fig. 4(a). Using a
half-wave plate, quarter-wave plate, and a fixed polarizer, polarization
states were prepared from a variable-wavelength laser operating in the
$1500-1575$ nm telecommunications bandwidth. The prepared polarization state
was then focused to a spotsize of approximately $15\,\mu$m normally incident
on the metasurface. The diffraction orders were then collected using an
objective lens with a high numerical aperture, and imaged onto an infrared CCD
camera using a convex lens. As a separate measurement, the camera was replaced
with a calibrated power meter in order to determine the total power that was
transmitted through the metasurface.

We measured the intensities of the diffraction orders over varying input
polarization orders by capturing images on a camera, such as shown in Fig.
4(c). First, the diffraction order spot locations were determined, and then
the individual intensities were extracted by integrating over the area of each
diffraction spot.

We performed measurements for a calibration set of 360 distinct input
polarization states and calculated the instrument matrix $\mathbf{B}$ for
diffraction orders $(\pm 2,\pm 1)$ by fitting its parameters according to Eq.
(13) using the known input states and measured output powers. We note that
each diffraction output intensity is defined by a row of the matrix
$\mathbf{B}$ that contains four elements, and therefore four or more
calibration measurements are needed to determine uniquely all the $\mathbf{B}$
elements. We intentionally use a large size of the calibration set to reduce
the effect of random noise in individual measurements through the averaging.

Figure 5:  (a) Poincaré sphere representation of the basis states as
calculated from the $(\pm 2,\pm 1)$ diffraction orders of the fabricated
metasurface. These are, as predicted, partially polarized states which deviate
from the original numerical design, shown as correspondingly coloured dots on
the Poincaré sphere. (b) The experimentally characterized inverse condition
numbers of the metasurface across a wavelength range are plotted using the
solid green line. The blue, dashed line represents the numerical result
through CST Studio simulations of the optimal metasurface design, and the
orange, dotted line represents the theoretical maximum. (c) The minimum (blue)
and maximum (red) power directed to each diffraction order, as determined by
experiment. Diffraction orders $\pm 3$ were measured but not utilized in the
calculations, and are shown for comparison. (d) Poincaré sphere showing a
comparison of input states versus the reconstructed states using the
experimentally characterized metasurface.

We compute the polarization bases of the measured instrument matrix and plot a
representative example in Fig. 5(a), corresponding to the wavelength of $1560$
nm. Here, we see that the metasurface operation has deviated from the designed
values as shown in Fig. 2(a), which may be ascribed to fabrication errors.
Nevertheless, these do not prevent the polarization characterization as we
discuss below.

We use the measured instrument matrix to compute the inverse condition number
dependence on the wavelengths, see Fig. 5(b). Importantly, it is only slightly
lower (by up to 15%) than the numerical design values indicated with the
dashed line, across the whole target operating bandwidth. For comparison, the
dotted line marks the theoretical maximum of $1/\sqrt{3}\simeq 0.58$ 18. Such
experimentally achieved performance is close to the value of $1/2.08\simeq
0.48$ for an interleaved metasurface 9, yet providing a fundamentally better
beam quality without distortions in the vertical direction.

We use the experimental data to visualize the minimum and maximum powers,
across all possible input polarization states and normalized to the input
power, at the selected output ports in Fig. 5(c). We note that the ratios of
these minimum to maximum powers define the polarization extinction, which
appears to be of the order of $10^{-1}$. Despite such a low extinction, an
accurate state reconstruction is still possible using the developed approach
based on the framework of generalized quantum measurements.

We demonstrate the reconstruction of the input states from the measured powers
at the selected diffraction orders $(\pm 2,\pm 1)$ after the metasurface. From
experimental data, a random set of angles for the half- and quarter-waveplates
were chosen. From these angles, and taking into account that the first
polarizer was positioned vertically, we thus define the input polarization
states $\rho_{in}^{(1)}$ before incidence on the metasurface. We also
independently reconstruct the input states from the measured intensities
$\vec{P}_{out}$ and the experimentally characterized instrument matrix
$\mathbf{B}$, as:

$\vec{S}=\mathbf{B}^{-1}\cdot\vec{P}_{out},$ (17)

and then determine the input density matrix using Eq. (12). We plot the
prepared and reconstructed states on the same Poincaré sphere for comparison,
as shown in Fig. 5(d).

In this plot, we show 44 data points, of which only 4 overlap with the
calibration dataset, thereby providing a valid test of the reconstruction
accuracy for unknown input states. We quantify the difference between the
normalized input ($\vec{S}^{(in)}$) and reconstructed ($\vec{S}^{(rec)}$)
Stokes vectors as
$\delta=1-|\langle\vec{S}^{(in)}\rvert\vec{S}^{(rec)}\rangle|^{2}/(\langle\vec{S}^{(in)}\rvert\vec{S}^{(in)}\rangle\langle\vec{S}^{(rec)}\rvert\vec{S}^{(rec)}\rangle)$,
such that $\delta=0$ only if the vectors are the same, up to a total
intensity. Then, we find that for the states in Fig. 5(d), the error in
reconstruction is $\delta\leq 2.9\%$.

Whereas multi-photon quantum experiments are beyond the scope of the current
work, we use the experimentally characterized instrument matrix of the
metasurface to calculate an inverse condition number for two-photon
reconstruction and find that $\kappa^{-1}\simeq 0.15$ can be achieved at a
wavelength of $1520nm$ (see Supplementary Information Section S2 for a
wavelength dependence), close to the theoretically optimized value of $0.21$.

## Conclusions

We have presented a general approach for the complete measurement of
polarization in both quantum multi-photon states and classical light using
nanostructured metasurfaces. Our method nontrivially combines the principles
of single-shot polarimetry and generalized quantum measurements with POVM
formalism, enabling accurate polarization reconstruction even in the presence
of significant fabrication errors or environmental changes, by performing a
simple device calibration. Our experimental measurements demonstrate a
complete reconstruction of classical polarization states with a maximum error
of $2.9\%$, while maintaining a high optical beam quality after the
metasurface for efficient detection. We anticipate that our new concept will
facilitate diverse quantum and classical applications, from laboratories to
satellite imaging systems, benefiting from extremely compact meta-devices
providing real-time polarimetric measurements.

## Acknowledgments

We thank Jihua Zhang and Dragomir Neshev for the helpful discussions.

## Funding sources

This work was supported by the Australian Research Council (NI210100072,
CE200100010), US AOARD (19IOA053), the Thuringian Ministry for Economy,
Science, and Digital Society (2021 FGI 0043), the German Federal Ministry of
Education and Research (FKZ 13N14877), the German Research Foundation (IRTG
2675), the European Union through the ERASMUS+ program, the German Academic
Exchange Service (grant 57388353), and UA-DAAD exchange scheme. This work used
the ACT node of the NCRIS-enabled Australian National Fabrication Facility
(ANFF-ACT).

## Supporting information

Table S1; numerically optimized values of phase shifts corresponding to the
metasurface renders shown in Figs. 2(a,c), Fig. S1; figure plotting
experimentally calculated two-photon inverse condition number.

## References

  * Martinez 2018 Martinez, A. Polarimetry enabled by nanophotonics. _Science_ 2018, _362_ , 750–751, DOI: 10.1126/science.aau7494
  * Intaravanne and Chen 2020 Intaravanne, Y.; Chen, X. Z. Recent advances in optical metasurfaces for polarization detection and engineered polarization profiles. _Nanophotonics_ 2020, _9_ , 1003–1014, DOI: 10.1515/nanoph-2019-0479
  * Rubin et al. 2021 Rubin, N. A.; Shi, Z. J.; Capasso, F. Polarization in diffractive optics and metasurfaces. _Adv. Opt. Photon._ 2021, _13_ , 836–970, DOI: 10.1364/AOP.439986
  * Pors et al. 2015 Pors, A.; Nielsen, M. G.; Bozhevolnyi, S. I. Plasmonic metagratings for simultaneous determination of Stokes parameters. _Optica_ 2015, _2_ , 716–723, DOI: 10.1364/OPTICA.2.000716
  * Zhang et al. 2019 Zhang, X. Q.; Yang, S. M.; Yue, W. S.; Xu, Q.; Tian, C. X.; Zhang, X. X.; Plum, E.; Zhang, S.; Han, J. G.; Zhang, W. L. Direct polarization measurement using a multiplexed Pancharatnam-Berry metahologram. _Optica_ 2019, _6_ , 1190–1198, DOI: 10.1364/OPTICA.6.001190
  * Arbabi et al. 2018 Arbabi, E.; Kamali, S. M.; Arbabi, A.; Faraon, A. Full-Stokes Imaging Polarimetry Using Dielectric Metasurfaces. _ACS Photonics_ 2018, _5_ , 3132–3140, DOI: 10.1021/acsphotonics.8b00362
  * Wei et al. 2017 Wei, S. W.; Yang, Z. Y.; Zhao, M. Design of ultracompact polarimeters based on dielectric metasurfaces. _Opt. Lett._ 2017, _42_ , 1580–1583, DOI: 10.1364/OL.42.001580
  * Yang et al. 2018 Yang, Z. Y.; Wang, Z. K.; Wang, Y. X.; Feng, X.; Zhao, M.; Wan, Z. J.; Zhu, L. Q.; Liu, J.; Huang, Y.; Xia, J. S.; Wegener, M. Generalized Hartmann-Shack array of dielectric metalens sub-arrays for polarimetric beam profiling. _Nat. Commun._ 2018, _9_ , 4607, DOI: 10.1038/s41467-018-07056-6
  * Wang et al. 2018 Wang, K.; Titchener, J. G.; Kruk, S. S.; Xu, L.; Chung, H. P.; Parry, M.; Kravchenko, I. I.; Chen, Y. H.; Solntsev, A. S.; Kivshar, Y. S.; Neshev, D. N.; Sukhorukov, A. A. Quantum metasurface for multiphoton interference and state reconstruction. _Science_ 2018, _361_ , 1104–1107, DOI: 10.1126/science.aat8196
  * Chekhova and Banzer 2021 Chekhova, M.; Banzer, P. _Polarization of Light In Classical, Quantum, and Nonlinear Optics_ ; De Gruyter: Berlin, 2021; DOI: 10.1515/9783110668025
  * Rubin et al. 2019 Rubin, N. A.; D’Aversa, G.; Chevalier, P.; Shi, Z. J.; Chen, W. T.; Capasso, F. Matrix Fourier optics enables a compact full-Stokes polarization camera. _Science_ 2019, _365_ , eaax1839, DOI: 10.1126/science.aax1839
  * Rubin et al. 2022 Rubin, N. A.; Chevalier, P.; Juhl, M.; Tamagnone, M.; Chipman, R.; Capasso, F. Imaging polarimetry through metasurface polarization gratings. _Opt. Express_ 2022, _30_ , 9389–9412, DOI: 10.1364/OE.450941
  * Chen et al. 2016 Chen, W. T.; Torok, P.; Foreman, M. R.; Liao, C. Y.; Tsai, W. Y.; Wu, P. R.; Tsai, D. P. Integrated plasmonic metasurfaces for spectropolarimetry. _Nanotechnology_ 2016, _27_ , 224002, DOI: 10.1088/0957-4484/27/22/224002
  * Ding et al. 2017 Ding, F.; Pors, A.; Chen, Y. T.; Zenin, V. A.; Bozhevolnyi, S. I. Beam-Size-Invariant Spectropolarimeters Using Gap-Plasmon Metasurfaces. _ACS Photonics_ 2017, _4_ , 943–949, DOI: 10.1021/acsphotonics.6b01046
  * Sun et al. 2022 Sun, T.; Hu, J. P.; Zhu, X. J.; Xu, F.; Wang, C. H. Broadband Single-Chip Full Stokes Polarization-Spectral Imaging Based on All-Dielectric Spatial Multiplexing Metalens. _Laser Photon. Rev._ 2022, _16_ , 2100650, DOI: 10.1002/lpor.202100650
  * Camayd-Munoz et al. 2020 Camayd-Munoz, P.; Ballew, C.; Roberts, G.; Faraon, A. Multifunctional volumetric meta-optics for color and polarization image sensors. _Optica_ 2020, _7_ , 280–283, DOI: 10.1364/OPTICA.384228
  * James et al. 2001 James, D. F. V.; Kwiat, P. G.; Munro, W. J.; White, A. G. Measurement of qubits. _Phys. Rev. A_ 2001, _64_ , 052312, DOI: 10.1103/PhysRevA.64.052312
  * Foreman et al. 2015 Foreman, M. R.; Favaro, A.; Aiello, A. Optimal Frames for Polarization State Reconstruction. _Phys. Rev. Lett._ 2015, _115_ , 263901, DOI: 10.1103/PhysRevLett.115.263901
  * Pittman et al. 1995 Pittman, T. B.; Shih, Y. H.; Strekalov, D. V.; Sergienko, A. V. Optical imaging by means of 2-photon quantum entanglement. _Phys. Rev. A_ 1995, _52_ , R3429–R3432, DOI: 10.1103/PhysRevA.52.R3429
  * Altuzarra et al. 2019 Altuzarra, C.; Lyons, A.; Yuan, G. H.; Simpson, C.; Roger, T.; Ben-Benjamin, J. S.; Faccio, D. Imaging of polarization-sensitive metasurfaces with quantum entanglement. _Phys. Rev. A_ 2019, _99_ , 020101, DOI: 10.1103/PhysRevA.99.020101
  * Vega et al. 2021 Vega, A.; Pertsch, T.; Setzpfandt, F.; Sukhorukov, A. A. Metasurface-Assisted Quantum Ghost Discrimination of Polarization Objects. _Phys. Rev. Appl._ 2021, _21_ , 064032, DOI: 10.1103/PhysRevApplied.16.064032
  * Brandt 1999 Brandt, H. E. Positive operator valued measure in quantum information processing. _Am. J. Phys._ 1999, _67_ , 434–439, DOI: 10.1119/1.19280
  * Hamieh et al. 2004 Hamieh, S.; Kobes, R.; Zaraket, H. Positive-operator-valued measure optimization of classical correlations. _Phys. Rev. A_ 2004, _70_ , 052325, DOI: 10.1103/PhysRevA.70.052325
  * Renes et al. 2004 Renes, J. M.; Blume-Kohout, R.; Scott, A. J.; Caves, C. M. Symmetric informationally complete quantum measurements. _J. Math. Phys._ 2004, _45_ , 2171–2180, DOI: 10.1063/1.1737053
  * Ziman 2008 Ziman, M. Process positive-operator-valued measure: A mathematical framework for the description of process tomography experiments. _Phys. Rev. A_ 2008, _77_ , 062112, DOI: 10.1103/PhysRevA.77.062112
  * Shapiro 2008 Shapiro, J. H. Continuous positive-operator-valued measurement of photon polarization. _Phys. Rev. A_ 2008, _77_ , 052330, DOI: 10.1103/PhysRevA.77.052330
  * Al Khafaji et al. 2022 Al Khafaji, M. A.; Cisowski, C. M.; Jimbrown, H.; Croke, S.; Padua, S.; Franke-Arnold, S. Single-shot characterization of vector beams by generalized measurements. _Opt. Express_ 2022, _30_ , 458352, DOI: 10.1364/OE.458352
  * Martinez et al. 2023 Martinez, D.; Gomez, E. S.; Carine, J.; Pereira, L.; Delgado, A.; Walborn, S. P.; Tavakoli, A.; Lima, G. Certification of a non-projective qudit measurement using multiport beamsplitters. _Nat. Phys._ 2023, _19_ , 190–197, DOI: 10.1038/s41567-022-01845-z
  * Barnett 2022 Barnett, S. M. On single-photon and classical interference. _Phys. Scr._ 2022, _97_ , 114004, DOI: 10.1088/1402-4896/ac971a
  * Titchener et al. 2016 Titchener, J. G.; Solntsev, A. S.; Sukhorukov, A. A. Two-photon tomography using on-chip quantum walks. _Opt. Lett._ 2016, _41_ , 4079–4082, DOI: 10.1364/OL.41.004079
  * Oren et al. 2017 Oren, D.; Mutzafi, M.; Eldar, Y. C.; Segev, M. Quantum state tomography with a single measurement setup. _Optica_ 2017, _4_ , 993–999, DOI: 10.1364/OPTICA.4.000993
  * Titchener et al. 2018 Titchener, J. G.; Grafe, M.; Heilmann, R.; Solntsev, A. S.; Szameit, A.; Sukhorukov, A. A. Scalable on-chip quantum state tomography. _npj Quant. Inform._ 2018, _4_ , 19, DOI: 10.1038/s41534-018-0063-5
  * Bayraktar et al. 2016 Bayraktar, O.; Swillo, M.; Canalias, C.; Bjork, G. Quantum-polarization state tomography. _Phys. Rev. A_ 2016, _94_ , 020105, DOI: 10.1103/PhysRevA.94.020105
  * Kok and Lovett 2010 Kok, P.; Lovett, B. W. _Introduction to Optical Quantum Information Processing_ ; Cambridge University Press: Cambridge, 2010; DOI: 10.1017/CBO9781139193658
  * Zhu et al. 2021 Zhu, Y.; Raj, V.; Li, Z. Y.; Tan, H. H.; Jagadish, C.; Fu, L. Self-Powered InP Nanowire Photodetector for Single-Photon Level Detection at Room Temperature. _Adv. Mater._ 2021, _33_ , 2105729, DOI: 10.1002/adma.202105729
  * Park and Crozier 2015 Park, H.; Crozier, K. B. Elliptical silicon nanowire photodetectors for polarization-resolved imaging. _Opt. Express_ 2015, _23_ , 7209–7216, DOI: 10.1364/OE.23.007209
  * Wang et al. 2019 Wang, K.; Suchkov, S. V.; Titchener, J. G.; Szameit, A.; Sukhorukov, A. A. Inline detection and reconstruction of multiphoton quantum states. _Optica_ 2019, _6_ , 41–44, DOI: 10.1364/OPTICA.6.000041
  * Filippov and Man’ko 2010 Filippov, S. N.; Man’ko, V. I. Inverse spin-s portrait and representation of qudit states by single probability vectors. _J. Russ. Laser Res._ 2010, _31_ , 32–54, DOI: 10.1007/s10946-010-9122-x
  * Bogdanov et al. 2010 Bogdanov, Y. I.; Kulik, S. P.; Moreva, E. V.; Tikhonov, I. V.; Gavrichenko, A. K. Optimization of a quantum tomography protocol for polarization qubits. _JETP Lett._ 2010, _91_ , 686–692, DOI: 10.1134/S0021364010120143
  * Miranowicz et al. 2014 Miranowicz, A.; Bartkiewicz, K.; Perina, J.; Koashi, M.; Imoto, N.; Nori, F. Optimal two-qubit tomography based on local and global measurements: Maximal robustness against errors as described by condition numbers. _Phys. Rev. A_ 2014, _90_ , 062123–12, DOI: 10.1103/PhysRevA.90.062123
  * Arbabi et al. 2015 Arbabi, A.; Horie, Y.; Bagheri, M.; Faraon, A. Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission. _Nat. Nanotech._ 2015, _10_ , 937–U190, DOI: 10.1038/NNANO.2015.186
  * Lalanne 1997 Lalanne, P. Improved formulation of the coupled-wave method for two-dimensional gratings. _J. Opt. Soc. Am. A_ 1997, _14_ , 1592–1598, DOI: 10.1364/JOSAA.14.001592
  * Lalanne and Morris 1996 Lalanne, P.; Morris, G. M. Highly improved convergence of the coupled-wave method for TM polarization. _J. Opt. Soc. Am. A_ 1996, _13_ , 779–784, DOI: 10.1364/JOSAA.13.000779
  * Hugonin and Lalanne 2021 Hugonin, J. P.; Lalanne, P. Reticolo software for grating analysis. _arXiv preprint arXiv:2101.00901_ 2021, DOI: 10.48550/arXiv.2101.00901, (Date accessed: 2020-03-01)
  * Fano 1949 Fano, U. Remarks on the classical and quantum-mechanical treatment of partial polarization. _J. Opt. Soc. Am._ 1949, _39_ , 859–863, DOI: 10.1364/JOSA.39.000859
  * Goldberg et al. 2021 Goldberg, A. Z.; De La Hoz, P.; Bjork, G.; Klimov, A. B.; Grassl, M.; Leuchs, G.; Sanchez-Soro, L. L. Quantum concepts in optical polarization. _Adv. Opt. Photon._ 2021, _13_ , 1–73, DOI: 10.1364/AOP.404175

Supplementary Materials

##  Appendix A Optimized metasurface parameters

Table S1 contains the numerically optimized values of phase shifts along the
$x$ and $y$ directions as well as the angle of rotations of each
nanoresonator, corresponding to the metasurface renders shown in Figs. 2(a,c).

Resonator no. | One-photon polarimetry | Two-photon polarimetry
---|---|---
| $\phi_{x}$ | $\phi_{y}$ | $\theta$ | $\phi_{x}$ | $\phi_{y}$ | $\theta$
1 | -0.29 | -0.15 | 0.43 | -1.19 | -1.5 | 0.72
2 | -0.45 | 0.11 | -1.34 | -1.08 | 0.85 | 0.94
3 | -0.73 | 0.33 | -1.48 | -1.09 | 0.078 | -0.092
4 | -1.48 | -0.30 | 1.1 | -1.31 | -1.12 | -0.28
5 | -1.06 | -0.22 | 0.86 | -1.39 | 1.44 | 0.76
6 | 1.08 | -0.38 | 0.75 | 1.11 | 1.38 | 0.53
7 | -1.51 | -0.42 | 0.83 | 1.47 | -1.42 | -1.29
8 | 1.52 | -0.71 | 0.86 | -1.4 | 1.39 | 0.128
9 | 1.36 | -0.72 | 1.01 | -0.40 | 1.19 | -0.0016
10 | 1.54 | -0.26 | 1.0 | 0.7 | 1.22 | -0.362
11 | 1.18 | -1.13 | -1.25 | 1.06 | 0.072 | 0.26
12 | 1.05 | -0.48 | 1.52 | -0.77 | 1.11 | 1.42
13 | 0.29 | -0.41 | 1.56 | 1.2 | 1.52 | -0.29
14 | -0.76 | -0.34 | 1.37 | 1.46 | 1.34 | 1.17
15 | -1.22 | -0.41 | 1.02 | -1.55 | -1.07 | -0.72
16 | -1.29 | -0.33 | 0.77 | 0.288 | -1.25 | 0.66
Table S1: Phases and rotations in radians of nanoresonators forming
metasurfaces optimized for one- and two-photon polarimetry.

##  Appendix B 2-photon condition number

Figure S1:  (a) Poincaré sphere representation of the basis states as
calculated from the $(\pm 2,\pm 1,0)$ diffraction orders of the fabricated
metasurface. (b) The two-photon state reconstruction condition numbers vs. the
wavelength determined from the experimentally characterized metasurface
instrument matrix (green line) and from the finite-difference numerical
simulations with CST studio (blue line). (c) The minimum (blue) and maximum
(red) power directed to each diffraction order, as determined by experiment.