# Schwinger particle production: rapid switch off of the external field versus
dynamical assistance

I. A. Aleksandrov  Department of Physics, Saint Petersburg State University,
Universitetskaya Naberezhnaya 7/9, Saint Petersburg 199034, Russia  Ioffe
Institute, Politekhnicheskaya Street 26, Saint Petersburg 194021, Russia  D.
G. Sevostyanov  Department of Physics, Saint Petersburg State University,
Universitetskaya Naberezhnaya 7/9, Saint Petersburg 199034, Russia  V. M.
Shabaev  Department of Physics, Saint Petersburg State University,
Universitetskaya Naberezhnaya 7/9, Saint Petersburg 199034, Russia

###### Abstract

We consider the process of electron-positron pair production in the presence
of strong electric backgrounds being rapidly switched on and off and examine
the total particle yield. For sufficiently sharp field profiles, the particle
number can be substantially enhanced. It is demonstrated that this enhancement
is quite similar to the phenomenon of dynamical assistance by a weak high-
frequency field superimposed on a strong background. Both these mechanisms are
analyzed by means of exact numerical computations and the locally-constant
field approximation. We identify the time scale of the switching profile
leading to the pair-production enhancement and argue that the particle yield
is highly unlikely to be increased by shaping the switch off of realistic
laser pulses. On the other hand, we confirm that it is feasible to observe the
dynamically assisted Schwinger effect by adding a rapidly oscillating field to
a strong electric background.

##  I Introduction

One of the most striking phenomena predicted by quantum electrodynamics (QED)
is the process of electron-positron pair production from vacuum in intense
electromagnetic fields sauter_1931 ; heisenberg_euler ; schwinger_1951 . In
the regime of strong coupling, where the external background cannot be treated
within perturbation theory (PT), this effect (Schwinger mechanism) still
awaits experimental observations. By continually increasing the pulse
intensity that can be achieved at the modern laser facilities,
experimentalists have stimulated great interest in strong-field QED over the
past decades, so it is also extremely important for theoretical groups to
explore the possibilities for enhancing the particle yield and propose most
favorable practical scenarios. Here we will focus on two general ideas.

First, it is known that evaluating the adiabatic number of particles at
intermediate times, when the external field is on, one may obtain values that
are orders of magnitude larger than the residual (asymptotic) pair yield (see,
e.g., Refs. fedotov_prd_2011 ; blaschke_prd_2013 ; blinne_prd_2014 ; zahn_2015
; otto_jpp_2016 ). Although it is the latter which is to be observed in the
experiment, one can attempt to maximize it by choosing a more advantageous
field profile. Recently, it was suggested that by properly switching the
external field off, one can indeed make the large adiabatic particle number
observable ilderton_prd_2022 . The switching profile should be sufficiently
sharp, and its location should be accurately adjusted. This issue was
addressed in Ref. ilderton_prd_2022  based on the analysis of the electron
(positron) number density evaluated nonperturbatively and the total particle
yield computed within PT. In the present study, we examine the total number of
pairs without perturbative expansions.

Second, in Ref. schuetzhold_prl_2008  it was demonstrated that the number of
pairs created from vacuum can be drastically enhanced by superimposing a weak
rapidly oscillating pulse on a strong quasistatic external background.
Although the former has a small amplitude, it provides a high-frequency
harmonic which efficiently stimulates the pair production process. As this so-
called dynamically assisted Schwinger effect (DASE) appears to be very
promising from the experimental viewpoint, it has been explored in numerous
investigations (see, e.g., Refs. orthaber_2011 ; fey ; kohlfuerst_prd_2013 ;
linder_prd_2015 ; akal_prd_2014 ; hebenstreit_plb_2014 ; otto_plb_2015 ;
panferov_epjd_2016 ; aleksandrov_prd_2018 ; otto_epja_2018 ;
sitiwaldi_plb_2018 ; torgrimsson_prd_2018 ; huang_prd_2019 ; chavez_prd_2019 ;
olugh_plb_2020 ; taya_prr_2020 ; li_prd_2021 ).

Here we discuss how the two aforementioned phenomena enhancing the particle
yield are related. The adiabatic number of $e^{+}e^{-}$ pairs evaluated at
intermediate times is the number of physical pairs which would be measured if
the external field were suddenly turned off at that time ilderton_prd_2022 .
This suggests that it is the high-frequency mode induced by the rapid
switching profile which multiplies the particle yield. As will be demonstrated
in this paper, this mechanism is closely related to the DASE. We will
scrutinize this relationship and identify the common patterns of the two
effects. In order to provide a connection with the experimental setups, we
will find out which values of the field parameters lead to the enhancement of
the particle number and assess the feasibility of the corresponding practical
scenarios.

To calculate the number of $e^{+}e^{-}$ pairs, we will first employ our
numerical procedures based on the Furry-picture quantization in momentum space
aleksandrov_prd_2016 ; aleksandrov_prd_2017_2 ; aleksandrov_kohlfuerst  and
quantum kinetic equations aleksandrov_qke_2020 ; sevostyanov_prd_2021 . These
nonperturbative approaches neglect the radiative corrections but fully
incorporate the interaction between the quantized Dirac field and the
classical background. For simplicity, we assume the external field to depend
solely on time (the corresponding temporal profile is treated without any
additional approximations). Second, it turns out that great insights can be
obtained by means of the so-called locally-constant field approximation (LCFA)
bunkin_tugov ; narozhny_pla_2004 ; dunne_prd_2006 ; hebenstreit_prdr_2008 ;
bulanov_prl_2010 ; gavrilov_prd_2017 ; sevostyanov_prd_2021 ;
aleksandrov_prd_2019_1 . The validity of this approach to computing the total
particle yield was recently examined in Ref. sevostyanov_prd_2021  (see also
Ref. aleksandrov_prd_2019_1 ). Finally, to illustrate the main patterns of the
enhancement mechanisms, we will perform the analysis of the Fourier transforms
of the external electric backgrounds.

The paper is organized as follows. In Sec. II we describe the field
configurations to be investigated. In Sec. III we present the general
expressions obtained within the LCFA. Section IV contains our analysis of the
enhancement mechanisms. The experimental prospects regarding particle
enhancement are discussed in Sec. V. We conclude in Sec. VI.

Throughout the text, we employ the units where $\hbar=c=1$ and the electric
charge $e<0$ obeys $\alpha=e^{2}/(4\pi)$ ($\alpha\approx 1/137$ is the fine-
structure constant).

##  II Setup. External electric field

We choose the classical external background in the form of a linearly
polarized electric field $E(t)$ vanishing at $t\to\pm\infty$. The main
patterns of the DASE will be discussed in the case of a combination of two
Sauter pulses:

$E_{\text{I}}(t)=\frac{E}{\cosh^{2}(t/T)}+\frac{\varepsilon}{\cosh^{2}(t/\tau)}.$
(1)

The strong pulse parameters obey $|e|ET\gtrsim m$, so it cannot be treated
within PT. Varying the weak pulse duration $\tau$, we will identify the onset
of dynamical assistance and thoroughly examine the $\tau$ dependence of the
total number of pairs.

The role of the switching profile and possible enhancement due to the
switching effects will be analyzed for the following external field
configuration:

$E_{\text{II}}(t)=\begin{cases}E&\text{if}\leavevmode\nobreak\
\leavevmode\nobreak\ |t|\leq T/2,\\\
E\cos^{2}\big{[}\frac{\pi}{2\tau}(|t|-T/2)\big{]}&\text{if}\leavevmode\nobreak\
\leavevmode\nobreak\ T/2<|t|\leq T/2+\tau,\\\ 0&\text{otherwise}.\end{cases}$
(2)

Here $T$ is the temporal extent of the plateau region, where the external
field is constant, and $\tau$ governs the duration of the smooth switch-on and
switch-off parts.

Since the external field is assumed to be spatially homogeneous, the total
number of pairs will be proportional to the volume of the system, $N\sim V$,
so in what follows we will present the results in terms of the pair yield per
unit volume $\nu=N/V$. The corresponding numbers will also carry the
subscripts “I” and “II”.

##  III Locally-constant field approximation

The idea of the LCFA is to employ a closed-form expression for the total
number of pairs produced per unit volume and time in the presence of a
constant electromagnetic field. This reads nikishov_constant

$\frac{dN}{dtd\boldsymbol{x}}\,[\mathcal{E},\mathcal{H}]=\frac{e^{2}\mathcal{E}\mathcal{H}}{4\pi^{2}}\coth{\frac{\pi\mathcal{H}}{\mathcal{E}}}\,\mathrm{exp}\bigg{(}\\!-\frac{\pi
m^{2}}{|e|\mathcal{E}}\bigg{)},$ (3)

where

$\displaystyle\mathcal{E}$ $\displaystyle=$
$\displaystyle\sqrt{\sqrt{\mathcal{F}^{2}+\mathcal{G}^{2}}+\mathcal{F}},$ (4)
$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\sqrt{\sqrt{\mathcal{F}^{2}+\mathcal{G}^{2}}-\mathcal{F}}.$ (5)

Here the invariant quantities $\mathcal{F}$ and $\mathcal{G}$ are defined in
terms of the electric and magnetic field vectors via
$\mathcal{F}=(\boldsymbol{E}^{2}-\boldsymbol{H}^{2})/2$ and
$\mathcal{G}=\boldsymbol{E}\cdot\boldsymbol{H}$. Accordingly, $\mathcal{E}$
and $\mathcal{H}$ are the magnitudes of the electric and magnetic field
components, respectively, in the frame where they are parallel to each other.
The Schwinger critical value of the electric field strength is
$E_{\text{c}}=m^{2}/|e|$.

Within the LCFA, one plugs the actual spatiotemporal dependence of the
external inhomogeneous field into Eq. (3) and integrates it then over $t$ and
$\boldsymbol{x}$ (see Refs. bunkin_tugov ; narozhny_pla_2004 ; dunne_prd_2006
; hebenstreit_prdr_2008 ; bulanov_prl_2010 ; gavrilov_prd_2017 ;
sevostyanov_prd_2021 ; aleksandrov_prd_2019_1 ). In the case of a spatially
uniform electric field $E(t)$, this procedure yields the following expression
for the total number of pairs per unit volume:

$\nu^{\text{(LCFA)}}=\frac{e^{2}}{4\pi^{3}}\int\\!dt\,E^{2}(t)\,\mathrm{exp}\left(\\!-\frac{\pi
m^{2}}{|eE(t)|}\right).$ (6)

In what follows, we will perform the calculations using the specific forms (1)
and (2). Although in the case of two Sauter pulses (1), the LCFA prediction
(6) is nonlinear with respect to the field parameters, employing the function
$E_{\text{II}}(t)$ leads to the following result:

$\nu^{\text{(LCFA)}}_{\text{II}}=[T+f(E/E_{\text{c}})\tau]g(E),$ (7)

where

$g(E)=\frac{e^{2}E^{2}}{4\pi^{3}}\mathrm{exp}\left(\\!-\frac{\pi
m^{2}}{|eE|}\right),$ (8)

which represents the number of pairs created in a constant electric field $E$
per unit volume and time. We also introduced

$f(\xi)=\frac{4}{\pi}\int\limits_{0}^{\pi/2}dx\cos^{4}x\
\mathrm{exp}\left(\\!-\frac{\pi}{\xi}\tan^{2}x\right).$ (9)

