# Inverse-designed metaphotonics for hypersensitive detection

Maxim S. Elizarov1, Yuri S. Kivshar2,3, and A. Fratalocchi1
<EMAIL_ADDRESS>www.primalight.org 1PRIMALIGHT, Faculty of
Electrical Engineering; Applied Mathematics and Computational Science, KAUST,
Thuwal 23955-6900, Saudi Arabia  2Australian National University, Canberra ACT
2601, Australia  3ITMO University, St. Petersburg 197101, Russia

###### Abstract

Controlling the flow of broadband electromagnetic energy at the nanoscale
remains a critical challenge in optoelectronics. Surface plasmon polaritons
(or plasmons) provide subwavelength localization of light, but are affected by
significant losses. On the contrary, dielectrics lack a sufficiently robust
response in the visible to trap photons similar to metallic structures.
Overcoming these limitations appears elusive, as it implies devising a path to
circumvent causality in the quantum-mechanical form of matter. Here we
demonstrate that addressing this problem is possible if we employ a novel
approach based on suitably deformed reflective metaphotonic structures. The
complex geometrical shape engineered in these reflectors emulates
nondispersive index responses, which can be inverse-designed following
arbitrary form factors. We discuss the realization of essential components
such as resonators with an ultra-high refractive index of $n=100$ in diverse
profiles. These structures support localization of light in the form of bound
states in the continuum (BIC), fully localized in air, in a platform in which
all refractive index regions are physically accessible. We discuss our
approach to sensing applications, designing a class of sensors where the
analyte directly contacts areas of ultra-high refractive index. Leveraging
this feature, we report differential sensitivities up to $350$ nm/RIU in
structures with footprints of approximately one micron. These performances are
two times better than the closest competitor with a similar form factor.
Inversely designed reflective metaphotonics offers a flexible technology for
controlling broadband light, supporting optoelectronics’ integration with
large bandwidths in circuitry with miniaturized footprints.

## Introduction

In applications that exploit the propagation of electromagnetic beams, it is
crucial to control broadband fields in small volumes of matter. At GHz
frequencies, this reduces chip footprint and increases device speed. At
optical wavelengths, strongly localized fields provide the key to many
significant effects, including low-threshold lasing, efficient nonlinear
harmonic generations, high-resolution imaging, highly sensitive detection,
low-power communications, enhanced security, and complex wavefront engineering
lee2019ultrahigh ; koshelev2019nonradiating ; Koshelev288 ; wu2021bound ;
huttenhofer2021anapole ; wang_maier_crescent_met_2021 ; Nasir_Zayats_2021 ;
smirnova2020nonlinear ; chen2021topologically . However, despite impressive
recent progress in this field, photonics integration still lags behind
electronics, which employs now pioneering devices with characteristic scales
of 10 nm or less el0 ; el1 .
To trap light at subwavelength scales, one employs a traditional approach
based on collective excitations of electrons and electromagnetic waves in
metallic resonators supporting hybrid modes described by surface plasmon
polaritons. Currently, metals provide the highest degree of light confinement
via plasmonic resonance maier . However, substantial losses at optical
wavelengths create a significant challenge to exploiting these materials for
advanced applications franano ; dar0 ; koshelev2019nonradiating .
The recently emerged field of dielectric resonant metaphotonics provides an
alternative physical mechanism of light localization via low-order dipole and
multipole Mie resonances that support highly efficient nanoresonators and all-
dielectric metasurfaces koshelev2020dielectric . This so-called “Mie-tronics”
approach mie_won  successfully demonstrated the suppression of radiative
losses of individual dielectric resonators while engineering unique optical
modes with high-quality factors facilitating nonlinear effects at the
nanoscale Koshelev288 .
However, the refractive index $n$ attainable in the visible frequency range is
relatively small (usually $n<5$), especially compared to the effective
refractive index of metals or dielectrics at longer wavelengths. The currently
available materials do not allow realizing a large manifold of concepts
initially developed for plasmonics and microwave metamaterials, due to either
large metallic losses or unavailability of materials with sufficiently large
refractive index response at optical wavelengths. This issue presently
confines conceptual demonstrations of many metaphotonics ideas at microwave
frequencies shen2013planar ; bonache2016controlling ; yermakov2018experimental
.
Thus, a novel approach to break this limit is highly desirable for
implementing many revolutionary ideas with a novel generation of low-loss
optical materials and for applications of deep subwavelength optics. Here we
introduce a novel platform to implement many of those ideas in lossless
reflective materials that can help devise many applications. Our approach
exploits an inverse design of deformed reflecting structures with user-defined
nondispersive refractive indices and spatial form profiles. We discuss basic
optical circuitry such as resonators with an ultra-high refractive index of
$n=100$ and design a new class of hypersensitive sensors that showed a record
value of sensitivity. We believe our approach will provide novel means of
strong energy confinement, opening the door to integrated optoelectronic
devices with simultaneous control of electrons and photons of large bandwidths
in channels of comparable size.

## Results

### Inverse-designed materials via geometrical deformations

The crucial challenge in controlling light localization at the nanoscale is
the quantum mechanical structure of matter combined with reciprocity,
representing the intuitive condition that a physical system cannot anticipate
the future doi:10.1002/9783527618156 . Quantum mechanics dictates that the
susceptibility response $\chi(\omega)=n(\omega)^{2}-1$ of any accessible
material is a frequency $\omega$ dependent function combining a series of
Lorentz-type resonances (see Chapter 3 of 10.5555/1817101 ):

$\chi(\omega)=\sum_{m}\left(\frac{a_{m}}{\omega_{m}-\omega}+\frac{a_{m}}{\omega_{m}^{*}+\omega}\right),$
(1)

with complex oscillator frequencies $\omega_{m}$ and amplitudes $a_{m}$. The
causality principle manifests in Eq. (1) as a rigid relationship, expressed by
Kramers-Kronig (KK) relations doi:10.1002/9783527618156 , between the real and
imaginary part of $\chi(\omega)$. Due to the direct relationship between
$\chi(\omega)$ and $n(\omega)$, KK relations state the impossibility to
independently control the real part of $n(\omega)$, representing the effective
refractive index of the material, from the material losses, defined by the
imaginary part of the refractive index $n(\omega)$.
Metallic media possess atomic resonances $\omega_{m}$ at visible wavelengths
maier , and acquire a strong light localization ability at visible
frequencies. However, KK relations constraint high-localization frequencies to
points of high absorption, blending the effects of energy localization and
losses irreversibly. Dielectric materials, conversely, possess electronic
resonances $\omega_{m}$ in the deep-blue or ultra-violet frequency range and
exhibit transparency over large optical bandwidths. However, the accessible
refractive index of these materials in the lossless window is inevitably
modest, compromising the resulting localization power, especially if compared
to plasmonic media.
Figure 1 shows an alternative material platform that could overcome the issue
mentioned above. The starting configuration is a universal basic structure
composed by a reflective substrate lying in the $(x,z)$ plane, with a user-
defined, semi-infinite material defined on top (Fig. 1a). The medium possesses
a nondispersive and lossless anisotropic refractive index
$\mathbf{n}(\mathbf{r})=\sqrt{\bm{\epsilon}(\mathbf{r})\cdot\bm{\mu}(\mathbf{r})}$,
with $\bm{\epsilon}(\mathbf{r})$ and $\bm{\mu}(\mathbf{r})$ dielectric
permittivity and magnetic permeability tensors, respectively, and
$\mathbf{r}=(x,y,z)$ position coordinates.
While a material defined as such does not directly exist in nature, we here
show that it is possible to engineer this structure by an inverse application
of transformation optics Pendry549 ; tra2 ; McCall_2018 ; So_2019 . Inverse
design is a powerful method for creating a material with desired properties
pilozzi2018machine ; so2020deep  while transformation optics establishes a
relationship between coordinate transformations and geometric materials
equivalent to light propagation. Following this correspondence, we wish to
implement the user-defined medium $\mathbf{n}(\mathbf{r})$ from a geometrical
deformation of coordinates $\mathbf{r}^{\prime}=\bm{\Omega}(\mathbf{r})$ with:

$\mathbf{r}^{\prime}=\begin{bmatrix}x^{\prime}\\\ y^{\prime}\\\
z^{\prime}\end{bmatrix}=\begin{bmatrix}x^{\prime}(x,y,z)\\\
y^{\prime}(x,y,z)\\\ z^{\prime}(x,y,z)\end{bmatrix}=\bm{\Omega}(\mathbf{r}),$
(2)

with $\mathbf{r}^{\prime}$ transformed coordinates. In the transformed space
$\mathbf{r}^{\prime}$, Maxwell equations remain identically the same with the
introduction of a new material with dielectric permittivity
$\bm{\epsilon}^{\prime}(\mathbf{r}^{\prime})=\frac{\nabla\bm{\Omega}\cdot\bm{\epsilon}\cdot\nabla\bm{\Omega}^{\dagger}}{|\nabla\bm{\Omega}^{\dagger}|}$
and magnetic permeability
$\bm{\mu}^{\prime}(\mathbf{r}^{\prime})=\frac{\nabla\bm{\Omega}\cdot\bm{\mu}\cdot\nabla\bm{\Omega}^{\dagger}}{|\nabla\bm{\Omega}^{\dagger}|}$
with
$\nabla\bm{\Omega}=\frac{\partial\mathbf{r}^{\prime}}{\partial\mathbf{r}}$ the
Jacobian matrix tra2 . Light cannot differentiate between the materials
$\mathbf{n}(\mathbf{r})$ and
$\mathbf{n}^{\prime}(\mathbf{r}^{\prime})=\sqrt{\bm{\epsilon}^{\prime}(\mathbf{r}^{\prime})\cdot\bm{\mu}^{\prime}(\mathbf{r}^{\prime})}$,
experiencing an identical evolution in the spaces $\mathbf{r}$ and
$\mathbf{r}^{\prime}$. For any response pair $n(\mathbf{r})$ and
$n^{\prime}(\mathbf{r}^{\prime})$ created via $\bm{\Omega}(\mathbf{r})$,
transformation optics guarantees the causality of the resulting material in
both spaces $\mathbf{r}$ and $\mathbf{r}^{\prime}$ (see Methods for a detailed
demonstration).
The main idea is to inverse design the geometrical deformation (2) so that the
transformed refractive index is that of a vacuum, with
$\bm{n}^{\prime}(\mathbf{r}^{\prime})=\mathbb{1}$:

$\nabla\bm{\Omega}\cdot\bm{\epsilon}\cdot\nabla\bm{\Omega}^{\dagger}\cdot\nabla\bm{\Omega}\cdot\bm{\mu}\cdot\nabla\bm{\Omega}^{\dagger}=|\nabla\bm{\Omega}^{\dagger}|^{2}\cdot\mathbb{1}.$
(3)

When this special condition occurs, the transformed medium
$\bm{n}^{\prime}(\mathbf{r}^{\prime})$ appears as a deformed reflector
immersed in a vacuum, and the deformed surface defined from (2) acquires the
ability to emulate the material with the user-defined refractive index
$\bm{n}$ (Fig. 1b-c).
Equation (3) comprises a system of nine nonlinear partial differential
equations, to be solved for the unknown transformation
$\bm{\Omega}(\mathbf{r})$ for a given desired material response in
$\bm{\epsilon}(\mathbf{r})$ and $\bm{\mu}(\mathbf{r})$. As a boundary
condition, we assume that for $\mathbf{r}\rightarrow\infty$ the transformation
(2) tends to the identity $\bm{\Omega}(\mathbf{r})=\mathbf{r}$, with
$\nabla\bm{\Omega}=\mathbb{1}$. This condition ensures that Eq. (3) represents
physical materials localized in a finite area of space, while relaxing to a
vacuum
$\bm{\epsilon}(\mathbf{r})=\bm{\mu}(\mathbf{r})=\mathbf{n}(\mathbf{r})=\mathbb{1}$
in the far-field for $\mathbf{r}\rightarrow\infty$.
In this article, we begin to study the solution of Eqs. (3) for scalar
dielectric or magnetic materials defined by a single quantity
$\epsilon(\mathbf{r})\equiv\epsilon_{zz}(\mathbf{r})$ and
$\mu(\mathbf{r})\equiv\mu_{zz}(\mathbf{r})$, assumed without loss of
generality on the $z$ axis:

$\displaystyle\bm{\epsilon}=\begin{bmatrix}\mathbb{1}&0\\\
0&\epsilon(\mathbf{r})\end{bmatrix},$
$\displaystyle\bm{\mu}=\begin{bmatrix}\mathbb{1}&0\\\
0&\mu(\mathbf{r})\end{bmatrix}.$ (4)

The resulting refractive index $\mathbf{n}=\sqrt{\bm{\epsilon}\cdot\bm{\mu}}$
acquires a traditional expression depending on a single scalar quantity
$n(\mathbf{r})=\sqrt{\epsilon(\mathbf{r})\mu(\mathbf{r})}$. Supplementary Note
I shows that a necessary condition for the solution of Eq. (3) is that the
deformation (2) is $z-$invariant with $\nabla=[\nabla_{\perp},0]$ and
$\nabla_{\perp}=[\frac{\partial}{\partial x},\frac{\partial}{\partial y}]$. In
this condition, Eqs. (3) reduce to:

$\begin{bmatrix}\left(\nabla_{\perp}\Omega_{\perp}\cdot\nabla_{\perp}\Omega^{\dagger}_{\perp}\right)^{2}&0\\\
0&n^{2}(\mathbf{r})\end{bmatrix}=|\nabla_{\perp}\Omega^{\dagger}_{\perp}|^{2}\cdot\mathbb{1},$
(5)

with $\bm{\Omega}_{\perp}=[x^{\prime}(x,y),y^{\prime}(x,y)]$. Supplementary
Note II shows that the solution to Eqs. (5) is an inverse conformal mapping
schinzinger1991conformal , with $\bm{\Omega}=[\rm I\\!R(\Omega),\rm
I\\!I(\Omega)]$ defined from the real $\rm I\\!R$ and imaginary $\rm I\\!I$
part of a single analytic function $\Omega(u)$ of complex coordinate $u=x+iy$
satisfying the scalar equation
$n^{2}(\mathbf{r})=\left|\frac{d\Omega(u)}{du}\right|^{2}$. We solve this
problem with a pseudospectral approach fornberg_1996  based on suitably
defined rational Chebyshev polynomials doi:10.1002/nme.392  in the complex
domain. We expand the unknown deformation $\Omega(u)$ as follows:

$\Omega\left(\frac{u}{a}\right)=\frac{u}{a}+\beta_{0}+\sum_{m=1}^{\infty}\beta_{m}T_{m}\left(\frac{u-ia}{u+ia}\right),$
(6)

with $T_{m}(v)$ the Chebyshev polynomial of order $m$, $a$ an arbitrary
spatial scaling constant, and $\beta_{m}$ unknown coefficients. Chebyshev
polynomials $T_{m}(v)$ are a complete basis that represents any complex
function on the unit circle $-1\leq|v|\leq 1$ fornberg_1996 . In (6), we use a
a Cayley transform $v=\frac{u-ia}{u+ia}$ kreyszig11  that maps the unit circle
to the positive semi-infinite space $y\geq 0$ and $-\infty\leq x\leq\infty$,
providing a rational series expansion that represents any analytic function
$\Omega(u)$ in the complex domain. Additionally, we impose the condition
$\beta_{0}=-\sum_{m=1}^{\infty}\beta_{m}$ to reduce Eq. (6) in the far-field
to the identity transformation $\Omega(u/a)=u/a$, as the user can verify by
direct substitution and as requested for the inverse solution of Eq. (3).
By expressing the derivative $\frac{dT_{m}(u)}{du}=m\cdot U_{m-1}(u)$ with
Chebyshev polynomials of second kind $U_{m}$, we obtain the nonlinear equation
for the inverse design of the transformation $\Omega(u)$:

$n^{2}\left(\frac{\mathbf{r}}{a}\right)=\left|1+\frac{2ia}{(u+ia)^{2}}\cdot\sum_{m=0}^{\infty}\beta_{m+1}(m+1)\cdot
U_{m}\left(\frac{u-ia}{u+ia}\right)\right|^{2}.$ (7)

In the scalar limit of Eqs. (4), the inverse design problem is equivalent to
finding the set of coefficients $\beta_{m}$ that satisfy Eqs. (7) for a user-
defined refractive index $n(\mathbf{r})$ given at the input. We solve Eq. (7)
by statistical learning via nonlinear regression hastie_09_elements-
of.statistical-learning . Given a user-defined distribution of refractive
index $n(\mathbf{r})$, we create a training dataset composed of a discrete
number $n=0,1,...$ of refractive index values $n(\mathbf{r}_{n})$ sampled on a
Chebyshev grid $\mathbf{r}_{n}$ fornberg_1996 . Each dataset input-output
couple $[n(\mathbf{r}_{n}),\mathbf{r}_{n}]$, when substituted in Eq. (7)
originates a nonlinear equation for the regression coefficients $\beta_{m+1}$.
These are calculated via nonlinear least-square, by using a trust-region
convex minimization routine ceres-solver . After solving for the coefficients
$\beta_{m+1}$, Eq. (6) predicts the deformation of coordinates $\Omega(u)$ for
any point $\mathbf{r}^{\prime}$ of the space and completes the solution of the
problem.
In the family of materials designed via (7), the refractive index
$n(\mathbf{r})$ is solely defined from the spatial curvature of $\Omega(u)$
arising from $\left|\frac{d\Omega(u)}{du}\right|^{2}$. As such,
$n(\mathbf{r})$ does not show any theoretical limit between zero and infinity,
with accessible values limited to implementing the required curvature arising
from the solution of (7).

### Nanoresonators with ultra-high refractive index

At visible frequencies, high refractive index available materials with $n\leq
5$ are typically III-IV semiconductors including Ge, Si, GaAs, GaSb. In these
media, it is traditionally challenging to exploit their highest values of the
refractive index because of high losses kasap2017springer ; baranov2017all ,
and applications typically rely on $3\leq n\leq 4$
miroshnichenko2015nonradiating ; Koshelev288 ; Liu2016AOM ; Schmid_2017 ;
YANG20171 ; PhysRevLett.119.243901 ; fran ; melik2021fano ; yurirev ;
zhizhchenko2019single . To the best of the authors’ knowledge, the reported
highest refractive index in the visible is $n=5.1$ for a block copolymer self-
assembly metasurface kim2016highly .
We here illustrate how it is possible via inverse design to engineer materials
with ultra-high refractive index $(n\approx 100)$ and user-defined profile. We
set a Gaussian distribution of refractive index $n(\mathbf{r})$ with amplitude
$n(\mathbf{r}=0)=100$ and planar width $x/a=0.1$. Figure 2a shows the
corresponding training dataset (gray circle markers). In the solution of (6),
we evaluate the minimum number $M$ of coefficients $\beta_{m+1}$ required to
inverse solve the problem by increasing $M$ and computing at each stage the
prediction mean square error $MSE=||n(\mathbf{r})-n_{pred}(\mathbf{r})||$,
with $n_{pred}(\mathbf{r})$ the refractive index predicted from the
transformation $\Omega(u)$. Few modes with $M=30$ yield a solution with
accuracy $10^{-7}$ (Fig. 2b ), and prediction values along the axis $x$
matching the required Gaussian profile (Fig. 2a solid red line). Figure 2c
illustrates the refractive index distribution in the design space
$\mathbf{r}$. The index relaxes exponentially in the $y-$axis, providing a
diffusive-like nanoresonator profile with an ultra-high refractive index
region defined in a narrow length along $y$.
Figure 2d (solid gray area) shows the reflective deformation in the accessible
space $\mathbf{r}^{\prime}$ that emulates the material in Fig. 2c. The
deformation corresponds to the mapping
$\mathbf{r^{\prime}}=\bm{\Omega}(\mathbf{r})$ for $y=0$ and
$-\infty<x<\infty$. The solid blue lines in Fig. 2d report the coordinate
deformation in the remaining area of space for $y\neq 0$. A single coordinate
transformation $\Omega(u)$ defines an entire family of materials with
different distributions of refractive indices. If we consider a coordinate
line at varying $x$ and constant $y_{0}$ in Fig. 2c, and set $y_{0}$ as a new
origin of coordinate along $y$ via $y\rightarrow y-y_{0}$, this deformation
generates a new reflective material in the accessible space $\mathbf{r}$ with
the index distribution defined by $n(\mathbf{r})$ evaluated at the new origin
$y=0$ and varying $x$. The corresponding refractive index maintains the same
spatial profile along $x$, but acquires a different amplitude due to the
relaxation of the index values towards infinity (Fig. 2c). Figures 2e-h
illustrate this possibility, showing the engineering of different reflective
deformations yielding Gaussian index profiles with maxima of $n=20$ and
$n=70$.
Supplementary Fig. 1 illustrate examples of ultra-high index resonators
realized with different form factors. Panels a-c show the case of a flat-top
resonator with a super-Gaussian spatial index profile. Supplementary Fig. 1d-f
shows the example of two adjacent resonators with different refractive
indices, demonstrating the flexibility of this technique in designing complex
index modulations that are not possible to observe at visible frequencies with
conventional materials.

### Resonator modes

Resonators designed by complex geometrical deformations support singular
energy localization states appearing as Bound States in the Continuum (BIC)
BogdanovFratalocchiKivshar . We studied singular resonances from the user-
defined space $\mathbf{r}$, by decomposing the electric field
$\mathbf{E}=\mathbf{E}_{t}+\mathbf{E}_{y}$ into a transverse $\mathbf{E}_{t}$
component, lying on the $(x,z)$ plane, and a normal contribution
$\mathbf{E}_{y}$ along $y$. We then define the following eigenvalue problem
for the transverse field:

$\frac{\partial}{\partial
y}\mathbf{E}_{t}(\mathbf{r})=\gamma(y)\cdot\mathbf{E}_{t}(\mathbf{r}),$ (8)

