# Complex Poynting vector in gyromagnetic media and its impact on power flow
in guided modes

Rajarshi Sen  Department of Electronics and Electrical Communication
Engineering, Indian Institute of Technology Kharagpur, Kharagpur, West Bengal,
India  Sarang Pendharker  Department of Electronics and Electrical
Communication Engineering, Indian Institute of Technology Kharagpur,
Kharagpur, West Bengal, India<EMAIL_ADDRESS>
###### Abstract

In this paper, we show the relation between the time-varying spinning nature
of the instantaneous Poynting vector and the phasor form of the complex
Poynting vector in the bulk of the gyromagnetic medium. We show the presence
of a transverse reactive power component in the bulk of the gyrotropic medium,
even for plane wave propagation. We use a simple quantification technique of
Poynting vector spin to analyze the rotation of the instantaneous Poynting
vector using the complex phasor form of the time-averaged Poynting vector. For
a transverse electric mode, we show the similarity between the Poynting vector
spin and the photonic spin of the magnetic field. The Poynting vector spin is
then used to represent the transverse power transfer across a ferrite-air
interface supporting TE surface wave modes. Considering a gyromagnetic
ferrite-filled rectangular waveguide following the TE mode profile, we show
the correspondence between the Poynting vector spin and the overall positive
and backward power flow. We analytically propose a mechanism to engineer the
region of the waveguide supporting backward and forward power propagation.

##  I Introduction

Gyrotropic materials form an important class of materials encompassing
naturally occurring and engineered materials. Gyrotropic materials can be
broadly classified into gyroelectric and gyromagnetic materials. Among these
two classes, gyromagnetic materials have seen widespread applications due to
their availability in the microwave frequency region. Microwave ferrites have
been extensively used to realize isolators[1, 2], circulators[3, 4, 5], and
phase shifters[6, 7, 8]. Apart from these conventional applications, they have
been recently used to realize nonreciprocity in numerous applications, such as
topological insulators[9, 10], thin films[11, 12], and mode-converting
waveguides[13, 14]. Despite the extensive use of gyromagnetic media in various
microwave and millimeter wave applications, an understanding of important
physical phenomena corresponding to such media was still missing. Researchers
consider isofrequency surfaces helpful in investigating the wave propagation
phenomena in complex media. These isofrequency surfaces have been used to
investigate wave propagation characteristics in both natural[15, 16, 17] as
well as engineered materials[18, 19, 20, 21].

We recently performed a thorough investigation of the nonreciprocal photonic
spin profile in the bulk of the media using isofrequency surfaces for
gyrotropic media[22]. This investigation widened our understanding of the
impact of gyrotropy over manipulating the plane wave propagation in the bulk
of the media and the existence of asymmetrical photonic spin profiles for
guided modes. However, a comprehensive understanding of power propagation in
gyrotropic media is still missing. We believe a detailed investigation of the
Poynting vector investigation will help us reveal key insights into the wave
propagation behavior in gyrotropic materials and provide governing principles
for engineering novel applications.

Poynting vector, which corresponds to the power flux density, has been used
extensively in different media and structures [23, 24, 25, 26, 27, 28, 29, 30,
31, 32] to observe the propagation of power. Conventionally, the Poynting
vector is represented using the instantaneous and time-averaged forms. The
time-averaged Poynting vector gives us an idea about the overall power flow in
the medium. The most general form of this time-averaged Poynting vector is its
complex form, consisting of both real and imaginary components, where the
imaginary component corresponds to the reactive power component. However, we
mostly consider only the real power component because, for most applications,
the investigation of real power itself is sufficient, such as reflection of
power from metasurfaces[33, 34], coupling of power for guided modes[35, 36,
37], and trapping of light for zero and negative Poynting vector[38, 39].
However, neglecting the imaginary (reactive) part of the complex Poynting
vector leads to an incomplete picture of the power flow in the medium.

It has been shown that the complex Poynting vector, which includes both real
and reactive power component plays an important role in the investigation of
propagating beams[40], the near field of Hertzian dipoles[41], evanescent
waves[42, 43] and optical forces over particles [40, 44, 45]. Recently, we
have seen research efforts being made to link this complex time-averaged
Poynting vector with the time-varying instantaneous Poynting vector. Litvin
has demonstrated the relationship between complex time-averaged Poynting
vector and the elliptical rotation of the instantaneous Poynting vector in the
case of Gaussian and Bessel beams[46]. Kim et al. made similar observations in
[47], where the longitudinal and transverse components of the time-averaged
Poynting vector are used to explain the spinning nature of the instantaneous
Poynting vector across an interface-supporting surface wave. This transverse
reactive power is crucial as it can explain the source of backward power
propagation observed in guided modes. Such investigation into the nature of
the Poynting vector in gyrotropic materials is missing.

In this paper, we show that gyrotropic material supports a complex Poynting
vector in the bulk for plane wave propagation. We further analyze the complex
Poynting vector with the help of isofrequency surfaces and link it with the
spinning nature of the instantaneous Poynting vector for plane waves in the
bulk of the medium. We quantify this rotation of the instantaneous Poynting
vector by introducing a Poynting vector spin. In the special case of
transverse electric mode, we show that the Poynting vector spin is directly
proportional to the photonic spin of the magnetic field
(Im($\vec{H}^{*}\times\vec{H}$)). We then use this Poynting vector spin to
explain the instantaneous transverse power flow across an interface between
gyromagnetic ferrite and air. Further, considering a ferrite-filled
rectangular waveguide, we show the correspondence between the backward
propagation of power density and the opposite sense of Poynting vector spin.
We provide guiding principles, using which the zero-power density point in the
cross-section of the waveguide can be engineered. The results presented in
this paper play a crucial role in understanding the nature of the
instantaneous Poynting vector in gyrotropic media and the corresponding
relation between the reactive and real components of the time-averaged
Poynting vector.

This paper is organized as follows. In section II, we analyze the power flow
phenomenon in the bulk of the gyromagnetic medium. Considering the transverse
electric mode, in section III, we quantify the temporal rotation of the
instantaneous Poynting vector. In section IV, using the understanding of the
temporal behavior of the Poynting vector for TE mode in gyromagnetic media, we
consider two different cases of surface wave propagation. We reveal the role
of transverse reactive power in supporting the unidirectional surface wave
propagation in gyromagnetic ferrite. Further, in section V, we report the
presence of negative power flow in a ferrite-filled waveguide, which is
validated using full-wave simulations. Our analytical expressions help in
engineering the Poynting vector distribution across the cross-section of the
waveguide. Finally, we conclude our results in section VI.

##  II Poynting vector in bulk of gyromagnetic media

Figure 1: (a) Graphical illustration of the temporal variation of the instantaneous Poynting vector (blue conical trace) in an elliptical manner around the time-averaged Poynting vector (red arrow). (b), (c) Isofrequency curves corresponding to the elliptical and hyperbolic regimes at frequencies 6 and 11 GHz, respectively. The isofrequency curves are color-mapped with the transverse reactive power. The arrows represent the time-averaged Poynting vector normal to the isofrequency surfaces. Table 1: YIG ferrite material specifications Parameter | Specification
---|---
Magnetic saturation (4$\pi M_{s}$) |  $1800$ Gauss
Magnetic bias ($H_{0}$) |  $3570$ Oersted
Dielectric permeability ($\epsilon_{f}$) | 14

For a comprehensive investigation of power flow over gyrotropic material, we
select YIG ferrite as our material of choice. The material specifications are
given in Table 1. Using the YIG ferrite material parameters, we can find the
permeability and gyrotropy terms, $\mu^{\prime}$ and $\kappa^{\prime}$,
respectively, corresponding to the operating frequency. We can express the
anisotropic permeability tensor for the $\hat{z}$-biased ferrite using these
anisotropy terms as

$\overset{\leftrightarrow}{\mu}_{r}=\begin{bmatrix}\mu^{\prime}&-j\kappa^{\prime}&0\\\
j\kappa^{\prime}&\mu^{\prime}&0\\\ 0&0&1\\\ \end{bmatrix}.$ (1)

Considering the uniaxial nature of gyromagnetic ferrite, we can simplify the
investigation by analyzing the isofrequency surfaces corresponding to a 2D
propagation scenario. However, the selected 2D plane should be parallel to the
bias-axis, which in this case is the $\hat{z}$-axis. Here, we select the $X-Z$
plane as the plane of wave propagation, and correspondingly, we will have wave
vector component $k_{y}=0$ and non-zero value of in-plane wave-vector
components $k_{x}$ and $k_{z}$. For this work, we exclude all the losses for
simplicity. Thus, the material parameters $\mu^{\prime}$ and
$\kappa^{\prime}$, and the wave-vector components $k_{x}$ and $k_{z}$ along
the isofrequency curves will be purely real quantities.

Once we have computed the isofrequency contour, we can find the field vector
intensities for the bulk of the medium using the method followed in [22].
Using the field intensities, we can finally compute both the instantaneous
pointing vector $\vec{P}_{ins}=\text{Re}(\vec{E})\times\text{Re}(\vec{H})$ and
time-averaged Poynting vector $\vec{P}_{avg}=0.5\vec{E}^{*}\times\vec{H}$ in
its Phasor form (see Appendix A).

$\vec{P}_{ins}=H_{0}^{2}k_{0}^{2}\Bigg{(}\Bigg{.}\Bigg{(}\Bigg{.}\frac{k_{x}(\epsilon_{f}k_{0}^{2}\kappa^{\prime
2}(k_{x}^{2}-\epsilon_{f}k_{0}^{2})^{2}+k_{z}^{2}(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})^{2})}{2\omega\epsilon_{0}(k_{x}^{2}-\epsilon_{f}k_{0}^{2})^{2}(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})^{2}}\\\
+\frac{k_{x}(-\epsilon_{f}k_{0}^{2}\kappa^{\prime
2}(k_{x}^{2}-\epsilon_{f}k_{0}^{2})^{2}+k_{z}^{2}(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})^{2})}{2\omega\epsilon_{0}(k_{x}^{2}-\epsilon_{f}k_{0}^{2})^{2}(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})^{2}}\cos(2\omega
t)\Bigg{.}\Bigg{)}\hat{x}-\frac{\kappa^{\prime}k_{x}(\epsilon_{f}k_{0}^{2}-k_{x}^{2}-k_{z}^{2})}{2\omega\epsilon_{0}(\epsilon_{f}k_{0}^{2}-k_{x}^{2})(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})}\sin(2\omega
t)\hat{y}\\\
+\Bigg{(}\Bigg{.}\frac{k_{z}((k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})^{2}+\epsilon_{f}k_{0}^{2}\kappa^{\prime
2}(\epsilon_{f}k_{0}^{2}-k_{x}^{2}))}{2\omega\epsilon_{0}(\epsilon_{f}k_{0}^{2}-k_{x}^{2})(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})^{2}}+\frac{k_{z}((k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})^{2}-\epsilon_{f}k_{0}^{2}\kappa^{\prime
2}(\epsilon_{f}k_{0}^{2}-k_{x}^{2}))}{2\omega\epsilon_{0}(\epsilon_{f}k_{0}^{2}-k_{x}^{2})(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})^{2}}\cos(2\omega
t)\Bigg{.}\Bigg{)}\hat{z}\Bigg{.}\Bigg{)}$ (2)
$\vec{P}_{avg}=H_{0}^{2}k_{0}^{2}\Bigg{(}\Bigg{.}\frac{k_{x}(\epsilon_{f}k_{0}^{2}\kappa^{\prime
2}(k_{x}^{2}-\epsilon_{f}k_{0}^{2})^{2}+k_{z}^{2}(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})^{2})}{2\omega\epsilon_{0}(k_{x}^{2}-\epsilon_{f}k_{0}^{2})^{2}(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})^{2}}\hat{x}\\\
+j\frac{\kappa^{\prime}k_{x}(\epsilon_{f}k_{0}^{2}-k_{x}^{2}-k_{z}^{2})}{2\omega\epsilon_{0}(\epsilon_{f}k_{0}^{2}-k_{x}^{2})(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})}\hat{y}+\frac{k_{z}(\epsilon_{f}k_{0}^{2}\kappa^{\prime
2}(\epsilon_{f}k_{0}^{2}-k_{x}^{2})+(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})^{2})}{2\omega\epsilon_{0}(\epsilon_{f}k_{0}^{2}-k_{x}^{2})(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime})^{2}}\hat{z}\Bigg{.}\Bigg{)}$
(3)

Here, $k_{0}$ denotes the free-space propagation constant. We observe that the
instantaneous Poynting vector for real power $\vec{P}_{inst}$ in the bulk YIG
ferrite contains both time-varying and time-invariant terms, which we can see
in eq. (2). Interestingly, the time-invariant terms exist only along the plane
of wave propagation, i.e., having $x$\- and $z$-components only. On the
contrary, the time-varying component of $\vec{P}_{inst}$ exists along the
$X-Z$ plane as well as transverse $\hat{y}$-direction. The transverse
component, however, is $90^{\circ}$ out of phase with respect to the in-plane
component, as indicated by their sine and cosine terms, respectively. This
phase difference between the orthogonal components leads to an elliptical
sweep of $\vec{P}_{inst}$ with respect to time. This spatio-temporal sweep of
$\vec{P}_{inst}$ is illustrated in Fig. 1(a). For simplicity, we have
considered only one direction of propagation having positive values of $k_{x}$
and $k_{z}$. The rotational phenomenon is depicted in this figure with the
help of a dashed green curve around the elliptical trace. Also, we have taken
four time-samples $t=0$, $t=T/8$, $t=T/4$, and $t=3T/8$, which correspond to
the phases $0$, $\pi/4$, $\pi/2$, and $3\pi/4$, respectively. The
$\vec{P}_{inst}$ at these four time-instances are shown along the elliptical
trace.

