]https://polly.phys.msu.ru/en/labs/Tribelsky/

# The Poynting vector field singularities: Effects of symmetry and energy
conservation

Michael I. Tribelsky  E-mail: mitribel_at_gmail.com (replace “_at_” by @) [ M.
V. Lomonosov Moscow State University, Moscow, 119991, Russia  Center for
Photonics and 2D Materials, Moscow Institute of Physics and Technology,
Dolgoprudny 141700, Russia

###### Abstract

The study’s primary goal is to reveal the generic effects of the problem
symmetry, its violation, and energy conservation law on the Poynting vector
field singularities based on the study of resonant scattering of a linearly
polarized plane electromagnetic wave by an infinite right cylinder. The
polarization plane of the incident wave has an arbitrary orientation against
the cylinder axis ($z$-axis), and the wave vector is antiparallel to the
$x$-axis. The angle between the polarization plane and the $z$-axis plays the
role of a bifurcation parameter. We show that any deviation of the incident
wave from the pure TE or TM orientations makes the pattern of the Poynting
vector field lines three-dimensional. Meanwhile, the translational symmetry
along the $z$-axis remains. Accordingly, all singular points of the Poynting
vector field, but the ones lying on the $x$-axis, become “false”
singularities. They are singular in the projection of the field lines on the
plane $xy$. However, in the three-dimensional space, these points are regular
owing to the finiteness of the Poynting vector $z$-component. In contrast, the
singularities belonging to the $x$-axis remain the actual singular points at
any angle between the $z$-axis and the polarization plane since the
$z$-component of the Poynting vector for them vanishes owing to the problem
symmetry. We study the bifurcations related to the creation (annihilation) of
the false and actual singular points due to their splitting (merger) because
of the bifurcation parameter variations. In all inspected cases, a pitchfork
bifurcation occurs: the distance between the diverging (converging)
singularities as well as the corresponding roots of the characteristic
equation vary as the square root of a normalized deviation of the bifurcation
parameter from its critical value. We present a phenomenological theory,
explaining all observed peculiarities of the bifurcations. We unveil the
qualitative impact of dissipation on the symmetry of the bifurcation
scenarios. Remarkably, the bifurcation type is robust against the specific
choice of the bifurcation parameter. We verify this by comparing the
bifurcations caused by variations of the angle between the polarization plane
and the cylinder axis at a fixed radius of its cross section with the ones
caused by changes of the radius for the pure TM incidence of the scattering
wave.

##  I Introduction

The fundamental role of singularities in determining the topological
properties of any field is well-known. In the case of the electromagnetic
field, the study of singularities constitutes a discipline known as singular
optics. Remarkable achievements of singular optics are discussed in numerous
publications; see, e.g., [1, 2, 3, 4, 5, 6, 7, 8] and reference therein.
Singularities of the electromagnetic field may be related to those in the
fields $\mathbf{E}$ (electric), $\mathbf{H}$ (magnetic), and $\mathbf{S}$ (the
Poynting vector). Though all these fields are interconnected, each has its own
peculiarities, and the corresponding singularities should be discussed
individually.

Regarding singularities of the Poynting vector field, recently, their study
has been shifted to complicated configurations, such as, e.g., those in a non-
diffractive tractor beam [9], Bessel beam [10], etc. However, as a rule, the
developed approaches do not differ from those for an abstract vector field.
Meanwhile, the fields $\mathbf{E}$, $\mathbf{H}$, and $\mathbf{S}$ satisfy
Maxwell’s equations and the energy conservation law, respectively. Though
everybody knows it, the fact that these conditions impose certain constraints
on the features of the singular points has not been either noticed or given
sufficient attention. As a result, even in the most straightforward cases of
scattering of a linearly polarized plane electromagnetic wave, some important
generic properties of the singularities have not yet been revealed. In our
previous publications, we discussed some of them related to the effects of the
spatial dimension [11] and dissipative properties of the scatterer [12]. The
consideration was restricted by the highly symmetric light scattering problems
by a sphere and infinite right circular cylinder. In the latter case, only the
pure TE and TM polarizations have been considered.

In the present paper, we discuss deep connections between the properties of
the singularities and the problem symmetry. We develop a phenomenological
theory of the bifurcation phenomena. With its help, we unveil the effects of
symmetry violation and classify the types of bifurcations related to the
creation (annihilation) of singular points. We illustrate the general
discussion with a detailed study of light scattering by an infinite right
circular cylinder. We consider the normal incidence and an arbitrary
orientation of the incident wave polarization plane against the cylinder axis.
Variations of the angle between the vector $\mathbf{E}$ of the incident wave
and the cylinder axis make it possible to observe a smooth violation of the
symmetry from the highly symmetric pure TE or TM polarizations to the maximal
range of the violation at the angle equal to $45^{\circ}$. To verify the
universality of the obtained bifurcation features, we also inspect
bifurcations initiated by the change of the cylinder radius for the pure TM
polarization of the incident wave. Comparing the two cases (the fixed radius
and variable orientation of the polarization plane vs. pure TM polarization
and variable radius) confirms the results’ generic nature.

In all specific examples illustrating the general discussion, the calculations
are made with the help of the exact solutions to the corresponding problem.
The accuracy of the calculations is sufficient to resolve the finest details
of the discussed phenomena.

The paper has the following structure. In Sec. II, we give the problem
formulation and introduce our method to attack it. In Sec. III, we reveal the
symmetry of the fields $\mathbf{E}$, $\mathbf{H}$, and $\mathbf{S}$ at light
scattering by a cylinder, under the conditions specified above. In Sec. IV, we
discuss the general features of the Poynting vector field lines (streamlines)
and singularities. Sec. V is the key one: there, we develop a phenomenological
theory of the Poynting vector field singularities and inspect how the symmetry
and energy conservation law affect them. In Sec. VI, we present specific
examples of light scattering by a Germanium cylinder, illustrating the
developed phenomenology. In Conclusions, we briefly formulate the main results
of the study.

##  II Problem formulation and Method

Figure 1: Mutual orientations of the cylinder, the coordinate system, and the
vectors $\mathbf{k}$, $\mathbf{E}$, and $\mathbf{H}$ in the incident wave;
$\alpha$ is the angle between the $z$-axis and vector $\mathbf{E}$ equal to
the one between the $y$-axis and vector $\mathbf{H}$.

We consider an infinite right circular cylinder with a complex permittivity
$\varepsilon=\varepsilon^{\prime}+i\varepsilon^{\prime\prime}$ and
permeability $\mu=1$ (which is common in the optical frequency range) in a
vacuum irradiated by a linearly polarized plane monochromatic electromagnetic
wave. The temporal dependence of the electric and magnetic fields is selected
in the form $\exp(-i\omega t)$, so that in what follows $\mathbf{E}$ and
$\mathbf{H}$ stand for the corresponding complex amplitudes. We suppose
$\varepsilon^{\prime\prime}\geq 0$, i.e., we do not discuss a scatterer with a
population inversion, though the obtained results can be readily generalized
to this case by the formal sign change of $\varepsilon^{\prime\prime}$. For
the sake of simplicity, we consider only the normal incidence. We choose the
conventional orientation of the coordinate system, whose $z$-axis coincides
with the axis of the cylinder and the wave vector $\mathbf{k}$ of the incident
wave is antiparallel to the $x$-axis [13]. The mutual orientations of the
cylinder, the coordinate frame, and the vectors $\mathbf{k}$, $\mathbf{E}$,
and $\mathbf{H}$ in the incident wave are shown in Fig. 1. In this case, the
incident wave fields $\mathbf{E}$ and $\mathbf{H}$ have the following
components in the Cartesian coordinate system:
$\mathbf{E}=E_{0}(0,\sin\alpha,\cos\alpha)$;
$\mathbf{H}=H_{0}(0,\cos\alpha,-\sin\alpha)$;

We discuss the conventional real Poynting vector, whose definition (in
Gaussian units) is

$\mathbf{S}=\frac{c}{16\pi}(\mathbf{E}^{*}\\!\times\mathbf{H}+\mathbf{E}\times\mathbf{H}^{*}),$
(1)

where the asterisk stands for the complex conjugation.

However, sometimes the complex Poynting vector
$\hat{\mathbf{S}}=\frac{c}{16\pi}\mathbf{E}^{*}\\!\times\mathbf{H}$ is also
meaningful. Its imaginary part characterizes the alternating flow of the so-
called “stored energy” [14]. This imaginary part plays an important role in
some problems of light-matter interaction [15, 16, 17, 18, 19, 20, 21]. In
most cases discussed below, the generalization of the obtained results to the
complex $\hat{\mathbf{S}}$ is quite straightforward.

To make the results invariant against the specific choice of the system of
units (Gaussian, SI, etc.), we normalize $\mathbf{E}$, $\mathbf{H}$,
$\mathbf{S}$, and radius vector $\mathbf{r}$ on the corresponding values in
the incident wave and the radius of the cylinder cross section ($R$),
respectively. In what follows, we employ only these dimensionless quantities.
We keep for them the same notations since it cannot give rise to confusion.

All specific examples are based on the exact solution of the problem. This
solution exists at an arbitrary orientation of the vector $\mathbf{k}$ and
polarization plane against the cylinder axis. However, in the general case, it
is very cumbersome [13]. Therefore, for the problem in question, it is more
convenient, employing its linearity, to present the incident wave fields
$\mathbf{E}_{\rm inc}$ and $\mathbf{H}_{\rm inc}$ as a sum of the fields of
the pure TM- and TE-polarized waves, namely $\mathbf{E}_{\rm
inc}=\mathbf{E^{\text{\tiny TM}}_{\rm inc}}+\mathbf{E^{\text{\tiny TE}}_{\rm
inc}}$; $\mathbf{H}_{\rm inc}=\mathbf{H^{\text{\tiny TM}}_{\rm
inc}}+\mathbf{H^{\text{\tiny TE}}_{\rm inc}}$ and use the more simple exact
solutions for the pure TM and TE polarizations [13]. Here the vectors
$\mathbf{E^{\text{\tiny TM,\/TE}}_{\rm inc}}$ and $\mathbf{H^{\text{\tiny
TM,\/TE}}_{\rm inc}}$ have the following components in the Cartesian
coordinate system; see Fig. 1:

$\mathbf{E^{\text{\tiny TM}}_{\rm
inc}}=(0,0,\cos\alpha),\;\;\mathbf{H^{\text{\tiny TM}}_{\rm
inc}}=(0,\cos\alpha,0),$ (2)

and

$\mathbf{E^{\text{\tiny TE}}_{\rm
inc}}=(0,\sin\alpha,0),\;\;\mathbf{H^{\text{\tiny TE}}_{\rm
inc}}=(0,0,-\sin\alpha).$ (3)

Then, the fields $\mathbf{E}$ and $\mathbf{H}$, both outside and inside the
cylinder, are equal to the sum of those given by the solutions for the pure
TM- and TE-polarized incident waves, whose amplitudes are defined above. Next,
knowing the fields $\mathbf{E}$ and $\mathbf{H}$, we calculate $\mathbf{S}$
according to Eq. (1).

##  III Symmetry

For the pure TM or TE polarizations, the problem is symmetric against the
transformation $z\rightarrow-z$. Accordingly, the plane $z=0$ is invariant for
the vector Poynting field lines. An oblique orientation of the polarization
plane with respect to the $z$-axis breaks this symmetry, even if the vector
$\mathbf{k}$ still remains perpendicular to the $z$-axis. As a result, the
plane $z=0$ is not invariant anymore: generally speaking, for the points
belonging to this plane, the $z$-component of the Poynting vector ($S_{z}$)
does not vanish; see also the last paragraph of this section. However, the
problem still possesses the translational symmetry along the $z$-axis. It
means, that the fields $\mathbf{E}$, $\mathbf{H}$, and $\mathbf{S}$ depend on
the two spatial variables $x$ and $y$ solely.

To understand more subtle properties of the symmetry, note that owing to the
linearity of the problem, different components of fields $\mathbf{E}$ and
$\mathbf{H}$ are coupled only through the boundary conditions at the surface
of the cylinder. These conditions stipulate the continuity of the tangential
components of the fields. How does it affect the scattered fields outside the
cylinder and the ones inside it for the pure TM- or TE-polarized incident
wave? Consider, for definiteness, the pure TM case. Selecting the unit vector
$\mathbf{e}_{z}$ as one of the two independent tangential vectors to the
cylinder surface and taking the tangential unit vector belonging to the $xy$
plane as the other (actually, this is the vector $\mathbf{e}_{\varphi}$ of the
corresponding cylindrical coordinate system; see Fig. 2), we obtain that
$\mathbf{E^{\text{\tiny TM}}_{\rm inc}}$ is tangential to the cylinder surface
at any its point and does not have any projection onto $\mathbf{e}_{\varphi}$.
Therefore, the $x$ and $y$ components of the corresponding scattered field and
the field within the cylinder are not coupled to the incident wave field,
i.e., they are not excited. Thus, for all these fields, only $z$-component
does not vanish.

Figure 2: Unit vectors $\mathbf{e}_{r}$ and $\mathbf{e}_{\varphi}$ for an
arbitrary angle $\phi$.

In contrast, $\mathbf{H^{\text{\tiny TM}}_{\rm inc}}$ does not have any
projection onto $\mathbf{e}_{z}$ but at any $y\neq 0$ has non-zero projections
on the $x$ and $y$ axes. Accordingly, only these components of $\mathbf{H}$
are excited in this case. The conclusion is that for the pure TM-polarized
incident wave both the scattered electromagnetic field and the one inside the
cylinder occur the TM polarized too. Obviously, the analogous conclusion is
valid for the pure TE polarization of the incident wave.

Then, the solution of Maxwell’s equations for the given problem should be as
follows:

$\mathbf{E}=(E^{\text{\tiny TE}}_{x},E^{\text{\tiny TE}}_{y},E^{\text{\tiny
TM}}_{z}),\;\;\mathbf{H}=(H^{\text{\tiny TM}}_{x},H^{\text{\tiny
TM}}_{y},H^{\text{\tiny TE}}_{z}),$ (4)

where the superscripts TM and TE designate the solutions for the corresponding
pure-polarized incident waves, whose amplitudes are given by Eqs. (2), (3),
respectively. Importantly, all components of $\mathbf{E}$ and $\mathbf{H}$ in
Eq. (4) are functions of $x$ and $y$ and do not depend on $z$. Is there any
additional symmetry in this coordinate dependence?

According to the made problem formulation, the incident wave propagates
against the $x$-axis. It violates the $x\rightarrow-x$ symmetry. What is about
the symmetry against the $y\rightarrow-y$ transformation? To answer the
question, note that in the cylindrical coordinate system, the transformation
$y\rightarrow-y$ is equivalent to $\varphi\rightarrow-\varphi$; see Fig. 2.
This transformation does not affect the $y$-component of the
$\mathbf{e}_{\varphi}$ equal to $\cos\varphi$ and changes the sign of
$(e_{\varphi})_{x}=\sin\varphi$ to the opposite; see Fig. 2. Accordingly, for
the pure TM or TE polarizations both, the $y$-components of the fields
$\mathbf{E}$ and $\mathbf{H}$ (if any) should be invariant against this
transformation, while their $x$-components should change the signs. Applying
this rule to Eq. (4) we obtain

$\displaystyle(E^{\text{\tiny TE}}_{x}(x,y),E^{\text{\tiny
TE}}_{y}(x,y),E^{\text{\tiny TM}}_{z}(x,y))=$ $\displaystyle(-E^{\text{\tiny
TE}}_{x}(x,-y),E^{\text{\tiny TE}}_{y}(x,-y),E^{\text{\tiny TM}}_{z}(x,-y)),$
(5) $\displaystyle(H^{\text{\tiny TM}}_{x}(x,y),H^{\text{\tiny
TM}}_{y}(x,y),H^{\text{\tiny TE}}_{z}(x,y))=$ $\displaystyle(-H^{\text{\tiny
TM}}_{x}(x,-y),H^{\text{\tiny TM}}_{y}(x,-y),H^{\text{\tiny TE}}_{z}(x,-y)).$
(6)

Finally, according to Eq. (1), it gives rise to the following rules for the
Poynting vector components:

$\displaystyle(S_{x}(x,y),S_{y}(x,y),S_{z}(x,y))=$
$\displaystyle(S_{x}(x,-y),-S_{y}(x,-y),-S_{z}(x,-y)).$ (7)

Identities (5)–(7) also may be verified constructively based on the explicit
form of the exact solution for the problem in question.

Note that the obtained results are in full agreement with the conclusions
presented at the beginning of this section based on the problem symmetry.
Indeed, recall that, according to Eq. (1), vector $\mathbf{S}$ is
perpendicular to both $\mathbf{E}$ and $\mathbf{H}$. On the other hand, as it
follows from Eq. (4), for the pure TM polarization of the incident wave,
vector $\mathbf{E}$ of the solution of Maxwell’s equations is always parallel
to the $z$-axis. The same is true for the pure TE polarization for vector
$\mathbf{H}$. Then, in both these cases, the $z$-component of vector
$\mathbf{S}$ must equal zero identically. In contrast, for a mixed
polarization, $S_{z}$ may vanish only at certain specific values of $x$ and
$y$. This point will be important for the subsequent analysis.

##  IV Streamlines, true and false singularities

The pattern of the Poynting vector streamlines determines the topological
structure of this field. By definition, at any regular field point, the
Poynting vector is tangential to the streamline passing through this point.
Then, it is convenient to represent the streamlines in a parametric form,
considering them as “trajectories” of points moving in a three-dimensional
(3D) space: $\mathbf{r}=\mathbf{r}(t)$. In this case, the “velocity”,
$d\mathbf{r}/dt$ is tangential to the corresponding trajectory, i.e., parallel
to $\mathbf{S}$ at the same point of it. It means that the vector
$d\mathbf{r}/dt$ is proportional to $\mathbf{S}(\mathbf{r})$. We should stress
that $t$ here is just a parameter. It is not the actual time, and
$d\mathbf{r}/dt$ is not the actual velocity describing the propagation of the
electromagnetic energy along the streamlines: multiplication of $t$ by any
constant does not change the streamline shape. Therefore, the proper choice of
this constant always makes it possible to convert to unity the proportionality
coefficient between $d\mathbf{r}/dt$ and $\mathbf{S}(\mathbf{r})$. Then,
eventually, the equation determining the streamlines takes the form

$\frac{d\mathbf{r}}{dt}=\mathbf{S}(\mathbf{r}).$ (8)

Now we have to make an important comment. It has been stressed above that for
the problem in question, the Poynting vector field does not depend on $z$.
However, it is not the case for the streamlines. Indeed, generally speaking,
all three components of $\mathbf{S}$ do not equal zero. Then, according to Eq.
(8), for a given streamline $dz/dt=S_{z}\neq 0$. Since the coordinates $x$ and
$y$ are also $t$-dependent for the same reasons, it means that the streamlines
are essentially 3D. The independence $\mathbf{S}$ of $z$ means only that any
translation of a given 3D streamline along the $z$-axis transforms it to the
identical streamline, which is 3D too.

As mentioned above, the types and positions of singular points of Eq. (8)
determine the global topological structure of the field
$\mathbf{S}(\mathbf{r})$. Therefore these points are of special interest to
us. According to the general rules, the singularities are located at the
points of intersections of the three null isoclines: $S_{x,y,z}=0$. However,
for the problem under consideration, the field $\mathbf{S}$ depends only on
$x$ and $y$ and does not depend on $z$. In this case we have the three
equations $S_{x,y,z}(x,y)=0$ for two variables $x$ and $y$. This problem is
over-determined and does not have any solution. The only exception is the case
when the condition $S_{w}=0$ holds identically owing to the problem symmetry.
Here $w$ denotes any of the coordinates $x,y,z$.

As pointed out above, we do have such a case for the pure TM or TE
polarization when $S_{z}\equiv 0$ at any $x$ and $y$. What is about mixed
polarization? In this case, we have a singled out plane $y=0$. Indeed,
according to Eq. (7), $S_{y}(x,y)=-S_{y}(x,-y)$ and $S_{z}(x,y)=-S_{z}(x,-y)$.
Then, $S_{y}(x,0)=S_{z}(x,0)=0$, and for all points belonging to the plane
$y=0$ the equations $S_{x,y,z}=0$ transforms into the single equation
$S_{x}(x,0)=0$. Its roots, $x_{\rm sin}$ (if any), give the position of
singularities. The conclusion is that, for the problem in question, the only
true singularities may be observed either for the pure TM- or TE-polarized
incident wave or (at a mixed polarization) must belong to the plane $y=0$.

Now recall that the field $\mathbf{S}$ is $z$-independent. It means that, if
at a singular point $\mathbf{S}(x_{\rm sin},y_{\rm sin})=0$, this condition
holds at any value of $z$. Therefore, instead of singular points in all cases,
we have continuous straight singular lines parallel to the $z$-axis. Moreover,
the invariance of the problem against any translation along the $z$-axis makes
one of the roots of the characteristic equation for these singularities
identically equal to zero.

In addition to the discussed true singularities, there are the ones that we
call “false.” The point is that we can consider projections of the actual 3D
streamlines onto the $xy$ plane. These projections also play an essential role
in understanding the $\mathbf{S}$-field topological structure. The streamline
pattern in these projections is described by the two-dimensional (2D) version
of Eq. (8), where the vectors $\mathbf{r}$ and $\mathbf{S}$ have only $x$ and
$y$ components. The corresponding 2D pattern also may have singularities whose
coordinates are the solutions of the equations $S_{x,y}(x,y)=0$. However, we
must remember that, as shown above, if the $y$ coordinate of these points is
not equal to zero, $S_{z}$ does not vanish there. Thus, in 3D, these points
remain regular.

To illustrate this general discussion, we present an example of the true and
false singular points; see Fig. 3. The figure corresponds to the irradiation
of a circular cylinder with the permittivity $\varepsilon=17.775+0.024i$. It
is the permittivity of Germanium at $\lambda=1590$ nm [22];
$\alpha=45.403^{\circ}$, the size parameter $kR=1.62$. Note the sharp drop of
$|\mathbf{S}|$ in the vicinity of the $x$-axis owing to the vanishing of
$S_{y,z}$ at $y=0$.

##  V Phenomenology

Figure 3: The Poynting vector field within a Germanium circular cylinder;
$\alpha=45.403^{\circ}$. (a) 2D projection of the streamline pattern exhibits
three singularities indicated as black dots: two saddles ($\mathbf{r}_{1,2}$)
lying on the $x$-axis and a stable focus ($\mathbf{r}_{3}$) which does not
belong to this axis. While the saddles are the true singularities, the focus
is not. It is seen in (b) showing a 3D image of the separatrix emerging from
saddle $\mathbf{r}_{1}$ and winding around the vertical straight line
(dashed), whose projection onto the $xy$-plane corresponds to
$\mathbf{r}_{3}$. See text for details. Pay attention to the difference in the
scales of the axes.

Here we elucidate how the problem symmetry affects bifurcations in the
Poynting vector field. To make the consideration the most general, we develop
a phenomenological approach, taking into account only the restrictions imposed
by the symmetry, spatial dimension, and energy conservation law. We use the
explicit form of the solution only to illustrate  the general results.

The main idea is to expand the field $\mathbf{S}(x,y)$ in powers of the
deviations of $x$ and $y$ from a certain singled-out point (see below). Higher
symmetry of the problem corresponds to a larger number of constraints.
Therefore, in what follows, we consider the most complicated cases of a
creation (annihilation) of singularities belonging to the plane of the mirror
symmetry, $y=0$, when the number of the constraints is maximal. To describe
bifurcations occurring in less symmetric points, one has to remove the
constraints which do not arise in these cases.

As shown above, for the problem in question, the only true singularities are
parallel to the $z$-axis straight lines belonging to the plane $y=0$. The
other “singularities” are false since, for them, $S_{z}\neq 0$. However, if we
consider the projection of the streamlines onto the plane $z=0$ for the
corresponding 2D pattern, the false singularities transform into true 2D
singular points. We discuss the creation (annihilation) of singular points of
this 2D pattern occurring in the $x$-axis or close to it.

Let $\alpha$ be a control parameter whose smooth variations determine the
topological structure of the pattern. In the preceding sections, it was the
angle between the vector $\mathbf{E}$ of the incident wave and the $z$-axis.
However, it could be any other quantity, e.g., the size parameter, the real or
imaginary part of permittivity, etc. Its specific definition is not essential
for the developed approach. It is important only that there exists a
bifurcation value of this parameter $\alpha_{c}$ at which two singularities
merge and disappear or are born and separate.

We are interested in the values of $\alpha$ in the vicinity of $\alpha_{c}$.
In this case, it is convenient to introduce the normalized control parameter

$\eta=\pm\frac{\alpha-\alpha_{c}}{\alpha_{c}},$ (9)

where the sign in Eq. (9) is selected in such a way that the singularities
exist for positive values of $\eta$.

For the reduced 2D problem, the coordinates of singular points are the
solutions of the set of equations $S_{x}(x,y)=0$ and $S_{y}(x,y)=0$. In a
generic case, $S_{x,y}(x,y)$ are smooth, differentiable functions. Then, the
condition $S_{x}(x,y)=S_{x}(x,-y)$, see Eq. (7), gives rise to the equality
$\partial S_{x}/\partial y=0$ at $y=0$.

Next, as shown above, $S_{y}(x,0)=0$. Since the singularities are the
intersection points of two null isoclines, close to the bifurcation “moment”
the line $S_{x}(x,y)=0$ intersects the $x$-axis (the other null isocline) in
two points situated close to each other. That is to say, the functions
$S_{x}(x,0)$ vanishes at two close points. Then, in between, there is its
extremum, where $\partial S_{x}/\partial x=0$. Thus, we have obtained that,
close to the bifurcation, there is a point on the $x$-axis where $\partial
S_{x}/\partial x=\partial S_{x}/\partial y=0$. It is precisely the point where
the expansion of $\mathbf{S}(x,y)$ should be performed. In the subsequent
analysis, it is convenient to introduce a local coordinate frame with origin
coinciding with this point. The point itself may be either an extremum of the
function $S_{x}(x,y)$ or its saddle point. We discuss these two cases
separately.

###  V.1 $S_{x}$ extrema

The expansion of a smooth function about the point of its extremum in the
origin of the local coordinate frame gives rise to the following generic
expression for $S_{x}(x,y)$:

$S_{x}\approx S_{x0}+\frac{1}{2}\left[\left(\frac{\partial^{2}S_{x}}{\partial
x^{2}}\right)_{\\!\\!0}\\!\\!x^{2}+\left(\frac{\partial^{2}S_{x}}{\partial
y^{2}}\right)_{\\!\\!0}\\!\\!y^{2}\right]+\left(\frac{\partial^{2}S_{x}}{\partial
x\partial y}\right)_{\\!\\!0}\\!\\!xy,$

where the subscript “0” means that the corresponding quantities are calculated
at the origin of the local coordinate frame. The quantities
$\left({\partial^{2}S_{x}}/{\partial x^{2}}\right)_{\\!0}$ and
$\left({\partial^{2}S_{x}}/{\partial y^{2}}\right)_{\\!0}$ are negative, when
$S_{x0}$ is a point of a local maximum of $S_{x}(x,0)$, and positive, if it is
a local minimum. Note that since close to the point $x=0$, there are two
points where $S_{x}(x,0)=0$, in the case of a local maximum, $S_{x0}>0$,
while, in the case of a local minimum, it is negative.

Our case is not generic since $S_{x}(x,y)$ must satisfy the condition
$S_{x}(x,y)=S_{x}(x,-y)$. It stipulates the vanishing of
$\left({\partial^{2}S_{x}}/{\partial x\partial y}\right)_{\\!0}$. Eventually,

$S_{x}\approx\pm\left(a_{0}^{2}-a_{xx}^{2}x^{2}-a_{yy}^{2}y^{2}\right),$ (10)

where the signs plus and minus correspond to the vicinity of a local maximum
and minimum of $S_{x}(x,y)$, respectively. Here

$a_{0}^{2}\equiv|S_{x0}|,\;\\!a_{xx}^{2}\equiv\frac{1}{2}\left|\left(\frac{\partial^{2}S_{x}}{\partial
x^{2}}\right)_{\\!\\!0}\right|,\;\\!a_{yy}^{2}\equiv\frac{1}{2}\left|\left(\frac{\partial^{2}S_{x}}{\partial
y^{2}}\right)_{\\!\\!0}\right|,$ (11)

and all $a$’s are supposed to be non-negative. In this case, the equation of
the null isocline $S_{x}(x,y)=0$ describes an ellipse with the diameters
$d_{x}=2a_{0}/a_{xx}$ and $d_{y}=2a_{0}/a_{yy}$.

What is about the other null isocline, $S_{y}(x,y)=0$? Once again, we have to
begin with the generic expansion of a smooth function of $x$ and $y$, taking
into account all terms up to the second order in small $x$ and $y$ and
imposing the constraints stipulated by the problem symmetry. Note that,
according to this constraint, $S_{y}(x,0)=0$ at any $x$ (see above), i.e., one
of the branches of the null isocline $S_{y}(x,y)=0$ is the straight line $y=0$
(the $x$-axis). To satisfy this condition, all terms in the expansion of
$S_{y}(x,y)$ must have $y$ as a common factor. The most general type of this
expression is

$S_{y}\approx y\left[\left(\frac{\partial S_{y}}{\partial
y}\right)_{\\!\\!0}+\left(\frac{\partial^{2}S_{x}}{\partial x\partial
y}\right)_{\\!\\!0}\\!\\!x+\frac{1}{2}\left(\frac{\partial^{2}S_{x}}{\partial
y^{2}}\right)_{\\!\\!0}\\!\\!y\right]$ (12)

However, we have another constraint: $S_{y}(x,y)=-S_{y}(x,-y)$. It requires
the vanishing of $\left({\partial^{2}S_{x}}/{\partial y^{2}}\right)_{\\!0}$ in
the above expression. Finally, we obtain

$S_{y}\approx y(b_{y}+xb_{xy}),$ (13)

where the meaning of $b_{y,xy}$ follows from the comparison of Eq. (13) with
Eq. (12). Then, according to Eq. (13) the null isocline $S_{y}=0$ consists of
the two perpendicular straight lines: $y=0$ and $x=-b_{y}/b_{xy}$. In this
case the null isoclines $S_{x,y}=0$ have four intersection points with the
coordinates

$\displaystyle\mathbf{r}_{1,2}=\left(\pm\frac{a_{0}}{a_{xx}},0\right),$ (14)
$\displaystyle\mathbf{r}_{3,4}=\left(-\frac{b_{y}}{b_{xy}},\pm\sqrt{\frac{a_{0}^{2}b_{xy}^{2}-a_{xx}^{2}b_{y}^{2}}{a_{yy}^{2}b_{xy}^{2}}}\right),$
(15)

see Fig. 4

Figure 4: Null isocline $S_{x}=0$ (blue ellipse, solid line) and two branches
of null isocline $S_{y}=0$ (red dashed lines). Singular points
$\mathbf{r}_{1,2,3,4}$ are marked with crosses. When the control parameter
$\eta$ decreases, the ellipse contracts, while the dashed vertical line
remains at the same position. As a result, point $\mathbf{r}_{1}$ moves along
the $x$-axis toward the vertical branch of null isocline $S_{y}=0$, while
points $\mathbf{r}_{3,4}$ move toward the $x$-axis. Eventually, these three
points merge. For the smaller values of $\eta$ only $\mathbf{r}_{1,2}$ remain.
With a further decrease in $\eta$ and the corresponding contraction of the
ellipse, $\mathbf{r}_{1,2}$ approach each other and merge at $\eta=0$. See
text for details.

Importantly, all coefficients $a$ and $b$ in the above expressions depend on
the normalized control parameter, $\eta$. Moreover, at $\eta=0$, two points
$\mathbf{r}_{1,2}$ merge at the origin of the coordinate frame (we remind that
the origin is at the point of the local extremum, which lies in between
$\mathbf{r}_{1,2}$). At this “moment,” the origin inherits the singular nature
of $\mathbf{r}_{1,2}$ and becomes a singular point of the $\mathbf{S}$-field.
It means that $\mathbf{S}(0,0)=0$ at $\eta=0$. Then, expanding $S_{x0}$ in
powers of small $\eta$, we obtain that the expansion of
$a_{0}^{2}\equiv|S_{x0}|$ (not the one for $a_{0}(\eta)$!) begins from a
linear in $\eta$ term.

Regarding all other $a$’s and $b$’s, in a generic case, their expansions begin
with constants. Then, at small $\eta$ tending to zero, the diameters of the
ellipse in Fig. 4 contracts as $\sqrt{\eta}$, see Eq. (11). At the same time,
the vertical line of the null isocline $S_{y}=0$ remains practically fixed. If
so, at $\eta\rightarrow 0$ the points $\mathbf{r}_{3,4}$ move toward
$\mathbf{r}_{1}$ (or $\mathbf{r}_{2}$, depending on the sign of
$-b_{y}/b_{xy}$); eventually, merge with it and annihilate.

Thus, a decrease in $\eta$ results in two sequential bifurcations
corresponding to the fusion of the singularities. The first bifurcation is
$3\rightarrow 1$ (coalescence of $\mathbf{r}_{3,4}$ with either
$\mathbf{r}_{1}$ or $\mathbf{r}_{2}$), the second is $2\rightarrow 0$ (fusion
of $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$).

However, this scenario has an exception, namely the case $b_{y}=0$. In this
case, the points $\mathbf{r}_{3,4}$ are situated on the $y$-axis. Then,
instead of the two sequential bifurcations, we have a single bifurcation of
the type $4\rightarrow 0$.

To unveil the physical grounds for vanishing $b_{y}$, we have to recall the
energy conservation law. For the problem under consideration, its differential
form reads [23]:

${\rm div}\,\mathbf{S}=-q,$ (16)

where $q$ is the power dissipated in unit of volume. In our case, ${\rm
div}\,\mathbf{S}=\partial S_{x}/\partial x+\partial S_{y}/\partial y$ since
$\partial S_{z}/\partial z=0$. Then, according to Eqs. (10), (13) we obtain
that in the given approximation

$b_{y}+\left(b_{xy}\mp 2a_{xx}^{2}\right)x=-q.$ (17)

In non-dissipative media $q\equiv 0$, i.e., the l.h.s. of Eq. (17) must vanish
identically. It results in the conditions

$b_{y}=0,\;\;b_{xy}=\pm 2a_{xx}^{2}.$ (18)

Thus, based just on the general requirements following from the problem
symmetry and energy conservation, we have obtained that the sequence of the
bifurcations $4\rightarrow 2$ and $2\rightarrow 0$ may take place in
dissipative media only, while non-dissipative ones may exhibit the degenerate
bifurcation $4\rightarrow 0$ solely.

To understand the types of the fusing singular points we have to consider the
values of the Jacobian matrix of the system of equations $S_{x,y}(x,y)=0$ at
every singular point $\mathbf{r}_{n}$ and find its eigenvalues, i.e., the
roots $\kappa_{n}$ of the corresponding characteristic equation.

The degenerate non-dissipative case, when $q=0$ and constraint (18) is
applicable, is more straightforward for the analysis, so we begin with it. In
this case, applying the above procedure to Eqs. (10)–(15) gives rise to the
following expressions valid in the proximity of both local maxima and minima
of $S_{x}(x,y)$:

$\kappa^{\pm}_{1}=\kappa^{\pm}_{2}=\pm
2a_{0}a_{xx},\;\;\kappa^{\pm}_{3}=\kappa^{\pm}_{4}=\pm 2ia_{0}a_{xx},$ (19)

where the subscripts have the same meaning as those in $\mathbf{r}_{1,2,3,4}$
(see Fig. 4). Thus, points $\mathbf{r}_{1,2}$ are saddles, while
$\mathbf{r}_{3,4}$ are centers.

Note, $\mathbf{r}_{i}$ and $\kappa_{i}^{\pm}$ ($i=1,2,3,4$) all are scaled as
$\sqrt{\eta}$ owing to Eqs. (14), (15), (18), (19), and the mentioned above
$\eta$-dependence of coefficients $a$’s, $b$’s. In other words, here we have a
typical pitchfork bifurcation.

Next, for the discussed 2D streamline pattern the conditions Eq. (7) mean that
the pattern possesses a mirror symmetry with respect to the $x$-axis. Then,
the reflection against this axis transforms the stable and unstable
separatrixes of the saddles $\mathbf{r}_{1,2}$ into themselves. However, there
are just two types of lines remaining invariant for this transformation,
namely the $x$-axis itself and any line parallel to the $y$-axis. Thus, for
the saddles belonging to the $x$-axis, one of the separatrixes is the
$x$-axis, while the other is normal to it. Moreover, since a certain bunch of
streamlines, emerging from the vicinity of one saddle along the unstable
direction, approaches the proximity of the other along its stable direction,
it means that if for one saddle the separatrix parallel to the $y$-axis is
unstable (the corresponding $\kappa$ is positive), for the other it is stable
($\kappa<0$). The same conclusion is valid for the separatrixes coinciding
with the $x$-axis. This result may also be verified by directly calculating
the eigenvectors at these singular points.

What happens in dissipative media? To answer the question, we have to discuss
the explicit expression for $q$ in Eq. (16). Bearing in mind that, for the
problem under consideration, permeability $\mu\equiv 1$, in the chosen
dimensionless variables one has

$q=\varepsilon^{\prime\prime}kR|\mathbf{E}|^{2},$ (20)

where $k=\omega/c$ is the wavenumber of the incident wave and $c$ stands for
the speed of light in a vacuum [23, 12].

Now recall that Eqs. (10), (13) for $S_{x}$ and $S_{y}$ are approximate: only
the terms quadratic in both $x$ and $y$ are taken into account. Accordingly,
the accuracy of div$\,\mathbf{S}$, calculated based on these expressions, is
up to linear terms in $x$ and $y$. To maintain the same accuracy in the
expansion for $q(x,y)$, we must drop all terms higher than linear ones.

On the other hand, as it follows from Eq. (5), $E_{x}(x,0)=0$. Then, the
expansion of $E_{x}(x,y)$ in powers of small $x$ and $y$ begins with a linear
term, and hence $|E_{x}|^{2}$ is quadratic in small $y$ and should be dropped.

Next, the same Eq. (5) requires that $\partial E_{y}/\partial y=\partial
E_{z}/\partial y=0$ at $y=0$. It means that linear in $y$ terms do not enter
the expansions of $E_{y,z}(x,y)$. Regarding the quadratic terms, they should
be neglected according to the specified accuracy. Thus, eventually, we obtain

$q\approx\varepsilon^{\prime\prime}kR\left(|E_{y}(0,0)|^{2}+|E_{z}(0,0)|^{2}\right)=const>0.$
(21)

In this case, to keep the equality in Eq. (17) valid at any $x$ we must change
Eq. (18) to the following

$b_{y}\approx const<0,\;\;b_{xy}\approx\pm 2a_{xx}^{2}.$ (22)

Note that the condition $b_{xy}\mp 2a_{xx}^{2}\approx 0$ removes the problem
of the non-physical negativeness of $q$ at $x<b_{y}/(2a_{xx}^{2}\mp b_{xy})$,
which otherwise would arise according to Eq. (17).

In what follows, for definiteness, we consider the vicinity of a local maximum
of $S_{x}(x,y)$, choosing the plus sign on the r.h.s. of Eq. (10) (the
vicinity of a local minimum is considered in the same manner). In this case,
$b_{xy}\approx 2a_{xx}^{2}>0$ and the points $\mathbf{r}_{3,4}$ belong to the
right semiplane, as shown in Fig. 4.

Then, the same approach as that described above for the non-dissipative case
gives rise to the following results

$\displaystyle\kappa^{+}_{1,2}=\pm 2a_{0}a_{xx}+b_{y},\;\;\kappa^{-}_{1,2}=\mp
2a_{0}a_{xx},$ (23)
$\displaystyle\kappa^{\pm}_{3}=\kappa^{\pm}_{4}=\left(b_{y}\pm\sqrt{5b_{y}^{2}-(4a_{0}a_{xx})^{2}}\;\right)/2.$
(24)

In Eq. (23), the upper and lower signs correspond to $\mathbf{r}_{1}$ and
$\mathbf{r}_{2}$, respectively. Note, as follows from Eqs. (23)–(23),
$\kappa_{i}^{+}+\kappa_{i}^{-}=b_{y}$; $i=1,2,3,4.$

Consider the obtained eigenvalues. We remind that, for the case under
consideration, $b_{y}<0$, while all other $a$’s and $b$’s are nonnegative;
and, for the approximation made, the only coefficient depending on the control
parameter $\eta$ is $a_{0}\propto\sqrt{\eta}$. All other coefficients are
supposed to be constants. We will move from a certain initial positive value
of $\eta$ toward $\eta=0$ so that $a_{0}(\eta)$ decreases monotonically.

As for point $\mathbf{r}_{2}$, the change of $\eta$ does not affect its type.
Indeed, according to Eq. (23), $\kappa^{+}_{2}=-2a_{0}a_{xx}+b_{y}<0$;
$\kappa^{-}_{2}=2a_{0}a_{xx}>0$. These $\kappa$’s keep the same sign at any
value of $a_{0}(\eta)>0$, so the point $\mathbf{r}_{2}$ always is a saddle.

Regarding point $\mathbf{r}_{1}$, the eigenvalue
$\kappa^{-}_{1}=-2a_{0}a_{xx}$ remains a real negative quantity at any
$a_{0}$. In contrast, $\kappa^{+}_{1}=2a_{0}a_{xx}+b_{y}$ changes its sign as
$a_{0}$ decreases. As for eigenvalues $\kappa_{3,4}$, characterizing
singularities $\mathbf{r}_{3,4}$, they also change their type, when $a_{0}$
varies. Specifically, there are three critical values of $a_{0}(\eta)$, namely
(i) $a_{0}=-b_{y}\sqrt{5}/(4a_{xx})$, (ii) $a_{0}=-b_{y}/(2a_{xx})$ and (iii)
$a_{0}=0$. Let us discuss the streamline pattern in the proximity of each of
them.

(V.1.1) At $S_{x0}\equiv a_{0}^{2}>5b_{y}^{2}/(4a_{xx})^{2}$, the expression
under the square root in Eq. (24) is negative. Then, $\kappa_{3,4}$ are
complex quantities with the same negative real part, i.e., the points
$\mathbf{r}_{3,4}$ are stable foci. Regarding $\kappa^{+}_{1}$, for these
values of $a_{0}$ it is positive, so here the point $\mathbf{r}_{1}$ is a
saddle.

(V.1.2) At $S_{x0}\equiv a_{0}^{2}=5b_{y}^{2}/(4a_{xx})^{2}$, the square root
in Eq. (24) vanishes, $\kappa^{\pm}_{3,4}$ become real negative quantities,
and foci $\mathbf{r}_{3,4}$ transform into stable nodes. As for singularity
$\mathbf{r}_{1}$, at this value of $a_{0}$ we have
$\kappa^{+}_{1}=2a_{0}a_{xx}(1-2/\sqrt{5})>0$. Thus, $\mathbf{r}_{1}$ still is
a saddle.

(V.1.3) At $5b_{y}^{2}/(4a_{xx})^{2}>S_{x0}\equiv
a_{0}>b_{y}^{2}/(2a_{xx})^{2}$, we have $\kappa^{\pm}_{3,4}$ so points
$\mathbf{r}_{3,4}$ remain stable nodes. Regarding, $\mathbf{r}_{1}$, the
quantity $\kappa^{+}_{1}$ is positive in this domain, and the singularity type
for $\mathbf{r}_{1}$ does not change: it is a saddle.

(V.1.4) $S_{x0}\equiv a_{0}^{2}=b_{y}^{2}/(2a_{xx})^{2}$. At this value of
$a_{0}$ the three points, $\mathbf{r}_{1,3,4}$ merge; see Eqs. (15), (22), and
Fig. 4: $\kappa^{+}_{1,3,4}=0$, $\kappa^{-}_{1,3,4}=-2a_{0}a_{xx}<0$.

(V.1.5) $b_{y}^{2}/(2a_{xx})^{2}>S_{x0}\equiv a_{0}^{2}>0$. In this case
points $\mathbf{r}_{3,4}$ disappear; see Fig. 4. Point $\mathbf{r}_{2}$
remains a saddle, while point $\mathbf{r}_{1}$ inherits the properties of
absorbed singularities $\mathbf{r}_{3,4}$ and becomes a stable node.

(V.1.6) At $S_{x0}\equiv a_{0}^{2}=0$, points $\mathbf{r}_{1,2}$ merge. At
$\eta<0$, there are no singularities in the vicinity of the origin of the
local coordinate system.

It is convenient to introduce $\eta_{1,2}$ so that

$a_{0}(\eta_{1})=-\frac{b_{y}\sqrt{5}}{4a_{xx}},\;\;a_{0}(\eta_{2})=-\frac{b_{y}}{2a_{xx}},$
(25)

(we remind, $b_{y}<0$) and the corresponding local control parameters
$\tilde{\eta}_{n}=(\eta-\eta_{n})/\eta_{n}$, ($n=1,2$). Then, expanding
$a_{0}^{2}(\eta)$ in powers of $\eta$, we obtain the following full
bifurcation picture.

(V.1.7) Proximity of $\eta=\eta_{1}$. The number of the singularities does not
change. At $\tilde{\eta}_{1}>0$, points $\mathbf{r}_{3,4}$ are stable foci. At
$\tilde{\eta}_{1}\rightarrow 0$ from the side of positive $\tilde{\eta}_{1}$,
the imaginary parts of $\kappa^{\pm}_{3,4}$ tend to zero as
$\sqrt{\tilde{\eta}_{1}}$. At $\tilde{\eta}_{1}<0$, points $\mathbf{r}_{3,4}$
are stable nodes.

(V.1.8) Proximity of $\eta=\eta_{2}$. The stable nodes $\mathbf{r}_{3,4}$,
existing at $\tilde{\eta}_{2}>0$, merge with point $\mathbf{r}_{1}$ and
disappear. At vanishing $\tilde{\eta}_{2}$, $\kappa^{+}_{1,3,4}$ also vanish
as linear functions of $\tilde{\eta}_{2}$. These linear functions remain the
same for both $\tilde{\eta}_{2}>0$ and $\tilde{\eta}_{2}<0$. In contrast, the
distances between $\mathbf{r}_{1,3,4}$ are meaningful only at
$\tilde{\eta}_{2}>0$ and vanish as $\sqrt{\tilde{\eta}_{2}}$ at
$\tilde{\eta}_{2}\rightarrow 0$; see Eqs. (14), (15).

(V.1.9) Proximity of $\eta=0$. Stable node $\mathbf{r}_{1}$ merges with saddle
$\mathbf{r}_{2}$ at $\eta=0$ and both annihilate. Before the annihilation
$\kappa^{-}_{1,2}$ have the same modulus but opposite signs and vary as
$\sqrt{\eta}$; see Eq. (23). The distance between the two singularities also
obeys the same square root law; see Eq. (14). Note also that
$\kappa^{-}_{2}>0$, see Eq. (23) and the corresponding unstable separatrix of
saddle $\mathbf{r}_{2}$ is directed to the stable node $\mathbf{r}_{1}$ along
the $x$-axis.

###  V.2 $S_{x}$ saddles

In the vicinity of a saddle point of the surface $S_{x}(x,y)$, the
coefficients of $x^{2}$ and $y^{2}$ in Eq. (10) have opposite signs. Except
for this feature, Eq. (13) remains the same. However, this seemingly minor
change in the governing equations gives rise to drastic changes in the
bifurcation scenario. Let us discuss it in detail. For definiteness, we
suppose that instead of Eq. (10) now we have

$S_{x}\approx S_{x0}-a_{xx}^{2}x^{2}+a_{yy}^{2}y^{2},$ (26)

where, as before, $a_{xx,yy}>0$, however, now the sign of $S_{x0}$ is not
fixed. The case

$S_{x}\approx S_{x0}+a_{xx}^{2}x^{2}-a_{yy}^{2}y^{2},$

is quite analogous and will not be discussed.

For $S_{x}$ given by Eq. (26), Eq. (22) keeps the same form with $b_{y}\approx
const<0$ and $b_{xy}\approx 2a_{xx}>0$. The null isocline $S_{y}=0$ also does
not change. However, for the null isocline $S_{x}=0$, instead of an ellipse,
we have hyperbolas; see Fig. 5. It is convenient to discuss cases $S_{x0}>0$
and $S_{x0}<0$ separately.

$S_{x0}\equiv a_{0}^{2}>0$. In this case, depending on the $a_{0}$ value,
there are either only two singular points $\mathbf{r}_{1,2}$ on the $x$-axis,
or two more singularities $\mathbf{r}_{3,4}$ lying beyond the $x$-axis; see
Fig. 5a. Their coordinates are given by the same Eqs. (14), (15) after the
formal replacement $a_{yy}^{2}\rightarrow-a_{yy}^{2}$, cf. Eqs (10), (26). The
points $\mathbf{r}_{3,4}$ exist at $a_{xx}^{2}b_{y}^{2}\geq
a_{0}^{2}b_{xy}^{2}$.

$S_{x0}\equiv-a_{0}^{2}<0$. In this case, only singularities
$\mathbf{r}_{3,4}$ remain; see Fig. 5b. Their coordinates are given by Eq.
(15) after the replacement $a_{yy}^{2}\rightarrow-a_{yy}^{2}$,
$a_{0}^{2}\rightarrow-a_{0}^{2}$.

Once again, we begin the analysis with the degenerate non-dissipative case
when $b_{y}=0$ and the vertical branch of the null isocline $S_{y}=0$
coincides with the $y$-axis. In this case, at positive $S_{x0}$, points
$\mathbf{r}_{3,4}$ do not exist. Regarding points $\mathbf{r}_{1,2}$, the
corresponding $\kappa^{(1,2)}_{1,2}$ are given by the same Eq. (19). Thus,
these singular points are saddles. As $S_{x0}$ decreases, the points
$\mathbf{r}_{1,2}$ move toward the origin of the coordinate frame. It is
accompanied by a decrease in $\kappa^{\pm}_{1,2}$. At $S_{x0}=0$,
$\mathbf{r}_{1,2}$ merge at the origin of the coordinate system, and
$\kappa^{\pm}_{1,2}$ vanish.

A further decrease in $S_{x0}$ makes it negative, and two new singular points
$\mathbf{r}_{3,4}$ emerge from the origin and separate symmetrically moving
along the $y$-axis in opposite directions.

Figure 5: Null isoclines $S_{x}(x,y)=0$ (solid) and $S_{y}(x,y)=0$ (dashed) in
the vicinity of a saddle point of the surface $S_{x}(x,y)$, at different
values of $S_{x0}$. (a) $S_{x0}>0$; (b) $S_{x0}<0$. In the discussed
approximation, $S_{x0}$ changes do not affect the branches of the null
isocline $S_{y}(x,y)=0$ but do the ones of $S_{x}(x,y)=0$. Numbers 1,2,3 mark
three typical positions of the latter at three different values of $S_{x0}$.
An increase in the numbers corresponds to growth in $|S_{x0}|$. See text for
details.

Notably, in contrast to the cases discussed in the previous subsection, now
these points are saddles, whose $\kappa^{\pm}_{3,4}$ still equal to $\pm
2a_{0}a_{xx}$ (we remind that for $S_{x0}<0$ the quantity $a_{0}$ is defined
as $\sqrt{|S_{x0}|}$).

Collecting all together, we can say that, in this case, we always have two
saddles situated either on the $x$\- or $y$-axis. Their $\kappa^{\pm}=\pm
2a_{0}a_{xx}$, and are scaled as $\sqrt{|\eta|}$. The distance between the
singularities also obeys the same square root law.

This analysis shows that all equations derived in the previous subsections may
be applied to the present case after the replacement $a_{yy}^{2}\rightarrow-
a_{yy}^{2}$, at $S_{x0}\geq 0$, supplemented by the transformation
$a_{0}^{2}\rightarrow-a_{0}^{2}$, at $S_{x0}<0$. For this reason, we can
proceed to the inspection of the general bifurcation picture for the
dissipative case ($b_{y}\neq 0$) based on the results obtained in subsection
V.1 in the proximity of an extremum of $S_{x}(x,y)$. It is convenient to begin
the discussion from positive values of $S_{x0}$, gradually decreasing it.

(V.2.1) $S_{x0}\equiv a_{0}^{2}>b_{y}^{2}/(2a_{xx})^{2}$. In this case (curves
3 in Fig. 5a) only two singularities, $\mathbf{r}_{1,2}$ exist. They lie on
the $x$-axis symmetric against the point $x=0$, and their coordinates are
given by Eq. (14). As the control parameter varies, the distance between these
points changes as $a_{0}$, i.e., as $\sqrt{\eta}$, where $a_{0}=0$ at
$\eta=0$. The eigenvalues of these singularities are given by Eq. (23), which
is not affected by the replacement rules specified above since it does not
contain $a_{yy}$ and the sign of $a_{0}$ remains the same as that in Eq. (23),
because $S_{x0}>0$. Then, we readily obtain that, in the discussed domain of
the values of $S_{x0}$, $\kappa^{\pm}_{1,2}$ are purely real and, for a given
$i=1,2$, $\kappa^{\pm}_{i}$ have opposite signs, that is to say, both
singularities are saddles.

(V.2.2) $S_{x0}\equiv a_{0}^{2}=b_{y}^{2}/(2a_{xx})^{2}>0$. At this value of
$S_{x0}$ two more singular points, $\mathbf{r}_{3,4}$ begin to detach from
$\mathbf{r}_{1}$; see right curve 2 in Fig. 5a. The eigenvalues of the
resulting degenerate singularity are $\kappa^{+}=0$ and $\kappa^{-}=b_{y}<0$.

(V.2.3) $b_{y}^{2}/(2a_{xx})^{2}>S_{x0}\equiv a_{0}^{2}>0$. According to Eq.
(23), we have $\kappa^{+}_{1,2}<0$, $\kappa^{-}_{1}<0$, $\kappa^{-}_{2}>0$.
Thus, $\mathbf{r}_{1}$ is a stable node, while $\mathbf{r}_{2}$ is a saddle.
As for the new-born singularities, $\mathbf{r}_{3,4}$, their eigenvalues are
given by Eq. (24), which also can be applied without any modification. For the
case in question, $5b_{y}^{2}-16a_{0}^{2}a_{xx}^{2}>b_{y}^{2}>0$. Hence,
$\kappa^{\pm}_{3,4}$ are real quantities with the opposite signs, i.e., the
points $\mathbf{r}_{3,4}$ are saddles. Their coordinates are given by Eq.
(15), where $a_{0}^{2}$ should be replace by $-a_{0}^{2}$, supplemented by the
condition $b_{xy}^{2}=4a_{xx}^{2}$; see Eq. (22).

As $S_{x0}$ decreases, the points $\mathbf{r}_{3,4}$ symmetrically diverge
from the $x$-axis along the vertical branch of the null isocline $S_{y}=0$,
and the distance between them increases as
$\sqrt{b_{y}^{2}-(2a_{0}a_{xx})^{2}}/(a_{yy}a_{xx})$; see Eqs. (15), (22).
Introducing $\eta_{2}$, according to Eq. (25), and
$\tilde{\eta}\equiv(\eta-\eta_{2})/\eta_{2}$ we obtain that in proximity of
$\eta=\eta_{2}$ this distance varies as $\sqrt{\tilde{\eta}}$. In contrast,
$\kappa^{\pm}_{3,4}$ remain linear functions of $\tilde{\eta}$.

(V.2.4) $S_{x0}=0$. At this value of $S_{x0}$, the two branches of the null
isocline $S_{x}(x,y)=0$ degenerate into the two straight lines
$y=\pm(a_{xx}/a_{yy})x$ intersecting with the horizontal branch of the null
isocline $S_{y}(x,y)=0$ at the origin of the coordinate frame, and points
$\mathbf{r}_{1,2}$ merge. According to Eq. (23), the corresponding $\kappa$’s
are as follows: $\kappa^{+}_{1,2}=b_{y}<0$ and $\kappa^{-}_{1,2}=0$.

(V.2.5) $S_{x0}\equiv a_{0}^{2}<0$. Only singular points $\mathbf{r}_{3,4}$
remain; see Fig. 5b. The corresponding eigenvalues equal each other and are
$\kappa^{\pm}_{3,4}=\left(b_{y}\pm\sqrt{5b_{y}^{2}+16a_{0}^{2}a_{xx}^{2}}\,\right)/{2}$,
cf. Eq. (24). Both terms under the square root sign are positive, so these
singularities are saddles.

To conclude this section, note that in the vicinity of an extremum of $S_{x}$,
points $\mathbf{r}_{3,4}$ can be either foci or nodes; see cases (V.1.1),
(V.1.2)–(V.1.4). This is not the case in the vicinity of its saddle. It is
explained by the fact that the eigenvalues $\kappa_{3,4}$ could be complex at
$a_{0}^{2}>5b_{y}^{2}/(4a_{xx})^{2}$; see Eq. (24). On the other hand, these
singularities exist, provided $a_{0}^{2}<b_{y}^{2}/(2a_{xx})^{2}$; see Eq.
(15), supplemented with the replacements $a_{yy}^{2}\rightarrow-a_{yy}^{2}$
and $b_{xy}^{2}\rightarrow 4a_{xx}^{4}$; see Eq. (22). Then, since $1/4<5/16$,
$\kappa_{3,4}$ remain purely real in the entire range of the existence of
these singularities. That is to say, in the vicinity of a saddle point of
$S_{x}$, singular points $\mathbf{r}_{3,4}$ exist only as saddles.

##  VI Examples

Figure 6: Projection of streamlines onto the $xy$-plane. (a) Case (V.1.3);
$\alpha=45.404023^{\circ}$. (b) Case (V.1.5); $\alpha=45.40405^{\circ}$. See
text for details. Figure 7: Projection of the streamlines onto the $xy$-plane.
(a) Case (V.2.1); $kR=1.57311$. (b) Case (V.2.3); $kR=1.573091$. (c) Case
(V.2.5); $kR=1.573000$ See text for details.

Here we present the results of calculations, illustrating the above general
discussion, in the case of a Germanium right circular cylinder irradiated in a
vacuum by a monochromatic linearly polarized plane electromagnetic wave with
$\lambda=1590$ nm. According to Ref. [22], the permittivity of Germanium is
taken as follows: $\varepsilon=17.775+0.024i$. In the vicinity of a local
maximum of $S_{x}$, we consider an oblique orientation of the polarization
plane with respect to the cylinder axis. Thus, in this case, the control
parameter is the angle $\alpha$; see Fig. 1, while the radius of the cylinder
cross section is fixed so that the size parameter $kR=1.62$. In contrast, to
present the case of a saddle point of $S_{x}$, we consider the pure TM
orientation of the incident wave, varying the radius of the cylinder base so
that, in this case, the size parameter is the bifurcation one.

We have to stress that the smallness of the domain of the control parameters
where the bifurcation values corresponding to all discussed above bifurcations
are situated makes the experimental observation of the bifurcation presented
below a challenging problem. This fact should not disappoint the reader. The
main goal of the presented examples is to verify the correctness of the
developed phenomenological approach. As for the possible experimental
evidence, it may be obtained for other materials and/or range of parameters
where the domain is broad enough. The choice of the corresponding materials
and parameters is a separate problem beyond the scope of the present paper.

The vicinity of a maximum of $S_{x}$. The oblique orientation of the
polarization plane against the cylinder axis. The size parameter is fixed and
equal to 1.62.

  * •

$\alpha=45.403^{\circ}$. This value of $\alpha$ corresponds to case (V.1.1).
The projection of streamlines onto the $xy$-plane (the upper semiplane only)
are shown in Fig. 3(a).

  * •

$\alpha=45.404023^{\circ}$. This case corresponds to (V.1.3). The stable foci
$\mathbf{r}_{3,4}$ have transformed into stable nodes; see Fig. 6(a).

  * •

$\alpha=45.40405^{\circ}$. This is case (V.1.5): points $\mathbf{r}_{3,4}$
have merged and annihilated, transferring their stable node type to point
$\mathbf{r}_{1}$, while point $\mathbf{r}_{2}$ remains a saddle; see Fig.
6(b).

The vicinity of a saddle point of $S_{x}$. The pure TM orientation of the
incident wave. The size parameter varies.

  * •

$kR=1.57311$. This value of the size parameter corresponds to case (V.2.1):
only two singularities $\mathbf{r}_{1,2}$ exist. Both are saddles; see Fig.
7(a).

  * •

$kR=1.573091$. At this value of $kR$ we have case (V.2.3): four singular
points; $\mathbf{r}_{1}$ is a stable node, $\mathbf{r}_{2,3,4}$ are saddles;
see Fig. 7(b).

  * •

$kR=1.573000$. Case (V.2.5): only $\mathbf{r}_{3,4}$ remains. Both are
saddles; see Fig. 7(c).

Regarding the dependence of various quantities on the corresponding control
parameters discussed in this section, we have checked them in the proximity of
both a maximum of $S_{x}$ and its saddle point. The calculations show that all
predicted scalings (see (V.1.7)–(V.1.9), (V.2.1), (V.2.3), and the discussion
in the text) hold within the employed numerical accuracy.

##  VII Conclusions

The obtained results have revealed the deep connection between the topological
structure of the Poynting vector field, the problem symmetry, and energy
conservation law. Specifically, in the scattering of a monochromatic linearly
polarized plane wave by an infinite right circular cylinder, a deviation of
the orientation of the polarization plane from the pure TE or TM configuration
makes the Poynting vector streamlines essentially three-dimensional. In this
case, while the projection of the streamline pattern onto the plane of the
cylinder base may exhibit various types of singularities, in the three-
dimensional space, all of them, but the ones belonging to the $x$-axis (see
Fig. 1), are regular points.

We have also shown that the bifurcation in the two-dimensional projection of
streamline patterns may be explained and classified within the framework of
the developed phenomenological theory. The latter is based on expanding the
Poynting vector components in powers of the deviations of the coordinates from
the specific characteristic points, supplemented by the constraints imposed by
the symmetry and energy conservation law. All predictions made by this
phenomenological theory have been illustrated and verified by analyzing the
exact solution of the scattering problem for a Germanium cylinder. The
obtained results unravel essential features of the electromagnetic energy
circulation at light scattering by particles (including nanoparticles) and
shed new light on this old problem.

###### Acknowledgements.

The author is grateful to Boris Y. Rubinstein for valuable discussions and
help in the visualization of the field patterns. This research was funded by
the Russian Science Foundation (Project No. 21-12-00151) for the analytical
study and the Ministry of Science and Higher Education of the Russian
Federation (Agreement No. № 075-15-2022-1150) for computer calculations.

## References

  * Wang  _et al._ [2004] Z. Wang, B. Luk’yanchuk, M. Hong, Y. Lin, and T. Chong, Physical Review B 70, 035418 (2004).
  * Bashevoy  _et al._ [2005] M. Bashevoy, V. Fedotov, and N. Zheludev, Optics express 13, 8372 (2005).
  * Mokhov  _et al._ [2006] S. Mokhov, R. El-Ganainy, and D. N. Christodoulides, Optics Express 14, 3255 (2006).
  * Mokhun  _et al._ [2007] I. Mokhun, R. Khrobatin, A. Mokhun, and J. Viktorovskaya, Optica Applicata 37, 261 (2007).
  * Dennis  _et al._ [2009] M. R. Dennis, K. O’Holleran, and M. J. Padgett (Elsevier, 2009) pp. 293–363.
  * Novitsky and Barkovsky [2009] A. V. Novitsky and L. M. Barkovsky, Physical Review A 79, 033821 (2009).
  * Yue  _et al._ [2019] L. Yue, B. Yan, J. N. Monks, R. Dhama, C. Jiang, O. V. Minin, I. V. Minin, and Z. Wang, Scientific Reports 9, 20224 (2019).
  * Angelsky  _et al._ [2021] O. V. Angelsky, A. Y. Bekshaev, S. G. Hanson, I. I. Mokhun, M. V. Vasnetsov, and W. Wang, Frontiers in Physics 9, 10.3389/fphy.2021.651964 (2021).
  * Gao  _et al._ [2014] D. Gao, A. Novitsky, T. Zhang, F. C. Cheong, L. Gao, C. T. Lim, B. Luk’yanchuk, and C.-W. Qiu, Laser & Photonics Reviews 9, 75 (2014).
  * Klimov  _et al._ [2021] V. Klimov, R. Heydarian, and C. Simovski, Journal of the Optical Society of America B 38, C84 (2021).
  * Tribelsky and Rubinstein [2022a] M. I. Tribelsky and B. Y. Rubinstein, Nanomaterials 12, 1878 (2022a).
  * Tribelsky and Rubinstein [2022b] M. I. Tribelsky and B. Y. Rubinstein, Nanomaterials 12, 3164 (2022b).
  * Bohren and Huffman [1998] C. F. Bohren and D. R. Huffman,  _Absorption and Scattering of Light by Small Particles_ (WILEY-VCH Verlag, 1998).
  * Jackson [1998] J. D. Jackson,  _Classical Electrodynamics_ (Wiley, 1998) p. 832.
  * Bliokh  _et al._ [2014a] K. Y. Bliokh, Y. S. Kivshar, and F. Nori, Physical review letters 113, 033601 (2014a).
  * Bliokh  _et al._ [2014b] K. Y. Bliokh, A. Y. Bekshaev, and F. Nori, Nature Communications 5, 10.1038/ncomms4300 (2014b).
  * Bekshaev  _et al._ [2015] A. Y. Bekshaev, K. Y. Bliokh, and F. Nori, Physical Review X 5, 011039 (2015).
  * Xu and Nieto-Vesperinas [2019] X. Xu and M. Nieto-Vesperinas, Physical review letters 123, 233902 (2019).
  * Khonina  _et al._ [2021] S. N. Khonina, S. A. Degtyarev, A. V. Ustinov, and A. P. Porfirev, Optics Express 29, 18634 (2021).
  * Tang and Cohen [2010] Y. Tang and A. E. Cohen, Physical Review Letters 104, 163901 (2010).
  * Lininger  _et al._ [2022] A. Lininger, G. Palermo, A. Guglielmelli, G. Nicoletta, M. Goel, M. Hinczewski, and G. Strangi, Advanced Materials , 2107325 (2022).
  * [22] M. Polyanskiy, Refractive index database, http://refractiveindex.info/, accessed: Nov. 01, 2022\.
  * Landau  _et al._ [2013] L. D. Landau, J. Bell, M. Kearsley, L. Pitaevskii, E. Lifshitz, and J. Sykes,  _Electrodynamics of continuous media_ , Vol. 8 (Elsevier, 2013).