# a q-EW-TOPSIS model of grey correlation for supply capacity evaluation

Jia-Ming Liao  Yu-Jie Huang  Ke-Ming Shen<EMAIL_ADDRESS>School of Science, East China University of Technology, Nanchang, China

###### Abstract

The paper describes a new supply capacity evaluation model based on the non-
extensive statistical entropy. The traditional EW-TOPSIS model is selected as
baseline and the GRA method is used to modify it. The correction results in
the non-extensive parameter $q$ which leads to the so-called $q$-EW-TOPSIS
model. This new model has advantages over the traditional EW-TOPSIS model,
including the ability to accurately evaluate indicator weights with smaller
sample sizes and weaker rules, and a more stable and closer-to-complete
structure due to the use of entropy evaluation and mutual restriction between
indicators. This study provides a more reliable and universal modified EW
model. It is proved to be a more compatible model with systems and own greater
credibility.

###### pacs:

02.50.-r; 02.70.Rr; 05.10.-a; 89.75.-k

##  I Introduction

Being an advanced way of organization and management technology, modern
logistics has attracted more and more attentions of both enterprises and
governments. Nowadays, the methods of logistics capacity evaluation are well
studied in various areas ref1 ; ref2 ; ref3 ; ref4 . It is also becoming a
vital issue to quantitatively analyze the supply capacity, comprehensively
examine the advantages and disadvantages from different levels, and
comprehensively weigh the advantages and disadvantages, in order to provide a
scientific basis for the decision-making of the whole system.

The conventional evaluation method is to evaluate and make decisions through
AHP (Analytic Hierarchy Process) AHP ; AHP1 . This method conducts a
quantitative analysis of qualitative problems and makes the various factors in
complex problems orderly by dividing them into interrelated levels. The
importance of elements in one or two levels is then quantitatively described.
Finally, the weight factor, which reflects the relevant importance order of
elements in each level, is calculated mathematically. Relatively simple,
flexible, and practical, AHP has been widely used as soon as proposed.
However, under the indicators of similar importance, this subjective analysis
method cannot solve the defects of large accidental errors existing in
different multiple scoring. In addition, it is prone to errors when there are
many evaluation factors and a large scale, and the evaluation results are not
accurate and reliable enough. There are also many classical objective
evaluation models, such as the combined evaluation method cem  and data
envelopment method data , both of which do not meet the needs of researchers
because of their inaccuracy. It is often necessary to use a variety of
classical algorithms together, and then select the optimal model from them. In
addition, there are a lot of strong adaptive algorithms that can avoid the
subjectivity, such as the comprehensive evaluation method of BP (Back-
Propagation) neural network BP ; BP1 ; BP2 . It carries on the contribution
analysis of variables to eliminate insignificant factors, and gives out the
current best weight factors by considering the space changed over time. While
BP needs a lot of training samples and long training time as well. Moreover,
the complexity of this algorithm is quite high. On the other hand, for most
traditional algorithms, being a class of static models, it is also difficult
to maintain the quality of the model’s accuracy as data evolves. Furthermore,
when the data size is not sufficient and the impact of the loss of information
during the acquisition channel cannot be ignored, there is no corresponding
method to compensate for this loss of information, which mostly leads to large
changes in the evaluated weights with data increasing.

TOPSIS (Technique for Order Preference by Similarity to Ideal Solution) TOPSIS
; TOPSIS0  is an effective evaluation method based on multi-objective
decision-making. It normalizes the original data, calculates the weight of
limited evaluation objects by weighting, and evaluates the quality of the
evaluation objects. EW-TOPSIS method (Entropy weight TOPSIS method) EWM , as a
modified TOPSIS method with weights determined by information entropy, ensures
the accuracy of evaluation. However, considering the uncertainty of
information with data changing, the complexity of the model is moderate. It
somehow makes the determined weights change greatly, which leads to the fact
that stable weights should be reached around the EW-TOPSIS results rather than
themselves. Thus the data need to be perturbed correspondingly many times and
this model is also to be synthesized many times to find a metastable point,
which is then cumbersome.

In order to further optimize the TOPSIS model, we hereby introduce the non-
extensive statistical theory. It starts by generalizing the standard Maxwell-
Boltzmann statistics and covers lots of entropic forms such as the ones
proposed by Renyi Renyi , Tsallis TA , Abe Abe , Kaniadakis Kaniadakis  among
others. In this work, the Tsallis entropy is applied which obtains a non-
extensive parameter, $q$. When $q\to 1$ it recovers the classical theory.
Based on the non-extensive Tsallis entropy weight, a $q$-EW-TOPSIS model is
then applied into the investigations on the supply capacity evaluation.

##  II theoretical framework

We have created an algorithmic flowchart to facilitate reader comprehension,
as illustrated in Figure II-1.

Figure II-1: The algorithm flow chart of the q-EW TOPSIS model.

Subsequently, I will delve into the specific details of the algorithm.

In this paper, we utilize partially available data provided by the
Mathematical Modeling Suppliers and Customers of the Chinese Society of
Industrial and Applied Mathematics to verify the advantages of the model and
algorithm, as the model and algorithm do not strictly require specific data
size and form dataset .

###  II.1 EW-TOPSIS model

TOPSIS method was firstly proposed by C. L. Hwang and K. Yoon in 1981 ref6 .
It ranks a limited number of evaluation objects according to their proximity
to the idealized target, and evaluates the relative merits of the existed
objects. It is a multi-index comprehensive evaluation method, which can
synthesize multiple indicators into a single matrix for investigation. In our
model, the EW-TOPSIS model is applied to calculate out a reference vector,
which serves as a set of evaluation sequences to be evaluated; the scores are
calculated from original weights under ideal conditions, namely, the scores of
different suppliers in this work.

First of all, all the relevant parameters are specifically classified as
follows. The supply stability is defined as a minimum indicator which is also
called the cost-type indicator, namely the smaller value of it leads to the
better result. On the contrary, there are three maximum indicators, the supply
quantity (i.e., efficiency-type indicator), the supply continuity and the
ambiguous supply capacity. All of them should be forward normalized in
matrixes by using the range method next.

For the minimum indicators,

$x^{\prime}_{ij}=\frac{x_{\max}-x_{j}}{x_{\max}-x_{\min}},$ (1)

where $x^{\prime}_{ij}$ stands for the positively normalized indicators,
$x_{min(max)}$ for its minimum(maximum) value of $j$-items. Similarly, for
maximum indicators,

$x^{\prime}_{ij}=\frac{x_{j}-x_{\min}}{x_{\max}-x_{\min}}.$ (2)

Thus can we obtain the probability matrix as,

$p_{ij}=\frac{x^{\prime}_{ij}}{\sum_{i}^{m}x^{\prime}_{ij}}\quad,\quad 0\leq
p_{ij}\leq 1,$ (3)

with $m$ denoting the total number of microstates of system.

To better classify the EW-TOPSIS model, we firstly consider the classical
statistics. It is well known that Shannon entropy goes like,

$S_{Sh}[p_{i}]=E\left[\ln\frac{1}{p_{i}}\right]=-\sum_{i=1}^{m}p_{i}\ln
p_{i},$ (4)

where $E[x]$ represents the expected value of a variable $x$. Worthy to note
that hereby the Boltzmann constant $k_{B}=1$ is applied without losing
generality. The information entropy $e_{j}$ of the $j$th-parameter is then
defined as

$e_{j}=S_{j}[p_{ij}]=-\sum_{i=1}^{m}p_{ij}\ln\left(p_{ij}\right).$ (5)

The information utility value $d_{j}$ is calculated as follows:

$d_{j}=1-e_{j},$ (6)

together with the weight factor for each evaluation indicator:

$W_{S}(j)=\frac{d_{j}}{\sum_{j=1}^{n}d_{j}},$ (7)

where $n$ denotes the number of $j$th-evaluation indicator.

In order to eliminate the influence of dimensions of different indexes, the
matrix needs to be further standardized. Next, we normalize the forward matrix
X, with its element $z_{ij}$

$z_{ij}=\frac{x_{ij}}{\sqrt{\sum_{i=1}^{m}x_{ij}^{2}}}$ (8)

and $z^{*}_{ij}$

$z_{ij}^{*}=z_{ij}\cdot W_{S}(j).$ (9)

The maximum or optimized scheme is then defined as:

$\displaystyle Z^{*+}$
$\displaystyle=(Z_{1}^{*+},~{}Z_{2}^{*+},~{}\cdots,~{}Z_{m}^{*+})$
$\displaystyle=\big{(}max\\{z_{11},z_{21},\cdots,z_{n1}\\},~{}max\\{z_{12},z_{22},\cdots,z_{n2}\\},~{}\cdots,~{}max\\{z_{1m},z_{2m},\cdots,z_{nm}\\}\big{)}.$
(10)

Similarly we have the minimum or deteriorated scheme,

$\displaystyle Z^{*-}$
$\displaystyle=(Z_{1}^{*-},~{}Z_{2}^{*-},~{}\cdots,~{}Z_{m}^{*-})$
$\displaystyle=\big{(}min\\{z_{11},z_{21},\cdots,z_{n1}\\},~{}min\\{z_{12},z_{22},\cdots,z_{n2}\\},~{}\cdots,~{}min\\{z_{1m},z_{2m},\cdots,z_{nm}\\}\big{)}.$
(11)

Using them the distance between the $i$th-evaluation subject
($i=1,2,\cdots,n$) and the maximum (minimum) value is well clarified,

$D_{i}^{\pm}=\sqrt{\sum_{j=1}^{m}\left(Z_{j}^{*\pm}-z_{ij}^{*}\right)^{2}}.$
(12)

The normalized score $S_{i}$ for the $i$th subject is finally obtained as:

$S_{i}=\frac{D_{i}^{-}}{D_{i}^{+}+D_{i}^{-}}.$ (13)

It is obvious that $S_{i}$ belongs to [0,1], and larger values of $S_{i}$ lead
to the larger $D_{i}^{+}$, which gets closer to the maximum scheme, otherwise
the opposite. Thus can we obtain a set of reference vectors calculated using
the original EW-TOPSIS model as the scores.

###  II.2 Grey Relational Analysis

Before utilizing the non-extensive theory to generalize the classical EW-
TOPSIS model, it is crucial to thoroughly investigate the correlation degree
among the evaluation sequences and the reference scores, in order to
iteratively refine the data. The concept of grey systems was firstly proposed
by Professor Deng Jilong GRA , a control scientist and engineer, in comparison
to white systems and black systems. According to control theory, the color
usually represents the amount of information known about a system. White
systems represent systems with sufficient information, such as a mechanical
system, where the relationships between elements are determinable. On the
other hand, black systems represent systems where the structure is not clear,
often referred to as black boxes. Grey systems fall in between, representing
systems where only partial understanding is available. In our model, due to
the lack of information and the impact of the loss of information, the
correlation degree between the evaluation sequences and the reference scores
can vary, and thus the use of grey systems theory to correct the deviation is
appropriate.

The Grey Relational Analysis (GRA) is a very active branch of the grey systems
theory. Its basic principle is to determine the proximity of the relationship
between different sequences based on the geometric shape of the corresponding
sequence curves. Many researchers have already combined this method with other
algorithms, validating its effectiveness in supplier selection and evaluation
GRA1 ; GRA2 . Therefore, GRA will be applied in this study to calculate the
correlation degree among the evaluation sequences and the reference scores.

GRA is to calculate the degree in this work since in the grey system grading
curves of each parameter could be quantitatively analyzed in it. Nevertheless,
it is only adopted to obtain the correction amplitude to modify the idealized
EW-TOPSIS model because of the inevitable errors existed in the evaluation of
weight factors.