At this point, we make an important observation from eq. (2) and Fig. 1(a). On
time-averaging the instantaneous Poynting vector, the net real power component
along the transverse direction vanishes, and only a real component remains
along the plane of wave propagation. This non-zero real component of time-
averaged power is represented using a red arrow in Fig. 1(a). However, on
considering the complex time-averaged $\vec{P}_{avg}$ of eq. (3), we find that
there exists a purely imaginary transverse reactive component along
$\hat{y}$-axis. The magnitude of this transverse reactive power component of
$\vec{P}_{avg}$ in eq. (3) is similar to the magnitude of the transverse
component of $\vec{P}_{inst}$ in eq. (2), and is important in defining the
axis of the elliptical trace along $\hat{y}$-axis. Thus, we can say that, by
considering both the real and reactive power components of the complex
$\vec{P}_{avg}$, the complete time-varying behavior of $\vec{P}_{inst}$ can be
understood. Considering the chosen ferrite medium, we now compute the time-
averaged Poynting vector $\vec{P}_{avg}$ for two frequencies $f=6$ GHz and
$11$ GHz using eq. (3) while considering $H_{0}=1$ for simplicity. These two
frequencies correspond to the elliptical and hyperbolic regimes with positive
and negative values of $\mu^{\prime}$, respectively. The ($\mu^{\prime}$,
$\kappa^{\prime}$) values for these two frequencies are ($1.79$, $0.47$), and
($-1.38$, $-2.62$), respectively. In these two regimes, the values of
$\mu^{\prime}$ and $\kappa^{\prime}$ are suitable for supporting two
independent isofrequency surfaces[22]. Hence, with the help of these two
cases, we can completely observe the power flow mechanism in the bulk of
gyromagnetic media.

The real part of $\vec{P}_{avg}$ denoting the real power in the medium is
shown using arrows over the isofrequency curves, and the imaginary component
corresponding to the reactive power is represented using colormap in Fig. 1(b)
and (c) for these two frequencies, respectively. We can observe the direction
of the real component of $\vec{P}_{avg}$ is always normal to the isofrequency
curve. This geometrical relation of the time-average real-power propagation
being normal to the isofrequency surface can be verified throughout the
isofrequency surface covering all propagation directions. Further, we observe
that the time-average Poynting vector for closed isofrequency curves such as
ellipses diverges outward. Whereas, for open curves such as hyperbolas, the
Real time-average Poynting vector tends to converge along the transverse axis
of the hyperbola.

This phenomenon of time-averaged Poynting vector (which only has a real
component) being normal to the isofrequency surface has already been
investigated for non-gyrotropic anisotropic media [48, 49, 50]. However, for
gyrotropic material, the scenario is more complex with the presence of the
imaginary component of $\vec{P}_{avg}$, which accounts for the overall
reactive power in the bulk of the medium. The transverse reactive power
(aligned along $\hat{y}$-axis) for propagation along $X-Z$ plane changes its
sign for forward and backward propagation corresponding to $+\hat{x}-$ and
$-\hat{x}-$direction, respectively. In uniaxial media such as gyromagnetic
ferrite, considering the symmetry of the isofrequency surface around the bias
axis, the 2D analogous isofrequency curves portrayed the complete information.
Now, since the transverse reactive power follows opposite sign along the
opposite sides of the 2D isofrequency curve, the rotation of the isofrequency
curve to form a complete isofrequency surface will give us an impression of
the reactive power flux to encircle the bias axis of the gyromagnetic medium.
This intrinsic reactive power in plane waves in the bulk of the medium results
from the gyrotropy of the material. The imaginary reactive power component of
eq. (3) is linearly proportional to the gyrotropy term $\kappa^{\prime}$.
Thus, we can infer that in the absence of gyrotropy, the reactive power in the
bulk of the medium vanishes. In the next section, we shall quantify the
rotation of the instantaneous Poynting vector $\vec{P}_{inst}$.

##  III Quantification of Poynting vector rotation

The temporal behavior of $\vec{P}_{inst}$ has a strong relation with the
phasor form of complex time-averaged Poynting vector $\vec{P}_{avg}$.
Understanding the time-varying nature of $\vec{P}_{inst}$ in complex media
such as gyrotropic media is a challenge. To solve this challenge, we quantify
the rotation of $\vec{P}_{inst}$ using the definition of Poynting vector spin
$\vec{S}_{p}=\text{Im}(\vec{P}_{avg}^{*}\times\vec{P}_{avg})$. This method of
defining the Poynting vector spin is similar to the way photonic spin is
defined for the electric and magnetic fields as
$S_{e}=\text{Im}(\vec{E}^{*}\times\vec{E})$ and
$S_{h}=\text{Im}(\vec{H}^{*}\times\vec{H})$, respectively. Note that
$\vec{S}_{p}$ shall be non-zero only when $\vec{P}_{avg}$ is complex, which
signifies the quantification of the rotation of the instantaneous Poynting
vector $\vec{P}_{inst}$ using the real and reactive power components.

Let us now consider a special condition of wave propagation in the bulk of YIG
ferrite when the direction of propagation is along $\pm\hat{x}$-axis. Under
this condition, one of the propagation modes corresponds to the transverse
electric (TE) mode, supporting transverse electric $E_{z}$, transverse
magnetic $H_{y}$, and longitudinal magnetic field components $H_{x}$. In this
special case of TE propagation, the $\vec{S}_{p}$ becomes proportional to
$\vec{S}_{h}$ (see Appendix B), as electric field spin in this mode is absent.
This absence of electric field spin simplifies our investigation as the
computation of photonic spin corresponding to the magnetic field shall give us
the Poynting vector spin $\vec{S}_{p}$. We will now use this approach of
quantification of rotation of the Poynting vector to investigate the real and
reactive power propagation in ferrite-associated surface waves and waveguide
modes.

##  IV Poynting vector spin across ferrite-air interface

Figure 2: Surface wave propagation along (a-e) NPM-air and (f-j) ferrite-air
interfaces. (a,f) Schematic representation of the surface wave propagation,
(b,g) dispersion curves, (c,h) Re($E_{z}$) as colormap showing transverse
electric field component for the surface wave and traces of $\vec{P}_{inst}$,
(d,i) Real power component Re($\vec{{P}_{avg}}$), and (e,j) Poynting vector
spin $\vec{S}_{p}$

Surface wave modes propagate along the interface between two media while being
evanescent in the direction transverse to the interface. This surface wave
confinement is typically realized in an interface between a dielectric and a
medium supporting either negative permittivity[51, 52, 53] or permeability
value[54, 16]. In [47], Kim et al. considered surface wave propagation along
an interface between a dielectric and a negative permittivity medium. They
showed the relation between the opposite sense of rotation of $\vec{P}_{inst}$
and the direction of real power flow Re($\vec{P}_{avg}$) across the interface.
However, such negative permittivity materials are reciprocal. Gyromagnetic
ferrite supports negative $\mu_{eff}=(\mu^{\prime 2}-\kappa^{\prime
2})/\mu^{\prime}$ corresponding to the TE propagation mode for a particular
frequency region, which depends on the material parameters and applied bias.
We can facilitate highly nonreciprocal surface wave propagation along a
ferrite-air interface in this frequency range.

We now consider two separate interfaces for investigating power flow for
surface wave propagation: (a) negative permeability medium (NPM)-air and (b)
ferrite-air interface. The medium parameters considered for the NPM-air
interface are $\epsilon_{np}=1$, $\mu_{np}=-3$. Whereas, for the YIG ferrite,
we consider material parameters given in Table 1, except the $4\pi M_{s}$,
which is taken as $3750$ Gauss. The magnetic bias applied in the ferrite
material is along the $-\hat{z}$-axis. For this material configuration of the
YIG ferrite, the effective permeability in ferrite for TE mode of propagation
$\mu_{eff}$ becomes negative within a particular frequency band. We shall only
consider this special frequency band in this investigation of surface wave
modes along the ferrite-air interface.

The schematic representation for the surface wave propagation along the NPM-
air and ferrite-air interface is shown in 2(a) and (f), respectively. The
$X-Z$ plane is the plane of the interface at $y=0$. The waves will encounter
evanescence conditions along $|y|>0$. Next, we see the dispersion relation in
these two cases. First we find the field components $E_{z}$, $H_{x}$, and
$H_{z}$ corresponding to TE mode, The field components are assigned
$e^{\alpha_{yd}y+j(k_{x}x-\omega t)}$ and $e^{-\alpha_{yu}y+j(k_{x}x-\omega
t)}$ for $y\leq 0$ and $y\geq 0$, respectively. Here, $k_{x}$ is the
longitudinal propagation constant, $\alpha_{yd}$ is the attenuation constant
in the air, and $\alpha_{yu}$ is the attenuation constant for the NPM and
ferrite regions.

Next, we find the dispersion relations by applying the boundary conditions.
For the NPM-air interface, the dispersion curve is simpler, and we get a
closed-form expression for $k_{x}$ as

$k_{x}=\pm
k_{0}\sqrt{\frac{\mu_{a}\mu_{np}(\epsilon_{np}\mu_{a}-\epsilon_{a}\mu_{np})}{\mu_{a}^{2}-\mu_{np}^{2}}}$
(4)

whereas for the ferrite-air interface, we get the following characteristic
equation

$H_{0}k_{0}(k_{x}(\alpha_{yu}+\alpha_{yd}\mu^{\prime})+\kappa^{\prime}(\alpha_{yd}\alpha_{yu}-\epsilon_{f}k_{0}^{2}))=0$
(5)

where $\alpha_{yd}$ and $\alpha_{yu}$ are functions of $k_{x}$, $k_{0}$ and
the medium parameters (Appendix D). $k_{0}$ is the free space wavelength. We
now find the wave-vector value $k_{x}$ for the propagation mode along the NPM-
air interface and ferrite-air interface from eq. (4) and (5), respectively.
The dispersion curves for these two cases are shown in Fig. 2(b) and (g),
respectively. We observe that the dispersion relation for the NPM-air
interface is linear and reciprocal, whereas there is significant
nonreciprocity for the ferrite-air interface. Backward propagation ($k_{x}<0$)
exists only for a small frequency region, and most of this special frequency
region supports unidirectional surface wave propagation. We now select
frequencies $f_{1}=5$ GHz and $f_{2}=15.31$ GHz to investigate further the
Poynting vector corresponding to the NPM-air and ferrite-air interface,
respectively. These frequencies are marked in the dispersion curves using
dotted lines.

The electric field component $E_{z}$ for the NPM-air and ferrite-air
interfaces are shown in Fig. 2(c) and (h), respectively. We can observe strong
confinement along the interface for the surface mode in both cases. To
investigate the time-varying nature of $\vec{P}_{inst}$, we select five time-
steps $t=0$, $0.08T$, $0.16T$, $0.34T$, and $0.42T$. $\vec{P}_{inst}$
corresponding to these time-steps are overlayed over the electric field plot.
Interestingly, we observe a contrasting behavior for the sense of rotation of
$\vec{P}_{inst}$. The sense of rotation of the instantaneous Poynting vector
is opposite across the NPM-air interface. Whereas, for the ferrite-air
interface, we observe that $\vec{P}_{inst}$ rotation follows the same sense.

Using the expression for the field intensities, we find the analytical
expressions for real and reactive components of the time-averaged Poynting
vector. These real and imaginary power components for the NPN-air interface
are (Appendix C)

$\text{Re}(\vec{P}_{avg})=\begin{cases}\frac{k_{0}^{2}\mu_{a}(\alpha_{yd}^{2}+\epsilon_{a}k_{0}^{2}\mu_{a})e^{2\alpha_{yd}y}}{2\omega\epsilon_{0}k_{x}\alpha_{yd}^{2}}\hat{x}&\text{if
$y\leq 0$}\\\
\frac{k_{0}^{2}\mu_{np}(\alpha_{yu}^{2}+\epsilon_{np}k_{0}^{2}\mu_{np})e^{-2\alpha_{yu}y}}{2\omega\epsilon_{0}k_{x}\alpha_{yu}^{2}}\hat{x}&\text{if
$y\geq 0$}\end{cases}$ (6)
$\text{Im}(\vec{P}_{avg})=\begin{cases}-\frac{k_{0}^{2}\mu_{a}e^{2\alpha_{yd}y}}{2\omega\epsilon_{0}\alpha_{yd}}\hat{y}&\text{if
$y\leq 0$}\\\
\frac{k_{0}^{2}\mu_{np}e^{-2\alpha_{yu}y}}{2\omega\epsilon_{0}\alpha_{yu}}\hat{y}&\text{if
$y\geq 0$}\end{cases}$ (7)