defining the spectral decomposition of the operator $\frac{\partial}{\partial
y}$ along the normal $y$ axis. Of particular interest are trapped modes
arising from the singularities in the eigenvalue amplitude
$|\gamma(0)|\rightarrow\infty$ arising at $y=0$. For the corresponding
eigenvector $\mathbf{E}_{t}(\mathbf{r})$ to be bounded and physical, the
contribution $\frac{\partial}{\partial y}\mathbf{E}_{t}$ is required to stay
finite. This condition implies that $\mathbf{E}_{t}\rightarrow 0$ at $y=0$,
satisfying the boundary condition $\mathbf{E}_{t}=0$ imposed by a perfect
reflector over the entire $(x,z)$ plane for transverse-electric (TE) polarized
modes. Each singularity of Eq. (8) therefore defines a TE-polarized,
nonradiating state that exists without any impinging source, by autonomously
satisfying the boundary conditions on the reflector plane at $y=0$. The
corresponding eigenvector $\mathbf{E}_{t}(\mathbf{r})$ furnishes the spatial
profile of the resulting TE-polarized BIC. As the singularity originating the
BIC arises at $y=0$, the electromagnetic field composing the BIC mode
typically localizes in the proximity of the reflective plane at $y=0$. In the
accessible space $\mathbf{r}^{\prime}$, such energy localizations appear in
the vacuum-gaps areas created by geometric mapping $\bm{\Omega}(\mathbf{r})$
(Fig. 1b). The case for TM polarization follows the same analysis, with the
only difference of decomposing the magnetic field
$\mathbf{H}=\mathbf{H}_{t}+\mathbf{H}_{y}$ into transverse $\mathbf{H}_{t}$
and normal $\mathbf{H}_{y}$ contributions, and formulating Eq. (8) for
$\mathbf{H}_{t}$. The boundary condition of TM-polarized localizations,
$\mathbf{H}_{t}=0$ at $y=0$, implies that the energy maxima of these modes
localize on the reflector surface. TE localizations, conversely, possess an
electric field that completely localizes in the vacuum. As the latter
configuration is more advantageous for the applications discussed in this
work, we focus on TE modes while deferring the analysis of TM-polarized BIC to
future work.
We illustrate these results by solving the eigenvalue Eq. (8) with an approach
that generalizes the work of PAGNEUX20101834 . The idea is to project the
eigenvalue equation over a basis of scattering modes computed at
$\mathbf{r}\rightarrow\infty$, and resolve the resulting discrete equations by
a spectral method. Supplementary section III discusses implementation details
of this approach. Figure 3 illustrates calculation examples of BIC supported
by the ultra-high refractive index nanoresonator designed in Fig. 2a. Figure
3a plots the corresponding eigenvalue amplitude $|\gamma(0)|$ at $y=0$ for a
varying frequency measured in wavevector units $k=2\pi\omega/c$. As the
resonator structure is symmetric in the $(x,z)$ plane, we perform computations
for even (solid blue lines) and odd (solid orange lines) BIC modes separately.
The eigenvalue $\gamma$ shows the formation of clear singular states with
amplitudes reaching up to $10^{5}$ at $k=0.3429$. Each of these singularities
defines a BIC state with a different energy profile. Figure 3b reports the
spatial energy distribution in the accessible space, while panel c shows the
corresponding profile in the design space. In the design space, odd and even
TE polarized BIC modes confine electromagnetic energy inside the ultra-high
refractive region in proximity to the reflective surface. In the accessible
space, these modes manifest as narrow light localization in the vacuum gaps
defined by the complex reflector spatial profile.
An advantage of this form of energy localization is the possibility to access
the entire spatial distribution of high refractive index, as observed from the
mapping between the design and accessible spaces in Fig. 3b-c. This property
results from the user-defined refractive index region created by the surface
profile of the deformed reflector, which provides trapping areas for singular
states with a high degree of energy localization. In the next Section, we
leverage this feature to design a new family of compact hypersensitive
devices.

### Sensing application

Figure 4a summarizes the state-of-the-art in refractive index (RI) sensing in
sensitivity, measured in nm/RIU, vs. sensor footprint, reported in units of
microns yalcin2006optical ; kim2008integrated ; girault2015integrated ;
schmidt2014improving ; liu2016simultaneous ; wei2016direct ;
goykhman2010ultrathin ; tang2016packaged ; hu2008planar ; wang2013suspended ;
kim2019chip ; grist2013silicon ; hamed2013photonic ; di2009chemical ;
garcia2010label . Figure 4b order these results in terms of a single FOM
measuring the differential sensitivity per sensor pathlength, obtained as the
ratio between the differential sensitivity and the sensor footprint. Sensors
created with single resonator structures have typically lower performances
than arrays of periodic structures such as, e.g., photonic crystals (PhC),
which can reach sensitivities up to 1500 nm/RIU. Addressing sensor
miniaturization without compromising performances is challenging: single
resonator structures cannot benefit from nearest-neighbor constructive
interactions and leverage only the quality-$Q$ factor of the resonant modes
supported by the refractive index of the resonator. While in microdisk
resonators with footprints of few microns, the $Q$ factor can reach $Q=12000$
wang2013suspended , this value decreases to $Q=680$ in the presence of
analytes, compromising the overall performances. A detrimental issue is also
the spatial distribution of resonator modes, which usually localize their
energy in high refractive index regions. These areas are the portion of
resonator space inaccessible from the environment and that analytes cannot
probe. As a result of these shortcomings, sensitivities up to $300$ nm/RIU
requires single resonators with large footprints of tens of microns.
The inverse design platform developed in this work can help to address these
issues by designing compact resonators with ultra-high refractive index
profiles and exploiting the high-quality BIC resonances supported by these
structures. Figure 5a-b illustrates the main design idea. We consider a single
resonator in the design space possessing an ultra-high refractive index of
$n=100$ and Gaussian form factor, as shown in Fig. 2c, and place an analyte on
top (Fig. 5a blue area). In this configuration, the analyte permeates the
entire refractive index region, including the areas of $n=100$. In the
accessible space (Fig. 5b), the structure appears as a deformed reflector with
a conformal layer of material that fills up the cap space created by the
modulated surface.
Figure 5c presents the density of states (DOS) of the resonator. We perform
this calculation with Finite-Difference Time-Domain (FDTD) simulations. We
extract quality $Q$ factors via a robust Prony method obtained by augmenting
the original Prony scheme with statistical regression (see details in
Methods). The structure supports diverse BIC states with quality factors up to
$Q=5200$, which we choose for our analysis. By setting the unit length
$a=1~{}\mu$m, the BIC with the highest quality factor occurs at the resonant
wavelength is $\lambda_{BIC}=307~{}$nm, in the ultra-violet (UV) region.
Figure 5d illustrates the spatial energy distribution of the BIC in the
accessible space. The BIC provides robust electromagnetic energy localization
inside the deformed reflector, showing a complex distribution of multiple
energy hot spots. In this sensor configuration, the analyte fills up the
entire gap area (Fig. 5b) and experiences a strong light-matter interaction by
probing the high energy spots supported by the BIC state.
Figure 6 analyzes this mechanism in detail by computing the differential
sensitivity from the wavelength shift of the BIC resonance. We computed these
results from FDTD simulation by changing the refractive index of the analyte
with steps $\Delta n=5\cdot 10^{-3}$ RIU. Figure 6a reports the variation of
the density of states for analytes with refractive indices between $n=1.005$
and $n=1.035$. While we do not observe any appreciable variation of the
$Q-$factor, the system shows a perfect linear shift of the BIC resonant
frequency (Fig. 6b solid red line). The resulting differential sensitivity
stays between $S=350$ nm/RIU and $S=310$ nm/RIU, with an average value of
$S=330$ nm/RIU (Fig. 6b solid blue line). This design reports a record value
of sensitivity per unit pathlength (Fig. 4b red circle), providing a single
resonant structure with similar performances of periodic assembly of
resonators with footprints of one order magnitude larger (Fig. 4a).

## Conclusion

We have presented an inverse-design framework for implementing structures with
user-defined nondispersive refractive indices and spatial form profiles in
suitably deformed reflectors. We have discussed the basic optical circuitry,
such as resonators with an ultra-high refractive index of $n=100$ in diverse
user-defined shapes. Leveraging the high refractive index response of these
systems, we have designed a new class of hypersensitive sensors that
demonstrate record values of differential sensitivity per unit pathlength
$S_{pl}=290$ nm$\cdot\mu$m${}^{-1}\cdot\mathrm{RIU}^{-1}$.
Notably, such inverse metaphotonic reflectors can overcome two challenging
problems. They allow (i) to design structures with similar trapping
performances of plasmonic structures but with no losses and (ii) to implement
dielectric materials with user-defined refractive indices that go beyond the
values available in physical materials. While these inverse-designed
reflectors do not eliminate all issues, they replace physical limitations
arising from the points (i) and (ii) with a more straightforward problem of
engineering a surface with a user-implemented profile. Some existing
nanofabrication techniques are already relevant to address these issues. While
an exhaustive analysis is beyond the scope of this conceptual work, here we
can mention two-photon lithography, which can already reach resolutions in
below the micron range geng2019ultrafast  and is successively covered by a
reflective layer, or powder bed fusion, which is essentially a version of
three-dimensional printing for metallic structures elder2020nanomaterial .

## Acknowledgements

The authors acknowledge the use of resources of the Supercomputing Laboratory
at KAUST. Y.K. acknowledges a support from the Australian Research Council
(grant DP210101292), Russian Science Foundation (grant 21-72-30018), and the
US Army International Office (grant FA520921P0034).

## Methods

Causality. We here prove that a material arising from the solution of Eq. (3)
is physical and obeys the causality principle. Causality in optics is observed
via Kramers-Kronig relationships doi:10.1002/9783527618156 , which relate the
real and imaginary part of the dielectric susceptibility
$\bm{\epsilon}(\mathbf{r},\omega)=\bm{\epsilon}_{r}(\mathbf{r},\omega)+i\bm{\epsilon}_{i}(\mathbf{r},\omega)$
and magnetic permeability
$\bm{\mu}(\mathbf{r},\omega)=\bm{\mu}_{r}(\mathbf{r},\omega)+i\bm{\mu}_{i}(\mathbf{r},\omega)$
expressed in the frequency domain $\omega$:

$\displaystyle\bm{\sigma}(\mathbf{r},\omega)=\frac{1}{i\pi}\mathrm{PV}\int_{\infty}^{\infty}\frac{\bm{\sigma}(\mathbf{r},\omega)}{\omega^{\prime}-\omega}d\omega^{\prime},$
(9)

with $\bm{\sigma}(\mathbf{r},\omega)=\bm{\epsilon}(\mathbf{r},\omega)$ or
$\bm{\sigma}(\mathbf{r},\omega)=\bm{\mu}(\mathbf{r},\omega)$, and
$\mathrm{PV}$ indicating the principal value. Multiplying Eq. (9) by
$\frac{\nabla\Omega^{\dagger}}{\left|\nabla\Omega^{\dagger}\right|^{\frac{1}{2}}}$
from the right and by
$\frac{\nabla\Omega}{\left|\nabla\Omega\right|^{\frac{1}{2}}}$ from the left,
it becomes:

$\displaystyle\bm{\sigma}^{\prime}(\mathbf{r}^{\prime},\omega)=\frac{1}{i\pi}\mathrm{PV}\int_{\infty}^{\infty}\frac{\nabla\Omega(\mathbf{r})\cdot\bm{\sigma}(\mathbf{r},\omega)\cdot\nabla\Omega^{\dagger}(\mathbf{r})}{|\nabla\Omega^{\dagger}(\mathbf{r})|\cdot\left(\omega^{\prime}-\omega\right)}d\omega^{\prime},$
(10)

which represents Kramers-Kronig relationships for the optical transformed
material properties
$\bm{\sigma}^{\prime}(\mathbf{r}^{\prime},\omega)=\bm{\epsilon}^{\prime}(\mathbf{r}^{\prime},\omega),\bm{\mu}(\mathbf{r}^{\prime},\omega)$
in the transformed physical coordinates
$\mathbf{r}^{\prime}=\bm{\Omega}(\mathbf{r})$. As discussed in the previous
section, the optical permittivity $\bm{\epsilon}=\bm{\epsilon}(\mathbf{r})$
and permeability $\bm{\mu}=\bm{\mu}(\mathbf{r})$ emulated via Eq. (3) are
nondispersive and depend only on the spatial coordinates $\mathbf{r}$. Under
this condition, the causality condition Eq. (10):

$\displaystyle\bm{\sigma}^{\prime}(\mathbf{r}^{\prime})=\frac{\nabla\Omega(\mathbf{r})\cdot\bm{\sigma}(\mathbf{r})\cdot\nabla\Omega^{\dagger}(\mathbf{r})}{|\nabla\Omega^{\dagger}(\mathbf{r})|}\cdot\mathrm{PV}\int_{\infty}^{\infty}\frac{d\omega^{\prime}}{i\pi(\omega^{\prime}-\omega)}=\frac{\nabla\Omega(\mathbf{r})\cdot\bm{\sigma}(\mathbf{r})\cdot\nabla\Omega^{\dagger}(\mathbf{r})}{|\nabla\Omega^{\dagger}(\mathbf{r})|},$
(11)

