# Frustrated hops in Ring Polymer Surface Hopping: Real-time dynamics and
detailed balance

Dil K. Limbu  Farnaz A. Shakib<EMAIL_ADDRESS>Department of Chemistry
and Environmental Science, New Jersey Institute of Technology, Newark 07102,
NJ United States

###### Abstract

Ring Polymer Surface-Hopping (RPSH) has been recently introduced as a well-
tailored method for incorporating nuclear quantum effects (NQEs), such as
zero-point energy and tunneling, into non-adiabatic molecular dynamics
simulations. The practical widespread usage of RPSH demands a comprehensive
benchmarking of different reaction regimes and conditions with equal emphasis
on demonstrating both the cons and pros of the method. Here, we investigate
the fundamental questions related to the conservation of energy and detailed
balance in the context of RPSH. Using Tully’s avoided crossing model as well
as a 2-level system coupled to a classical bath undergoing Langevin dynamics,
we probe the critical problem of the proper treatment of the classically
forbidden transitions stemming from the surface hopping algorithm. We show
that proper treatment of these frustrated hops is key to the accurate
description of real-time dynamics as well as reproducing the exact quantum
Boltzmann population.

††preprint: AIP/123-QED

##  I Introduction

Trajectory-based mixed quantum-classical dynamics methods have been a
mainstream approach in describing electronic and/or nuclear quantum effects in
condensed-phase dynamics.Brown and Shakib (2021); Crespo-Otero and Barbatti
(2018); Wang, Akimov, and Prezhdo (2016); Tully (2012); Barbatti (2011) For
more than three decades, Tully’s Fewest Switches Surface Hopping (FSSH)Tully
(1990) has served as ground zero for developing a myriad of such methods. In
the original FSSH, the classical degrees of freedom (DOFs) evolve on the
adiabatic surfaces interrupted by instantaneous electronic transitions. For
each nuclear trajectory, the associated electronic density matrix is
determined by coherently integrating the electronic time-dependent Schrödinger
equation along the classical trajectory. The transition probabilities between
adiabatic states are then determined from the density matrix elements and
their time derivatives. Hammes-Schiffer and Tully extended this methodology to
non-adiabatic transitions between vibrational states as well, which allowed
simulation of proton-transfer (PT)Hammes-Schiffer and Tully (1994) and proton-
coupled electron transfer (PCET)Fang and Hammes-Schiffer (1997a, b, 1999)
reactions. Efforts to remedy the primary limitation of the FSSH ansatz,
commonly known as overcoherence problem or lack of decoherence, gave birth to
more sophisticated trajectory-based methods such as augmented surface
hopping,Jain, Alguire, and Subotnik (2016) decoherence-induced surface
hopping,Jaeger, Fischer, and Prezhdo (2012) simultaneous-trajectory surface
hopping,Shenvi, Subotnik, and Yang (2011) and coherent fewest-switches
surface-hopping,Tempelaar and Reichman (2018) to name a few. The lack of a
formal derivation of surface hopping was met by the advent of a first-
principle trajectory-based methodology utilizing quantum-classical Liouville
equation (QCLE)Kapral and Ciccotti (1999); Nielsen, Kapral, and Ciccotti
(2000) which was proved to be very accurate for PTHanna and Kapral (2005)
reactions. Shakib and Hanna extended QCLE to the study of the rate and
mechanism of PCETShakib and Hanna (2014, 2016a, 2016b) reactions yielding
exceptionally accurate results, albeit at the expense of a huge computational
cost. The number of classical trajectories to reach convergence in the QCLE
approach for PCET reactions turned out to be at the order of millions,Shakib
and Hanna (2016b) which is 1000 times more than the FSSH requirement.
Nevertheless, later, Kapral showed that by applying two approximations,
surface-hopping could be derived from QCLEKapral (2016) hence putting the
former on much firmer ground than the previously phenomenologically perceived
one. Very recently, developments aimed at dealing with the limitations of FSSH
via a mapping approach to surface hopping (MASH)Mannouch and Rechardson
(2023), which imposes internal consistency between nuclear and electronic
DOFs. Still, MASH performs its best with the addition of decoherence; however,
this is not an ad hoc correction scheme but the result of the rigorous
derivation of the method at the QCLE limit.Mannouch and Rechardson (2023) In
2012, Shushkov et al. tried to remedy another inherent limitation of surface-
hopping,Shushkov, Li, and Tully (2012) i.e., lack of nuclear quantum effects
(NQEs), by marriage between surface hopping and ring polymer molecular
dynamics (RPMD)Craig and Manolopoulos (2004). In this method, the non-
adiabatic electronic transitions are described by the FSSH algorithm. At the
same time, NQEs are incorporated with Feynman’s imaginary-time path-integral
formalism,Feynman and Hibbs (1965) giving birth to ring polymer surface
hopping (RPSH). Thus, RPSH is deemed a well-tailored method for investigating
multi-electron/multi-proton transfer dynamics in condensed phases. Shakib and
Huo applied RPSH with centroid approximation to investigate electronic non-
adiabatic dynamics in three infamous Tully models with explicit nuclear
quantization.Shakib and Huo (2017) They showed that RPSH can qualitatively
capture the correct branching probabilities, especially at low-temperature
limits, due to the inclusion of nuclear tunneling and zero-point energy via
the extended phase-space of the classical ring polymer. Interestingly, RPSH
was also capable of quantitatively reproducing the branching probabilities in
the model of "extended coupling with reflections", which is specifically
designed for the study of coherence/decoherence.

Despite all the promises that RPSH offers, it is still an approximate method
whose accuracy and validity limits should be carefully examined. Specifically,
while we know that it combines the powerful capabilities of FSSH and RPMD, we
should also ask about the shortcomings of its constituent methodologies. Are
they magnified or quenched within RPSH? And to what degree? In the current
work, we address the critical issues related to surface hopping phenomenon. We
will focus on the conservation of energy, which leads to frustrated hops, and
their effect on preserving detailed balance. We show how these frustrated hops
are manifested in RPSH and offer different remedies for dealing with them. The
effect of such remedies is discussed in terms of different model systems to
provide valuable insights into the choice of proper recipe for treating
frustrated hops with respect to the objectives of a study. While not a single
treatment of frustrated hops seems applicable in all models, RPSH
systematically operates better than FSSH. This is on top of the inclusion of
NQEs into the dynamics that separates RPSH from the ever-expanding pool of
surface hopping methods. Further, this paper portrays a clear picture of the
functionality of RPSH, which is crucial for future developments via a unified
community effort. In the remainder of this paper, we will first give a summary
of the RPSH algorithm at the centroid level as well as methodological
considerations regarding frustrated hops. The studied models will be
introduced at the beginning of each subsection in the Results and Discussions.
Concluding Remarks and our Outlook toward future research directions in this
field will follow.

##  II Methodological Considerations

###  II.1 Ring Polymer Surface Hopping

We present a brief introduction to the RPSH ansatz with centroid
approximation, interested reader is referred to Refs. 23 and 26 for further
details. In RPSH, every nuclear DOF is represented by a ring polymer which is
comprised of n copies of the nuclear DOF, known as beads, connected by
harmonic forces. The corresponding extended Hamiltonian is described with:

$H_{n}=\sum^{n}_{i=1}\left[\frac{\textbf{P}_{i}^{2}}{2M}+\frac{M\omega_{n}^{2}}{2}(\textbf{R}_{i}-\textbf{R}_{i-1})^{2}+V_{\alpha}(\textbf{R}_{i})\right].$
(1)

Here, $n$ is the total number of beads and $\omega_{n}=n/\beta\hbar$ where
$\beta=1/k_{\rm{B}}T$ is the reciprocal temperature. $\textbf{P}_{i}$ and
$\textbf{R}_{i}$ represent the momentum and position of each bead of the ring
polymer which moves on a single adiabatic surface
$|\alpha;\textbf{R}_{i}\rangle$ corresponding to the potential energy
$V_{\alpha}(\textbf{R}_{i})=\langle\alpha;\textbf{R}_{i}|\hat{V}|\alpha;\textbf{R}_{i}\rangle$.
At every time step of the nuclear dynamics, the position and momentum of the
centroid of the ring polymer is updated as:

$\bar{\textbf{R}}=\frac{1}{n}\sum_{i=1}^{n}\textbf{R}_{i}\;\;\;\;\;\;\;\;\;\;\;\;\;\bar{\textbf{P}}=\frac{1}{n}\sum_{i=1}^{n}\textbf{P}_{i}$
(2)

At this point, according to the surface hopping algorithm, the time-dependent
Shrödinger equation is numerically integrated as

$i\hbar\dot{c}_{\alpha}(t)=V_{\alpha}(\bar{\textbf{R}})\,c_{\alpha}(t)-i\hbar\sum_{\beta}\dot{\bar{\textbf{R}}}\cdot\textbf{d}_{\alpha\beta}(\bar{\textbf{R}})\,c_{\beta}(t)$
(3)

to confer the electronic coefficients $c_{\alpha}$ associated with each
adiabatic surface. The important difference with the original FSSH is that
this integration is carried out along the motion of the centroid. Both the
energy of the adiabatic surfaces, $V_{\alpha}(\bar{\textbf{R}})$ and the non-
adiabatic coupling vector between surfaces,
$\textbf{d}_{\alpha\beta}(\bar{\textbf{R}})=\langle\alpha;\bar{\textbf{R}}|\nabla_{\bar{\textbf{R}}}|\beta;\bar{\textbf{R}}\rangle$,
are evaluated at the centroid level. The same goes for the probability of
transition, i.e., switching between surfaces at each time step $\Delta t$,
which is defined based on density matrix elements,
$\rho_{\alpha\beta}=c_{\alpha}c_{\beta}^{*}$, as

$g_{\alpha\beta}=\frac{-2\mathrm{Re}(\rho_{\beta\alpha}^{*}\dot{\bar{\textbf{R}}}\cdot\textbf{d}_{\beta\alpha}(\bar{\textbf{R}}))\,\Delta
t}{\rho_{\alpha\alpha}}.$ (4)

The nonadiabatic transition from the current surface $\alpha$ to the next
surface $\beta$ occurs if $g_{\alpha\beta}$ is greater than a randomly-
generated number between 0 and 1. If a transition occurs, the entire ring
polymer hops to the next adiabatic surface while the velocity of each bead in
the ring polymer is re-scaled in order to conserve the total energy of the
quantum plus classical subsystems.

###  II.2 Conservation of energy and frustrated hops