Accordingly, the particle number evaluated within the LCFA linearly depends on
the parameters $T$ and $\tau$. We also note that
$f(\xi)\approx(2/\pi)\sqrt{\xi}$ for $\xi\ll 1$ (see Appendix A).

Equation (7) contains $g(E)$, which reflects the fact that a constant electric
background generates $e^{+}e^{-}$ pairs in the tunneling (nonperturbative)
regime. The LCFA is valid only in the domain of sufficiently large field
amplitudes and pulse durations. The approximate predictions and the results of
exact numerical computations will be discussed next.

##  IV Results. Enhancement mechanisms

Here we will first evaluate the total number of pairs produced for the two
scenarios described above as a function of $\tau$, which represents there the
weak-pulse duration and the duration of the switch-on and switch-off parts,
respectively. To gain additional insights into the enhancement mechanisms, we
will also analyze the Fourier spectra of the corresponding electric
backgrounds (1) and (2).

###  IV.1 Total particle yield. Dynamical assistance

In Fig. 1 we present the total number of pairs as a function of the weak-pulse
duration $\tau$ in the case of the external field (1). This example was
obtained for $E=0.2E_{\text{c}}$, $T=20m^{-1}$, and
$\varepsilon=0.05E_{\text{c}}$. The dashed line corresponds to the LCFA
prediction (6). The results of our exact numerical computations are
represented by the grey solid line. For sufficiently large $\tau$, the exact
values of $\nu_{\text{I}}$ almost coincide with
$\nu_{\text{I}}^{\text{(LCFA)}}$, but there is still a subtle discrepancy even
for larger $\tau$ since the LCFA is not completely justified even if the weak
pulse is absent: for the strong-pulse parameters chosen here and $\tau=0$, the
LCFA yields $4.0\times 10^{-10}m^{3}$, while the exact number of pairs
produced by the strong pulse amounts to $5.4\times 10^{-10}m^{3}$ (see the
horizontal light-blue line in Fig. 1). According to Ref. sevostyanov_prd_2021
, in the case of an individual Sauter pulse, the LCFA accurately predicts the
particle yield if the parameter $|eE|^{3/2}T$ is much larger than $m^{2}$, but
in our case it amounts to $|eE|^{3/2}T\approx 1.8m^{2}$. Note that for large
$\tau$ the LCFA is quite accurate in the case of two pulses although it
completely fails to describe the pair production process in the field of the
weak pulse alone (e.g., for $\tau=10m^{-1}$ the LCFA underestimates the
particle yield by approximately 14 orders of magnitude).

Figure 1: Total number of $e^{+}e^{-}$ pairs per unit volume created by the
external field (1) as a function of the weak-pulse duration $\tau$
($E=0.2E_{\text{c}}$, $T=20m^{-1}$, $\varepsilon=0.05E_{\text{c}}$). The
results are obtained within the LCFA and by means of our exact numerical
techniques. The dash-dotted line and thick dots correspond to the particle
yield produced by the weak pulse only [second term in Eq. (1)].

With decreasing $\tau$, the exact number of pairs grows rapidly, while the
LCFA prediction decreases. The latter point is evident already from the
general LCFA expression (6) as for smaller $\tau$ the local values $E(t)$
decrease. Whereas the LCFA incorporates only the tunneling pair-production
mechanism, it is the dynamical effect which enhances the exact particle
number. It turns out that the number of pairs in the combination of two Sauter
pulses can substantially exceed each of the two individual contributions. This
result constitutes the DASE.

Nevertheless, once the weak-pulse duration $\tau$ becomes too small, the
individual weak pulse may produce almost the same number of pairs as the two
pulses generate together. In this domain, the weak pulse does not assist the
strong one but rather governs the pair-production process unaided. In order to
simplify the analysis and benchmark the nonperturbative methods, we also
employ PT in the case of the individual weak pulse (dash-dotted line in Fig.
1). The total number of pairs per unit volume to the second order in the
external-field amplitude $\varepsilon$ reads (see Appendix B)

$\nu_{\text{weak}}^{\text{(PT)}}=\frac{e^{2}\varepsilon^{2}\tau^{4}}{4\pi}\int
d\boldsymbol{p}\,\frac{m^{2}+\boldsymbol{p}_{\perp}^{2}}{p_{0}^{2}}\,\frac{1}{\sinh^{2}(\pi\tau
p_{0})},$ (10)

where $\boldsymbol{p}_{\perp}$ is the transverse momentum component and
$p_{0}=\sqrt{m^{2}+\boldsymbol{p}^{2}}$. The expression (10) is valid only for
large values of the Keldysh parameter $\gamma_{\varepsilon}\equiv
m/|e\varepsilon\tau|\gg 1$. For instance, for $\tau=5m^{-1}$ the PT approach
underestimates the exact yield by a factor of $\sim 20$ and for
$\tau=10m^{-1}$ the PT result is already about ten orders of magnitude
smaller.

Finally, for very small $\tau$ ($m\tau\ll 1$), one can expand Eq. (10) in
powers of $\tau$ and obtain the following leading-order contribution (see
Appendix B):

$\nu_{\text{weak}}^{\text{(PT)}}\approx\frac{1}{9\pi}\,\frac{\varepsilon^{2}}{E_{\text{c}}^{2}}\,(m\tau)m^{3}.$
(11)

Accordingly, in the regime $m\tau\ll 1$ the number of pairs produced by the
weak field linearly depends on the pulse duration and vanishes at $\tau=0$.

Figure 2: Enhancement factor defined in Eq. (12) as a function of $\tau$ for
$E=0.2E_{\text{c}}$, $T=20m^{-1}$, and $\varepsilon=0.05E_{\text{c}}$.

In order to quantitatively estimate the efficiency of dynamical assistance,
one can first calculate the ratio
$\nu_{\text{I}}/(\nu_{\text{strong}}+\nu_{\text{weak}})$, which shows the
relative increase in the particle yield arising when the strong and weak
pulses are combined. However, this quantity becomes irrelevant in the domain
of large $\tau$ since in this case the number of pairs can be evaluated simply
by means of the LCFA, i.e., by taking into account the tunneling mechanism
[note that in the tunneling regime for $T=\tau$, the ratio is always large due
to the highly nonlinear dependence of the particle number on the field
amplitude, see Eq. (6)]. Accordingly, we suggest that one relies on the
enhancement factor defined via

$\mathcal{K}=\mathrm{min}\leavevmode\nobreak\
\Bigg{\\{}\frac{\nu_{\text{I}}}{\nu_{\text{strong}}+\nu_{\text{weak}}},\leavevmode\nobreak\
\frac{\nu_{\text{I}}}{\nu^{\text{(LCFA)}}_{\text{I}}}\Bigg{\\}}.$ (12)

In Fig. 2 we depict the enhancement factor $\mathcal{K}$ as a function of
$\tau$. In our example, the presence of the short weak pulse can enhance the
total number of particles by one order of magnitude (in principle, by varying
the field parameters, one can achieve a stronger enhancement, see, e.g., Ref.
aleksandrov_prd_2018 ). The most important point for the assessment of the
experimental prospects is that the onset of the DASE corresponds to
sufficiently small values of $\tau$: for $\varepsilon\ll E$, the combined
Keldysh parameter $\gamma_{\text{c}}\equiv m/|eE\tau|$ should be sufficiently
large, $\gamma_{\text{c}}\gtrsim 1$ schuetzhold_prl_2008 . The peak in Fig. 2
is located at $\tau_{*}\approx 2.3m^{-1}$. Including the weak-pulse amplitude
$\varepsilon$, which in our case not that small compared to $E$, we obtain
$m/|e(E+\varepsilon)\tau_{*}|\approx 1.7$. Let us underline that the condition
$\gamma_{\text{c}}\gtrsim 1$ reflects the onset of the DASE, but one has to
make sure that the nonperturbative nature of the pair production process is
preserved. Since the strong (weak) pulse operates in the tunneling
(perturbative) regime, the weak pulse should only assist the particle
production mechanism.

Dynamical assistance can be viewed as a process of additional photon
absorption galvanizing the tunneling mechanism. For this to happen, the energy
spectrum of the weak pulse should be sufficiently broad, so that it contains
large $\omega$, i.e., $\tau$ should be sufficiently small (note that the
energy spectrum of a Sauter pulse is centered at $\omega=0$). This will be
examined in more detail in Sec. IV.3. Instead of two Sauter pulses, one can
also analyze a combination of two oscillating pulses with a low and high
frequency, respectively (see, e.g., Refs. linder_prd_2015 ; akal_prd_2014 ;
hebenstreit_plb_2014 ; otto_plb_2015 ; panferov_epjd_2016 ;
aleksandrov_prd_2018 ; otto_epja_2018 ; sitiwaldi_plb_2018 ; olugh_plb_2020 ;
li_prd_2021 ). Although this configuration appears to be more relevant from
the experimental point of view, the main patterns of the DASE are the same.
For instance, the number of particles as a function of the weak-pulse
frequency $\omega$ has also two maxima as does the $\tau$-dependence revealed
in Fig. 1 once one maps $\omega\sim 1/\tau$ kohlfuerst_prd_2013 ;
kohlfuerst_arxiv_2022 . In the case of oscillating pulses, the onset of the
DASE is governed by the combined Keldysh parameter
$\gamma_{\text{c}}=m\omega/|eE|$.

Figure 3: Total number of $e^{+}e^{-}$ pairs in the case of the external field
(2) as a function of $\tau$ for various values of $T$ ($E=0.5E_{\text{c}}$).
The dashed lines correspond to the LCFA predictions. The total number of pairs
for $T=0$ is also evaluated within PT (dash-dotted line).

###  IV.2 Total particle yield. Switch-off effects

We now turn to the study of the enhancement mechanism due to a sharp switch
off of the external field. The main idea is to make use of the fact that the
adiabatic particle number at intermediate times can be orders of magnitude
larger than the final particle yield. To maximize the number of pairs, one can
rapidly turn the field off, so that the large adiabatic number transforms to a
physically observable quantity ilderton_prd_2022 . This enhancement mechanism
can be efficient only if the corresponding smooth switch off is sufficiently
rapid. Here our main goal is to analyze the time scale of the switching
profile required for the pair-production enhancement and to assess then the
possible experimental implementations of this mechanism. It turns out that one
can substantially deepen understanding of this scenario by comparing it with
the DASE. As will be seen below, the two phenomena are closely related.

We will examine the external-field configuration given in Eq. (2), which
represents a rectangular-like pulse with smooth switch-on and switch-off parts
of duration $\tau$ each. First, we note that in the case of a sharp
rectangular profile with parameters $T$ and $\tau=0$, the adiabatic number of
particles at time instant $t$ ($|t|<T/2$) obviously coincides with the final
pair yield for the pulse with duration $T^{\prime}=t+T/2$. The adiabatic
particle number is equally nonphysical as the sudden switch on and off of the
rectangular pulse for $\tau=0$. To make the field configuration more
realistic, we consider nonzero values of $\tau$.