Firstly consider the characteristic sequence of system:

$R=\left[\begin{array}[]{ccc}x_{1}(1),x_{1}(2),\cdots,x_{1}(n),S_{1}\\\
x_{2}(1),x_{2}(2),\cdots,x_{2}(n),S_{2}\\\
\vdots\qquad\vdots\qquad\ddots\qquad\vdots\\\
x_{m}(1),x_{m}(2),\cdots,x_{m}(n),S_{m}\end{array}\right],$ (14)

where the subsequences $x_{i}(j),~{}(i=1,2,\cdots,m;j=1,2,\cdots,n)$
correspond to the values of $j$th-indicators with $i$th-ID. The score $S_{i}$
is following the definition as in Eq.(13).

The correlation coefficients of the subsequence are then calculated for each
sequence. The correlation coefficient on $S_{i}$ with $x_{i}$ at $j$ point is
defined as:

$\gamma\left(S_{i},x_{i}(j)\right)=\frac{\min_{i}\min_{j}\left|S_{i}-x_{i}(j)\right|+\xi\max_{i}\max_{j}\left|S_{i}-x_{i}(j)\right|}{\left|S_{i}-x_{i}(j)\right|+\xi\max_{i}\max_{j}\left|S_{i}-x_{i}(j)\right|}.$
(15)

Here $\min_{i}\min_{j}\left|S_{i}-x_{i}(j)\right|$ denotes the absolute value
of the secondary minimum difference and
$\max_{i}\max_{j}\left|S_{i}-x_{i}(j)\right|$ of the maximum one,
respectively. On the other hand, $\xi$ is introduced as the grey distinguish
coefficient with its range of $[0,1]$; it is shown better when $\xi<0.5468$.

Using Eq.(15) another coefficient between $X_{i}$ and $X$ could be written as:

$\gamma\left(S(j),X_{i}\right)=\frac{1}{m}\sum_{i=1}^{m}\gamma\left(S_{i},x_{i}(j)\right).$
(16)

Since the grey correlation coefficients are dimensionless and the contribution
of a group of vectors to the grey system could be quantitatively analyzed
under the same grey resolution coefficient, they are normalized to obtain the
contribution weight of each index to the systemref8 ,

$W_{r}(j)=\frac{\gamma_{j}}{\sum_{j=1}^{n}\gamma_{j}}.$ (17)

In the analysis of the supply ability score in this work, $m=402$ indicates
the number of suppliers and $n=4$ is for four indicators after the forward
transformation. The scoring matrix $S$ obtained by the EW-TOPSIS model is used
as a reference sequence of $S_{i}$ in this work and other indicators as
comparison sequences of $x_{i}(j)$. Furthermore, the correlation degree is
calculated through four indicators, i.e., the supply stability, quantity,
continuity and the ambiguous supply capacity. The weight curve then gives a
guide to how much the model should change (up or down).

###  II.3 $q$-EW-TOPSIS model

It is inevitable that there is a significant amount of noise in the process of
information acquisition, which deviates significantly from the ideal case.
When the data volume is not large enough to ignore it, all simulations are no
longer exactly done by classical statistical models. In this subsection, a
non-extensive statistical theory is nicely applied as a useful tool for
further simulating the entropy weight method above.

The Tsallis non-extensive statistics, as a widely-used generalization of
classical Boltzmann-Gibbs statistics, introduces a non-extensive quantity $q$
to modify it with a noise loss. When $q\to 1$ it returns to the classical
situation.

Specifically, in 1988 C. Tsallis firstly proposed the non-extensive statistics
TA . Although it is one of generalizations of the classical entropy, c.f.
Eq.(4), the Tsallis entropy

$S_{q}=\ln_{q}W\equiv\frac{W^{1-q}-1}{1-q}$ (18)

has been nicely applied into various kinds of fields. It is also often seen as
the form of probabilities,

$\displaystyle
S_{Ts}=E[\ln_{q}(\frac{1}{p_{i}})]=\frac{\sum_{i}^{m}p_{i}^{q}-1}{1-q}.$ (19)

Note that here the $q$-logarithm is introduced,

$\displaystyle\ln_{q}(x):=\frac{x^{1-q}-1}{1-q}\quad(x>0,q\neq 1),$ (20)

with its inverse function, the $q$-exponential function,

$\displaystyle\exp_{q}(x):=[1+(1-q)x]^{\frac{1}{1-q}}\quad(q\neq 1).$ (21)

Both of them will go back to the normal case when the limit of $q\to 1$ is
taken. Moreover, the most obvious feature of non-extensive entropy is that the
entropy of a system consisting of two subsystems is not a direct addition of
the entropic functions of the two subsystems, namely the extensive property of
the entropy function is no longer preserved. With this non-extensive
parameter, this theory could be more widely and successfully used in solving
problems where the traditional statistical models could not be applied. For
example, it has been proved to be successful in the fields of self-gravitation
systems, neural networks, turbulence, etc..

In this paper, we introduce the non-extensive theory in order to modify the
values of entropy weight factors by re-fitting the entropy weight curves
within kinds of similar algorithms. A non-idealized model will be then
constructed, which eventually covers the complex connections in systems. This
will also greatly improve the model’s noise immunity, increase the model’s
accuracy and reliability during data iteration, and make the model more
generalizable.

Similar to the classical case in Eq.(5), the generalized formula of
information entropy value can be deduced as

$\tilde{e}_{j}=S_{Ts}\left(p_{ij}\right)=\frac{\sum_{i}^{m}p_{ij}^{q}-1}{1-q}.$
(22)

By solving the constraint equation

$\tilde{e}_{j}=W_{r}(j)$ (23)

we could give out the selected $q$ values of different evaluation sequence in
the $q$-entropy weight method. Note that $W_{r}(j)$, shown in
Eq.(LABEL:equa-W), represents the grey correlation correction weight. At the
same time, for consistency it is considered that $q$ will be unity only when
the grey correlation sequence cannot give out an exact $q$. Finally, the
expected $q$-value is re-calculated as:

$q=\frac{\sum_{j=1}^{n}q_{j}}{n},$ (24)

and the normalized weights with the $q$ above are given below:

$W_{T}(j)=\frac{\tilde{e}_{j}}{\sum_{j=1}^{n}\tilde{e}_{j}}$ (25)

##  III Results and Discussions

###  III.1 Data Processing

Data were taken within five years when the number of suppliers is the most
stable for processing. And the four evaluation indicators of supply capacity
are listed hereafter.

  1. 1.

Supply stability

Supply stability is the variance value $D(X)$ of the total raw materials
provided by each supplier within five years,

$D\left(X_{i}\right)=E\left[X_{i}^{2}\right]-\mu_{i}^{2},$ (26)

where $E[X_{i}^{2}]$ is the mean square value of raw materials provided by a
single supplier in the past five years, and $\mu_{i}^{2}$ is the square of the
mean value of raw materials provided by a single supplier in the past five
years.

Data matrix X is obtained by the imported data sets and then transposed.
$D(X_{i})$ of each column of this generating matrix is calculated to get the
single-row matrix EV. It concludes nearly all the information on the supply
stability from every supplier, so that it will be defined as one indicator of
the supply stability correspondingly. Last not least, an empty matrix W is
also constructed and denoted as the original evaluation matrix with loading
the matrix EV.

  2. 2.

Supply quantity

Summing the matrix X row by row, we can obtain a single-row vector S about the
supply quantity of suppliers,

$S=\sum_{j=1}^{n}X_{ij}.$ (27)

As another evaluation indicator, S will be also loaded into the matrix W.

  3. 3.

Supply continuity

For indicating the supply continuity, we re-build one empty matrix comparing
with the matrix X. The number of 0 is for no supply of the current month while
1 is the opposite. Then do the re-summation column by column and obtain the
single-row matrix L. Being the indicator of evaluating the supply continuity,
it is then loaded into W.

  4. 4.

Ambiguous supply capacity

Furthermore, two different matrixes will be generated by MATLAB by subtracting
the data sets of both the enterprise’s order quantity and supplier’s supply
quantity. For the indicator matrix of ambiguous supply capacity, the
difference of these two should be made and re-summed column by column next.
The final single-row matrix T, as our last indicator, will be loaded into the
matrix W as well.

###  III.2 Data analysis

Relevant results are graphically shown in this part. Firstly the original EW-
TOPSIS model was considered and its normalized score table of $S_{i}$ (seen as
in Eq.(13)) for the $i$th ID of suppliers was given as shown in Table 1. For
simplicity only the first 48 of all 402 suppliers are listed, which will be as
our reference standards for comparisons between the original and non-extensive
EW-TOPSIS models.

Table 1: Normalized Scores in Eq.(13) of the first 48 of all 402 suppliers:

S/N | ID | Score | S/N | ID | Score | S/N | ID | Score | S/N | ID | Score
---|---|---|---|---|---|---|---|---|---|---|---
1 | 229 | 99.997 | 13 | 330 | 41.883 | 25 | 201 | 22.650 | 37 | 037 | 19.030
2 | 361 | 92.635 | 14 | 356 | 40.331 | 26 | 365 | 21.991 | 38 | 218 | 18.963
3 | 140 | 85.343 | 15 | 268 | 40.204 | 27 | 031 | 21.931 | 39 | 007 | 18.353
4 | 108 | 68.966 | 16 | 306 | 39.311 | 28 | 395 | 21.848 | 40 | 266 | 18.328
5 | 151 | 56.638 | 17 | 194 | 33.519 | 29 | 040 | 20.674 | 41 | 123 | 18.323
6 | 340 | 50.647 | 18 | 352 | 30.792 | 30 | 364 | 20.294 | 42 | 338 | 18.278
7 | 282 | 50.111 | 19 | 348 | 29.833 | 31 | 367 | 20.017 | 43 | 150 | 18.080
8 | 275 | 47.359 | 20 | 143 | 29.464 | 32 | 055 | 19.769 | 44 | 114 | 17.923
9 | 329 | 46.843 | 21 | 307 | 25.659 | 33 | 346 | 19.686 | 45 | 314 | 17.763
10 | 139 | 45.278 | 22 | 247 | 24.418 | 34 | 080 | 19.293 | 46 | 291 | 17.087
11 | 131 | 42.095 | 23 | 374 | 23.162 | 35 | 294 | 19.257 | 47 | 086 | 17.028
12 | 308 | 41.968 | 24 | 284 | 22.745 | 36 | 244 | 19.041 | 48 | 098 | 15.855

Next we go on looking into the characteristic sequence of system as well as
its contribution weight. Figure III-1 describes the weight distributions with
the four indicators above. Being the baseline, the result of classical EW
model, as in Eq.(7), is firstly present as the blue dotted line. All the other
weights are derived by Eq.(17) with different values of the gray
discrimination coefficient $\xi$.

In this iteration, we will start with a gray discrimination coefficient of
0.001, with a step size of 0.0001, up to 0.005. We will then select the gray
discrimination coefficient that is most similar to the original entropy weight
as the optimal coefficient for this iteration and incorporate it into future
revisions.

Figure III-1: The graph of grey correlation coefficient weight and original
weight of each index under different gray discrimination coefficient ($\xi$)

It is obvious to see that the GRA curve with $\xi=0.0013$, namely the red
dashed line, best fits the classical EW curve among all. Corresponding
correlation coefficients $\gamma(S_{i},x_{i}(j))$ are then calculated out,
based on Eq.(15) and Eq.(16). Finally the GRA weight factors with this $\xi$
is given out.