and for the ferrite-air interface (Appendix D)

$\text{Re}(\vec{P}_{avg})=\begin{cases}\frac{k_{0}^{2}(\alpha_{ya}^{2}+\epsilon_{a}k_{0}^{2})e^{2\alpha_{ya}y}}{2\omega\epsilon_{0}k_{x}\alpha_{ya}^{2}}\hat{x}&\text{if
$y\leq 0$}\\\
\frac{k_{0}^{2}(\alpha_{yn}^{2}+\epsilon_{f}k_{0}^{2}\mu^{\prime})(\alpha_{yn}\kappa^{\prime}+k_{x}\mu^{\prime})e^{-2\alpha_{yn}y}}{2\omega\epsilon_{0}(\epsilon_{f}k_{0}^{2}\kappa^{\prime}-\alpha_{yn}k_{x})^{2}}\hat{x}\\\
&\text{if $y\geq 0$}\end{cases}$ (8)
$\text{Im}(\vec{P}_{avg})=\begin{cases}-\frac{k_{0}^{2}e^{2\alpha_{ya}y}}{2\omega\epsilon_{0}\alpha_{ya}}\hat{y}&\text{if
$y\leq 0$}\\\
-\frac{k_{0}^{2}(\alpha_{yn}\kappa^{\prime}+k_{x}\mu^{\prime})e^{-2\alpha_{yn}y}}{2\omega\epsilon_{0}(\epsilon_{f}k_{0}^{2}\kappa^{\prime}-\alpha_{yn}k_{x})^{2}}\hat{y}&\text{if
$y\geq 0$}\end{cases}$ (9)

From eq. (6) and (8), the real component of $\vec{P}_{avg}$ is completely
along $\hat{x}$-axis. This real component of $\vec{P}_{avg}$ is shown in Fig.
2(d) and (i) corresponding to the NPM-air and ferrite-air interfaces,
respectively. Additionally, for the NPM-air interface, the analytical results
are verified using CST Microwave Studio simulations. We observe that for the
NPM-air interface, the real power propagates along opposite directions on
either side of the interface. In contrast, the direction of real power
propagation remains consistent for the ferrite-air interface. This behavior of
real power flow being opposite in the two opposite sides across the NPM-air
interface and in the same direction for the ferrite-air interface is in
agreement with the corresponding sense of rotation of $\vec{P}_{inst}$. This
matching behavior of propagation direction of real power and rotational sense
of $\vec{P}_{inst}$ requires the computation of Poynting vector spin across
the NPM-air interface ($\vec{S}_{p,np}$) and ferrite-air interface
($\vec{S}_{p,f}$).

The analytical equations for $\vec{S}_{p,np}$ and $\vec{S}_{p,f}$ after
normalizing the magnetic field spatially are

$\vec{S}_{p,np}=\begin{cases}-\frac{2\alpha_{yd}k_{x}(\alpha_{yd}^{2}+\epsilon_{a}\mu_{a}k_{0}^{2})}{\alpha_{yd}^{4}+\epsilon_{a}^{2}\mu_{a}^{2}k_{0}^{4}+\alpha_{yd}^{2}(k_{x}^{2}+2\epsilon_{a}\mu_{a}k_{0}^{2})}\hat{z}&\text{if
$y\leq 0$}\\\
\frac{2\alpha_{yu}k_{x}(\alpha_{yu}^{2}+\epsilon_{np}\mu_{np}k_{0}^{2})}{\alpha_{yu}^{4}+\epsilon_{np}^{2}\mu_{np}^{2}k_{0}^{4}+\alpha_{yu}^{2}(k_{x}^{2}+2\epsilon_{np}\mu_{np}k_{0}^{2})}\hat{z}&\text{if
$y\geq 0$}\end{cases}$ (10)

$\vec{S}_{p,f}=\begin{cases}-\frac{2\alpha_{yd}k_{x}(\alpha_{yd}^{2}+\epsilon_{a}k_{0}^{2})}{\alpha_{yd}^{4}+\epsilon_{a}^{2}k_{0}^{4}+\alpha_{yd}^{2}(2\epsilon_{a}k_{0}^{2}+k_{x}^{2})}\hat{z}&\text{if
$y\leq 0$}\\\
-\frac{2(\epsilon_{f}k_{0}^{2}\kappa^{\prime}-\alpha_{yu}k_{x})(\alpha_{yu}^{2}+\epsilon_{f}k_{0}^{2}\mu^{\prime})}{\alpha_{yu}^{4}-2\alpha_{yu}\epsilon_{f}k_{0}^{2}\kappa^{\prime}k_{x}+\alpha_{yu}^{2}(k_{x}^{2}+2\epsilon_{f}k_{0}^{2}\mu^{\prime})+\epsilon_{f}^{2}k_{0}^{4}(\mu^{\prime
2}+\kappa^{\prime 2})}\hat{z}\\\ &\text{if $y\geq 0$}\end{cases}$ (11)

The subsequent plots of the Poynting vector spin $\vec{S}_{p,np}$ and
$\vec{S}_{p,f}$ across the NPM-air and ferrite-air interfaces are shown in
Fig. 2(e) and (j), respectively. As expected from the sense of rotation of
$\vec{P}_{inst}$, we observe the correspondence between the sign of
$\vec{S}_{p}$ and the direction of Re($\vec{P}_{avg}$).

For the ferrite-air interface in the specific frequency range having negative
$\mu_{eff}$, $\vec{S}_{p}$ maintains the same sign but a different magnitude
across either side of the interface. The inequality of the magnitude of
$\vec{S}_{p}$ corresponds to the discontinuity of the real power propagation
at the interface. Further, the $\vec{S}_{p}$ will also indicate the
propagation of imaginary reactive power across the interface. Previously, it
was observed that the material-induced photonic spin is locked to ferrite
material for bidirectional wave propagation[22]. Similarly, the Poynting
vector spin also remains locked in the ferrite medium irrespective of the
direction of propagation. In the next section, we consider a ferrite-filled
rectangular waveguide and observe this phenomenon of opposite signs of
$\vec{S}_{p}$ for forward and backward real power propagation in the
waveguide.

Figure 3: (a) Illustration of a $\hat{z}$-biased ferrite waveguide with cross-
section of dimension $a\times b$. The homogeneous bias $H_{0}$ is represented
using red arrows. The direction of propagation is along $\hat{x}$-axis. (b)
and (c) shows the $\vec{S}_{p}$ and Re($\vec{P}_{avg}$) along the cross-
section plane for TE${}_{\text{10}}$ mode. The plots include both analytical
and simulation results.

##  V Power flow dynamics in ferrite-filled waveguide

Figure 4: (a) 2D cross-section view of the waveguide depicting the arrangement
of the real power and the null-power point $y_{c}$. The variation of the null-
power point $y_{c}$ as a ratio with the half broad-side width $a/2$,
corresponding to (b) the variation of broad-side width $a$ and (c) frequency
of operation $f$. Both analytical and simulation results are shown.

Previously, we have performed a detailed study of the role of gyrotropy in
imposing a dominant photonic spin over the cross-section of a YIG ferrite-
filled waveguide (section II.C of [22]). The dominant photonic spin results
from the interaction of the material-induced and structure-induced photonic
spins. Here, we understand the resulting effect of the gyrotropy over the
power flow in the waveguide. Let us consider a rectangular waveguide of cross-
section $a=5$, $b=3$ mm, and length $l=30$ mm as shown in Fig. 3(a), which is
filled with $\hat{z}$-biased YIG ferrite (specifications similar to Table 1).
The operating frequency is selected as $7$ GHz. The waveguide walls are
assigned PEC boundary conditions to avoid losses. In this waveguide,
$\text{TE}_{\text{10}}$ is the dominant mode, and we shall select this mode
for our power-flow analysis. Correspondingly, the field variations for the
$\text{TE}_{\text{10}}$ mode shall be along the broad dimension of the
waveguide, i.e., along $\hat{y}$-axis. The field equations for this YIG
ferrite-filled waveguide in $\text{TE}_{\text{10}}$ mode of operation can be
written as

$E_{z}=\frac{A\omega\cos{k_{y}y}}{k_{y}}e^{jk_{x}x}$ (12)
$H_{x}=\frac{jA(\kappa^{\prime}k_{x}\cos{k_{y}y}-\mu^{\prime}k_{y}\sin{k_{y}y})}{k_{y}\mu_{0}(\kappa^{\prime
2}-\mu^{\prime 2})}e^{jk_{x}x}$ (13)
$H_{y}=\frac{A(\mu^{\prime}k_{x}\cos{k_{y}y}-\kappa^{\prime}k_{y}\sin{k_{y}y})}{k_{y}\mu_{0}(\kappa^{\prime
2}-\mu^{\prime 2})}e^{jk_{x}x}$ (14)

Similar to a dielectric-filled waveguide, we have $k_{y}=\pi/a$ in this
ferrite-filled waveguide. The normal component of magnetic flux density is
zero on the surface of the PEC wall, i.e., $\hat{y}\cdot\vec{B}=0$. This
condition would also make the normal component of the magnetic flux density
equal to zero along the waveguide wall in dielectric-filled waveguides.
However, in a gyromagnetic ferrite medium, due to the off-diagonal terms of
the permeability tensor, we have a non-zero normal component of magnetic field
intensity along the waveguide wall at $y=\pm a/2$. This non-zero component can
be observed in eq. (14) by putting $y=\pm a/2$. At $y=\pm a/2$ we have
$\sin{k_{y}y}=\pm 1$ and $\cos{k_{y}y}=0$. Hence, the non-zero $H_{x}$ and
$H_{y}$ at the waveguide wall correspond to non-zero Poynting vector spin
$\vec{S}_{p}$ along the waveguide wall. The mathematical expression for the
$\vec{S}_{p}$ along the $\hat{y}$-axis of the waveguide is

$\begin{split}\vec{S}_{p}=&\frac{2A^{2}}{k_{y}^{2}\mu_{0}^{2}(\mu^{\prime
2}-\kappa^{\prime 2})^{2}}(k_{y}\mu^{\prime}\sin k_{y}y\\\
&-k_{x}\kappa^{\prime}\cos k_{y}y)(k_{x}\mu^{\prime}\cos
k_{y}y-k_{y}\kappa^{\prime}\sin k_{y}y)\hat{z}\end{split}$ (15)

Using this equation of $\vec{S}_{p}$, we analytically compute the Poynting
vector spin across the waveguide cross-section, which is further validated
using CST Microwave Studio in Fig. 3(b). We observe a contrasting behavior of
the Poynting vector spin in this ferrite-filled waveguide with respect to the
conventional dielectric waveguide. In the latter case, an equal amount of
positive and negative $\vec{S}_{p}$ exists across the center of the cross-
section of the waveguide. Whereas, in this ferrite-filled waveguide, the
$\vec{S}_{p}$ crosses the zero-line in two places. The first cross-over point
is now offset from the center of the waveguide, and the second cross-over
point (non-existent in the case of a dielectric-filled waveguide) exists in
the other half of the waveguide. The analytical expressions for the first
($y_{a}$) and second cross-over ($y_{c}$) points along the cross-section are

$y_{a}=\frac{a}{\pi}\tan^{-1}\left(\frac{\kappa^{\prime}}{\mu^{\prime}}\sqrt{\frac{\epsilon_{f}\mu_{eff}k_{0}^{2}a^{2}}{\pi^{2}}-1}\right)$
(16)
$y_{c}=\frac{a}{\pi}\tan^{-1}\left(\frac{\mu^{\prime}}{\kappa^{\prime}}\sqrt{\frac{\epsilon_{f}\mu_{eff}k_{0}^{2}a^{2}}{\pi^{2}}-1}\right)$
(17)

From eq. (16) and (17), we can understand the role of gyrotropy over this
shift of the cross-over points. For a non-gyrotropic material, the gyrotropic
parameter $\kappa$ would be $0$, which will lead to $y_{a}=0$ and $y_{c}=a/2$,
similar to the case of a dielectric-filled waveguide.

These cross-over points having $\vec{S}_{p}=0$ exist due to the absence of
transverse reactive power and longitudinal real power at $y_{a}$ and $y_{c}$,
having $H_{x}=0$ and $H_{y}=0$, respectively. We are more interested in the
second cross-over point $y_{c}$ where the real power is absent, and beyond
which $\vec{S}_{p}$ makes a second transition into the negative region. This
point is marked using dotted lines in Fig. 3(b). Correspondingly, we plot the
real power propagation (Re$\vec{P}_{avg}$) along the forward $\hat{x}$-axis of
the waveguide using analytical equation and CST Microwave Studio in Fig. 3(c).
The cross-over point $y_{c}$ is marked using a dotted line. We observe that
the region supporting additional negative $\vec{S}_{p}$ beyond $y_{c}$
contains a backward propagation of the real power. This backward power is
absent for a dielectric-filled waveguide, as $y_{c}$ itself shall lie in the
outermost wall at $a/2$.