In Fig. 3 we present the total number of pairs produced as a function of
$\tau$ for $E=0.5E_{\text{c}}$ and various durations of the plateau region
$T$. Here we observe that the $\tau$ dependence is very similar to that
displayed in Fig. 1. According to Eq. (7), the LCFA prediction is linear and
decreases with decreasing $\tau$ (the dashed lines have a log-like shape
although they may seem straight in the graph). Within the LCFA, for $\tau=0$
we obtain the constant-field result (8) multiplied by the rectangular-pulse
duration $T$. On the other hand, the exact number of pairs proves to be
considerably larger than the LCFA estimates in the domain of sufficiently
small $\tau$. For large values of $\tau$, the exact particle yield can be
accurately obtained within the LCFA. In the case of $T=0$, the field profile
turns to a $\cos^{2}$ pulse, whose shape looks similar to that of the Sauter
pulse. In this case, we also employ PT and reproduce the exact result in the
region $|eE|\tau\ll m$, i.e., the dash-dotted line in Fig. 3 coincides with
the solid light-blue one for small $\tau$. (We note that the PT approach is
much less accurate here than in the context of the DASE because the amplitude
$E$ of the rectangular pulse is ten times larger than $\varepsilon$ in Fig.
1.) If the pulse duration $\tau$ becomes extremely small, $m\tau\ll 1$, the PT
result exhibits a linear behavior (see Appendix C),

$\nu_{T=0}^{\text{(PT)}}\approx\frac{1}{16\pi}\,\frac{E^{2}}{E_{\text{c}}^{2}}\,(m\tau)m^{3}.$
(13)

The $\tau$ dependences for nonzero $T$ in Fig. 3 have obviously nonzero values
at $\tau=0$ and cannot be described by means of PT since the parameter
$|eE|T/m$ is not small.

Our primary goal is to identify the onset of the enhancement mechanism under
consideration. First, let us note that nonperturbative pair production takes
place in the domain $|eE|T\gtrsim m$. It is convenient to introduce
$T_{0}=m/|eE|$, which represents a “formation time”: within this time scale, a
virtual $e^{+}e^{-}$ pair can gain enough energy to become real. We are
interested in the regime $T\gtrsim T_{0}$. According to Ref.
sevostyanov_prd_2021 , the switch-on and switch-off parts can be accurately
treated within the LCFA if $|eE|^{3/2}\tau\gg m^{2}$. In this domain, the
particle yield cannot be enhanced since the LCFA predictions monotonically
depend on $\tau$. Suppose now that for given $\tau$ the LCFA is applicable and
consider the contributions from the switch-on and -off intervals to the
particle yield according to Eq. (6). As was stated above [see Eq. (7)], they
can be estimated as follows:

$\nu^{\text{(LCFA)}}_{\tau}=g(E)f(E/E_{\text{c}})\tau\approx\frac{2}{\pi}g(E)\sqrt{\frac{E}{E_{\text{c}}}}\tau.$
(14)

If one decreases $\tau$, the LCFA becomes inapplicable at $\tau\sim
m^{2}/|eE|^{3/2}=(E_{\text{c}}/E)^{3/2}m^{-1}$. For this value of $\tau$, Eq.
(14) yields $\nu^{\text{(LCFA)}}_{\tau}\sim g(E)T_{0}$. The LCFA term arising
from the plateau region amounts to $g(E)T$, so the relative contribution of
the switch-on and -off parts is $T_{0}/T\lesssim 1$. By increasing $T$, one
can further decrease the relative number of particles produced due to the
switching intervals at the threshold of the LCFA applicability. It means that
for larger values of $T$, the LCFA can provide quite accurate results for
smaller $\tau$, i.e., even if the condition $|eE|^{3/2}\tau\gg m^{2}$ no
longer holds. This issue is clearly demonstrated in Fig. 3: larger values of
$T$ improves the accuracy of the LCFA. These results were also confirmed by
performing calculations for other values of the field amplitude $E$.

In what follows, we will use the inequality $|eE|^{3/2}\tau\lesssim m^{2}$ as
a necessary condition for the pair-production enhancement. This represents a
sufficient condition only for $T=0$, but instead of specifying $T>0$, we will
obtain an upper bound on $\tau$ necessary for the enhancement of the particle
yield. It can be easily seen that the DASE criterion $\gamma_{\text{c}}\gtrsim
1$ governing the onset of the pair-production enhancement is stronger than
$|eE|^{3/2}\tau\lesssim m^{2}$ as it should be (the LCFA is not able to
capture the DASE).

We also point out that in Ref. ilderton_prd_2022  it was suggested that the
final particle number can be enhanced by rapidly switching off a Sauter pulse
at the vicinity of its maximum. In Ref. ilderton_prd_2022  the total number of
pairs was evaluated only be means of PT. We also performed exact
nonperturbative calculations using the same switching profiles as in Ref.
ilderton_prd_2022  and confirmed the main findings reported there. Moreover,
we found that in order to examine the main qualitative patterns of the
enhancement mechanism and provide the necessary quantitative estimates, it
suffices to explore the simple field configuration (2) although by accurately
shaping the external pulse, one may further increase the particle yield.

Our findings can also be illustrated by analyzing the Fourier transforms of
the corresponding external field configurations, which is the subject of the
next section.

Figure 4: Normalized Fourier transforms: (a), (b)
$\tilde{E}_{\text{I}}(\omega)$ [see Eq. (16)] for $E/\varepsilon=4$,
$T=10m^{-1}$, and two different values of $\tau$; (c), (d)
$\tilde{E}_{\text{II}}(\omega)$ [see Eq. (17)] for $T=40m^{-1}$ and two values
of $\tau$.

###  IV.3 Fourier analysis

Here we will examine the Fourier transforms of the external backgrounds (1)
and (2) defined via

$\tilde{E}(\omega)=\int\limits_{-\infty}^{\infty}E(t)\mathrm{e}^{i\omega
t}dt.$ (15)

We find

$\displaystyle\tilde{E}_{\text{I}}(\omega)$ $\displaystyle=$
$\displaystyle\pi\omega\bigg{[}\frac{ET^{2}}{\sinh\pi\omega
T/2}+\frac{\varepsilon\tau^{2}}{\sinh\pi\omega\tau/2}\bigg{]},$ (16)
$\displaystyle\tilde{E}_{\text{II}}(\omega)$ $\displaystyle=$
$\displaystyle\frac{2\pi^{2}E}{\omega(\pi^{2}-\omega^{2}\tau^{2})}\,\sin\frac{\omega(T+\tau)}{2}\cos\frac{\omega\tau}{2}.$
(17)

In Fig. 4 we present the functions (16) and (17) divided by
$\tilde{E}_{\text{I}}(0)=2(ET+\varepsilon\tau)$ and
$\tilde{E}_{\text{II}}(0)=E(T+\tau)$, respectively. All of the Fourier
transforms are centered at $\omega=0$ and quite narrow for relatively large
$\tau$. With decreasing $\tau$, the spectra gain high-frequency contributions,
which may stimulate the pair production process. For instance, in Fig. 4(b)
the presence of a short and weak Sauter pulse alters the large-$\omega$ part
of the spectrum providing it with broad wings. Although the relative
contribution of the high-frequency intervals seems tiny, these harmonics lead
to a significant enhancement of the pair number. In the scenario concerning a
rapid switch off of the external field, the qualitative structure of the
Fourier spectra is the same apart from the presence of an oscillatory pattern.
For smaller values of $\tau$, one observes an increase in the high-frequency
contributions. More pronounced oscillations in this case appear due to the
sharp edges of the rectangular-like pulse [for $\tau=0$ one obtains a well-
known Fourier transform $\sim\sin(\omega T/2)/\omega$, which is very similar
to the function displayed in Fig. 4(d)].

As was demonstrated above, in the case of weak and rapidly varying fields, one
can accurately describe the pair production process by means of PT although
the main contribution in the spectrum may correspond to $\omega=0$. This means
that depending on the field parameters, the large-$\omega$ harmonics may
completely govern particle production. However, within the enhancement
mechanisms investigated in our study the zero-energy contribution plays also
an essential role. The similarity between the Fourier transforms (16) and (17)
and between their dependence on the external field parameters indicates again
that the two effects of the enhancement are indeed related. In the next
section, we will discuss the experimental feasibility of the two mechanisms.

$E/E_{\text{c}}$ |  $I$ ($\text{W}/\text{cm}^{2}$) |  $\lambda_{\text{I}}$ ($\mu$m) |  $\lambda_{\text{II}}$ ($\mu$m)
---|---|---|---
0.1 | $2.3\times 10^{27}$ | $8\times 10^{-6}$ | $2\times 10^{-5}$
0.05 | $5.8\times 10^{26}$ | $2\times 10^{-5}$ | $7\times 10^{-5}$
0.02 | $9.2\times 10^{25}$ | $4\times 10^{-5}$ | $3\times 10^{-4}$
0.01 | $2.3\times 10^{25}$ | $8\times 10^{-5}$ | $8\times 10^{-4}$
0.005 | $5.8\times 10^{24}$ | $2\times 10^{-4}$ | $2\times 10^{-3}$
0.002 | $9.2\times 10^{23}$ | $4\times 10^{-4}$ | $9\times 10^{-3}$
0.001 | $2.3\times 10^{23}$ | $8\times 10^{-4}$ | $2\times 10^{-2}$
Table 1: Characteristic scale (wavelength) of the weak pulse
$\lambda_{\text{I}}$ and the switching profile $\lambda_{\text{II}}$ necessary
for the pair-production enhancement as a function of the electric field
amplitude $E$ ($I$ is the corresponding peak intensity).

##  V Experimental prospects

To find out whether the enhancement mechanisms examined in this study can be
implemented experimentally, we will now evaluate the spatiotemporal scale of
the rapid weak field and switch-off parts, respectively.

The onset of the pair-production enhancement in the context of the DASE
corresponds to $\gamma_{\text{c}}\equiv m/|eE\tau|=\gamma_{\text{crit}}\sim
1$. We note that depending on the external field profile, the critical value
of the Keldysh parameter can be different. For instance,
$\gamma_{\text{crit}}=\pi/2$ schuetzhold_prl_2008  or
$\gamma_{\text{crit}}=\sqrt{\ln(E/\varepsilon)}$ linder_prd_2015 .
Nevertheless, these values are basically of the order of unity, so we will
assume $\gamma_{\text{crit}}=1$. Let us now estimate the characteristic
wavelength of the rapid weak pulse as $\lambda_{\text{I}}=2\tau$. Then the
onset of dynamical assistance is approximately given by
$\lambda_{\text{I}}=2(E_{\text{c}}/E)m^{-1}$.

As was demonstrated above, by rapidly switching off the external field, one
can enhance the particle yield only if the LCFA is not justified. Although the
validity of the LCFA depends on the duration of the plateau region $T$, we
will employ the necessary condition $|eE|^{3/2}\tau\lesssim m^{2}$ and obtain
the following upper bound:
$\lambda_{\text{II}}=2(E_{\text{c}}/E)^{3/2}m^{-1}$.

It is known that the presence of the temporal oscillations and large volume
factor $V$ effectively lower the threshold for $e^{+}e^{-}$ pair production,
so that the field amplitude $E=0.1E_{\text{c}}$ seems already sufficient for
observing the Schwinger mechanism (see, e.g., Ref. aleksandrov_prd_2022 ). We
vary $E$ from $0.001E_{\text{c}}$ to $0.1E_{\text{c}}$ and display our
estimates in Table 1.

Since modern laser facilities generating most intense laser pulses operate in
the regime of $\lambda\sim 1\leavevmode\nobreak\ \mu\text{m}$, shaping the
laser field on the scale of $\lambda_{\text{II}}$ appears to be completely
unrealistic. On the other hand, it is possible to generate relatively weak
electromagnetic pulses having a much shorter wavelength. While such pulses
cannot significantly change the strong-field profile, they can be superimposed
on the strong component according to the idea of the DASE. This scenario
should be feasible as the values of $\lambda_{\text{I}}$ correspond to
realistic energy scales: the data presented in Table 1 yields the photon
energy of 1–150 keV. Finally, we point out that the parameters of the weak
high-frequency pulse assisting the Schwinger mechanism might not allow one to
treat it as a classical background. In this case, one should introduce quanta
of the electromagnetic field (photons) in order to describe the DASE (see,
e.g., Ref. dunne_prd_2009 ).