In order to clarify the difference between the classical and the $q$-EW-TOPSIS
models, the non-extensive normalized weights (seen in Eq.(25)) should be
studied too. GRA weight factors above for all four indicators are substituted
into the relation of Eq.(23) to get the non-extensive parameters $q_{i}$,
respectively, which are clearly collected in Table 2. Note that from now on
the average non-extensive parameter $q=\sum_{j}^{n}q_{j}/n$ is introduced
without losing of generality. In virtue of the validity of weight, the weight
factor in the $q$-EW-TOPSIS model is eventually solved, following Eq.(25).

Table 2: Non-extensive parameters of $q_{j}$ calculated by Eq.(23) for each
indicator:

Indicators | SS | SQ | SC | ASC | Average value
---|---|---|---|---|---
Parameter $q$ | 0.3216 | 0.1346 | 1.0001 | 0.0993 | 0.3888

Figure III-2: QEWM-TOPSIS score of each index Table 3: Non-Extensive Scores of
the first 48 of all 402 suppliers in $q$-EW-TOPSIS model:

S/N | ID | Score | S/N | ID | Score | S/N | ID | Score | S/N | ID | Score
---|---|---|---|---|---|---|---|---|---|---|---
1 | 229 | 99.994 | 13 | 330 | 41.553 | 25 | 284 | 22.016 | 37 | 244 | 18.204
2 | 361 | 92.615 | 14 | 356 | 39.980 | 26 | 395 | 21.800 | 38 | 218 | 18.125
3 | 140 | 85.316 | 15 | 268 | 39.852 | 27 | 365 | 21.240 | 39 | 338 | 17.582
4 | 108 | 68.855 | 16 | 306 | 38.946 | 28 | 031 | 21.177 | 40 | 007 | 17.509
5 | 151 | 56.451 | 17 | 194 | 33.059 | 29 | 040 | 19.880 | 41 | 266 | 17.484
6 | 340 | 50.411 | 18 | 352 | 30.277 | 30 | 364 | 19.489 | 42 | 123 | 17.480
7 | 282 | 49.870 | 19 | 348 | 29.470 | 31 | 367 | 19.204 | 43 | 150 | 17.239
8 | 275 | 47.091 | 20 | 143 | 28.919 | 32 | 55 | 18.949 | 44 | 114 | 17.107
9 | 329 | 46.570 | 21 | 307 | 25.313 | 33 | 346 | 18.863 | 45 | 314 | 16.934
10 | 139 | 45.028 | 22 | 247 | 23.741 | 34 | 037 | 18.607 | 46 | 291 | 16.294
11 | 131 | 41.767 | 23 | 201 | 22.697 | 35 | 080 | 18.461 | 47 | 086 | 16.283
12 | 308 | 41.639 | 24 | 374 | 22.446 | 36 | 294 | 18.424 | 48 | 098 | 15.100

By incorporating the non-extensive normalized weights into Eq.(22, 25), the
new model’s score sequence can be obtained. We will then order it and present
the top 48 in the table.

It is of great importance to research on the comparisons between these two
different models. As shown in Figure III-2, the scores of all suppliers within
these two different models are clearly compared. There are no big differences
among them all, especially the sequence ordered by scores. For instance, the
best 10 suppliers share the same ID numbers within both of models although
they own different scores. This is especially indicated when comparing the
data results in Table 1 and Table 3. For deeper clarification, Table 3 still
exhibits the top 48 suppliers with their scores.

In the above figure, the upper part shows the distribution of scores for two
evaluation sequences, and the lower figure shows the error analysis of the
original score and the $q$-EW TOPSIS score (measured by the ratio method). It
can be clearly seen that after reordering the supplier evaluation ID according
to the original EW-TOPSIS model score, the old model and the new model score
curve are compared. It can be seen that a certain decrease in the score of the
middle part of the sequence is made, and a certain compensation is made for
the score of the tail part of the sequence. Overall, the score fluctuates but
is not inconsistent with the facts, which proves that the $q$-EW-TOPSIS model
does not deviate too much from the original model’s score distribution and
verifies the reliability of the model. On the other hand, using the new model
for scoring, 99 suppliers’ rankings have changed, with an impact rate of
24.63$\%$. Therefore, in terms of supplier selection, this method is also
different from traditional evaluation methods. This ensures the correctness
and reliability of the model itself.

In order to better review this model, the difference proportion $\delta$ is
established as

$\delta(j)=\frac{\left|W_{S}(j)-W_{T}(j)\right|}{W_{T}(j)}.$ (28)

In this paper we have the values of it as 0.5352, 0.0096, 0.0606 , and 0.4741,
for four evaluation indicators, respectively. Normally the smaller value of it
leads to the greater influence of non-extensivity on the system. Its average
value is then given out, $\bar{\delta}=0.2396$.

In order to validate the model’s robustness and generalization capabilities,
we need to artificially manipulate the data to simulate the impact of data
changes on model weights under real conditions, in order to evaluate the
quality of our model.

This further demonstrates the robustness and generalizability of the model, as
it is able to maintain a stable weight distribution even when the amount of
data changes. Additionally, the small variance in the weight averages suggests
that the model is able to provide consistent results under different
conditions. Overall, this experiment shows that the model is able to perform
well in dynamic and unpredictable environments.

Table 4: The mean and variance of the weight of the original weight w0 and the
weight of this model w1 during the data iteration process.

Indicators | SS | SQ | SC | ASC
---|---|---|---|---
$w0_{ave}$ | 0.0057 | 0.7733 | 0.2135 | 0.0075
$w1_{ave}$ | 0.0064 | 0.7724 | 0.2130 | 0.0082
$w0_{var}$ | 0.0000 | 0.0031 | 0.0032 | 0.0075
$w1_{var}$ | 0.0000 | 0.0024 | 0.0025 | 0.0000

Figure III-3: Radar chart of the weights of the corresponding evaluation index
by taking the number of different suppliers as the step size.

Table 4 illustrates the average and variance of the original weight $w0$ and
the weight of the proposed model $w1$ throughout the entire data iteration
process. As shown, the mean of the weights for the proposed model is
relatively similar to that of the original model, however, the variance of the
weights for the proposed model is slightly better than the original model. To
better visualize the improvements of the proposed model, next the relative
error of each evaluation indicator is shown at each data update.

Next we set the average weight of each model after each iteration as the
standard weight, and plot the relative error of the weight (current period
weight - standard weight / standard weight) as a radar chart, seen in Figure
III-3. As can be seen from the chart, the weights determined by the new method
have errors that are generally smaller than those of the original model under
general conditions.

In summary, we have demonstrated the superiority of the proposed non-extensive
statistical theory in simulating the entropy weight method. By comparing the
relative errors of the weights calculated by the original and present models,
we have shown that the new-proposed model has better robustness and
generalization ability. The results of our study have verified the reliability
and effectiveness of the proposed model in mitigating the impact of noise and
incomplete information in data iteration, making it able to maintain high
accuracy and credibility under different data conditions.

It is constructive to plot the distributions of scores obtained in our $q$-EW-
TOPSIS model together with the others. As shown in Figure III-4, we compare
the scores among five different TOPSIS models, namely the classical and non-
extensive EW-TOPSIS models, IW-TOPSIS model, CRITIC-TOPSIS model and CV-TOPSIS
one.

In details, CV-TOPSIS (Coefficient of Variation TOPSIS), IW-TOPSIS
(independent weight coefficient method TOPSIS), and CRITIC-TOPSIS are variants
of the TOPSIS method. These variants introduce modifications to the original
TOPSIS methodology to address specific issues and improve its performance in
decision-making problems.

The CV-TOPSIS (Coefficient of Variation TOPSIS) method introduces the
coefficient of variation to consider the stability and relative difference of
each indicator data. In the CV-TOPSIS method, the coefficient of variation is
calculated for each decision criterion to evaluate the stability of the
criterion. The coefficient of variation represents the ratio of the standard
deviation of the data set of the criterion to its mean, and the smaller it is,
the more stable the criterion. Then, the coefficient of variation is used as a
weight and multiplied by each decision criterion’s score. By calculating the
distance between each decision alternative and the ideal solution, and
converting it into a composite score, the alternative is rankedCV1 ; CV2 ; CV3
.

The IW-TOPSIS (Independent Weight Coefficient Method TOPSIS) method is another
variation of the TOPSIS method that introduces an independent weight
coefficient to consider the importance and relative difference of each
criterion data. In the IW-TOPSIS method, the decision maker specifies the
independent weight of each decision criterion. The independent weight
represents the relative importance of the criterion and determines how the
weight is allocated to the decision alternatives. Then, by calculating the
distance of each decision alternative from the ideal solution and converting
it into a composite score, the alternatives are rankedIW1 ; IW2 ; IW3 .

The CRITIC-TOPSIS (Criteria Weights Determined by Interval Comparison of
Rankings) method is an improvement over the TOPSIS method. It introduces the
scoring method and the interval comparison method to determine the importance
and relative gap of each index data. In the CRITIC-TOPSIS method, the decision
maker first specifies the ideal solution for each decision criterion, which is
the best score the criterion can reach. The decision maker also specifies the
key factors of each decision criterion, which reflect the critical importance
of the criterion to the decision outcome. The key factors can be positive or
negative and they can respectively enhance or decrease the score of the
decision alternative. Then, by calculating the distance between each decision
alternative and the ideal solution, converting it into a composite score, the
alternatives are rankedCRITIC1 ; CRITIC2 .

In general, these three methods are extensions of TOPSIS that improve the
accuracy and efficiency of evaluation by expanding the evaluation criteria.
They are currently the three most widely used objective evaluation models.

Figure III-4: TOPSIS score distribution curves for the five objective
evaluation methods

Overall, it is clear from Figure III-4 that our model owns a higher
discrimination capability when evaluating suppliers with higher ranking, which
is in line with the current industry’s demand. Additionally, when compared to
CV-TOPSIS (green line), it also shows better discrimination in the later
stages of the score range. Furthermore, the score distribution of our model is
not significantly different from the standard entropy weight model, indicating
the accuracy of it.

In order to verify the performance of our model when data loss is negligible,
we need to examine whether the evaluation method of the gray correlation model
is superior to that of the $q$-EW-TOPSIS model. Additionally, we want to test
whether the model will revert to the weights that can be described by the
classical model when data is abundant. Therefore, we have chosen a dataset on
employee performance evaluation (for instance, the data in this work comes
from the open source data of Tempe dataset2 ), with 2,243 participants and 72
evaluation criteria. Our goal is to observe the model’s performance when data
is abundant and to verify its universality in both evaluation and decision-
making models.

Figure III-5: Error-bar plot of QEWM versus EWM when the sample size and the
number of evaluation indicators are large.

As depicted in Figure III-5, our proposed model, $q$-EW-TOPSIS, demonstrates
better differentiation among the evaluation indices when applied to a large
number of indices. This indicates that the model maintains a high level of
performance even when the data volume and number of evaluation indices are
substantial. The mean value of absolute error is found to be $5.02\%$, with a
range of $0.0795\%\sim 7.73\%$, all of which fall within an acceptable range.
Notably, the evaluation results exhibit a higher degree of fluctuation in the
middle portion of the dataset but remain relatively stable at the beginning
and end. This further confirms that $q$-EW-TOPSIS model can effectively
address the oversimplification issues encountered in classical models when
assessing and ranking the characteristics of suppliers. Additionally, this
modified method maintains a consistent order among all suppliers.

Furthermore, when the number of evaluation indices, as well as the size of
subjects, gets large enough, our model will approach to the normal one, which
means $q\to 1$.

##  IV Summary and Discussion