By construction, surface hopping methodologies conserve total quantum plus
classical energy by enforcing the hopping trajectories to have enough kinetic
energy to compensate for the potential energy difference between surfaces, or
states. Transition attempts that do not fulfill this requirement are deemed as
frustrated hops and return to the original states.Tully (1990) Proper
treatment of these frustrated hops is an essential step in carrying out
surface hopping simulations. Tully demonstrated frustrated hops as
trajectories hitting a wall and coming back, hence the component of their
velocity in the direction of nonadiabatic coupling vector needed to be
reversed.Tully (1990); Hammes-Schiffer and Tully (1994) Later, Müller and
Stock showed that not reversing the velocity of frustrated hops leads to
significantly improved results in photoinduced relaxation dynamics in
comparison to reversing the velocityMüller and Stock (1997) hinting that the
choice of treatment can be system-dependent. In our first implementation of
RPSH,Shakib and Huo (2017) we opted to not reverse the velocity of the beads
of the ring polymer in the case of frustrated hops. To assess the accuracy of
the method and shed light on manifestation of frustrated hops in RPSH, here,
we revisit Tully’s single avoided crossing model and establish the effect of
reversing or not reversing the velocity of frustrated hops on the adiabatic
population transfer profile in the low to medium momentum regions where the
nuclear quantum effects are paramount. Furthermore, we employ two different
remedies suggested by TruhlarJasper and Truhlar (2003) and SubotnikJain and
Subotnik (2015) for frustrated hop treatment. Jasper and Truhlar combined the
features of velocity "reversing" (VR) and "not reversing" (NR) in a
prescription called $\Delta V$ which, whenever a frustrated hop occurs, allows
the trajectory to feel the target adiabatic state $\beta$. Accordingly, in the
case of a frustrated hop in RPSH, we compare
$\dot{\bar{\textbf{R}}}\cdot\textbf{d}_{\alpha\beta}$ to the component of the
force in the direction of the nonadiabatic coupling vector, defined as

$F_{\beta}=-\nabla
V_{\beta}({\bar{\textbf{R}}})\cdot\textbf{d}_{\alpha\beta}(\bar{\textbf{R}}).$
(5)

If these two quantities have the same sign, we do not reverse the velocity but
will do so if they have opposite signs. On the other hand, Jain and Subotnik
suggested a comparison between the components of forces from both the current
and the target adiabatic states. Accordingly, in RPSH, we calculate:

$\bigg{(}\textbf{d}_{\alpha\beta}(\bar{\textbf{R}})\cdot\big{(}-\nabla
V_{\alpha}({\bar{\textbf{R}}})\big{)}\bigg{)}\bigg{(}\textbf{d}_{\alpha\beta}(\bar{\textbf{R}})\cdot\big{(}-\nabla
V_{\beta}({\bar{\textbf{R}}})\big{)}\bigg{)}$ (6)

and reverse the velocity if it is smaller than zero. Henceforth, we shall call
this approach $\Delta V^{2}$ as it bears resemblance to the $\Delta V$
approach.

###  II.3 Frustrated hops and detailed balance

The requirement to conserve energy leads to disruption of internal consistency
despite the original promise of surface hopping algorithm, i.e., the fraction
of trajectories ended up on each state is not consistent with the average
quantum probabilities anymore. Even though, frustrated hops are not the
Achilles heel of the surface hopping methodologies, note that transitions from
higher-energy to lower-energy states are always permitted and it is the
reverse that encounters rejection. As a result, presence of frustrated hops
leads to an almost correct apportionment of nuclear DOFs according to the
Boltzmann distribution. It was shown early-on that in a 2-state system coupled
to a chain of N classical DOFs, the original FSSH algorithm approximately, in
some reaction regimes exactly, preserves detailed balance.Parandekar and Tully
(2005); Schmidt, Parandekar, and Tully (2008) Recently, Prezhdo and coworkers
revisited this study with a focus on the relationship between detailed balance
and proper treatment of frustrated hops.Sifain, Wang, and Prezhdo (2016) They
showed that the VR approach slightly improves the detailed balance in this
model system compared to the NR approach over a high-temperature region of 350
K-2500 K. This conclusion is, of course, opposite to the earlier findings of
Müller and StockMüller and Stock (1997) albeit in a different system. Here, we
will employ a slightly modified version of the chain model to investigate the
capability of RPSH in preserving detailed balance, with different treatments
of frustrated hops, compared to the original FSSH. Furthermore, we will
demonstrate our results at the significantly important region of low
temperatures.

##  III Results and Discussion

###  III.1 Single-avoided crossing model: Branching probabilities

Here, we focus on Tully’s single-avoided crossing (SAC) model, shown in the
inset of Fig. 1, which is defined by a $2\times 2$ diabatic matrix with the
details being provided in the Supporting Information (SI). Diagonalizing this
matrix yields two adiabatic states with a single avoided crossing and a
corresponding non-adiabatic coupling vector. This 2-state system is then
coupled to one nuclear DOF whose wavefunction is initialized on state
$|1\rangle$ as a Gaussian wavepacket in the form of

$G(R)=\bigg{(}\frac{2\alpha}{\pi}\bigg{)}^{\frac{1}{4}}\mathrm{exp}(-\alpha(R-R_{0})^{2}+ik(R-R_{0})).$
(7)

This is manifested in a Gaussian distribution of nuclear position around
$R_{0}=-15$ a.u. with the width of $\sigma_{R}=1/\sqrt{2\alpha}$. Here,
$\alpha=0.25$ a.u. and $k$ is the deterministic incoming momentum of the
nuclear DOF with the mass $M=2000$ a.u. Note that, here, the classical DOF is
expanded to a ring polymer with four beads whose initial positions are
randomly selected from this Gaussian distribution. The choice of four beads
for the ring polymer is based on the convergence of the adiabatic
populations.Shakib and Huo (2017) We first focus on the branching
probabilities as a function of the incoming momentum obtained from an ensemble
of 10,000 RPSH trajectories ran for 0.5 ps with a time-step of 1 a.u.
($\sim$0.024 fs). Branching probabilities show the final adiabatic populations
of states $|1\rangle$ and $|2\rangle$ with differentiating between
trajectories that upon reaching at avoided crossing transmit from the reactant
to the product side or reflect back to the reactant side. Fig. 1 demonstrate
the branching probabilities, namely transmission and reflection on state
$|1\rangle$ ($T_{1}$ and $R_{1}$, respectively), in a low-momentum region,
i.e. $k=2-6.5$ a.u.

Figure 1: Branching probabilities in the Tully’s single avoided crossing (SAC)
model at the low momentum region; Transmission (a) and reflection (b) on state
$|1\rangle$, T1 and R1, respectively. The model is shown in the in-set with
the red line representing the non-adiabatic coupling vector.

According to the exact results shown in dashed lines, in this region, the
incoming momentum is not sufficient for non-adiabatic transition and the
wavepacket moves adiabatically on state $|1\rangle$. This dynamics is governed
by nuclear tunneling where exact results show a smooth increase of $T_{1}$ in
the expense of a decrease in $R_{1}$.

Figure 2: Branching probabilities in the SAC model at the intermediate
momentum region; transmission on state $|1\rangle$, $T_{1}$, (a), transmission
on state $|2\rangle$, $T_{2}$, (b), and reflection on state $|1\rangle$,
$R_{1}$ (c).

We have previously shownShakib and Huo (2017) how FSSH fails to capture this
quantum-mechanical phenomenon but RPSH recovers the correct physical
behaviour, i.e. a smooth transition between $R_{1}$ and $T_{1}$, through the
extended phase-space of a classical ring polymer. Here, we evaluate the
effects of the proper treatment of frustrated hops on branching probabilities
applying four different schemes for adjusting their velocities, namely
velocity reversing (VR) or not reversing (NR) along with $\Delta V$ and
$\Delta V^{2}$ approachesJasper and Truhlar (2003); Jain and Subotnik (2015)
as explained in the previous section. As can be seen in Fig 1, the four
schemes show similar behaviour in capturing the smooth transition in the very
low momentum region but start deviating around $k=3.3$ a.u. where only the NR
scheme follows the same behaviour of the exact results in 100% reversal of
$R_{1}$ to $T_{1}$. The other three schemes underestimate the magnitude of
adiabatic transmission at higher values of momentum. The reason goes back to
the fact that trajectories facing frustrated hops in NR scheme still continue
their motion with positive momentum. They tunnel through the energy barrier on
state $|1\rangle$ and end up in the product side. On the other hand, fully or
partially reversing the velocity in the other three schemes forces such a
trajectory to go back to the reactant side negating the possibility of
tunneling through the energy barrier. This results in an artificial increase
of $R_{1}$ in the expense of decrease in $T_{1}$. As we increase the momentum
beyond $k=6.5$ a.u. we reach an intermediate reaction regime where non-
adiabatic trajectories appear alongside adiabatic trajectories. This higher
momentum results on successful transitions from state $|1\rangle$ to
$|2\rangle$ and a smooth increase of $T_{2}$ as shown in Fig. 2b. As can be
seen, all four schemes retrieve the correct physical behaviour due to the
inclusion of ZPE in the dynamics by RPSH. We have already shown that FSSH
fails to do so.Shakib and Huo (2017) On the other hand, there is still
significant deviation between NR and the other three schemes in retrieving the
correct physical behaviour in the case of $T_{1}$ and $R_{1}$ with only the
former being successful. To shed more light on this finding, we take a look at
the number of frustrated hops averaged over the number of trajectories
($N_{FHs}$) as a function of the incoming momentum in the RPSH method, Fig. 3.
A comparison between these results and Fig. 1 and Fig. 2 show that the
deviation between the four schemes starts as soon as we have frustrated hops
in our simulations, i.e. around $k=3.3$ a.u., and persists till around $k=8.9$
a.u. where $N_{FHs}$ goes to zero.

Figure 3: The number of frustrated hops at different velocity re-scaling
schemes in RPSH in the low and intermediate momentum region of the SAC model.

Clearly, in the adiabatic or intermediate reaction regimes RPSH trajectories
encounter frustrated hops and even partially reversing the velocity can
prevent an otherwise successful trajectory dynamics from taking place. Another
point seen in Fig. 3 is the lower number of frustrated hops in the NR scheme
compared to the other three schemes. The reason can be traced back to the fact
that in the NR scheme, successful trajectories continue moving to the product
side regardless of how many times they encounter frustrated hops or whether
this encounter happens in the vicinity of an energy barrier on the reactant or
the product side. For the other three schemes, on the other hand, a trajectory
on the product side might encounter a frustrated hop and due to reversing the
velocity tunnels back to the reactant side. Now, if that trajectory encounters
another reversing of velocity it will go back to the product side. Overall,
reversing the velocity allows the trajectories to spend more time on the
avoided crossing region and hence $N_{FHs}$ increases. This can enhance the
deviation of the results of velocity re-scaling schemes from the exact
results.

###  III.2 N-particle chain model: Detailed balance

The original RPMD, by construction, yields real-time MD trajectories that
preserve exact quantum Boltzmann distribution.Craig and Manolopoulos (2004)
Its expansion to multi-state non-adiabatic dynamics, however, did not always
lead to preserving the detailed balance mainly due to the zero-point energy
leakage. Nevertheless, multi-variable (MV)-RPMD, where both electronic states
and nuclear DOFs are represented by continuous Cartesian coordinates,Ananth
(2013) demonstrated that ring polymer isomorphism can exactly preserve
detailed balance in multi-state systems. On the other hand, as was explained
before, FSSH approximately preserves detailed balance due to the presence of
frustrated hops. Now, the question is whether RPSH can borrow these desired
characteristics and improve upon. To investigate this matter, here, we employ
a 2-state system, with constant values of the energy gap
($\Delta=E_{2}-E_{1}$) and non-adiabatic coupling vector ($d_{12}$), coupled
to a chain of $N$ nuclear DOFs.Parandekar and Tully (2005) The quantum
subsystem is coupled to only the first particle in the chain and the nearest-
neighbor potential energy between the particles is a quartic Morse potential
as

$V(\mathbf{R})=\sum_{i=1}^{N}V_{M}(R_{i}-R_{i+1})$ (8)

where

$V_{M}(R)=V_{0}(a^{2}R^{2}-a^{3}R^{3}+0.58a^{4}R^{4}).$ (9)