reduces to transformation optics relationships of the material response,
demonstrating that a material designed with this approach is perfectly causal.
This result relies on the geometrical nature of transformation
$\bm{\Omega}(\mathbf{r})$ used to create the material in (3), which does not
depend on the frequency $\omega$. In this case,
$\mathrm{PV}\int_{\infty}^{\infty}\frac{1}{i\pi(\omega^{\prime}-\omega)}d\omega^{\prime}=1$
and the real and imaginary part of Eq. (11) are decoupled, allowing to
engineer a material that is fully causal.

BIC Q-factor computations. We perform numerical simulations with our parallel
NANOCPP FDTD solver. In the FDTD analysis, we studied the response of the
structure on broadband TE-polarized light from the TFSF source by probing the
electric field inside the nanoscale resonator and then by constructing the
density of states (DOS). To resolve the $Q$-factor of high-Q modes within the
DOS, we analyze the mode’s energy decay by augmenting the Prony algorithm with
a cross-validation model, as explained in the following. The Prony algorithm
decomposes the signal into a sum of complex exponential:

$\hat{f}(t)=\sum_{i=1}^{M}A_{i}e^{-\sigma_{i}t}\cos(\omega_{i}t+\phi_{i})$
(12)

with amplitudes $A_{i}$, decaying constants $\sigma_{i}$, frequencies
$\omega_{i}$ and phases $\phi_{i}$. To compute these terms, we first divide
the FDTD computed over time into two distinct sets: a training set from
$t=t_{0}$ to $t=t_{1}$, and a testing set from $t=t_{1}$ to the simulation
time end $t=t_{2}$. We chose the intervals $t_{1}-t_{0}$ and $t_{2}-t_{1}$ to
be approximately equal.
For a given input number of modes and sampling rate, we first solve linear
prediction model for the training set with least-squares algorithm, find roots
of a characteristic polynomial of the model, and then solve the original set
of linear equations to obtain amplitudes and phases of an individual mode.
Next, we apply this solution to the testing test and calculate the mean
squared error between the FDTD data and the Prony model prediction. We then
construct a surface of MSE values for different values of number of modes $M$
and sampling rates $T_{s}$, and select the values of $M$ and $T_{s}$ that
yield the best prediction. For each mode, we then compute the Q-factors as:

$Q=\frac{\omega_{i}}{2\sigma_{i}}$ (13)

## References

  * (1) June Sang Lee, Ji Yeon Park, Yong Hwan Kim, Seokwoo Jeon, Olivier Ouellette, Edward H Sargent, Dong Ha Kim, and Jerome K Hyun.  Ultrahigh resolution and color gamut with scattering-reducing transmissive pixels.  Nature communications, 10(1):1–9, 2019.
  * (2) Kirill Koshelev, Gael Favraud, Andrey Bogdanov, Yuri Kivshar, and Andrea Fratalocchi.  Nonradiating photonics with resonant dielectric nanostructures.  Nanophotonics, 8(5):725–745, 2019.
  * (3) Kirill Koshelev, Sergey Kruk, Elizaveta Melik-Gaykazyan, Jae-Hyuck Choi, Andrey Bogdanov, Hong-Gyu Park, and Yuri Kivshar.  Subwavelength dielectric resonators for nonlinear nanophotonics.  Science, 367(6475):288–292, 2020.
  * (4) Mengfei Wu, Lu Ding, Randy P. Sabatini, Laxmi Kishore Sagar, Golam Bappi, Ramón Paniagua-Domínguez, Edward H. Sargent, and Arseniy I. Kuznetsov.  Bound state in the continuum in nanoantenna-coupled slab waveguide enables low-threshold quantum-dot lasing.  Nano Letters, 21(22):9754, Nov 2021.
  * (5) Ludwig Huttenhofer, Andreas Tittl, Lucca Kuhner, Emiliano Cortés, and Stefan A Maier.  Anapole-assisted absorption engineering in arrays of coupled amorphous gallium phosphide nanodisks.  ACS Photonics, 8(5):1469–1476, 2021.
  * (6) Juan Wang, Julius Kühne, Theodosios Karamanos, Carsten Rockstuhl, Stefan A. Maier, and Andreas Tittl.  All-dielectric crescent metasurface sensor driven by bound states in the continuum.  Advanced Functional Materials, 31(46):2104652, 2021.
  * (7) Mazhar E. Nasir, Alexey V. Krasavin, R. Margoth Córdova-Castro, Cillian P. T. McPolin, Jean-Sebastien G. Bouillard, Pan Wang, and Anatoly V. Zayats.  Mode engineering in large arrays of coupled plasmonic–dielectric nanoantennas.  Advanced Optical Materials, 9(11):2001467, 2021.
  * (8) Daria Smirnova, Daniel Leykam, Yidong Chong, and Yuri Kivshar.  Nonlinear topological photonics.  Applied Physics Reviews, 7(2):021306, 2020.
  * (9) Yafeng Chen, Zhihao Lan, Jensen Li, and Jie Zhu.  Topologically protected second harmonic generation via doubly resonant high-order photonic modes.  Physical Review B, 104(15):155421, 2021.
  * (10) M. Fuechsle, J. A. Miwa, S. Mahapatra, H. Ryu, S. Lee, O. Warschkow, L. C. L. Hollenberg, G. Klimeck, and M. Y. Simmons.  A single-atom transistor.  Nat. Nanotech., 7:242–246, 2012.
  * (11) S. B. Desai, S. R. Madhvapathy, A. B. Sachid, J. P. Llinas, Q. Wang, G. H. Ahn, G. Pitner, M. J. Kim, J. Bokor, C. Hu, H.-S. P. Wong, and A. Javey.  Mos2 transistors with 1-nanometer gate lengths.  Science, 354:99–102, 2016.
  * (12) Stefan Alexander Maier.  Plasmonics: Fundamentals and Applications.  Springer, London, 2007.
  * (13) Jianfeng Huang, Changxu Liu, Yihan Zhu, Silvia Masala, Erkki Alarousu, Yu Han, and Andrea Fratalocchi.  Harnessing structural darkness in the visible and infrared wavelengths for a new source of light.  Nat. Nanotech., 11:60–66, 2016.
  * (14) Yi-Fan Huang, Surojit Chattopadhyay, Yi-Jun Jen, Cheng-Yu Peng, Tze-An Liu, Yu-Kuei Hsu, Ci-Ling Pan, Hung-Chun Lo, Chih-Hsun Hsu, Yuan-Huei Chang, Chih-Shan Lee, Kuei-Hsien Chen, and Li-Chyong Chen.  Improved broadband and quasi-omnidirectional anti-reflection properties with biomimetic silicon nanostructures.  Nat. Nanotech., 2:770–774, 2007.
  * (15) Kirill Koshelev and Yuri Kivshar.  Dielectric resonant metaphotonics.  ACS Photonics, 8(1):102–112, 2020.
  * (16) R Won.  Into the ‘mie-tronic’ era.  Nature Photonics, 13:585–587, 2019.
  * (17) Xiaopeng Shen and Tie Jun Cui.  Planar plasmonic metamaterial on a thin film with nearly zero thickness.  Applied Physics Letters, 102(21):211909, 2013.
  * (18) Jordi Bonache, Gerard Zamora, Ferran Paredes, Simone Zuffanelli, Pau Aguilà, and Ferran Martín.  Controlling the electromagnetic field confinement with metamaterials.  Scientific reports, 6(1):1–6, 2016.
  * (19) Y Yermakov, Anna A Hurshkainen, Dmitry A Dobrykh, Polina V Kapitanova, Ivan V Iorsh, Stanislav B Glybovski, and Andrey A Bogdanov.  Experimental observation of hybrid te-tm polarized surface waves supported by a hyperbolic metasurface.  Physical Review B, 98(19):195404, 2018.
  * (20) Craig F. Bohren Donald R. Huffman.  Electromagnetic Theory.  John Wiley & Sons, Ltd, 2007.
  * (21) Robert W. Boyd.  Nonlinear Optics, Third Edition.  Academic Press, Inc., USA, 3rd edition, 2008.
  * (22) J. B. Pendry, A. Aubry, D. R. Smith, and S. A. Maier.  Transformation optics and subwavelength control of light.  Science, 337(6094):549–552, 2012.
  * (23) D. Felbacq and G. Bouchitté.  Metamaterials Modelling and Design.  Pan Stanford, New York, 2017.
  * (24) Martin McCall, John B Pendry, Vincenzo Galdi, Yun Lai, S A R Horsley, Jensen Li, Jian Zhu, Rhiannon C Mitchell-Thomas, Oscar Quevedo-Teruel, Philippe Tassin, Vincent Ginis, Enrica Martini, Gabriele Minatti, Stefano Maci, Mahsa Ebrahimpouri, Yang Hao, Paul Kinsler, Jonathan Gratus, Joseph M Lukens, Andrew M Weiner, Ulf Leonhardt, Igor I Smolyaninov, Vera N Smolyaninova, Robert T Thompson, Martin Wegener, Muamer Kadic, and Steven A Cummer.  Roadmap on transformation optics.  Journal of Optics, 20(6):063001, may 2018.
  * (25) Sunae So and Junsuk Rho.  Geometrically flat hyperlens designed by transformation optics.  Journal of Physics D: Applied Physics, 52(19):194003, mar 2019.
  * (26) Laura Pilozzi, Francis A Farrelly, Giulia Marcucci, and Claudio Conti.  Machine learning inverse problem for topological photonics.  Communications Physics, 1(1):1–7, 2018.
  * (27) Sunae So, Trevon Badloe, Jaebum Noh, Jorge Bravo-Abad, and Junsuk Rho.  Deep learning enabled inverse design in nanophotonics.  Nanophotonics, 9(5):1041–1057, 2020.
  * (28) R. Schinzinger and P.A.A. Laura.  Conformal Mapping: Methods and Applications.  Elsevier, 1991.
  * (29) Bengt Fornberg.  A Practical Guide to Pseudospectral Methods.  Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 1996.
  * (30) Ben-yu Guo, Jie Shen, and Zhong-qing Wang.  Chebyshev rational spectral and pseudospectral methods on a semi-infinite interval.  International Journal for Numerical Methods in Engineering, 53(1):65–84, 2002.
  * (31) Erwin Kreyszig, Herbert Kreyszig, and E. J. Norminton.  Advanced Engineering Mathematics.  Wiley, Hoboken, NJ, tenth edition, 2011.
  * (32) Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  The elements of statistical learning: data mining, inference and prediction.  Springer, 2 edition, 2009.
  * (33) Sameer Agarwal, Keir Mierle, and Others.  Ceres solver.  http://ceres-solver.org.
  * (34) Safa Kasap and Peter Capper.  Springer handbook of electronic and photonic materials.  Springer, 2017.
  * (35) Denis G Baranov, Dmitry A Zuev, Sergey I Lepeshov, Oleg V Kotov, Alexander E Krasnok, Andrey B Evlyukhin, and Boris N Chichkov.  All-dielectric nanophotonics: the quest for better materials and fabrication techniques.  Optica, 4(7):814–825, 2017.
  * (36) Andrey E Miroshnichenko, Andrey B Evlyukhin, Ye Feng Yu, Reuben M Bakker, Arkadi Chipouline, Arseniy I Kuznetsov, Boris Luk’yanchuk, Boris N Chichkov, and Yuri S Kivshar.  Nonradiating anapole modes in dielectric nanoparticles.  Nature communications, 6(1):1–8, 2015.
  * (37) Sheng Liu, Gordon A. Keeler, John L. Reno, Michael B. Sinclair, and Igal Brener.  Iii–v semiconductor nanoresonators—a new strategy for passive, active, and nonlinear all-dielectric metamaterials.  Advanced Optical Materials, 4(10):1457–1462, 2016.
  * (38) M Schmid.  Review on light management by nanostructures in chalcopyrite solar cells.  Semiconductor Science and Technology, 32(4):043003, mar 2017.
  * (39) Zhong-Jian Yang, Ruibin Jiang, Xiaolu Zhuo, Ya-Ming Xie, Jianfang Wang, and Hai-Qing Lin.  Dielectric nanoresonators for light manipulation.  Physics Reports, 701:1 – 50, 2017.
  * (40) Mikhail V. Rybin, Kirill L. Koshelev, Zarina F. Sadrieva, Kirill B. Samusev, Andrey A. Bogdanov, Mikhail F. Limonov, and Yuri S. Kivshar.  High-$q$ supercavity modes in subwavelength dielectric resonators.  Phys. Rev. Lett., 119:243901, Dec 2017.
  * (41) Juan S. Totero Gongora, Andrey E. Miroshnichenko, Yuri S. Kivshar, and Andrea Fratalocchi.  Anapole nanolasers for mode-locking and ultrafast pulse generation.  Nature Communications, 8:15535, 2017.
  * (42) Elizaveta Melik-Gaykazyan, Kirill Koshelev, Jae-Hyuck Choi, Sergey S Kruk, Andrey Bogdanov, Hong-Gyu Park, and Yuri Kivshar.  From fano to quasi-bic resonances in individual dielectric nanoantennas.  Nano Letters, 21(4):1765–1771, 2021.
  * (43) Yuri Kivshar.  All-dielectric meta-optics and non-linear nanophotonics.  Natl Sci Rev, 5:144–148, 2018.
  * (44) Alexey Zhizhchenko, Sergey Syubaev, Alexander Berestennikov, Alexey V Yulin, Alexey Porfirev, Anatoly Pushkarev, Ivan Shishkin, Kirill Golokhvast, Andrey A Bogdanov, Anvar A Zakhidov, et al.  Single-mode lasing from imprinted halide-perovskite microdisks.  ACS nano, 13(4):4140–4147, 2019.
  * (45) Ju Young Kim, Hyowook Kim, Bong Hoon Kim, Taeyong Chang, Joonwon Lim, Hyeong Min Jin, Jeong Ho Mun, Young Joo Choi, Kyungjae Chung, Jonghwa Shin, et al.  Highly tunable refractive index visible-light metasurface from block copolymer self-assembly.  Nature communications, 7(1):1–9, 2016.
  * (46) Andrey A. Bogdanov, Andrea Fratalocchi, and Yuri Kivshar.  The science of harnessing light’s darkness.  Nanophotonics, 10(17):4171–4173, 2021.
  * (47) V. Pagneux.  Multimodal admittance method in waveguides and singularity behavior at high frequencies.  Journal of Computational and Applied Mathematics, 234(6):1834–1841, 2010.  Eighth International Conference on Mathematical and Numerical Aspects of Waves (Waves 2007).
  * (48) Ayca Yalcin, Ketul C Popat, John C Aldridge, Tejal A Desai, John Hryniewicz, Nabil Chbouki, Brent E Little, Oliver King, Vien Van, Sai Chu, et al.  Optical sensing of biomolecules using microring resonators.  IEEE Journal of Selected Topics in Quantum Electronics, 12(1):148–155, 2006.
  * (49) Gun-Duk Kim, Geun-Sik Son, Hak-Soon Lee, Ki-Do Kim, and Sang-Shin Lee.  Integrated photonic glucose biosensor using a vertically coupled microring resonator in polymers.  Optics Communications, 281(18):4644–4647, 2008.
  * (50) Pauline Girault, Nathalie Lorrain, Luiz Poffo, Mohammed Guendouz, Jonathan Lemaitre, Christiane Carré, Michel Gadonna, Dominique Bosc, and Guillaume Vignaud.  Integrated polymer micro-ring resonators for optical sensing applications.  Journal of applied Physics, 117(10):104504, 2015.
  * (51) Shon Schmidt, Jonas Flueckiger, WenXuan Wu, Samantha M Grist, Sahba Talebi Fard, Valentina Donzella, Pakapreud Khumwan, Emily R Thompson, Qian Wang, Pavel Kulik, et al.  Improving the performance of silicon photonic rings, disks, and bragg gratings for use in label-free biosensing.  In Biosensing and Nanomedicine VII, volume 9166, page 91660M, 2014\.
  * (52) Penghao Liu and Yaocheng Shi.  Simultaneous measurement of refractive index and temperature using a dual polarization ring.  Applied optics, 55(13):3537–3541, 2016.
  * (53) Heming Wei and Sridhar Krishnaswamy.  Direct laser writing polymer micro-resonators for refractive index sensors.  IEEE Photonics Technology Letters, 28(24):2819–2822, 2016.
  * (54) Ilya Goykhman, Boris Desiatov, and Uriel Levy.  Ultrathin silicon nitride microring resonator for biophotonic applications at 970 nm wavelength.  Applied Physics Letters, 97(8):081108, 2010.
  * (55) Ting Tang, Xiang Wu, Liying Liu, and Lei Xu.  Packaged optofluidic microbubble resonators for optical sensing.  Applied Optics, 55(2):395–399, 2016.
  * (56) Juejun Hu, Nathan Carlie, Ning-Ning Feng, Laeticia Petit, Anu Agarwal, Kathleen Richardson, and Lionel Kimerling.  Planar waveguide-coupled, high-index-contrast, high-q resonators in chalcogenide glass for sensing.  Optics letters, 33(21):2500–2502, 2008.
  * (57) Xiaokun Wang, Xiaowei Guan, Qiangsheng Huang, Jiajiu Zheng, Yaocheng Shi, and Daoxin Dai.  Suspended ultra-small disk resonator on silicon for optical sensing.  Optics letters, 38(24):5405–5408, 2013.
  * (58) Yeseul Kim and Hansuek Lee.  On-chip label-free biosensing based on active whispering gallery mode resonators pumped by a light-emitting diode.  Optics express, 27(23):34405–34415, 2019.
  * (59) Samantha M Grist, Shon A Schmidt, Jonas Flueckiger, Valentina Donzella, Wei Shi, Sahba Talebi Fard, James T Kirk, Daniel M Ratner, Karen C Cheung, and Lukas Chrostowski.  Silicon photonic micro-disk resonators for label-free biosensing.  Optics express, 21(7):7994–8006, 2013.
  * (60) S Hamed Mirsadeghi, Ellen Schelew, and Jeff F Young.  Photonic crystal slot-microcavity circuit implemented in silicon-on-insulator: High q operation in solvent without undercutting.  Applied Physics Letters, 102(13):131115, 2013.
  * (61) A Di Falco, L O’faolain, and TF Krauss.  Chemical sensing in slotted photonic crystal heterostructure cavities.  Applied physics letters, 94(6):063503, 2009.
  * (62) Jaime García-Rupérez, Veronica Toccafondo, María José Bañuls, Javier García Castelló, Amadeu Griol, Sergio Peransi-Llopis, and Ángel Maquieira.  Label-free antibody detection using band edge fringes in soi planar photonic crystal waveguides in the slow-light regime.  Optics Express, 18(23):24276–24286, 2010.
  * (63) Qiang Geng, Dien Wang, Pengfei Chen, and Shih-Chi Chen.  Ultrafast multi-focus 3-d nano-fabrication based on two-photon polymerization.  Nature communications, 10(1):1–7, 2019.
  * (64) Brian Elder, Rajan Neupane, Eric Tokita, Udayan Ghosh, Samuel Hales, and Yong Lin Kong.  Nanomaterial patterning in 3d printing.  Advanced Materials, 32(17):1907142, 2020.