Summarizing, we have generalized the supply capacity evaluation model with
respect to the non-extensive statistical entropy. The traditional EW-TOPSIS
model is selected as our baseline and GRA method is used as the correlation
feedback to provide the amplitude and direction of model correction, and
finally brings out the non-extensive parameter $q$, which will determine the
$q$-EW-TOPSIS model in turn.

As discussed above, scores weighted by the classical EW-TOPSIS model
distribute too idealized and concentrated, which somehow fails describing the
characteristics of suppliers or loses some information accordingly. For the
GRA method, it succeeds in small sample demands or the system that covers few
regularities of data request. Thus it can help deduce the value of $q$, the
non-extensive parameter. With it the $q$-EW-TOPSIS model could be perfectly
established.

Compared to the results of the traditional EW-TOPSIS method, the $q$-EW-TOPSIS
model maintains its advantages by preserving the sorting order of suppliers,
while also expanding the scoring range for increased universality. This new
model is characterized by its ability to accurately evaluate the weights of
indicators, even with smaller sample sizes and weaker rules. Unlike
traditional models that require larger sample sizes for optimal performance,
the $q$-EW-TOPSIS model exhibits superior performance and minimal deviation as
a whole.

In this work, it is argued that the use of different indicator weights can
capture the ever-changing information in an incomplete system. The entropy
evaluation method is employed to perform numerical iterations in the system,
leading to a more stable and closer-to-complete structure. The mutual
restriction between indicators is used to overcome the limitations of the
traditional EW-TOPSIS model in the evaluation of supply capacity, resulting in
increased compatibility with the system and greater credibility of the
results. By using the entropy evaluation method, it can perform numerical
iterations in the system in order to achieve a more stable structure, which is
quite closer to the complete one. In this process, the mutual restriction
between different indicators is used to make up for the deficiency of the
traditional EW-TOPSIS model itself in the evaluation of supply capacity,
making it more compatible with the system and increasing the credibility of
the results.

Acknowledgments

This work has been partly supported by the funding for the Doctoral Research
of ECUT Nos. DHBK2019211. The author Ke-Ming Shen also wants to thank his Tuer
for kindly discussions.

## References

  * (1) Lima, J. F. R., Osiro, L., $\&$ Carpinetti, L. C. R., “Supplier selection and order allocation using a multi-objective simulated annealing algorithm,” Applied soft computing, vol. 21, pp. 194-209, 2014.
  * (2) George A. Zsidisin, et al., “Supplier selection: review and analysis of literature from 1994-2004,” International Journal of Physical Distribution and Logistics Management, 2004.
  * (3) Vijay R. Kannan, $\&$ Keah Choon Tan, “Supplier selection: An analytical hierarchy process (AHP) approach,” Journal of supply chain management, vol. 38, no. 3, pp. 11-21, 2002.
  * (4) Narges Banaeian, et al., “Supplier selection using a hybrid MCDM method based on grey relational analysis and DEMATEL,” Computers $\&$ Operations Research, vol. 89, pp. 337-347, 2018.
  * (5) Thomas L. Saaty, “The Analytic Hierarchy Process,” Springer, Berlin, Heidelberg, pp. 109-121, 1988.
  * (6) Thomas L. Saaty, “The Hierarchon: A Dictionary of Hierarchies,” International Journal of Services Sciences, vol. 1, no. 1, pp. 83-98, 2008.
  * (7) Mithat Zeydan, Cuneyt Colpan, and Cemal Cobanoglu, “A novel MCDM approach for supplier selection: Interval-valued 2-tuple linguistic grey relational analysis,” Expert systems with applications, vol. 38, no. 3, pp. 2741-2751, 2011.
  * (8) Abdullah Aldamak and Saeed Zolfaghari, “A new hybrid MCDM model for supplier selection using the combination of fuzzy TOPSIS and entropy,” Measurement, vol. 106, pp. 161-172, 2017.
  * (9) Celebi, Dilay, and Demet Bayraktar. “An integrated neural network and data envelopment analysis for supplier evaluation under incomplete information.” Expert Systems with Applications 35.4 : 1698-1710, 2008.
  * (10) Quan, Quan, and Zhongqiang Zhang. “Supply capability evaluation of intelligent manufacturing enterprises based on improved BP neural network.” Journal of Mathematics 2022: 1-8, 2022.
  * (11) Liu, Li, and Wenxue Ran. “Research on supply chain partner selection method based on BP neural network.” Neural Computing and Applications 32: 1543-1553, 2020.
  * (12) Velasquez, M. $\&$ Hester, P.T., “International journal of operations research.” International Journal of Operations Research, 10(2), 56-66, 2013.
  * (13) Behzadian, Majid, et al. “A state-of the-art survey of TOPSIS applications.” Expert Systems with applications 39.17: 13051-13069, 2012.
  * (14) Xing, Yahong, et al. “Evaluation System of Distribution Network Admission to Roof Distributed Photovoltaic Based on AHP-EW-TOPSIS.” 2021 IEEE 5th Conference on Energy Internet and Energy System Integration (EI2). IEEE, 2021.
  * (15) A. Renyi. Probability Theory. North Holland, Amsterdam, 1970.
  * (16) Tsallis, C., “Possible generalization of Boltzmann-Gibbs statistics,” J. Stat. Phys., vol. 52, pp. 479, 1988.
  * (17) Abe, Sumiyoshi. “A note on the q-deformation-theoretic aspect of the generalized entropies in nonextensive physics.” Physics Letters A 224.6: 326-330, 1997.
  * (18) G. Kaniadakis, “Non-linear kinetics underlying generalized statistics,” Physica A, vol. 296, pp. 405, 2001.
  * (19) Hwang Ching-Lai and Kwangsun Yoon, Multiple attribute decision making. Springer, Berlin, Heidelberg.
  * (20) Zhou, L., $\&$ Li, X. (n.d.). “Supplier selection based on grey system theory and improved entropy weight method.” The 8th China Intelligent Computing Conference and the 3rd Annual Conference of the China Branch of the International Electronic Commerce Federation.
  * (21) Deng, J.L., “The Introduction to Grey System Theory,” The Journal of Grey System, Vol. 1, No. 1, pp.1-24, 1989
  * (22) Luthra, Sunil, et al. “An integrated framework for sustainable supplier selection and evaluation in supply chains.” Journal of cleaner production 140: 1686-1698, 2017.
  * (23) Rajesh, R., and V. Ravi. “Supplier selection in resilient supply chains: a grey relational analysis approach.” Journal of Cleaner Production 86: 343-359, 2015.
  * (24) https://data.tempe.gov/datasets/city-of-tempe-2016-employee-survey-data/explore
  * (25) T. Sun, J. Wang and C. Wang, “TOPSIS threat assessment based on the coefficient of variation,” 2020 39th Chinese Control Conference (CCC), Shenyang, China, pp. 175-179, 2020.
  * (26) Vavrek, Roman, and Jana Chovancová. “Assessment of economic and environmental energy performance of EU countries using CV-TOPSIS technique.” Ecological Indicators 106: 105519, 2019.
  * (27) H. Su and D. Z. Du, “The Improved Coefficient of Variation TOPSIS Model and its Application,” Measurement, vol. 72, pp. 174-183, 2015.
  * (28) Hai-dong, W. U., et al. “Evaluation of Driving Behavior Economy Based on Big Data of New Energy Bus.” China Journal of Highway and Transport 35.3 : 177, 2022.
  * (29) Scientific Platform Serving for Statistics Professional 2021. SPSSPRO
  * (30) Jiang, Tingting, et al. “Strategy of Energy Conservation and Emission Reduction in Residential Building Sector: A Case Study of Jiangsu Province, China.” Journal of Environmental and Public Health 2023, 2023.
  * (31) Mohamadghasemi, A., HadiVencheh, A., $\&$ Hosseinzadeh Lotfi, F., The multiobjective stochastic CRITICTOPSIS approach for solving the shipboard crane selection problem. International Journal of Intelligent Systems, 35(10), 1570-1598,2020.
  * (32) Rostamzadeh, Reza, et al. “Evaluation of sustainable supply chain risk management using an integrated fuzzy TOPSIS-CRITIC approach.” Journal of Cleaner Production 175: 651-669, 2018.
  * (33) Chinese Society for Industrial and Applied Mathematics.,The ordering and transportation of raw materials for a manufacturing enterprise.question C, 2021. $http://www.mcm.edu.cn/html\\_cn/node/4d73a36cc88b35bd4883c276afe39d89.html$


# Visible metalenses with high focusing efficiency fabricated using
nanoimprint lithography

Andrew McClung  Department of Electrical and Computer Engineering, University
of Massachusetts Amherst, Amherst, MA 01003, USA  Mahsa Torfeh  Department of
Electrical and Computer Engineering, University of Massachusetts Amherst,
Amherst, MA 01003, USA  Vincent J. Einck  Department of Polymer Science and
Engineering,
University of Massachusetts Amherst, Amherst, MA 01003, USA  James J. Watkins
Department of Polymer Science and Engineering,
University of Massachusetts Amherst, Amherst, MA 01003, USA  Amir Arbabi
<EMAIL_ADDRESS>Department of Electrical and Computer Engineering, University
of Massachusetts Amherst, Amherst, MA 01003, USA

###### Abstract

Metasurfaces enable precise control over the properties of light and hold
promise for commercial applications. However, fabricating visible metasurfaces
suitable for high-volume production is challenging and requires scalable
processes. Nanoimprint lithography is a cost-effective and high-throughput
technique that can meet this scalability requirement. This work presents a
mask-templating nanoimprint lithography process for fabricating metasurfaces
with varying fill factors and negligible wavefront aberrations using composite
stamps. As a proof-of-concept, a 6 mm diameter metalens formed of silicon
nitride nano-posts with a numerical aperture of 0.2 that operates at 550 nm is
demonstrated. The nanoimprinted metalens achieves a peak focusing efficiency
of (81 $\pm$ 1)%, comparable to the control metalens made with electron beam
lithography with a focusing efficiency of (89 $\pm$ 1)%. Spatially resolved
deflection efficiency and wavefront data, which informs design and process
optimization, is also presented. These results highlight nanoimprint
lithography as a cost-effective, scalable method for visible metasurface
fabrication that has potential for widespread adoption in consumer electronics
and imaging systems.

## Introduction

Metasurfaces are optical elements created by arranging subwavelength
scatterers (known as ‘meta-atoms’) on a surface. Metasurfaces have several
desirable features: They can control the wavefront, polarization, and
amplitude of light with subwavelength resolution [1]. Compared to traditional
refractive optical elements, metasurfaces are significantly lighter and their
planar form factor simplifies integration into compact optical systems.
Furthermore, their design allows for cascading, enabling the realization of
optical metasystems [2, 3, 4, 5, 6, 7, 8]. These advantages have sparked
interest in the utilization of metasurfaces for commercial applications [9,
10]. However, widespread adoption will require metasurface designs and
fabrication processes suitable for high-volume production.

For metasurfaces to perform well, the dimensions and positions of their meta-
atoms must be precisely controlled. Subwavelength critical dimensions and
spacing are needed to achieve high efficiency, which presents challenges in
patterning visible and shorter wavelengths metasurfaces. Commonly, metasurface
patterns for academic research are defined using electron-beam lithography
(EBL), a technique that meets these precise positioning and resolution
requirements. Unfortunately, EBL has low throughput and is unsuitable for mass
production. While a fabrication process for metasurface lenses (metalenses)
using deep-ultraviolet photolithography, which is a scalable patterning
approach, has been demonstrated [11], the metalens pattern in that study was
constrained to minimum feature sizes equal to the wavelength used in
lithography (248 nm), limiting both numerical aperture (NA) and focusing
efficiency. Smaller features can be achieved by using shorter wavelength
illumination, multiple patterning, or other advanced photolithographic
techniques. However, the use of these technologies increases the cost of
manufacture substantially [12].