The atom farthest from the quantum subsystem ($R_{N}$) is connected to
Langevin dynamics. All the particles in the chain are represented by ring
polymers comprised of four beads based on our convergence tests, see the SI
Figure S1. Subsequently, the EOM for the last particle would be slightly
different from the conventional classical Langevin dynamics. Each mode of a
free ring polymer ($k$) exhibits uncoupled harmonic oscillator dynamics when
represented in the normal mode basis.Ceriotti  _et al._ (2010) Accordingly,
the corresponding EOM for the last ring polymer is written as

$\dot{P}_{N}^{(k)}=-\frac{\partial V}{\partial
R_{N}^{(k)}}-\gamma^{(k)}P_{N}^{(k)}+\sqrt{\frac{2M\gamma^{(k)}}{\beta_{n}\Delta
t}}.$ (10)

where $\beta_{n}=\beta/n$.

Table 1:  Simulation parameters used for the $N$-particle chain model. Parmeter | Value | Unit
---|---|---
$N$ | 20 |
$m$ | 12.0 | amu
$V_{0}$ | 175.0 | kJ/mol
$a$ | 4.0 | Å-1
$\gamma$ | $10^{14}$ | $s^{-1}$
$\Delta$ | 8.0 | kJ/mol
$d_{12}$ |  $-$6.0 | Å-1

The second term includes the Langevin friction constant $\gamma^{(k)}$. For
the excited modes of the ring polymer $(k>0)$ $\gamma^{(k)}=2\omega_{k}$ where
$\omega_{k}=2\sin(2\pi/n)/\beta_{n}\hbar$ and for the centroid mode
$\gamma^{(0)}=\gamma$. The last term in Eq. 10 represents a random force the
width of which depends both on the number of beads within $\beta_{n}$ and the
time step of the simulation $\Delta t$.

Figure 4: Equilibrium populations on the first state $P_{1}$ (a), and on the
second state $P_{2}$ (b) of the 2-state quantum subsystem in the chain model
obtained from RPSH ($\bullet$) and FSSH ($\Diamondblack$) methods with NR
scheme. Exact Boltzmann populations are shown in solid lines. Panel (c) shows
the relative errors of RPSH and FSSH in producing $P_{2}$ with respect to the
Boltzmann population.

All parameters of our 2-state chain model are listed in Table 1 with some
variations compared to the original model.Parandekar and Tully (2005)
Specifically, we changed the energy gap to 8.0 kJ/mol and the nonadiabatic
coupling vector to -6.0 Å-1 to expand our dynamics simulations to a wider
range of temperatures at 200 K$-$2500 K, with a special emphasis on the low
temperature region. Fig. 4a and 4b show the populations of the two states
($P_{1}$ and $P_{2}$) over a wide range of temperatures with the solid line
representing the exact Boltzmann populations ($P_{\rm{B}}$). The RPSH
populations are obtained from 1000 trajectories for 50 ps with $\Delta t=0.01$
fs and taking the average of the last 20 ps to avoid any dependence on the
initial conditions. FSSH results are obtained similarly and are shown for
comparison. The RPSH results are very satisfactory but at the same time are
very close to the FSSH data. For a better comparison, Fig. 4c demonstrates the
relative error of the RPSH and FSSH results for the second state defined as
$(P_{2}-P_{\rm{B}})/P_{\rm{B}}$. As can be seen, at high temperatures RPSH and
FSSH show identical results with the lowest relative error. As the temperature
is decreased, the error of both methods increases but RPSH acts much better
than FSSH in preserving detailed balance. Here, we used the NR scheme to
obtain these results following the observations in the previous section.
However, it should be noted that Prezhdo and coworkers earlier used a similar
model to show that VR scheme has a slightly improving effect on preserving
detailed balance with FSSH in a temperature range of 350 K$-$2500 K.Sifain,
Wang, and Prezhdo (2016) In Fig. 5, we demonstrate the relative error of the
four different velocity re-scaling schemes in preserving detailed balance in
RPSH. Among the four schemes, the VR scheme indeed leads to almost perfect
preservation of detailed balance in both high and low temperatures.
Nevertheless, regardless of the employed velocity re-scaling scheme, RPSH
always operates better than FSSH at low temperature limit. This can be traced
back to the natural preservation of the quantum Boltzman distribution in RPMD
which is manifested even in the multi-state dynamics here.

Figure 5: Relative errors as defined before in RPSH results with four
different velocity re-scaling schemes.

###  III.3 Differentiating between the two models

We showed in the case of the SAC model that reversing the velocity of the
adiabatic RPSH trajectories prevents an otherwise successful trajectory from
going from the reactant to the product side. This observation can be
generalized to any situation where tunneling through the energy barrier is a
prominent factor in the success of adiabatic reactions. It should be noted
that one can enforce the correct dynamics even with total or partial velocity
reversing in such systems by not allowing the trajectories to attempt a
transition. However, that needs a priori information about the reaction
regimes of different systems. On the other hand, regardless of the velocity
re-scaling scheme, RPSH is more successful than FSSH in preserving detailed
balance in the case of the 2-state quantum subsystem coupled to a linear chain
of 20 nuclear DOFs. The reason can be traced back to two factors. First, RPMD
technically is a phase-space representation for the quantum Boltzmann
distribution. Second, this extended phase-space, with a wider distribution of
momentum, leads generally to a higher number of frustrated hops in RPSH
compared to FSSH, see Fig. 6a. And, as we know the existence of frustrated
hops is the reason behind approximate preservation of the detailed balance in
surface hopping methods. It is noted that $N_{FHs}$ of the RPSH method is
similarly higher than FSSH in the case of the SAC model, Fig. 6b.

Figure 6: Number of frustrated hops averaged over number of trajectories in
the chain model over a range of temperatures (a) and in the SAC model over a
range of initial momentum (b) obtained from RPSH and FSSH methods.

Then, why the real-time dynamics in the SAC model is more accurately recovered
by not reversing the velocity of the frustrated hops in RPSH trajectories
while reversing it leads to exactly preserving the detailed balance in the
chain model, as was already shown for the FSSH case as well ?Sifain, Wang, and
Prezhdo (2016) One should note the essential differences between the two
models with the lack of an avoided crossing in the chain model being the most
important one. The chain model is designed to be a representative of a quantum
system in thermal equilibrium with its surrounding. Furthermore, it is not
suitable for investigating scattering events. Hence, the important observation
here is the superiority of RPSH in preserving detailed balance with the
additional insight that one can fine tune the results by choosing an
appropriate treatment for frustrated hops. On the other hand, real-time
dynamics of multi-state systems is best recovered by not interfering with the
role of NQEs throughout the dynamics. This is of course not the case for other
surface hopping methods that neglect NQEs. Nevertheless, care should be taken
in the choice of an appropriate treatment of frustrated hops with respect to
the system of interest.

##  IV Conclusions and Outlooks

In this work, Ring Polymer Surface Hopping (RPSH) non-adiabatic molecular
dynamics methodology is critically revisited in order to shed light on its
accuracy and validity limits. The most exciting aspect of RPSH, i.e., the
inclusion of nuclear quantum effects into molecular dynamics, can be further
harvested if the cons and pros of the method are exposed and documented in
different systems and reaction regimes. RPSH, similar to other surface hopping
algorithms, conserves the total quantum plus classical energy at the level of
each individual trajectory. This leads to the creation of frustrated hops
that, although they disrupt the internal consistency of the method but lead to
the conservation of exact quantum Boltzmann distribution. We showed that RPSH
yields better results in preserving detailed balance than the original FSSH
method in the low-temperature region using a 2-state system coupled to
Langevin dynamics via a linear chain of 20 nuclear DOFs. We also showed that
different schemes of re-scaling velocity have a profound effect on preserving
correct quantum Boltzmann populations where reversing the velocity after each
frustrated hop leads to almost exact conservation of detailed balance in a
range of low to high temperatures. This was opposite to what we observed in
calculating the branching probabilities in Tully’s avoided crossing model.
While we have more frustrated hops in RPSH than FSSH for this model, we also
have an interplay between frustrated hops and nuclear tunneling in the
adiabatic reaction regime. Reversing the velocity of frustrated hops can
disrupts the motion of an adiabatic trajectory that otherwise would tunnel
through a reaction barrier from the reactant to the product side. The correct
real-time dynamics in the case of this avoided crossing model is preserved
better by not reversing the velocity of frustrated hops. Hence, we advise
caution in choosing the proper treatment of frustrated hops in RPSH where one
needs to consider the nature of the problem to be studied.

To further bring RPSH to the mainstream, research should be directed on (i) a
systematic implementation of decoherence correction and (ii) ,if necessary,
reformulating the algorithm for the proper treatment of internal consistency.
Surface-hopping methods take care of electronic coherence in an elegant way.
However, they have always been plagued with the over-coherence problem or
neglect of decoherence.Subotnik and Shenvi (2011) We have already shown that,
in contrast to FSSH, RPSH can capture distribution-dependent decoherence in
Tully’s "extended coupling with reflection" model even with a deterministic
initial momentum.Shakib and Huo (2017) This improvement is because, despite
initial deterministic momentum, each bead’s momentum gradually differs through
the dynamical propagation, resulting in a broader centroid momentum
distribution compared to the distribution in FSSH. Therefore, different
reflected ring polymer trajectories gain different phases, and over an
ensemble of the ring polymer trajectories, the high-frequency oscillations in
reflection coefficients cancel out. However, we cannot yet claim this is a
general effect. In fact, the original FSSH can also enjoy similar improvements
if the initial conditions of the MD simulations involve a Maxwell-Boltzmann
distribution instead of deterministic momenta. Application of RPSH in more
sophisticated models designed for investigation of decoherence as well as
equipping it with different decoherence correction schemes is currently being
pursued in our group. On the other hand, as briefly mentioned in the
methodological considerations, surface hopping methods violate internal
consistency and as a result, populations of electronic states can be different
if they are calculated based on the number of nuclear trajectories on each
surface or the electronic amplitudes of those states or a mixture of
both.Landry, Falk, and Subotnik (2013) Recently, there has been a growing
argument on the importance of dealing with this issue rather than "decoherence
correction“ to systematically improve the surface hopping results.Mannouch and
Rechardson (2023) A part of the research in our group focuses on this issue in
the context of the real-time population transfer dynamics of a couple of
three-state Morse potential model systems that are designed for investigating
photo-excited relaxation. As can be seen in the SI Fig. S2, RPSH yields time-
dependent diabatic state populations that are in good agreement with the exact
quantum dynamics simulations for all three models. However, the three
approaches explained earlier give different population profiles. On the other
hand, according to our observations, no frustrated hops are recorded for these
simulations. Hence, these models provide a unique case for disentangling the
concepts of "internal consistency" and "decoherence correction" in the context
of the RPSH methodology.

###### Acknowledgements.

The authors acknowledge support from the New Jersey Institute of Technology
(NJIT). This research has been enabled by the use of computing resources and
technical support provided by the HPC center at NJIT. This work partially used
Bridges2 at Pittsburgh Supercomputing Center through allocation CHE200007 from
the Extreme Science and Engineering Discovery Environment (XSEDE), which was
supported by National Science Foundation Grant no. 1548562.Towns  _et al._
(2014)

## Conflict of Interest

The authors have no conflicts to disclose.

## Data Availability Statement

The data that support the findings of this study are available within the
article and its supplementary material.