Figure 1:  Designing materials via geometrical deformations: general
principles. a) Ideal structure for controlling light, characterized by a
reflective substrate with a material possessing an arbitrary refractive index
distribution on top. b) Equivalent deformed structure found by inverse design
and composed of a reflective substrate in air that possesses the same
electromagnetic behavior of (a). c) Mapping
$\mathbf{r^{\prime}}=\mathbf{\Omega}(\mathbf{r})$ function connecting both
domains.

Figure 2:  Ultra-high refractive index nanoresonators. a) Target Refractive
index profile (circle markers) and inverse design refractive index (solid red
line). b) Corresponding deformation in the accesible space. c) Convergence of
nonlinear regression for a varying number of modes $M$. d) Spatial
distribution of inverse designed refractive index in the design space. e-g)
Refractive index distribution with amplitudes $n=20$ and $n=70$, respectively,
and f-h corresponding deformations.  Figure 3:  BIC states of inverse designed
ultra-high index resonators. a) Eigenvalue amplitude $|\gamma(0)|$ vs.
frequency $k$. b)-g) Electromagnetic energy spatial distribution for singular
states in the accessible $\mathbf{r^{\prime}}$ and design $\mathbf{r}$ spaces,
respectively.  Figure 4:  State-of-the-art in RI optical sensors. a)
Sensitivity vs. sensor footprint. b) Single figure of merit displaying the
sensitivity per sensor footprint.  Figure 5:  RI sensor based on inversely
designed ultra-high index material a) Sensor idea in the design space
comprising a single resonator with ultra-high refractive index of $n=100$ and
an analyte placed on top. b) Corresponding configuration in the accessible
space. c) Density of states (DOS) of the resonator (solid blue line) with BIC
states possessing $Q=5200$). d) Electromagnetic energy spatial distribution of
the BIC state in the accessible space.  Figure 6:  Performance of the RI
sensing device. a) Spectral shift of the BIC resonant wavelength for a varying
refractive index $n$ of the analyte. b) Wavelength shift $\Delta\lambda$ and
differential sensitivity $S$ vs. refractive index change of the analyte


# Origin of negative electrocaloric effect in Pnma-type antiferroelectric
perovskites

Ningbo Fan  Institute of Theoretical and Applied Physics, Soochow University,
Suzhou 215006, China  School of Physical Science and Technology, Soochow
University, Suzhou 215006, China  Jorge Íñiguez  Materials Research and
Technology Department, Luxembourg Institute of Science and Technology (LIST),
Avenue des Hauts-Fourneaux 5, L-4362 Esch/Alzette, Luxembourg  Department of
Physics and Materials Science, University of Luxembourg, Rue du Brill 41,
L-4422 Belvaux, Luxembourg  L. Bellaiche  Physics Department and Institute for
Nanoscience and Engineering, University of Arkansas, Fayetteville, Arkansas
72701, USA  Bin Xu<EMAIL_ADDRESS>Institute of Theoretical and Applied
Physics, Soochow University, Suzhou 215006, China  School of Physical Science
and Technology, Soochow University, Suzhou 215006, China

###### Abstract

Anomalous electrocaloric effect (ECE) with decreasing temperature upon
application of an electric field is known to occur in antiferroelectrics
(AFEs), and previous understanding refers to the field-induced canting of
electric dipoles if there is no phase transitions. Here, we use a first-
principle-based method to study the ECE in Nd-substituted BiFeO3 (BNFO)
perovskite solid solutions, which has the Pnma-type AFE ground state. We
demonstrate another scenario to achieve and explain anomalous ECE, emphasizing
that explicit consideration of octahedral tiltings is indispensable for a
correct understanding. This mechanism may be general for AFEs for which the
antipolar mode is not the primary order parameter. We also find that the
negative ECE can reach a large magnitude in BNFO.

Electrocaloric effect (ECE) can make temperature change via adiabatic
application (or removal) of an electric field, providing an efficient approach
for cooling or heating Correlia and Zhang (2014); Fähler  _et al._ (2012);
Moya  _et al._ (2014). While ferroelectric (FE) or relaxor materials typically
have “normal” positive sign of ECE, i.e., the temperature increases by
applying a voltage, antiferroelectrics (AFEs) are known to have anomalous ECE
that can yield an opposite sign Pirc  _et al._ (2014); Geng  _et al._ (2015);
Vales-Castro  _et al._ (2021) . These two types of ECE can be utilized in
combination to improve the performance of cooling/heating devices.