Nanoimprint lithography (NIL) is a promising alternative to photolithography
for applications where costs of nanoscale patterning dominate [12]. NIL is a
low-cost, high-throughput manufacturing method in which a stamp with nanoscale
features mechanically deforms an imprint material to create a desired pattern
[13]. While it is not yet competitive with photolithography for semiconductor
manufacture, current NIL technology is viable for devices with less stringent
demands in terms of defect density, overlay accuracy, and throughput. This
makes it an attractive manufacturing method for optical elements [14] like
polarizers [15, 16], microring resonators [17], waveguides [18] and photonic
crystals [19, 20].

NIL is also a promising candidate for scalable metasurface fabrication, as
evidenced by several recent demonstrations. NIL metasurface fabrication
processes fall into two main categories: direct imprint [21, 22, 23, 24] and
mask templating [25, 26, 27]. In direct imprint, the stamp shapes an ink that
is subsequently cured to form the metasurface. This approach requires very few
processing steps, which minimizes both fabrication time and cost. However,
only a few materials can be made into curable inks suitable for metasurfaces.

In mask templating, the stamp shapes a temporary layer that serves as an etch
or lift-off mask. Although it requires more processing steps than direct
imprint, mask templating offers the advantage of using material platforms
proven in EBL workflows, as all other processing steps (e.g., material
deposition or dry etching) remain the same.

Here, we present a mask-templating NIL process using tri-layer composite
stamps [28] for patterning metasurfaces with spatially varying fill factors
that operate in the visible wavelength range. The ability to pattern
structures with varying fill factors enables the realization of polarization-
insensitive metasurfaces operating based on propagation phase, in contrast to
geometric phase metasurfaces that are limited to circularly polarized light
[21, 23, 25, 26]. We evaluate the fabrication process by designing and
patterning a polarization-insensitive silicon nitride metalens with a 6 mm
diameter and an NA of 0.2 at a design wavelength of 550 nm. Silicon nitride is
an established material platform for visible wavelength metasurfaces [29, 30,
31]. The metalens produced by our process exhibits a nearly diffraction-
limited spot and a focusing efficiency of (81 $\pm$ 1)%. A control metalens
made using EBL exhibited a focusing efficiency of (89 $\pm$ 1)%, indicating
that our NIL metalens was fabricated with high fidelity. Furthermore, to
provide insights into the factors contributing to the slightly lower focusing
efficiency on NIL metalenses, we present spatially resolved deflection
efficiency by local probing and wavefront data obtained through a stable
holographic technique.

## Results and discussion

Figure 1:  Schematic of scalable NIL fabrication process. (a) A PDMS elastomer
stamp is cast from a master mold comprising holes in thin film amorphous
silicon on fused silica. The stamp is used to imprint a pattern in PMMA
resist. The pattern is transferred to a silicon nitride layer using a hard
mask and dry etch process. (b) A single master mold is used to make many
stamps, each of which is used to make many metasurfaces.

Our NIL patterning process, illustrated schematically in Figure 1a, involves
using a hybrid silicone elastomer stamp that is cast on a master mold to
imprint a resist on a metalens substrate. The thermal NIL process that we
employ is a well-studied process with proven repeatability and reliability
[32]. To make metasurfaces, we evaporate an etch mask into the resist, lift it
off, and use the etch mask to transfer the pattern into the substrate. Through
this approach, a single master mold can be used to make several stamps, each
of which can pattern several wafers by imprinting, enabling high-volume
replication of a pattern (Figure 1b) [33, 34].

To showcase the process and assess the performance of the fabricated
metasurfaces, we opted to design, fabricate, and characterize metalenses with
a diameter of 6 mm and NA of 0.2, resulting in a focal length of 14.7 mm, for
operation at a wavelength of 550 nm (Figure 2a). Metalenses that are corrected
for spherical aberrations have practical applications and serve as benchmarks
for evaluating the performance of different metasurface design and fabrication
techniques. The choice of a large aperture was deliberate, as the magnitude of
aberrations scale with the diameter of the aperture. We anticipated that the
NIL process might introduce more substantial aberrations than EBL due to
potential errors in the placement of the zone boundaries caused by stamp
stretching and deformation during the stamp fabrication and the imprinting
process. Therefore, a large aperture metalens, which is more sensitive to such
aberrations, is suitable for evaluating the performance of the NIL process.

As shown in Figure 2a, the metalens is composed of silicon nitride nano-posts
that are 650 nm in height. These nano-posts are positioned on a fused silica
substrate and have hexagonal cross-sections. They are arranged on a triangular
lattice with a lattice constant of 430 nm, and their widths vary between 100
nm and 310 nm. The minimum nano-post width and gap of 100 nm and 120 nm were
not imposed by the NIL process but were selected to facilitate potential
wafer-scale fabrication of the master mold using a cost-effective process such
as immersion deep ultraviolet (DUV) lithography. Figure 2b shows the simulated
transmittance and phase characteristics for periodic nano-posts arrays as a
function of their width. To determine the design parameters and estimate the
metalens focusing efficiency, we employed the grating averaging technique [1].
This technique allowed us to establish a design, depicted in Figure 2c, that
correlates a desired phase with the width of the nano-post. We obtained a
simulated focusing efficiency of 90.2% for the metalens shown in Figure 2,
which has an NA of 0.2. We also designed and simulated silicon nitride
metalenses with larger NAs while limiting the minimum nano-post widths to 100
nm and the minimum gaps between them to 120 nm. The results are presented in
Figure S1 and indicate that a focusing efficiency of more than 75% can be
achieved for metalenses having an NA smaller than 0.7.

Figure 2: Metalens design. (a) A metalens with a diameter of 6 mm and a focal
length of 14.7mm (NA of 0.2) is designed. The metalens is formed of hexagonal
silicon nitride nano-posts arranged on a triangular lattice and resting on a
fused silica substrate. (b) Phase and transmittance as a function of nano-post
width. (c) The design curve that maps the desired phase to the nano-post
width.

To create the metalens pattern, we fabricated a master mold by EBL and dry
etching holes in a 200-nm-thick film of amorphous silicon deposited on a fused
silica substrate (see Supporting Information for details). In the fluorine
plasma chemistry we use, the etch rate of amorphous silicon is substantially
higher than that of fused silica, and the use of amorphous silicon on fused
silica substrate ensures that all the holes are etched to nearly the same
depth. After fabrication of the master mold, we fluorinated it to reduce
adhesion during stamp fabrication.

To make a stamp, we spun a hard PDMS (hPDMS) film onto the master mold and
cured it by baking it in an oven. We then placed the hPDMS side of the mold in
an uncured puddle of commercial PDMS that was resting on a flexible glass
backing, as illustrated in Figure 3a. Cured PDMS spacers with a thickness of
140 $\mu$m were employed to maintain a fixed distance between the master mold
and the glass backing. The use of bilayer PDMS with optimized thicknesses and
the flexible yet inelastic glass backing ensures that the stamp does not
stretch or deform during the imprint process because stamp deformation can
lead to substantial metasurface wavefront aberrations. The entire assembly was
then baked to cure the PDMS, and the mold was released to yield a stamp
comprising nano-posts that are approximately 200 nm in height. Because the
stamp is used to template a mask, it does not need to have the same feature
height as the final device: all features in the stamp have an aspect ratio of
2 or smaller. The small aspect ratio of the stamp features enables patterning
with a resolution down to 10 nm [36].

Figure 3: Fabrication process of NIL metalenses. (a) Stamp assembly: hPDMS is
spin-coated on the master mold, partially cured, and then adhered to a
flexible glass substrate by a commercial PDMS (Sylgard 184). Spacers set the
thickness of the stamp. (b) Thermal NIL: The PDMS stamp is used to imprint
features in PMMA. (c) Hard mask deposition: an Al2O3 mask is evaporated into
the resist template and lifted off. (d) Etching: the hard mask is used to etch
metalens patterns into silicon nitride. (e) Photograph and SEM image of the
master mold, (f) a photograph of the stamp, and (g) a photograph and an SEM
image of an imprinted metalens. The roughness observed in the SEM image is due
to the gold coating that was done to mitigate charging during imaging.

To fabricate the metalenses, we first spun a layer of polymethyl methacrylate
(PMMA) on top of a silicon nitride film that was deposited on a fused silica
substrate. We then used a commercial thermal NIL tool to imprint the metalens
patterns into the PMMA layer. The imprinting process is depicted in Figure 3b.
To remove the thin residual layer of PMMA at the bottoms of the features, we
etched the imprint in oxygen plasma. The optical and scanning electron
microscope (SEM) images of the imprinted PMMA displayed in Figure S2 in the
Supporting Information demonstrate that the patterns were transferred with
high fidelity and minimal defects. Following the successful imprinting of the
patterns, we evaporated an aluminum oxide etch mask into the PMMA (Figure 3c)
and lifted it off. Finally, we etched the metalens patterns into the silicon
nitride layer, resulting in the fabrication of metalenses (Figure 3d).
Photographs and SEM images of the master mold, a stamp and an imprinted
metalens are shown in Figure 3e–g. For performance comparison, we also
produced a control metalens with the same design as the NIL metalens but using
EBL. In addition, the experimental results for a second NIL metalens that was
fabricated using the same NIL process but with a sub-optimal silicon nitride
etch process are presented in Figure S3. A more comprehensive description of
the fabrication processes for the master mold, stamps, metalenses, and the EBL
metalens can be found in the Supporting Information.

We characterized the performance of the NIL metalens by measuring its focal
spot and focusing efficiency. The focal spot shape mostly provides information
about wavefront aberrations caused by potential errors in the zone boundary
locations (i.e., positions of the nano-posts) that could arise from the
deformation of the stamp during the imprinting process. Aberrations of this
nature are typically not observed in metalenses patterned using EBL or
photolithography due to their high positioning accuracy [37, 38]. On the other
hand, the focusing efficiency, which represents the proportion of incident
power that is focused by the metalens, is used to evaluate the accuracy of the
nanostructure pattern generated through the NIL process. The focusing
efficiency is sensitive to the shapes and dimensions of the meta-atoms,
allowing for an assessment of the fidelity of the nanostructure pattern
created using NIL.

We measured the focal spots of the NIL and EBL metalenses using a custom-built
microscope, as shown schematically in Figure 4a. We passed a laser beam with
wavelength of 552 nm through a beam expander and illuminated the metalenses,
which had a diameter of 6 mm and an NA of 0.2 at 550 nm. An objective lens and
a tube lens magnified and imaged the focal spots of metalenses onto a camera
(see Supporting Information for details). Figure 4b presents the measured
focal spot images of the NIL and EBL metalenses, and the simulated focal spot
of an ideal metalens with the same NA. An ideal metalens is a local metalens
that directs all the normally incident power to its focal spot without any
spherical aberrations [39]. The focal spots of ideal metalenses with an NA
smaller than 0.8 can be well approximated with an Airy disk [39]. We used the
focal spot data in Figure 4b to calculate the modulation transfer function
(MTF) of the metalenses. The MTFs of the NIL and EBL metalenses, as well as
the simulated MTF of a diffraction-limited ideal metalens, are shown in Figure
4c. The similarity between the aberrations observed in the focal spots of the
NIL and EBL metalenses suggests that the slight deviations from the
diffraction limit may be due to the nonuniform intensity distribution and
wavefront aberrations of the incident expanded beam. Nonetheless, the nearly
diffraction-limited response of the NIL metalens confirms the negligible error
in the position of the nano-posts. It is worth noting that when the process
parameters during the imprinting process were suboptimal, we observed
significant aberrations in the focal spots of the NIL metalenses. Hence,
optimizing the NIL process is essential for achieving high-quality
metasurfaces.