## References

  * Brown and Shakib (2021) S. E. Brown and F. A. Shakib, “Recent progress in approximate quantum dynamics methods for the study of proton-coupled electron transfer reactions,” Phys. Chem. Chem. Phys. 23, 2535–2556 (2021).
  * Crespo-Otero and Barbatti (2018) R. Crespo-Otero and M. Barbatti, “Recent advances and perspectives on nonadiabatic mixed quantum–classical dynamics,” Chem. Rev. 118, 7026–7068 (2018).
  * Wang, Akimov, and Prezhdo (2016) L. Wang, A. Akimov, and O. V. Prezhdo, “Recent progress in surface hopping: 2011–2015,” J. Phys. Chem. Lett. 7, 2100–2112 (2016).
  * Tully (2012) J. C. Tully, “Perspective: Nonadiabatic dynamics theory,” J. Chem. Phys. 137, 22A301 (2012).
  * Barbatti (2011) M. Barbatti, “Nonadiabatic dynamics with trajectory surface hopping method,” WIREs Comput. Mol. Sci. 1, 620–633 (2011).
  * Tully (1990) J. C. Tully, “Molecular dynamics with electronic transitions,” J. Chem. Phys. 93, 1061–1071 (1990).
  * Hammes-Schiffer and Tully (1994) S. Hammes-Schiffer and J. C. Tully, “Proton transfer in solution: Molecular dynamics with quantum transitions,” J. Chem. Phys. 101, 4657–4667 (1994).
  * Fang and Hammes-Schiffer (1997a) J.-Y. Fang and S. Hammes-Schiffer, “Proton-coupled electron transfer reactions in solution: Molecular dynamics with quantum transitions for model systems,” J. Chem. Phys. 106, 8442–8454 (1997a).
  * Fang and Hammes-Schiffer (1997b) J.-Y. Fang and S. Hammes-Schiffer, “Nonadiabatic dynamics for processes involving multiple avoided curve crossings: Double proton transfer and proton-coupled electron transfer reactions,” J. Chem. Phys. 107, 8933–8939 (1997b).
  * Fang and Hammes-Schiffer (1999) J.-Y. Fang and S. Hammes-Schiffer, “Improvement of the internal consistency in trajectory surface hopping,” J. Phys. Chem. A 103, 9399–9407 (1999).
  * Jain, Alguire, and Subotnik (2016) A. Jain, E. Alguire, and J. E. Subotnik, “An efficient, augmented surface hopping algorithm that includes decoherence for use in large-scale simulations,” J. Chem. Theory Comput. 12, 5256–5268 (2016).
  * Jaeger, Fischer, and Prezhdo (2012) H. M. Jaeger, S. Fischer, and O. V. Prezhdo, “Decoherence-induced surface hopping,” J. Chem. Phys. 137, 22A545 (2012).
  * Shenvi, Subotnik, and Yang (2011) N. Shenvi, J. E. Subotnik, and W. Yang, “Simultaneous-trajectory surface hopping: A parameter-free algorithm for implementing decoherence in nonadiabatic dynamics,” J. Chem. Phys. 134, 144102 (2011).
  * Tempelaar and Reichman (2018) R. Tempelaar and D. R. Reichman, “Generalization of fewest-switches surface hopping for coherences,” J. Chem. Phys. 148, 102309 (2018).
  * Kapral and Ciccotti (1999) R. Kapral and G. Ciccotti, “Mixed quantum-classical dynamics,” J. Chem. Phys. 110, 8919 (1999).
  * Nielsen, Kapral, and Ciccotti (2000) S. Nielsen, R. Kapral, and G. Ciccotti, “Non-adiabatic dynamics in mixed quantum-classical systems,” J. Stat. Phys. 101, 225–242 (2000).
  * Hanna and Kapral (2005) G. Hanna and R. Kapral, “Quantum-classical liouville dynamics of nonadiabatic proton transfer,” J. Chem. Phys. 122, 244505 (2005).
  * Shakib and Hanna (2014) F. A. Shakib and G. Hanna, “An analysis of model proton-coupled electron transfer reactions via the mixed quantum-classical Liouville approach,” J. Chem. Phys. 141, 044122 (2014).
  * Shakib and Hanna (2016a) F. A. Shakib and G. Hanna, “New insights into the nonadiabatic state population dynamics of model proton-coupled electron transfer reactions from the mixed quantum-classical Liouville approach,” J. Chem. Phys. 144, 024110 (2016a).
  * Shakib and Hanna (2016b) F. Shakib and G. Hanna, “Mixed quantum-classical Liouville approach for calculating proton-coupled electron-transfer rate constants,” J. Chem. Theory Comput. 12, 3020–3029 (2016b).
  * Kapral (2016) R. Kapral, “Surface hopping from the perspective of quantum–classical liouville dynamics,” Chem. Phys. 481, 77–83 (2016).
  * Mannouch and Rechardson (2023) J. R. Mannouch and J. O. Rechardson, “A mapping approach to surface hopping,” J. Chem. Phys. 158, 104111 (2023).
  * Shushkov, Li, and Tully (2012) P. Shushkov, R. Li, and J. C. Tully, “Ring polymer molecular dynamics with surface hopping,” J. Chem. Phys. 137, 22A549 (2012).
  * Craig and Manolopoulos (2004) I. R. Craig and D. E. Manolopoulos, “Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics,” J. Chem. Phys. 121, 3368–3373 (2004).
  * Feynman and Hibbs (1965) R. P. Feynman and A. R. Hibbs,  _Quantum mechanics and path integrals_ (McGraw-Hill, New York, 1965).
  * Shakib and Huo (2017) F. A. Shakib and P. Huo, “Ring polymer surface hopping: Incorporating nuclear quantum effects into nonadiabatic molecular dynamics simulations,” J. Phys. Chem. Lett. 8, 3073–3080 (2017).
  * Müller and Stock (1997) U. Müller and G. Stock, “Surface-hopping modeling of photoinduced relaxation dynamics on coupled potential-energy surfaces,” J. Chem. Phys. 107, 6230–6245 (1997).
  * Jasper and Truhlar (2003) A. W. Jasper and D. G. Truhlar, “Improved treatment of momentum at classically forbidden electronic transitions in trajectory surface hopping calculations,” Chem. Phys. Lett. 369, 60–67 (2003).
  * Jain and Subotnik (2015) A. Jain and J. E. Subotnik, “Surface hopping, transition state theory, and decoherence. ii. thermal rate constants and detailed balance,” J. Chem. Phys. 143, 134107 (2015).
  * Parandekar and Tully (2005) P. V. Parandekar and J. C. Tully, “Mixed quantum-classical equilibrium,” J. Chem. Phys. 122, 094102 (2005).
  * Schmidt, Parandekar, and Tully (2008) J. R. Schmidt, P. V. Parandekar, and J. C. Tully, “Mixed quantum-classical equilibrium: Surface hopping,” J. Chem. Phys. 129, 044104 (2008).
  * Sifain, Wang, and Prezhdo (2016) A. E. Sifain, L. Wang, and O. V. Prezhdo, “Communication: Proper treatment of classically forbidden electronic transitions significantly improves detailed balance in surface hopping,” J. Chem. Phys. 144, 211102 (2016).
  * Ananth (2013) N. Ananth, “Mapping variable ring polymer molecular dynamics: A path-integral based method for nonadiabatic processes,” J. Chem. Phys. 139, 124102 (2013).
  * Ceriotti  _et al._ (2010) M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopoulos, “Efficient stochastic thermostatting of path integral molecular dynamics,” J. Chem. Phys. 133, 124104 (2010).
  * Subotnik and Shenvi (2011) J. E. Subotnik and N. Shenvi, “Decoherence and surface hopping: When can averaging over initial conditions help capture the effects of wave packet separation?” J. Chem. Phys. 134, 244114 (2011).
  * Landry, Falk, and Subotnik (2013) R. Landry, M. J. Falk, and J. E. Subotnik, “Communication: The correct interpretation of surface hopping trajectories: How to calculate electronic properties,” J. Chem. Phys. 139, 211101 (2013).
  * Towns  _et al._ (2014) J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D. Peterson, R. Roskies, and J. R. Scott, “Xsede: Accelerating scientific discovery,” Comput. Sci. Eng. 16, 62–74 (2014).
  * Abedi, Maitra, and Gross (2010) A. Abedi, N. T. Maitra, and E. K. U. Gross, “Exact factorization of the time-dependent electron-nuclear wave function,” Phys. Rev. Lett. 105, 123002 (2010).
  * Abedi, Maitra, and Gross (2012) A. Abedi, N. T. Maitra, and E. K. U. Gross, “Correlated electron-nuclear dynamics: Exact factorization of the molecular wavefunction,” J. Chem. Phys. 137, 22A530 (2012).
  * Ha and Min (2022) J.-K. Ha and S. K. Min, “Independent trajectory mixed quantum-classical approaches based on the exact factorization,” J. Chem. Phys. 156, 174109 (2022).


# Hyper Resolute Ultra thin Low Cost All-Dielectric Broadband Achromatic
Metalenses

Giuseppe Emanuele Lio  CNR Nanotec-Institute of Nanotechnology, UOS Cosenza,
87036 Rende (CS), 87036, Italy  Physics Department, University of Calabria,
I-87036 Arcavacata di Rende (CS), Italy  Antonio Ferraro  E-mail:
<EMAIL_ADDRESS>Physics Department, University of Calabria, I-87036
Arcavacata di Rende (CS), Italy  CNR Nanotec-Institute of Nanotechnology, UOS
Cosenza, 87036 Rende (CS), 87036, Italy  Tiziana Ritacco  CNR Nanotec-
Institute of Nanotechnology, UOS Cosenza, 87036 Rende (CS), 87036, Italy
Physics Department, University of Calabria, I-87036 Arcavacata di Rende (CS),
Italy  Dante Maria Aceti  Physics Department, University of Calabria, I-87036
Arcavacata di Rende (CS), Italy  Antonio De Luca  Physics Department,
University of Calabria, I-87036 Arcavacata di Rende (CS), Italy  CNR Nanotec-
Institute of Nanotechnology, UOS Cosenza, 87036 Rende (CS), 87036, Italy
Roberto Caputo  E-mail<EMAIL_ADDRESS>Physics Department,
University of Calabria, I-87036 Arcavacata di Rende (CS), Italy  CNR Nanotec-
Institute of Nanotechnology, UOS Cosenza, 87036 Rende (CS), 87036, Italy
Institute of Fundamental and Frontier Sciences, University of Electronic
Science and Technology of China, Chengdu 610054  Michele Giocondo  CNR
Nanotec-Institute of Nanotechnology, UOS Cosenza, 87036 Rende (CS), 87036,
Italy

###### Abstract

Metalenses offer the ground-breaking opportunity to realize highly performing
low-weight, flat and ultrathin, optical elements which substantially reduce
size and complexity of imaging systems. Today, a major challenge in metalenses
design is still the realization of achromatic optical elements ideally
focussing a broad wavelength spectrum at a single focal length. Here we
present a fast and effective way to design and fabricate extremely thin all-
dielectric metalenses, optimally solving achromaticity issues by means of
machine learning codes. The enabling technology for fabrication is a recently
developed hyper resolute two-photon direct laser writing lithography
equipment. The fabricated metalenses, based on a completely flat and ultrathin
design, show intriguing optical features. Overall, achromatic behavior, focal
length of 1.14 mm, depth of focus of hundreds of microns and thickness of only
few nanometers allow considering the design of novel and efficient imaging
systems in a completely new perspective.