Using the expression of $y_{c}$ in eq. (17), we have the means to engineer
this cross-over point and alter the ratio of negative to positive propagation
of real power. In Fig. 4(a), the movement of the cross-over point $y_{c}$
along $\hat{y}$-axis is graphically illustrated. The variation of the broad-
width dimension $a$ and the frequency of operation $f$ changes $y_{c}$. At
this point, we have to note that the second cross-over point can be moved
leftward maximum up to the center of the waveguide $y_{c}=0$. At this point,
there is an equal amount of forward and backward propagation of real power,
and a cut-off condition is enforced in the waveguide. This gyrotropy-enforced
cut-off condition was first reported in Section IV(C) of [22]. Figure 4(b) and
(c) show the analytical and CST Microwave Studio plots of the ratio of $y_{c}$
with half the broad width dimension $a/2$ for variation of $a$ and $f$,
respectively. We can see significant movement of $y_{c}$ using variations of
both parameters. Note that this zero cross-over point will be located on the
other side of the waveguide for a backward propagating wave. The tunability of
the cross-over point $y_{c}$ will be beneficial for engineering nonreciprocal
gyrotropic guiding structures.

##  VI Conclusion

This paper establishes the strong relationship between the phasor form of the
time-averaged Poynting vector and its instantaneous form, which has an
intrinsic sense of rotation in a bulk gyrotropic media. We capture the sense
of rotation of the instantaneous Poynting vector using the Poynting vector
spin, which can be computed using the complex time-averaged Poynting vector.
Our results reveal the existence of reactive power components in the bulk of
gyrotropic material for plane waves, which can be represented aptly by the
Poynting vector spin. For the specific case of transverse electric mode
propagation, the Poynting vector spin is shown to be identical to the photonic
spin of the magnetic field component. The photonic spin of the magnetic field
is induced by the gyrotropy. Further, we investigate the TE-guided modes (i)
at the ferrite-air interface and (ii) in the ferrite-filled waveguide to
reveal the impact of Poynting vector spin on the positive and negative power
flow in these guiding structures. We show that the gyrotropy-induced asymmetry
and nonreciprocity in Poynting vector spin and the resultant negative power
flow in the ferrite-filled waveguide can be controlled by the waveguide and
gyrotropic parameters. This work provides the analysis using which new avenues
of designing gyrotropic material-based guiding-wave structures can open up.

##  Appendix A Derivation of Poynting vector for bulk ferrite

Using the mechanism used in [22], we find the magnetic field in the bulk of
ferrite

$\vec{H}_{bf}=H_{0}(\hat{x}+((j\epsilon_{f}k_{0}^{2}\kappa^{\prime})/(k_{x}^{2}+k_{z}^{2}-\epsilon_{f}k_{0}^{2}\mu^{\prime}))\hat{y}\\\
+((k_{x}k_{z})/(k_{x}^{2}-\epsilon_{f}k_{0}^{2}))\hat{z})e^{j(k_{x}x+k_{z}z-\omega
t)}$ (18)

The corresponding electric field vector can be found using Maxwell’s equation
as
$\vec{E}_{bf}=-(1/(\omega\epsilon_{0}\epsilon_{f})\overset{\leftrightarrow}{k}\cdot\vec{H}_{bf})$

Now we find the time-averaged Poynting vector $\vec{P}_{avg}$ as

$\vec{P}_{avg}=\frac{\vec{E}_{bf}^{*}\times\vec{H}_{bf}}{2}$ (19)

Using the instantaneous form of the magnetic field and electric field
component, we can find the instantaneous Poynting vector $\vec{P}_{inst}$ as

$\vec{P}_{inst}=\text{Re}(\vec{E}_{bf}e^{-j\omega
t})\times\text{Re}(\vec{H}_{bf}e^{-j\omega t})$ (20)

The real and imaginary components of equation 19 give us the time average real
and reactive power in the medium. In the instantaneous form, the time-varying
nature of real power flow in the medium can be analyzed using eq. 20.

##  Appendix B Poynting Vector and Photonic spin

The rotation of the field intensity vector ($\vec{H}$ or $\vec{E}$) can be
properly measured using the definition of Photonic spin and mathematically
expressed as $\vec{S}_{h}=\Im(\vec{H}^{*}\times\vec{H})$ for magnetic field
and $\vec{S}_{e}=\Im(\vec{E}^{*}\times\vec{E})$ for the electric field. This
definition of Photonic spin is analogous to computing the third Stoke’s
Parameter. Considering a transverse electric (TE) mode having field components
$E_{z}$, $H_{x}$, and $H_{y}$ propagating along $\hat{x}$ axis, the photonic
spin for the electric field is absent, and only $\vec{S}_{h}$ remains.
Further, the components of the time-averaged Poynting vector $\vec{P}_{avg}$
for this TE mode will exist only in the $\hat{x}$ $\hat{y}$ axis, the latter
being the reactive power direction. The time-averaged Poynting vector in
phasor form can be expressed as

$\vec{P}_{avg}=\frac{1}{2}(\vec{E}^{*}\times\vec{H})=P_{x}\hat{x}+P_{y}\hat{y}$
(21)

Now, following an approach similar to computing the photonic spin, we can
quantify the rotation of the instantaneous Poynting vector and call it
Poynting vector spin ($\vec{S}_{p}$). As we have real and imaginary component
of $\vec{P}_{avg}$ along $\hat{x}$ and $\hat{y}$-directions, respectively,
$\vec{S}_{p}$ can be defined as

$\begin{split}\vec{S}_{p}&=\Im(\vec{P}_{avg}^{*}\times\vec{P}_{avg})\\\
&=\Im(P_{x}^{*}P_{y}-P_{y}^{*}P_{x})\\\
&=0.25\Im(-E_{z}^{*}H_{y}^{*}E_{z}H_{x}+E_{z}^{*}H_{x}^{*}E_{z}H_{y})\\\
&=0.25|E_{z}|^{2}\Im(H_{x}^{*}H_{y}-H_{y}^{*}H_{x})\\\
&=0.25|E_{z}|^{2}\Im(\vec{H}^{*}\times\vec{H})\\\
&=0.25|E_{z}|^{2}\vec{S}_{h}\end{split}$ (22)

Since we only have non-zero $H_{x}$ and $H_{y}$ components, the photonic spin
$\vec{S}_{h}$ is along $\hat{z}$-axis, similarly, the Poynting vector spin
$\vec{S}_{p}$ is also along $\hat{z}$-axis.

We have the Poynting vector spin $\vec{S}_{p}$ related with the photonic spin
$\vec{S}_{h}$ as $\vec{S}_{p}=0.25|E_{z}|^{2}\vec{S}_{h}$. As the photonic
spin is normalized spatially, the term $|E_{z}|$ does not hold any significant
information, and we can safely relate $\vec{S}_{p}\sim\vec{S}_{h}$.

##  Appendix C Derivations for air-NPM interface

To support the wave confinement in the interface $y=0$, we define the TE mode
field components as

$E_{z}=\begin{cases}\frac{jH_{0}k_{0}^{2}\mu_{a}e^{\alpha_{yd}y}}{\omega\epsilon_{0}\alpha_{yd}}&\text{if
$y\leq 0$}\\\
-\frac{jH_{0}k_{0}^{2}\mu_{np}e^{-\alpha_{yu}y}}{\omega\epsilon_{0}\alpha_{yu}}&\text{if
$y\geq 0$}\end{cases}$ (23)

$H_{x}=\begin{cases}H_{0}e^{\alpha_{yd}y}&\text{if $y\leq 0$}\\\
H_{0}e^{-\alpha_{yu}y}&\text{if $y\geq 0$}\end{cases}$ (24)

$H_{y}=\begin{cases}-\frac{jH_{0}(\alpha_{yd}^{2}+\epsilon_{a}k_{0}^{2}\mu_{a})e^{\alpha_{yd}y}}{\alpha_{yd}k_{x}}&\text{if
$y\leq 0$}\\\
\frac{jH_{0}(\alpha_{yu}^{2}+\epsilon_{np}k_{0}^{2}\mu_{np})e^{-\alpha_{yu}y}}{\alpha_{yu}k_{x}}&\text{if
$y\geq 0$}\end{cases}$ (25)

For simplicity, we omitted the basis function $e^{k_{x}x-\omega t}$ from the
above field equations. On application of the boundary condition of $H_{x}$
being continuous at the interface $y=0$, we obtain the dispersion relation as

$k_{x}=\pm
k_{0}\sqrt{\frac{\mu_{a}\mu_{np}(\epsilon_{np}\mu_{a}-\epsilon_{a}\mu_{np})}{\mu_{a}^{2}-\mu_{np}^{2}}}$
(26)

Next, we compute the real ($0.5\Re(\vec{E}^{*}\times\vec{H})$) and reactive
power ($0.5\Im(\vec{E}^{*}\times\vec{H})$) using the field equations (eq.
23-23) as

$0.5\Re(\vec{E}^{*}\times\vec{H})=\begin{cases}\frac{H_{0}^{2}k_{0}^{2}\mu_{a}(\alpha_{yd}^{2}+\epsilon_{a}k_{0}^{2}\mu_{a})e^{2\alpha_{yd}y}}{2\omega\epsilon_{0}k_{x}\alpha_{yd}^{2}}\hat{x}&\text{if
$y\leq 0$}\\\
\frac{H_{0}^{2}k_{0}^{2}\mu_{np}(\alpha_{yu}^{2}+\epsilon_{np}k_{0}^{2}\mu_{np})e^{-2\alpha_{yu}y}}{2\omega\epsilon_{0}k_{x}\alpha_{yu}^{2}}\hat{x}&\text{if
$y\geq 0$}\end{cases}$ (27)
$0.5\Im(\vec{E}^{*}\times\vec{H})=\begin{cases}-\frac{H_{0}^{2}k_{0}^{2}\mu_{a}e^{2\alpha_{yd}y}}{2\omega\epsilon_{0}\alpha_{yd}}\hat{y}&\text{if
$y\leq 0$}\\\
\frac{H_{0}^{2}k_{0}^{2}\mu_{np}e^{-2\alpha_{yu}y}}{2\omega\epsilon_{0}\alpha_{yu}}\hat{y}&\text{if
$y\geq 0$}\end{cases}$ (28)

We can substitute $H_{0}=1$ A/m to simplify these equations of the time-
averaged Poynting vector. Further, we normalize the field component spatially
for the computation of photonic spin $\vec{S}_{h}$. Correspondingly, we can
compute the Poynting vector spin as
$\vec{S}_{p,np}=\Im(\vec{P}_{avg}^{*}\times\vec{P}_{avg})$.

##  Appendix D Derivation for ferrite-air interface

The field intensities corresponding to the TE mode of wave propagation along
the interface are defined as

$E_{z}=\begin{cases}\frac{jH_{0}k_{0}^{2}e^{\alpha_{yd}y}}{\omega\alpha_{yd}\epsilon_{0}}&\text{if
$y\leq 0$}\\\
\frac{jH_{0}k_{0}^{2}(\alpha_{yu}\kappa^{\prime}+k_{x}\mu^{\prime})e^{-\alpha_{yu}y}}{\omega\epsilon_{0}(\epsilon_{f}k_{0}^{2}\kappa^{\prime}-\alpha_{yu}k_{x})}&\text{if
$y\geq 0$}\end{cases}$ (29)

$H_{x}=\begin{cases}H_{0}e^{\alpha_{yd}y}&\text{if $y\leq 0$}\\\
H_{0}e^{-\alpha_{yu}y}&\text{if $y\geq 0$}\end{cases}$ (30)
$H_{y}=\begin{cases}-\frac{jH_{0}(\alpha_{yd}^{2}+\epsilon_{a}k_{0}^{2})e^{\alpha_{yd}y}}{\alpha_{yd}k_{x}}&\text{if
$y\leq 0$}\\\
-\frac{jH_{0}(\alpha_{yu}^{2}+\epsilon_{f}k_{0}^{2}\mu^{\prime})e^{-\alpha_{yu}y}}{\epsilon_{f}k_{0}^{2}\kappa^{\prime}-\alpha_{yu}k_{x}}&\text{if
$y\geq 0$}\end{cases}$ (31)

Using the boundary conditions of the tangential electric and magnetic field to
be continuous at the interface, we obtain the characteristic equation for the
solution of the propagation constant $k_{x}$ as

$H_{0}k_{0}(k_{x}(\alpha_{yu}+\alpha_{yd}\mu^{\prime})+\kappa^{\prime}(\alpha_{yd}\alpha_{yu}-\epsilon_{f}k_{0}^{2}))=0$
(32)

The terms $\alpha_{yd}$ and $\alpha_{yu}$ denotes the decay of the field
components along $-\hat{y}$ and $\hat{y}$ directions, respectively. These
attenuation constants can be mathematically defined as
$\alpha_{yd}=\sqrt{k_{x}^{2}-\epsilon_{a}k_{0}^{2}}$ and
$\alpha_{yu}=\sqrt{k_{x}^{2}-\epsilon_{f}\mu_{eff}k_{0}^{2}}$, where
$\mu_{eff}=(\mu^{\prime 2}-\kappa^{\prime 2})/\mu^{\prime}$. Further, we can
find the real and reactive power components using the field equations below.