Figure 4: Metalens characterization. (a) Schematic of image relay used for
measuring focal spots. (b) Focal spot images of NIL and EBL metalenses and
that of an ideal metalens. (c) Modulation transfer functions of NIL, EBL, and
ideal metalenses. (d) Schematic of measurement setups used for determining
metalens focusing efficiency. (e) Measured focusing efficiency values for NIL
and EBL metalenses at a few wavelengths between 532 and 633 nm.

We also conducted measurements to determine the focusing efficiency of the
metalenses. The focusing efficiency is defined as the fraction of power
incident on the metalens that is focused by the metalens. In Figure 4d, we
schematically depict the setup used for these measurements. To quantify the
optical power focused by the metalens, we employed a 6-mm-diameter aperture
positioned behind the metalens, along with a 40-$\mu$m-diameter pinhole placed
at the focal plane of the metalens, and measured the optical power that passed
through the pinhole. The purpose of the aperture was to reduce the size of the
incident beam, matching it to the aperture size of the metalens. On the other
hand, the pinhole served to block the light that passed through the metalens
but did not converge to the desired focal spot, including the zeroth and other
transmission diffraction orders apart from the first order. To ensure that the
optical power passing through the pinhole predominantly originated from the
focused light and had minimal contributions from unfocused light, we conducted
additional tests. Specifically, we obtained a high dynamic range image of the
intensity distribution at the focal plane of the NIL metalens and compared it
to the simulated intensity distribution of an ideal metalens. Detailed
information about these procedures can be found in the Supporting Information.
Figure S4 in the Supporting Information presents the resulting intensity
distributions, along with the encircled power percentage as a function of the
pinhole diameter. The simulated encircled power percentage indicates that
98.6% of the power focused by an ideal metalens passes through the
40-$\mu$m-diameter pinhole.

We measured the incident power by removing the metalens and the pinhole and
measuring the optical power that passed through the 6-mm-diameter aperture. To
determine the focusing efficiency, we divided the power transmitted through
the pinhole by the incident power and accounted for corrections due to the
Fresnel reflection at the fused silica–air interface at the backside of the
metalens substrate (3.5%) and the percentage of power focused by an ideal
metalens that is blocked by the pinhole (1.4%). Figure 4e shows the measured
values of the focusing efficiency for NIL and EBL metalenses at wavelengths
between 532 nm and 633 nm. The NIL metalens exhibits a peak focusing
efficiency of (81 $\pm$ 1)% at the wavelength of 552 nm which is close to the
(89 $\pm$ 1)% obtained for the EBL metalens at the same wavelength. The
focusing efficiency data are also presented in Table S1 in the Supporting
Information. The high focusing efficiency observed for the NIL metalens
validates the high accuracy of the nanostructure patterns produced using the
NIL technique. This outcome is consistent with the SEM inspection results
depicted in Figure 3g, affirming the suitability of NIL for fabricating high-
quality visible wavelength metasurfaces.

The focusing efficiency, which is a coherent average of efficiency over the
metalens aperture, does not capture the potential efficiency variation across
the aperture or provide insights into its underlying causes. Identifying and
understanding potential causes of efficiency variation are crucial for
maximizing the metalens efficiency by adjusting the metalens design or
fabrication parameters such as the NIL process or etching parameters. To
detect such variations and gain insights into their potential causes, we
conducted two additional measurements: local deflection efficiency and high-
resolution holographic complex field measurements. We conducted deflection
efficiency measurements at various locations on the NIL metalens by
illuminating small regions of the aperture with a loosely focused beam (ca. 50
$\mu$m diameter). The purpose of these measurements was to determine the
percentage of power directed toward the metalens’ focal point. Figure S5 in
the Supporting Information presents the results, which demonstrate that the
deflection efficiency rises as the distance from the center of the metalens
increases.

Figure 5: Holographic characterization. (a) Schematic of the optical setup
used for holographic complex field measurement. (b) Measured phase profile
over a small region on the surface of metalens. The blue dot indicates the
location of (x, y) = (1000 $\mu$m, 0 $\mu$m) assuming the center of the
metalens is at (x, y) = (0 $\mu$m, 0 $\mu$m). (c) The Difference between the
measured and ideal phase over the region shown in (b). (d) Left: Diagram of
metalens showing regions of interest. Regions i, ii and iii lie at distances
of 500 $\mu$m, 1000 $\mu$m and 2000 $\mu$m from the metalens center. Right:
Line cuts of the measured and ideal phase difference in regions i, ii, and
iii.

The spatially resolved deflection efficiency measurement provides valuable
information and aligns with the focusing efficiency values obtained using the
pinhole. However, it does not reveal the underlying causes of the lower
efficiency observed near the center of the metalens. To further investigate
this issue, we employed a high-resolution digital holography technique [40,
41, 2] to measure the complex-valued optical field on the output aperture of
the metalens. Figure 5a illustrates the schematic of the holographic
microscope used in the experiment. A coherent illumination at 552 nm was
projected onto a diffraction grating, and the resulting image of the metalens
surface was captured after passing through a relay. By filtering the first
diffraction order of the grating in the Fourier plane of the relay, a coherent
reference beam was created, which interfered with the relayed image. We
captured images of the metalens surface and utilized a Fourier technique to
obtain a high-resolution map of the complex-valued optical field on the
metalens surface. The holographic imaging setup and measurements are described
in greater detail in the Supporting Information.

The extracted phase of the field at a region approximately 1 mm from the lens
center is shown in Figure 5b, and the phase error, which is the difference
between the measured phase and the desired phase profile, are shown in Figure
5c. Line cuts showing the phase error at different distances from the metalens
center are presented in Figure 5d. Notably, larger phase errors are noticeable
near the center of the metalens, but they decrease in magnitude as the radius
increases, in line with the trend observed in the local deflection efficiency
results. Moreover, the phase error indicates that the realized phase coverage
over a Fresnel zone exceeds $2\pi$. This could be attributed to factors such
such as nano-post widths being larger than the design values. To further
quantify the effect of such errors, we used the grating averaging technique
[1], and found the focusing efficiency of the metalens as a function of
systematic errors in the nano-post height and width, and random errors in the
nano-post widths. The results of this study are presented in Figure S6 and
indicate that the focusing efficiency is more sensitive to errors in the nano-
post widths than their height.

The NIL metalens was primarily designed for optimal performance at a
wavelength of 550 nm. However, the focusing efficiency values depicted in
Figure 4e indicate that it can also produce high-contrast images at other
wavelengths within the visible spectrum. To demonstrate its imaging
properties, we assembled a camera using the metalens and an RGB image sensor,
shown schematically in Figure 6a: we illuminated the object with different
light sources—two narrowband LED sources, whose spectra are shown in Figure
6b, and the fluorescent lamps installed in our lab—placed the metalens camera
assembly approximately 30 cm from the object, and adjusted the distance
between metalens and image sensor to obtain the best focus. Figure 6c shows
images captured with our camera. Clear images are formed when the object is
illuminated with narrowband green and red light sources, and the image formed
under broadband illumination exhibits the chromatic aberration expected for
metalenses that employ phase wrapping. To provide a color reference and for
comparison, we imaged the object using the same image sensor and an aspheric
lens with similar NA (Figure 6c, right). The varying focusing efficiency of
the NIL metalens across the visible spectrum can alter the apparent color of
objects within the scene. However, this effect can be readily corrected by
adjusting the white balance of the image, which is a digital post-processing
step that has not been performed on photos shown in Figure 6c.

Figure 6: Imaging demonstration. (a) Schematic of the imaging setup. Light
emanating from a distant object is focused by the metalens onto a color CMOS
image sensor. (b) Spectra of narrowband green and red LED illumination sources
used in imaging. (c) Comparison of photos taken with a 6-mm-diameter metalens
under narrowband green, narrowband red, and white light illumination, and a
refractive aspheric lens. The photos represent the raw data captured by the
image sensor, and no processing has been performed on them.

## Conclusion and outlook

The scalable fabrication of high-efficiency metasurfaces operating in the
visible wavelength range opens up opportunities for their widespread adoption
in consumer applications. Unlike some of the previous NIL techniques developed
for visible metasurfaces, the approach described in this study offers distinct
advantages by not imposing limitations on the metasurface materials or their
patterns, including fill factors and shapes. Our methodology employs an
amorphous silicon on fused silica master mold, which ensures consistent height
of the stamp features, along with a thin hard mask that allows for the
utilization of a thin PMMA layer. This combination enables precise patterning
of nanostructures with spatially varying fill factors, achieving high fidelity
in the process. The technique’s high resolution facilitated the patterning of
small features for visible metasurfaces, and it could potentially be extended
to pattern metasurfaces even at UV wavelengths.

At scale, this fabrication method could significantly reduce the cost of
producing visible-wavelength metalenses compared to equivalent metalenses
created using photolithography. Furthermore, the nano-post meta-atoms utilized
in this study enable the implementation of metasurfaces that not only
manipulate the phase of light but also control its polarization. As a result,
the NIL technique holds promise for the cost-effective fabrication of such
elements, expanding the range of applications and functionalities achievable
in metasurface devices. In fact, while preparing this work for publication
[43], we learned that a mask templating process for manufacturing niobium
oxide metalenses, though with a lower efficiency and a smaller diameter, has
been recently developed by a commercial company [44], further validating the
advantages of the NIL technique for the manufacturing of metalenses.

In numerous applications, the efficiency of metasurfaces is of utmost
importance, making its accurate measurement crucial. The approach discussed
here, which involves direct focusing efficiency measurement by measuring the
optical power passing through a small-diameter pinhole, combined with high
dynamic range imaging for precise capture of the transmitted light, can serve
as the standard technique for comparing and benchmarking different design and
fabrication methods for metalenses. It is worth noting that when measuring
high numerical aperture (NA) metalenses, a power meter equipped with an
integrating sphere may be required to ensure accurate measurement of the
focused power.

While the results of focusing efficiency measurements can be validated through
local deflection efficiency measurements, it is important to note that
deflection efficiency measurement cannot replace the focusing efficiency
measurement due to the loss of phase information in the former approach [44].
However, local deflection efficiency measurements can provide additional
insights and complement the focusing efficiency data. Furthermore, the
holographic complex-field measurement technique introduced in this study has
valuable applications in troubleshooting and optimizing the design and
fabrication process. It offers direct information about the local optical
field, which cannot be obtained through intensity-based characterization
techniques alone. Therefore, this holographic complex-field measurement
technique can be employed for thorough analysis, fine-tuning, and optimization
of metasurface design and fabrication methods.

The high focusing efficiency demonstrated here indicates the feasibility of
cascaded multilayer meta-optics using metasurfaces fabricated through the NIL
technique. In such systems, the total efficiency is the product of the
efficiency of each individual metasurface element. Notably, the fabrication of
separate metasurface wafers through NIL allows for wafer-level alignment. This
alignment can be achieved using 3D metasurface alignment marks [45], which
offer sub-micron accuracy. As a result, fully-integrated optical metasystems
can be realized, eliminating the need for post-fabrication alignment and
assembly processes. This advancement enables the development of highly
efficient and compact optical systems that can be seamlessly integrated into
various applications.