##  I Introduction

The development of optical elements with nanoscale features has recently
raised a very huge and wide interest in materials science Genevet  _et al._
(2017). These research efforts are motivated by fundamental physical insights
demonstrating the possibility to realize perfectly flat lenses with eventual
application in a plethora of optical systems including integral imaging, phase
filters, discrete lenses for sensing, high numerical aperture and achromatic
optical devices, optical tunability or reconfigurability, phase dispersion
correction and correlation and augmented-reality (AR) applications Chen  _et
al._ (2018, 2019); Fan  _et al._ (2019); Idjadi and Aflatouni (2020); Gao  _et
al._ (2013); Yu and Capasso (2014); Byrnes  _et al._ (2016); Khorasaninejad
_et al._ (2017); Khorasaninejad and Capasso (2017a); Cui, Bai, and Sun (2019);
Chen, Zhu, and Capasso (2020). Initially, studies considered metasurfaces
where linear gradients of phase discontinuities lead to planar reflected and
refracted wavefronts. These metasurfaces, mainly working in the IR spectrum,
were characterized by discrete patterns involved to mime the functionality of
classical refractive lenses thus reconstructing the phase shift typically
introduced by the surface curvature Khorasaninejad  _et al._ (2015, 2016);
Khorasaninejad and Capasso (2017b). The replacement of bulky refractive optics
with diffractive planar elements represents an important leap forward for the
miniaturization of optical systems. However, by their nature, diffractive
optics suffer of chromatic aberrations due to the phase dispersion ($\theta$)
accumulated by light during its propagation within the optical element. It has
been largely demonstrated that this limitation can be overcome when a suitable
wavelength-dependent phase shift is directly imparted to the metasurface
distribution. However, this option requires considering focal length ($f$) and
beam deflection under a new point of view. In fact, period ($\Lambda$), number
of the discrete elements ($N$) involved in the final structure and radius of
the area occupied by the metalenses ($R_{N}$) become critical specifications
Aieta  _et al._ (2015). Similar metalenses, realized by Capasso et al., were
initially constituted by metal sub-elements replaced later on by dielectric
ones: an example is given by an array of $SiO_{2}$ nanoscopic conical bulges
with base diameter of few hundreds of nanometers and height of around $2\mu
m$Park  _et al._ (2019). While showing noticeable optical features, these
metalenses require advanced lithography techniques for their fabrication that
makes manufacturing on a large area extremely slow and time-consuming She  _et
al._ (2018). In the last years, much effort has been devoted to improve
metalenses and metasurfaces for several applications, thus involving them as
diffractive elements Ni  _et al._ (2013); Arbabi  _et al._ (2015a, b);
Naserpour, Hashemi, and Zapata-Rodríguez (2017); Banerji  _et al._ (2019);
Shrestha  _et al._ (2018), using various focal systems to obtain multiple
focal length hence correcting the chromaticity Wang  _et al._ (2017); Afridi
_et al._ (2018), for imaging Wang  _et al._ (2018); Zhou  _et al._ (2020) or
image reconstruction through metalenses array Zhu  _et al._ (2020), or
implementation in quantum electrodynamics systems Huang  _et al._ (2019). In
this work, we propose a novel and effective approach for developing ultrathin
all dielectric flat metalenses. The design involves generative adversarial
networks (GAN) whereas, for the fabrication, a two photon direct laser writing
(TP-DLW) process, upgraded to hyper resolution performance, is employed Lio
_et al._ (2020a). The latter is achieved by exploiting a plasmonic nano-cavity
metamaterial (Metal-Insulator-Metal-Insulator, MIMI) Lio  _et al._ (2020b)
that considerably improves the resolution of the TP-DLW process. The
fabrication of 2D and 3D nanostructures with typical size of few tens of
nanometers is thus enabled. The involved nano-cavity is an epsilon near zero
($\varepsilon_{NZ}$) metamaterial designed to support two resonances, centered
at $\lambda=780nm$ and $\lambda=390$ Lio  _et al._ (2020a), more details in
Method Section. In the following, it is detailed how to engineer the design of
ultra-thin flat metalenses, with thickness as low as 30-40nm, and their
fabrication through an upgraded TP-DLW process able to pattern structures
within less than a minute time. The latter has great significance for
practical and industrial applications. Despite its ultra-thin thickness, the
proposed metalenses show exciting features like apochromatic behavior, focal
length of 1.14mm, numerical aperture of 0.87, FWHM of 0.9-1.6 $\mu m$ and
depth of focus of 100-150 $\mu m$ in the visible range. These properties pave
the way for the application of our metalenses in several technological fields
like miniaturized photocameras for mobile phones, new concept microscopies,
anticounterfeiting or telecommunication.

##  II Metalenses: numerical design and TP-DLW fabrication

The proposed metalenses reproduce the optical behavior of a classic lens by
means of low-loss dielectric metasurfaces. However, their functionalities go
beyond those of classical lenses by introducing a variety of optical modes
enabling dispersive phase compensation and chromatic aberrations suppression
Park  _et al._ (2019). In order to efficiently design the metalenses, their
geometrical parameters are retrieved through an inverse design method mediated
by deep machine learning (DML). The latter exploits the ”GLOnets” Jiang and
Fan (2019) and generative adversarial networks (GANs) to realize a target
refractive index pattern. In detail, the algorithm iteratively modifies
geometrical features of polymeric nanoridges like height $H$, width $W$,
period $\Lambda$ and maximum transverse dimension $u$ (Figure S1a) until the
metalens performs the desired optical function. The range of sizes, within
which these parameters can vary, are derived from the study detailing the
hyper resolution TP-DLW technique Lio  _et al._ (2020a). In particular, height
and width for each element respectively vary from $5\pm 2$ to $40\pm 2nm$ and
from $200\pm 2$ to $500\pm 2nm$. At each iteration, the developed GAN
algorithm modifies the geometrical parameters to let the nanoridges design
converge to desired focal length $f=(1:1.2mm)$ and $u=(0:100\mu m)$. At this
point of the design, the spatial distribution of the nanoridges is 1D. The
metalens with cylindrical symmetry is obtained by revolving the nanoridge
pattern around its propagation axis. The GAN uses the actual geometrical
parameters and the considered wavelength to calculate the produced electric
near field $E(x_{0},y)$ and compare it with the one that would be produced by
an equivalent (same focus distance) classic lens. For each considered
wavelength, namely $405nm,532nm,633nm,780nm$ and $1000nm$, the calculated
geometrical parameters that satisfy the comparison within a certain error are
stored. Details about the adapted ”GLOnets” algorithm are reported in the
supplementary Figure S1b, some examples of the developed element arrangement
are reported in Figure S2. By exploiting the following Fresnel diffraction
formula, the calculated near field $E(x_{0},y)$ can be used to obtain the
envelope $E(f,u)$ in far field of the wave propagating from the metalens.

$E(f,u)=\frac{1}{\sqrt{\lambda}}\int_{-\infty}^{\infty}E(x_{0},y)e^{-i\pi
y^{2}/(\lambda f)}e^{-i\pi uy/(\lambda f)}dy$ (1)

where $u/(\lambda f)$ is the spatial frequency, $E(x_{0},y)$ and $E(f,u)$ are
the electric field amplitude distributions at the exit plane $x=x_{0}$ and the
focal plane $x=f$, respectively. From the above analysis, the metalens
structure has the following retrieved parameters: total radius $u=100\mu m$,
period $\Lambda=1\mu m$, height and width of the nanoridge element $H=30nm$
and $W=200nm$ respectively. In order to avoid hole diffraction, the inner
center of the metalens is filled with lines. By using the previous formula
with the retrieved geometrical parameters, the beam profile at the common
focal length of 1.14 mm is evaluated for each considered wavelength; the
results are reported in Figure 1a.

Figure 1: a) Beam envelope at the focal plane (common for five different
wavelengths) calculated through the Fresnel diffraction formula by using the
near field $E(x_{0},y)$ distribution provided by the GAN algorithm. b)
Electric field evaluated by using Finite Element Method (FEM) beam envelope
method for five different wavelengths and c) Beam profile retrieved at the
common focal plane (FP) by calculating the Point Spread Function (PSF); d)
Schematic view of numerically designed achromatic metalens considering the
hyper resolution achieved through the TP-DLW technique upgraded by a metal-
insulator-metal-insulator (MIMI) nano-cavity.

Figure 1b reports a further analysis conducted calculating the electric far
field for each considered wavelength by means of a Finite Element Method (FEM)
beam envelope model. The performance of an imaging system can also be
quantified by calculating its point spread function (PSF). The amplitude PSF
of a lens is defined as the transverse spatial variation of the amplitude of
the image received at the focal plane (FP) when the lens is illuminated by a
perfect plane source Kino and Corle (1996) (Figure 1c). It is worth noting
that the agreement between the results obtained using the different methods is
quite impressive: each wavelength profile converges to the focal plane with a
Full Width Half Maximum (FWHM) that is less than the diameter of the whole
metalens, further confirming its achromatic behavior. Thereafter, the depth of
focus (DOF), or the distance from the focal point to the position where the
intensity profile drops to half its maximum, is evaluated from the electric
far field maps. This results close to $50\mu m$ for shorter wavelength and
more than $200\mu m$ for longer ones, as indicated by dashed green lines drawn
on Figure 1b. According with the performed numerical simulation and past
studies Khorasaninejad and Capasso (2017a); Decker  _et al._ (2019), the
designed metalenses present also a zero phase shift at the same focal length
for all considered wavelengths, see Figure S3. This feature covers a
fundamental role in the broadband achromatic lens projection. Once designed,
the metalenses are fabricated exploiting an upgraded two-photon direct laser
writing technique which adopts a metal-insulator-metal-insulator (MIMI) nano-
cavity for the hyper resolution performance (Figure 1d). In the last years,
the TP-DLW technique has emerged as a new frontier for micro- and nano-
fabrication often replacing electron beam and UV lithography. It is actually
used to realize miniaturized efficient bulk optics Gissibl  _et al._ (2016).
However, the present technology performance needs to be pushed to its limits
to achieve the required patterning of all dielectric metasurfaces like ultra-
thin diffractive optical elements Chen  _et al._ (2013, 2018) or complex
circular structures made by single pillars Arbabi  _et al._ (2020). A work
from our group has recently demonstrated that a substantial improvement of the
TP-DLW capabilities is obtained when the process is upgraded through the
exploitation of a MIMI substrate. The introduction of the MIMI metamaterial in
the TP-DLW process represents a ground-breaking technological leap: typical
sizes of few tens and few hundreds of nanometers in height and slit width
respectively Lio  _et al._ (2020a) are easily achievable and the fabrication
of all dielectric ultra-thin nano-structures as flat metalenses is obtained in
a fraction of the time required by other standard processes. The array of
($5X4$) metalenses considered in this work has been fabricated in less than
five minutes thus confirming the fast fabrication procedure as extremely
useful for mass production and application in end-user imaging devicesZhu  _et
al._ (2020). After writing and development steps, a morphological
characterization on the fabricated metalenses has been performed using an
Atomic Force Microscope (AFM), the acquired image are reported in Figure 2a.
By showing a width $W\sim 200nm$ and a height of only $H\sim 30nm$, the AFM
profile of each circle (Figure 2b) confirms the ultra-thin characteristic of
the proposed metalens .