$\Re(\vec{E}^{*}\times\vec{H})=\begin{cases}\frac{H_{0}^{2}k_{0}^{2}(\alpha_{yd}^{2}+\epsilon_{a}k_{0}^{2})e^{2\alpha_{yd}y}}{2\omega\epsilon_{0}k_{x}\alpha_{yd}^{2}}\hat{x}&\text{if
$y\leq 0$}\\\
\frac{H_{0}^{2}k_{0}^{2}(\alpha_{yu}^{2}+\epsilon_{f}k_{0}^{2}\mu^{\prime})(\alpha_{yu}\kappa^{\prime}+k_{x}\mu^{\prime})e^{-2\alpha_{yu}y}}{2\omega\epsilon_{0}(\epsilon_{f}k_{0}^{2}\kappa^{\prime}-\alpha_{yu}k_{x})^{2}}\hat{x}\\\
&\text{if $y\geq 0$}\end{cases}$ (33)
$\Im(\vec{E}^{*}\times\vec{H})=\begin{cases}-\frac{H_{0}^{2}k_{0}^{2}e^{2\alpha_{yd}y}}{2\omega\epsilon_{0}\alpha_{yd}}\hat{y}&\text{if
$y\leq 0$}\\\
-\frac{H_{0}^{2}k_{0}^{2}(\alpha_{yu}\kappa^{\prime}+k_{x}\mu^{\prime})e^{-2\alpha_{yu}y}}{2\omega\epsilon_{0}(\epsilon_{f}k_{0}^{2}\kappa^{\prime}-\alpha_{yu}k_{x})^{2}}\hat{x}&\text{if
$y\geq 0$}\end{cases}$ (34)

We can simplify the time-averaged Poynting vector equation by substituting
$H_{0}=1$ A/m. Subsequently, we compute the Poynting vector spin
$\vec{S}_{p,f}=\Im(\vec{P}_{avg}^{*}\times\vec{P}_{avg})$ after normalizing
the magnetic field vector $\vec{H}$ at every spatial coordinate.

## References

  * Ghaffar  _et al._ [2019] F. A. Ghaffar, J. R. Bray, M. Vaseem, L. Roy, and A. Shamim, IEEE Transactions on Magnetics 55, 1 (2019).
  * Noferesti and Djerafi [2021] M. Noferesti and T. Djerafi, IEEE Transactions on Magnetics 57, 1 (2021).
  * Marzall  _et al._ [2021] L. Marzall, D. Psychogiou, and Z. Popović, IEEE Transactions on Microwave Theory and Techniques 69, 1217 (2021).
  * Olivier  _et al._ [2020] V. Olivier, L. Huitema, B. Lenoir, H. Turki, C. Breuil, P. Pouliguen, and T. Monediere, IEEE Transactions on Microwave Theory and Techniques 68, 2521 (2020).
  * Chou  _et al._ [2018] H.-T. Chou, C.-H. Chang, and Y.-T. Chen, IEEE Transactions on Antennas and Propagation 66, 5934 (2018).
  * Ghaffar and Shamim [2015] F. A. Ghaffar and A. Shamim, IEEE Transactions on Magnetics 51, 1 (2015).
  * Nafe and Shamim [2015] A. Nafe and A. Shamim, IEEE Transactions on Microwave Theory and Techniques 63, 2264 (2015).
  * Kagita  _et al._ [2017] S. Kagita, A. Basu, and S. K. Koul, IEEE Transactions on Magnetics 53, 1 (2017).
  * Zhou  _et al._ [2020] P. Zhou, G.-G. Liu, X. Ren, Y. Yang, H. Xue, L. Bi, L. Deng, Y. Chong, and B. Zhang, Light: Science & Applications 9, 133 (2020).
  * Xi  _et al._ [2021] X. Xi, X.-M. Li, K.-P. Ye, H.-B. Wu, J. Chen, and R.-X. Wu, New Journal of Physics 23, 083042 (2021).
  * Tikhonov and Litvinenko [2019] V. V. Tikhonov and A. N. Litvinenko, Applied Physics Letters 115, 072410 (2019).
  * Tang  _et al._ [2021] S. Tang, M. S. Sarker, K. Ma, H. Yamahara, H. Tabata, and M. Seki, Applied Physics Letters 119, 082402 (2021).
  * Afshani and Wu [2019] A. Afshani and K. Wu, IEEE Transactions on Microwave Theory and Techniques 67, 3350 (2019).
  * Afshani and Wu [2021] A. Afshani and K. Wu, IEEE Transactions on Microwave Theory and Techniques 69, 4406 (2021).
  * Fu and Qin [2022] Y. Fu and H. Qin, Journal of Plasma Physics 88, 835880401 (2022).
  * Macêdo and Camley [2019] R. Macêdo and R. E. Camley, Phys. Rev. B 99, 014437 (2019).
  * Lokk [2017] E. G. Lokk, Journal of Communications Technology and Electronics 62, 251 (2017).
  * Sakhno  _et al._ [2022] D. Sakhno, E. Koreshin, and P. Belov, Opt. Lett. 47, 4451 (2022).
  * Yu  _et al._ [2016] K. Yu, Z. Guo, H. Jiang, and H. Chen, Journal of Applied Physics 119, 203102 (2016), https://pubs.aip.org/aip/jap/article-pdf/doi/10.1063/1.4952378/13046001/203102_1_online.pdf .
  * Papadakis  _et al._ [2018] G. T. Papadakis, D. Fleischman, A. Davoyan, P. Yeh, and H. A. Atwater, Nature Communications 9, 296 (2018).
  * Fedorin [2023] I. Fedorin, Scientific Reports 13, 16151 (2023).
  * Sen and Pendharker [2022] R. Sen and S. Pendharker, Phys. Rev. A 105, 023528 (2022).
  * Warmbier  _et al._ [2012] R. Warmbier, G. S. Manyali, and A. Quandt, Phys. Rev. B 85, 085442 (2012).
  * Liu  _et al._ [2004] Z. Liu, Z. Lin, and S. T. Chui, Phys. Rev. B 69, 115402 (2004).
  * Belov [2003] P. A. Belov, Microwave and Optical Technology Letters 37, 259 (2003), https://onlinelibrary.wiley.com/doi/pdf/10.1002/mop.10887 .
  * Cheng and Cui [2006] Q. Cheng and T. J. Cui, Phys. Rev. B 73, 113104 (2006).
  * Costa  _et al._ [2011] J. a. T. Costa, M. G. Silveirinha, and A. Alù, Phys. Rev. B 83, 165120 (2011).
  * Popov  _et al._ [2006] A. K. Popov, V. V. Slabko, and V. M. Shalaev, Laser Physics Letters 3, 293 (2006).
  * Mackay and Lakhtakia [2009] T. G. Mackay and A. Lakhtakia, Phys. Rev. B 79, 235121 (2009).
  * Rondón-Ojeda and Soto-Eguibar [2018] I. Rondón-Ojeda and F. Soto-Eguibar, Wave Motion 78, 176 (2018).
  * Halterman  _et al._ [2007] K. Halterman, S. Feng, and P. L. Overfelt, Phys. Rev. A 76, 013834 (2007).
  * Edwards  _et al._ [2008] B. Edwards, A. Alù, M. E. Young, M. Silveirinha, and N. Engheta, Phys. Rev. Lett. 100, 033903 (2008).
  * Taghvaee  _et al._ [2021] H. Taghvaee, F. Liu, A. Díaz-Rubio, and S. Tretyakov, Phys. Rev. B 104, 235409 (2021).
  * Kwon and Tretyakov [2017] D.-H. Kwon and S. A. Tretyakov, Phys. Rev. B 96, 085438 (2017).
  * Fang  _et al._ [2015] Y.-t. Fang, H.-Q. He, J.-x. Hu, L.-k. Chen, and Z. Wen, Phys. Rev. A 91, 033827 (2015).
  * Kodera  _et al._ [2013] T. Kodera, D. L. Sounas, and C. Caloz, IEEE Transactions on Microwave Theory and Techniques 61, 1030 (2013).
  * Suzuki  _et al._ [2008] E. Suzuki, S. Arakawa, M. Takahashi, H. Ota, K. I. Arai, and R. Sato, IEEE Transactions on Instrumentation and Measurement 57, 1014 (2008).
  * Huang  _et al._ [2008] Y. J. Huang, W. T. Lu, and S. Sridhar, Phys. Rev. A 77, 063836 (2008).
  * Peacock and Broderick [2003] A. C. Peacock and N. G. R. Broderick, Opt. Express 11, 2502 (2003).
  * Xu and Nieto-Vesperinas [2019] X. Xu and M. Nieto-Vesperinas, Phys. Rev. Lett. 123, 233902 (2019).
  * Baumgartner  _et al._ [2022] P. Baumgartner, A. Masiero, C. Riener, and T. Bauernfeind, Electronics 11, 10.3390/electronics11131967 (2022).
  * Nieto-Vesperinas and Xu [2021] M. Nieto-Vesperinas and X. Xu, Phys. Rev. Res. 3, 043080 (2021).
  * Picardi  _et al._ [2019] M. F. Picardi, M. Neugebauer, J. S. Eismann, G. Leuchs, P. Banzer, F. J. Rodríguez-Fortuño, and A. V. Zayats, Light: Science & Applications 8, 52 (2019).
  * Antognozzi  _et al._ [2016] M. Antognozzi, C. R. Bermingham, R. L. Harniman, S. Simpson, J. Senior, R. Hayward, H. Hoerber, M. R. Dennis, A. Y. Bekshaev, K. Y. Bliokh, and F. Nori, Nature Physics 12, 731 (2016).
  * Nieto-Vesperinas and Xu [2022] M. Nieto-Vesperinas and X. Xu, Light: Science & Applications 11, 297 (2022).
  * Litvin [2012] I. A. Litvin, J. Opt. Soc. Am. A 29, 901 (2012).
  * Kim  _et al._ [2012] K.-Y. Kim, J. Kim, I.-M. Lee, and B. Lee, Phys. Rev. A 85, 023840 (2012).
  * Kong [2008] J. A. Kong,  _Electromagnetic wave theory_ (EMW Publishing, 2008).
  * Schilling [2006] J. Schilling, Phys. Rev. E 74, 046618 (2006).
  * Cortes and Jacob [2013] C. L. Cortes and Z. Jacob, Phys. Rev. B 88, 045407 (2013).
  * Shalin  _et al._ [2015] A. S. Shalin, S. V. Sukhov, A. A. Bogdanov, P. A. Belov, and P. Ginzburg, Phys. Rev. A 91, 063830 (2015).
  * Gao  _et al._ [2015] W.-L. Gao, F.-Z. Fang, Y.-M. Liu, and S. Zhang, Light: Science & Applications 4, e328 (2015).
  * Kotov and Lozovik [2016] O. V. Kotov and Y. E. Lozovik, Phys. Rev. B 93, 235417 (2016).
  * Li  _et al._ [2021] M. Li, J. Song, and Y. Jiang, Phys. Rev. B 103, 045307 (2021).


# A Good Snowman is Hard to Plan

Miquel Bofill,1 Cristina Borralleras, 2 Joan Espasa, 3 Gerard Martín,1 Gustavo
Patow,1 Mateu Villaret1

###### Abstract

In this work we face a challenging puzzle video game: _A Good Snowman is Hard
to Build_. The objective of the game is to build snowmen by moving and
stacking snowballs on a discrete grid. For the sake of player engagement with
the game, it is interesting to avoid that a player finds a much easier
solution than the one the designer expected. Therefore, having tools that are
able to certify the optimality of solutions is crucial.

Although the game can be stated as a planning problem and can be naturally
modelled in PDDL, we show that a direct translation to SAT clearly outperforms
off-the-shelf state-of-the-art planners. As we show, this is mainly due to the
fact that reachability properties can be easily modelled in SAT, allowing for
shorter plans, whereas using axioms to express a reachability derived
predicate in PDDL does not result in any significant reduction of solving time
with the considered planners. We deal with a set of 51 levels, both original
and crafted, solving 43 and with 8 challenging instances still remaining to be
solved.

## Introduction

During the design of video games, one of the most important aspects to
consider is the difficulty of each level. A common problem while designing the
levels is the possibility of ignoring a possible solution to the level that is
much easier than the ones the designer had initially in mind (Silli 2010).
These unplanned solutions often are not desired because they break the
difficulty slope, making the game uninteresting. To this end, it is of crucial
interest to have a tool able to provide (certified) optimal solutions for
given scenarios of the game.

In this work we focus on finding optimal solutions for a discrete time and
space puzzle. In particular we consider the PSPACE-complete game _A Good
Snowman is Hard to Build_ (Davis and Hazelden 2015; He, Liu, and Yang 2017),
where a character has to move snowballs to build snowmen. The game has some
resemblance with the widely studied Sokoban (Culberson 1997) game, because the
character pushes balls (instead of boxes) in a matrix-like maze with some
obstacles. However, it generally gets more involved because not only scenario
characteristics (positions of balls and snow) may change when moving
snowballs, but the final positions of snowmen are not fixed, hence the goal is
conceptually disjunctive, making the search space much bigger in comparison.
Therefore, both the modelling and the increasing computational complexity
makes the considered game challenging. In these kind of games, the difficulty
of a level is sometimes measured using the number of times the character needs
to move an object. There are two main reasons for this: first, moving the
character from one location to another uses to be trivial because it consists
on identifying if there is a path from one location to another; second,
pushing an object may close free paths.

