# Shape optimization for high efficiency metasurfaces: theory and
implementation

P. Dainese  Corning Research and Development Corporation, 184 Science Center
Dr, Painted Post, NY 14870, USA<EMAIL_ADDRESS>L. Marra  Corning
Research and Development Corporation, 184 Science Center Dr, Painted Post, NY
14870, USA  D. Cassara  Harvard John A. Paulson School of Engineering and
Applied Sciences 1, Cambridge, MA 02138, USA  A. Portes  Harvard John A.
Paulson School of Engineering and Applied Sciences 1, Cambridge, MA 02138, USA
J. Oh  Harvard John A. Paulson School of Engineering and Applied Sciences 1,
Cambridge, MA 02138, USA  J. Yang  Corning Research and Development
Corporation, 184 Science Center Dr, Painted Post, NY 14870, USA  A. Palmieri
Harvard John A. Paulson School of Engineering and Applied Sciences 1,
Cambridge, MA 02138, USA  J.R. Rodrigues  Corning Research and Development
Corporation, 184 Science Center Dr, Painted Post, NY 14870, USA  A.H. Dorrah
Harvard John A. Paulson School of Engineering and Applied Sciences 1,
Cambridge, MA 02138, USA  F. Capasso  Harvard John A. Paulson School of
Engineering and Applied Sciences 1, Cambridge, MA 02138, USA

###### Abstract

Complex non-local behavior makes designing high efficiency and multifunctional
metasurfaces a significant challenge. While using libraries of meta-atoms
provide a simple and fast implementation methodology, pillar to pillar
interaction often imposes performance limitations. On the other extreme,
inverse design based on topology optimization leverages non-local coupling to
achieve high efficiency, but leads to complex and difficult to fabricate
structures. In this paper, we demonstrate numerically and experimentally a
shape optimization method that enables high efficiency metasurfaces while
providing direct control of the structure complexity. The proposed method
provides a path towards manufacturability of inverse-designed high efficiency
metasurfaces.

##  1 Introduction