Such negative (or inverse) caloric effect is also known to occur in other
occasions, e.g., magnetic Heusler alloys, transitions between FE phases of
different polarization directions, and application of an electric field
against the polarization of a FE phase without switching Gottschall  _et al._
(2018); Marathe  _et al._ (2017); Grünebohm  _et al._ (2018); however, its
origin in AFEs is less well understood. AFEs materials are characterized by
anti-polar atomic distortions that can be switched to a FE state under an
electric field, and two mechanisms to explain their negative ECE with no AFE-
FE transition have been proposed: 1) the “dipole-canting” model that dipolar
entropy increases by misaligning the anti-parallel dipoles upon application of
the field Geng  _et al._ (2015); 2) the perturbative theory based on the
Maxwell relation that only temperature and electric field dependencies of
polarization need to be considered Graf and Íñiguez (2021). Interestingly, all
these mechanisms only take the electric degrees of freedom explicitly into
account. In contrast, most of the known AFEs are neither proper type (that is,
the AFE phase is rarely driven by an AFE soft mode Milesi-Brault  _et al._
(2020)), nor systems with the anti-polar mode being the only significant order
parameter. In fact, quite often, the AFE mode is secondary and coupled to
other degrees of freedom, such as the octahedral tiltings in perovskites. For
instance, PZO has a strong instability of anti-phase octahedral tilting Xu
_et al._ (2019), while in Pnma-type perovskite, such as rare earth
orthoferrites and CsPbI3, the anti-polar distortion arises from the
condensation of both the in-plane anti-phase ($\omega_{R_{x,y}}$) and out-of-
plane in-phase tiltings ($\omega_{M_{z}}$) via trilinear coupling Bellaiche
and Íñiguez (2013). Although these tilting modes are non-polar, they couple
strongly with the polar and anti-polar modes so that they can be influenced by
the electric field as well, and in consequence contribute to the ECE.

To get a deeper understanding of the (negative) ECE in AFE, analysis based on
Landau models involving the most relevant degrees of freedom have been proved
to be very useful Planes  _et al._ (2014); Zhang  _et al._ (2011); Ma  _et
al._ (2016); Edström and Ederer (2020), and it may thus be necessary that all
the important order parameters are taken into consideration. Furthermore, some
previous phenomenological models are often over simplified, since only one
dimension is assumed Pirc  _et al._ (2014). In reality, the direction of the
applied field with respect to the crystallographic axis should have different
effects regarding ECE. In this Letter, we take the antiferroelectric Nd
substituted BiFeO3 (BNFO) solid solution as an example and demonstrate that
the octahedral tiltings can have very important effect on the sign and
magnitude of the ECE. We also construct a phenomenological model that allows
us to rationalize the contributions of each degree of freedom. In particular,
the dipoles alone are found to be insufficient to explain the negative ECE,
while contributions from the in-phase and anti-phase tilting modes are
indispensible. Moreover, BNFO is predicted to yield rather large negative ECE
close to the AFE-to-FE transition.

BiFeO3 (BFO) stabilizes in a R3c ground state, but rare-earth doping with
composition larger than 20$-$30% is sufficient to alter it to the Pnma
structure Karimi  _et al._ (2009a). Here, we adopt the effective Hamiltonian
scheme of Ref. Xu  _et al._ (2015a) to study the Bi0.6Nd0.4FeO3 solid solution
under electric field at finite temperatures. With this composition, BNFO is
stabilized in the AFE Pnma phase at room temperature and can transform to a FE
state under an electric field Karimi  _et al._ (2009a, b); Kan  _et al._
(2010); Xu  _et al._ (2015a, 2017). The solid solutions are simulated by a
12$\times$12$\times$12 supercell (containing 8,640 atoms) using Monte-Carlo
(MC) simulations, in which the Bi and Nd atoms are randomly distributed (see
Supplemental Material (SM) Section S1 for details sm ).

Figure 1: Effect of $\mathcal{E}$ field on the temperature dependence of the
order parameters. (a),(b) $[001]$ field of 0.87 MV/cm (solid lines), in
comparison with zero field (dashed lines), (c),(d) ibid for $[1\bar{1}0]$
field of 0.87 MV/cm, (e),(f) ibid for $[110]$ field of 0.61 MV/cm. The four
relevant order parameters are anti-phase octahedral tilting $\omega_{R}$, in-
phase tilting $\omega_{M}$, polarization $P$, and antiferroelectric vector
$X$. Note that $P_{x}$=$-P_{y}$ for $[1\bar{1}0]$ field. The vertical dotted
(or dashed) lines delimit different phases under finite (or zero) field.

First, let us check how the order parameters are influenced by the electric
field $\mathcal{E}$ 1no , using the effective Hamiltomian scheme. Note that
the initial Pnma structure has zero polarization, an in-plane anti-polar
vector $X$ along the $[110]$ direction, and $a^{-}a^{-}c^{+}$ tilting in
Glazer’s notation (which corresponds to finite $x$\- and $y$-components of the
antiphase tilting vector, $\omega_{R,x}=\omega_{R,y}$, and finite
$z$-component of the in-phase tilting vector, $\omega_{M,z}$). Three
representative field directions are investigated (Fig. 1), together with the
zero-field data (dashed curves) for comparison (for other $\mathcal{E}$
magnitudes, see SM Section S2 sm ). With no applied field, the Pnma phase
transforms to the paraelectric (PE) cubic phase at 1400 K (dashed lines in
Fig. 1).

Under $[001]$ field, we consider a representative case with $\mathcal{E}=$
0.87 MV/cm, with which the AFE state transforms to the ferroelectric P4mm
phase at 880 K, characterized by a large polarization $P_{z}$ along $[001]$
and no octahedral tiltings. In the AFE state, one can see that the temperature
dependence of $\omega_{M,z}$ does not differ much from the zero-field case,
whereas $P$ and $\omega_{R,(x,y)}$ ($\omega_{R,x}=\omega_{R,y}$) show apparent
changes. The moderate field-induced change of $X$ can be understood to a good
approximation via the change in $\omega_{R,(x,y)}$ since $X_{(x,y)}$ should be
proportional to the product $\omega_{R,(x,y)}\omega_{M,z}$ – as a result of a
trilinear coupling between $X_{(x,y)}$, $\omega_{R,(x,y)}$ and $\omega_{M,z}$
Xu  _et al._ (2015b). Finite $P_{z}$ is induced by the field, whereas the
field-induced suppression of $\omega_{R,(x,y)}$ is due to the competitive
coupling with $P_{z}$.

Moreover, if the field is applied along $[1\bar{1}0]$ (Figs. 1c and 1d), $P$
is induced along the same direction, i.e., a finite $P_{x}$=$-P_{y}$ first
develops, and a transition to a FE Cc phase occurs at 790 K, a structure
characterized by a polarization in $\langle uuv\rangle$ direction ($u>v$) and
$a^{-}a^{-}c^{-}$ octahedral tiltings ($a>c$). The third $\mathcal{E}$
direction is along $[110]$, along which $P$ develops, while $\omega_{M,z}$ is
much suppressed and $\omega_{R,(x,y)}$ is more or less unchanged within the
AFE-based state 2no . Such AFE phase then transforms also into a Cc phase at
780 K (Fig. 1e and 1f).

Let us now concentrate on the ECE coefficient $\alpha=\frac{\partial
T}{\partial\mathcal{E}}\rvert_{S}$ with $T$ being temperature and $S$ being
entropy, which can be calculated from the cumulant formula using outputs of
the MC simulationsJiang  _et al._ (2017, 2018); Bin-Omran  _et al._ (2016);
Jiang  _et al._ (2021a).

$\displaystyle\alpha_{\text{MC}}=-Z^{*}a_{\text{latt}}T\left\\{\frac{\left<\left|\bm{u}\right|E_{\text{tot}}\right>-\left<\left|\bm{u}\right|\right>\left<E_{\text{tot}}\right>}{\left<E_{\text{tot}}^{2}\right>-\left<E_{\text{tot}}\right>^{2}+\frac{21(k_{\text{B}}T)^{2}}{2N}}\right\\}\;,$
(1)

where $Z^{*}$ is the Born effective charge associated with the local mode,
$a_{\text{latt}}$ represents the lattice constant of the five-atom pseudo-
cubic perovskite cell, $T$ is the simulation temperature,
$\left|\bm{u}\right|$ is the supercell average of the magnitude of the local
mode, $E_{\text{tot}}$ is the total energy given by the effective Hamiltonian,
$k_{\text{B}}$ is the Boltzmann constant, $N$ is the number of sites in the
supercell, and $\left<\;\right>$ denotes average over the MC sweeps at a given
temperature.

For fields applied along the $[001]$ direction (Fig. 2a), similar to the case
in PZO-based AFE Geng  _et al._ (2015); Vales-Castro  _et al._ (2021),
$\alpha$ is negative in the AFE-based state, and its magnitude increases with
temperature, which maximizes at the transition point where the AFE-based state
disappears. Across the phase transition, $\alpha$ jumps to be positive in the
FE state, then (slightly) increases with temperature, as such qualitative
temperature dependence is known for ferroelectrics Jiang  _et al._ (2021b).

In order to have insightful analysis of the ECE, Figs. 2b and 2c also report
two other quantities related to electro-caloric response, namely the total
isothermal change in entropy $\Delta S$, and the adiabatic temperature change
$\Delta T$, as well as their individual contributions. Practically, in order
to be able to compute the total entropy change, we consider the following
Landau model for Pnma-type AFE Xu  _et al._ (2015b).

$\displaystyle F=$ $\displaystyle
F_{0}+\frac{1}{2}a_{\omega_{R}}(T)\omega_{R}^{2}+\frac{1}{4}b_{\omega_{R}}\omega_{R}^{4}+\frac{1}{2}a_{\omega_{M}}(T)\omega_{M}^{2}+\frac{1}{4}b_{\omega_{M}}\omega_{M}^{4}$
$\displaystyle+\frac{1}{2}a_{X}(T)X^{2}+\frac{1}{4}b_{X}X^{4}+\frac{1}{2}a_{P}(T)P^{2}+\frac{1}{4}b_{P}P^{4}$
$\displaystyle-\mathcal{E}P-cX\omega_{R}\omega_{M}+\frac{1}{2}d_{1}P^{2}\omega_{R}^{2}+\frac{1}{2}d_{2}P^{2}\omega_{M}^{2}\;.$
(2)

Note that $F_{0}$ is the field-independent part of the free energy, only the
quadratic coefficients $a_{\text{op}}=A_{\text{op}}(T-T_{0}^{\text{op}})$ have
explicit temperature dependence, with “op” being the order parameter of
$\omega_{R}$, $\omega_{M}$, $P$, or $X$, respectively, and $T_{0}^{\text{op}}$
is a transition temperature for each order parameter. Also note that the
tiltings are known to be the primary order parameters in many Pnma-systems,
while $X$ is secondary arising from the $cX\omega_{R}\omega_{M}$ trilinear
coupling involving in-phase and anti-phase tiltings Zhao  _et al._ (2014).
Recalling that $S=-\frac{\partial F}{\partial
T}\rvert_{\mathcal{E},\text{op}}$, the change of entropy associated with a
change of $\mathcal{E}$ can be written as:

$\displaystyle\Delta S(T)=$ $\displaystyle-
A_{\omega_{R}}\left[\omega_{R}^{2}(T,\mathcal{E}_{2})-\omega_{R}^{2}(T,\mathcal{E}_{1})\right]$
$\displaystyle-
A_{\omega_{M}}\left[\omega_{M}^{2}(T,\mathcal{E}_{2})-\omega_{M}^{2}(T,\mathcal{E}_{1})\right]$
$\displaystyle-
A_{P}\left[P^{2}(T,\mathcal{E}_{2})-P^{2}(T,\mathcal{E}_{1})\right]$
$\displaystyle-
A_{X}\left[X^{2}(T,\mathcal{E}_{2})-X^{2}(T,\mathcal{E}_{1})\right]$
$\displaystyle=$ $\displaystyle\Delta S_{\omega_{R}}+\Delta
S_{\omega_{M}}+\Delta S_{P}+\Delta S_{X}\;,$ (3)