Figure 2: a) AFM morphological characterization of the metalenses with b) the
related height and width profile of the circle constituting it. c) Optical
image acquired with microscope with magnification 20x of the metalenses
matrix, in the inset a zoom at 50x of one structure. Particular of one
metalens at 50x observed using different band pass filter on the incoming
light.

Finally, an optical microscopy analysis has been conducted highlighting the
metalenses matrix and details of each one. Figure 2d shows an image acquired
with a magnification 20x while the inset shows the zoom (magnification 50x) of
a single structure. Successively, according with the peculiar feature of the
nano-cavity Lio  _et al._ (2020b), a metalens has been observed including a
band pass filter (BPF) centered at three different wavelengths, two close to
the double $\varepsilon_{NZ}$ resonances and one far from them. The first test
has been performed by using a filter centered at $\lambda_{BPF}=400nm$. In
this condition, the structure is very clear and it is also possible to
distinguish horizontal lines in the inner circle of the metalens (Figure 2c,
top right). These lines fill the inner center of the metalens to avoid hole
diffraction. The same test has been repeated with a BPF centered at
$\lambda_{BPF}=460nm$, now the structure appears blur and the circles broad
(Figure 2c, bottom left). Finally, the frame reported in Figure 2c bottom
right has been collected by using a BPF centered at $\lambda_{BPF}\sim 780nm$.
As expected, the system appears again very clear and the features of each
structure are highly defined. The high resolution at blue and infrared
wavelengths is an effect related to the double resonances and a consequence of
the zero light spreading property of the $\varepsilon_{NZ}$ MIMI metamaterial
exploited during the fabrication process Lio  _et al._ (2020a). Further
detailed about the unusual optical behaviour of $\varepsilon_{NZ}$
metamaterials can be found in past studiesFang and Zhang (2002); Casse  _et
al._ (2010); Zhao  _et al._ (2011); Newman, Cortes, and Jacob (2013); Shekhar,
Atkinson, and Jacob (2014).

##  III Experimental Results:
Focal Length of all-dielectric metalenses

Figure 3: a) Schematic representation of the setup used to characterize the
metalens. Stacks of the beam divergence from the focal length ($1.14mm$)
indicated as z=0 to $\pm 200\mu m$ at b) $\lambda=405nm$ c)$\lambda=532nm$
d)$\lambda=633nm$ e) using the previous three lasers at the same time and f)
using a Xenon lamp source. g) Full Width Half Maximum (FWHM) for all
considered cases. It is $0.9\mu m$ for each single laser (BL, GL, RL), $1.5\mu
m$ for the three lasers (3L) and $1.65\mu m$ for the white lamp (WL). h) CIE
chromaticity xy-coordinate plot for each laser wavelength, combination of them
(3L) and white lamp (WL); for laser with $\lambda=445nm$ and the combination
with the GL and RL producing RGB spot (grey circle).

A complete optical characterization was performed on the fabricated metalenses
exploiting the confocal home built setup shown in Figure 3a. More details are
available in materials and methods. The use of a beam profiler (BP) permits to
collect, at the same time, the intensity of the signal, the beam profile, high
quality 2D pictures of the sample and a 3D reconstruction made by the
intensity, as largely done in past studies She  _et al._ (2018); Decker  _et
al._ (2019); Lee  _et al._ (2020). The focal length, the beam divergence and
the point spread function (PSF) analysis has been performed using three lasers
(3L), namely a blue one at $\lambda=405nm$ (BL), a green one at
$\lambda=532nm$ (GL) and red one at $\lambda=633$ (RL), and a white light
source originating from a Xenon lamp (WL) with intensity of $50\mu W/mm^{2}$.
The response of the metalens at different propagation distances using each
laser at one-by-one wavelength is reported in Figure 3b,c and d respectively.
Then, the three lasers beams impinge on the sample at the same time to deduce
the metalens behavior when enlightened by a white RGB coherent source (Figure
3e). Finally, the metalenses are studied under the effect of a broad range
white light ($250nm$ to $1100nm$) as reported in Figure 3f. All the
measurements validate the metalens broadband achromatic behavior with a focal
length of 1.14 mm. The latter has been experimentally measured using the
procedure reported in the supplementary Figure S4. In order to compare the
numerical prediction and the experimental results, the full width half maximum
(FWHM) at the focal plane ($f$) for each considered case has been evaluated
and reported in Figure 3g. For each single laser wavelength (BL, GL and RL)
the PSF presents a value of FWHM of only $0.9\mu m$, for 3L the FWHM is
$1.5\mu m$ and for the WL it is $1.6\mu m$. The experimental depth of focus
(DOF) of the metalenses, which depends on the impinging wavelength, agrees
with the numerical studies showing average values about $|100-150|\mu m$
around the common focal distance for all wavelengths, see Figure S5 for more
details. As described in the numerical section, the DOF is shorter in the
visible spectrum than in the IR one. Similar results have been demonstrated
for grating-based metasurfaces Bayati  _et al._ (2020). The collected results
lead to evaluate the numerical aperture (NA) of the proposed metalens to
$NA\simeq nD/2f=0.87$, details are reported in the supplementary information.
Such high NA value is extremely useful for realizing very precise lens systems
for imaging and confocal instruments. Finally, in order to further confirm the
achromatic behavior of the proposed metalens, a CIE chromaticity xy-coordinate
graph has been evaluated for all the considered cases and reported in Figure
3h. The use of a $\lambda=405nm$ produces a blue/purple point and the three
collinear lasers (3L) results shifted in a blur white on the tri-stimulus
region (warm light)Broadbent (2004). As expected, the white lamp is centered
on the pure white point (cold light). In order to further confirm the RGB
behavior, another test has been done using a laser with a different blue
wavelength, $\lambda=445nm$, obtaining this time a practically blue spot. Now,
if this blue laser is combined with previous green (GL) and red (RL) lasers, a
white spot comparable with the WL is obtained, reported in the gamut plot as
the RGB grey circle. These results corroborate the broadband achromatic
behavior of the proposed ultra-thin dielectric metalenses fabricated with a
hyper resolute TP-DLW lithography process. The whole characterization confirms
extremely exciting properties of the proposed metalens that can be exploited,
due to their miniaturization, as optical components for micro-nano imaging,
smart devices, in situ investigation as a local spectroscopic probe, color
sensitivity and so on Suwannasopon  _et al._ (2019); Huang  _et al._ (2019);
Fan  _et al._ (2019); Wei  _et al._ (2020); Lee  _et al._ (2020); Malekovic
_et al._ (2020); Kim  _et al._ (2020).

##  IV Conclusions and Perspectives

In this work, ultra-thin flat metalenses with vertical sizes in the order of
few tens of nanometers have been reported. Design through GANs algorithms,
fabrication by hyper-resolute TP-DLW process and optical characterization have
bbeen demonstrated. The metalenses are able to work in a broad visible range
focusing light at 1.14 mm focal lengtt, DOF of about $|100-150|\mu m$ and high
numerical aperture of 0.87. These metalenses present ground-breaking
advantages as low-cost single step and time-saving fabrication, nanometric
sizes, all dielectric composition. Moreover, they are engineered to be
achromatic across the entire visible spectrum. The achieved results represent
a significant advance in the state-of-the-art of research on metalenses, that
are typically limited in their application by bandwidth and chromatic
aberration. These achromatic ultra-thin flat metalenses can find vast
implementation across industry and scientific research, such as in
miniaturized optical devices, advanced microscopy, local imaging and nano-
lithography.

## Materials and methods

Numerical model The metalens design has been performed using numerical models
implemented in Python and in COMSOL Multiphysics. The first step is using
Python to code a generative adversarial networks (GANs) able to take in
account all element variables. The developed GAN algorithms is based on
”GLOnets” Jiang and Fan (2019). Then, by a training based on a deep machine
learning (DML) process, the neural network returns the best parameters that
satisfy the fixed conditions, in particular a common focal length at multiple
wavelengths. Successively, the parameters retrieved using GAN code have been
used to verify that the near field $E(x_{0},y)$ produced by the metasurface
satisfy the Fresnel diffraction formula ($E(f,u)$) and that the beam profile
at the fixed focal length is maximum. Finally, in order to further validate
the results in far-field numerical simulations, the beam envelop module of
COMSOL Multiphysics has been used returning the electric far-field along the
propagation direction. By placing a cut line at the focal length in the
electric field map, the point spread function (PSF) of the propagated beam is
retrieved.
Metalenses fabrication The nano-cavity has been deposited by DC sputtering on
glass substrate with thickness of $\sim 100\mu m$) and it is constituted
starting from the glass as follow: $30nm$ Ag layer, $160nm$ ZnO cavity, $30nm$
Ag layer and $30nm$ ZnO layer. The system is placed in front of the two-photon
source. This lithography apparatus is used to define the metasurface lens
pattern in a drop of photo-sensitive resist placed on top of the nano-cavity.
The laser power used for the process is equal to $20mW$ and the scan speed is
$4000\mu m$. After the writing process, the whole system (nano-cavity with
cured and uncured photo-resin) is soaked in a bath of propylene glycol methyl
ether acetate (PGMEA) for 15 minutes and then cleaned in a bath of isopropanol
alcohol (IPA) for other 5 minutes. Lio  _et al._ (2020a)
AFM characterization AFM has been performed using a Zeiss LSM 780 with high
resolute tips with precision of $\pm 2nm$. Each scan has been collected at
high resolution 1024x1024 px in order to reduce any background noise.
Optical setup and characterization A home-built setup composed of various
optical elements has been used. A collimation line for white light (WLS) that
is coming from a Xenon lamp is realized with a 40X (Obj. 40X) objective and a
fiber coupler (FC) mounted on a 3-axis stage (XYZmM) in order to collect the
WLS by a fiber and send it collimated on the main confocal line. Then, on the
same path it has been placed another 3-axis stage with an objective of 10X
(Obj. 10X) used only to detect the metalenses on the sample. The objective is
moved backward during the measurements in order to maintain the collimation of
light impinging on the metalens. Beyond the metalens, there is an other
objective 40X (Obj. 40X) that collects the signal from the sample and
collimates it into the detector represented by a Thorlabs beam profiler (BP),
Thorlabs BC106N-VIS spectral range from $350nm$ to $1100nm$, mounted on a
double-axis stage (ZXmM). Another lens (L) collects the signal from the
objective and sent it on the BP ccd. For acquire the spectrum needed to
calculate the CIE chromaticity xy-coordinate, a flip mirror (FM) is placed
before the collection lens aming to send the collimated signal on an iris that
select only the focal spot. Then, the spectrum is acquired using a Flame Ocean
Optics spectrometer. The metalens as moved forward and backward along Z in
order to collect the beam divergence. Between the fiber collimator and the
first objective it has been placed a beam splitter (BS2) that sends the light
from lamp or from the laser (line) in the main path. The laser line is
composed by a blue laser that sends the light through a dichroic mirror (DM1)
that reflects the green and leads both lasers on a beam splitter (BS1) that
collects the red laser and it sent all lasers on the second beam splitter
(BS2).

## Contributions