##  VI Conclusion

In the present study, we investigated the dynamically assisted Schwinger
effect and pair production in a strong external background being rapidly
switched on and off. It was shown that both scenarios may lead to a
substantial increase in the total number of pairs produced. It turns out that
the two mechanisms are quite similar to each other as they both require the
presence of a strong background and high-frequency harmonic. By varying the
field parameters, we analyzed the onset of the pair-production enhancement and
discussed the experimental prospects of the two scenarios.

Even in the case of a constant background, the adiabatic number of pairs
always incorporates dynamical contribution due to a sharp switch off of the
external field at the corresponding time instant. However, this initially
nonphysical contribution can be used to maximize the particle yield by
introducing a smooth switching function. The crucial property of this
enhancement mechanism is the characteristic time scale of the switching
profile. It was shown that this scale is at least two-three orders of
magnitude smaller than the period of electromagnetic radiation having a
wavelength of 1 $\mu\text{m}$. This suggests that shaping laser fields in
order to take advantage of a sharp switch off is unfeasible from the
experimental viewpoint.

On the other hand, one can attempt to design a setup where the pair production
process occurs in a combination of a strong electric background and a rapidly
oscillating field. According to our estimates, this scenario should involve
high-frequency radiation with photon energy of 1–150 keV depending on the
strong-field amplitude. As we were interested only in the onset of dynamical
assistance, photons of higher energy may, in fact, be even more advantageous,
provided the pair production process remains nonperturbative.

###### Acknowledgements.

The study was funded by RFBR and ROSATOM, project No. 20-21-00098. I.A.A.
acknowledges the support from the Foundation for the advancement of
theoretical physics and mathematics “BASIS”.

##  Appendix A Asymptotic expansion of the function (9)

Substituting $z=\tan x$ in Eq. (9), one obtains

$f(\xi)=\frac{4}{\pi}\int\limits_{0}^{\infty}\frac{\mathrm{e}^{-\pi
z^{2}/\xi}}{(1+z^{2})^{3}}\,dz.$ (18)

This integral can be evaluated analytically in terms of the Gauss error
function, but to find the small-$\xi$ expansion it suffices to integrate the
Taylor series of $1/(1+z^{2})^{3}$ multiplied by the exponential function
containing $\xi$. This immediately brings us to

$f(\xi)=\frac{2}{\pi}\sqrt{\xi}\Big{[}1-\frac{3}{2\pi}\,\xi+\frac{9}{2\pi^{2}}\,\xi^{2}+\mathcal{O}(\xi^{3})\Big{]}.$
(19)

Figure 5: Leading-order Feynman diagram describing electron-positron pair
production. The line with the cross corresponds to the interaction with the
external classical background. The particle number is proportional to the
amplitude squared and integrated over the electron (positron) momentum.

##  Appendix B Perturbation theory. Individual Sauter pulse

Let us evaluate the number of pairs to the leading (second) order in the
external field amplitude. The corresponding Feynman diagram is depicted in
Fig. 5, where the cross represents the interaction with the classical
background. A general expression in the case of a spatially homogeneous
external field has the following form:

$\nu^{\text{(PT)}}=\frac{e^{2}}{2(2\pi)^{3}}\int
d\boldsymbol{p}\,\frac{m^{2}+\boldsymbol{p}_{\perp}^{2}}{p_{0}^{4}}\,\Bigg{|}\int\limits_{-\infty}^{\infty}E(t)\mathrm{e}^{2ip_{0}t}\,dt\Bigg{|}^{2},$
(20)

where
$p_{0}=\sqrt{m^{2}+\boldsymbol{p}^{2}}=\sqrt{m^{2}+p_{\parallel}^{2}+\boldsymbol{p}_{\perp}^{2}}$.
To perform the calculations for a Sauter pulse
$E(t)=\varepsilon/\cosh^{2}(t/\tau)$, we use the Fourier transform

$\int\limits_{-\infty}^{\infty}\frac{\mathrm{e}^{i\omega
t}}{\cosh^{2}(t/\tau)}\,dt=\frac{\pi\omega\tau^{2}}{\sinh(\pi\omega\tau/2)},$
(21)

which can be obtained, e.g., by means of the residue theorem [it was also used
in Eq. (16)]. It brings us to Eq. (10). The PT approach accurately predicts
the particle yield if $\gamma_{\varepsilon}\equiv m/|e\varepsilon\tau|\gg 1$.
Note that the integrand represents the number density of the particles
produced and coincides with Eq. (44) from Ref. gavrilov_prd_1996  (one should
take into account the spin degeneracy).

In order to derive the small-$\tau$ expression (11), we will perform the
integration in terms of dimensionless $\xi=\tau p_{\parallel}$ and
$\eta=\tau^{2}|\boldsymbol{p}_{\perp}|^{2}$. One obtains

$\frac{\nu^{\text{(PT)}}}{m\tau}=\frac{\varepsilon^{2}}{2E_{\text{c}}^{2}}\,m^{3}\,\mathcal{I}(m^{2}\tau^{2}),$
(22)

where

$\mathcal{I}(\zeta)=\int\limits_{0}^{\infty}d\xi\int\limits_{0}^{\infty}d\eta\,\frac{\eta+\zeta}{\xi^{2}+\eta+\zeta}\,\frac{1}{\sinh^{2}(\pi\sqrt{\xi^{2}+\eta+\zeta})}.$
(23)

Substituting then $x=\xi^{2}+\eta+\zeta$, $y=\xi^{2}$ and integrating over $y$
from $0$ to $x-\zeta$, we arrive at

$\mathcal{I}(\zeta)=\frac{1}{3}\int\limits_{\zeta}^{\infty}\bigg{(}2+\frac{\zeta}{x}\bigg{)}\,\frac{\sqrt{x-\zeta}}{\sinh^{2}(\pi\sqrt{x})}\,dx.$
(24)

The contribution which comes from the second term in the parentheses vanishes
for $\zeta\to 0$ as it approximately equals $\sqrt{\zeta}/(6\pi)$. The value
$\mathcal{I}(0)$ is nonzero due to the first term:

$\mathcal{I}(0)=\frac{2}{3}\int\limits_{0}^{\infty}\frac{\sqrt{x}}{\sinh^{2}(\pi\sqrt{x})}\,dx=\frac{2}{9\pi}.$
(25)

Together with Eq. (22), it yields Eq. (11).

Finally, we point out that in the case of $1+1$ QED, there are no
perpendicular components $\boldsymbol{p}_{\perp}$, so we have a one-
dimensional integral in Eq. (10). Instead of Eq. (11), one obtains here

$\nu^{\text{(PT,
1+1)}}\approx\frac{1}{8\pi^{2}}\,\frac{\varepsilon^{2}}{E_{\text{c}}^{2}}\,(m\tau)^{2}m.$
(26)

##  Appendix C Perturbation theory for the $\cos^{2}$ pulse

The external background (2) for $T=0$ reads $E(t)=E\cos^{2}[\pi
t/(2\tau)]\theta(\tau-|t|)$. First, we find

$\int\limits_{-\tau}^{\tau}\cos^{2}\Big{(}\frac{\pi
t}{2\tau}\Big{)}\,\mathrm{e}^{i\omega
t}\,dt=\frac{\pi^{2}\sin\omega\tau}{\omega(\pi^{2}-\omega^{2}\tau^{2})}.$ (27)

This expression can be obtained from Eq. (17) by setting $T=0$. Then the PT
expression (20) can be evaluated by means of the same substitutions as those
employed in Appendix B. We obtain

$\frac{\nu^{\text{(PT)}}}{m\tau}=\frac{E^{2}}{2E_{\text{c}}^{2}}\,m^{3}\,\mathcal{J}(m^{2}\tau^{2}),$
(28)

where

$\mathcal{J}(\zeta)=\frac{\pi^{2}}{768}\int\limits_{\zeta}^{\infty}\frac{\sqrt{x-\zeta}(2x+\zeta)\sin^{2}2\sqrt{x}}{x^{3}(x-\pi^{2}/4)^{2}}\,dx.$
(29)

This function has a finite nonzero limit as $\zeta$ tends to zero:

$\mathcal{J}(0)=\frac{\pi^{2}}{384}\int\limits_{0}^{\infty}\frac{\sin^{2}2\sqrt{x}}{x^{3/2}(x-\pi^{2}/4)^{2}}\,dx.$
(30)

Substituting $y=2\sqrt{x}$, we recast this expression into

$\mathcal{J}(0)=\frac{\pi^{2}}{6}\int\limits_{0}^{\infty}\frac{\sin^{2}y}{y^{2}(y^{2}-\pi^{2})^{2}}\,dy=\frac{1}{8\pi}.$
(31)

The last integral can be computed by means of the residue theorem. Equations
(28) and (31) lead to Eq. (13).

In the case of $1+1$ QED, one obtains

$\nu^{\text{(PT,
1+1)}}\approx\frac{1}{32\pi^{2}}\,\frac{E^{2}}{E_{\text{c}}^{2}}\,(m\tau)^{2}m.$
(32)