## References

  * [1] S. M. Kamali, E. Arbabi, A. Arbabi, A. Faraon,  _Nanophotonics_ 2018, _7_ , 6 1041.
  * [2] A. Arbabi, E. Arbabi, S. M. Kamali, Y. Horie, S. Han, A. Faraon,  _Nat. Commun._ 2016, _7_ , 1 13682.
  * [3] B. Groever, W. T. Chen, F. Capasso,  _Nano Lett._ 2017, _17_ , 8 4902.
  * [4] M. Faraji-Dana, E. Arbabi, A. Arbabi, S. M. Kamali, H. Kwon, A. Faraon,  _Nat. Commun._ 2018, _9_ , 1 4196.
  * [5] H. Kwon, E. Arbabi, S. M. Kamali, M. Faraji-Dana, A. Faraon,  _Nat. Photonics_ 2020, _14_ , 2 109.
  * [6] A. McClung, S. Samudrala, M. Torfeh, M. Mansouree, A. Arbabi,  _Sci. Adv._ 2020, _6_ , 38 eabc7646.
  * [7] A. Arbabi, A. Faraon,  _Nat. Photonics_ 2023, _17_ , 1 16.
  * [8] D. N. Neshev, A. E. Miroshnichenko,  _Nat. Photonics_ 2023, _17_ , 1 26.
  * [9] STMicroelectronics and Metalenz Partner to Transform Optical Sensing for Consumer, Automotive, and Industrial Applications (Press Release), 2021.
  * [10] NIL Technology launches flat NIR metalenses ready for mass production (Press Release), 2021.
  * [11] J.-S. Park, S. Zhang, A. She, W. T. Chen, P. Lin, K. M. Yousef, J.-X. Cheng, F. Capasso,  _Nano Lett._ 2019, _19_ , 12 8673.
  * [12] International Roadmap for Devices and Systems,  Technical report, Institute of Electrical and Electronics Engineers, 2022.
  * [13] H. Schift,  _J. Vac. Sci. Technol., B: Microelectron. Nanometer Struct.–Process., Meas., Phenom._ 2008, _26_ , 2 458.
  * [14] S. Barcelo, Z. Li,  _Nano Convergence_ 2016, _3_ , 1 1.
  * [15] Z. Yu, P. Deshpande, W. Wu, J. Wang, S. Y. Chou,  _Appl. Phys. Lett._ 2000, _77_ , 7 927.
  * [16] S.-W. Ahn, K.-D. Lee, J.-S. Kim, S. H. Kim, J.-D. Park, S.-H. Lee, P.-W. Yoon,  _Nanotechnology_ 2005, _16_ , 9 1874.
  * [17] C.-y. Chao, L. J. Guo,  _J. Vac. Sci. Technol., B: Microelectron. Nanometer Struct.–Process., Meas., Phenom._ 2002, _20_ , 6 2862.
  * [18] N. Koo, J. Wuk Kim, M. Otto, C. Moormann, H. Kurz,  _J. Vac. Sci. Technol., B: Nanotechnol. Microelectron.: Mater., Process., Meas., Phenom._ 2011, _29_ , 6 06FC12.
  * [19] R. Ji, M. Hornung, M. A. Verschuuren, R. van de Laar, J. van Eekelen, U. Plachetka, M. Moeller, C. Moormann,  _Microelectron. Eng._ 2010, _87_ , 5-8 963.
  * [20] M. R. Beaulieu, N. R. Hendricks, J. J. Watkins,  _ACS Photonics_ 2014, _1_ , 9 799.
  * [21] G. Yoon, K. Kim, D. Huh, H. Lee, J. Rho,  _Nat. Commun._ 2020, _11_ , 1 2268.
  * [22] V. J. Einck, M. Torfeh, A. McClung, D. E. Jung, M. Mansouree, A. Arbabi, J. J. Watkins,  _ACS Photonics_ 2021, _8_ , 8 2400.
  * [23] M. Miyata, A. Kohno, T. Komatsu, T. Hashimoto,  _Opt. Mater. Express_ 2022, _12_ , 10 4169.
  * [24] D. E. Jung, I. R. Howell, V. J. Einck, F. D. Arisoy, L. D. Verrastro, A. McClung, A. Arbabi, J. J. Watkins,  _ACS Appl. Nano Mater._ 2023, _6_ , 3 2009.
  * [25] G.-Y. Lee, J.-Y. Hong, S. Hwang, S. Moon, H. Kang, S. Jeon, H. Kim, J.-H. Jeong, B. Lee,  _Nat. Commun._ 2018, _9_ , 1 4562.
  * [26] C. A. Dirdal, G. U. Jensen, H. Angelskår, P. C. V. Thrane, J. Gjessing, D. A. Ordnung,  _Opt. Express_ 2020, _28_ , 10 15542.
  * [27] U. Quaade, V. E. Johansen, J. M. Llinas, J. Fly, A. Samadi, M. Schmidlin, N. A. Mortensen,  In _Photonic and Phononic Properties of Engineered Nanostructures XII_. SPIE, 2022 PC1201012.
  * [28] T. W. Odom, J. C. Love, D. B. Wolfe, K. E. Paul, G. M. Whitesides,  _Langmuir_ 2002, _18_ , 13 5314.
  * [29] A. Zhan, S. Colburn, R. Trivedi, T. K. Fryett, C. M. Dodson, A. Majumdar,  _ACS Photonics_ 2016, _3_ , 2 209.
  * [30] S. Colburn, A. Zhan, E. Bayati, J. Whitehead, A. Ryou, L. Huang, A. Majumdar,  _Opt. Mater. Express_ 2018, _8_ , 8 2330.
  * [31] Z.-B. Fan, Z.-K. Shao, M.-Y. Xie, X.-N. Pang, W.-S. Ruan, F.-L. Zhao, Y.-J. Chen, S.-Y. Yu, J.-W. Dong,  _Phys. Rev. Appl._ 2018, _10_ , 1 014005.
  * [32] S. Y. Chou, P. R. Krauss, P. J. Renstrom,  _J. Vac. Sci. Technol., B: Nanotechnol. Microelectron.: Mater., Process., Meas., Phenom._ 1996, _14_ , 6 4129.
  * [33] H. Schmitt, P. Duempelmann, R. Fader, M. Rommel, A. Bauer, L. Frey, M. Brehm, A. Kraft,  _Microelectron. Eng._ 2012, _98_ 275\.
  * [34] N. Tucher, O. Höhn, H. Hauser, C. Müller, B. Bläsi,  _Microelectronic Engineering_ 2017, _180_ 40\.
  * [35] A. Arbabi, E. Arbabi, M. Mansouree, S. Han, S. M. Kamali, Y. Horie, A. Faraon,  _Sci. Rep._ 2020, _10_ , 1 1.
  * [36] A. Pandey, S. Tzadka, D. Yehuda, M. Schvartzman,  _Soft Matter_ 2019, _15_ , 13 2897.
  * [37] P. Fischer, S. Y. Chou,  _Appl. Phys. Lett._ 1993, _62_ , 23 2989.
  * [38] G. Wallraff, W. D. Hinsberg,  _Chem. Rev._ 1999, _99_ 1801\.
  * [39] A. McClung, M. Mansouree, S. Samudrala, A. Arbabi,  In _Conference on Lasers and Electro-Optics_ , OSA Technical Digest. Optical Society of America, Washington, DC, 2020 FM1R.1.
  * [40] C. Edwards, A. Arbabi, G. Popescu, L. L. Goddard,  _Light: Sci. Appl._ 2012, _1_ , 9 e30.
  * [41] G. Popescu, T. Ikeda, R. R. Dasari, M. S. Feld,  _Opt. Lett._ 2006, _31_ , 6 775.
  * [42] B. Bhaduri, C. Edwards, H. Pham, R. Zhou, T. H. Nguyen, L. L. Goddard, G. Popescu,  _Adv. Opt. Photonics_ 2014, _6_ , 1 57.
  * [43] A. McClung, M. Torfeh, V. J. Einck, J. J. Watkins, A. Arbabi,  In _High Contrast Metastructures XI_. SPIE, 2022 PC120110V.
  * [44] B. R. Williams, D. Bacon-Brown, R. Edwards, S. Nielson, A. Farhang, S. Johnson, R. Petrova, M. C. George, S. Ogden,  In _High Contrast Metastructures XII_ , volume 12432. SPIE, 2023 37–43.
  * [45] A. McClung, B. Mirzapourbeinekalaye, A. Arbabi,  In _High Contrast Metastructures XII_. SPIE, 2023 PC124320W.

## Acknowledgements

This work was supported by the DARPA Extreme Optics and Imaging program.
Fabrication was performed at the core facilities of the Institute for Applied
Life Sciences (IALS) at the University of Massachusetts Amherst and at the
Center for Nanoscale Systems (CNS) at Harvard University, which is part of the
National Nanotechnology Coordinated Infrastructure Network (NNCI). The NNCI is
supported by the National Science Foundation under NSF Award No. 1541959.

## Supporting Information

Master design. We determined the design curve (Figure 2c) using a method
described in a previous publication [1]. In the design, we assumed refractive
indices of 2.02 for silicon nitride and 1.46 for fused silica. We explored
different designs with varying nano-post heights and lattice constants. Using
the grating averaging technique, we evaluated the focusing efficiency of
metalenses with a numerical aperture (NA) of 0.2 for each design. Among the
designs considered, we selected the one shown in Figure 2a, which exhibited a
simulated focusing efficiency of 95%.

Master mold fabrication. To make the master mold, we first deposited a 200-nm-
thick amorphous silicon film on a 0.5-mm-thick fused silica substrate using a
PECVD (PECVD, STS). We then spun on a positive tone electron beam resist (ZEP
520A, Zeon Specialty Materials) and a charge dissipation layer (ARPC 5090,
Allresist) and exposed the metalens pattern using an electron beam lithography
system (ELS-HS50, Elionix). We removed the charge dissipation layer, developed
the resist, and transferred the pattern, which consisted of holes, to the
amorphous silicon film by dry etching (ICP RIE, STS). We removed the resist
using a solvent (Remover PG, Kayaku Advanced Materials) and treated the
surface using an oxygen plasma asher (SCE-104, Anatech USA). We fluorinated
the master mold by exposing it to
heptadecafluoro-1,1,2,2-tetrahydrodecyltrichlorosilane (Gelest) in a
desiccator pumped to rough vacuum and left overnight. Finally, we glued the
master mold to a glass microscope slide using a UV-curable adhesive (NOA 84,
Norland Products Inc.). Figure 3e shows an image of the master mold at this
stage.

Stamp fabrication. To cast a stamp, we first spun a layer of PDMS (Sylgard
184, Dow) on a flexible, 200-$\upmu$m-thick glass substrate (AF-32, Schott) to
a thickness of ca. 140 $\upmu$m, and baked the sample overnight at 80 ºC.
After the PDMS was cured, we cut a channel approximately the width of the
metalens pattern into the PDMS layer using a razor blade. We then prepared an
hPDMS solution consisting of 1.7 g of 7-8% vinyl methysiloxane-
dimethylsiloxane (Gelest), 5 $\upmu$L of 2,4,6,8-tetramethyl-2,4,6,8
tetravinylcyclotetrasiloxane (Sigma Aldrich), 9 $\upmu$L of platinum-
divinyltetramethyldisiloxane (Gelest) and 5 g of toluene. We vortex mixed this
solution, added 0.5 mL of 25-30% methylhydrosiloxane-dimethylsiloxane
(Gelest), and mixed again. We spun this solution onto the master mold at 3
kRPM for 40 seconds and baked the sample in an oven at 100 ∘C for 45 minutes.
We then poured a puddle of uncured PDMS (Sylgard 184, Dow), into the channel
cut into the PDMS on glass workpiece, taking care to avoid forming bubbles,
and placed the master mold, hPDMS side down, in this puddle, pressing to
ensure the edges of the master mold were flush with the spun and cured PDMS.
We baked this assembly overnight in an oven at 80 ∘C. To release the stamp, we
let the assembly cool, submerged it in a dish of isopropyl alcohol, and gently
flexed the glass backing away from the master mold, holding the master mold
against the bottom of the dish using the slide. We removed the stamp from the
isopropyl alcohol, rinsed it with deionized water, and blew it dry with
nitrogen. Figure 3f shows an image of a stamp at this stage.