G.E.L. conceived the idea, designed and numerically simulated the metalenses.
G.E.L. improved the GLOnets code using a suitable GAN algorithms to design the
metalenses. G.E.L., A.F. and T.R. fabricated and fully characterized the
structures. G.E.L., A.F. and D.M.A. performed the optical measurements and
analyzed the data. R.C, A.DL. and M.G. devised the experiments and supervised
the work. G.E.L. and A.F. prepared and wrote the manuscript with input from
all authors.

## Acknowledgements

The authors thank the Infrastructure “BeyondNano” (PONa3-00362) of CNR-Nanotec
for the access to research instruments. They also thank the “Area della
Ricerca di Roma 2”, Tor Vergata, for the access to the ICT Services (ARToV-
CNR) for the use of the COMSOL Multiphysics Platform, Origin Lab, and Matlab.

## References

  * Genevet  _et al._ (2017) P. Genevet, F. Capasso, F. Aieta, M. Khorasaninejad, and R. Devlin, “Recent advances in planar optics: from plasmonic to dielectric metasurfaces,” Optica 4, 139–152 (2017).
  * Chen  _et al._ (2018) W. T. Chen, A. Y. Zhu, V. Sanjeev, M. Khorasaninejad, Z. Shi, E. Lee, and F. Capasso, “A broadband achromatic metalens for focusing and imaging in the visible,” Nature nanotechnology 13, 220–226 (2018).
  * Chen  _et al._ (2019) W. T. Chen, A. Y. Zhu, J. Sisler, Z. Bharwani, and F. Capasso, “A broadband achromatic polarization-insensitive metalens consisting of anisotropic nanostructures,” Nature communications 10, 1–7 (2019).
  * Fan  _et al._ (2019) Z.-B. Fan, H.-Y. Qiu, H.-L. Zhang, X.-N. Pang, L.-D. Zhou, L. Liu, H. Ren, Q.-H. Wang, and J.-W. Dong, “A broadband achromatic metalens array for integral imaging in the visible,” Light: Science & Applications 8, 1–10 (2019).
  * Idjadi and Aflatouni (2020) M. H. Idjadi and F. Aflatouni, “Nanophotonic phase noise filter in silicon,” Nature Photonics 14, 234–239 (2020).
  * Gao  _et al._ (2013) Y. Gao, Z. Xin, B. Zeng, Q. Gan, X. Cheng, and F. J. Bartoli, “Plasmonic interferometric sensor arrays for high-performance label-free biomolecular detection,” Lab on a Chip 13, 4755–4764 (2013).
  * Yu and Capasso (2014) N. Yu and F. Capasso, “Flat optics with designer metasurfaces,” Nature materials 13, 139–150 (2014).
  * Byrnes  _et al._ (2016) S. J. Byrnes, A. Lenef, F. Aieta, and F. Capasso, “Designing large, high-efficiency, high-numerical-aperture, transmissive meta-lenses for visible light,” Optics express 24, 5110–5124 (2016).
  * Khorasaninejad  _et al._ (2017) M. Khorasaninejad, Z. Shi, A. Y. Zhu, W.-T. Chen, V. Sanjeev, A. Zaidi, and F. Capasso, “Achromatic metalens over 60 nm bandwidth in the visible and metalens with reverse chromatic dispersion,” Nano letters 17, 1819–1824 (2017).
  * Khorasaninejad and Capasso (2017a) M. Khorasaninejad and F. Capasso, “Metalenses: Versatile multifunctional photonic components,” Science 358, eaam8100 (2017a).
  * Cui, Bai, and Sun (2019) T. Cui, B. Bai, and H.-B. Sun, “Tunable metasurfaces based on active materials,” Advanced Functional Materials 29, 1806692 (2019).
  * Chen, Zhu, and Capasso (2020) W. T. Chen, A. Y. Zhu, and F. Capasso, “Flat optics with dispersion-engineered metasurfaces,” Nature Reviews Materials , 1–17 (2020).
  * Khorasaninejad  _et al._ (2015) M. Khorasaninejad, F. Aieta, P. Kanhaiya, M. A. Kats, P. Genevet, D. Rousso, and F. Capasso, “Achromatic metasurface lens at telecommunication wavelengths,” Nano letters 15, 5358–5362 (2015).
  * Khorasaninejad  _et al._ (2016) M. Khorasaninejad, W. T. Chen, R. C. Devlin, J. Oh, A. Y. Zhu, and F. Capasso, “Metalenses at visible wavelengths: Diffraction-limited focusing and subwavelength resolution imaging,” Science 352, 1190–1194 (2016).
  * Khorasaninejad and Capasso (2017b) M. Khorasaninejad and F. Capasso, “Metalenses: Versatile multifunctional photonic components,” Science 358, eaam8100 (2017b).
  * Aieta  _et al._ (2015) F. Aieta, M. A. Kats, P. Genevet, and F. Capasso, “Multiwavelength achromatic metasurfaces by dispersive phase compensation,” Science 347, 1342–1345 (2015).
  * Park  _et al._ (2019) J.-S. Park, S. Zhang, A. She, W. T. Chen, P. Lin, K. M. Yousef, J.-X. Cheng, and F. Capasso, “All-glass, large metalens at visible wavelength using deep-ultraviolet projection lithography,” Nano letters 19, 8673–8682 (2019).
  * She  _et al._ (2018) A. She, S. Zhang, S. Shian, D. R. Clarke, and F. Capasso, “Large area metalenses: design, characterization, and mass manufacturing,” Optics express 26, 1573–1585 (2018).
  * Ni  _et al._ (2013) X. Ni, S. Ishii, A. V. Kildishev, and V. M. Shalaev, “Ultra-thin, planar, babinet-inverted plasmonic metalenses,” Light: Science & Applications 2, e72–e72 (2013).
  * Arbabi  _et al._ (2015a) A. Arbabi, Y. Horie, M. Bagheri, and A. Faraon, “Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission,” Nature nanotechnology 10, 937 (2015a).
  * Arbabi  _et al._ (2015b) A. Arbabi, Y. Horie, A. J. Ball, M. Bagheri, and A. Faraon, “Subwavelength-thick lenses with high numerical apertures and large efficiency based on high-contrast transmitarrays,” Nature communications 6, 1–6 (2015b).
  * Naserpour, Hashemi, and Zapata-Rodríguez (2017) M. Naserpour, M. Hashemi, and C. J. Zapata-Rodríguez, “Recent progress in far-field optical metalenses,” Metamaterials-Devices and Applications (2017).
  * Banerji  _et al._ (2019) S. Banerji, M. Meem, A. Majumder, F. G. Vasquez, B. Sensale-Rodriguez, and R. Menon, “Imaging with flat optics: metalenses or diffractive lenses?” Optica 6, 805–810 (2019).
  * Shrestha  _et al._ (2018) S. Shrestha, A. C. Overvig, M. Lu, A. Stein, and N. Yu, “Broadband achromatic dielectric metalenses,” Light: Science & Applications 7, 1–11 (2018).
  * Wang  _et al._ (2017) S. Wang, P. C. Wu, V.-C. Su, Y.-C. Lai, C. H. Chu, J.-W. Chen, S.-H. Lu, J. Chen, B. Xu, C.-H. Kuan, _et al._ , “Broadband achromatic optical metasurface devices,” Nature communications 8, 1–9 (2017).
  * Afridi  _et al._ (2018) A. Afridi, J. Canet-Ferrer, L. Philippet, J. Osmond, P. Berto, and R. Quidant, “Electrically driven varifocal silicon metalens,” Acs Photonics 5, 4497–4503 (2018).
  * Wang  _et al._ (2018) S. Wang, P. C. Wu, V.-C. Su, Y.-C. Lai, M.-K. Chen, H. Y. Kuo, B. H. Chen, Y. H. Chen, T.-T. Huang, J.-H. Wang, _et al._ , “A broadband achromatic metalens in the visible,” Nature nanotechnology 13, 227–232 (2018).
  * Zhou  _et al._ (2020) Y. Zhou, H. Zheng, I. I. Kravchenko, and J. Valentine, “Flat optics for image differentiation,” Nature Photonics , 1–8 (2020).
  * Zhu  _et al._ (2020) X. Zhu, W. Fang, J. Lei, Z. Li, F. Xie, Y. Cao, Y. Zhang, F. Qin, and X. Li, “Supercritical lens array in a centimeter scale patterned with maskless uv lithography,” Optics Letters 45, 1798–1801 (2020).
  * Huang  _et al._ (2019) T.-Y. Huang, R. R. Grote, S. A. Mann, D. A. Hopper, A. L. Exarhos, G. G. Lopez, G. R. Kaighn, E. C. Garnett, and L. C. Bassett, “A monolithic immersion metalens for imaging solid-state quantum emitters,” Nature communications 10, 1–8 (2019).
  * Lio  _et al._ (2020a) G. E. Lio, T. Ritacco, A. Ferraro, A. D. Luca, R. Caputo, and M. Giocondo, “Hyper resolution two photon direct laser writing using enz nano-cavity,” (2020a), arXiv:2007.13509 [physics.app-ph] .
  * Lio  _et al._ (2020b) G. E. Lio, A. Ferraro, M. Giocondo, R. Caputo, and A. De Luca, “Color gamut behavior in epsilon near-zero nanocavities during propagation of gap surface plasmons,” Advanced Optical Materials , 2000487 (2020b).
  * Jiang and Fan (2019) J. Jiang and J. A. Fan, “Global optimization of dielectric metasurfaces using a physics-driven neural network,” Nano letters 19, 5366–5372 (2019).
  * Kino and Corle (1996) G. S. Kino and T. R. Corle,  _Confocal scanning optical microscopy and related imaging systems_ (Academic Press, 1996).
  * Decker  _et al._ (2019) M. Decker, W. T. Chen, T. Nobis, A. Y. Zhu, M. Khorasaninejad, Z. Bharwani, F. Capasso, and J. Petschulat, “Imaging performance of polarization-insensitive metalenses,” ACS Photonics 6, 1493–1499 (2019).
  * Gissibl  _et al._ (2016) T. Gissibl, S. Thiele, A. Herkommer, and H. Giessen, “Two-photon direct laser writing of ultracompact multi-lens objectives,” Nature Photonics 10, 554 (2016).
  * Chen  _et al._ (2013) X. Chen, L. Huang, H. Mühlenbernd, G. Li, B. Bai, Q. Tan, G. Jin, C.-W. Qiu, T. Zentgraf, and S. Zhang, “Reversible three-dimensional focusing of visible light with ultrathin plasmonic flat lens,” Advanced Optical Materials 1, 517–521 (2013).
  * Arbabi  _et al._ (2020) A. Arbabi, E. Arbabi, M. Mansouree, S. Han, S. M. Kamali, Y. Horie, and A. Faraon, “Increasing efficiency of high numerical aperture metasurfaces using the grating averaging technique,” Scientific Reports 10, 1–10 (2020).
  * Fang and Zhang (2002) N. Fang and X. Zhang, “Imaging properties of a metamaterial superlens,” in  _Proceedings of the 2nd IEEE Conference on Nanotechnology_ (IEEE, 2002) pp. 225–228.
  * Casse  _et al._ (2010) B. Casse, W. Lu, Y. Huang, E. Gultepe, L. Menon, and S. Sridhar, “Super-resolution imaging using a three-dimensional metamaterials nanolens,” Applied Physics Letters 96, 023114 (2010).
  * Zhao  _et al._ (2011) Y. Zhao, A. A. Nawaz, S.-C. S. Lin, Q. Hao, B. Kiraly, and T. J. Huang, “Nanoscale super-resolution imaging via a metal–dielectric metamaterial lens system,” Journal of Physics D: Applied Physics 44, 415101 (2011).
  * Newman, Cortes, and Jacob (2013) W. D. Newman, C. L. Cortes, and Z. Jacob, “Enhanced and directional single-photon emission in hyperbolic metamaterials,” JOSA B 30, 766–775 (2013).
  * Shekhar, Atkinson, and Jacob (2014) P. Shekhar, J. Atkinson, and Z. Jacob, “Hyperbolic metamaterials: fundamentals and applications,” Nano convergence 1, 14 (2014).
  * Lee  _et al._ (2020) D. Lee, J. Gwak, T. Badloe, S. Palomba, and J. Rho, “Metasurfaces-based imaging and applications: from miniaturized optical components to functional imaging platforms,” Nanoscale Advances (2020).
  * Bayati  _et al._ (2020) E. Bayati, R. Pestourie, S. Colburn, Z. Lin, S. G. Johnson, and A. Majumdar, “Inverse designed metalenses with extended depth of focus,” ACS Photonics 7, 873–878 (2020).
  * Broadbent (2004) A. D. Broadbent, “A critical review of the development of the cie1931 rgb color-matching functions,” Color Research & Application: Endorsed by Inter-Society Color Council, The Colour Group (Great Britain), Canadian Society for Color, Color Science Association of Japan, Dutch Society for the Study of Color, The Swedish Colour Centre Foundation, Colour Society of Australia, Centre Français de la Couleur 29, 267–272 (2004).
  * Suwannasopon  _et al._ (2019) S. Suwannasopon, F. Meyer, C. Schlickriede, P. Chaisakul, J. Limtrakul, T. Zentgraf, N. Chattham, _et al._ , “Miniaturized metalens based optical tweezers on liquid crystal droplets for lab-on-a-chip optical motors,” Crystals 9, 515 (2019).
  * Wei  _et al._ (2020) Y. Wei, Y. Wang, X. Feng, S. Xiao, Z. Wang, T. Hu, M. Hu, J. Song, M. Wegener, M. Zhao, _et al._ , “Compact optical polarization-insensitive zoom metalens doublet,” Advanced Optical Materials , 2000142 (2020).
  * Malekovic  _et al._ (2020) M. Malekovic, M. Urann, U. Steiner, B. D. Wilts, and M. Kolle, “Soft photonic fibers for colorimetric solvent vapor sensing,” Advanced Optical Materials , 2000165 (2020).
  * Kim  _et al._ (2020) J. A. Kim, D. J. Wales, A. J. Thompson, and G.-Z. Yang, “Fiber-optic sers probes fabricated using two-photon polymerization for rapid detection of bacteria,” Advanced Optical Materials 8, 1901934 (2020).