where $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ are the initial and final
electric field, respectively, and the values of the order parameters are
obtained from MC simulations. The adiabatic temperature change due to ECE can
be further obtained via $\Delta T=-(T/C_{E})\Delta S$, where $C_{E}$ is the
phonon specific heat Kutnjak  _et al._ (1999) – which is calculated here by
density functional theory (see SM Section S3 sm ). We thus have:

$\displaystyle\Delta T=-\frac{T}{C_{E}}\Delta S=\Delta T_{\omega_{R}}+\Delta
T_{\omega_{M}}+\Delta T_{P}+\Delta T_{X}\;.$ (4)

Equations (Origin of negative electrocaloric effect in Pnma-type
antiferroelectric perovskites) and (4) are used to evaluate the entropy and
temperature change (as well as their individual contributions associated with
$\omega_{R}$, $\omega_{M}$, $P$ and $X$), in which each order parameter is
directly obtained from MC simulations while the coefficients
$A_{\text{op}}=-a_{\text{op}}^{\text{0-K}}/T_{0}^{\text{op}}$ are extracted by
the following two steps: (1) finding $a_{\text{op}}^{\text{0-K}}$ from fitting
the 0 K double-well energies based on our effective Hamiltonian to the Landau
model (Eq. Origin of negative electrocaloric effect in Pnma-type
antiferroelectric perovskites); (2) determining $T_{0}^{\text{op}}$ from
fitting $\alpha$ calculated from the cumulant formula to that based on the
Landau model with $\alpha_{\text{fit}}=\Delta T/(2\Delta\mathcal{E})$ for
$\Delta\mathcal{E}$=0.02 MV/cm (see SM Section S4 for detailssm ). Note that
$\alpha_{\text{fit}}$ is shown in Fig. 2a and agrees rather well (over a large
temperature range) with the $\alpha$ directly obtained from the MC simulations
in the AFE-based state, therefore demonstrating the relevance of the
aforementioned Landau model. Also note that the obtained $T_{0}^{\text{op}}$
reflects the depth of the energy wells.

Figure 2: Temperature dependence of the calculated ECE coefficient $\alpha$,
isothermal entropy change, and temperature change under (a)-(c) $[001]$ (0.87
MV/cm), (d)-(f) $[1\bar{1}0]$ (0.87 MV/cm), and (g),(i) $[110]$ (0.61 MV/cm)
electric field. Colors in $\Delta S$ and $\Delta T$ denote the total and
individual contributions from each order parameter. $\alpha_{\text{MC}}$ is
computed by the cumulant approach, and $\alpha_{\text{fit}}$ is a fitted curve
based on the Landau model.

One can then realize from Eq. (4) two important features for each $\Delta
T_{\text{op}}$ associated with a specific order parameter (op), that are: (1)
such change is directly proportional to the $A_{\text{op}}$ coefficient
related to that order parameter; and (2) such change involves the difference
in the square of the order parameter between two electric fields, e.g.,
$\left[\omega_{R}^{2}(T,\mathcal{E}_{2})-\omega_{R}^{2}(T,\mathcal{E}_{1})\right]$
for the anti-phase tilting. These four different individual changes in
temperature are reported in Fig. 2c for field along $[001]$.

Here, we only focus on the ECE in the AFE state, since latent heat needs to be
considered when first-order AFE-to-FE transition occurs – which is not
included in Eq. Origin of negative electrocaloric effect in Pnma-type
antiferroelectric perovskites. Therefore, ECE involving the FE and PE states
are not investigated within the scope of the present study. $\Delta S$ and
$\Delta T$ in Figs. 2b and 2c, respectively, are practically calculated from
Eqs. (Origin of negative electrocaloric effect in Pnma-type antiferroelectric
perovskites) and (4) with $\mathcal{E}_{1}$=0 and $\mathcal{E}_{2}$=0.87 MV/cm
being along $[001]$. We numerically find that $\Delta T$ is negative and bears
qualitatively similar $T$ dependence as that of $\alpha$, i.e., its magnitude
increases with increasing temperature and maximizes at the transition point.
On the other hand, $\Delta S$ has an opposite sign to that of $\Delta T$, as
dictated by their relationship $\Delta T=-(T/C_{E})\Delta S$. To understand
how ECE is contributed by each order parameter, we can focus on the different
$\Delta T_{\text{op}}$ displayed in Fig. 2c. For this $[001]$ field direction,
$\omega_{R}$, $\omega_{M}$, and $X$ show negative contributions, with the
magnitude from $\omega_{R}$ being the largest and that from $\omega_{M}$ being
negligible, together with the rather large positive contribution from $P$.
According to Eq. 4, the negative (or positive) signs of the four different
$\Delta T_{\text{op}}$ can be understood by the fact that the corresponding
order parameter is suppressed (or enhanced) with the application of the field.
In fact, as suggested by Eqs. (Origin of negative electrocaloric effect in
Pnma-type antiferroelectric perovskites) and (4), a large magnitude of $\Delta
T_{\text{op}}$ relies on two factors, either those with large field-induced
changes in their order parameter such as $\omega_{R}$ and $P$ here, or those
with large coefficient $A_{\text{op}}$ such as $X$ here 3no , whereas $\Delta
T_{\omega_{M}}$ for field along $[001]$ is negligible since both
$A_{\omega_{M}}$ and the field-induced change in $\omega_{M}$ are very small.

Figure 3: The ECE temperature change as a function of temperature under
various electric field magnitude and directions. (a) $[001]$, (b)
$[1\bar{1}0]$, and (c) $[110]$ field.

The same procedure is adopted to study the field directions along
$[1\bar{1}0]$ and $[110]$. For the $[1\bar{1}0]$ field, as shown in Figs.
2d-f, a qualitatively similar ECE behavior is found compared with those from
the $[001]$ field (namely, a negative $\alpha$ having a magnitude increasing
when approaching the AFE-to-FE transition), except that the predicted
$\alpha$, $\Delta S$ and $\Delta T$ are smaller in magnitude, e.g., the
largest negative $\alpha$ = $-$1.0 mK$\cdot$kV/cm, $\Delta T$ = $-$9.2 K for
the $[1\bar{1}0]$ field, to be compared with $-$7.8 mK$\cdot$kV/cm and $-$22.8
K for the $[001]$ field, respectively; however, the relative contributions of
the order parameters share qualitatively similar characteristics, i.e.,
$\omega_{R}$ and $X$ are associated to negative $\Delta T_{\omega_{R}}$ and
$\Delta T_{X}$, respectively, with $\omega_{R}$ contributing the most,
together with the opposite effect from $P$ (positive $\Delta T_{P}$) and
negligible contribution from $\omega_{M}$ (i.e., weak $\Delta
T_{\omega_{M}}$). Note that $\alpha$ and $\Delta T$ are also negative and
increase in magnitude with $T$ within the AFE-based state for the $[110]$
field (Fig. 2g-i), and their magnitude is about $-$2.5 mK$\cdot$kV/cm and
$-$3.9 K, respectively, at the transition point.

However, albeit the qualitative similarities from the sign of $\Delta T_{P}$
and $\Delta T_{X}$ for all three $\mathcal{E}$ directions considered in Figs.
2, $\Delta T_{\omega_{R}}$ becomes positive and $\Delta T_{\omega_{M}}$ shows
large negative contribution only for $\mathcal{E}$ along $[110]$. Such unique
features are associated to the small enhancement of $\omega_{R}$ and large
suppression of $\omega_{M}$ under $[110]$ fields (see Fig. 1e). Further, the
small ECE for $[110]$ field can be well understood by the fact that the
negative $\Delta T_{\omega_{M}}$ and $\Delta T_{X}$ contributions are largely
cancelled by the positive $\Delta T_{P}$ and $\Delta T_{\omega_{R}}$
contributions. It is also worth mentioning that $\Delta T_{P}$ is always
positive and $\Delta T_{X}$ is always negative, as $P$ increases while $X$
decreases upon application of the $\mathcal{E}$ field (see Eq. 4). In other
words, being consistent with Eq. Origin of negative electrocaloric effect in
Pnma-type antiferroelectric perovskites, the applied field makes $P$ more
ordered and $X$ more disordered.

It is interesting to note that the contrast of different ECE behaviors from
various field directions are related to the relative orientation between
$\mathcal{E}$ and $X$; for instance, $[001]$ and $[1\bar{1}0]$ fields are
perpendicular to the antipolar vector, whereas the $[110]$ field is parallel
to it. As detailed in SM Section S5 sm , via the calculated energy landscape,
we show numerically that minimization of the total energy dictates that
increasing $P_{z}$ or $P_{(x,\bar{y})}$ in the Pnma state would strongly
suppress $\omega_{R}$, whereas such suppression is much less for $\omega_{M}$;
and in contrast, the opposite is found for $P_{(x,y)}$, i.e., $\omega_{M}$ is
strongly suppressed but $\omega_{R}$ is less correlated when $P_{(x,y)}$
increases.

In fact, a continuous variation of $\Delta T_{\text{op}}$ can be realized if
the $\mathcal{E}$ direction changes from $[110]$ towards $[001]$. In SM
Section S6 sm , we include results for $\mathcal{E}$ lying along $[111]$,
$[112]$, $[113]$, and $[114]$ directions. While the overall $\Delta T$ remains
negative and $\Delta T_{P}$ is always positive (because $P$ always increases
with the field in the AFE-based phase), $\Delta T_{\omega_{R}}$ becomes more
negative while $\Delta T_{\omega_{M}}$ becomes less negative, as the angle
between $\mathcal{E}$ and $X$ changes from 0 to 90 degrees. Nonetheless, it is
worth to note that the contributions associate to $X$ and $P$ also show
opposite trend with such rotation of the field, i.e., $\Delta T_{P}$ decreases
(less positive) and $\Delta T_{X}$ increases (less negative).

For each field direction, we can further change the field strength to tune the
AFE transition temperature, where the maximum ECE occurs, thus manipulating
the negative ECE. As depicted in Fig. 3, for $\mathcal{E}$ along $[001]$,
$[1\bar{1}0]$ and $[110]$, the maximized negative $\Delta T$ is given by
intermediate field, which allows the AFE-to-FE transition to occur at rather
high temperature. More precisely, if $\mathcal{E}$ is too small, the FE state
is absent; if $\mathcal{E}$ is too large, the AFE-to-FE transition happens at
low temperature and $\Delta T$ is relatively small. Remarkably, $\Delta T$ is
predicted to reach $-$29.0 K if $\mathcal{E}$ is along $[001]$ with a
magnitude of 0.70 MV/cm at 1100 K, therefore being very large compared with
the previously reported value of $-$5.8 K in La-doped Pb(Zr,Ti)O3 Geng  _et
al._ (2015). However, if one is interested in a large negative ECE at room
temperature (RT), a specific field direction and strength should be chosen,
e.g., the largest negative $\Delta T$ of $-$4.7 K is realized at RT with a
$[001]$ field having a magnitude of 1.13 MV/cm, with which the AFE-to-FE
transition is close to RT. This value is comparable to that of PZO-based AFEs
at RT Geng  _et al._ (2015); Vales-Castro  _et al._ (2021).