## References

  * (1) F. Sauter, Z. Phys. 69, 742 (1931).
  * (2) W. Heisenberg and H. Euler, Z. Phys. 98, 714 (1936).
  * (3) J. Schwinger, Phys. Rev. 82, 664 (1951).
  * (4) A. M. Fedotov, E. G. Gelfer, K. Yu. Korolev, and S. A. Smolyansky, Phys. Rev. D 83, 025011 (2011).
  * (5) D. B. Blaschke, B. Kämpfer, S. M. Schmidt, A. D. Panferov, A. V. Prozorkevich, and S. A. Smolyansky, Phys. Rev. D 88, 045017 (2013).
  * (6) A. Blinne and H. Gies, Phys. Rev. D 89, 085001 (2014).
  * (7) J. Zahn, J. Phys. A 48, 475402 (2015).
  * (8) A. Otto, T. Nousch, D. Seipt, B. Kämpfer, D. Blaschke, A. D. Panferov, S. A. Smolyansky, and A. I. Titov, J. Plasma Phys. 82, 655820301 (2016).
  * (9) A. Ilderton, Phys. Rev. D 105, 016021 (2022).
  * (10) R. Schützhold, H. Gies, and G. Dunne, Phys. Rev. Lett. 101, 130404 (2008).
  * (11) M. Orthaber, F. Hebenstreit, and R. Alkofer, Phys. Lett. B 698, 80 (2011).
  * (12) C. Fey and R. Schützhold, Phys. Rev. D 85, 025004 (2012).
  * (13) C. Kohlfürst, M. Mitter, G. von Winckel, F. Hebenstreit, and R. Alkofer, Phys. Rev. D 88, 045028 (2013).
  * (14) I. Akal, S. Villalba-Chávez, and C. Müller, Phys. Rev. D 90, 113004 (2014).
  * (15) F. Hebenstreit and F. Fillion-Gourdeau, Phys. Lett. B 739, 189 (2014).
  * (16) M. F. Linder, C. Schneider, J. Sicking, N. Szpak, and R. Schützhold, Phys. Rev. D 92, 085009 (2015).
  * (17) A. Otto, D. Seipt, D. Blaschke, B. Kämpfer, and S. A. Smolyansky, Phys. Lett. B 740, 335 (2015).
  * (18) A. D. Panferov, S. A. Smolyansky, A. Otto, B. K. Kämpfer, D. B. Blaschke, and L. Juchnowski, Eur. Phys. J. D 70, 56 (2016).
  * (19) I. A. Aleksandrov, G. Plunien, and V. M. Shabaev, Phys. Rev. D 97, 116001 (2018).
  * (20) A. Otto, H. Oppitz, and B. Kämpfer, Eur. Phys. J. A 54, 23 (2018).
  * (21) I. Sitiwaldi and B.-S. Xie, Phys. Lett. B 777, 406 (2018).
  * (22) G. Torgrimsson, C. Schneider, and R. Schützhold, Phys. Rev. D 97, 096004 (2018).
  * (23) X.-G. Huang and H. Taya, Phys. Rev. D 100, 016013 (2019).
  * (24) S. Villalba-Chávez and C. Müller, Phys. Rev. D 100, 116018 (2019).
  * (25) O. Olugh, Z. L. Li, B.-S. Xie, Phys. Lett. B 802, 406 (2020).
  * (26) H. Taya, Phys. Rev. Research 2, 023257 (2020).
  * (27) L.-J. Li, M. Mohamedsedik, and B.-S. Xie, Phys. Rev. D 104, 036015 (2021).
  * (28) I. A. Aleksandrov, G. Plunien, and V. M. Shabaev, Phys. Rev. D 94, 065024 (2016).
  * (29) I. A. Aleksandrov, G. Plunien, and V. M. Shabaev, Phys. Rev. D 96, 076006 (2017).
  * (30) I. A. Aleksandrov and C. Kohlfürst, Phys. Rev. D 101, 096009 (2020).
  * (31) I. A. Aleksandrov, V. V. Dmitriev, D. G. Sevostyanov, and S. A. Smolyansky, Eur. Phys. J. Special Topics 229, 3469 (2020).
  * (32) D. G. Sevostyanov, I. A. Aleksandrov, G. Plunien, and V. M. Shabaev, Phys. Rev. D 104, 076014 (2021).
  * (33) F. V. Bunkin and I. I. Tugov, Dokl. Akad. Nauk SSSR 187, 541 (1969) [Sov. Phys. Dokl. 14, 678 (1970)].
  * (34) N. B. Narozhny, S. S. Bulanov, V. D. Mur, and V. S. Popov, Phys. Lett. A 330, 1 (2004).
  * (35) G. V. Dunne, Q. Wang, H. Gies, and C. Schubert, Phys. Rev. D 73, 065028 (2006).
  * (36) F. Hebenstreit, R. Alkofer, and H. Gies, Phys. Rev. D 78, 061701(R) (2008).
  * (37) S. S. Bulanov, V. D. Mur, N. B. Narozhny, J. Nees, and V. S. Popov, Phys. Rev. Lett. 104, 220404 (2010).
  * (38) S. P. Gavrilov and D. M. Gitman, Phys. Rev. D 95, 076013 (2017).
  * (39) I. A. Aleksandrov, G. Plunien, and V. M. Shabaev, Phys. Rev. D 99, 016020 (2019).
  * (40) A. I. Nikishov, Zh. Eksp. Teor. Fiz. 57, 1210 (1969) [Sov. Phys. JETP 30, 660 (1970)].
  * (41) C. Kohlfürst, N. Ahmadiniaz, J. Oertel, and R. Schützhold, arXiv:2107.08741.
  * (42) I. A. Aleksandrov, A. Di Piazza, G. Plunien, and V. M. Shabaev, Phys. Rev. D 105, 116005 (2022).
  * (43) G. Dunne, H. Gies, and R. Schützhold, Phys. Rev. D 80, 111301(R) (2009).
  * (44) S. P. Gavrilov and D. M. Gitman, Phys. Rev. D 53, 7162 (1996).


# Inverse design and experimental realization of plasma metamaterials

Jesse A. Rodríguez<EMAIL_ADDRESS>Department of Mechanical
Engineering, Stanford University, Stanford, CA 94305  Mark A. Cappelli
Department of Mechanical Engineering, Stanford University, Stanford, CA 94305

###### Abstract

We apply inverse design methods to produce two-dimensional triangular-lattice
plasma metamaterial (PMM) devices which are then constructed and demonstrated
experimentally. Finite difference frequency domain simulations are used along
with forward-mode automatic differentiation to optimize the plasma densities
of each of the plasma elements in the PMM to perform beam steering and
demultiplexing under transverse magnetic polarization. The optimal device
parameters are then used to assign plasma density values to elements that make
up an experimental version of the device. Device performance is evaluated
against both the simulated results and human-designed alternatives, showing
the benefits and disadvantages of in-silico inverse design and paving the way
for future fully in-situ optimization.

††preprint: APS/123-QED

##  I Introduction

The inverse design technique is aptly named; it consists of the process of
optimizing a device given an a-priori set of performance metrics. In
electromagnetically active systems, the technique becomes a type of
constrained optimization problem with Maxwell’s equations serving as the
constraints. The performance metrics are encoded in an objective function that
increases with favorable performance, e.g. $L_{2}$ inner products of the
electromagnetic (EM) field solutions for the current design iteration and the
desired EM fields in regions of interest. This objective is a function of
user-defined parameters that encode the domain/configuration of the device.
This technique has been very fruitfully applied in EM systems over the last
decade, leading to novel devices that perform functions so complex and at such
high efficiencies that they would be impossible to realize through expert-
level human design [1, 2, 3, 4, 5, 6, 7]. Some examples of such devices are
high performance metalenses [8, 9], metasurfaces [10, 11], photonic crystals
[12, 13, 14, 15, 16], and even optical computing devices that can perform
functions such as boolean logic operations [17] and matrix-vector
multiplication [18].

As mentioned above, in order to construct an objective that enables
optimization of the device domain, the designer must come up with a
parameterization scheme. Unfortunately, the configuration space of physically
realizable devices is much smaller than the configuration space of devices
that can be modeled, so the domain is parameterized in a manner that is in
agreement with experimental facilities. Typically, given the limits of working
with existing dielectric materials which have a fixed dielectric response, the
domain is binarized and the parameters simply correspond to whether or not
material with a given dielectric constant is present at a location in the
domain [19, 20]. Of course, this limits the space of possible devices, but
nevertheless this approach has yielded startlingly strong results. There is
also some exciting recent work that involves the use of additive manufacturing
to produce two-dimensional devices with a finite number of grayscale
permittivity values within a certain interval through fine layering [21],
greatly expanding the configuration space for inverse design devices.

In this and previous studies [17, 22], we choose to encode our device domain
as a plasma metamaterial (PMM), that is, as a domain filled with tunable
plasma elements that can be reconfigured within a continuous range of
permittivity values. A diagram that illustrates the inverse design algorithm
in our context is shown in Figure 1. In a PMM, each of the device parameters
corresponds to the plasma frequency (proportional to the square root of the
plasma density) of one of the elements in the device. The plasma density is
varied in practice by tuning the power delivery to the plasma elements; in our
case by increasing power supply voltage for a low-temperature gas discharge,
but possible by many methods in practice such as tuning the power delivered to
a resonator cavity or of the source of a laser-generated plasma. In-silico,
the domain is encoded as such and then a finite difference frequency domain
simulation (FDFD) is computed using an autograd-compliant Maxwell solver
called Ceviche [23]. ’Autograd-compliant’ here means that by defining an
objective that depends on the EM fields produced by the solver ($L_{2}$ inner
products as mentioned before), numerical gradients can easily be computed with
respect to the parameters that encode the domain using a widely-known python
package. Ceviche supplies the gradient of the fields with respect to the input
parameters via forward-mode differentiation, a technique where the Jacobian is
computed during the forward pass, e.g. as the solver progresses to produce the
field solutions. The resulting gradient is then used to iteratively optimize
the device.

Figure 1: (a) Flow chart describing the algorithmic design of our PMM array.
$J$ represents the modal source for the FDFD simulation. Examples of the
training parameter vector $\rho$, permittivity matrix $\varepsilon_{r}$ and
simulated E fields for a simple PMM device are provided in (b).

The choice to encode the device as a PMM has a number of advantages that are a
direct result of the unique dielectric response of the plasma. Figure 2 shows
the way the dielectric permittivity of a non-magnetized plasma varies with
incident wave frequency. The permittivity (neglecting the ion’s response) is
of the following form, known as the Drude permittivity:

$\displaystyle\varepsilon=1-\frac{\omega_{p}^{2}}{\omega^{2}+i\omega\gamma}$

where $\gamma$ is the collision/damping rate,
$\omega_{p}^{2}=\frac{n_{e}e^{2}}{\varepsilon_{0}m_{e}}$ is the plasma
frequency squared, $n_{e}$ is the electron density, $e$ is the electron
charge, $m_{e}$ is the electron mass, and $\varepsilon_{0}$ is the free-space
permittivity. We can see from this functional form that in the collisionless
limit, the plasma permittivity can be tuned continuously in the interval
$(-\infty,1)$ by varying the plasma density for a given operating frequency
$\omega$. By choosing an operating frequency that is close to the plasma
frequency of the elements, we can easily make drastic changes to their
electromagnetic response; i.e. that of a metallic element or a near-zero-index
medium. This dispersion relation not only makes the configuration space
infinite as opposed to a binarized parameterization scheme, but much more
importantly, this allows a single array of plasma elements to perform several
functions, and since the plasma elements can be quickly deactivated, you also
can switch the device on and off very quickly $\sim 10$ ms or less. More
detail about switching times is given in the methods section below.

Figure 2: Permittivity of non-magnetized plasma vs incident wave frequency
normalized by the plasma frequency.

Most importantly, since all the parameters can be varied precisely in the
physical device, it is possible to perform the inverse design algorithm
entirely in-situ. This expedites the process of designing these devices by
eliminating the need for expensive numerical simulations and also removes the
errors associated with inaccurate modeling of the device elements. Fully in-
situ inverse design is just now becoming experimentally feasible [24], and the
use of plasma in inverse design devices offers high-precision, continuous
reconfigurability of device elements, attributes which are ideally suited for
in-situ optimization. In this study, we will examine the merits and
disadvantages of creating inverse design devices the conventional way, by
optimizing them in-silico and then seeing how they perform experimentally. The
development of the experimental platform which is detailed in this study will
enable the transition to in-situ inverse design in short order, but first we
must examine how our best in-silico optimization efforts translate to
experimental results. In the following sections we will numerically optimize
and experimentally demonstrate a two-dimensional PMM that can perform both
beam steering and demultiplexing.

##  II Methods

###  II.1 Computational Inverse Design

Figure 3(a) provides a diagram of our PMM device. The device has a triangular
lattice with lattice spacing $a=24$mm and is in the shape of a hexagon with
six elements on each side, yielding 91 elements total. The outer diameter of
the $\varepsilon=3.8$ quartz tubes is 15mm and the inner diameter is 13mm. In
figure 3(b), the simulation domain is shown with an arbitrary average plasma
density distribution among the elements, where each domain parameter
corresponds to the space-averaged plasma frequency of one of the elements. The
elements are modeled as either spatially uniform or as having a density
profile that is axisymmetric and a zeroth-order Bessel function of the first
kind of radius, with the first root at $r=6.5$mm. In the collisional model
(which is also non-uniform), the simulated elements have a collision frequency
of $\gamma=1$ GHz. The blue line indicates where the source is located in the
domain, and the red lines indicate the locations in the domain where the
objective function is evaluated; we refer to these as ’probe slices’. The
probe slices surrounding the crystal are present to discourage leaking out of
the device, while the probe slices within the exit horns either promote
propagation into the correct horn or discourage propagation into the incorrect
horn. The relative permittivity of the input and output horn walls was set to
$\varepsilon=-1000$ to serve as a lossless metal, and the horn geometry was
chosen to match our experimental microwave horns (A INFO LB-20180 2-18 GHz).