## Supplementary material

## Generative Adversarial Networks (GANs)

For the metalens design we used a conditional Generative Adversarial Networks
(GANs), based on GLOnets (Metanet code) Jiang and Fan (2019), that have
wavelength, geometrical parameters, focal length and, if required, deflection
angle as inputs. We designed and optimized an ensemble of polymer nano ridges
that operate across a wide range of wavelengths in the UV-NIR range. This code
builds on analysis of conditional GAN that can optimize only a single device
in a training session. This means that focal length and wavelength are fixed,
while the other parameters are continuously changed from the input array
values. The thickness of the nano ridges is varied in the range from
$H\sim(5\pm 2:30\pm 2)nm$ and the incident light is TM/TE-polarized. For each
device, the nanoridge pattern period ($\Lambda$) is divided into segments
$N=\Lambda/W$ where $W\sim(200\pm 2:500\pm 2)nm$ represents the nanoridge
width. For each segment, the refractive index of polymer or air is considered.
The nanoridge pattern deflects normally incident light to a direction meeting
the focal point (Figure S1 a). The optimization scope is to maximize the
deflection efficiency of the metalens given an operating wavelength ranging
from 400 to 1000 nm and an outgoing angle ranging from $40^{\circ}$ to
$80^{\circ}$. A schematic representation of our conditional GAN flow chart is
presented in Figure S1b. The input are the operating wavelength $\lambda$, the
geometrical parameters $H$, $W$ and $\Lambda$, the desired focal length ($f$),
the maximum transversal dimension ($u$) and the near field equation. The
output is the refractive index profile n(u) of the device. The weights of the
neurons are parameterized as $W_{n}$. The generator, conditioned on
($\lambda$, $H$, $W$ and $\Lambda$, $f$ , $u$) produces results which are
compared with the near field profile of a classical lens. Due to the circular
symmetry of a classical lens, the nanoridge pattern is then revolved around
the propagation direction producing a circular pattern with the features
designed by the mediated deep machine learning algorithm. This comparison is
used to discard or store the obtained parameters. Finally, the parameters that
allow a common focal length for all the considered wavelengths are used for
the the numerical simulation using Finite Element Method models and the
experimental fabrication/implementation.

Figure S1: Global optimization based on a generative neural network (GAN)
opportunely modified. a) Schematic view of nano ridges that constitutes the
metalens. b) Schematic representation of the conditional GANs code for
metalens generation. The generator is built on fully connected layers (FC),
deconvolution layers (dconv), and a Gaussian filter. An identity shortcut
connection is also used by the electric near field equation. The input of GANs
code are the device’s geometrical parameters as height ($H$), width ($W$),
period ($\Lambda$) and finally wavelength ($\lambda$). During each iteration
of training, a batch of devices is generated and for each device efficiency
gradients (g) are calculated using forward electromagnetic simulations. These
gradients are propagated through the network to update the weights of neurons
($W_{n}$). The GAN verifies that the near field E(x0; y) produced by the
metasurface satisfy the Fresnel diffraction formula (E(f; u)). Then, the GAN
algorithm compares the electrical near field of the actual geometrical
parameters and the considered wavelength with the one of a classical lens
having the same focal length. Only the parameters that respect the fixed
condition over the focal length and the maximum lens transversal dimension at
multiple wavelengths are stored.

Details of the GAN algorithms are reported in Figure S2a depicting the working
principle where each constitutive element is changed until the right values
are reached. The algorithm has to satisfy an high value of efficiency for each
considered wavelength for different deflection angles $\theta$. Figure S2b
reports the more efficient evaluated element disposition patterns for some
cases $\lambda=405nm$, $\lambda=532nm$, $\lambda=633nm$, $\lambda=1000nm$. For
a bunch of iterations, the code returns an histogram with the efficiency
reached for each disposition, as reported in Figure S2c. Finally all the best
data are stored and matched to obtain the element sizes that satisfy the
initial condition as already discussed above. The enormous advantage of design
algorithms coupled with the TP-DLW cad software is represented by the direct
conversion of the produced design figure into TP-DLW writing file for a quick
and immediate metalens fabrication.

Figure S2: a) Sketch of the GAN algorithm working flow with the variable
parameters. b) Examples of produced element disposition in function of
wavelength and deflection angle, varying the height, width and period of the
constitutive metalens element. c) Histogram efficiency used to chose the
element disposition.

## Phase shift at different wavelengths

The phase shift of the proposed metalenses has been numerically analyzed, see
Figure S3. It is close to 0 rad for all wavelengths at the focal length of
$f=1.14mm$; similar results have been already presented for different
achromatic metasurfacesKhorasaninejad and Capasso (2017a); Decker  _et al._
(2019). The phase is calculated as the argument of the electric field
component along the propagation axis. The maximum transverse dimension of the
nanoridges pattern ”u”, equal to $100\mu m$, is the center of the metalens and
the phase has been evaluated using a cut line at the calculate focal distance
common for all wavelengths. The phase shift, according with past studies,
results equal to $\phi(x)\sim 2\pi n_{eff}(x)t/\lambda$. Where x represents
the dimension along the lens diameter with a maximum extension equal to
$x=2*u$. When this equation is imposed equal to zero for each wavelength, the
broadband achromatic lens design process starts Aieta  _et al._ (2015);
Genevet  _et al._ (2017); Chen, Zhu, and Capasso (2020).

Figure S3: The phase shift analyzed for different wavelengths at the same
focal length and along the lens diameter. The analyzed wavelengths are: a)
$450nm$, b) $532nm$, c) $633nm$, d) $780nm$ and finally e) $1000nm$.

In order to clarify the connection between the already studied metalens
systems and our GAN implementation, it is necessary to consider the phase
equation written as a function of the radial coordinate $r$ resulting equal to
$\phi(x)=-2\pi/\lambda_{d}(\sqrt{(}r^{2}-f^{2})-f)$, with
$r=\sqrt{(}x^{2}+y^{2})$. Fixing $f$ and the maximum lens radius
$r=u=(0:100)\mu m$, it is possible to run the GAN code until all conditions
are respected and obtaining the element sizes ($H$,$W$, $\Lambda$) as already
explained above. A comprehensive overview about the phase calculation and
metalens design with different numerical tools are reported in Chen, Zhu, and
Capasso (2020).

## Metalens focal length

In this section, it is reported the procedure to measure the metalens focal
length. A confocal homebuilt setup has been used allowing to remove the
optical elements that are not necessary for the complete measurement. The
first step consists to place the sample containing the metalens in the optical
path of the confocal configuration (Figure S4a). Then, once the focus is found
on the metalens, the input objective (10X) is removed. Finally, the light
impinges on the metalens and the output objective (32X) is moved along z until
the new focal spot is established, Figure S4b. The distance between the first
configuration and the second one corresponds to the focal length, which in our
case is $f=1.14mm$.

Figure S4: Sketch of the setup used to measure the focal length. a) Initial
confocal configuration. b) The metalens focal length measurements moving out
the input objective and moving along z the output objective until the produced
spot is again the smallest one.

## Numerical Aperture (NA)

The numerical aperture is evaluated using the following equation:

$NA=nsin\theta$ (2)

In the condition of a low lens magnification the equation (1) is reduced to

${\text{NA}}_{\text{i}}=n\sin\theta=n\sin\left[\arctan\left({\frac{D}{2f}}\right)\right]\approx
n{\frac{D}{2f}}$ (3)

Where $D$ is the diameter of the lens $D=2*R_{N}$, $n$ is the refractive index
of the metalenses $n\sim 1.56$ and $f$ is the focal length.

## Depth of Focus (DOF)

The experimental DOF of the proposed metalens has been reported for each
single light source and the combination of them, see Figure S5. It is
collected in the propagation range along $z$ from $-200\mu m$ to $200\mu m$.
As predicted by the numerical simulations, the DOF is short in the blu-green
spectral range, while it starts to increase when the wavelength is close to
the red-IR spectral range. The combination of lasers (3L) presents a DOF
extended as the red one. Instead, the withe light (WL) produces a depth of
focus very short due to the incoherent nature of the propagated light and also
because it is a continuous spectrum and the metalenses probably do not focus
all wavelengths emitted by the Xenon lamp.

Figure S5:  Full Width Half Maximum (FWHM) of the produced spot is analyzed
along the propagation direction $z$ in the range from $-200\mu m$ to $200\mu
m$. It evidences the short DOF for the blu laser (blue squares), increasing a
little bit for the green (green circles) and more for the red lasers (red up-
triangles). The three lasers (3L) at the same time presents the more extended
DOF that the previous cases (black down-triangles). The withe lamp (WL) due to
the non coherence along the propagation presents a DOF reduced that the
aspected one (dark green diamonds).