Metalens fabrication. To make the NIL metalenses, we first deposited a ca.
650-nm-thick silicon nitride film on a 0.5-mm-thick fused silica substrate
using a PECVD (PECVD, STS). We used a spectroscopic ellipsometer (RC2, J.A.
Woollam Ellipsometry Solutions) to measure the film’s thickness (667 nm) and
its refractive index, which is shown in Figure S7. We then spun on a 90-nm-
thick layer of PMMA A2 (Kayaku Advanced Materials) and soft baked the sample
on a hotplate at 180 ∘C for 3 minutes. We used a thermal imprinting tool
(NX-2600BA, Nanonex) to transfer the pattern on the stamp to the PMMA layer.
The imprint process consists of a 2 minute pump; a 50 second pre-imprint phase
at 180 ∘C and 100 psi; a 1 minute, 30 second imprint at 210 ∘C and 200 psi;
and a cooldown phase, after which the pressure is released. Pressures and
temperatures during a typical imprint are shown in Figure S8. We carefully
separated the stamp and workpiece after removing them from the tool. We then
descummed the PMMA using oxygen plasma (Vision 320, Advanced Vacuum),
deposited a 30-nm-thick aluminum oxide hard mask using electron beam
evaporation, lifted off the PMMA in a solvent (Remover PG, Kayaku Advanced
Materials) and sonicated the sample in acetone. Subsequently, we transferred
the metalens pattern to the silicon nitride layer by dry etching in a mixture
of SF_6 and C_4F_8 gases (PlasmaPro 100 Cobra, Oxford Instruments). Finally,
we removed the aluminum oxide mask using a TMAH-based developer (AZ 300 MIF,
AZ Electronic Materials). We fabricated the control metalens using an electron
beam lithography process similar to the one for the master mold described
above, using PECVD silicon nitride substrates instead of the amorphous silicon
substrate. After developing the electron beam resist, we evaporated a 30-nm-
thick aluminum oxide hard mask, lifted it off, and etched the control metalens
similarly to the NIL metalenses.

Focal spot characterization. We obtained images of the focal spots using the
setup shown schematically in Figure 4a. The light from a 552 nm laser (OBIS
552-20LS, Coherent) was passed through a 10$\times$ beam expander (Melles
Griot, 09lbx003) and illuminated the metalenses. A custom-built microscope
formed of a 50$\times$ objective (LMPlanFL N, Olympus) and 200 mm focal length
tube lens (AC254-200-A-ML, Thorlabs) relayed a magnified image onto a CCD
camera (CoolSnap K4, Photometrics).

Focusing and deflection efficiency measurements. We used five lasers (532 nm:
COMPASS 215M-15, Coherent; 552 nm: OBIS 552-20LS, Coherent; 559 nm: DPSS,
Laserland Australia, 544 nm: HeNe 1676, JDS Uniphase; 633 nm: HeNe 1107P-2194,
JDS Uniphase) in our focusing efficiency measurement. The laser beams were
expanded using a 10$\times$ beam expander (Melles Griot, 09lbx003) and passed
through a 6-mm-diameter aperture to match the metalens aperture. To measure
the optical power focused by the metalens, we placed a 40-$\upmu$m-diameter
pinhole (P40D, Thorlabs) in the focal plane of the metalens and measured the
power passing through the pinhole using a power meter (PM100D with S120C head,
Thorlabs). We corrected measured power for the percentage of the power on an
ideal metalens blocked by the aperture (i.e., divided the measured value by
0.986). To determine the optical power incident on the metalens aperture, we
remove the metalens and the pinhole, measured the power passing through the
6-mm-diamter aperture divide the result by 0.965 to correct for the reflection
off the back side of the metalens substrate. Finally, the focusing efficiency
was obtained by dividing the power of focused light by the incident power. To
find the distribution of power inside the 40-$\upmu$m-diameter pinhole (Figure
S4a), we constructed a high dynamic range composite image of the focal spot by
combining many different exposures. We also measured the deflection efficiency
(ratio of power in the first diffraction order to the incident power) of the
metalens using the setup shown in Figure S5a: the light from a laser (559 nm,
DPSS, Laserland Australia) was polarized and was focused on the plane of the
metalens by passing through a 200 mm focal length lens (AC254-200-A-ML,
Thorlabs). The loosely focused light ensured that the diffraction orders were
well separated. The power of the first diffraction order and the incident beam
were measured using a power meter (PM100D with S120C head, Thorlabs) and the
incident beam power was divided by 0.965 to correct for the reflection off the
back side of the metalens substrate. Finally, the deflection efficiency was
obtained by diving the power of the first diffraction order by the incident
beam power.

Holographic measurement. We determined the phase of light at the metalens
surface using the interferometric imaging setup shown schematically in Figure
5a. We illuminated the metalens with laser light (552 nm, OBIS 552-20LS,
Coherent) and, using a custom-built microscope formed of a 10$\times$
objective (MPlanFL N, Olympus) and 200 mm focal length tube lens
(AC254-200-A-ML, Thorlabs), relayed an image of its surface onto a
transmissive diffraction grating with 300 grooves/mm (GT25-03, Thorlabs). A
second image relay comprising a 75 mm focal length lens (LA1145, Thorlabs) and
500 mm focal length lens (LA1380, Thorlabs) projects this image onto a camera
(CoolSnap K4, Photometrics). A 30-$\upmu$m-diameter pinhole (P30D, Thorlabs)
in the Fourier plane of the relay acts as a low pass filter on the image from
the first order, creating a clean reference beam for the interferogram. A
screen blocks all other orders of the grating except for the zeroth order. We
use a computational Fourier-domain filtering technique to retrieve the
complex-valued image [2]. We used a freely available image processing library
[3] to unwrap the phase of the complex image. This unwrapped phase is shown in
Figure 5b. The target phase for the metalens is given by
$\phi(r)=k_{0}\left(f^{2}+r^{2}\right)^{1/2}+\phi_{0}$, where $k_{0}$ is the
free-space wavenumber, $f$ is the focal length of the metalens at its design
wavelength, $r$ is the radial distance from the metalens center and $\phi_{0}$
is the phase offset at the center of the metalens. The data in Figures 5c and
5d show the difference of this target phase and the measured phase data.

Imaging. We measured the spectra of the LEDs used as illumination sources in
the imaging experiment (Figure 6b) using an optical spectrum analyzer (Ando
AQ-6315a). The lens used in the reference image in Figure 6c was an aspheric
lens with 5 mm clear aperture and an NA of 0.16 (C260TM, Thorlabs). All the
images were captured on a CMOS camera with a color sensor (DCC1240C,
Thorlabs). Figure S9 shows a photograph of the metalens on its camera mount.

## References

  * [1] Amir Arbabi, Ehsan Arbabi, Mahdad Mansouree, Seunghoon Han, Seyedeh Mahsa Kamali, Yu Horie, and Andrei Faraon.  Increasing efficiency of high numerical aperture metasurfaces using the grating averaging technique.  Sci. Rep., 10(1):1–10, 2020.
  * [2] Basanta Bhaduri, Chris Edwards, Hoa Pham, Renjie Zhou, Tan H Nguyen, Lynford L Goddard, and Gabriel Popescu.  Diffraction phase microscopy: principles and applications in materials and life sciences.  Adv. Opt. Photonics, 6(1):57–119, 2014.
  * [3] Stefan Van der Walt, Johannes L Schönberger, Juan Nunez-Iglesias, François Boulogne, Joshua D Warner, Neil Yager, Emmanuelle Gouillart, and Tony Yu.  scikit-image: image processing in python.  PeerJ, 2:e453, 2014.

Table S1: Focusing efficiencies for the NIL and EBL metalenses. Errors represent uncertainty due to fluctuation in laser power. Wavelength (nm) | Metalens | Focusing efficiency (%)
---|---|---
532 | NIL | 70$\pm$1
EBL | 77$\pm$1
543 | NIL | 73$\pm$1
EBL | 82$\pm$2
552 | NIL | 81$\pm$1
EBL | 89$\pm$1
559 | NIL | 80$\pm$2
EBL | 89$\pm$2
633 | NIL | 68$\pm$1
EBL | 74$\pm$1

Figure S1: High NA silicon nitride metalenses with high focusing efficiency.
(a) Schematic of the metalenses and the meta-atoms comprising the metalenses.
(b) The design curve that maps the desired phase to the nano-post width. The
minimum nano-post width and the gap between the nano-posts are 100 nm and 120
nm, respectively. (c) Simulated focusing efficiency values for metalenses with
different NAs that are designed using the meta-atoms shown in (a) and the
design curve shown in (b). The simulation is performed using the grating
averaging technique [1]. The dashed line is a fit to simulated data points
(blue circles).

Figure S2: Images of imprinted PMMA. (a) Optical microscope, and (b) SEM
images of imprinted PMMA.

Figure S3: (a) A photo of the master mold used for fabricating the NIL
metalenses. (b) Photos of the EBL metalens, the NIL metalens, and a second NIL
metalens (NIL metalens 2). The NIL metalens 2 was fabricated prior to fully
optimizing the silicon nitride etch process. (c) A scanning electron
microscope image of the nano-posts on the NIL metalens 2 that shows a sub-
optimal etch process. (d) Measured focal spot for the NIL metalens 2 and an
ideal metalens. The negligible aberrations of the focal spot indicate that the
stamp was not deformed significantly during the imprinting process. The
measured focusing efficiency of the NIL metalens 2 was 73%, indicating the
effect of the sub-optimal etch.

Figure S4: High dynamic range spot image. (a) Logarithmic scale plot of the
measured focal spot at $\lambda$ = 559 nm for NIL metalens and the simulated
focal spot for an ideal metalens with the same NA. Circle (black, dashed)
shows the size of the 40-$\mu$m-diameter pinhole. (b) Relative encircled power
for the NIL and ideal metalenses normalized to the power incident on the
metalenses.

Figure S5: Deflection efficiency. (a) Schematic illustration and a photograph
of the deflection efficiency measurement setup. (b) Optical microscope image
of the NIL metalens used for deflection efficiency measurement. The radial
distance and the locations on the metalens corresponding to deflection angles
of 5∘ and 10∘ are marked. (c) Measured deflection efficiency for p and s
polarized light. (d) Schematic illustration and a photograph of the focusing
efficiency measurement setup. Efficiency values are presented in Figure 4e.
QWP: Quarter-wave plate.

Figure S6: Sensitivity of metalens focusing efficiency to fabrication errors.
(a) The focusing efficiency of the metalens shown in Fig. 2a as a function of
systematic (deterministic) errors in the height, and (b) width of all nano-
posts. (c) The focusing efficiency of the metalens shown in Fig. 2a to random
errors in nano-post widths. The error distribution is assumed to have a
Gaussian distribution and the abscissa represents the standard deviation (STD)
of the distribution. The simulations are performed using the grating averaging
technique [1]. The dashed line are fits to simulated data points (blue
circles).

Figure S7: Refractive index of silicon nitride. The refractive index values
are measured using a spectroscopic ellipsometer.

Figure S8: Thermal imprint. Pressure (top) and temperature (bottom) as a
function of time.

Figure S9: Photograph of metalens camera.