Figure 3: (a) Schematic of PMM and (b) relative permittivity of the simulation
domain with source and objective evaluation locations denoted in blue and red
respectively. The domain is shown with an arbitrary device configuration as an
example.

The basic idea is that a modal source is introduced at the input horn and
allowed to scatter through the training region before collecting at the
desired output horn. In the computational optimization phase, the fields are
propagated through the domain via Ceviche. As mentioned earlier, Ceviche
solves Maxwell’s equations using a FDFD solver, and therefore all computed
devices described here represent the steady state solution achieved after some
characteristic time. The simulation domain was discretized using a resolution
of 50 pixels per lattice constant $a$ and in each case presented dimensions of
$N_{x}=1000$ and $N_{y}=1200$ pixels. The final device configurations were
tested at the higher resolution of 150 pixels$/a$ and yielded results that
were qualitatively identical. A perfectly-matched boundary layer (PML) $2a$ in
width was applied along the domain boundaries. The simulated $\mathbf{E}$
fields are masked to compute integrals along the probe slices within the
problem geometry. These integrals are used to calculate the objective
$\mathcal{L}(\rho)$. Our simulation tool, Ceviche, is then used to compute
numerical gradients of the objective with respect to the training parameters
via forward-mode differentiation [23]. We then use the Adam optimization
algorithm [25] (gradient ascent with momentum, in essence) to iteratively
adjust $\rho$ and thereby maximize $\mathcal{L}$. Optimization was conducted
with learning rates ranging from $0.01-0.05$, and the default Adam
hyperparameters $\beta_{1}=0.9$, and $\beta_{2}=0.999$.

The polarization of the input source has a strong effect in devices of this
nature, either $E_{z}$ ($E$ out of the page), which we call the transverse
magnetic (TM) polarization, or $E_{x}$ ($H$ out of the page), which we call
the transverse electric (TE) polarization, with the latter case benefiting
(and suffering, depending on the desired device functionality) from the
presence of localized surface plasmon (LSP) modes [26, 27, 28, 29], while the
former makes more direct use of dispersive and refractive effects [27, 30,
31]. In our previous studies [17, 22], we see that for devices which seek to
preserve an input signal like those we present here, the TE polarization leads
to poor performance, so we choose to focus on the TM mode alone in this study.
Thus, we will not expect to see any small, sub-wavelength-scale field
structures around our plasma elements or the large transmission losses which
are indicative of the presence of LSP modes.

To summarize, the PMM optimization problem can be expressed as

$\displaystyle\max_{\rho}\quad\mathcal{L}(\mathbf{E})$
$\displaystyle\text{given}\quad\nabla\times\frac{1}{\mu_{0}}\nabla\times\mathbf{E}-\omega^{2}\varepsilon(\omega,\rho)\mathbf{E}=-i\omega
J.$

where $\mathcal{L}$ is the objective composed of a set of $L_{2}$ inner
product integrals of the simulated field with the desired propagation mode, E
is the electric field, $\rho$ is an $n$-dimensional vector (where $n$ is the
number of elements) that contains the parameters that set the plasma frequency
of each of the PMM elements, $\mu_{0}$ is the vacuum permeability, $\omega$ is
the field frequency, $\varepsilon(\omega,\rho)$ is the spatially-dependent
permittivity that is encoded by $\rho$, and $J$ is a current density used to
define a fundamental modal source at the input horn. In practice, the
permittivity distribution among the plasma elements is controlled by varying
the discharge current (which therein alters the plasma density) in each of the
PMM elements according to the Drude model (discussed above in the
introduction), which is where the dependence on $\omega$ arises.

The beam steering objective is the $L_{2}$ inner product of the simulated
field with an $m=1$ propagation mode (in order to preserve the input mode)
evaluated at the probe in the desired exit minus the sum of the integrated
field intensity at the incorrect exit and the other probe slices to discourage
loss:

$\mathcal{L}_{steer}=\int E\cdot E_{m=1}^{*}dl_{\text{desired
exit}}-\sum_{\text{prb}\neq\text{desired exit}}\int|E|^{2}dl_{\text{prb}}.$

The optimization objective for the demultiplexer is quite different from the
beam steering case since the objective must take into account two separate
simulations as the permittivity of the elements is frequency-dependent. With
this in mind, the demultiplexer objective is:

$\displaystyle\mathcal{L}_{mp}=$ $\displaystyle\left(\int E_{\omega_{1}}\cdot
E_{m=1}^{*}dl_{\omega_{1}\text{ exit}}\right)\left(\int E_{\omega_{2}}\cdot
E_{m=1}^{*}dl_{\omega_{2}\text{ exit}}\right)$
$\displaystyle-\int|E_{\omega_{1}}|^{2}dl_{\omega_{2}\text{
exit}}-\int|E_{\omega_{2}}|^{2}dl_{\omega_{1}\text{ exit}}$
$\displaystyle-\sum_{\text{prb}\neq\text{wvg
exits}}\left[\int|E_{\omega_{1}}|^{2}dl_{\text{prb}}+\int|E_{\omega_{2}}|^{2}dl_{\text{prb}}\right],$

where $E_{\omega_{1}}$ is the simulated field for the first frequency and
$E_{\omega_{2}}$ is the simulated field for the second frequency. The first
term in the objective is a product so that it ensures that the fields reach
the exit horn for both operating frequencies (therein avoiding local minima
where one gets good performance for just one of the two frequencies), the
second two terms discourage leakage into the incorrect exits for each
frequency, and the final term discourages leakage out of the device for both
operating frequencies. For both objectives, each term is normalized by its
value in the initial design iteration.

The maximum plasma frequency of the elements is set to be 8 GHz in the
computational optimization procedure as that was our conservative estimate for
the nominal maximum plasma density of our experimental sources. The optimal
parameters from the computational inverse design process are then used to set
the configuration of the experimental device according to the mapping defined
below in Section II.2

###  II.2 Experimental Realization

The experimental plasma elements in this study are custom-manufactured UV
germicidal discharge tubes made of high-purity quartz ($\varepsilon=3.8$) with
an inner tube diameter of 13mm and an outer diameter of 15mm. The tubes are
filled with argon gas at a fill pressure between 200 and 300 Pa, with a small
amount of Mercury present in each discharge to provide the UV radiation. The
use of UV discharge tubes is necessary because other types of off-the-shelf
discharge tubes (for lighting applications) often have a fluorescent coating
applied to the glass that is opaque to microwave-range radiation. When
operating, the tubes are estimated to have a gas temperature between 315 and
330 K according to infrared thermometer measurements. The discharges are
designed to operate with alternating current (AC) at a frequency between 20
and 50 kHz depending on the ballast circuit that is used. The nominal
operating current and voltage is 24V and 0.6A, for a power of $\approx$ 15W.
To push the discharges to higher plasma densities, we drive them with 24W DC
ballasts (Beasun RL15-425-18D24), each having its own dedicated programmable
DC power supply (Longwei LW-3010). The DC power supplies are connected to a
RS485 serial communication bus (with maximum 32 power supplies per RS485 bus)
and accept MODBUS protocol serial commands from a custom python library via a
USB to RS485 converter. A schematic and photograph of the experimental scheme
is included below in Figure 4. In practice, the discharges can be tuned
through a fairly large range of plasma frequencies by setting the current of
the DC power supply to its maximum value and then limiting the output voltage.
The discharges typically operate from a maximum DC power supply voltage of 32
V down to about 6 V and then fail to maintain the discharge at a plasma
frequency of about 2.4-4 GHz. We refer to this as the ’voltage tuning’ mode.
If the power supply is set to maximum voltage and the current is instead
limited, a different operating mode can be reached, allowing us to access
lower plasma frequency values all the way down to 0.4 GHz before the discharge
extinguishes. For the high power ballasts used in this study, we achieve
steady operation in the current tuning mode from approximately 60 mA to 200
mA.

Figure 4: (a) Schematic of experimental PMM system and (b) photographs of the
different components of the apparatus including the DC power supplies (left),
the ballasts (middle), and the PMM itself activated in an arbitrary
configuration (right).

In order to use the parameters from the in-silico optimization which
correspond to the plasma frequency of each of the plasma elements, we need to
generate a mapping between the plasma frequency and the power supply settings
that produce that frequency. This is quite challenging as we do not know the
precise manufacturer fill pressure or gas temperature and the plasma density
in our operating range is too small to use optical diagnostics for direct
measurements. To combat this and obtain an approximate mapping from the power
supply settings to plasma frequency, we use a 0-D kinetic simulation of the
electron energy distribution function (EEDF) within a plasma discharge,
BOLSIG+ [32, 33], assuming a known reduced electric field, $E/n$, where $n$ is
the estimated neutral gas number density. Our estimate for the neutral gas
number density is obtained using the ideal gas law along with estimates of the
fill pressure and neutral gas temperature. We measure the time-varying voltage
across and current through the plasma discharge for a range of operating
conditions. The RMS voltage, estimated neutral gas number density, and
estimates of the length of the discharge enables estimates of $E/n$. The EEDF
determined using BOLSIG+ yields the electron mobility. The product of the
mobility and electric field provides the space-averaged electron drift
velocity $u_{d}=\mu_{e}E$ within the discharge. Finally, we can estimate the
electron density using the measured current; $n_{e}=I_{RMS}/eu_{d}A_{tube,i}$,
yielding the plasma frequency. Here, $A_{tube,i}$ is the internal cross
sectional area of the discharge tube. We note that this is complex path to
determine the plasma frequency, where several sources of uncertainty may be
encountered.

Although the discharge frequency is much higher than the plasma recombination
rate and as such the plasma density is quasi-steady, we have confirmed via
high-speed photography that at the peak of the voltage cycle the plasma has
the structure of a DC glow discharge with cathode and Faraday dark spaces
clearly visible. At low-plasma density operating conditions, we also observe
striations in the positive column of the discharge. This means that when
considering the effective electric field applied to the plasma that we feed to
the BOLSIG+ kinetic model, we have to take into account the cathode fall
voltage and be sure that we use the correct discharge length. For our
discharges, the discharge length varies between 190mm in the current-limited
mode and 230mm in the voltage limited mode, and we estimate the cathode fall
voltage to be somewhere between 6 and 11 V [34]. Taking all of these
considerations into account, we perform a sweep of all parameters (fill
pressure, temperature, cathode fall voltage) and perform a one-term polynomial
fit that allows us to sweep between the highest and lowest-density cases by
varying a parameter $k\in[0,1]$. Since the plasma frequency values appeared to
be slightly too large for these types of discharges (owing to the many
uncertainties in determining plasma density), we also added a scale factor $S$
that would allow us to adjust the fit to lower values over the entire tuning
range. The quasi-experimental fit is included below in Figure 5.

Figure 5: The quasi-experimental mapping for plasma tuning (bottom) and the
plots including all the cases from the parameter sweep in BOLSIG+ that form
the basis for the mapping (top).