Wavefront shaping using metasurfaces has attracted significant scientific and
technological interest in recent years [1, 2, 3, 4, 5, 6]. The ability to
control many degrees of freedom of an incoming beam, including phase,
amplitude, polarization, and dispersion has led to a number of demonstrations
in areas ranging from imaging, polarization optics, communications, atom
trapping, optical computing and image processing, nonlinear optics and others
[7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. This versatility stems from a large
design space enabled by engineering meta-atoms with different geometrical
shapes and materials readily available in nano-fabrication, which allows
tapping into different mechanisms such as Mie scattering, waveguide
propagation phase, Pancharatnam-Berry phase, guided mode resonances, and more
broadly Bloch-mode engineering [17, 18]. Despite successful demonstrations,
designing metasurfaces is still a significant challenge due to the complex
nature of such physical mechanisms, particularly when sub-wavelength and often
non-periodic arrangement leads to strong non-local coupling between meta-atoms
[5, 19, 6].

The most common design approach is based on libraries of meta-atoms (Figure
1a), using either parametrized shapes or free-form meta-atoms [20, 21, 22, 23,
3, 7, 24]. The simplicity of this method relies on the assumption that each
meta-atom is placed on a unit cell with periodic boundary conditions, allowing
fast computation of its response using numerical methods. Metagratings,
metalenses and more complex devices have been created using libraries.
However, except for truly periodic arrays, this assumption is necessarily
broken in real devices and pillar-to-pillar interaction imposes a fundamental
limitation to performance by perturbing the realized phase and intensity
profile, making it difficult to achieve high efficiency in general. Another
limitation of the library method is that the unit cell response is computed
for a specific angle of incidence, typically normal, and for a specific
polarization. In general, however, an incoming beam may contain a spectrum of
incidence angles, as is the case even for a simple metalens illuminated with
non-collimated light, and in more general holograms used for vector beam
generation, spatial mode multiplexing and many others.

Figure 1:  Metasurface design methodologies. In the library method (a), the
phase and amplitude response of pre-determined meta-atoms is calculated as a
function of its parameters (e.g. diameter for a circular pillar). The response
is then used to create devices by dividing the surface in unit cells, each
containing a specific element of the library. In topology optimization (b),
the design domain is typically initialized with a randomized refractive index
distribution, and the adjoint algorithm is used to converge towards a final
binarized geometry. In shape optimization (c), only smooth variations in
existing geometrical features are allowed.

To overcome these challenges, topology optimization based on the adjoint
method has been proposed [25, 26, 18, 27, 28]. In this approach, an initial
random refractive index distribution iteratively converges to a final binary,
i.e. manufacturable, structure (Figure 1b). Rigorous electromagnetic
simulation of the full device ensures that the interaction between meta-
structures is fully considered. Furthermore, the adjoint formulation enables
computing gradients of a given figure of merit with respect to the refractive
index at every pixel in the domain using only two electromagnetic simulations
(forward and adjoint simulations), critical for optimization problems with
such high dimensionality. Despite successful designs with very high
efficiencies, controlling the complexity of the final structure is a
significant challenge, particularly for large scale manufacturing. During the
optimization process, the device structure evolves according to the
topological derivative at each iteration, and often leads to features that are
difficult to control in patterning and etching processes, such as the
appearance of sharp boundaries and small features such as islands of
materials, small holes and small gaps between structures. These challenges
have stimulated recent research to incorporate manufacturing robustness in the
design process [25, 29, 28, 6], either external to the adjoint formulation
such as applying structural blurring every so often during optimization, or
including penalty terms in the figure of merit to guide the topological
derivative in a way to minimize these issues.

The adjoint formulation can also be used to extract gradients of the figure of
merit with respect to boundary shifts at the interface between two materials.
In this case, the topology of the structure is unchanged during the
optimization, and only deformations of the existing boundaries occur (Figure
1c). This approach has been applied to optimize the shape of planar photonic
devices, such as waveguide splitters, crossing and bends in photonic crystal
waveguides [30, 31, 32]. Leveraging boundary gradients for metasurface design
has had limited investigation, with only a few examples where the basic meta-
atom shape remained unchanged but their sizes were optimized. This was applied
to design metalenses by optmizing the side lengths of rectangular pillars
[33], the radii of circular pillars [33], as well as to design metagratings
where the semi-axis of elliptical pillars were optimized [34, 35, 36]. In this
paper, we generalize this approach and investigate a shape optimization method
that achieves efficiencies higher than the library method by fully considering
pillar-to-pillar interaction, while providing greater control over the
structure complexity compared to topology optimization. The initial shape of
each meta-atom in the device is smoothly deformed throughout the process, with
the shape complexity controlled through a Fourier decomposition of the adjoint
boundary gradients. Direct control of all boundaries naturally incorporates
fabrication constraints such as minimum feature size and minimum gap, and by
construction, excludes the appearance of holes. Similarly to topology
optimization, shape optimization can be applied to any kind of metasurface
devices, can handle any input and target field distributions, as well as
include multiple objectives. The paper is organized as follows: we outline the
formulation in Section 2, and then apply the shape optimization method to
design several high efficiency metagratings and metalenses in Section 3.
Experimental results are presented in Section 4, and we draw conclusions in
Section 5.

##  2 Shape optimization formulation

Topology and shape optimization in photonics are both inverse design
techniques based on the adjoint method [25, 26]. A flowchart of the shape
optimization method is shown in Figure 2. The design domain is initialized
with a given set of meta-atoms, for example using a uniform array of circular
pillars, a library based device, or even a random distribution of pillars. At
each iteration, two electromagnetic simulations are performed, a forward and
an adjoint simulation, from which shape gradients for all pillars are computed
for a given figure of merit. If at any iteration the figure of merit of the
device reaches a desired target or if it converges to a local maxima, the
optimization stops. Otherwise, the gradients are used to update the shapes and
continue on to the next iteration. Before updating the shapes, the gradients
are first decomposed in Fourier basis and further fabrication constraints can
be implemented.

Figure 2: Flow chart of the optimization process. The design domain is
initialized with a given set of meta-atoms, for example with a library-based
metasurface. Forward and adjoint simulations are performed and the figure of
merit computed. If the figure of merit satisfies the target or if it has
converged the optimization stops. If not, the shifting boundary gradients are
calculated from the simulated fields, followed by Fourier decomposition and
possibly additional fabrication constraints. The shape is then deformed
according to the gradients and a new simulation is performed. The process
repeats until convergence.

The mathematical formulation of the adjoint shape gradient and subsequent
Fourier decomposition is outlined here, and more details can be found in
Appendix A. As illustrated in Figure 3a, in the forward simulation the
incident field propagates through the metasurface, while in the adjoint
simulation the target field propagates backwards. In this example, the forward
field is simply a plane wave incident at normal direction, and the target
field is a plane wave propagating at a certain deflection angle. Through a
reciprocity argument [37], the forward $(E)$ and adjoint $(E_{a})$ fields at
the surface of each pillar can be used to compute variations in the figure of
merit $(\eta)$ due to an arbitrary (but small) boundary shift $u_{\perp}$:

$\delta\eta=\frac{\omega\delta\epsilon}{2P}Re\left[jF^{*}\int{u_{\perp}(E_{\parallel}\cdot
E_{a,\parallel}+\frac{1}{\epsilon_{1}\epsilon_{2}}D_{\perp}\cdot D_{a,\perp})\
da}\right],$ (1)

where $\omega$ is the optical frequency,
$\delta\epsilon=\epsilon_{2}-\epsilon_{1}$ is the difference in dielectric
permittivity between the meta-atom and the surrounding medium, and $P$ is a
normalization power. The subscripts $\parallel$ and $\perp$ represent the
tangential and normal components of the fields ($F$ is related to the device
efficiency, $\eta=|F|^{2}$). The integration in equation 1 is performed on
every pillar’s surface. Since we wish to retain vertical pillars for top-down
fabrication, we enforce the boundary displacement $u_{\perp}$ to be uniform
along $z$, and allow only variation along the cross-sectional boundary. With
that, we can explicitly write the efficiency change as an integral along the
pillar closed cross-sectional boundary as:

$\delta\eta=\frac{\omega\delta\epsilon}{2P}\oint{u_{\perp}g\ ds},$ (2)

where the gradient function $g$ is defined accordingly as:

$g=Re\left[jF^{*}\int{(E_{\parallel}\cdot
E_{a,\parallel}+\frac{1}{\epsilon_{1}\epsilon_{2}}D_{\perp}\cdot D_{a,\perp})\
dz}\right].$ (3)

Figure 3: Details of the shape optimization algorithm. The forward and adjoint
fields are computed in the design domain, and are then used to calculate a
gradient function $g$ on the surface of all existing meta-atoms. (a) and (b)
show the fields and the derived boundary gradients for TM and TE
polarizations. The gradient function determines the magnitude and direction
that each point on the surface must be displaced to increase the device
efficiency. In (c), the resulting shape gradient is plotted (blue) along with
its Fourier decomposition (black) for the pillar highlighted in (b). One can
observe that the gradient containts significant contributions up to order
$m=6$.

From the expression of $\delta\eta$ in equation 2, it is clear that the
efficiency change is always positive if we choose the boundary deformation
$u_{\perp}$ to follow the gradient function $g$, i.e., by choosing
$u_{\perp}=hsign(\delta\epsilon)g$, where $h$ is a scaling factor.
Mathematically, this is simply the projection of the gradient function on
itself, as the integral represents the inner product in the closed boundary
domain. The term $(E_{\parallel}\cdot
E_{a,\parallel}+\frac{1}{\epsilon_{1}\epsilon_{2}}D_{\perp}\cdot D_{a,\perp})$
in equation 3 explicitly determines how the forward and adjoint fields
determine the shape gradient $u_{\perp}$. This term is plotted on the pillars’
surface in Figure 2b, with an arbitrary -1 to 1 color scale. There are several
aspects worth discussing. First, clearly some pillars tend to increase in size
(red colors) while others tend to shrink (blue colors), creating a size
gradient that eventually will form the desired metagrating. Second, the
gradients are not symmetric along the circular boundary, indicating that the
shapes will tend to deviate from simple circles. Physically, this asymmetry
arises from the field discontinuities at regions where the linearly polarized
incident field is normal to the pillar surface. Finally, as a consequence of
the last point, TE and TM polarizations create different gradients, and tend
to deform the pillars differently. This is quite clearly seen in the second
pillar from the left, where the gradients seem somewhat rotated by $90\deg$.
This difference in gradients means that it is more challenging to optimize a
structure to simultaneously diffract both TE and TM polarizations with high
efficiencies. This argument can be extended to problems with multiple
objectives, for example optimizing a multi-wavelength device, where each
wavelength tends to create different forward and adjoint field distributions,
and therefore different shape gradients.

The gradients discussed in Figure 3b continuously deform the shape at each
iteration. As discussed, it is essential for manufacturability that certain
constraints are observed. The deformation function $u_{\perp}=g$ may in
general be very complicated, leading to complex shapes after many iterations.
Since we are dealing with a closed boundary, any function can be expanded in
terms of an appropriate basis defined in such domain. We chose Fourier as it
easily allows restriction to smooth round structures. Decomposing the gradient
function $g$ in a Fourier series, the boundary displacement is expressed as:

$u_{\perp}=h\,\text{sign}(\delta\epsilon)\left(\frac{a_{0}}{2}+\sum_{m=1}^{\infty}a_{m}\text{cos}\,m\theta+b_{m}\text{sin}\,m\theta\right),$
(4)

where $(a_{m},b_{m})$ are the expansion coefficients. With that, the
efficiency variation in equation 2 is then written as:

$\delta\eta=\frac{\omega\delta\epsilon}{4P}\,h\,s\,\text{sign}(\delta\epsilon)\left[\frac{a_{0}^{2}}{2}+\sum_{m=1}^{\infty}(a_{m}^{2}+b_{m}^{2})\right],$
(5)

Note that because we use a set of orthogonal basis (Fourier in this example),
each coefficient appears squared in the brackets. This means that each term of
the expansion independently contributes to increasing the metasurface
efficiency. We can freely choose to drop certain coefficients and still, the
remaining ones will always push the efficiency upwards (or remain unchanged if
a local optimum has been reached, but never reduce the efficiency). For
example, one might choose to keep only the zero-order ($m=0$) coefficient,
ensuring that the pillars remain circular throughout the optimization process
(of course their diameters are allowed to change). Controlling the Fourier
order $m$ gives explicit control of the trade-off between performance and
complexity at every iteration. To illustrate this point, we show in Figure 3c
the shape gradient and its decomposition for one of the pillars (highlighted
in Figure 3b). As one can see, restricting to zero-order $m=0$ simply
increases the pillar diameter. Adding the first-order term $m=1$ introduces a
shift in the center position of the pillar, while the second-order $m=2$
creates a certain ellipticity. Finally, adding up to the sixth $m=6$ order is
sufficient to closely represent the full gradient function for this particular
case. In the numerical examples presented in Section 3, we show various
devices designed with different Fourier orders, illustrating the ability
control the device complexity.

Another important constraint for fabrication is to ensure that the gap between
particles is not too small, which can severely intensify the issue of gap-
dependent etch rate, leading to under- or over-etching regions and sidewall
angle variations. In the shape optimization method, at every iteration we have
the explicit boundary coordinates $\bm{r_{i}}$ before deformation and
$\bm{r_{i+1}}=\bm{r_{i}}+u_{\perp}\bm{n_{i}}$ after deformation, where
$\bm{n_{i}}$ is the normal unit vector at a given point on the boundary. In
the simplest form, we can limit the scaling factor $h$ so that the deformed
boundary always respects a target $gap$ to the unit cell boundary. Finally, it
is also important to limit the minimum feature size to avoid challenges with
patterning, etching and pillars falling over, and again, this becomes
straightforward given direct knowledge of the boundary coordinates.

##  3 Numerical results

To illustrate the method, we applied shape optimization to design several
metagratings and metalenses at $1.55~{}\mu m$ operating wavelength, all based
on $1~{}\mu m$ tall amorphous silicon (aSi) pillars on a glass substrate. In
every example discussed here, we imposed minimum gap between pillars of 90 nm
and minimum feature size of 80 nm. Furthermore, we targeted polarization
insensitive operation and therefore simultaneously optimized for TE and TM
incident polarizations. In the first example, a meta-grating was designed to
deflect light at $51^{\circ}$. The initial structure was created based on the
library approach, and contains 4 circular pillars, each on a 500 nm unit cell.
In the library method, the pillars’ diameters are chosen so that they impart a
phase profile with steps of $2\pi/N$ ($N=4$ in this example). The phase
imparted by the first unit cell is, however, arbitrary, as different values
simply represent different global phases. One can therefore freely choose the
diameter of the first pillar as long as subsequent diameters correspond to a
$2\pi/N$ phase shift. Despite nominally imparting equivalent phase profiles,
these meta-gratings with different pillar sizes exhibit distinct non-local
interaction, and may perform very differently. This is exactly the case
observed here for the $51^{\circ}$, in which certain choices of the first
pillar diameter create a resonance that severely reduces the first order
diffraction efficiency for TE light. We applied shape optimization to two
metagratings differing only by the choice of the first pillar diameter, one
exhibiting low efficiency due to the appearance of a resonance for TE, and
another with the highest efficiency for this library (which was found after
sweeping the first pillar diameter for all values available in the library
from 100 nm to 410 nm).

Figure 4: Optimization of a dual-polarization 51-degree beam deflector. Shape
optimization was applied to two library designs differing only by the choice
of the first pillar diameter: in (a), the initial design exhibits low
efficiency due to the appearance of a resonance for TE polarization while in
(b) the initial structure exhibited the highest efficiency possible for such
library. The curves show the evolution of the absolute efficiency at each
iteration, where dashed curves represent restricting the gradient function to
$m=0$ Fourier order (i.e. maintaining circular shape), while solid line allows
up to $m=6$. The final structures for both Fourier orders are shown in the
insets. In (c), shape optimization was applied to design meta-gratings from 0
to 70 degrees (final geometries shown as insets). The color legend on top of
the figure applies to all plots.

The initial library design for the resonant metagrating has first order
diffraction efficiency of only about 30% for TE polarization (see iteration =
0 in Figure 4a), while it is as high as 85% for TM polarization. Such low
performance for TE cannot be predicted from the performance of the individual
meta-atoms, as all values of pillar diameters in the library are highly
transmissive (above 80%), and of course their respsonse is polarization
independent by symmetry. We then applied shape optimization and observed that
the performance is significantly improved as shown in Figure 4a. The
enhancement was observed either when the Fourier order is restricted to $m=0$
(dashed lines, final geometry in upper inset) or when we allowed up to $m=6$
(solid lines, final geometry in lower inset). Note that in both cases, at the
beginning of the optimization over the first few iterations, the improvement
in TE polarization comes at the expense of the TM efficiency. This initial
degradation in the TM efficiency could not be recovered with the $m=0$ case,
demonstrating a clear limitation in restricting the shapes to circles. With
higher Fourier orders, eventually the TM efficiency recovers leading to a
device with efficiencies of 82% and 84% for TE and TM, respectively. Despite
allowing higher Fourier order, the final shapes are still smooth and respect
the gap and minimum feature constraints. In the second case shown in Figure
4b, the initial structure was the highest average TE/TM efficiency obtained
with this library. Despite that, the initial efficiency of 76% for TE was
still significantly improved to 85%, with no degradation on the TM efficiency
(which actually slightly improved from 88% to 90%). This example illustrates
that the shape optimization can be used to eliminate resonances that are not
possible to predict from the library alone, as well as to push the performance
beyond the highest achievable with such library. Degradation in efficiency due
to non-local effects is not specific to this $51^{\circ}$ metagrating, and was
observed in other deflection angles from 10 to 70 degrees, in some cases TE
and in other cases for TM polarizations. Furthermore, a choice of initial
diameter that leads to relatively good performance for one deflection angle is
not necessarily the best choice for another angle, illustrating another
challenge of relying solely on the library approach. We designed various other
metagratings using shape optimzation from angles varying between 10 and 70
degrees, with the final efficiencies and shapes shown in Figure 4c. For
example, the 70 degrees metagrating showed average TE/TM efficiency of 74%.

As mentioned, the shape optimization method is not limited to metagratings and
can be applied to any general metasurface. To illustrate this, we designed a
metalens operating at $1.55~{}\mu m$ and the results are shown in Figure 5.
The optimization domain was again initialized with a structure based on a
library of aSi circular pillars on a glass substrate. The library unit cell is
$500~{}nm$ and the pillars range from $100$ to $410~{}nm$ in diameter, with
height of $1~{}\mu m$. The lens was designed to impart a phase transformation
$\phi=k(\sqrt{r^{2}+f^{2}}-f)$ with focus of $5.3~{}\mu m$ over a total
diameter of $2r=23~{}\mu m$, nominally resulting in a numerical aperture of
$NA=sin[tan^{-1}(r/f)]=0.92$. However, the lens was illuminated with a
Gaussian beam with waist radius $w=5.2~{}\mu m$, which focused at $f=5.3~{}\mu
m$ limits the NA of the system to approximately $NA=0.83$. We optimized two
metalenses that differ only by a global phase, represented by the diameter of
the innermost pillars of $d=250~{}nm$ and $d=300~{}nm$. As can be seen on
Figure 5a, these nominally equivalent metalenses have different efficiencies
of 75% and 81% (values at iteration 0), respectively. Shape optimization was
applied with Fourier decomposition up to $m=6$, and again with minimum gap of
$90~{}nm$ and minimum feature of $80~{}nm$. Given the circular symmetry of the
device, we only simulated one-quarter of the structure, and only one
polarization (TE). Polarization insensitivity is enforced by computing the
gradient for TM polarization from the TE gradient rotated by $90\deg$. The
optimization results in Figure 5a show that the efficiency is substantially
improved for both initial structures, reaching 88% and 90%.

Figure 5: (a) evolution of the absolute efficiency for two different initial
structures; cuts of the z-component of Poynting vector in the xz plane in (b)
and in the xy plane in (c). The color scale shown in (c) also applies for (b).
The full metalens structure is shown in (d), and a zoomed region at the
metalens center before and after shape optimizaiton is shown in (e).

These results reinforce two aspects already observed for metagratings: first,
shape optimization can improve the efficiency of library structures,
regardless of the specific initial choice. This is not to say that every
choice leads to the same final efficiency though (as is clear in Figure 5a).
Second, the efficiency is improved beyond what could be achieved with such
library. Figure 5a shows cuts of the z-component of Poynting vector in the xz
plane (z being the propagation direction), starting at $z=0$ (base of the
pillars) and in the free-space region (above the $1~{}\mu m$ pillar height)
with very little undesired diffraction observed. The beam at the focus plane
is shown in Figure 5c. The full metalens structure for the $d=300~{}nm$ case
is shown Figure 5d, and Figure 5e shows a zoomed region at the metalens center
before and after shape optimization. As can be seen, all features have smooth
boundaries, the minimum feature size is approximately $130~{}\mu m$, and
respecting the $80~{}nm$ minimum gap. As outlined in Appendix A, all values
reported here represent the absolute efficiency, defined as the projection of
the output field onto the target field divided by the power incident on the
metasurface.

##  4 Fabrication and characterization

To validate the method, various metagratings were fabricated and characterized
experimentally. The samples were fabricated using a conventional top-down
approach, allowing for high-throughput and reproducibility. Choosing amorphous
silicon as the nanopillars’ material allows for low losses and high contrast
refractive index across the telecommunications c-band. The metasurfaces were
fabricated on a $500~{}\mu m$ thick fused silica wafer, coated with $1~{}\mu
m$ aSi using PECVD. An adhesive layer (HDMS) is spun on the wafer to promote
adhesion of the negative e-beam resist (ma-N 2400) and the latter is then
baked and coated with a charging dissipating solution (e-spacer). The
nanopillars pattern is then exposed using e-beam lithography and developed in
AZ 726 developer, while the e-spacer layer is removed in water. Exploiting a
RIE (SPTS Rapier) technique, with simultaneous injection of etcher (SF6) and
passivation (C4F8) gases in the chamber, the pattern is transferred onto the
a-Si, using the resist as etching mask. The residual resist is then removed
with oxygen plasma asher (Matrix Plasma Asher). After fabrication, the samples
were characterized using a tuneable laser to measure the diffraction
efficiency over a wide spectral range for both TE and TM polarizations. All
values reported here represent the absolute diffraction efficiencies, defined
as the ratio between the output power in the first order diffraction and the
power incident on the metasurface, exactly how it is defined in the
simulations (often in the literature, relative diffraction efficiencies are
reported, and some times, even absolute efficiencies values don’t properly
consider Fresnel reflections at the glass-air interface). A detailed
description of the experimental setup and procedure is provided in Appendix B.

Figure 6: Simulated (a) and measured (b) TE spectral responses for the
resonant 51-degree beam deflector. Different colors show the first order
diffraction efficiency for the initial and shape optimized designs with
Fourier order $m=0$ and $m=6$ (the color legend applies to both (a) and (b)
plots). A scanning electron microscope image of the fabricated structures for
$m=6$ is shown in the inset (scale bar is 200 nm). In (c), the efficiency at
the design wavelength is plotted for TE and TM polarizations.

Figures 6a and b show the simulated and experimental spectral responses for
the resonant $51^{\circ}$ metagrating design from Figure 4a. We only show the
TE polarization as the efficiency for TM is relatively flat in this wavelength
region. The initial library design (in blue) exhibits a clear resonance at
$1.55~{}\mu m$. In the experiment, this resonance was slightly red-shifted to
$1.58~{}\mu m$, which we attribute to small variations in the fabrication
parameters such as refractive index, small sidewall angle and pillar sizes.
For the sake of comparison between experiment and simulations, all
efficiencies were then measured at the resonance wavelength, as indicated by
the solid dots in the figures. In both, simulations and experiments, the shape
optimized metagratings exhibited improved efficiencies. Comparing their
spectra, we observe that shape optimization improves the diffraction
efficiency by shifting by shifting the resonance away from the design
wavelength. Furthermore, one can see that the shape optimization induces a
tilt in the spectrum and the efficiency at the design wavelength, to the right
of the resonance, is further increased. Remarkably, these spectral signatures,
shift and tilt, are clearly observed in the experimental results in Figure 6b.
A comparison of the efficiencies predicted in simulations and observed
experimentally is shown in Figure 6c for both TE and TM, with reasonably good
agreement. It is also remarkable that the experimental results show the trade-
off predicted in the TE-TM efficiency for the two different Fourier orders:
simple circles ($m=0$) increase the TE at the expense of TM while $m=6$ breaks
this trade-off.

Figure 7: Simulated (a) and measured (b) TE spectral responses for the non-
resonant 51-degree beam deflector. Different colors show the first order
diffraction efficiency for the initial and shape optimized designs with
Fourier order $m=0$ and $m=6$ (the color legend applies to both (a) and (b)
plots). An scanning electron microscope image of the fabricated structures for
$m=6$ is shown in the inset (scale bar is 300 nm). In (c), the efficiency at
the design wavelength is plotted for TE and TM polarizations.

We also fabricated the non-resonant $51^{\circ}$ meta-grating design from
Figure 4b, and the results for TE are shown in Figure 7. As mentioned, this
was the highest efficiency obtained with this library and one can see that
differently from the previous case, it pushed away the resonance approximately
34 nm above the design wavelength. Again the spectrum is slightly shifted in
the experiment, and therefore all measurements were shifted accordingly.
Differently than in the previous resonant case, the shape optimization now
does not lead to significant shifts in the resonance, and mostly push the
spectrum upwards at the design wavelength. The experimental results exhibit
the same behavior as predicted in simulations, and a quantitative comparison
in the efficiencies is shown in Figure 7c. Similarly to the resonant case,
here we also observe that the higher Fourier order $m=6$ exhibits high
efficiency for both TE and TM, reaching about 83% and 84% absolute
efficiencies (higher than both the library design and the $m=0$ device).

Figure 8: Simulated (a) and measured (b) spectral responses for both TE and TM
polarizations for a meta-gratings at 70-degree deflection angle. Different
colors show the first order diffraction efficiency for the initial and shape
optimized designs with Fourier order $m=6$ (the color legend applies to both
(a) and (b) plots). An scanning electron microscope image of the fabricated
structures is shown in the inset (scale bar is 200 nm). In (c), the efficiency
at the design wavelength is plotted for TE and TM polarizations.

Finally, we then designed and fabricated a $70^{\circ}$ metagrating, showing
broadband and polarization insensitive operation. The spectral efficiencies
for the initial library design and for the shape optimized metagratings are
shown in Figure 8 in dashed and solid lines, respectively. The experimental
results in (b) reproduce qualitatively the simulations in (a), with measured
peak efficiency above 70%. The shape optimized meta-gratings with Fourier
order $m=6$ exhibit a relatively flat spectrum centered at $1.55~{}\mu m$,
with substantially higher efficiencies than the library for both TE and TM. A
quantitative comparison between the simulation and experimental efficiencies
is shown in Figure 8c.

##  5 Discussion

The shape optimization method discussed in this paper presents an alternative
to conventional library approach in metasurface design and to the more general
topology optimization. As discussed, despite being simple and computationally
fast, the library method falls short in many aspects, most notably by not
capturing non-local coupling effects and not being amenable to general
waveform incidences. Our results demonstrate that library structures
exhibiting low efficiency due to the aforementioned limitations can be
substantially improved with shape optimization. We show results where a
resonant library structure with only 30% TE efficiency was enhanced to 82%,
while still maintaining high TM efficiency at 84%, smooth boundaries and
respecting minimum features size and minimum gaps. Furthermore, we showed that
even the highest efficiency library is also significantly improved going from
76% to 85% for TE, while maintaining high efficiency for TM at 90%, and again
respecting fabrication constraints. The general shape optimization explores
more degrees of freedom than parametrized structures, and it is not therefore
surprising that it leads to higher efficiencies. For example, parametrized aSi
elliptical pillars were used to design 50 deg gratings at 900 nm wavelength
and obtained absolute efficiency of 64% for a single polarization [36].
Comparatively, shape optimization leads to TE and TM efficiency of 85% and
90%. Using rectangular pillars, a metalens with 0.78 numerical aperture
operating at 850 nm was designed using parametrized rectangular aSi pillars,
and obtained efficiency of 78% (no minimum feature size or gap constraints).
In contrast, shape optimized metalens with slightly higher NA of 0.83 resulted
in efficiency of 90%, including fabrication constraints.

Compared to topology optimization, shape optimization explores a reduced
design space. In topology, the output structure can contain many more detailed
features, and it most likely leads to higher overall efficiencies than shape
optimization. For example, topology optimized gratings showed theoretical
efficiencies for a 70 degrees deflection angle around 96% (average TE/TM)
[18]. Comparatively, shape optimization obtained lower values of 74% average
TE/TM as shown in Figure 4c. However, as discussed such topology optimized
structures are more difficult to fabricate than shape optimized structures.
Another topology optimized metagrating at 75 degrees deflection angle that was
actually fabricated showed experimental efficiencies of 74% and 75%, for TE
and TM, respectively [18]. The definition of efficiency here was the power in
the deflected beam normalized to the power transmitted through the bare
silicon dioxide substrate. This definition differs from ours in the sense that
we compute efficiency as the power in the deflected beam normalized by the
power incident on the metasurface. The difference is the Fresnel reflection,
which if adjusted brings the measured values to 71% and 72%. Comparatively,
the measured efficiencies for the 70 degrees shape optimized structures from
Figure 8 were on a similar level at 70%.

As a final note, the shape optimization algorithm was demonstrated here for
two objective functions, i.e., TE and TM polarization, but it can be easily
extended to more objectives such as multiple wavelengths or multiple angles of
incidence etc. The computational cost of this optimization method is
expensive, as it requires solving Maxwell equations for both forward and
adjoint fields. In our optimization, we used a commercial finite-element
solver in the frequency domain. Multi-wavelength optimization can greatly
benefit from finite-difference time domain solvers, as the fields for all
wavelengths can be obtained from only one forward and one adjoint simulations.
Finally, combining this approach with GPU-accelerated FDTD can greatly speed
up solving the fields [38, 39]. In general, the adjoint formulation tends to
converge to a local optimum, with different initial conditions converging to
different local optima, as was illustrated in the examples discussed here.
Coupling the optimization method with parallelization is therefore beneficial
to allow many starting geometries to be explored efficiently.

In summary, we presented a general shape optimization method that enables
optimization of high efficiency metasurfaces. Coupled with a Fourier
decomposition of the boundary gradient, the method enables high performance
devices while providing greater control over the structure complexity and
rigorously enforcing minimum feature size and minimum gaps. Various
metagratings and metalenses were simulated, with experimental results
validating the expected efficiencies. We believe these results provide a path
towards manufacturability of inverse-designed high efficiency metasurfaces.

## References

  * [1] Nanfang Yu and Federico Capasso  “Flat optics with designer metasurfaces”  In _Nature materials_ 13.2 Nature Publishing Group UK London, 2014, pp. 139–150
  * [2] Alexander V Kildishev, Alexandra Boltasseva and Vladimir M Shalaev  “Planar photonics with metasurfaces”  In _Science_ 339.6125 American Association for the Advancement of Science, 2013, pp. 1232009
  * [3] Amir Arbabi, Yu Horie, Mahmood Bagheri and Andrei Faraon  “Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission”  In _Nature nanotechnology_ 10.11 Nature Publishing Group UK London, 2015, pp. 937–943
  * [4] Arseniy I Kuznetsov, Andrey E Miroshnichenko, Mark L Brongersma, Yuri S Kivshar and Boris Luk’yanchuk  “Optically resonant dielectric nanostructures”  In _Science_ 354.6314 American Association for the Advancement of Science, 2016, pp. aag2472
  * [5] Philippe Lalanne and Pierre Chavel  “Metalenses at visible wavelengths: past, present, perspectives”  In _Laser & Photonics Reviews_ 11.3 Wiley Online Library, 2017, pp. 1600295
  * [6] Arseniy I Kuznetsov, Mark L Brongersma, Jin Yao, Mu Ku Chen, Uriel Levy, Din Ping Tsai, Nikolay I Zheludev, Andrei Faraon, Amir Arbabi and Nanfang Yu  “Roadmap for optical metasurfaces”  In _ACS Photonics_ ACS Publications, 2024
  * [7] Wei Ting Chen, Alexander Y Zhu, Vyshakh Sanjeev, Mohammadreza Khorasaninejad, Zhujun Shi, Eric Lee and Federico Capasso  “A broadband achromatic metalens for focusing and imaging in the visible”  In _Nature nanotechnology_ 13.3 Nature Publishing Group UK London, 2018, pp. 220–226
  * [8] Ehsan Arbabi, Seyedeh Mahsa Kamali, Amir Arbabi and Andrei Faraon  “Full-Stokes imaging polarimetry using dielectric metasurfaces”  In _Acs Photonics_ 5.8 ACS Publications, 2018, pp. 3132–3140
  * [9] Noah A Rubin, Gabriele D’Aversa, Paul Chevalier, Zhujun Shi, Wei Ting Chen and Federico Capasso  “Matrix Fourier optics enables a compact full-Stokes polarization camera”  In _Science_ 365.6448 American Association for the Advancement of Science, 2019, pp. eaax1839
  * [10] Ahmed H Dorrah and Federico Capasso  “Tunable structured light with flat optics”  In _Science_ 376.6591 American Association for the Advancement of Science, 2022, pp. eabi6860
  * [11] Jaewon Oh, Kangmei Li, Jun Yang, Wei Ting Chen, Ming-Jun Li, Paulo Dainese and Federico Capasso  “Adjoint-optimized metasurfaces for compact mode-division multiplexing”  In _ACS photonics_ 9.3 ACS Publications, 2022, pp. 929–937
  * [12] Jaewon Oh, Jun Yang, Louis Marra, Ahmed H Dorrah, Alfonso Palmieri, Paulo Dainese and Federico Capasso  “Metasurfaces for free-space coupling to multicore fibers”  In _Journal of Lightwave Technology_ IEEE, 2023
  * [13] Mu Ku Chen, Yue Yan, Xiaoyuan Liu, Yongfeng Wu, Jingcheng Zhang, Jiaqi Yuan, Zhengnan Zhang and Din Ping Tsai  “Edge detection with meta-lens: from one dimension to three dimensions”  In _Nanophotonics_ 10.14 De Gruyter, 2021, pp. 3709–3715
  * [14] Andrea Cordaro, Brian Edwards, Vahid Nikkhah, Andrea Alù, Nader Engheta and Albert Polman  “Solving integral equations in free space with inverse-designed ultrathin optical metagratings”  In _Nature Nanotechnology_ 18.4 Nature Publishing Group UK London, 2023, pp. 365–372
  * [15] Sindhu Jammi, Andrew R Ferdinand, Zheng Luo, Zachary L Newman, Gregory Spektor, Junyeob Song, Okan Koksal, Akash V Rakholia, William Lunden and Daniel Sheredy  “Three-dimensional, multi-wavelength beam formation with integrated metasurface optics for Sr laser cooling”  In _arXiv preprint arXiv:2402.08885_ , 2024
  * [16] Ming Lun Tseng, Michael Semmlinger, Ming Zhang, Catherine Arndt, Tzu-Ting Huang, Jian Yang, Hsin Yu Kuo, Vin-Cent Su, Mu Ku Chen and Cheng Hung Chu  “Vacuum ultraviolet nonlinear metalens”  In _Science advances_ 8.16 American Association for the Advancement of Science, 2022, pp. eabn5644
  * [17] Philippe Lalanne, Jean Paul Hugonin and Pierre Chavel  “Optical properties of deep lamellar gratings: a coupled Bloch-mode insight”  In _Journal of Lightwave Technology_ 24.6 IEEE, 2006, pp. 2442
  * [18] David Sell, Jianji Yang, Sage Doshay, Rui Yang and Jonathan A Fan  “Large-angle, multifunctional metagratings based on freeform multimode geometries”  In _Nano letters_ 17.6 ACS Publications, 2017, pp. 3752–3757
  * [19] Sawyer D Campbell, David Sell, Ronald P Jenkins, Eric B Whiting, Jonathan A Fan and Douglas H Werner  “Review of numerical optimization techniques for meta-device design”  In _Optical Materials Express_ 9.4 Optica Publishing Group, 2019, pp. 1842–1863
  * [20] Wilhelm Stork, Norbert Streibl, H Haidner and P Kipfer  “Artificial distributed-index media fabricated by zero-order gratings”  In _Optics letters_ 16.24 Optica Publishing Group, 1991, pp. 1921–1923
  * [21] Frederick T Chen and Harold G Craighead  “Diffractive lens fabricated with mostly zeroth-order gratings”  In _Optics letters_ 21.3 Optica Publishing Group, 1996, pp. 177–179
  * [22] Philippe Lalanne, Simion Astilean, Pierre Chavel, Edmond Cambril and Huguette Launois  “Blazed binary subwavelength gratings with efficiencies larger than those of conventional échelette gratings”  In _Optics letters_ 23.14 Optica Publishing Group, 1998, pp. 1081–1083
  * [23] Nanfang Yu, Patrice Genevet, Mikhail A Kats, Francesco Aieta, Jean-Philippe Tetienne, Federico Capasso and Zeno Gaburro  “Light propagation with phase discontinuities: generalized laws of reflection and refraction”  In _Science_ 334.6054 American Association for the Advancement of Science, 2011, pp. 333–337
  * [24] Eric B Whiting, Sawyer D Campbell, Lei Kang and Douglas H Werner  “Meta-atom library generation via an efficient multi-objective shape optimization method”  In _Optics Express_ 28.16 Optica Publishing Group, 2020, pp. 24229–24242
  * [25] Jakob Søndergaard Jensen and Ole Sigmund  “Topology optimization for nano-photonics”  In _Laser & Photonics Reviews_ 5.2 Wiley Online Library, 2011, pp. 308–321
  * [26] Sean Molesky, Zin Lin, Alexander Y Piggott, Weiliang Jin, Jelena Vucković and Alejandro W Rodriguez  “Inverse design in nanophotonics”  In _Nature Photonics_ 12.11 Nature Publishing Group UK London, 2018, pp. 659–670
  * [27] Haejun Chung and Owen D Miller  “High-NA achromatic metalenses by inverse design”  In _Optics Express_ 28.5 Optica Publishing Group, 2020, pp. 6945–6965
  * [28] Di Sang, Mingfeng Xu, Mingbo Pu, Fei Zhang, Yinghui Guo, Xiong Li, Xiaoliang Ma, Yunqi Fu and Xiangang Luo  “Toward high-efficiency ultrahigh numerical aperture freeform metalens: from vector diffraction theory to topology optimization”  In _Laser & Photonics Reviews_ 16.10 Wiley Online Library, 2022, pp. 2200265
  * [29] Evan W Wang, David Sell, Thaibao Phan and Jonathan A Fan  “Robust design of topology-optimized metasurfaces”  In _Optical Materials Express_ 9.2 Optica Publishing Group, 2019, pp. 469–482
  * [30] Victor Liu and Shanhui Fan  “Compact bends for multi-mode photonic crystal waveguides with high transmission and suppressed modal crosstalk”  In _Optics express_ 21.7 Optica Publishing Group, 2013, pp. 8069–8075
  * [31] Christopher M Lalau-Keraly, Samarth Bhargava, Owen D Miller and Eli Yablonovitch  “Adjoint shape optimization applied to electromagnetic design”  In _Optics express_ 21.18 Optica Publishing Group, 2013, pp. 21693–21701
  * [32] Nicolas Lebbe, Charles Dapogny, Edouard Oudet, Karim Hassan and Alain Gliere  “Robust shape and topology optimization of nanophotonic devices using the level set method”  In _Journal of Computational Physics_ 395 Elsevier, 2019, pp. 710–746
  * [33] Mahdad Mansouree, Andrew McClung, Sarath Samudrala and Amir Arbabi  “Large-scale parametrized metasurface design using adjoint optimization”  In _ACS Photonics_ 8.2 ACS Publications, 2021, pp. 455–463
  * [34] Erez Gershnabel, Mingkun Chen, Chenkai Mao, Evan W Wang, Philippe Lalanne and Jonathan A Fan  “Reparameterization Approach to Gradient-Based Inverse Design of Three-Dimensional Nanophotonic Devices”  In _ACS Photonics_ 10.4 ACS Publications, 2022, pp. 815–823
  * [35] You Zhou, Chenkai Mao, Erez Gershnabel, Mingkun Chen and Jonathan A Fan  “Large-Area, High-Numerical-Aperture, Freeform Metasurfaces”  In _Laser & Photonics Reviews_ Wiley Online Library, 2024, pp. 2300988
  * [36] You Zhou, Yixuan Shao, Chenkai Mao and Jonathan A Fan  “Inverse-designed metasurfaces with facile fabrication parameters”  In _Journal of Optics_ 26.5 IOP Publishing, 2024, pp. 055101
  * [37] Lev Davidovich Landau, John Stewart Bell, MJ Kearsley, LP Pitaevskii, EM Lifshitz and JB Sykes  “Electrodynamics of continuous media”  elsevier, 2013
  * [38] Tyler W Hughes, Momchil Minkov, Victor Liu, Zongfu Yu and Shanhui Fan  “A perspective on the pathway toward full wave simulation of large area metalenses”  In _Applied Physics Letters_ 119.15 AIP Publishing, 2021
  * [39] Jinhie Skarda, Rahul Trivedi, Logan Su, Diego Ahmad-Stein, Hyounghan Kwon, Seunghoon Han, Shanhui Fan and Jelena Vučković  “Low-overhead distribution strategy for simulation and optimization of large-area metasurfaces”  In _npj Computational Materials_ 8.1 Nature Publishing Group UK London, 2022, pp. 78
  * [40] Steven G Johnson, Mihai Ibanescu, MA Skorobogatiy, Ori Weisberg, JD Joannopoulos and Yoel Fink  “Perturbation theory for Maxwell’s equations with shifting material boundaries”  In _Physical review E_ 65.6 APS, 2002, pp. 066611
  * [41] K-M Chen  “A mathematical formulation of the equivalence principle”  In _IEEE Transactions on Microwave Theory and Techniques_ 37.10 IEEE, 1989, pp. 1576–1581

##  Appendix A Mathematical formulation

The adjoint formulation is discussed in A.1, including the definition of the
boundary gradient. A validation of the gradient with two numerical examples is
provided in A.2, and the Fourier decomposition is discussed in A.3. Finally, a
generalization to an arbitrary figure of merit is discussed in A.4.

###  A.1 Adjoint-based boundary gradient

The objective when designing a metasurface is to determine a specific
configuration in the constituent materials, meaning their geometrical shapes
and dielectric permittivity, so that an incoming field is transformed into a
target field $E_{t}$ at the output surface $S$ (Figure 9a). Suppose we start
with an arbitrary configuration represented by a dielectric permittivity
distribution $\epsilon$, which in general varies in space, and then change it
from $\epsilon^{\prime}=\epsilon+\delta\epsilon$ by modifying a certain region
$\chi$ within the optimization domain (Figure 9b). For example, $\epsilon$ may
represent an arrangement of meta-atoms (with permittivity $\epsilon_{2}$)
embedded in a background material (with permittivity $\epsilon_{1}$).
Enlarging one of the meta-atoms changes the permittivity in the region $\chi$
from the background value $\epsilon_{1}$ to the meta-atom’s permittivity
$\epsilon_{2}$. In general, a set of perturbations create a new configuration
represented by $\epsilon^{\prime}$. The question we want to answer is where
and how to change the structure so that the field changes from $E$ to
$E^{\prime}=E+\delta E$, moving in a direction closer to our target field
$E_{t}$ at the output surface $S$.

Figure 9: (a) an incident field interacts with a material configuration with
permittivity $\epsilon$ to generate a field $E$ in space; (b) a new
configuration $\epsilon^{\prime}=\epsilon+\delta\epsilon$ generates a field
$E^{\prime}=E+\delta E$ (note that the region where the permittivity changed
is denoted by $\chi$). In (c), the same field correction $\delta E$ can be
produced in the initial configuration $\epsilon$ when a proper current
$J_{\delta}$ is added in the perturbed region $\chi$.

Typically, a figure of merit is used to ‘measure’ how close the field $E$ is
to the target $E_{t}$. For example, $E_{t}$ might be the mode of a waveguide
positioned at the output of the domain or simply a desired field in free-space
(such as a particular diffraction order in a grating or a focused Gaussian
beam and so on). We can express the projection of the field $E$ into the
target field $E_{t}$ as:

$F=\frac{1}{4\ P}\int_{S}{\left(E\times H_{t}^{\ast}+E_{t}^{\ast}\times
H\right)\cdot n}da,$ (6)

where $n$ is the normal unit vector at the output surface S, pointing outwards
from the simulation domain. If we normalize both the incoming and target
fields ($P$=1 W), then the efficiency of our device is directly obtained by
$\eta=\left|F\right|^{2}$. Any modification in the metasurface structure that
changes the fields by $\left(\delta E,\delta H\right)$ impacts the efficiency
as:

$\delta\eta=2\ \mathrm{\mathrm{Re}}\left[F^{\ast}\delta F\right],$ (7)

where

$\displaystyle\delta F$ $\displaystyle=\frac{1}{4\ P_{t}}\int_{S}{\left(\delta
E\times H_{t}^{\ast}+E_{t}^{\ast}\times\delta H\right)\cdot n}da$ (8)
$\displaystyle=\frac{1}{4\ P_{t}}\int_{S}\left[-\left(n\times
H_{t}^{\ast}\right)\cdot\delta E+\left(n\times E_{t}^{\ast}\right)\cdot\delta
H\right]da.$

The challenge is to determine how to modify the metasurface geometry in each
iteration so $\delta\eta$ increases, as opposed to brute force trial and error
every possible modification. The derivation goes as follows: we first perform
two simulations to extract information in the existing medium configuration,
so called forward and adjoint simulations. In the forward simulation, we
calculate the field $E$ for the unperturbed medium $\epsilon$ excited by our
input source. This allows us to estimate any new polarization currents
$J_{\delta}$ that might be generated when the medium is modified by
$\delta\epsilon$ (i.e., any small amount of material placed on the surface
that generates a new polarization current), as illustrated in Figure 9c. In
the adjoint simulation, we place specific currents sources $J_{a}$ on the
output surface $S$ that excites our medium backwards, generating a field in
our domain denoted by $E_{a}$. We choose the currents such the adjoint field
is equivalent to propagating our target field backwards into the medium. The
next step is to relate these currents $J_{a}$ and $J_{\delta}$ using the
reciprocity theorem. This will give us a recipe for choosing where to place
$J_{\delta}$ (i.e. where to place $\delta\epsilon$) that ensures our
efficiency increases. We now follow these steps mathematically.

The field correction $\delta E$ obtained in the perturbed medium
$\epsilon^{\prime}$ can be generated in the ‘known’ configuration $\epsilon$
if a proper current $J_{\delta}$ is placed in the perturbed regions $\chi$,
where:

$J_{\delta}=j\omega\delta\epsilon\ E^{\prime},$ (9)

This is illustrated in Figure 9c. Physically, this current is the time
derivative of the additional polarization in the medium
$J_{\delta}=\partial_{t}\delta P=\partial_{t}\delta\epsilon
E^{\prime}=j\omega\delta\epsilon E^{\prime}$, where we assume all fields are
harmonics with dependency $e^{j\omega t}$. In other words, $J_{\delta}$ are
the polarization currents in the medium. This can be seen by directly writing
down Ampère’s law in a source free region in both perturbed and unperturbed
medium:

$\nabla\times\ H=j\omega\epsilon\ E,\text{ and }\nabla\times\
H^{\prime}=j\omega\epsilon^{\prime}E^{\prime},$ (10)

where again $\epsilon^{\prime}=\epsilon+\delta\epsilon$ and
$E^{\prime}=E+\delta E$. Subtracting one from the other we obtain directly:

$\nabla\times\delta\ H=j\omega\epsilon\delta\ E+J_{\delta}.$ (11)

The last equation states that the fields $\left(\delta E,\delta H\right)$ are
generated in the unperturbed medium $\epsilon$ by the current source
$J_{\delta}$. At first this expression does not seem to be very useful because
we do not yet know the perturbed field $E^{\prime}$ (and therefore we can’t
determine $J_{\delta}=j\omega\delta\epsilon E^{\prime}$). However, note that
we do not need to know $E^{\prime}$ everywhere in space, we only need to find
an approximation for $E^{\prime}$ in the perturbation regions $\chi$, i.e.,
where $\delta\epsilon$ is non-zero. This is illustrated in Figure 10. The
dashed line represents the original boundary of a given feature (a circular
pillar in this case) that is then distorted by an amount $u_{\bot}$, which of
course can vary along the perimeter of the shape. The tangential component of
the electric field is continuous across the boundary and therefore can be
considered approximately unchanged if the boundary displacement is
infinitesimal (i.e. $E_{\parallel}^{\prime}=E_{\parallel}$ for
$u_{\bot}\rightarrow 0$). This is however not true for the normal electric
field. Since it is discontinuous, the field change is finite no matter how
small the boundary displacement is [40]. We can circumvent this problem by
noting that the normal displacement field $D_{\bot}$ is continuous, and
therefore:

$E^{\prime}\cong\ E_{\parallel}+\frac{1}{\epsilon^{\prime}}D_{\bot},$ (12)

and so

$J_{\delta}\cong\
j\omega\delta\epsilon\left(E_{\parallel}+\frac{1}{\epsilon\prime}D_{\bot}\right),$
(13)

Note that this approximation treats correctly abrupt material boundaries,
i.e., with discontinuous permittivity across the boundary [40]. In most
topological optimization algorithms, the medium is represented by an
approximate continuous permittivity, in which case the field is always
continuous and therefore $E^{\prime}\cong E$. This means that in topology
optimization, the material slowly evolves from a continuous distribution,
eventually converging to a binary index distribution (i.e. manufacturable
medium). This evolution however makes it difficult to control appearances of
new boundaries, which can be dealt directly with shape optimization as
discontinuous boundaries are naturally treated.

Figure 10: an initial boundary (dashed line) is deformed by a normal
displacement field $u_{\perp}$. In regions where $u_{\perp}$ is positive the
permittivity changes by $\delta\epsilon=\epsilon_{2}-\epsilon_{1}$, while
where $u_{\perp}$ is negative $\delta\epsilon=\epsilon_{1}-\epsilon_{2}$. In
the perturbed regions (i.e. near the boundaries), the corrected electric field
can be estimated by considering the continuous tangential electric
$E_{\parallel}$ and normal displacement $D_{\perp}$ fields, obtained in the
original configuration.

For the adjoint simulation, we illuminate the structure with an incident field
$\left(-E_{t}^{\ast},H_{t}^{\ast}\right)$, where the minus reflects the
reversed propagation direction (i.e. propagating into the domain). This can be
done through the Principle of Equivalence [41], which states that any incident
field into a domain can be produced by equivalent electric and magnetic
currents $J_{a}$ and $K_{a}$ on the surface $S$ enclosing this domain, where:

$J_{a}=-n\times H_{t}^{\ast}\text{, and }K_{a}=-n\times E_{t}^{\ast},$ (14)

where again these surface currents produce a field in our unperturbed domain
denoted by $E_{a}$. The interpretation of this is clear, it is just the
response of the metasurface when excited by our target field instead of our
original incident field.

These considerations result in two current-field pairs with a meaningful
physical interpretation: the first pair comes from the forward simulation,
where we determine that polarization currents $J_{\delta}$ located in the ‘to-
be-perturbed’ regions $\chi$ produces the field correction $\left(\delta
E,\delta H\right)$. The second pair comes from the adjoint simulation, where
the currents $\left(J_{a},K_{a}\right)$ located on the domain surface $S$
produce the adjoint field $\left(E_{a},H_{a}\right)$. The goal is now to find
a relation between these two pairs that tells us where the currents
$J_{\delta}$ should be placed so that $E+\delta E\rightarrow E_{t}$. This is
accomplished by evoking the Lorentz reciprocity theorem [37], which states
that in a medium characterized by symmetric $(\epsilon,\mu)$, the fields
$(E_{1},H_{1})$ and $(E_{2},H_{2})$, respectively produced by two different
sets of localized currents $(J_{1},K_{1})$ and $(J_{2},K_{2})$, satisfy the
relationship:

$\int_{\Omega}\left(J_{1}\cdot E_{2}-K_{1}\cdot
H_{2}\right)dv=\int_{\Omega}\left(J_{2}\cdot E_{1}-K_{2}\cdot H_{1}\right)dv.$
(15)

Using the Dirac’s delta to represent our surface current as volume currents
$\left(J_{a},K_{a}\right)\ast\delta S$, where $\delta S$ is surface
differential, the volume integrals involving the adjoint currents become
surface integrals, and therefore the Lorentz reciprocity applied to our pairs
of current-field becomes:

$\int_{S}\left(J_{a}\cdot\delta E-K_{a}\cdot\delta
H\right)da=\int_{\chi}(J_{\delta}\cdot E_{a})dv.$ (16)

Note that the volume integrals involving $J_{\delta}$ is performed only in the
region $\chi$, because only there $J_{\delta}$ is non-zero. Furthermore,
because of our choice for the adjoint currents, we can express the
differential in the figure of merit as:

$\displaystyle\delta F$ $\displaystyle=\frac{1}{4\
P}\int_{S}\left[-\left(n\times H_{t}^{\ast}\right)\cdot\delta E+\left(n\times
E_{t}^{\ast}\right)\cdot\delta H\right]da$ (17) $\displaystyle=\frac{1}{4\
P}\int_{S}\left[J_{a}\cdot\delta E-K_{a}\cdot\delta H\right]da$
$\displaystyle=\frac{1}{4\ P}\int_{\chi}{(J_{\delta}\cdot E_{a})}dv.$

Using the expression for $J_{\delta}$ and writing
$E_{a}=E_{a,\parallel}+\frac{1}{\epsilon_{2}}D_{a,\bot}$(which is exact), we
obtain:

$\delta F=\frac{j\omega\delta\epsilon}{4\
P}\int_{\chi}{\left(E_{\parallel}\cdot E_{a,\parallel}+\frac{1}{\epsilon\
\epsilon\prime}D_{\bot}\cdot D_{a,\bot}\right)}dv.$ (18)

The reciprocity theorem allows us to predict the change in the figure of merit
by performing an integration involving quantities directly at the metasurface
region, and not on the output surface S. This also means that maximizing this
integral in $\chi$ directly maximizes our figure of merit. The integration in
the volume $\chi$ can be mapped to an integration on the surface of the
original shape. Say $s$ denotes the arc length along the closed boundary
$\Gamma$ encircling the nanopillar and $u_{\bot}$ is the normal component of
the boundary displacement (Figure 10). The volume element is therefore
$dv=u_{\bot}dsdz$. Assuming the shape only changes in the x,y plane (i.e. all
pillars are uniform along the propagation direction z), then:

$\delta F=\frac{j\omega\delta\epsilon}{4\
P}\int{u_{\bot}\left(E_{\parallel}\cdot E_{a,\parallel}+\frac{1}{\epsilon\
\epsilon\prime}D_{\bot}\cdot D_{a,\bot}\right)}dsdz,$ (19)

and the change in efficiency can be written finally as:

$\delta\eta=\frac{\omega\delta\epsilon}{2\ P}\ \int_{\Gamma}{(u_{\bot}g)}ds,$
(20)

where we define a gradient function
$g=Re\left[jF^{\ast}\int{\left(E_{\parallel}\cdot
E_{a,\parallel}+\frac{1}{\epsilon\ \epsilon\prime}D_{\bot}\cdot
D_{a,\bot}\right)}dz\right]$, which includes an integration along the pillar
height. The efficiency change is always positive if we choose $u_{\bot}=h\
sign\left(\delta\epsilon\right)\ g$, where $h$ is a scaling factor.

With that, the algorithm is summarized as follows:

  * •

Initialize the structure with a set of known geometrical features (meta-
atoms), with any shape. For example, a set of uniform circular pillars or a
metasurface previously design with a library approach;

  * •

Compute the forward $E$ and adjoint $E_{a}$ fields;

  * •

With these fields, calculate the gradient function $g$ along the boundary of
every pillar;

  * •

Update the boundary shape by displacing every point in its perimeter by
$u_{\bot}=h\ sign\left(\delta\epsilon\right)\ g$. Go back to step 2 until the
efficiency reaches a local maximum;

The scaling factor $h$ does not alter the deformation function $u_{\bot}$ but
simply scales how much the overall shape is deformed. Obviously if $h$ is too
large, the approximation for $E^{\prime}$ and thus for $J_{\delta}$ eventually
becomes a poor one. There are many ways in which $h$ could be chosen. A simple
method is to first calculate $h$ so that the efficiency change is bounded (say
no more than 1% or no less than 0.01% or any other value). Then, one might
impose the maximum of $u_{\bot}$ (say less than 5 nm) and re-scale $h$. A more
general method to deform the boundary is using level-set functions.

###  A.2 Gradient validation

To validate the gradients obtained using the adjoint formulation, we computed
how the first order diffraction efficiency change as the parameters of the
geometries are varied. In Figure 11a, the diffraction efficiency is plotted as
a function of the width w of a rectangular ridge for both TE (y-pol) and TM
(x-pol) polarized light. The derivatives were then numerically calculated and
plotted in solid lines in Figure 11c. For comparison, the derivatives obtained
directly from the adjoint gradient are shown in circles, in excellent
agreement with the exact ones. In Figure 11b, the efficiency of a metasurface
that contains two circular pillars is calculated as a function of $d_{2}$,
diameter of the second pillar (while $d_{1}$ is fixed). Again, the derivates
computed from the efficiency curves (solid lines) match the one obtained from
the adjoint gradient, as shown in Figure 11d.

Figure 11: Validation of the boundary gradient calculation. In (a), the
efficiency of a metagrating containing only one rectangular nanopillar is
numerically calculated as a function of the width w, and in (b) the efficiency
of a metasurface with two circular pillars is calculated as a function of the
diameter of one of them ($d_{2}$). Here, blue and orange represent efficiency
for x-polarized (along larger domain direction) and y-polarized light (along
smaller domain direction), respectively. From the efficiency, the exact
gradient (solid line) is calculated for the two metasurfaces and compared with
the gradient obtained from the boundary shift expression (open dots). Top and
bottom figures share the same horizontal axis.

###  A.3 Fourier decomposition

It is often desirable to restrict the shape to simple smooth geometries, for
example to simplify fabrication. The deformation function $u_{\perp}$
calculated as outlined here may be very complicated, with sharp peaks or
valleys. Since we are dealing with a closed boundary, any function can be
expanded in terms of an appropriate basis. We chose Fourier as it easily
allows restriction to smooth round structures (the arguments presented here
can be generalized to any basis). The gradient function $g$ can be expanded in
terms of its Fourier components:

$g=\frac{a_{0}}{2}+\sum_{m=1}^{\infty}{a_{m}\cos{m\theta}+b_{m}\sin{m\theta}},$
(21)

where

$a_{m}=\frac{1}{\pi}\int_{-\pi}^{\pi}{g\cos{m\theta}d\theta},\\\
b_{m}=\frac{1}{\pi}\int_{-\pi}^{\pi}{g\sin{m\theta}d\theta},\text{ and}\\\
\theta=2\pi\frac{s}{s_{t}}.$ (22)

Here, $s_{t}$ is the total boundary length of a given meta-atom. A circular
shape is maintained if we retain only the zero-order term $a_{0}$:
$u_{\bot}=h\ sign\left(\delta\epsilon\right)\ \frac{a_{0}}{2}$. Retaining
first order terms represent slight displacements, second order terms represent
elliptical distortions and so on. In general:

$u_{\bot}=h\ sign\left(\delta\epsilon\right)\
\left(\frac{a_{0}}{2}+\sum_{m=1}^{\infty}{a_{m}\cos{m\theta}+b_{m}\sin{m\theta}}\right)$
(23)

and therefore, the change in efficiency is:

$\delta\eta=\frac{\omega\delta\epsilon}{4\ P}\
hs_{t}sign\left(\delta\epsilon\right)\left[\frac{a_{0}^{2}}{2}+\sum_{m=1}^{\infty}\left(a_{m}^{2}+b_{m}^{2}\right)\right]$
(24)

Note that due to the orthogonality of the Fourier basis, each term contributes
independently to an increase in efficiency. Therefore, we can choose to retain
a finite number of terms. Finally, if we want to optimize for multiple
parameters, for example TE and TM polarizations, one way is to define multiple
figures of merit, each leading to its own gradient functions:

$\delta\eta_{TE}=\frac{\omega\delta\epsilon}{2\
P}\int_{\Gamma}{u_{\bot}g_{TE}}ds,\text{ and }\\\
\delta\eta_{TM}=\frac{\omega\delta\epsilon}{2\
P}\int_{\Gamma}{u_{\bot}g_{TM}}ds.$ (25)

They can be combined in different ways, for example by using the inverse of
efficiency as a weighing function $(w_{TE,TM}=1/\eta_{TE,TM})$ so that the
optimization tends to balance well the performance of both polarizations,
$g=\left(w_{TM}g_{TM}+w_{TE}g_{TE}\right)/\left(w_{TM}+w_{TE}\right)$. One may
also define a single figure of merit as the sum of multiple figures of merit,
for example to take into account different wavelengths.

###  A.4 General figure of merit

The formulation here was based on a specific definition of the figure of
merit, expressed as the projection of the output field into the target field,
i.e., $\eta=\left|F\right|^{2}$ where $F=\frac{1}{4\
P_{t}}\int_{S}{\left(E\times H_{t}^{\ast}+E_{t}^{\ast}\times H\right)\cdot
n}da.$ This can be generalized as follows:

$F=\int_{V}f\left(E,H\right)dv,$ (26)

where for convenience, we normalize f so that it has units of inverse of
volume. The change in efficiency is then expressed as:

$\delta F=\int_{V}\left(\frac{\partial f}{\partial E}\cdot\delta
E+\frac{\partial f}{\partial H}\cdot\delta H\right)dv,\mathrm{\ and}\\\ \delta
F=\frac{1}{4\ P_{t}}\int_{S}\left(J_{a}\cdot\delta E-K_{a}\cdot\delta
H\right)da,$ (27)

where care must be taken to deal with derivatives with respect to complex
vectorial fields. The formulation then follows identically if we choose the
adjoint currents as:

$\left(J_{a},K_{a}\right)\delta S=\left(\frac{\partial f}{\partial
E},-\frac{\partial f}{\partial H}\right).$ (28)

##  Appendix B Experimental setup

A schematic of the experimental setup is shown above in Figure 12. The
linearly polarized output from a Santec TSL-570 tunable laser is collimated
and sent through a focusing lens to create a smaller sized beam at the focal
plane. A $\lambda/2$ plate is placed after the lens to control the orientation
of the linear polarization. The sample is located at the focal plane of the
lens, which is set to coincide with the image plane of the camera. A 10x
objective and a tube lens are used to image the sample onto the camera which
is used for alignment. A power meter is placed at different positions along
the beam path to input power and diffracted power.

Figure 12: Schematic of the experimental setup used to measure diffraction
efficiency.

Alignment: first, the beam is aligned onto the camera without the sample in
the system. Then, the position of the focusing lens is adjusted until the
focus of the beam coincides with the image plane of the camera. With a focal
length of 100 mm, the focusing lens produces a beam with a diameter of
approximately $100~{}\mu m$ at the focus for a wavelength of $\lambda$=1550
nm. This beam size is chosen to ensure that the beam is small enough to fit
entirely within one of the $250~{}\mu m$ x $250~{}\mu m$ grating patterns.
Next, the polarization is set by rotating the $\lambda/2$ plate to minimize
the power passing through a linear polarizer that is set to x-pol (horizontal)
or y-pol (vertical), relative to the table. Then, the sample is placed at the
focus of the beam in the image plane of the camera. With the sample positioned
so that the beam passes through the bare glass substrate, interference fringes
resulting from a tilt of the sample are visible on the camera. The sample is
aligned for normal incidence by adjusting the tilt angle until these
interference fringes disappear. Finally, the sample is positioned so that the
beam passes through one of the grating patterns and is rotated until the
height of the diffracted beam matches the height of the input beam.

Measurement Procedure: the input power is measured by placing the power meter
in the beam path right before the sample. The power is measured as the
tuneable laser is swept from 1480 nm – 1640 nm in steps of 0.1 nm. Then, the
beam is centered onto one of the grating patterns with the help of the imaging
system. The diffracted power is measured by placing the power meter in the
path of the 1st-order diffracted beam, positioned as close as possible to the
sample while ensuring that the 0-order beam is not detected by the sensor (2
cm). Again, the power is measured while the tuneable laser is swept over the
same wavelength range.

Figure 13: (a) The measured transmission through the bare substrate (blue
line) and the 1st-order diffraction of one of the gratings (orange line). (b)
Diagram of the Fabry-Perot model for both the bare substrate and one of the
grating patterns. The blue- and orange-colored arrows represent where the
transmission was measured for the bare substrate and grating, respectively.

Data Processing: in our simulations, the diffraction efficiency of the device
is defined relative to the power inside the substrate before interaction with
the metasurface. During measurement, however, we can only measure the power
before the substrate and use this as our reference. This results in many
oscillations in the measured transmission due to interference within the fused
silica substrate, as shown in Figure 13a. This can be understood by
recognizing that the substrate acts as a Fabry-Perot cavity between the
surfaces of the substrate and the 0th-order reflection from the metasurface. A
diagram is shown in Figure Figure 13b. The transmission of a Fabry-Perot
cavity for normal incidence is given by:

$T=\frac{T_{1}T_{2}}{1+R_{1}R_{2}-2\sqrt{R_{1}R_{2}}cos[2kn(\lambda)d]},$ (29)

where $d=500~{}\mu m$ is the thickness of the substrate, $n(\lambda)$ is the
wavelength-dependent refractive index of fused silica, $k=2\pi/\lambda$ is the
wavenumber, and $\lambda=1550~{}nm$. In this equation, the parameter $T_{2}$
represents the transmission of purely the second interface of the sample,
without any interference effects from the first interface. It is the
transmission of the sample relative to the power inside the substrate, exactly
how we define the diffraction efficiency. With knowledge of the transmission
and reflection coefficient for the first interface ($T_{1}$ and $R_{1}$), the
value of $T_{2}$ can be extracted from $T$, which we measure directly.
Consider the maximum and minimum values of the Fabry-Perot transmission
function. We find:

$T_{max}=\frac{T_{1}T_{2}}{\left(1-\sqrt{R_{1}R_{2}}\right)^{2}}\text{ and
}\\\ T_{min}=\frac{T_{1}T_{2}}{\left(1+\sqrt{R_{1}R_{2}}\right)^{2}}.$ (30)

The values of $T_{max}$ and $T_{min}$ can be extracted from the measurements.
$R_{1}$ and $T_{1}$ are the reflection and transmission coefficient for the
back surface of the substrate, with values given by the Fresnel equations. For
normal incidence from air to fused silica we have,

$R_{1}=\left|\frac{1-n\left(\lambda\right)}{1+n\left(\lambda\right)}\right|^{2}\text{
and }T_{1}=1-R_{1},$ (31)

where $n(\lambda)$ is given by the Sellmeier equation for fused silica. By
rearranging the equations for $T_{max}$ and $T_{min}$ we can find $T_{2}$ by
first solving for the unknown quantity $\sqrt{R_{1}R_{2}}$ using the equations
below:

$\sqrt{R_{1}R_{2}}=\frac{\left(1-\sqrt{\frac{T_{min}}{T_{max}}}\right)}{\left(1+\sqrt{\frac{T_{min}}{T_{max}}}\right)}\text{
and }\\\ T_{2}=\frac{T_{max}\left(1-\sqrt{R_{1}R_{2}}\right)^{2}}{T_{1}}.$
(32)

To verify this procedure for extracting $T_{2}$, we first consider the case of
a beam passing through the bare substrate for which $T_{2}$ can be readily
calculated according to the Fresnel equations. The results are shown below in
Figure 14.

Using the procedure outlined above for the case of transmission through the
bare substrate, we find an average value of $T_{2,measured}=0.9670079$ which
agrees very well with the value given by the Fresnel equations
$T_{2,Fresnel}=0.9670084$. The same processing is used for the 1st-order
transmission measurements of the meta-gratings to calculate the diffraction
efficiency. An example is shown below in Figure 15.

Figure 14: (a) The values of $T_{max}$ (green curve) and $T_{min}$ (red curve)
are interpolated from the measured transmission (blue curve) through bare
substrate. (b) The extracted value of $T_{2}$(blue curve) is compared with the
exact value given by the Fresnel equations (dashed black line). Figure 15: :
(a) The values of $T_{max}$ (green curve) and $T_{min}$ (red curve) are
interpolated from the measured 1st-order transmission (blue curve) of a meta-
grating. (b) The extracted value of $T_{2}$ (blue curve), defining the
measured diffraction efficiency.


# Are ultrafast inflows in AGN truly rare - or just much harder to see?

Ken Pounds, Physics and Astronomy, University of Leicester, Leicester, LE1
7RH, UK

(Accepted ; Submitted )

###### Abstract

Short-term variability and multiple velocity components in the powerful highly
ionized wind of the archetypal UFO PG1211+143 are indicative of inner disc
instabilities or short-lived accretion events. The detection of a high
velocity ($\sim$ 0.3c) inflow of highly ionized matter, located at 20 Rg,
offered the first direct observational support for the latter scenario, where
matter approaching at a high inclination to the black hole spin plane may
result in warping and tearing of the inner accretion disc, with subsequent
inter-ring collisions producing shocks, loss of rotational support and rapid
mass infall. Simultaneous soft x-ray spectra reveal a lower velocity ($\sim$
0.1c) inflow of less ionized matter, identified as ’upstream’ at 200 Rg, with
a line of sight through matter converging on the supermassive black hole. We
discuss here why ultrafast ionized winds are relatively common in luminous
Seyfert galaxies, while detection of the 0.3c inflow in PG1211+143 remains a
rare example.

###### keywords:

galaxies: active – galaxies: Seyfert: quasars: general – galaxies: individual:
PG1211+143 – X-ray: galaxies

††pagerange: Are ultrafast inflows in AGN truly rare - or just much harder to
see?–References††pubyear: 2005††pagerange: Are ultrafast inflows in AGN truly
rare - or just much harder to see?–References

##  1 Introduction

X-ray spectra from an XMM-Newton observation of the luminous Seyfert galaxy
PG1211+143 in 2001 first identified strongly blue-shifted absorption lines of
highly ionized gas, corresponding to a sub-relativistic outflow with a
velocity of 0.15$\pm$0.01$c$ (Pounds et al. 2003; Pounds & Page 2006). A
second example quickly followed with the detection of a powerful high velocity
wind in PDS 456 (Reeves et al. 2003), with archival data from XMM-Newton and
Suzaku subsequently showing that such ultra-fast, highly-ionized outflows
(UFOs) are relatively common in nearby, luminous AGN (Tombesi et al. 2010,
2011; Gofford et al. 2013). While those archival searches typically reported a
single velocity, in the few cases where an AGN was observed repeatedly, the
wind velocity was often different. Mkn509 is the best example from the XMM-
Newton data archive, with wind velocities of $\sim$0.173c, $\sim$0.139c and
$\sim$0.196c, separated by 5 years and 6 months respectively (Cappi et al.
2009).

Figure 1: rev2659 soft x-ray spectral fit with continuum and photoionized
emission features from the mean 2014 spectrum with the addition of: (top
panel) blue-shifted absorption corresponding to a prevailing outflow velocity
v$\sim$0.061c; (mid-panel) red-shifted absorption from inflowing matter at
v$\sim$0.096c; (lower panel) combined blue- and red-shifted absorption
components yielding an excellent overall spectral fit.

An extended 5-weeks XMM-Newton observation of PG1211+143 in 2014,covering 7
spacecraft orbits, showed a more complex velocity structure in the highly
ionized wind, with three primary (high column density) outflow velocities of
$v\sim 0.066c$, $v\sim 0.13c$ and $v\sim 0.19c$ (Pounds et al. 2016a), none
being consistent with the wind velocity observed in 2001. In a review of
continuum-driven or ’Eddington’ winds, King and Pounds (2015) furthermore
showed that the observability of an individual, short-lived wind ejection may
be of order months or less, as an expanding shell presents a lower absorbing
column to the nuclear x-ray source.

Considering the simultaneous observation of multiple expanding shells of
absorbing gas in the context of super-Eddington accretion, Pounds, Lobban and
Nixon (2017) noted that short-lived wind components might be a natural
consequence of the way in which matter accretes in an AGN, typically falling
from far outside the sphere of gravitational influence of the SMBH and with
essentially random orientation. The resulting accretion stream will in general
orbit in a plane misaligned to the spin of the central black hole (King &
Pringle 2006, 2007), with the inner disc subject to Lense-Thirring precession
around the spin vector, and orbits at smaller radii precessing sufficiently
fast to cause the tearing-away of independent rings of gas.

Computer simulations by Nixon et al. (2012) furthermore showed that - as each
torn-off ring precesses on its own timescale - two neighbouring rings will
eventually collide, with the shocked material losing rotational support and
falling inwards to a new radius defined by its residual angular momentum. A
substantial increase in the local accretion rate might then allow a previously
sub-Eddington flow to become briefly super-Eddington, with excess matter being
ejected as a wind at or above the local escape velocity, as predicted in the
classic paper on accretion disc theory by Shakura and Sunyaev (1973). For
typical AGN disc parameters, Nixon et al. (2012) show the expected ‘tearing
radius’ is of order a few hundred gravitational radii from the black hole -
with predicted (escape) wind speeds in the range observed (King and Pounds
2015; fig 4a).

The detection of red-shifted X-ray absorption spectra during part of the same
2014 XMM-Newton observation of PG1211+143, corresponding to a substantial
inflow of matter approaching the black hole at velocities up to $\sim$ 0.3c
(Pounds et al. 2018; hereafter Po18), offered the first observational support
for the above accretion scenario. Simultaneous hard and soft X-ray spectra
provided independent confirmation of the high inflow velocity, with the soft
X-ray absorption indicating less ionized (higher density) matter embedded in
the primary highly ionized flow. No further confirmed examples of an ultrafast
AGN inflow have been reported, at least to this author’s knowledge, which is
surprising in light of the above discussion. As noted in Po18, the probability
of detecting such an extreme velocity inflow, close to the black hole, will be
reduced both by a lower chance of a narrow converging flow being in the line
of sight to the continuum source, and the likely short duration of the infall.
Such considerations suggest a lower velocity, less highly ionized inflow, such
as that also seen in PG1211+143, could offer a better opportunity for
detection

In the present paper we report a new analysis of the soft x-ray data from the
Reflection Grating Spectrometer (RGS; den Herder et al.2001) in the 2014 XMM-
Newton observation of PG1211+143, separately identifying blue- and red-shifted
absorption components. In addition to the relatively weak high redshift
absorber (at $z=0.43$), we confirm a stronger soft x-ray absorber at a
redshift $z=0.189$, corresponding to an inflow velocity $\sim$0.1c, and (free-
fall) location at 200 Rg. We identify this absorption with matter lying
’upstream’ of the ultrafast, highly ionized inflow reported in Po18, in the
context of a line of sight along an accelerating flow converging on the black
hole.

We assume a cosmological redshift for PG1211+143 of $z=0.0809$ (Marziani et
al. 1996), with a black hole mass of $4\times 10^{7}$$\rm~{}M_{\odot}$ (Kaspi
et al. 2000) indicating the historical bolometric luminosity is close to
Eddington. All velocities are corrected for the relativistic Doppler effect.
The spectral analysis is based on the software package described in XSPEC
(Arnaud et al. 1996.

##  2 Soft x-ray absorption coincident with the hard x-ray inflow in
PG1211+143.

In order to quantify the strong red-shifted absorption unique to the second
orbit (rev2659) of the 2014 campaign, we start with the mean 2014 soft x-ray
spectral profile of PG1211+143 published in Pounds et al.2016b (Po16b), based
on the sum of RGS1 and RGS2 data from all 7 XMM-Newton orbits (for a combined
exposure of 1.27 Ms). Fitting the data over the 11-28 Å waveband, where the
data are of highest statistical quality, the continuum is modelled by power
law and black body components, both attenuated by the Galactic column. Strong
positive and negative residuals were matched in Po16b with photoionized
emission and absorption attributed to a circumstellar outflow. Details of the
multiplicative absorption and additive emission grids used are also given in
Po16b.

Here, we substitute the composite 7-orbit RGS data with that from rev2659, and
repeat the spectral fit, retaining continuum and photoionised emission, with
only the normalisation of each component left free. The blue-shifted
absorption reported in Po16b is initially excluded, though we expect it to be
present at some level throughout all 7 orbits. The fit is now quite poor, with
$\chi^{2}_{\nu}$ of 554/436 over the 11-28 Å waveband.

In a second spectral fit (Figure 1; top panel), blue-shifted absorption is
added, finding the primary (high column) outflow velocities of $\sim 0.06c$
and $\sim 0.13$ reported in Po16b, with column density and ionization
parameters listed in the Table 1. This addition recovered a much improved fit,
with $\chi^{2}_{\nu}$ of 485/432.

The data from rev2659 are next examined in the same way for red-shifted
absorption. All 3 components mentioned in Pounds et al. 2018 are found, with
column density, ionization parameter and redshifts listed in Table 1, together
with their individual significance. The fit statistic has $\chi^{2}_{\nu}$ of
525/427, weaker than for the blue-shifted absorption alone, but still highly
significant.

Finally, re-fitting the 7-orbit baseline emission model with both red- and
blue-shifted absorption yields an excellent spectral fit ($\chi^{2}_{\nu}$ of
455/421), with the combined flows provide a better match, in particular to the
complex of inner-shell Fe absorption lines (or UTA; eg Behar et al. 2001). For
visual clarity in Figure 1 we confine the blue- and red-shifted spectral plots
to the dominant flow component in each case.

Table 1: Parameters of 3 red- and 2 blue-shifted soft x-ray absorption components from a spectral fit to the RGS data from rev2659. Each absorber is further defined by its ionisation parameter $\xi$ (erg cm s-1) and column density NH (cm-2).  log$\xi$ | NH (1020) | redshift | vin | $\Delta\chi^{2}$
---|---|---|---|---
2.58$\pm$0.15 | 11$\pm$4 | 0.189$\pm$0.001 | 0.095$\pm$0.001 | 17/3
3.42$\pm$0.18 | 40$\pm$43 | 0.172$\pm$0.001 | 0.080$\pm$0.001 | 12/3
0.76$\pm$0.23 | 8$\pm$4 | 0.433$\pm$0.002 | 0.28$\pm$0.01 | 9/3
0.85$\pm$0.18 | 4.2$\pm$2.1 | -0.052$\pm$0.001 | -0.131$\pm$0.001 | 13/3
1.32$\pm$0.08 | 8$\pm$4 | -0.0168$\pm$0.0005 | -0.0611$\pm$0.0005 | 65/3

In the context of the present paper we speculate that previous analyses of the
soft x-ray spectra may well have missed a short-lived inflow in the presence
of a stronger, more persistent outflow. We acknowledge that was the case in
our analysis of the variable soft x-ray wind in PG1211+143 (Reeves et al.
2018), carried out before the highly ionised inflow detection at 0.3c was
made.

##  3 Discussion

The discovery of ultra-fast highly ionized outflows (UFOs) in early x-ray
spectra of the narrow line Seyfert galaxy PG1211+143 (Pounds et al. 2003) and
the luminous QSO PDS 456 (Reeves et al. 2003) opened up a new field of study
of AGN, well suited to the uniquely high throughput of x-ray spectrometers on
ESA’s XMM-Newton, launched in late 1999. King and Pounds (2003, 2015) noted
that such winds are a natural result of a high accretion ratio, with excess
matter being driven off by radiation pressure when the accretion rate exceeds
the local Eddington limit (Shakura and Sunyaev 1973).

While this picture provided a satisfactory explanation of most UFOs, where a
single detection yielded a unique wind velocity, the extended study of
PG1211+143 in 2014 found a more complex outflow profile, with velocities of
$\sim$0.06c, $\sim$0.13c and $\sim$0.18c detected in the stacked data set.
Such complexity was clearly inconsistent with a wind profile launched from a
flat axi-symmetric disc (SS73), suggesting some intrinsic disc instability or
rapidly variable accretion rate being introduced to the inner disc.

The detection of red-shifted X-ray absorption spectra during part of the same
2014 XMM-Newton observation of PG1211+143, where the extreme redshift
($\sim$0.48) in highly ionized matter indicated absorption in a substantial
inflow approaching the black hole at a velocity of $\sim$ 0.3c (Po18), offered
the first observational support for such rapidly rapidly variable accretion.
While the inflow observation was of high significance in only one of seven
individual observations, inflows of lower column density and smaller redshift
were detected in 5 of the other 6 EPIC pn observations, with redshifts ranging
from 0.20 to 0.36 (and line-of-sight inflow velocities from 0.1 – 0.23c).

As noted in the Introduction, periods of highly variable accretion to the
inner disc are a likely consequence of the way in which AGN accrete, where gas
initially falls towards the SMBH with essentially random orientation. Lense-
Thirring precession would cause misaligned orbits to precess around the black
hole spin vector, with the innermost ring(s) warping and potentially breaking
off. Collision between neighbouring rings, rotating at different rates, would
then shock, with loss of rotational support leading to matter falling freely
to a new radius defined by its residual angular momentum, where it may form a
new disc (Nixon et al. 2012).

Our new soft x-ray analysis, reported here, shows both blue- and red-shifted
photionized absorption strongly affecting the soft x-ray spectrum of
PG1211+143 during the 2014 XMM-Newton observations. The higher significance of
the outflow components is in part due to the 7-fold greater exposure for a
wind that persists throughout the XMM-Newton observation, with the outflow at
$\sim 0.061c$ particularly strong (see also Reeves et al.2018). In contrast,
red-shifted absorption is limited to rev2659 and much diluted in the overall
2014 data.

We identify the soft x-ray inflow component at $\sim 0.433c$ with higher
density matter surviving close to the primary, highly ionised inflow reported
in Po18. While providing independent support of the primary $\sim$0.3c inflow,
the accretion mass rate remains dominated by the highly ionised flow.

The strongest soft x-ray inflow component, at $\sim$0.095c, is the most
interesting and here we identify that slower moving and less highly ionized
absorber with matter ’upstream’ from the 0.3c inflow, with the correlation of
ionization parameter and velocity being also indicative of a converging and
accelerating inflow approaching the SMBH. In this impoortant respect we differ
from Reeves et al. (2018), where the absorption in rev2659 was discussed as a
random cloud or filament crossing the line of sight. Assuming, as before, a
free-fall velocity determined solely by the black hole’s gravity, the strong
soft x-ray inflow is located at 200 Rg.

Splitting the $\sim$100ks observation into 5 equal time intervals finds both
high and low ionisation inflows present at similar strength throughout
rev2659, showing the flow mass rate is maintained throughout the rev2659
observation. However, the high ionization, high column inflow is not detected
in either preceding (rev2652) or following (rev2661) spacecraft orbits.

Po18 estimate an accretion rate during rev2659 for the highly ionized inflow
of $\sim 10^{23}$ gm s-1, yielding a luminosity of $\sim 10^{43}$ erg s-1 for
a mass/energy conversion efficiency $\eta\sim$0.1, and lasting for $\sim
10^{5}$ s. In comparison, the mean bolometric luminosity of PG1211+143 is
$\sim 10^{45}$ erg s-1, with $\sim 10^{44}$ erg s-1 arising from the inner
disc region relevant here, suggesting a small number of infall events -
similar to that observed - would make a significant contribution to the inner
disc accretion rate, while remaing out of sight.

While the lower column density of a soft x-ray counterpart of the 0.3c inflow
is of less significance to the overall accretion rate, detection of co-moving,
higher density matter in the higher resolution RGS spectra provides valuable
confirmation of that extreme velocity. In addition, the strongest red-shifted
absorber provides an upstream measure of an accelerating stream of matter as
it converges and becomes more highly ionized on close approach to the SMBH in
PG1211+143. Most relevant in the context of the present paper, the slower and
broader upstream inflow is likely to be easier to detect. A search for such
soft x-ray inflows in new or archival XMM-Newton observations is highly
recommended.

##  4 Data availability

The data underlying this article are available in the XMM archive at
http://nxsa.esac.esa.int/nxsa-web.

##  Acknowledgements

XMM-Newton is a space science mission developed and operated by the European
Space Agency. We acknowledge in particular the excellent work of ESA staff in
Madrid in successfully planning and conducting the XMM-Newton observations. I
am particularly grateful for the collaboration of Andrew Lobban who provided
the initial data reduction on which this note is based, as well as being an
outstanding colleague in researching the x-ray properties of PG1211+143 over
many years.

## References

  * [1] Arnaud K.A. 1996, ASP Conf. Series, 101, 17
  * [2] Behar E., Sako M., Kahn S.M. 2001, ApJ, 563, 497
  * [3] Cappi M. et al. 2009, A&A, 504, 401
  * [4] den Herder J.W. et al. 2001, A&A, 365, L7
  * [5] Garcia J., Dauser T., Reynolds C.S., Kallman T.R., McClintock J.E., Wilms J., Eikmann W. 2013, ApJ, 768, 146
  * [6] George I.M., Fabian A.C. 1991, MNRAS, 249, 352
  * [7] Gofford J., Reeves J.N., Tombesi T., Braito V., Turner T.J., Miller L., Cappi M. 2013, MNRAS, 430, 60
  * [8] Kaspi S. Smith P.S., Netzer H., Maoz D., Jannuzi B.T., Giveon U. et al. 2000, ApJ, 533, 631
  * [9] King A.R., Pounds K.A. 2003, MNRAS, 345, 657
  * [10] King A.R., Pounds K.A. 2015, ARA&A, 53, 115
  * [11] King A.R. & Pringle J.E. 2006, MNRAS, 373, 90
  * [12] King A.R. & Pringle J.E. 2007, MNRAS, 377, 25
  * [13] Marziani P., Sulentic J.W., Dultzin-Hacyan D., Clavani M., Moles M. 1996, ApJS, 104, 37
  * [14] Murphy E.M., Lockman F.J., Laor A., Elvis M. 1996, ApJS, 105, 369
  * [15] Nandra K. and Pounds K.A. 1994, MNRAS, 268, 405
  * [16] Nixon C., King A.R., Price D., Frank J. 2012, ApJ, 757, 24
  * [17] Parker M.L. et al. 2015, MNRAS, 447, 72
  * [18] Pounds K.A., Reeves J.N., King A.R., Page K.L., O’Brien P.T., Turner M.J.L. 2003, MNRAS, 345, 705
  * [19] Pounds K.A., Reeves J.N., Page K.L., O’Brien P.T. 2004, Astrophys. J, 616, 696
  * [20] Pounds K.A., Page K.L. 2006, MNRAS, 360, 1123
  * [21] Pounds K.A., Lobban A., Reeves J.N., Vaughan S. 2016(a), MNRAS, 457, 2951 (Po16a)
  * [22] Pounds K.A., Lobban A., Reeves J.N., Vaughan S., Costa M. 2016(b), MNRAS, 459, 4389 (Po16b)
  * [23] Pounds K.A., Lobban A.P., Nixon C. 2017, Astron. Nachr., 338, 249
  * [24] Pounds K.A., Nixon C., Lobban A.P., King A.R. 2018, MNRAS, 481, 1832 (Po18)
  * [25] Reeves J.N., O’Brien P.T., Ward M.J. 2003, ApJ, 593, 65
  * [26] Reeves J.N., Lobban A.P., Pounds K.A. 2018, ApJ, 854, 1538
  * [27] Shakura N.I. & Sunyaev R.A. 1973, A&A, 24, 337 (SS73)
  * [28] Shull J.M and Van Steenberg M. 1982, Astrophys J. Supp., 48, 95
  * [29] Strueder L. et al. 2001, A&A, 365, L18
  * [30] Titarchuk L. 1994, Astrophys.J., 434, 570
  * [31] Tombesi F., Cappi M., Reeves J.N., Palumbo G.C., Yaqoob T., Braito V., Dadina M. 2010, ApJ, 742, 44
  * [32] Tombesi F., Cappi M., Reeves J.N., Palumbo G.C., Braito V., Dadina M. 2011, A&A, 521, A57
  * [33] Vaughan S. and Fabian A.C. 2004, MNRAS, 348, 1415