_A Good Snowman is Hard to Build_ is naturally characterised as a planning
problem, where we are asked to find a sequence of actions such as rolling
balls on snow, and stacking or unstacking them, so that balls of the
appropriate sizes are joined into snowmen. In the planning literature, most
models for Sokoban have two actions: _move_ , for changing the location of the
character, and _push_ for changing the location of a box by making the
character push it. Some previous works devise methods to focus the search
efforts on certain parts of the problem. In Ivankovic and Haslum (2015), the
authors show how a reachability derived predicate can be specified with axioms
and the _move_ actions can be completely avoided in the model. The idea is to
ensure with the reachability predicate that the pushing location is
_reachable_ from the current character location. This results in shorter plans
and in a significant reduction of the search space and hence in the time
required to solve the instances. Similarly, a proposal to automatise the
inference of axioms and to reformulate the problem accordingly is given in
Miura and Fukunaga (2017). This is done with success for Sokoban and
reachability, as shown by using axiom supporting model-based planners with
Integer programming and Answer Set Programming technologies.

In this paper, we start by providing a model for _A Good Snowman is Hard to
Build_ written in PDDL (Haslum et al. 2019), the de-facto standard language
for AI planners. This model contains actions to both _move_ the character and
to _push_ snowballs, and we evaluate the performance of some state-of-the-art
planners on it. We also provide a PDDL model using axioms for the reachability
predicate, hence, only considering _push_ actions, which preconditions require
that the location where the action is performed is reachable from the location
where the character was at the previous time step. Unfortunately, despite the
detailed instructions for the planner from Ivankovic and Haslum (2015) we have
not been able to install it, and the planner from Miura and Fukunaga (2017)
was not available. Therefore, we use a state-of-the-art planner supporting
derived predicates but without any specialized heuristics for them. To fully
assess the overhead of the _move_ actions, we also consider a “what if”
scenario, where the character can do unsound teleports. That is, teleports
that do not guarantee that there exists a valid path from source to
destination (a sort of cheating). Interestingly, experiments show that this
does not have a significant impact in solving time.

Next, we present an ad-hoc planning as SAT approach. Again, we then get rid of
_move_ actions by encoding reachability. This results in a significant
reduction of the time horizon and, accordingly, of the solving times, allowing
to solve many more instances to optimality. By using the proposed planning as
SAT approach with reachability, we are able to certify the minimum number of
ball movements needed to solve most of the original levels of _A Good Snowman
is Hard to Build_ and of the additional crafted instances that we also
provide.

The main contribution of this paper is a new case study with a set of
challenging PDDL instances (with and without derived predicates) and a
detailed comparison of planning tools against planning as SAT methodology. For
the planning as SAT approach we use the s-t-reachability encoding from Gebser,
Janhunen, and Rintanen (2014); Rankooh and Rintanen (2021) and the folklore
one based on neighbour counting. We adapt these encodings to deal with a
priori unknown source and target locations in a dynamic environment.

For the sake of reproducibility, the presented PDDL models and instances and
the SAT encoding generators can be found in the supplementary material at
https://github.com/udg-lai/KEPS2023. In addition, the full detailed
experimental results are also included.

## The Game

_A Good Snowman Is Hard To Build_ is a single-agent puzzle video game where
the goal is to push snowballs in a maze to build some snowmen by stacking
three balls of decreasing size.

The game elements are the agent (i.e., the black character controlled by the
player), the playable cells, which may or not contain snow, and the snowballs,
which are initially distributed on the playable cells. Snowballs have three
possible sizes: small, medium and large. As in Sokoban, the only allowed
action is moving the agent in one of four directions. The results of _moving_
depend on the cells in front:

  * •

_Move_ : When the agent walks into a free cell, he simply moves to that cell.

  * •

_Roll_ : When the agent moves into a cell with a single ball, and there is a
free cell in front of the ball, the ball gets pushed and the agent occupies
the cell previously occupied by the ball. If a ball is pushed into a snow
cell, the snow disappears and the ball increases in size, up to a maximum.
Still, the snow is always removed.

  * •

_Push_ : A ball can be pushed on a stack of balls if the size of the ball in
the top is bigger than the one being pushed. Then, the agent occupies the cell
previously occupied by the pushed ball like when rolling.

  * •

_Pop_ : Trying to move into a stack of balls will pop the topmost ball but
will not change the location of the agent. This action can only happen if the
ball falls directly into a cell without any ball.

Notice that, as in Sokoban, the agent is not allowed to pull balls.

The goal of the game is to build snowmen composed by a pile of three balls of
decreasing size. The scenarios of the game considered consist of three, six or
nine snowballs, hence, one, two or three snowmen. Snowmen can be built
anywhere. Figure 1 depicts an example of how to solve one of the levels of the
game.

Figure 1: _Andy_ level, showing the execution of the optimal solution
_lluRurDlldddrUluRuurrrdLulD_. Letters represent the direction of movement.
Uppercase letters indicate ball movements.

## PDDL Formulation

Since the problem is very similar to Sokoban, we consider a PDDL formulation
based on the Sokoban domain used in the 2008 and 2011 International Planning
Competitions. However, when encoding it in PDDL there are various key
differences, as the usage of quantifiers and conditional effects is essential.
This is due to elements such as: disappearing snow, ball stacking and the
disjunctive goal. In this section we partially describe this PDDL
formulation.111Full models and instances can be found at
https://github.com/udg-lai/KEPS2023.

#### Types, objects and fluents

The objects considered have the following types: loc, dir and ball. The loc
type represents a grid cell location, which may contain the character or a
stack of balls. In addition to these, a cell may also contain snow. The dir
type enumerates the four directions the character and balls can move (up,
down, left, right) and the ball type defines all possible balls of the
scenario, e.g., ball1, ball2 and ball3. Note that when building two snowmen
there will be six balls, whereas when building three snowmen there will be
nine balls. In order to model the current state, the following predicates are
defined:

  * •

(snow ?l - loc): location l is covered in snow.

  * •

(next ?from ?to - loc ?d - dir): locations from and to are are adjacent;
location to is in direction d relative to from.

  * •

(occ ?l - loc): location l contains at least one ball.

  * •

(char_at ?l - loc): character is at location l.

  * •

(ball_at ?b - ball ?l - loc): ball b is at location l.

  * •

(ball_size_s ?b - ball): ball b has small size.

  * •

(ball_size_m ?b - ball): ball b has medium size.

  * •

(ball_size_l ?b - ball): ball b has large size.

  * •

(goal): the goal is satisfied.

#### Actions

ppos from to

Figure 2: Location parameters of action move_ball.

Three actions are defined: move_character, move_ball and goal. As its name
indicates, action move_character moves the character from its current location
to some of its next adjacent locations, according to a direction.

Action move_ball moves a single ball. This action has as parameters a ball to
be moved (b) and three aligned and adjacent locations (ppos, from and to).
Figure 2 illustrates the role of the three location parameters. With this
action, ball b always is going to move, but the character may not.

⬇