In addition to the high-speed photography images mentioned above, we have also
conducted Schottky diode measurements of the transient dielectric response of
the discharges as they are activated and deactivated, showing that the rise
time of the dielectric response is $\sim 10$ ms and the fall time is $\sim 60$
ms. Due to the current configuration of the device where each RS485 bus hosts
up to 25 power supplies, the actual time to set and activate the array is much
longer, about 7 seconds, since each discharge is set and activated in series
to avoid cross-talk within each RS485 bus. In practice, should a user need
faster actuation, they need only parallelize the computational array setting
procedure and use more RS485 buses. The theoretical limit is to have one cpu
process communicating directly with a single power supply at a time with as
many processes running in parallel as necessary. Since the power supplies
operate at 9600 baud, each message is 8 bytes, and each actuation requires two
messages (set voltage/current and activate), the total time added by the
communication protocol is about 15 ms, so the time to switch the entire array
on could be lowered to $\sim 25$ ms and the deactivation time could be reduced
to $\sim 75$ ms. Temporal modulation of the plasma frequency of an already-
running discharge could occur at a maximum of about 75 Hz. Of course, more
sophisticated communication protocols and interconnects along with specially
designed power supplies and ballast circuits could push these actuation times
down further. Experimental data supporting these claims about the discharge
structure and switching time can be found in the appendix of reference [35].

To perform the measurements described in the Results section, we first
construct the apparatus as pictured in Fig. 4, place the array into a custom-
made anechoic chamber consisting of a box with microwave absorber panels
attached to the interior surfaces, and then carry out an initial warm up of
the device by turning every element in the device on for 15 seconds of every
minute over a period of 10 minutes. This warm-up procedure makes sure that the
temperature of the discharge tube will be steady throughout the measurements
and it also allows the 3D-printed scaffolding to reach its steady operating
temperature such that it won’t flex/deform between measurements. Once the
warm-up period is concluded, the PMM is activated according to the optimal
parameters from the in-silico optimization procedure for 15 seconds of each
minute for as many minutes as it takes for all measurements to be collected.
As stated before, the PMM device can activate and set the operating conditions
for 91 elements in about 7 seconds. Python scripts are constructed to ensure
that the PMM is not allowed to cool to room temperature by failing to maintain
the 15 sec/min duty cycle.

For each measurement, $S_{21}$ and $S_{31}$ are measured by a Rohde and
Schwarz ZNB40 4-port vector network analyzer that was calibrated using a Rohde
and Schwarz ZN-Z54 calibration unit. The $S$-parameters are collected in
10,000 points from 2-12 GHz with an averaging factor of 10. For each device
objective, we test the optimal parameters for three different plasma element
models (uniform density profile, Bessel function profile, and Bessel function
profile with collisions, see ref. [22] for more details regarding the plasma
model) and sweep the fitting parameters $k$ and $S$ to attempt to find the
correct plasma frequency mapping. This corresponds to about 30 measurements
for each objective, and the best / most informative case is displayed below in
Section III.

##  III Results

###  III.1 Beam Steering

The simulated and experimental results for the beam steering objective where a
source delivers power to one of two receiving ports are presented below in
Figure 6. During the in-silico optimization, the plasma density within each
discharge is modeled as being either uniform with no collisions, non-uniform
with no collisions, or non-uniform with collisions with $k$ and $S$ used to
parameterize the relationship between the plasma density and operating voltage
or current (see Figure 5).

Figure 6: Field simulations $|\textbf{E}_{z}|$ and optimal device domains
$\Re[\varepsilon]$ (left) along with experimental transmission spectra (right)
for (a) the PMM device with all the elements inactive, (b) the 5 GHz operating
frequency optimal beam steering design and, (c) the 6 GHz operating frequency
optimal beam steering design. The vertical dashed black lines indicate the
design frequency in (b) and (c), and the fit parameters $k$ and $S$ are given
along with the plasma element model used (see Methods for details).

As expected, when all elements are inactive (Fig.6(a)) the spectra of both
ports are nearly identical within about 5 dB of one another at all
frequencies. The small differences between the two measured spectra (right
frame in Fig.6(a)) are due to minor defects in the device scaffold and in the
alignment of the microwave horns. In both the 5 and 6 GHz active
configurations, the transmission in port 2 (the desired port) is nearly
maximal at the frequency for which the device was optimized and we have about
10-15 dB isolation between the ports. In the 5 GHz configuration in
particular, there is a drop in transmission in port 3 at almost precisely the
operating frequency predicted by the simulations assuming a uniform plasma
model with a discharge voltage mapping to plasma density defined by $k=0.1$
and $S=1.0$. In the 6 GHz plasma active case we see even better isolation at
another frequency point, with over 30 dB isolation at about 4.7 GHz with a
transmission in port 2 almost as high as at the operating frequency. We
discuss the possible cause(s) for this and its significance in Section IV. In
both cases, the power supply setting to plasma density fitting parameters that
led to the best performance ($k\approx 0.2$ and $S=1.0$), suggest that the
actual maximum plasma density of the discharges when driven by their power
supplies is is about 14 GHz (see Figure 5 for details) although prolonged
operation at this power level causes damage to the discharge.

###  III.2 Demultiplexer

The simulated and experimental results for the demultiplexer objective are
presented below in Figure 7. Here, we use the same device as that above, but
seek to separate different source frequencies. For example, in panel (b) of
Figure 7, we direct 5.5 GHz and 6 GHz content sourced at port 1, to ports 2
and 3 respectively.

Figure 7: Field simulations $|\textbf{E}_{z}|$ and optimal device domains
$\Re[\varepsilon]$ (left) along with experimental transmission spectra (right)
for (a) the PMM device with all the elements inactive at 5 and 7 GHz, (b) the
5.5 / 6 GHz optimal demultiplexer design and, (c) the 5 / 7 GHz optimal
demultiplexer design. The vertical dashed black lines indicate the design
frequencies in (b) and (c), and the fit parameters $k$ and $S$ are given along
with the plasma element model used (see Methods for details).

Figure 7(a) shows the baseline inactive plasma case, with experiments once
again indicating that ports 2 and 3 receive similar signals (within 3-5 dB)
over a wide range of frequencies. With active plasma cases, we find the best
performance for the non-uniform model using $k=0.0$ and $S=0.7$ (similar to
the beam steering objective) for the 5.5 / 6 GHz demultiplexer (Figure 7(b)),
but for the 5 / 7 GHz demultiplexer (Figure 7(c)) we see the transmission
characteristics we want for the collisional plasma element model at $k=1.0$,
$S=1.0$. In the 5.5 / 6 GHz case, the transmission spectrum appears as one
would expect in the region of the operating frequencies, with the transmission
of the two ports crossing at the midpoint between the operating frequencies,
where we see about 10 dB isolation in the correct direction. Interestingly,
the same configuration performs the demultiplexer objective at the frequencies
of $\sim$5.25 GHz and 8.05 GHz with up to 40 dB isolation, hinting at the
potential of the device if the optimization was to be carried out _in-situ_.
The 5 / 7 GHz demultiplexer does not perform its objective properly at the
design frequencies, but instead, performs very well at the slightly shifted
frequencies of about 5.6 and 7.4 GHz.

Figure 8: Field simulations $|\textbf{E}_{z}|$ and optimal device domains
$\Re[\varepsilon]$ (left) along with experimental transmission spectra (right)
for the PMM device with our best attempt at an a priori human-designed 5 / 7
GHz demultiplexer. The vertical dashed black lines indicate the design
frequencies, and the fit parameters $k$ and $S$ are chosen according to the
best values for the devices conceived via inverse design.

In figure 8, we see the simulated and experimental results when a 5 / 7 GHz
demultiplexer is designed by hand by creating what would function as low-index
bridges to the correct output port for either design frequency. The general
principle is that port 3 is protected by a region of sources that have an
average plasma frequency $\sim 6$ GHz, so those elements serve as a reflective
barrier to the 5 GHz source but are transparent to the 7 GHz source. The edges
of the domain are set to their maximum plasma density to discourage leakage.
$k=0.0$, and $S=0.7$ were chosen as the fitting parameters to limit the
voltage of the edge elements an avoid damaging the discharges. In both the
simulation and the experiment, the device doesn’t appear to steer the 7 GHz
source to either output, and there is only slight isolation between the ports
in the correct direction at 5 GHz. The experimental spectrum is fairly
unremarkable without obviously meaningful structures at any frequency.

##  IV Discussion & Conclusion

The results presented in the previous section show that although we can
reliably get somewhat strong performance by optimizing our devices in-silico,
there are ultimately some shortcomings. In most cases, the configurations
actually resulted in better performance at frequencies other than the design
frequency as seen with the 6 GHz beam steering device in Fig. 6 and the 5.5 /
6 GHz demultiplexer in Fig. 7. These cases show that up to $\sim$40 dB
isolation is achievable in practice at various frequencies throughout the
frequency domain, indicating that the device has a high degree of dynamic
range (possibly due to some manner of Fabry-Pérot resonance effect), but the
in-silico optimizer was only able to produce 10-15 dB isolation at the design
frequencies.

While there can be a few reasons for these shortcomings, we attribute them to
simplifications made in the simulations. The in-silico optimization procedure
uses a 2D domain with plasma elements of varying levels of complexity. The
most sophisticated plasma element model that includes both a non-uniform
plasma density profile and collisional damping failed to produce the best
performance in all but one case, illustrating that we still lack the detailed
knowledge of the discharge properties that are needed to produce robust
results at design conditions. Furthermore, the system is not entirely
electromagnetically isolated as assumed in the simulations. To avoid
interferences due to scattering from the surroundings of the PMM, experiments
are carried out in a custom-fabricated anechoic chamber, but even this isn’t
entirely sufficient. In addition, the finite size of the tubes brings the
electrodes and the wiring into play, both of which undoubtedly interact with
the source fields. Unfortunately, a 3D simulation that takes many of these
non-ideal conditions into account is far too costly to use within an inverse
design scheme.

We also found that the mapping between the power supply settings (voltage or
current) and the plasma frequency of the discharge tubes, facilitated by the
electron-energy distribution solver, BOLSIG+, tends to overestimate the plasma
density, where only the lowest density case according to the BOLSIG+
calculations is compatible with the experimental results ($k=0$). One
interesting departure from the other results regarding the experimental
mapping is the 5 / 7 GHz demultiplexer that performs best for the collisional
plasma element model at $k=1.0$, $S=1.0$, indicating that when using that
element model, the optimal plasma frequencies from the in-silico optimization
are necessarily larger. Because the elements are lossy, the optimizer may be
biased to higher plasma densities to encourage scattering off of the plasma
elements instead of propagation through them. It’s also possible that the
collisionality is stronger in the experimental device at higher power (the
collision frequency is assumed to be 1 GHz in the in-silico collisional plasma
model), therefore, to match the simulated result one has to decrease the
plasma frequency in the experiment. Also, since the device does not function
at precisely the intended frequencies, it is possible that this configuration
was simply stumbled upon by chance. Much of these challenges can be avoided by
in-situ optimization.

The attempt at designing the device intuitively for objectives such as the
demultiplexer is an illustration of how we cannot hope to outperform inverse
design, particularly so when our model of the physical device is not precise.
Despite the shortcomings discussed above, the devices still show strong
evidence that the in-silico optimization can result in performance that at
least addresses the design criteria, if not at a high degree of efficiency.
Designing by intuition is therefore not an option; but of course we have also
shown here that the in-silico inverse design procedure fails to fully utilize
the dynamic range of the system elements. Luckily, since our device is
composed entirely of plasma elements, we can perform the inverse design
process entirely in-situ.