Finally, let us examine the two previously proposed mechanisms. Regarding the
“dipole-canting” mechanism Geng  _et al._ (2015), in AFE BNFO, for all the
investigated $\mathcal{E}$ directions, if only the dipole degrees of freedom
(i.e., $P$ and $X$) are considered, canting of dipoles indeed occurs in
response to the field (Fig. 1). However, we numerically find that $\Delta
T_{P}$ and $\Delta T_{X}$ have opposite signs, but the response of $P$ to
$\mathcal{E}$ field is much greater than $X$ for all field directions. The
corresponding change of temperature ($\Delta T_{P}+\Delta T_{X}$) is thus
always positive, instead of negative, which therefore can not explain the
overall negative ECE. This observation strongly suggests that explicit
consideration of octahedral degrees of freedom is mandatory to correctly
account for the (negative) ECE in AFEs. On the other hand, the perturbative
approach solely based on polarization (as functions of $\mathcal{E}$ field and
temperature) Graf and Íñiguez (2021) agrees reasonably well with our
predictions at low fields (see SM Section S7 sm ), indicating that the Maxwell
relation $\frac{\partial S}{\partial\mathcal{E}}\rvert_{T}=\frac{\partial
P}{\partial T}\rvert_{\mathcal{E}}$ holds for negative ECE systems with
contributions from multiple degrees of freedom.

In summary, our quantitative method show that, as octahedral tiltings are
often primary order parameters in AFE perovskites, they have to be taken
explicitly into account to construe the negative electrocaloric effect, and
mere inclusion of dipolar degrees of freedom is insufficient. We also predict
that AFE BNFO solid solutions are very promising to achieve giant temperature
decrease near the AFE transition upon application of an electric field. Our
atomistic mechanism also suggests that enhanced negative ECE may be possible
upon freezing or suppression of order parameter(s) with positive $\Delta
T_{\text{op}}$, for instance via interfacial proximity effect.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
under Grant No. 12074277, Natural Science Foundation of Jiangsu Province
(BK20201404), the startup fund from Soochow University and the support from
Priority Academic Program Development (PAPD) of Jiangsu Higher Education
Institutions. L.B. thanks the Office of Naval Research for the support under
Grant No. N00014-21-1-2086.

## References

  * Correlia and Zhang (2014) T. Correlia and Q. Zhang,  _Electrocaloric materials: new generation of coolers_ (Springer, Berlin, 2014).
  * Fähler  _et al._ (2012) S. Fähler, U. K. Rößler, O. Kastner, J. Eckert, G. Eggeler, H. Emmerich, P. Entel, S. Müller, E. Quandt, and K. Albe, Adv. Eng. Mater. 14, 10 (2012).
  * Moya  _et al._ (2014) X. Moya, S. Kar-Narayan, and N. D. Mathur, Nat. Mater. 13, 439 (2014).
  * Pirc  _et al._ (2014) R. Pirc, B. Rožič, J. Koruza, B. Malič, and Z. Kutnjak, Europhys. Lett. 107, 17002 (2014).
  * Geng  _et al._ (2015) W. Geng, Y. Liu, X. Meng, L. Bellaiche, J. F. Scott, B. Dkhil, and A. Jiang, Adv. Mater. 27, 3165 (2015).
  * Vales-Castro  _et al._ (2021) P. Vales-Castro, R. Faye, M. Vellvehi, Y. Nouchokgwe, X. Perpiñà, J. M. Caicedo, X. Jordà, K. Roleder, D. Kajewski, A. Perez-Tomas, _et al._ , Phys. Rev. B 103, 054112 (2021).
  * Gottschall  _et al._ (2018) T. Gottschall, A. Gracia-Condal, M. Fries, A. Taubel, L. Pfeuffer, L. Manosa, A. Planes, K. P. Skokov, and O. Gutfleisch, Nat. Mater. 17, 929 (2018).
  * Marathe  _et al._ (2017) M. Marathe, D. Renggli, M. Sanlialp, M. O. Karabasov, V. V. Shvartsman, D. C. Lupascu, A. Grünebohm, and C. Ederer, Phys. Rev. B 96, 014102 (2017).
  * Grünebohm  _et al._ (2018) A. Grünebohm, Y.-B. Ma, M. Marathe, B.-X. Xu, K. Albe, C. Kalcher, K.-C. Meyer, V. V. Shvartsman, D. C. Lupascu, and C. Ederer, Energy Technol. 6, 1491 (2018).
  * Graf and Íñiguez (2021) M. Graf and J. Íñiguez, Commun. Mater 2, 1 (2021).
  * Milesi-Brault  _et al._ (2020) C. Milesi-Brault, C. Toulouse, E. Constable, H. Aramberri, V. Simonet, S. De Brion, H. Berger, L. Paolasini, A. Bosak, J. Íñiguez, _et al._ , Phys. Rev. Lett. 124, 097603 (2020).
  * Xu  _et al._ (2019) B. Xu, O. Hellman, and L. Bellaiche, Phys. Rev. B 100, 020102 (2019).
  * Bellaiche and Íñiguez (2013) L. Bellaiche and J. Íñiguez, Phys. Rev. B 88, 014104 (2013).
  * Planes  _et al._ (2014) A. Planes, T. Castan, and A. Saxena, Philos. Mag. 94, 1893 (2014).
  * Zhang  _et al._ (2011) J. Zhang, A. Heitmann, S. Alpay, and G. Rossetti Jr, Integr. Ferroelectr. 125, 168 (2011).
  * Ma  _et al._ (2016) Y.-B. Ma, N. Novak, J. Koruza, T. Yang, K. Albe, and B.-X. Xu, Phys. Rev. B 94, 100104 (2016).
  * Edström and Ederer (2020) A. Edström and C. Ederer, Phys. Rev. Lett. 124, 167201 (2020).
  * Karimi  _et al._ (2009a) S. Karimi, I. Reaney, Y. Han, J. Pokorny, and I. Sterianou, J. Mater. Sci. 44, 5102 (2009a).
  * Xu  _et al._ (2015a) B. Xu, D. Wang, J. Íñiguez, and L. Bellaiche, Adv. Funct. Mater. 25, 552 (2015a).
  * Karimi  _et al._ (2009b) S. Karimi, I. Reaney, I. Levin, and I. Sterianou, Appl. Phys. Lett. 94, 112903 (2009b).
  * Kan  _et al._ (2010) D. Kan, L. Pálová, V. Anbusathaiah, C. J. Cheng, S. Fujino, V. Nagarajan, K. M. Rabe, and I. Takeuchi, Adv. Funct. Mater. 20, 1108 (2010).
  * Xu  _et al._ (2017) B. Xu, J. Íñiguez, and L. Bellaiche, Nat. Commun. 8, 15682 (2017).
  * (23) See Supplemental Material for (1) Computational details, (2) Temperature dependencies of the order parameters for various fields, (3) DFT-calculated specific heat, (4) Determination of parameters in the Landau free-energy model, (5) Correlation between polarization and octahedral tiltings, (6) Additional $\mathcal{E}$-field directions from $[110]$ to $[001]$ and (7) Comparison with a perturbative approach, which additionally includes Refs. Prosandeev  _et al._ (2013); Albrecht  _et al._ (2010); Kornev  _et al._ (2007); Lisenkov  _et al._ (2009); Zhong  _et al._ (1994); Zhong and Vanderbilt (1995); Kornev  _et al._ (2006); Daumont  _et al._ (2012); Kresse and Furthmüller (1996a, b); Blöchl (1994); Perdew  _et al._ (1996); Baroni  _et al._ (1987, 2001); Gonze and Vigneron (1989); Gonze and Lee (1997); Togo and Tanaka (2015); Íñiguez  _et al._ (2014); Wojdeł and Íñiguez (2014).
  * (24) BNFO is multiferroic and in principle the magnetic moments can also respond to the applied electric field. Here it is not considered since we focus on the structural order parameters for a general AFE and numerically the magnetic contribution to $\alpha$ is found to the small (see Ref. Jiang  _et al._ (2021b)).
  * Xu  _et al._ (2015b) B. Xu, D. Wang, H. J. Zhao, J. Íñiguez, X. M. Chen, and L. Bellaiche, Adv. Funct. Mater. 25, 3626 (2015b).
  * (26) Note that by AFE-based state, we mean the AFE state with all its associated order parameters but on top of which polarization develops along the field direction.
  * Jiang  _et al._ (2017) Z. Jiang, S. Prokhorenko, S. Prosandeev, Y. Nahas, D. Wang, J. Íñiguez, E. Defay, and L. Bellaiche, Phys. Rev. B 96, 014114 (2017).
  * Jiang  _et al._ (2018) Z. Jiang, Y. Nahas, S. Prokhorenko, S. Prosandeev, D. Wang, J. Íñiguez, and L. Bellaiche, Phys. Rev. B 97, 104110 (2018).
  * Bin-Omran  _et al._ (2016) S. Bin-Omran, I. A. Kornev, and L. Bellaiche, Phys. Rev. B 93, 014104 (2016).
  * Jiang  _et al._ (2021a) Z. Jiang, B. Xu, S. Prosandeev, Y. Nahas, S. Prokhorenko, J. Íñiguez, and L. Bellaiche, Phys. Rev. B 103, L100102 (2021a).
  * Jiang  _et al._ (2021b) Z. Jiang, B. Xu, S. Prosandeev, Y. Nahas, S. Prokhorenko, J. Íñiguez, and L. Bellaiche, Phys. Rev. B 103, L100102 (2021b).
  * Zhao  _et al._ (2014) H. J. Zhao, J. Íñiguez, W. Ren, X. M. Chen, and L. Bellaiche, Phys. Rev. B 89, 174101 (2014).
  * Kutnjak  _et al._ (1999) Z. Kutnjak, B. Rožič, and R. Pirc, Wiley encyclopedia of electrical and electronics engineering , 1 (1999).
  * (34) op with large coeﬀicient $A_{\text{op}}$ is exptected to have strongly $T$-dependent change of the corresponding vibrational entropy.
  * Prosandeev  _et al._ (2013) S. Prosandeev, D. Wang, W. Ren, J. Íñiguez, and L. Bellaiche, Adv. Funct. Mater. 23, 234 (2013).
  * Albrecht  _et al._ (2010) D. Albrecht, S. Lisenkov, W. Ren, D. Rahmedov, I. A. Kornev, and L. Bellaiche, Phys. Rev. B 81, 140401 (2010).
  * Kornev  _et al._ (2007) I. A. Kornev, S. Lisenkov, R. Haumont, B. Dkhil, and L. Bellaiche, Phys. Rev. Lett. 99, 227602 (2007).
  * Lisenkov  _et al._ (2009) S. Lisenkov, D. Rahmedov, and L. Bellaiche, Phys. Rev. Lett. 103, 047204 (2009).
  * Zhong  _et al._ (1994) W. Zhong, D. Vanderbilt, and K. M. Rabe, Phys. Rev. Lett. 73, 1861 (1994).
  * Zhong and Vanderbilt (1995) W. Zhong and D. Vanderbilt, Phys. Rev. Lett. 74, 2587 (1995).
  * Kornev  _et al._ (2006) I. A. Kornev, L. Bellaiche, P.-E. Janolin, B. Dkhil, and E. Suard, Phys. Rev. Lett. 97, 157601 (2006).
  * Daumont  _et al._ (2012) C. Daumont, W. Ren, I. Infante, S. Lisenkov, J. Allibe, C. Carrétéro, S. Fusil, E. Jacquet, T. Bouvet, F. Bouamrane, _et al._ , J. Phys.: Condens. Matter 24, 162202 (2012).
  * Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996a).
  * Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Comput. Mater. Sci 6, 15 (1996b).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Baroni  _et al._ (1987) S. Baroni, P. Giannozzi, and A. Testa, Phys. Rev. Lett. 58, 1861 (1987).
  * Baroni  _et al._ (2001) S. Baroni, S. De Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  * Gonze and Vigneron (1989) X. Gonze and J.-P. Vigneron, Phys. Rev. B 39, 13120 (1989).
  * Gonze and Lee (1997) X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997).
  * Togo and Tanaka (2015) A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).
  * Íñiguez  _et al._ (2014) J. Íñiguez, M. Stengel, S. Prosandeev, and L. Bellaiche, Phys. Rev. B 90, 220103 (2014).
  * Wojdeł and Íñiguez (2014) J. C. Wojdeł and J. Íñiguez, Phys. Rev. B 90, 014105 (2014).