1(:action move_ball

2 :parameters

3 (?b - ball ?ppos ?from ?to - loc ?d - dir)

4 :precondition (and

5 (next ?ppos ?from ?d) (next ?from ?to ?d)

6 (ball_at ?b ?from)

7 (char_at ?ppos)

8 (forall (?o - ball)

9 (or

10 (= ?o ?b) ; Ball != from ball action

11 (or ; Implication

12 (not (ball_at ?o ?from))

13 (or ; Ball ?b smaller than ?o

14 (and

15 (ball_size_small ?b)

16 (ball_size_medium ?o))

17 (and

18 (ball_size_small ?b)

19 (ball_size_large ?o))

20 (and

21 (ball_size_medium ?b)

22 (ball_size_large ?o))))))

23 (or

24 (forall (?o - ball)

25 (or (= ?o ?b)

26 (not (ball_at ?o ?from))))

27 (forall (?o - ball)

28 (not (ball_at ?o ?to))))

29 (forall (?o - ball)

30 (or ; Implication

31 (not (ball_at ?o ?to))

32 (or ... ; Ball ?b smaller than ?o

33 ))

34 :effect (and

35 (occ ?to)

36 (not (ball_at ?b ?from)) (ball_at ?b ?to)

37 (when

38 (forall (?o - ball)

39 (or (= ?o ?b)

40 (not (ball_at ?o ?from))))

41 (and (not (char_at ?ppos))

42 (char_at ?from)

43 (not (occ ?from))))

44 (not (snow ?to))

45 (when

46 (and (snow ?to)

47 (ball_size_s ?b))

48 (and (not (ball_size_s ?b))

49 (ball_size_m ?b)))

50 (when

51 (and (snow ?to)

52 (ball_size_m ?b))

53 (and (not (ball_size_m ?b))

54 (ball_size_l ?b)))

55 (increase (total-cost) 1)))

To perform the action we first need to ensure that the given three locations
are aligned and adjacent (line 5), ball b is located at the from location
(line 6), and the character is located at the ppos location (line 7). As the
game rules states, we can only move a ball if it does not have any other ball
on top, i.e., if it is the smallest ball of the stack (lines 8–22). Also, if a
ball is on top of a stack of balls, it cannot be moved to another stack of
balls (lines 23–28). Finally, to maintain a correct game state, a ball can
only be pushed to a stack of balls if it is smaller than all other balls in
the stack (lines 29–33).

When the action is performed, since a ball is moved, it modifies the occupancy
of some location (lines 35 and 37–43). Note that occupancy is conditional
since the character might be rolling, pushing or popping a ball. We also have
to update the ball location (line 36). If the character rolls or pushes a
ball, its location changes (lines 37–43). We remove the snow from the ball
location unconditionally (line 44). If there is snow in the target location,
we consider two cases: increasing ball size from small to medium (lines 45-49)
or from medium to large (lines 50-54). At last, we declare that this action
has cost 1 (line 55).

The goal action (in the single snowman scenarios) has three balls and a
location as parameters. Via its precondition it enforces that the three balls
are different objects located at the same place. Note that the ball sizes
constraint is enforced in the possible state transitions of the move_ball
action.

⬇

1(:action goal

2 :parameters

3 (?b0 ?b1 ?b2 - ball ?p0 - loc)

4 :precondition (and

5 (not (= ?b0 ?b1)) (not (= ?b0 ?b2))

6 (not (= ?b1 ?b2))

7 (ball_at ?b0 ?p0) (ball_at ?b1 ?p0)

8 (ball_at ?b2 ?p0))

9 :effect (and (goal)) )

As said in the introduction, in Sokoban-like games, sometimes the notion of
optimal plans prioritizes box movements minimization. We propose to also
consider this notion of optimality here, hence, correspondingly, only snow
ball movements are considered. This is meaningful in Snowman because, for a
human player, it is trivial to identify if there is a path from the current
location of the character to the location of the next snow ball movement.
Therefore, we have the statement (:metric minimize (total-cost)) and only the
action move_ball has a cost.

#### Cheating variant

In order to evaluate how costly is having to find out all the movements of the
character, we also evaluate a modified formulation, where the move_character
action is removed and in the move_ball action it is assumed that the character
can teleport (see the empirical evaluation section). Note that this relaxation
of the formulation may lead to non-valid plans since sometimes balls will be
blocking some paths.

#### Reachability variant

PDDL offers logical axioms to specify derived predicates, a feature able to
cleanly express reachability. Similarly as done in Ivankovic and Haslum (2015)
with Sokoban, we declare a new predicate reachable with two parameters, the
source and target location. Given a pair of ?s (source) and ?t (target)
locations, the target will be reachable from the source if the target is not
occupied and, it is either the same location as the source or there exists a
middle location ?m such that it is next to the source and it is able to reach
the target.

⬇

1(:derived (reachable ?s ?t - location)

2 (and

3 (or ;; there is a path between ?s and ?t

4 (= ?s ?t) ;; either is the same location,

5 ;; or it is inductively reachable

6 (exists (?d - direction)

7 (exists (?m - location)

8 (and

9 (next ?s ?m ?d)

10 (not (occupancy ?m))

11 (reachable ?m ?t)))))

12 ;; and target location has no ball in it

13 (not (occupancy ?t))))

We can now remove the action that moves the character. We need an additional
parameter ?prevpos in the move_ball action, to refer to the position where the
character ended in previous state. We then require (reachable ?prevpos ?ppos)
in move_ball action precondition, enforcing that the position where the ball
push is going to be executed from is reachable from where the character was.
Finally, an additional effect will be required to update the character
position, i.e., change the position of the character from ?prevpos to ?ppos.

## Planning as SAT

Here we recall some basic ideas of Planning as SAT. In this approach (Kautz
and Selman 1992), a planning problem is encoded to a Boolean formula, with the
property that any model of this formula will correspond to a valid plan.

Since the length of a valid plan is not known a priori, the basic idea is to
encode the existence of a plan of $T$ steps with a formula $f(T)$. Then, the
method for finding the shortest plan consists in iteratively checking the
satisfiability of $f(T)$ for $T=0,1,2,\dots$ until a satisfiable formula is
found.

Variables need to be replicated for each time step. E.g., $a^{t}$ will denote
if action $a$ is executed at time $t$. Then, the general (standard) encoding
goes as follows. First of all, it is stated that the execution of an action
implies its preconditions, for each $t\in 0..T-1$:

$a^{t}\rightarrow\mathit{Pre}^{t}\text{ for every action
}a=\langle\mathit{Pre},\mathit{Eff}\rangle$

Also, if an action is executed, its effects must take place at the next time
step:

$a^{t}\rightarrow\mathit{Eff}^{t+1}\text{ for every action
}a=\langle\mathit{Pre},\mathit{Eff}\rangle$

Preconditions and effects are typically given as sets of literals. Here, by
$\mathit{Pre}^{t}$ and $\mathit{Eff}^{t+1}$ we denote the corresponding
formulas (usually, conjunctions of literals) on the time-indexed state
variables.

Moreover, a change in the value of every state variable $v$ can occur only if
an action that can change this value has been executed:

$(v^{t}\neq v^{t+1})\rightarrow\bigvee\\{a^{t}\mid
a=\langle\mathit{Pre},\mathit{Eff}\rangle,v\in\mathit{Eff}\\}$

These are called frame axioms. Finally, it is stated that exactly one action
is executed at each time step $t$, and that the goal holds at time $T$.

## Snowman SAT encoding

In this section we propose an encoding for solving the problem at hand,
following the planning as SAT approach.

For clarity and space limitations we do not provide the whole encoding, but
the viewpoint (state variables) and an excerpt of the formulas including the
goal and relevant transition constraints and frame axioms. Also, for
simplicity reasons, we only describe the case of a single snowman.

We represent states with Boolean variables stating, for each location, whether
(i) there is snow or not, (ii) there is a ball of a particular size or not,
and (iii) there is the character or not. The actions considered are only four,
corresponding to a ball movement per possible direction. This will eventually
result in rolling, pushing or popping a ball, depending on the current state.
In other words, we only need to know the direction of the action.

We consider $L$ to be the set of valid (non-wall) locations, $S$ the set of
locations with snow, and $T$ the number of time steps considered. Below we
detail the main variables and constraints.

#### Variables

For all $l$ in $L$, $t$ in $0..T$:

  * •

$s_{l}^{t}$: there is snow at location $l$ at time $t$

  * •

$bs_{l}^{t},bm_{l}^{t},bl_{l}^{t}$: there is a small, medium or large ball at
location $l$ at time $t$

  * •

$c_{l}^{t}$: character is at location $l$ at time $t$

For all $t$ in $0..T-1$:

  * •

$n^{t},s^{t},e^{t},w^{t}$: direction of action at time $t$ is north, south,
east or west

#### Constraints

The goal consists in requiring all balls to be at the same location at the end
of the plan:

$\displaystyle\forall l\in L\quad(bs_{l}^{T}\leftrightarrow bm_{l}^{T}\wedge
bm_{l}^{T}\leftrightarrow bl_{l}^{T})$

For each time step $t$ in $0..T-1$ we have the following constraints.

  * •

Exactly one action is executed per time step: we impose an exactly-one
constraint over $n^{t},s^{t},e^{t},w^{t}$ variables.

  * •

Action preconditions and effects (excerpt of the action to move the character
north):

Let $L_{n}$ be the set of valid locations with a wall at north and let
$L_{nn}$ be the set of valid locations with a wall two locations ahead at
north.

When the character is at any location in $L_{n}$, it cannot go north:

$\displaystyle\forall l\in L_{n}\quad c_{l}^{t}\to\neg n^{t}$

Otherwise, if the character moves north and it has a wall two locations ahead,
there cannot be any ball in front of him, and at the next time step the
location of the character has changed accordingly (here $l_{n}$ denotes the
location at north of $l$):

$\displaystyle\forall l\in L_{nn}\setminus L_{n}\quad c_{l}^{t}\land
n^{t}\to(move:$ $\displaystyle\neg c_{l}^{t+1}\land c_{l_{n}}^{t+1}\land\neg
bs_{l_{n}}^{t}\land\neg bm_{l_{n}}^{t}\land\neg bl_{l_{n}}^{t})$

Finally, if the character moves north without having a wall two locations
ahead, apart from moving it can also _roll_ a ball north, _push_ a ball into a
stack of balls or _pop_ a ball from a stack of balls. For the sake of brevity
we only describe the _push_ north action (here $l_{nn}$ denotes the location
two steps ahead at north of $l$):

$\displaystyle\forall l\in L\setminus\\{L_{n}\cup L_{nn}\\}\quad
c_{l}^{t}\land n^{t}\to(move:\dots$ $\displaystyle{}\lor push:(\neg
c_{l}^{t+1}\land c_{l_{n}}^{t+1}\land((\neg bs_{l_{n}}^{t+1}\land
bs_{l_{nn}}^{t+1}\land{}$ $\displaystyle bs_{l_{n}}^{t}\land\neg
bm_{l_{n}}^{t}\land\neg bl_{l_{n}}^{t}\land{}\neg
bs_{l_{nn}}^{t}\land(bm_{l_{nn}}^{t}\lor bl_{l_{nn}}^{t}))$
$\displaystyle{}\lor(\neg bm_{l_{n}}^{t+1}\land bm_{l_{nn}}^{t+1}\land{}$
$\displaystyle\neg bs_{l_{n}}^{t}\land bm_{l_{n}}^{t}\land\neg
bl_{l_{n}}^{t}\land\neg bs_{l_{nn}}^{t}\land\neg bm_{l_{nn}}^{t}\land
bl_{l_{nn}}^{t})))$ $\displaystyle{}\lor roll:\dots\lor pop:\dots)$

  * •

Frame axioms impose that the state cannot change without a reason: snow cannot
be created, or if it disappears from a location it must be because a ball
occupies that location, etc.:

$\displaystyle\forall l\in L\quad(\neg s_{l}^{t}\to\neg s_{l}^{t+1})\land{}$
$\displaystyle(s_{l}^{t}\land\neg s_{l}^{t+1}\to bm_{l}^{t+1}\lor
bl_{l}^{t+1})\land\dots$

## Optimal plans with respect to ball movements

As previously done with PDDL, now we also modify the basic encoding by
considering only ball movement actions, asking in their preconditions that the
location from where the action is performed is reachable from the last
location of the character. Additionally, a cheating variant encoding is
included in the experiments for performance comparison.

### Standard Graph Reachability

We need to characterize reachability from the location of the character to the
appropriate location next to the ball to be moved. Gebser, Janhunen, and
Rintanen (2014) describe an SMT encoding for _source-target reachability_ in
graphs, where sources and targets are fixed beforehand. In this section we
provide a similar encoding to SAT which is able to deal with unknown a priori
sources and targets. We also incorporate the fact that accessible locations
may change due to ball movements.

The basic idea of the encoding is to prove the existence of a _reachability
path_ from the source to the target location. We stress that the way of
characterising reachability is intrinsically different in SAT or SMT than in
PDDL. In particular, we cannot mimic the idea followed in PDDL, where a
location $t$ is reachable from $s$ either if $s$ equals $t$ or there is some
(unoccupied) neighbour $s^{\prime}$ of $s$ such that $t$ is reachable from
$s^{\prime}$. The reason is that, whereas planning adheres to the closed-world
assumption (i.e., all unknown values are assumed to be false), this is not the
case for SAT and SMT. In other words, the translation of the inductive
reachability axioms of PDDL into SAT would be satisfied by any model where all
corresponding reachability variables are set to true, hence not encoding
reachability at all.

To deal with the open-world assumption of SAT (and SMT) when encoding
reachability in graphs, we essentially need to break cyclic relations, i.e.,
paths from source to target must be encoded as a transitive and antisymmetric
relation. As described in Gebser, Janhunen, and Rintanen (2014), acyclicity
can be easily modelled with SMT, since an ordering on locations can be imposed
by assigning a numeric value to each of them. The details of the encoding are
given below and the way we impose the ordering constraints in SAT is given in
next subsection.

We consider a directed graph representation of the maze where the nodes
correspond to the valid locations of the maze, and where there is a directed
edge from $l$ to $l^{\prime}$ if and only if $l$ and $l^{\prime}$ are adjacent
locations in the maze. In other words, we define a directed graph $G=(V,E)$
where $V=L$ and $E=\\{(l,l^{\prime}),(l^{\prime},l)\mid l,l^{\prime}\text{ are
adjacent locations in }L\\}$. The encoding goes as follows.

#### Variables

We introduce for all $l$ in $L$, $t$ in $0..T-1$, the following Boolean
variable:

  * •

$r_{l}^{t}$: node $l$ is reachable by the character at time $t$

For a node $l$ to be reachable, either it must be the source or it must be
reached through an edge ending in $l$ and coming from a reachable node.
Therefore, we need variables to state if an edge is building the reachability
path or not. For all edges $(l,l^{\prime})$ in $E$ and time $t$ we introduce
the following Boolean variable:

  * •

$e_{l\,l^{\prime}}^{t}$: edge $(l,l^{\prime})$ is in the reachability path at
time $t$

Moreover, cycles in those paths must be avoided. Therefore the SMT encoding
from Gebser, Janhunen, and Rintanen (2014) would also introduce, for each
location $l$ in $L$ and time $t$ the following integer variable:

  * •

$a_{l}^{t}\in 1..|L|$: value associated to node $l$ at time $t$; this is used
to break cycles: a topological ordering will be enforced among variables $a$
of the nodes in the reachability path to avoid cycles.

#### Constraints

The reachability constraints are the following, where
$next(l)=\\{l^{\prime}\mid(l,l^{\prime})\in E\\}$.

$\displaystyle\forall_{l\in L}\quad$ $\displaystyle
r_{l}^{t}\rightarrow\big{(}c_{l}^{t}\bigvee_{l^{\prime}:l\in
next(l^{\prime})}~{}e_{l^{\prime}\,l}^{t}\big{)}$
$\displaystyle\forall_{(l,l^{\prime})\in E}\quad$ $\displaystyle
e_{l\,l^{\prime}}^{t}\rightarrow(r_{l}^{t}\wedge
a_{l}^{t}<a_{l^{\prime}}^{t})$

Notice that this last ordering constraint is an SMT constraint because it is
using the difference logic atom $a^{t}_{l}<a^{t}_{l^{\prime}}$. We show in
next subsection how to replace it in SAT.

If there is a ball in a location then it is not reachable:

$\displaystyle\forall_{l\in L}\quad(bs_{l}^{t}\lor bm_{l}^{t}\lor
bl_{l}^{t})\rightarrow\lnot r_{l}^{t}$

###### Example 1

Consider the graph below. To ease notation, we refer to locations as numbers.

$1$$3$$2$$4$$5$$6$

At each time-step, some location $l$ will be asked to be reachable in some
action’s precondition. In other words, $r_{l}$ will be asked to be true for
some $l$ in $1..6$.

The encoding for reachability would be the following (we use $b_{l}$ as a
shorthand for $bs_{l}\lor bm_{l}\lor bl_{l}$ and omit time superindexes).

$\begin{array}[]{lll}r_{1}\to c_{1}&e_{13}\to r_{1}\wedge
a_{1}<a_{3}&b_{1}\to\neg{r_{1}}\\\ r_{2}\to c_{2}\vee e_{12}\vee
e_{32}&e_{12}\to r_{1}\wedge a_{1}<a_{2}&b_{2}\to\neg{r_{2}}\\\ r_{3}\to
c_{3}\vee e_{13}\vee e_{53}&e_{32}\to r_{3}\wedge
a_{3}<a_{2}&b_{3}\to\neg{r_{3}}\\\ r_{4}\to c_{4}\vee
e_{34}\vee{e_{64}}&e_{34}\to r_{3}\wedge a_{3}<a_{4}&b_{4}\to\neg{r_{4}}\\\
r_{5}\to c_{5}&e_{46}\to r_{4}\wedge a_{4}<a_{6}&b_{5}\to\neg{r_{5}}\\\
r_{6}\to c_{6}\vee e_{46}&e_{53}\to r_{5}\wedge
a_{5}<a_{3}&b_{6}\to\neg{r_{6}}\\\ &e_{64}\to r_{6}\wedge
a_{6}<a_{4}\end{array}$

### Translation of Ordering Constraints to SAT

The previous ordering constraints on the values of the numeric variables $a$
associated to locations can be easily translated to SAT. We propose not to
encode the numbers to binary form, but to encode the acyclicity relation
directly to SAT.

A _strict partial order_ is a relation $<$ that is irreflexive and transitive
(which implies antisymmetry as well). This is all we need to ensure
acyclicity. We can encode such a relation by adding the constraints

$\displaystyle\neg p_{ii}^{t}$ (irreflexivity) $\displaystyle p_{ij}^{t}\land
p_{jk}^{t}\to p_{ik}^{t}$ (transitivity)

where $p_{ij}^{t}$ are Boolean variables, for locations $i,j$ and time $t$.
Then, to get rid of SMT constraints and obtain a full SAT encoding, we can
simply replace the constraints

$\forall_{(l,l^{\prime})\in E}\quad
e_{l\,l^{\prime}}^{t}\rightarrow(r_{l}^{t}\wedge
a_{l}^{t}<a_{l^{\prime}}^{t})$

by

$\forall_{(l,l^{\prime})\in E}\quad
e_{l\,l^{\prime}}^{t}\rightarrow(r_{l}^{t}\wedge p_{ll^{\prime}}^{t})$

and no numeric variables $a_{l}^{t}$ are needed at all.

Notice that transitivity constraints $p_{ij}^{t}\land p_{jk}^{t}\to
p_{ik}^{t}$ are only needed for neighbours $j$ of $i$, since transitivity
follows by induction. This allows to reduce the number of such constraints
from cubic to quadratic.

### Grid Graph Reachability

Grids are a particular case of graphs, and a simpler encoding for reachability
in them is possible. The idea is to build a single path from the source to the
target as follows:

define a Boolean variable for each location denoting if it is included in the
path, and impose that the source and destination are in the path and exactly
one of their neighbours is in the path, and that any valid (non-obstacle)
location different from the source and the destination has either zero or two
neighbours in the path.

Not surprisingly, this simple approach is the best performing in our
experiments.

### Invariants

To help pruning the search space, some simple invariants can be considered.
Since balls cannot decrease their size, the following constraints have to be
fulfilled at any time: the number of large balls cannot exceed the number of
snowmen, and there must be at least as many small balls as snowmen. Imposing
these constraints has shown to be useful only in the harder instances.

Other invariants have been also considered, like forbidding to move balls to
locations where they would get stuck without the possibility of building a
snowmen, but this does not showed significant improvement.

## Empirical Evaluation

For our experiments we used the planners Fast Downward (Helmert 2006) v22.12
and SymK (Speck et al. 2019) v3.0, and the SAT solver Kissat (Biere et al.
2020) v3.0.0. These two planners were chosen due to their well-known
performance and support for the required features present in the PDDL models.
In particular, for Fast Downward we used three configurations: (i) the
integrated Stone Soup portfolio (Seipp and Röger 2018) (ii) LAMA (Richter and
Westphal 2010) and (iii) the blind heuristic. The Fast Downward Stone Soup
portfolio won the satisficing and cost-bounded tracks of the 2018
International Planning Competition (IPC). When using Fast Downward in
preliminary experiments, we observed that all configurations of Stone Soup
always resorted to the blind heuristic in our instances, and therefore we only
report results for LAMA and the blind heuristic. As SymK natively supports a
variety of PDDL features that are rarely supported by other planners, such as
conditional effects and derived predicates with axioms, we used its default
bidirectional search. Finally, the Kissat SAT solver and its variations were
the winners of various tracks of the 2021 and 2022 SAT competitions (Froleyks
et al. 2021).

To test the efficiency of the proposed approaches we considered all 30
instances of the main part of the game (the real world).222The game has a
second part (the dream world), governed by a different set of rules. We also
included a set of 9 _hidden_ instances, that are present in the game files but
not used in the actual game. Finally, we crafted 12 instances for some edge
cases. More concretely, we considered much bigger scenarios,333The game
scenarios contain at most 40 locations whilst the crafted ones have up to 99
locations. instances with snow everywhere and some really easy to solve
instances with two snowmen. In summary, we considered a total of 51 instances
ranging from 1 to 3 snowmen and from 14 to 99 valid locations. We ran the
experiments on a cluster of compute nodes equipped with Intel Xeon E-2234 CPU
@ 3.60GHz processors, where each execution was given a time-limit of 1 hour
and 16 GB of memory.

We considered three variations of the PDDL model: basic, cheating and
reachability. _Basic_ denotes the model with actions able to move the
character as well as balls. _Cheating_ is a model where we simply remove the
action where the character moves, and the character can do “unsound”
teleports. That is, it does not need to consider if there exists a path
between a source and a target position before acting from the target position.
Finally, _reachability_ denotes a model where derived predicates are used to
require reachability between subsequent action positions in the game. With
these three models we aim to evaluate the effect of different encodings of
reachability in PDDL.

| Solved | PAR-2
---|---|---
SAT | 27 | 196942
SAT-cheating | 43(19) | 62414
SAT-R-order | 43 | 62818
SAT-R-count | 43 | 60858
SAT-R-count+ | 43 | 60947
basic-blind | 27 | 177082
basic-LAMA | 14 | 275961
basic-SymK | 29 | 166980
cheating-blind | 25(16) | 188086
cheating-LAMA | 20(11) | 234919
cheating-SymK | 29(14) | 162480
reachability-blind | 26 | 183732
reachability-LAMA | 16 | 263718
reachability-SymK | 23 | 205981
Table 1:  Number of solved instances (within the 1 hour timeout), and PAR-2
scores (the score of a solver is defined as the sum of all runtimes for solved
instances + $2\times$timeout for unsolved instances) for each approach. For
cheating models the number of valid solutions is given in parenthesis. The
total number of instances considered is 51. Figure 3: Cummulative number of
solved instances for all the considered approaches. Time in seconds in
logarithmic scale.

Table 1 summarises all results.444Detailed results are given in the
supplementary material at https://github.com/udg-lai/KEPS2023. The different
SAT configurations are as follows: SAT refers to the basic SAT encoding
without reachability constraints, SAT-cheating to a variation of the SAT
encoding with “unsound” teleports, and the -R versions implement the
reachability encodings with the ordering and counting variations, where +
denotes the addition of implied constraints (invariants). The rest of the rows
show the combined performance of each planner and PDDL variant considered. For
example, basic-blind considers the basic model solved with Fast Downward using
the blind heuristic.

Column Solved shows the number of instances solved by each system. For the
cheating variants, it shows in parenthesis how many solutions are valid out of
all the ones solved. Column PAR-2 shows the sum of all runtimes for solved
instances + 2$\times$3600 seconds per timed-out instances. The cheating and
SAT-R variants are clearly the strongest approaches, solving 43 instances out
of 51, and showing that the reachability encoding is critical to be able to
obtain more (sound) solutions within the given time limit. The cheating
variant is not faster than the reachability based ones. Notice that it has
less variables but it also forbids less partial assignments, and therefore it
may spend more time exploring possible moves. This effect may be amplified by
the planning as satisfiability approach taken, as for finding a length-optimal
plan it previously proves unfeasibility of all shorter plans.

When comparing the cheating and basic PDDL models, the cheating variant only
results into notable gains with LAMA. When comparing the cheating model with
the reachability model, there is a notable difference in both LAMA and SymK.
Although not visible in the summarized results, amongst the PDDL approaches
with two or three snowmen, only SymK was able to solve two of them within the
given time limit. The rest of the planners were not able to preprocess any of
them due to the increased complexity of the domain.

Following the PAR-2 scores, three main groups can be observed: the SAT-R
variants as the best performers, SAT, blind search and SymK as the second
group, and finally LAMA being notably slower. The implied constraints
(invariants) did not seem to help in reducing the solving time, at least for
the instances considered. Figure 3 depicts the cummulative number of solved
instances over time. Note the time axis is in logarithmic scale. This figure
allows to identify more clearly how fast are the three aforementioned groups.

## Conclusions and Future Work

In this work we contributed new challenging planning benchmarks. In the
framework of puzzle-like games, we believe that SAT models like the ones we
have presented could be used to assist in the (semi)automatic design of new
scenarios, for instance, by ensuring that randomly generated or crafted
scenarios satisfy the desired constraints in terms of difficulty of the game
levels.

We have shown how relatively simple encodings to SAT outperform by far state-
of-the-art planners in the game _A Good Snowman is Hard to Build_ , especially
when considering reachability constraints.

Some works have introduced SAT and SMT solvers with support for detecting
acyclicity and reachability in graphs (Gebser, Janhunen, and Rintanen 2014;
Bayless et al. 2015). Therefore, an alternative approach could be to use a SAT
or SMT solver with built-in support for reachability, such as MonoSAT (Bayless
et al. 2015).

As future work, we plan to address the simultaneous (parallel) execution of
actions as an additional way of reducing the time span, in order to overcome
scalability issues and be able to solve the harder instances in a reasonable
amount of time. This will require non-trivial checking of possible
interferences between actions.

## Acknowledgements

Work partially supported by grant PID2021-122274OB-I00 funded by
MCIN/AEI/10.13039/501100011033 and by ERDF A way of making Europe.

## References

  * Bayless et al. (2015) Bayless, S.; Bayless, N.; Hoos, H. H.; and Hu, A. J. 2015.  SAT Modulo Monotonic Theories.  In Bonet, B.; and Koenig, S., eds., _Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence, January 25-30, 2015, Austin, Texas, USA_ , 3702–3709. AAAI Press.
  * Biere et al. (2020) Biere, A.; Fazekas, K.; Fleury, M.; and Heisinger, M. 2020.  CaDiCaL, Kissat, Paracooba, Plingeling and Treengeling Entering the SAT Competition 2020.  In Balyo, T.; Froleyks, N.; Heule, M.; Iser, M.; Järvisalo, M.; and Suda, M., eds., _Proc. of SAT Competition 2020 – Solver and Benchmark Descriptions_ , volume B-2020-1 of _Department of Computer Science Report Series B_ , 51–53. University of Helsinki.
  * Culberson (1997) Culberson, J. 1997.  Sokoban is PSPACE-complete.  Technical report, 97-02, Department of Computer Science, University of Alberta.
  * Davis and Hazelden (2015) Davis, B.; and Hazelden, A. 2015.  A Good Snowman Is Hard To Build.  https://agoodsnowman.com/.  [Online; accessed 23-March-2023].
  * Froleyks et al. (2021) Froleyks, N.; Heule, M.; Iser, M.; Järvisalo, M.; and Suda, M. 2021.  SAT Competition 2020.  _Artificial Intelligence_ , 301: 103572.
  * Gebser, Janhunen, and Rintanen (2014) Gebser, M.; Janhunen, T.; and Rintanen, J. 2014.  SAT Modulo Graphs: Acyclicity.  In _Logics in Artificial Intelligence - 14th European Conference, JELIA 2014, Funchal, Madeira, Portugal, September 24-26, 2014. Proceedings_ , 137–151.
  * Haslum et al. (2019) Haslum, P.; Lipovetzky, N.; Magazzeni, D.; and Muise, C. 2019.  _An Introduction to the Planning Domain Definition Language_.  Synthesis Lectures on Artificial Intelligence and Machine Learning. Morgan & Claypool Publishers.
  * He, Liu, and Yang (2017) He, W.; Liu, Z.; and Yang, C. 2017.  Snowman is PSPACE-complete.  _Theoretical Computer Science_ , 677: 31 – 40.
  * Helmert (2006) Helmert, M. 2006.  The Fast Downward Planning System.  _J. Artif. Intell. Res._ , 26: 191–246.
  * Ivankovic and Haslum (2015) Ivankovic, F.; and Haslum, P. 2015.  Optimal Planning with Axioms.  In Yang, Q.; and Wooldridge, M. J., eds., _Proceedings of the Twenty-Fourth International Joint Conference on Artificial Intelligence, IJCAI 2015, Buenos Aires, Argentina, July 25-31, 2015_ , 1580–1586. AAAI Press.
  * Kautz and Selman (1992) Kautz, H. A.; and Selman, B. 1992.  Planning as Satisfiability.  In _10th European Conference on Artificial Intelligence, ECAI 92, Vienna, Austria, August 3-7, 1992. Proceedings_ , 359–363.
  * Miura and Fukunaga (2017) Miura, S.; and Fukunaga, A. 2017.  Automatic Extraction of Axioms for Planning.  In _Proceedings of the Twenty-Seventh International Conference on Automated Planning and Scheduling, ICAPS 2017, Pittsburgh, Pennsylvania, USA, June 18-23, 2017._ , 218–227.
  * Rankooh and Rintanen (2021) Rankooh, M. F.; and Rintanen, J. 2021.  Propositional Encodings of Acyclicity and Reachability by using Vertex Elimination.  _CoRR_ , abs/2105.12908.
  * Richter and Westphal (2010) Richter, S.; and Westphal, M. 2010.  The LAMA Planner: Guiding Cost-Based Anytime Planning with Landmarks.  _J. Artif. Intell. Res._ , 39: 127–177.
  * Seipp and Röger (2018) Seipp, J.; and Röger, G. 2018.  Fast downward stone soup 2018.  _IPC2018–Classical Tracks_ , 72–74.
  * Silli (2010) Silli, E. 2010.  Mirror’s Edge - Level Design Challenges & Solutions.  Paper presented at GDC Europe 2010, Cologne, Germany.
  * Speck et al. (2019) Speck, D.; Geißer, F.; Mattmüller, R.; and Torralba, Á. 2019.  Symbolic Planning with Axioms.  In Lipovetzky, N.; Onaindia, E.; and Smith, D. E., eds., _Proceedings of the Twenty-Ninth International Conference on Automated Planning and Scheduling (ICAPS 2019)_ , 464–472. AAAI Press.