Fully in-situ inverse design addresses almost all the shortcomings described
above and also opens up many exciting opportunities. We no longer have to be
concerned with the fidelity of the modeling since we will be using the
physical device to perform the optimization. Since we find that the device is
capable of $\sim$50 dB isolation at various frequencies, we would expect to be
able to achieve this and likely more since those cases were obtained
accidentally. This also presents an opportunity to refine our element model.
Once a high-performance design is achieved via in-situ inverse design, the
parameters can be fed into a simulation domain where we modify the plasma
model until we see the same performance, allowing us to more thoroughly
understand the physics of our plasma sources. Moving to fully in-situ
optimization also removes some of the limitations on our geometry. As
mentioned before, 3D geometry would be far too computationally costly to use
for in-silico inverse design, but it incurs no extra cost in in-situ inverse
design. Thus, we could use device structures like that of ref. [27]. Even for
this 2D configuration, each iteration in the optimization procedure can take
10 minutes with a 64-core workstation, while the in-situ iterations an be
limited to mere seconds in practice regardless of the geometry. We have also
shown in prior work that our plasma sources have a pronounced gyrotropic
response when magnetized [36, 37], but materials with anisotropic permittivity
tensors like magnetized plasma can’t be modeled with our simulation tool. By
using a large Helmholtz configuration, entire PMM devices can be magnetized to
take advantage of the very rich physics inherent to magnetized plasmas.

In conclusion, we show in this study that PMM devices can be optimized in-
silico to perform beam steering and demultiplexing at a level of performance
that is not possible to obtain using conventional intuitive design methods.
Ultimately, the in-silico inverse design process has shortcomings that lead to
a failure to utilize the full potential of the PMM configuration. Since the
PMM device can be reconfigured in seconds or less, the inverse design process
can be performed entirely in-situ, side-stepping the computational cost and
modeling inaccuracies associated with in-silico inverse design, thereby
promising better performance.

###### Acknowledgements.

The authors would like to thank Dr. Benjamin Wang for his assistance in
procuring the experimental components for this study. This research is
supported by the Air Force Office of Scientific Research through a Multi-
University Research Initiative (MURI), Grant No. FA9550-21-1-0244, with Dr.
Mitat Birkan and Dr. Arje Nachman as the Program Managers. J.A.R. acknowledges
support from the Charles H. Kruger Stanford Graduate Fellowship.

## References

  * Su  _et al._ [2020] L. Su, D. Vercruysse, J. Skarda, N. V. Sapra, J. A. Petykiewicz, and J. Vučković, Nanophotonic inverse design with SPINS: Software architecture and practical considerations, Applied Physics Reviews 7, 11407 (2020).
  * Miller [2013] O. D. Miller, Photonic design: From fundamental solar cell physics to computational inverse design, arXiv preprint arXiv:1308.0212 (2013).
  * Liu  _et al._ [2013] D. Liu, L. H. Gabrielli, M. Lipson, and S. G. Johnson, Transformation inverse design, Optics express 21, 14223 (2013).
  * Hughes  _et al._ [2018] T. W. Hughes, M. Minkov, I. A. Williamson, and S. Fan, Adjoint method and inverse design for nonlinear nanophotonic devices, ACS Photonics 5, 4781 (2018).
  * Molesky  _et al._ [2018] S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, Inverse design in nanophotonics, Nature Photonics 12, 659 (2018).
  * Christiansen and Sigmund [2021] R. E. Christiansen and O. Sigmund, Inverse design in photonics by topology optimization: tutorial, JOSA B 38, 496 (2021).
  * Andrade  _et al._ [2019] N. M. Andrade, S. Hooten, S. A. Fortuna, K. Han, E. Yablonovitch, and M. C. Wu, Inverse design optimization for efficient coupling of an electrically injected optical antenna-led to a single-mode waveguide, Opt. Express 27, 19802 (2019).
  * Christiansen  _et al._ [2020] R. E. Christiansen, Z. Lin, C. Roques-Carmes, Y. Salamin, S. E. Kooi, J. D. Joannopoulos, M. Soljačić, and S. G. Johnson, Fullwave maxwell inverse design of axisymmetric, tunable, and multi-scale multi-wavelength metalenses, Opt. Express 28, 33854 (2020).
  * Meem  _et al._ [2020] M. Meem, S. Banerji, C. Pies, T. Oberbiermann, A. Majumder, B. Sensale-Rodriguez, and R. Menon, Large-area, high-numerical-aperture multi-level diffractive lens via inverse design, Optica 7, 252 (2020).
  * Pestourie  _et al._ [2018] R. Pestourie, C. Pérez-Arancibia, Z. Lin, W. Shin, F. Capasso, and S. G. Johnson, Inverse design of large-area metasurfaces, Opt. Express 26, 33732 (2018).
  * Chung and Miller [2020] H. Chung and O. D. Miller, Tunable metasurface inverse design for 80% switching efficiencies and 144° angular deflection, ACS Photonics 7, 2236 (2020), https://doi.org/10.1021/acsphotonics.0c00787 .
  * Borel  _et al._ [2004] P. I. Borel, A. Harpøth, L. H. Frandsen, M. Kristensen, P. Shi, J. S. Jensen, and O. Sigmund, Topology optimization and fabrication of photonic crystal structures, Opt. Express 12, 1996 (2004).
  * Minkov  _et al._ [2020] M. Minkov, I. A. Williamson, L. C. Andreani, D. Gerace, B. Lou, A. Y. Song, T. W. Hughes, and S. Fan, Inverse design of photonic crystals through automatic differentiation, ACS Photonics 7, 1729 (2020).
  * Burger  _et al._ [2004] M. Burger, S. J. Osher, and E. Yablonovitch, Inverse problem techniques for the design of photonic crystals, IEICE transactions on electronics 87, 258 (2004).
  * Lin  _et al._ [2016] Z. Lin, A. Pick, M. Lončar, and A. W. Rodriguez, Enhanced spontaneous emission at third-order dirac exceptional points in inverse-designed photonic crystals, Phys. Rev. Lett. 117, 107402 (2016).
  * Jiang  _et al._ [2020] J. Jiang, M.-K. Chen, and J. A. Fan, Deep neural networks for the evaluation and design of photonic devices, ArXiv abs/2007.0 (2020).
  * Rodríguez  _et al._ [2021] J. A. Rodríguez, A. I. Abdalla, B. Wang, B. Lou, S. Fan, and M. A. Cappelli, Inverse design of plasma metamaterial devices for optical computing, Phys. Rev. Applied 16, 014023 (2021).
  * Nikkhah  _et al._ [2023] V. Nikkhah, A. Pirmoradi, F. Ashtiani, B. Edwards, F. Aflatouni, and N. Engheta, Inverse-designed low-index-contrast structures on silicon photonics platform for vector-matrix multiplication (2023), arXiv:2302.00793 [physics.optics] .
  * Piggott  _et al._ [2015] A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vučković, Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer, Nature Photonics 9, 374 (2015).
  * Su  _et al._ [2018] L. Su, A. Y. Piggott, N. V. Sapra, J. Petykiewicz, and J. Vuckovic, Inverse design and demonstration of a compact on-chip narrowband three-channel wavelength demultiplexer, Acs Photonics 5, 301 (2018).
  * Huang  _et al._ [2023] Q. Huang, L. T. Gan, and J. A. Fan, Conformal volumetric grayscale metamaterials, Advanced Materials 35, 2204688 (2023), https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.202204688 .
  * Rodríguez and Cappelli [2022] J. A. Rodríguez and M. A. Cappelli, Inverse design of plasma metamaterial devices with realistic elements, Journal of Physics D: Applied Physics 55, 465203 (2022).
  * Hughes  _et al._ [2019] T. W. Hughes, I. A. D. Williamson, M. Minkov, and S. Fan, Forward-Mode Differentiation of Maxwell’s Equations, ACS Photonics 6, 3010 (2019).
  * Pai  _et al._ [2023] S. Pai, Z. Sun, T. W. Hughes, T. Park, B. Bartlett, I. A. D. Williamson, M. Minkov, M. Milanizadeh, N. Abebe, F. Morichetti, A. Melloni, S. Fan, O. Solgaard, and D. A. B. Miller, Experimentally realized in situ backpropagation for deep learning in photonic neural networks, Science 380, 398 (2023), https://www.science.org/doi/pdf/10.1126/science.ade8450 .
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2014), arXiv:1412.6980 [cs.LG] .
  * Righetti  _et al._ [2018] F. Righetti, B. Wang, and M. Cappelli, Enhanced attenuation due to lattice resonances in a two-dimensional plasma photonic crystal, Physics of Plasmas 25, 124502 (2018).
  * Wang  _et al._ [2019] B. Wang, J. A. Rodríguez, and M. A. Cappelli, 3D woodpile structure tunable plasma photonic crystal, Plasma Sources Science and Technology 28, 02LT01 (2019).
  * Sakai  _et al._ [2018] O. Sakai, A. Iwai, Y. Omura, S. Iio, and T. Naito, Wave propagation in and around negative-dielectric-constant discharge plasma, Physics of Plasmas 25, 031901 (2018), https://doi.org/10.1063/1.5009413 .
  * Sakai  _et al._ [2013] O. Sakai, J. Maeda, T. Shimomura, and K. Urabe, Functional composites of plasmas and metamaterials: Flexible waveguides, and variable attenuators with controllable phase shift, Physics of Plasmas 20, 073506 (2013), https://doi.org/10.1063/1.4816698 .
  * Wang and Cappelli [2016a] B. Wang and M. A. Cappelli, Waveguiding and bending modes in a plasma photonic crystal bandgap device, AIP Advances 6, 65015 (2016a).
  * Wang and Cappelli [2016b] B. Wang and M. A. Cappelli, A plasma photonic crystal bandgap device, Applied Physics Letters 108, 161101 (2016b).
  * Hagelaar and Pitchford [2005] G. J. M. Hagelaar and L. C. Pitchford, Solving the boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models, Plasma Sources Science and Technology 14, 722 (2005).
  * [33] SIGLO database.
  * Hilscher [2002] A. Hilscher, Determination of the cathode fall voltage in fluorescent lamps by measurement of the operating voltage, Journal of Physics D: Applied Physics 35, 1707 (2002).
  * Rodríguez [2023] J. A. Rodríguez,  _Electromagnetic Wave Manipulation with Plasma Metamaterials_ , Ph.D. thesis, Stanford University (2023).
  * Houriez  _et al._ [2022] L. S. Houriez, H. Mehrpour Bernety, J. A. Rodríguez, B. Wang, and M. A. Cappelli, Experimental study of electromagnetic wave scattering from a gyrotropic gaseous plasma column, Applied Physics Letters 120, 10.1063/5.0095038 (2022), 223101, https://pubs.aip.org/aip/apl/article-pdf/doi/10.1063/5.0095038/16449035/223101_1_online.pdf .
  * Mehrpour Bernety  _et al._ [2022] H. Mehrpour Bernety, L. S. Houriez, J. A. Rodríguez, B. Wang, and M. A. Cappelli, A characterization of plasma properties of a heterogeneous magnetized low pressure discharge column, AIP Advances 12, 10.1063/5.0124845 (2022), 115220, https://pubs.aip.org/aip/adv/article-pdf/doi/10.1063/5.0124845/16476461/115220_1_online.pdf .