Gustavo E. Romero

*Gustavo E. Romero, Instituto Argentino de Radioastronomía (CONICET; CICPBA; UNLP), C.C. No. 5, 1894 Villa Elisa, Argentina.

# The content of astrophysical jets

Gustavo E. Romero  Instituto Argentino de Radioastronomía (IAR), CONICET;
CICPBA; UNLP, Buenos Aires, Argentina  Facultad de Ciencias Astronómicas y
Geofísicas, Universidad Nacional de La Plata, Buenos Aires, Argentina
<EMAIL_ADDRESS>
(3 November 2020; 4 November 2020; 5 November 2020)

###### Abstract

[Summary]Jets, collimated outflows of particles and fields, are observed in a
wide variety of astrophysical systems, including Active Galactic Nuclei of
various types, microquasars, gamma-ray bursts, and young stellar objects.
Despite intensive efforts along several decades of observational and
theoretical research, there are still many uncertainties and open questions
related to how jets are produced and what is their composition. In this review
I offer an outline of some current views on the content and basic properties
of astrophysical jets.

###### keywords:

Jets, black holes, magnetic fields, accretion

††articletype: Review00footnotetext: Abbreviations: AGN, Active Galactic
Nuclei; ADAF, Advection Dominated Accretion Flows; BH, Black Hole; DSA,
Diffusive Shock Acceleration; EHT, Event Horizon Telescope; FSRQ, Flat
Spectrum Radio Quasar; GRMHD, General relativistic magnetohydrodynamics; HD,
Hydrodynamics; ISCO, Innermost Stable Circular Orbit; MAD, Magnetically
Arrested Disk; MHD, Magnetohydrodynamics; RIAF, Radiatively Inefficient
Accretion Flow; SED, Spectral Energy Distribution; SMBH, Super Massive Black
Hole; YSO, Young Stellar Object.

##  1 Introduction

Astrophysical jets are collimated outflows of particles and fields observed on
various scales in several types of both Galactic and extragalactic sources.
These sources include Active Galactic Nuclei (AGNs), young stellar objects
(YSOs), X-ray binaries (microquasars), and gamma-ray bursts (GRBs). This is a
surprising variety of astrophysical systems involving many different scales
and numerous physical processes. The properties of jets also span astonishing
ranges. For instance, the jets of AGNs can remain well collimated along
distances of hundreds of kiloparsecs. In GRBs the jets seem to have initially
extreme velocities with Lorentz factors of several hundreds. The extent of
jets launched by YSOs can reach several parsecs. Regarding power, jets can go
from $10^{32}$ erg s-1 in YSOs up to $10^{52}$ erg s-1 in GRBs; i.e. they can
cover 20 or more orders of magnitude.

The first jet ever observed was discovered by Curtis 1 in the nearby galaxy
M87, a century ago. Despite intensive investigation across the entire
electromagnetic spectrum, our knowledge of astrophysical jets is still
fragmentary and incomplete. One important open issue is related to the
composition of these flows. Are they formed by an electron/proton plasma? Or
is it an electron/positron fluid? Or a mixture? Are some jets initially
Poynting-flux dominated? If so, how do they become matter loaded? What is the
origin of the relativistic particles that produce the non-thermal radiation we
observe from these jets? And what about the slow, non-relativistic jets we
observe in some protostars? The answers to these and other similar questions
are closely related to the problem of how the jets are launched and how they
interact with the environment. In this paper I shall discuss these problems,
offer some possible answers, and outline prospects for future research. I
shall focus mainly on the relativistic jets of AGNs, but I shall occasionally
comment on other types of jets. Because of limitations of space I decided to
focus on some general topics avoiding most details. This, I hope, will provide
the reader with a general framework to start with. In order to delve into the
subject, there are several useful monographs and books2, 3, 4, 5, 6.

##  2 Jet launching

The production of astrophysical jets seems to be always associated with the
presence of the same basic ingredients: a gravitational potential well,
angular momentum, accretion of matter, and magnetic fields. The potential well
is provided by a supermassive black hole in the case of AGNS, by a stellar
mass black hole or a neutron star in microquasars and GRBs, and by a protostar
in YSOs. The accreted matter with angular momentum comes from the central
region of a galaxy in AGNs and from a companion star in microquasars and short
gamma-ray bursts; in long GRBs what seems to be accreted is an imploding star;
finally in YSOs matter is accreted from the original molecular cloud where the
young star is being born. In all these processes an accretion disk is formed.
Some properties of these systems and the resulting jets are quantified in
Table 1.

Table 1: Properties of jets and central engines of various sources.

Property | GRBs | Microquasars | YSO | AGN
---|---|---|---|---
Mass of the accretor [M⊙] | $\sim 10$ | $\sim 10$ | $\sim 1-10$ | $\sim 10^{6}-10^{9}$
Size of the accretor [cm] | $\sim 10^{6}$ | $\sim 10^{6}$ | $\sim 10^{11}$ | $\sim 10^{11}-10^{15}$
Maximum magnetic field [G] | $\sim 10^{16}$ | $\sim 10^{7}$ | $\sim 10^{3}$ | $\sim 10^{3}-10^{5}$
Jet power [erg s-1] | $\sim 10^{50}-10^{52}$ | $\sim 10^{37}-10^{40}$ | $\sim 10^{32}-10^{36}$ | $\sim 10^{42}-10^{46}$
Lifetime [yr] | $\sim 10^{-6}$ | $\sim 10^{4}-10^{6}$ | $\sim 10^{4}-10^{5}$ | $\sim 10^{7}-10^{8}$
Jet size [pc] | $\sim 10$ | $\sim 10$ | $\sim 1$ | $\sim 10^{5}$

The structure and spectrum of the disks depend on the accretion rate $\dot{M}$
. Three basic kinds of accretion regimes can be distinguished depending on the
relation of the actual accretion rate to the critical rate defined as,

$\displaystyle M_{\rm crit}=\frac{L_{\rm Edd}}{c^{2}}\approx 1.4\times
10^{17}\frac{M}{M_{\odot}}\;\;\textrm{g\;s}^{-1},$ (1)

where

$\displaystyle L_{\rm Edd}=\frac{4\pi GMm_{p}c}{\sigma_{\rm T}}\approx
1.3\times 10^{38}\frac{M}{M_{\odot}}\;\;\textrm{erg\;s}^{-1},$ (2)

is the Eddington luminosity111This is the luminosity necessary to stop
spherical accretion by the radiation pressure generated by the accreted
thermal gas., $M$ is the mass of the accreting object, and the remaining
symbols have their usual meaning. At sub-critical rates $\dot{M}<\dot{M}_{\rm
crit}$ the disk is geometrically thin and optically thick; it can be described
using the standard model developed by Shakura and Sunyaev7. At very low
accretion rates, $\dot{M}<<\dot{M}_{\rm crit}$, optically-thin advection-
dominated regimes with very hot flows are possible8, 9. Finally, if the regime
is super-critical, $\dot{M}>>\dot{M}_{\rm crit}$, the disk becomes optically-
thick, geometrically slim or thick, and advection-dominated10, 11, 12.

The magnetic field in the disk is expected to be initially only toroidal. A
poloidal magnetic flux can be generated in situ through a turbulent dynamo
powered by magnetorotational instabilities13, 14. When the disk becomes a
magnetically arrested disk (MAD) in its inner region, powerful jets can be
launched with Lorentz factor $\Gamma_{\rm jet}\sim 10$ in the case of AGNs.

Powerful, magnetically dominated outflows can be driven either from the
accretion disk or the black hole magnetosphere222For YOSs clearly only the
first possibility is available.. Poloidal magnetic field lines inclined at
angles $i>30$ degrees to the disk rotation axis present an effective potential
that decreases along field lines. An outflow can be launched and driven away
by centrifugal forces in such a situation even for very low coronal
temperatures15. For angles $i<30$ degrees, the coronal plasma cannot be freely
driven by the centrifugal forces: instead, it must first overcome the
effective-potential barrier, which has a maximum far away from the disk.
Therefore, a strong thermal or radiative assistance is required to initiate
outflows with such geometry. This radiative assistance is naturally expected
in super-critical sources where the radiation pressure is dynamically
important. In any case, sub-relativistic, matter-loaded winds from the
accretion disks are expected in most cases. For YSOs this is the only
mechanism operating for launching jets.

The production of very strong and relativistic outflows requires a large
fraction of the gravitational energy of the accreting matter to be converted
to Poynting flux. This can only occur in the very central region of the
sources. The basic features of the magnetic model for the generation of
relativistic jets are the following16, 17, 18, 19, 20:

  * •

Relativistic jets are produced by rapidly rotating BHs fed by magnetized
accretion disks.

  * •

The magnetic field geometry in the innermost region is dominated by the
poloidal component.

  * •

The ultimate power source is the rotational energy of the black hole.

  * •

The energy is extracted via magnetic torque as Poynting flux.

  * •

Jet collimation is due to the external medium outside the MAD zone333Magnetic
self-collimation due to the toroidal component of the field is not effective
in relativistic plasma flows due to kink instabilities21, so confinement by
the pressure and inertia of an external medium seems to be quite essential22..

  * •

Jet acceleration is via conversion of the electromagnetic energy into bulk
kinetic energy.

A torsional Alfvén wave is generated by the rotational dragging of the
poloidal field lines in the ergosphere of the black hole. The field develops a
toroidal component producing a traveling perturbation. This wave transports
magnetic energy outward, causing the total energy of the plasma near the hole
to decrease to negative values. When this plasma with negative energy crosses
the horizon, the rotational energy of the black hole decreases. Through this
process, the energy of the spinning black hole is extracted by the magnetic
field. The result is that while plasma is carried into the hole only (and not
ejected), electromagnetic power is evacuated along the rotation axis as a
Poynting flux. The magnetic field is tied to the infalling plasma, not to the
horizon, as in the simplified split monopole model23. I emphasize that it is
the back-reaction of the magnetic field which accelerates the ergospheric
plasma to relativistic speeds counter to the hole’s rotation endowing the
plasma with negative energy. And is the accretion of this negative energy
plasma which spins down the hole. Since the energy output is due to the black
hole and not to the accreting matter (whose role is to produce and support the
poloidal field lines), the power of the jet can exceed the accretion power,
i.e.

$\displaystyle P_{\rm jet}=q\dot{M}c^{2}\;\;\;\textrm{with}\;q>1,$ (3)

as inferred for some blazars and GRBs24. A detailed calculation leads to the
following expression for the jet power23, 25:

$\displaystyle P\sim\frac{G^{2}}{8c^{3}}B^{2}M_{\rm BH}^{2}\;f(a),$ (4)

where $B$ is the field in the ergosphere, $M_{\rm BH}$ is the black hole mass,
and $a$ is the dimensionless spin parameter. The function $f(a)$ tends to
$a^{2}$ for low values of $a$ 19, yielding the well-known Blandford-Znajek
formula:

$P_{\rm BZ}\propto B^{2}M_{\rm BH}^{2}a^{2},$ (5)

or, in convenient units,

$P_{\rm BZ}\approx 10^{46}\,\left(\frac{B}{10^{4}\rm
G}\right)^{2}\left(\frac{M}{10^{9}M_{\odot}}\right)^{2}a^{2}\,\,\rm{erg\,s}^{-1}.$
(6)

This equation is accurate for $a\leq 0.5$. For higher spin, the formula under-
predicts the true jet power by a factor of $\approx 3$ and the 6th order
expansion of $f(a)$ is necessary25. In any case, Eq. 6 shows that a rapidly
rotating black hole in a MAD regime can produce a Poynting-flux-dominated jet
with an efficiency $q_{\rm BZ}=P_{\rm BZ}/\dot{M}c^{2}\gg 1$.

This Poynting outflow is created in a region where the magnetic field is $\sim
10^{4}$ G in AGNs (for other types of sources see Table 1). In such sources
the radius of the jet is of the order of the Schwarzschild radius of the black
hole, i.e. $\sim 10^{15}$ cm. This is much longer than the gyroradius of
relativistic particles:

$\displaystyle r_{\rm L}=\frac{\gamma_{p}m_{p}c^{2}}{eB}\simeq 3\times
10^{5}\gamma_{p}\left(\frac{B}{10^{3}\;\textrm{G}}\right)^{-1}\;\;\textrm{cm},$
(7)

where $\gamma_{p}$ is the Lorentz factor of protons. For electrons the
gyroradius is $\sim 10^{3}$ times smaller than for protons. Clearly, charged
particles cannot penetrate the outflow from outside: the jet is magnetically
shielded. The region around the black hole where the jet is launched is
illustrated in Figure 1 for a radiatively inefficient accretion flow.

Figure 1: Sketch of the region around a black hole with a radiatively
inefficient accretion flow. SSD means the Shakura-Sunyaev thin disk. Not drawn
to real scale. Adapted from Gutiérrez et al. 26.

Despite the fact that the evacuation funnel around the rotation axis of the
black hole is expected to be free of plasma, very-high-resolution images of
the central source of the nearby jet in M87 obtained with the Event Horizon
Telescope (EHT)27 at 230 GH and Very Long Baseline Array data at 43 and 86
GHz28 reveal the existence of radiation associated with the jet at distances
down to $\sim 10$ gravitational radii from the supermassive black hole29. This
radiation is produced by relativistic electrons or positrons injected and
accelerated somehow very close to the BH. The jet also accelerates from
$\approx 0.3c$ at 0.5 mas from the black hole to $\approx 2.7c$ at $20$ mas (1
mas $\approx 250\;r_{\rm g}$)30, showing an efficient conversion of magnetic
energy into bulk motion and internal energy of the emitting gas. What is the
origin of this matter so close to the black hole? A possible answer is that
charged particles are locally created by gamma-ray photon annihilation and
neutron decays. These neutral particles must be produced in the surrounding
disk20, 31, 32, 33.

##  3 Matter loading from the hot disk

The physics of the accretion of matter onto a black hole is complex since the
accretion process develops in many different ways depending on the ambient
conditions. To characterize the different regimes, it is useful to parametrize
the accretion rate as a function of the Eddington rate:

$\dot{m}=\dot{M}/\dot{M}_{\rm Edd},$ (8)

where $\dot{M}_{\rm Edd}=10L_{\rm Edd}/c^{2}$ (assuming an efficiency of
$10\%$).

At very high accretion rates $\dot{m}>>1$ the accretion is super-Eddignton and
photon-trapping in the disk becomes important34, 35: the photon diffusion
timescale exceeds the accretion timescale, so photons are advected towards the
compact object without being able to reach the surface of the disk and escape.
Advection of photons attenuates the disk luminosity, which remains of $\sim
L_{\rm Edd}$. Powerful winds are expelled from the disk by the intense
radiation pressure: practically all the accreted matter is ejected36, 37. This
regime occurs in GRBs, AGNs undergoing tidal disruption events, and some
microquasars such as SS 433.

At moderate to high accretion rates, $\dot{m}\gtrsim 0.05$, the usual picture
of the accretion flow onto supermassive black holes consists of a
geometrically thin ($h/r\lesssim 0.1$), optically thick, cold ($T\ll T_{\rm
vir}$) disk that cools efficiently, and a geometrically thick ($h/r\gtrsim
0.5$), optically thin, hot ($T\approx T_{\rm vir}$) corona above and below the
disk. The corona is thought to be responsible for the hard X-ray tail detected
from active galactic nuclei and microquasars38, 39, 40, 41. This emission is
produced by hot electrons via inverse Compton up-scattering of UV/X-ray
photons coming from the thin disk. A typical example is the well-known
microquasar Cygnus X-1.

At low accretion rates, $\dot{m}<0.05$, the thin disk does not extend down to
the innermost stable circular orbit (ISCO) and it is truncated at longer
distances from the hole. The plasma in the inner region is thought to be in
the form of a Radiatively Inefficient Accretion Flow (RIAF) 8, 9. A RIAF is
similar to a corona; namely, it is a hot, inflated, optically thin plasma. The
hot electrons there Compton up-scatter not only photons from the thin disk but
also low-energy (IR) photons produced by synchrotron radiation by the same
population of hot electrons ($T\sim 10^{9}$ K). The physical properties of the
optically thin, hot component in accretion flows, either a corona or a RIAF,
are less understood than those of a cold thin disk. In particular, this hot
plasma lies in the so-called collisionless regime, and thus a fraction of the
particle population can be far out of thermal equilibrium. The presence of a
non-thermal component in the coronal region leads to the emission of high-
energy photons and other particles42, 43, 44, 45.

Hot accretion flows are more efficient than thin disks in launching and
collimating relativistic jets since the conditions discussed in Sect. 2 are
more easily realized close to the black hole in such regime. The jet is then
initiated as a Poynting flux emerging from the BH ergosphere. This flow can be
loaded through neutral particles coming from the surrounding hot flow. More
specifically, the particles and involved processes are20:

  * •

High-energy photons. They are generated in the hot flow by several processes
(most notably proton-synchrotron and proton-photon interactions) and
annihilate in the evacuation funnel of the ergospheric jet against background
photons. The annihilation process, $\gamma+\gamma\rightarrow e^{-}+e^{+}$,
injects pairs at the base of the jet. Direct photon annihilation occurs in two
energy bands in the context of hot accretion flows: two photons with similar
energy ($\sim$ MeV) collide producing the low energy peak of the spectral
energy distribution, and a high-energy photon annihilates with a low-energy
photon (IR, from disk synchrotron radiation). This latter process generates
the second peak in the spectrum at higher energies. As an example I show in
Figure 2 the pair injection spectrum calculated by Romero and Gutiérrez20 for
a BH similar to that in M87. About $10^{3}$ pairs per second per cubic
centimeter are injected above the black hole in this case.

  * •

Protons injected by neutron decays. The neutrons are generated in the hot flow
by $pp$ and $p\gamma$ collisions. Energy distributions of neutron injection
for both processes in the same conditions as before are shown in Fig. 3. The
neutrons can decay inside the jet funnel yielding protons and electrons:
$n\rightarrow p+e^{-}+\bar{\nu_{e}}$. The protons, in turn, can produce
additional pairs by Bethe-Heitler mechanism: $p+\gamma\rightarrow
p+e^{-}+e^{+}$. More protons can be created by conversion of flying neutrons:
$n+\gamma\rightarrow p+\pi^{-}$.

Figure 2: Left panel: Number of pairs created per unit time and volume by
photon annihilation in a hot accretion flow as a function of the radius
(distance to the black hole) in units of the Schwarzschild radius. Parameters
similar to those of M87 (from Romero and Gutiérrez20). Right panel: Spectral
energy distribution of the same electron/positron pairs in the innermost
region of the accretion flow.

Figure 3: Energy distribution of neutron injection in the same accretion flow
of Fig. 1.

Altogether, most protons injected at the base of the jet are the result of
neutron decays, although direct neutron conversion into protons might
represent a significant fraction close to the black hole. The ratio of the
total number of pairs to the total number of protons is $R_{e^{\pm}/p}\sim
10^{5}$ and the associated energy ratio is $R_{u_{e}^{\pm}/u_{p}}\sim 10^{2}$.
Most of the jet mass load by neutral particles in the base of the jet seems to
be, then, under the form of pairs. The exact ratio, of course, will change
from one source to another depending on the acceleration efficiency and other
parameters, but the number of pairs will always exceed the number of protons.

##  4 Entrainment and obstacles

The production rate of pairs in the funnel is $\sim 10^{48}$ s-1 for an
accreting BH similar to the one in M87. This injection implies a power of
$L_{e^{\pm}}\sim 10^{42}$ erg s-1. If there is not further mass loading, the
bulk Lorentz factor of the jet should be $\Gamma\sim P_{\rm BZ}/L_{e^{\pm}}$.
If the spin of the black hole is $a\sim 1$, $B\sim 10^{4}$ G, and $M_{\rm
BH}\sim 10^{9}$ $M_{\odot}$, we have $P_{\rm BZ}\sim 10^{46}$ erg s-1 from Eq.
(6). Then, $\Gamma\sim 10^{4}$. This value is far from the typical value
$\Gamma\sim 10$ inferred from superluminal motions in AGNs. This discrepancy
is pointing out to the fact that somehow more baryonic matter should be
injected close to the base of the jet by some additional mechanism.

One possibility is the accretion of non-dipolar large scale magnetic fields
onto the central object. Numerical simulations show that such accretion
creates instabilities that mix disk matter with the outflow charging and
slowing down the jet46. The non-dipolar fields are the results of turbulence
in the accreting flow. A same source might perhaps switch from dipolar to
quadrupolar field accretion, so the jet would decelerate. Then, when the
conditions are propitious again, the jet might accelerate and collide with the
older outflow forming shocks. These shocks might be responsible for some of
the moving knots observed in several sources.

Baryons can be also supplied by interactions of the jet with different types
of obstacles found in its way out. Such obstacles include clouds of the broad
line region47, 48, clumps in a wind49, or stars50, 51, 52, 53. In microquasars
clumps and interactions with the stellar wind provided by the donor star are
likely the main source of protons for the jet54, 55, 56, 57.

In the case of the jets of microquasars, where the accretion rate is higher
than in most AGNs, baryon load can be important even close to the base of the
jet. As a consequence, some jets in these objects become sub-relativistic and
dark transporting energy to large distances without being detected unless they
find some obstacle that allows part of the kinetic energy to be reconverted
into internal energy and then radiated away58, 59, 60.

On the other extreme of the power range, the duration distribution of long
gamma-ray bursts reveals a plateau at durations shorter than $\sim 20$ s (in
the observer frame) and a power-law decline at longer durations. Such a
plateau arises naturally in the context of the collapsar model because the
engine has to operate long enough to push the jet out of the stellar envelope,
so the observed duration of the burst is the difference between the engine’s
operation time and the jet breakout time. The breakout time for the jet
inferred from the duration distribution is similar to the breakout time
obtained with analytic estimates and numerical simulations (both 2D and 3D) of
a hydrodynamic jet ($\sim 10$ s for typical parameters). This suggests that
jets from collapsars are not Poynting flux dominated61. The mass in these jets
might come initially from neutrino annihilation and, after the breakout, from
interactions with the external medium.

##  5 Signatures of matter in the jets

The power and actual composition of jets can be investigated through their
effects on the external medium and from modeling their broadband spectral
energy distribution (SED). Studies of the radio and X-ray lobes of radio
galaxies suggest that the proton content in the jets tends to be higher in
Fanaroff–Riley I (FRI) radio galaxies than in FR IIs62. Jet-environment
interactions also can inform about the jet power and content in the case of
microquasars and YSOs. Microquasars such as Cygnus X-1 or SS 433 produce heavy
hydrodynamical jets with a significant contribution from cold protons to the
total kinetic energy flux.60, 63.

The detection of gamma-ray emission that dominates the spectral energy
distribution in many blazars can be also used to probe the nature of the jets.
This radiation is usually thought to be the result of inverse Compton up
scattering of either synchrotron or thermal photons. If the jet is Poynting
flux dominated, the synchrotron component should be stronger than observed.
This suggests that the outflow is an electron-proton jet, or at most a pair
plasma plus an electron-proton fluid64.

An upper limit to the pair content of the initial jet can be imposed
considering that it should move through the external UV photon field of the
inner accretion disk. A large content of pairs in the relativistic flow would
Compton up scatter this radiation field producing a soft X-ray bump that is
not observed. The conclusion is that pure electron-positron jet models can be
excluded as over predicting soft X-ray radiation. Similarly, pure electron-
proton jets can also be ruled out because they under predict the soft X-ray
luminosity of most blazars65. What remains is a scenario where the flow starts
being Poynting flux dominated but soon is loaded with pairs and electron-
proton cold plasma65, 66, 20. The scales on which the magnetization of the jet
drops below unity depend on the initial matter load and the geometry of the
jet, but typically should not exceed $\sim 10^{3}-10^{4}$ Schwarzschild
radii20, at least for most AGNs.

Further insights on jet composition can be obtained from polarization
observations67. The observed radio polarization is dependent on plasma
composition due to Faraday conversion and rotation effects. In particular, the
contribution of the electrons is degenerate with that of the positrons and
their respective number densities appear additively in the radiative transfer
functions relating to emission and linear polarization. The degeneracy is
broken by the consideration of circular polarization, which cancels out for a
symmetric $e^{-}e^{+}$ plasma, but does not of $ep$ plasma. Observations of
circular polarization have been used to infer the presence of pairs in jets,
but the interpretation of the results depends on various assumptions if the
source is not adequately resolved68. Hence, more accurate determinations
should wait to detailed polarimetric imaging with the Event Horizon Telescope
of M8727. Among the observational signatures predicted for this source are a
slight bilateral asymmetry in the flux along the innermost part of the jet and
one independent linear polarization component, and bilateral anti-symmetry in
the other component. The clearest observable feature is an enhancement of
circular polarization with increasing mangetization and ion content69.

On large scales, by the other hand, the presence of hydrodynamical plasma is
revealed by the development of substantial back-flowing cocoons. Instead,
magnetic dominance of the jet creates a “nose cone” – a shaped head, which is
rarely observed, and might be produced also by some combinations of heavy jets
and ambient medium70, 71. Back-flowing cocoons have been recently observed in
microquasars as well, in accordance with the expected hydrodynamic nature of
their jets72, 73.

The cumulative evidence tends to support a picture of relativistic jets where
the outflow starts as a Poynting flux but soon becomes matter dominated. The
jet accelerates along distances of $10^{3}-10^{4}$ $R_{\rm Schw}$. The matter
is an electron-proton plasma plus a pair content, which can vary according to
the manifold conditions of different sources. The inertia is determined by the
protons, but radiative properties are mostly due to the leptons. The
relativistic jet is likely accompanied by a magneto-centrifugal wind driven by
the disk. Mixing and entrainment of matter from the wind is likely.

##  6 Conclusions

Despite several decades of intensive research on astrophysical jets, many
basic issues remain unclear. With the advent of gamma-ray instruments such as
Fermi, high-resolution X-ray telescopes as Chandra, and sub-millimeter
interferometry, a consistent picture is finally emerging. Both disks and
rotating black holes can launch jets if poloidal magnetic fields are accreted,
as predicted in the 1970s and early 1980s. Disks produce magnetohydrodynamic
outflows, that are mostly subrelativistic. Magneto-centrifugal outflows formed
by $ep$ gas are present in YSOs, microquasars, and some AGNs. Faster,
relativistic jets are launched by the ergosphere of Kerr black holes as a
Poynting flux but soon become matter dominated. The total power of some of
these jets might even exceed the accretion power, clearly indicating that
energy is being extracted from the rotation of the black hole. The kinetic
power beyond a fraction of parsec from the black hole is dominated by the
inertia of protons, but in general the radiative output is due to leptons
through synchrotron radiation and inverse Compton interactions. Hadronic
radiative processes, however, can be important for the baryon load of the
initial jet and during interactions with various kinds of obstacles. The
acceleration of hadrons up to ultra-high energies in vacuum electrostatic gaps
close to the black hole remains an open possibility that can trigger cascades
and high-energy gamma-ray emission in cases with very little accretion74, 75,
76.

The predominance of one type of composition or another can of course change
along jets, depending on the specific conditions. Forthcoming observations of
nearby jets with the EHT will put this picture to robust tests in the near
future.

## Acknowledgments

This work was supported by the Argentine agency CONICET (PIP 2014-00338), the
National Agency for Scientific and Technological Promotion (PICT 2017-2865)
and the Spanish Ministerio de Ciencia e Innovacion (MICINN) under grant
PID2019-105510GBC31 and though the “Center of Excellence Maria de Maeztu
2020-2023” award to the ICCUB (CEX2019-000918-M).

### Author contributions

This is a review article prepared wholly by the sole author.

### Financial disclosure

None reported.

### Conflict of interest

The author declares no potential conflict of interests.

## References

  * 1 Curtis H. D.. Descriptions of 762 Nebulae and Clusters Photographed with the Crossley Reflector. Publications of Lick Observatory. 1918;13:9-42.
  * 2 Begelman Mitchell C., Blandford Roger D., Rees Martin J.. Theory of extragalactic radio sources. Reviews of Modern Physics. 1984;56(2):255-351.
  * 3 Blandford Roger, Meier David, Readhead Anthony. Relativistic Jets from Active Galactic Nuclei. ARA&A. 2019;57:467-509.
  * 4 Beskin V. S.. MHD Flows in Compact Astrophysical Objects: Accretion, Winds and Jets.  2009\.
  * 5 Meier David L.. Black Hole Astrophysics: The Engine Paradigm.  2012\.
  * 6 Romero G. E., Vila G. S.. Introduction to Black Hole Astrophysics Lecture Notes in Physics, Berlin Springer Verlag, vol. 876: .  2014\.
  * 7 Shakura N. I., Sunyaev R. A.. Black holes in binary systems. Observational appearance.. A&A. 1973;500:33-51.
  * 8 Narayan Ramesh, Yi Insu. Advection-dominated Accretion: Self-Similarity and Bipolar Outflows. ApJ. 1995;444:231.
  * 9 Yuan Feng, Narayan Ramesh. Hot Accretion Flows Around Black Holes. ARA&A. 2014;52:529-588.
  * 10 Abramowicz M. A., Calvani M., Nobili L.. Thick accretion disks with super-Eddington luminosities. ApJ. 1980;242:772-788.
  * 11 Begelman M. C., Meier D. L.. Thick accretion disks - Self-similar, supercritical models. ApJ. 1982;253:873-896.
  * 12 Wiita P. J.. Physical properties of thick supercritical accretion disks. ApJ. 1982;256:666-680.
  * 13 Tout Christopher A., Pringle J. E.. Can a disc dynamo generate large-scale magnetic fields?. MNRAS. 1996;281(1):219-225.
  * 14 Liska M., Tchekhovskoy A., Quataert E.. Large-scale poloidal magnetic field dynamo leads to powerful jets in GRMHD simulations of black hole accretion with toroidal field. MNRAS. 2020;494(3):3656-3662.
  * 15 Blandford R. D., Payne D. G.. Hydromagnetic flows from accretion disks and the production of radio jets.. MNRAS. 1982;199:883-903.
  * 16 Koide Shinji, Meier David L., Shibata Kazunari, Kudoh Takahiro. General Relativistic Simulations of Early Jet Formation in a Rapidly Rotating Black Hole Magnetosphere. ApJ. 2000;536(2):668-674.
  * 17 Koide Shinji, Shibata Kazunari, Kudoh Takahiro, Meier David L.. Extraction of Black Hole Rotational Energy by a Magnetic Field and the Formation of Relativistic Jets. Science. 2002;295(5560):1688-1691.
  * 18 Meier David L., Koide Shinji, Uchida Yutaka. Magnetohydrodynamic Production of Relativistic Jets. Science. 2001;291(5501):84-92.
  * 19 Tchekhovskoy Alexander. Launching of Active Galactic Nuclei Jets Astrophysics and Space Science Library, vol. 414: :45.  2015\.
  * 20 Romero Gustavo, Gutiérrez Eduardo. The Origin of Matter at the Base of Relativistic Jets in Active Galactic Nuclei. Universe. 2020;6(7):99.
  * 21 Eichler David. Magnetic Confinement of Jets. ApJ. 1993;419:111.
  * 22 Globus N., Levinson A.. The collimation of magnetic jets by disc winds. MNRAS. 2016;461(3):2605-2615.
  * 23 Blandford R. D., Znajek R. L.. Electromagnetic extraction of energy from Kerr black holes.. MNRAS. 1977;179:433-456.
  * 24 Ghisellini G., Tavecchio F., Maraschi L., Celotti A., Sbarrato T.. The power of relativistic jets is larger than the luminosity of their accretion disks. Nature. 2014;515(7527):376-378.
  * 25 Tchekhovskoy Alexander, Narayan Ramesh, McKinney Jonathan C.. Black Hole Spin and The Radio Loud/Quiet Dichotomy of Active Galactic Nuclei. ApJ. 2010;711(1):50-63.
  * 26 Gutiérrez E. M., Vieyro F. L., Romero G. E.. Non-thermal processes in hot accretion flows onto supermassive black holes: an inhomogeneus model. A&A. 2021;in press.
  * 27 Event Horizon Telescope Collaboration , Akiyama Kazunori, Alberdi Antxon, et al. First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole. ApJ. 2019;875(1):L1.
  * 28 Hada Kazuhiro, Kino Motoki, Doi Akihiro, et al. The Innermost Collimation Structure of the M87 Jet Down to ~10 Schwarzschild Radii. ApJ. 2013;775(1):70.
  * 29 Hada Kazuhiro. The Structure and Propagation of the Misaligned Jet M87. Galaxies. 2017;5(1):2.
  * 30 Park Jongho, Hada Kazuhiro, Kino Motoki, et al. Kinematics of the M87 Jet in the Collimation Zone: Gradual Acceleration and Velocity Stratification. ApJ. 2019;887(2):147.
  * 31 Vila Gabriela S., Vieyro Florencia L., Romero Gustavo E.. Neutron Production in Black Hole Coronae and Proton Loading of Jets. In: International Journal of Modern Physics Conference Series, vol. 28: :1460191; 2014\.
  * 32 Toma K., Takahara F.. Baryon Loading of Active Galactic Nucleus Jets Mediated by Neutrons. ApJ. 2012;754(2):148.
  * 33 Kimura Shigeo S., Toma Kenji, Takahara Fumio. Effects of High-energy Particles on Accretion Flows onto a Supermassive Black Hole. ApJ. 2014;791(2):100.
  * 34 Ohsuga Ken, Mineshige Shin, Watarai Ken-ya. Spectral Energy Distribution in Supercritical Disk Accretion Flows through Photon-trapping Effects. ApJ. 2003;596(1):429-436.
  * 35 Ohsuga Ken, Mori Masao, Nakamoto Taishi, Mineshige Shin. Supercritical Accretion Flows around Black Holes: Two-dimensional, Radiation Pressure-dominated Disks with Photon Trapping. ApJ. 2005;628(1):368-381.
  * 36 Fukue Jun. Basic Properties of Supercritical Accretion Disks. PASJ. 2000;52:829-840.
  * 37 Fukue Jun. Critical Accretion Disk. PASJ. 2004;56:569-580.
  * 38 Bisnovatyi-Kogan G. S., Blinnikov S. I.. Disk accretion onto a black hole at subcritical luminosity.. A&A. 1977;59:111-125.
  * 39 Di Matteo T., Blackman E. G., Fabian A. C.. Two-temperature coronae in active galactic nuclei. MNRAS. 1997;291(1):L23-L27.
  * 40 Fabian A. C., Lohfink A., Kara E., Parker M. L., Vasudevan R., Reynolds C. S.. Properties of AGN coronae in the NuSTAR era. MNRAS. 2015;451(4):4375-4383.
  * 41 Vieyro F. L., Romero G. E.. Particle transport in magnetized media around black holes and associated radiation. A&A. 2012;542:A7.
  * 42 Romero G. E., Vieyro F. L., Vila G. S.. Non-thermal processes around accreting galactic black holes. A&A. 2010;519:A109.
  * 43 Vieyro F. L., Sestayo Y., Romero G. E., Paredes J. M.. Nonthermal processes and neutrino emission from the black hole GRO J0422+32 in a bursting state. A&A. 2012;546:A46.
  * 44 Kimura Shigeo S., Murase Kohta, Mészáros Peter. Multimessenger tests of cosmic-ray acceleration in radiatively inefficient accretion flows. Phys. Rev. D. 2019;100(8):083014.
  * 45 Inoue Yoshiyuki, Khangulyan Dmitry, Inoue Susumu, Doi Akihiro. On High-energy Particles in Accretion Disk Coronae of Supermassive Black Holes: Implications for MeV Gamma-rays and High-energy Neutrinos from AGN Cores. ApJ. 2019;880(1):40.
  * 46 McKinney Jonathan C., Blandford Roger D.. Stability of relativistic jets from rotating, accreting black holes via fully three-dimensional magnetohydrodynamic simulations. MNRAS. 2009;394(1):L126-L130.
  * 47 Araudo A. T., Bosch-Ramon V., Romero G. E.. Gamma rays from cloud penetration at the base of AGN jets. A&A. 2010;522:A97.
  * 48 del Palacio S., Bosch-Ramon V., Romero G. E.. Gamma rays from jets interacting with BLR clouds in blazars. A&A. 2019;623:A101.
  * 49 Araudo A. T., Bosch-Ramon V., Romero G. E.. High-energy emission from jet-clump interactions in microquasars. A&A. 2009;503:673-681.
  * 50 Komissarov S. S.. Mass-Loaded Relativistic Jets. MNRAS. 1994;269:394.
  * 51 Barkov M. V., Aharonian F. A., Bosch-Ramon V.. Gamma-ray Flares from Red Giant/Jet Interactions in Active Galactic Nuclei. ApJ. 2010;724:1517-1523.
  * 52 Bosch-Ramon V., Perucho M., Barkov M. V.. Clouds and red giants interacting with the base of AGN jets.. A&A. 2012;539:A69.
  * 53 Araudo A. T., Bosch-Ramon V., Romero G. E.. Gamma-ray emission from massive stars interacting with active galactic nuclei jets. MNRAS. 2013;436(4):3626-3639.
  * 54 Romero G. E., Torres D. F., Kaufman Bernadó M. M., Mirabel I. F.. Hadronic gamma-ray emission from windy microquasars. A&A. 2003;410:L1-L4.
  * 55 Owocki S. P., Romero G. E., Townsend R. H. D., Araudo A. T.. Gamma-Ray Variability from Wind Clumping in High-Mass X-Ray Binaries with Jets. ApJ. 2009;696(1):690-693.
  * 56 Perucho M., Bosch-Ramon V., Khangulyan D.. 3D simulations of wind-jet interaction in massive X-ray binaries. A&A. 2010;512:L4.
  * 57 Perucho M., Bosch-Ramon V.. 3D simulations of microquasar jets in clumpy stellar winds. A&A. 2012;539:A57.
  * 58 Fabrika S.. The jets and supercritical accretion disk in SS433. Astrophys. Space Phys. Res.. 2004;12:1-152.
  * 59 Gallo Elena, Fender Rob, Kaiser Christian, et al. A dark jet dominates the power output of the stellar black hole Cygnus X-1. Nature. 2005;436(7052):819-821.
  * 60 Heinz S.. Composition, Collimation, Contamination: The Jet of Cygnus X-1. ApJ. 2006;636(1):316-322.
  * 61 Bromberg Omer, Granot Jonathan, Piran Tsvi. On the composition of GRBs’ Collapsar jets. MNRAS. 2015;450(1):1077-1084.
  * 62 Croston J. H., Ineson J., Hardcastle M. J.. Particle content, radio-galaxy morphology, and jet power: all radio-loud AGN are not equal. MNRAS. 2018;476(2):1614-1623.
  * 63 Heinz Sebastian. Clues for the Composition of Relativistic Microquasar Jets. International Journal of Modern Physics D. 2008;17(10):1947-1952.
  * 64 Celotti A.. Composition and energetics of relativistic jets. Ap&SS. 2003;288(1):175-181.
  * 65 Sikora M., Madejski G.. On Pair Content and Variability of Subparsec Jets in Quasars. ApJ. 2000;534(1):109-113.
  * 66 Madejski Grzegorz (Greg), Sikora Marek. Gamma-Ray Observations of Active Galactic Nuclei. ARA&A. 2016;54:725-760.
  * 67 Wardle J. F. C., Homan D. C., Ojha R., Roberts D. H.. Electron-positron jets associated with the quasar 3C279. Nature. 1998;395(6701):457-461.
  * 68 Homan Daniel, Hovatta Talvikki, Kovalev Yuri, Lister Matthew, Pushkarev Alexander, Savolainen Tuomas. Constraints on Particles and Fields from Full Stokes Observations of AGN. Galaxies. 2018;6(1):17.
  * 69 Anantua Richard, Emami Razieh, Loeb Abraham, Chael Andrew. Determining the Composition of Relativistic Jets from Polarization Maps. ApJ. 2020;896(1):30.
  * 70 Komissarov S. S.. Numerical simulations of relativistic magnetized jets. MNRAS. 1999;308(4):1069-1076.
  * 71 Komissarov S. S., Falle S. A. E. G.. The evolution of slab jets and self-similar models of extragalactic radio sources. MNRAS. 2003;343(3):1045-1053.
  * 72 Martí Josep, Luque-Escamilla Pedro L., Romero Gustavo E., Sánchez-Sutil Juan R., Muñoz-Arjonilla Álvaro J.. Real-time evolution of a large-scale relativistic jet. A&A. 2015;578:L11.
  * 73 Martí Josep, Luque-Escamilla Pedro L., Bosch-Ramon Valentí, Paredes Josep M.. A galactic microquasar mimicking winged radio galaxies. Nature Communications. 2017;8:1757.
  * 74 Levinson Amir. Particle Acceleration and Curvature TeV Emission by Rotating, Supermassive Black Holes. Phys. Rev. Lett.. 2000;85(5):912-915.
  * 75 Levinson Amir, Rieger Frank. Variable TeV Emission as a Manifestation of Jet Formation in M87?. ApJ. 2011;730(2):123.
  * 76 Levinson Amir, Segev Noam. Existence of steady gap solutions in rotating black hole magnetospheres. Phys. Rev. D. 2017;96(12):123006.


∎

11institutetext:  I.A. Korneev, A.V. Slepnev, T.E. Vadivasova 22institutetext:
Saratov State University, Astrakhanskaya str., 83, 410012, Saratov, Russia
33institutetext: A.S. Zakharova 44institutetext: Technische Universität
Berlin, Hardenbergstraße 36, 10623, Berlin, Germany
55institutetext: V.V. Semenov 66institutetext: Technische Universität Berlin,
Hardenbergstraße 36, 10623, Berlin, Germany
Saratov State University, Astrakhanskaya str., 83, 410012, Saratov, Russia
66email<EMAIL_ADDRESS>
# Generalized model for steady-state bifurcations without parameters in
memristor-based oscillators with lines of equilibria

lvan A. Korneev  Andrei V. Slepnev  Anna S. Zakharova  Tatiana E. Vadivasova
Vladimir V. Semenov

(Received: date / Accepted: date)

###### Abstract

We demonstrate how the pitchfork, transcritical and saddle-node bifurcations
of steady states observed in dynamical systems with a finite number of
isolated equilibrium points occur in systems with lines of equilibria. The
exploration is carried out by using the numerical simulation and linear
stability analysis applied to a model of a memristor-based oscillator. First,
all the discussed bifurcation scenarios are considered in the context of
systems including Chua’s memristor with a piecewise-smooth characteristic.
Then the memristor characteristic is changed to a function that is smooth
everywhere. Finally, the action of the memristor forgetting effect is taken
into consideration. The presented results are obtained for electronic circuit
models, but the considered bifurcation phenomena can be exhibited by systems
with a line of equilibria of any nature.

###### Keywords:

memristor memristor-based oscillators line of equilibria bifurcations without
parameters pitchfork bifurcation transcritical bifurcation saddle-node
bifurcation

###### pacs:

05.10.-a 05.45.-a 84.30.-r

##  1 Introduction

The simultaneous and continuous dependence of the oscillatory dynamics both on
parameter values and initial conditions is a frequent occurrence in various
dynamical systems. In most cases, this peculiarity is associated with the
properties of conservative oscillators. However, such behaviour can be
exhibited by nonlinear dissipative oscillators with manifolds of equilibria.
The manifolds of equilibria consist of non-isolated equilibrium points and
takes different forms: line of equilibria fiedler2000-1 ; fiedler2000-2 ;
fiedler2000-3 ; liebscher2015 , surface of equilibria jafari2016 , circle
gotthans2015 ; gotthans2016  and square gotthans2016  equilibria, etc. Among
such manifolds one distinguishes $m$-dimensional normally hyperbolic manifolds
of equilibria characterized by $m$ pure imaginary eigenvalues or $m$
eigenvalues being equal to zero, whereas all the other eigenvalues have non-
zero real parts. In the simplest case these manifolds exist as a line of
equilibria. Systems with a line of equilibria have been mathematically
considered fiedler2000-1 ; fiedler2000-2 ; fiedler2000-3 ; liebscher2015 ;
riaza2012 ; riaza2016 ; corinto2020 . It has been shown that their significant
feature is the occurrence of so-called bifurcations without parameters, i.e.,
the bifurcations corresponding to fixed parameters when the condition of
normal hyperbolicity is violated at some points of the line of equilibria.

Memristor-based oscillators are widely represented in a variety of systems
with a line of equilibria. One of the frequently considered models is the
series RLC-circuit with negative resistance, where a flux-controlled memristor
is connected in parallel with the capacitor. Bifurcation mechanisms of the
periodic solution appearance in this system have been explored numerically and
analytically for different kinds of nonlinearity messias2010 ; botta2011 ;
riaza2012 ; semenov2015 ; korneev2017 ; korneev2017-2 ; korneev2021-2 .
Analytical solution of the model equations has enabled to recognize
significant features of the supercritical korneev2017 ; korneev2017-2  and
subcritical korneev2021-2  Andronov-Hopf bifurcation. In addition, the
recently published research korneev2021-2  has demonstrated how the saddle-
node bifurcation of limit cycles observed in classical self-oscillators with
hard self-oscillation excitation transforms in systems with a line of
equilibria. Thus, all the bifurcations of limit cycles exhibited by self-
oscillators with one degree of freedom have been described in terms of
oscillators with a line of equilibria. The similarity between classical self-
oscillators and the considered systems with a line of equilibria is emphasized
by the fact that undamped periodic oscillations in systems with a line of
equilibria can be synchronized by external periodic forcing korneev2020 .

In contrast to our previous publications focused on the self-oscillatory
dynamics, in the current paper we show how classical bifurcations of steady
states are manifested in systems with a line of equilibria. Three bifurcations
are studied: the pitchfork bifurcation, the transcritical bifurcation and the
saddle-node bifurcation of steady states. It is demonstrated that the
bifurcations under study become bifurcations without parameters due to the
memristor properties. Our purpose is to describe all the models and
bifurcations in a unified manner. For this reason, all the considered systems
belong to the same class of circuits, the RLC-circuit with feedback including
the memristor.

##  2 Model and methods

Figure 1: (a) Memristive circuit based on the memristive non-inverting
amplifier and analog divider; (b) Schematic circuit diagram of the studied
model (Eqs.(6)).

According to the idea proposed by Leon Chua chua1971  the memristor relates
the transferred electrical charge, $q(t)$, and the magnetic flux linkage,
$\varphi(t)$, by means of the linear relationship $dq=Wd\varphi$, whence it
follows that $W=W(\varphi)=\dfrac{dq}{d\varphi}$. By this way, using the
relationships $d\varphi=V_{\text{m}}dt$ and $dq=I_{\text{m}}dt$
($V_{\text{m}}$ is the voltage across the memristor, $I_{\text{m}}$ is the
current passing through the memristor), the memristor current-voltage
characteristic can be derived: $I_{\text{m}}=W(\varphi)V_{\text{m}}$. In such
a case, $W$ is the flux-controlled conductance (memductance) and depends on
the entire past history of $V_{\text{m}}(t)$:

$W(\varphi)=\dfrac{dq}{d\varphi}=q^{\prime}\left(\int\limits_{-\infty}^{t}{V_{\text{m}}(t)dt}\right).$
(1)

It must be noted that the term ’flux’ (or ’flux linkage’) is used below only
for denoting the memristor state variable being proportional to the integral
$\int\limits_{-\infty}^{t}{V_{\text{m}}(t)dt}$. The issues concerning the
physical realizability of the postulated relationship $dq=Wd\varphi$ are not
discussed in the current paper. Thus, the memristor is considered as a
resistive component which conductance is dictated by the state variable which
is not necessarily associated with magnetic phenomena. This approach reflects
the conception of ’memristive system’ chua1976  implying the mathematical
definition, which does not concern a physical sense of the dynamical variables
and their functional dependence. It allows for grouping a broad variety of
elements of different nature identified by a continuous functional dependence
of characteristics on previous states.

One of the simplest memristor model is Chua’s memristor described by the
piecewise-linear dependence $q(\varphi)$ which takes the following form for
the flux-controlled memristor:

$q(\varphi)=\begin{cases}(a-b)\varphi_{*}+b\varphi,&\varphi\geq\varphi_{*},\\\
a\varphi,&|\varphi|<\varphi_{*},\\\
-(a-b)\varphi_{*}+b\varphi,&\varphi\leq-\varphi_{*}.\end{cases}$ (2)

Then the memristor conductance $W(\varphi)$ becomes:

$W(\varphi)=\begin{cases}a,&|\varphi|<\varphi_{*},\\\
b,&|\varphi|\geq\varphi_{*}.\end{cases}$ (3)

Nonlinearity (3) can be approximated by the hyperbolic tangent function:

$W(\varphi)=\dfrac{b-a}{2}\tanh\left(k(\varphi^{2}-\varphi_{*})\right)+\dfrac{b+a}{2},$
(4)

where a parameter $k$ characterizes the sharpness of the transitions between
two memristor’s states. It has been shown in korneev2021-2  that changing the
memristor conductance function to the smooth one does not qualitatively modify
the memristor properties. The classical loop in the current-voltage
characteristic of the memristor driven by the external periodic influence
persists (see Fig. 2 in Ref. korneev2021-2 ). The memristor model including
tanh-nonlinearity is not the only smooth memristor model. There is a number of
smooth models describing various memristor properties tetzlaff2014 ; linn2014
; singh2019 ; ascoli2013-2 ; chang2011 ; chua2011 ; guseinov2021 .

Real memristive systems can ’forget’ the state history over time. In
particular, the ’forgetting’ effect in memristors based on metal oxides is
associated with the diffusion of charged particles chang2011 ; chen2013 ;
zhou2019  (however, the ’forgetting’ can happen very slowly). One of the
simplest form of the memristor state equation which implies the forgetting
effect is the following:

$\dfrac{dz}{dt}=g(x,z)=x-\delta z,$ (5)

where $z$ plays a role of a memristor state variable, $x$ is an input signal,
a parameter $\delta$ characterizes the forgetting effect strength.

Consider the circuit in Fig. 1 (a). It contains the operational-amplifier-
based non-inverting memristive amplifier which output signal is
$V_{a}=(1+R_{1}W(\varphi))V_{\text{in}}$, where $\varphi$ varies according to
memristor state equation (5) taken in the form
$\dfrac{d\varphi}{dt}=g(V_{\text{in}},\varphi)$. The second block is an analog
divider producing the output voltage such that the resulting circuit response
is $V_{\text{out}}=1+Z+R_{1}W(\varphi)$, where $Z$ is a summing input voltage.
The memristive circuit in Fig. 1 (a) is included as the memristive block $M$
into the schematic circuit diagram in Fig. 1 (b). The circuit in Fig. 1 (b)
represents the linear series RLC-circuit forced by the feedback signal
$V_{\text{f}}(t)$ being an output signal of the complex nonlinear amplifier
(the grey block in Fig. 1 (b)). The amplifier has two inputs: the voltage drop
across the resistor, $V_{R}=RI$ and the voltage across the capacitor, $V$. The
amplifier input currents are assumed to be zero. After the voltage $V_{R}$ is
reproduced by a differential amplifier, it becomes an input signal of the
memristive circuit $M$ which output $V_{\text{M}}$ is comes to the analog
multiplier and to the two-state switch. If the switch state is 1, the
corresponding output signal is $V_{\text{s}}=V_{\text{M}}V$. Otherwise, the
signal coming to the summing block $\sum$ is $V_{\text{s}}=V_{\text{M}}$. In
addition, the amplifier contains the nonlinear block for the square or cubic
transformation $sV^{n}$ of the input. Thus, the amplifier output signal is
$V_{\text{f}}=V_{\text{s}}+V+mV^{n}$. The presented in Fig. 1 (b) system is
described by the following dynamical variables: $V$ is the voltage across the
capacitor $C$, $I$ is the current through the inductor $L$ and $\varphi$ is
the magnetic flux linkage controlling the memristor. Using Kirchhoff’s laws,
one obtains differential equations for the considered system in physical time
$t^{\prime}$:

$\left\\{\begin{array}[]{l}C\dfrac{dV}{dt^{\prime}}=I,\\\ \\\
L\dfrac{dI}{dt^{\prime}}=-RI+V_{\text{s}}(V,V_{R},\varphi)+mV^{n},\\\ \\\
\dfrac{d\varphi}{dt^{\prime}}=RI-k\varphi,\end{array}\right.$ (6)

where the third equation reflects the memristor forgetting effect. In the
dimensionless variables $x=V/V_{0}$, $Y=I/I_{0}$ and
$z=\varphi/(L\varphi_{0})$ with $V_{0}=1$ V, $I_{0}=1$ A, $\varphi_{0}=\text{1
sec}\times V_{0}$ and the dimensionless time $t=[(V_{0}/(I_{0}L)]t^{\prime}$,
Eqs.(6) can be rewritten in the following form:

$\left\\{\begin{array}[]{l}\dfrac{dx}{dt}=\nu Y,\\\ \\\ \dfrac{dY}{dt}=-\gamma
Y+f(x,Y,z)+\mu x^{n},\\\ \\\ \dfrac{dz}{dt}=\gamma Y-\delta z,\\\
\end{array}\right.$ (7)

where $\nu=(L/C)(I_{0}/V_{0})^{2}$ is a dimensionless parameter being
numerically equal to $L/C$, $\gamma=R(I_{0}/V_{0})$ is a dissipation factor,
$f(x,Y,z)$ is a function reflecting the memristive properties of the block $M$
in Fig. 1 (d) and represents a dimensionless analog of the signal
$V_{\text{s}}$, $\mu$ is a dimensionless equivalent of the parameter $m$,
$\delta=kL(I_{0}/V_{0})$. The substitution $y=\gamma Y$ transforms model (7)
into the finalized form:

$\left\\{\begin{array}[]{l}\dfrac{dx}{dt}=\alpha y,\\\ \\\
\dfrac{1}{\gamma}\dfrac{dy}{dt}=-y+f(x,y,z)+\mu x^{n},\\\ \\\
\dfrac{dz}{dt}=y-\delta z,\\\ \end{array}\right.$ (8)

where $\alpha=\nu/\gamma$. The function $f(x,y,z)$ contains the functional
dependence on the memristor instantaneous state. In physical variables the
dependence is described as the expression for the memristor conductance given
either by Exps. (3) or (4). In the context of model (8), the function
$W_{M}(\varphi)$ corresponds to the equivalent dimensionless form $G_{M}(z)$.
Two options for $G_{M}(z)$ are under consideration:

$G_{M}(z)=\begin{cases}a,&|z|<z_{0},\\\ b,&|z|\geq z_{0},\end{cases}$ (9)

and

$G_{M}(z)=\dfrac{b-a}{2}\tanh\left(k(z^{2}-z_{0})\right)+\dfrac{b+a}{2}.$ (10)

System (8) is explored both theoretically by using the linear stability
analysis and numerically by means of integration methods. Numerical
simulations are carried out by integration of Eqs. (8) using the fourth-order
Runge-Kutta method with the time step $\Delta t=0.0001$ from different initial
conditions. Particular modifications of dynamical model (8) are considered
below in more details.

##  3 Pitchfork bifurcation

Suppose that the switch in Fig. 1 (b) operates in state 1 while the nonlinear
transformation is cubic. Then the expression for the feedback voltage takes
the form $V_{\text{f}}=(1+Z+R_{1}W(\varphi))V+V+mV^{3}$, where the state
variable varies according the equation $d\varphi/dt=RI-k\varphi$ and $m=-1$.
Adjusting the voltage $Z$ and the resistance $R_{1}$, one controls the
coefficient $(1+Z+R_{1}W(\varphi))$ which can possess both negative and
positive values. Then the physical model (see Eqs. (6)) corresponding to the
presence of piecewise-smooth memristor (9) or its continuous analog (10) is
transformed to the dimensionless form

$\left\\{\begin{array}[]{l}\dfrac{dx}{dt}=\alpha y,\\\ \\\
\dfrac{1}{\gamma}\dfrac{dy}{dt}=-y+\beta(z)x-x^{3},\\\ \\\
\dfrac{dz}{dt}=y-\delta z,\\\ \end{array}\right.$ (11)

for two options $\beta(z)$:

$\beta(z)=\begin{cases}\beta_{1},&|z|<1,\\\ \beta_{2},&|z|\geq 1,\end{cases}$
(12)
$\beta(z)=\dfrac{\beta_{2}-\beta_{1}}{2}\tanh\left(k(z^{2}-1)\right)+\dfrac{\beta_{2}+\beta_{1}}{2}.$
(13)

System (11) is studied for fixed parameters: $\alpha=\gamma=1$,
$\beta_{1}=-1$, $\beta_{2}=1$, $k=5$. First, the memristor forgetting effect
is excluded from the consideration: $\delta=0$.

Figure 2: Pitchfork bifurcation without parameters in the phase space of
system (11) for memristive nonlinearities (12) (panel (a)) and (13) (panel
(b)) and the absence of the forgetting effect, $\delta=0$. The attracting
lines of equilibria are marked by solid green and red lines, the unstable
branches of the central line of equilibria are marked by green dashed lines.
Phase trajectories are black arrowed curves. The right insets schematically
show the manifolds of equilibria as a projection on the plane $(x,z)$; (c)
Phase trajectories in system (11)-(13) for $\delta=10^{-3}$ tracing motions to
the stable steady state (green point) in the origin. Other system parameters
are: $\alpha=\gamma=1$, $\beta_{1}=-1$,$\beta_{2}=1$, $k=5$.

System (11) with piecewise-smooth function (12) does not contain the variable
$z$ in an explicit form. Then two cases $\beta=\beta_{1}$ and
$\beta=\beta_{2}$ are considered independently. In case $|z|<1$ the system
attractor is a manifold of steady states with the coordinates $x_{1}=0$,
$y_{1}=0$, $z\in(-1;1)$ (the green solid line in Fig. 2 (a)). All the
trajectories in the subspace $z\in(-1;1)$ are attracted to the manifold of
equilibria. Thus, the attractor represents a line of equilibria characterised
by the eigenvalues $\lambda_{1}=0$, $\lambda_{2,3}=(-1\pm i\sqrt{3})/2$
considered as stable fixed points since $\lambda_{2,3}$ have negative real
part. Each point of the line of equilibria is neutrally stable in the OZ-axis
direction. Hereinafter, using the terms ’stable’ or ’unstable’ point at the
line of equilibria, we mean the behaviour of trajectories (attraction or
repelling) in the neighbourhood of the equilibrium point, which is determined
by the eigenvalues $\lambda_{2,3}$.

For two subspaces $|z|\geq 1$, there exist three lines of equilibria. The
first one ($x_{1}=0$, $y_{1}=0$, $|z|\geq 1$) represent a line of saddle-like
equilibria (the green dashed line in Fig. 2 (a)) with eigenvalues
$\lambda_{1}=0$, $\lambda_{2,3}=(-1\pm\sqrt{5})/2$. The second and third lines
of equilibria (the red solid lines in Fig. 2 (a)) are two coexisting
attractors which consist of steady states ($x_{2,3}=\pm 1$, $y_{2,3}=0$,
$|z|\geq 1$) with the corresponding eigenvalues $\lambda_{1}=0$,
$\lambda_{2,3}=(-1\pm i\sqrt{7})/2$. Thus, for $|z|\geq 1$ the central line of
equilibria plays a role of an unstable fixed points between two attractors in
classical bistable oscillators.

In summary, changing $z$ and keeping parameter values to be fixed, one
implements the transition between the existence of a single attractor to the
regime of bistablility which consists in the coexistence of two stable lines
of equilibria. The described bifurcation transition corresponds to the
pitchfork bifurcation which occurs at $z=\pm 1$: the line of equilibria
becomes unstable for $|z|\geq 1$ and two coexisting attractors appear at the
moment of bifurcation. If system (11) involves tanh-function (13), the
similarity with the classical pitchfork bifurcation is more evident [Fig. 2
(b)]: there exists a single line of equilibria in the central subspace, which
loses stability when $z$ passes through the threshold value $z_{*}\approx\pm
1$. At the bifurcation moment two stable lines of equilibria begin in the
vicinity of the unstable line of equilibria.

For non-zero values of the parameter $\delta$ the lines of equilibria
disappear [Fig. 2 (c)]. In such a case, the steady state $x=y=z=0$ becomes a
single attractor in the phase space. Then the transition to the bistability
occurs in a classical way by varying parameter $\beta_{1}$: the pitchfork
bifurcation occurs at $\beta_{1}=0$. However, for small values $\delta$ the
phase trajectories describing motion to the stable fixed point can trace the
lines of equilibria initially existed at $\delta=0$ (see the trajectories in
Fig. 2 (c)).

##  4 Transcritical bifurcation

Suppose that the switch in Fig. 1 (b) operates in state 1 while the nonlinear
transformation is quadratic. Then the expression for the feedback voltage
becomes $V_{\text{f}}=(1+Z+R_{1}W(\varphi))V+V+mV^{2}$, where the state
variable varies according the equation $d\varphi/dt=RI-k\varphi$ and $m=-1$.
Following the same procedures as in the previous section, one obtains the
dimensionless model of the circuit:

$\left\\{\begin{array}[]{l}\dfrac{dx}{dt}=\alpha y,\\\ \\\
\dfrac{1}{\gamma}\dfrac{dy}{dt}=-y+\beta(z)x-x^{2},\\\ \\\
\dfrac{dz}{dt}=y-\delta z,\\\ \end{array}\right.$ (14)

where the function $\beta(z)$ has two options: Exp. (12) or Exp. (13). System
(14) is studied for fixed parameters: $\alpha=\gamma=1$, $\beta_{1}=-1$,
$\beta_{2}=1$, $k=5$.

First, system (14) is considered for piecewise-smooth function $\beta(z)$
(function (12)) and $\delta=0$. Two cases $\beta=\beta_{1}$ and
$\beta=\beta_{2}$ are considered independently. In case $|z|<1$ the system
attractor is a manifold of steady states with the coordinates $x_{1}=0$,
$y_{1}=0$, $z\in(-1;1)$ characterised by the eigenvalues $\lambda_{1}=0$,
$\lambda_{2,3}=(-1\pm i\sqrt{3})/2$ (the red solid line in Fig. 3 (a)). All
the phase trajectories in the subspace $z\in(-1;1)$ are attracted to the line
of equilibria $x_{1}=0$, $y_{1}=0$, $z\in(-1;1)$. In addition, there exists
the second line of equilibria ($x_{2}=-1$, $y_{2}=0$, $z\in(-1;1)$) (the green
dashed line in Fig. 3 (a)) exhibiting the properties of a saddle fixed point
(the corresponding eigenvalues are $\lambda_{1}=0$,
$\lambda_{2,3}=(-1\pm\sqrt{5})/2$).

Figure 3: Transcritical bifurcation without parameters in the phase space of
system (14) for memristive nonlinearities (12) (panel (a)) and (13) (panel
(b)) and the absence of the forgetting effect, $\delta=0$. The attracting
segments of the lines of equilibria are marked by solid red lines while the
repelling ones are marked by green dashed lines. Phase trajectories are black
arrowed curves. The right insets schematically show the manifolds of
equilibria as a projection on the plane $(x,z)$; (c) Phase trajectories in
system (14)-(13) for $\delta=10^{-3}$ tracing motions to the stable steady
state (green point). Other system parameters are: $\alpha=\gamma=1$,
$\beta_{1}=-1$,$\beta_{2}=1$, $k=5$.

For two subspaces $|z|\geq 1$, there also exist two lines of equilibria. The
first one ($x_{1}=0$, $y_{1}=0$, $|z|\geq 1$) represent a line of saddle-like
equilibria (the green dashed lines in Fig. 3 (a)) with eigenvalues
$\lambda_{1}=0$, $\lambda_{2,3}=(-1\pm\sqrt{5})/2$. The second line of
equilibria (the red solid lines in Fig. 3 (a)) is an attractor and consists of
steady states ($x_{2}=1$, $y_{2}=0$, $|z|\geq 1$) with the corresponding
eigenvalues $\lambda_{1}=0$, $\lambda_{2,3}=(-1\pm i\sqrt{3})/2$.

In conclusion, changing $z$ and keeping parameter values to be fixed one
implements the bifurcation transition at $z=\pm 1$ when the central segment of
the line of equilibria ($x_{1}=0$, $y_{1}=0$) changes the stability at the
same moment with the coexisting lines of equilibria ($x_{2}=\beta$,
$y_{2}=0$). If system (14) involves tanh-function (13), the described
transition is transformed into the transcritical bifurcation without
parameters [Fig. 3 (b)]: two lines of equilibria ($x_{1}=0$, $y_{1}=0$,
$z\in(-\infty;\infty)$) and ($x_{2}=\beta$, $y_{2}=0$,$z\in(-\infty;\infty)$)
intersect at $z=\pm 1$ and exchange the stability.

For non-zero values of the parameter $\delta$ the lines of equilibria
disappear [Fig. 3 (c)]. In such a case, the steady state $x=y=z=0$ becomes a
single attractor in the phase space and the transcritical bifurcation can be
induced only by varying the system parameter $\beta_{1}$. However, for small
values $\delta$ the phase trajectories describing motion to the stable fixed
point can trace the attracting segments of the lines of equilibria initially
existed at $\delta=0$ (see the trajectories in Fig. 3 (c)).

##  5 Saddle-node bifurcation

Suppose that the switch in Fig. 1 (b) operates in state 2 while the nonlinear
transformation is quadratic. Then the expression for the feedback voltage
becomes $V_{\text{f}}=(1+Z+R_{1}W(\varphi))+V+mV^{2}$, where the state
variable varies according the equation $d\varphi/dt=RI-k\varphi$ and $m=-1$.
Then the dimensionless model of the circuit takes the form:

$\left\\{\begin{array}[]{l}\dfrac{dx}{dt}=\alpha y,\\\ \\\
\dfrac{1}{\gamma}\dfrac{dy}{dt}=-y+\beta(z)-x^{2},\\\ \\\
\dfrac{dz}{dt}=y-\delta z,\\\ \end{array}\right.$ (15)

where the function $\beta(z)$ has two options: Exp. (12) or Exp. (13). System
(15) is studied for fixed parameters: $\alpha=\gamma=1$, $\beta_{1}=-1$,
$\beta_{2}=1$, $k=5$.

First, system (15) is considered for piecewise-smooth function $\beta(z)$
(function (12)) and $\delta=0$. Then in case $|z|<1$ the system has no
equilibria while two lines of equilibria ($x_{1,2}=\pm\sqrt{\beta_{2}}$,
$y_{1,2}=0$, $|z|\geq 1$) exist for two subspaces $|z|\geq 1$ [Fig. 4 (a)].
The first line of equilibria ($x_{1}=1$, $y_{1}=0$, $|z|\geq 1$) characterised
by the eigenvalues $\lambda_{1}=0$, $\lambda_{2,3}=(-1\pm i\sqrt{7})/2$ is the
system attractor. The second line of equilibria ($x_{1}=-1$, $y_{1}=0$,
$|z|\geq 1$) exhibits the saddle-equilibrium properties (the corresponding
eigenvalues are $\lambda_{1}=0$, $\lambda_{2,3}=(-1\pm 3)/2$).

Figure 4: Saddle-node bifurcation without parameters in the phase space of
system (15) for memristive nonlinearities (12) (panel (a)) and (13) (panel
(b)) and the absence of the forgetting effect, $\delta=0$. The attracting
segments of the lines of equilibria are marked by solid red lines while the
repelling ones are marked by green dashed lines. Phase trajectories are black
arrowed curves. The right insets schematically show the manifolds of
equilibria as a projection on the plane $(x,z)$; (c) Phase trajectory in
system (15)-(13) for $\delta=10^{-3}$ ending at infinity. Other system
parameters are: $\alpha=\gamma=1$, $\beta_{1}=-1$,$\beta_{2}=1$, $k=5$.

In summary, changing $z$ and keeping parameter values to be fixed one
implements the bifurcation transition at $z=\pm 1$ when two equilibria ($x=\pm
1$, $y=0$) appear: stable and unstable ones. If system (15) involves tanh-
function (13), the described transition is transformed into the saddle-node
bifurcation without parameters [Fig. 4 (b)]: two lines of equilibria ($x=\pm
1$, $y=0$, $|z|\geq 1$) appear at $z=\pm 1$. The first line of equilibria is
an attractor while the second one is unstable.

For non-zero values of the parameter $\delta$ the lines of equilibria
disappear. In such a case, the system has no equilibria. As demonstrated in
Fig. 4 (c), the phase trajectory can trace the attracting lines of equilibria
initially existed at $\delta=0$, but after that the trajectories come to
infinity.

## Conclusions

The circuit model proposed in the current paper represents a universal
dynamical system for the implementation of bifurcations without parameters
related to basic bifurcations of steady states: the pitchfork, transcritical
and saddle-node bifurcations. The demonstrated bifurcations are associated
with the existence of lines of equilibria. The transcritical and saddle-node
bifurcations imply the coexistence of a stable and unstable line of
equilibria. The pitchfork bifurcation realises the transition to the regime of
bistability which consists in the coexistence of two lines of equilibria
attracting the phase trajectories and an unstable line of equilibria between
them.

The impact of the memristor forgetting effect is manifested in the contraction
of limit sets along the OZ-axis. Then the resulting dynamics is determined by
the limit sets crossing the plane $z=0$ in the phase space. If any line of
equilibria intersect the plane $z=0$, it transforms into an isolated fixed
point with the same stability properties. Resultantly, a continuous dependence
of the oscillation characteristics on the initial condition $z_{0}$
disappears. A similar character of the forgetting effect action was reported
in earlier publicationskorneev2020 ; korneev2021 ; korneev2021-2 .

The transcritical bifurcation without parameter was discussed before riaza2018
. However, we are not aware of any models of dynamical systems exhibiting the
pitchfork and saddle-node bifurcations without parameters. Apparently, the
current paper is the first publication where the occurrence of the pitchfork
and saddle-node bifurcations without parameters is analysed.

###### Acknowledgements.

V.V.S. and A.S.Z. acknowledge support by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) – Projektnummer – 163436311-SFB-910.

## References

  * (1) Ascoli, A., Senger, V., Tetzlaff, R.: Memristor model comparison.  IEEE Circuits and Systems Magazine 13(2), 89–105 (2013)
  * (2) Botta, V., Néspoli, C., Messias, M.: Mathematical analysis of a third-order memristor-based Chua’s oscillator.  TEMA Tend. Mat. Apl. Comput. 12(2), 91–99 (2011)
  * (3) Chang, T., Jo, S., Kim, K., Sheridan, P., Gaba, S., Lu, W.: Synaptic behaviors and modeling of a metal oxide memristive device.  Applied Physics A 102(4), 857–863 (2011)
  * (4) Chen, L., Li, C., Huang, T., Chen, Y., Wen, S., Qi, J.: A synapse memristor model with forgetting effect.  Phys. Lett. A 377(45-48), 3260–3265 (2013)
  * (5) Chua, L.: Memristor – the missing circuit element.  IEEE Trans. on Circuit Theory CT-18(5), 507–519 (1971)
  * (6) Chua, L.: Resistance switching memories are memristors.  Appl. Phys. A 102(4), 765–783 (2011)
  * (7) Chua, L., Kang, S.: Memristive devices and systems.  Proceedings of the IEEE 64(2), 209–223 (1976)
  * (8) Corinto, F., Forti, M., Chua, L.: Nonlinear Circuits and Systems with Memristors.  Springer (2020)
  * (9) Fiedler, B., Liebscher, S.: Hopf bifurcation from lines of equilibria without parameters: II. systems of viscous hyperbolic balance laws.  SIAM J. Math. Anal. 31(6), 1396–1404 (2000)
  * (10) Fiedler, B., Liebscher, S., Alexander, J.: Generic Hopf bifurcation from lines of equilibria without parameters: I. theory.  Journal of Differential Equations 167(1), 16–35 (2000)
  * (11) Fiedler, B., Liebscher, S., Alexander, J.: Generic Hopf bifurcation from lines of equilibria without parameters: III. binary oscillators.  International Journal of Bifurcation and Chaos 10(7), 1613–1621 (2000)
  * (12) Gotthans, T., Petržela, J.: New class of chaotic systems with circular equilibrium.  Nonlinear Dynamics 81(3), 1143–1149 (2015)
  * (13) Gotthans, T., Sprott, J., Petrzela, J.: Simple chaotic flow with circle and square equilibrium.  International Journal of Bifurcation and Chaos 26(8), 1650,137 (2016)
  * (14) Guseinov, D., Matyushkin, I., Chernyaev, N., Mikhailov, A., Pershin, Y.: Capacitive effects can make memristors chaotic.  Chaos, Solitons and Fractals 144, 110,699 (2021)
  * (15) Jafari, S., Sprott, J., Pham, V.T., Volos, C., Li, C.: Simple chaotic 3d flows with surfaces of equilibria.  Nonlinear Dynamics 86(2), 1349–1358 (2016)
  * (16) Korneev, I., Semenov, V.: Andronov–Hopf bifurcation with and without parameter in a cubic memristor oscillator with a line of equilibria.  Chaos 27(8), 081,104 (2017)
  * (17) Korneev, I., Semenov, V., Slepnev, A., Vadivasova, T.: Complete synchronization of chaos in systems with nonlinear inertial coupling.  Chaos, Solitons and Fractals 142, 110,459 (2021)
  * (18) Korneev, I., Slepnev, A., Vadivasova, T., Semenov, V.: Forced synchronization of an oscillator with a line of equilibria.  Eur. Phys. J. Special Topics 229(12), 2215–2224 (2020)
  * (19) Korneev, I., Slepnev, A., Vadivasova, T., Semenov, V.: Subcritical Andronov–Hopf scenario for systems with a line of equilibria.  Chaos 31(7), 073,102 (2021)
  * (20) Korneev, I., Vadivasova, T., Semenov, V.: Hard and soft excitation of oscillations in memristor-based oscillators with a line of equilibria.  Nonlinear Dynamics 89(4), 2829–2843 (2017)
  * (21) Liebscher, S.: Bifurcation without Parameters, _Lectures Notes in Mathematics_ , vol. 2117.  Springer International Publishing (2015)
  * (22) Linn, E., Siemon, A., Waser, R., Menzel, S.: Applicability of well-established memristive models for simulations of resistive switching devices.  IEEE Transactions on Circuits and Systems I: Regular Papers 61(8), 2402–2410 (2014)
  * (23) Messias, M., Nespoli, C., Botta, V.: Hopf bifurcation from lines of equilibria without parameters in memristor oscillators.  International Journal of Bifurcation and Chaos 20(2), 437–450 (2010)
  * (24) Riaza, R.: Manifolds of equilibria and bifurcations without parameters in memristive circuits.  SIAM J. Appl. Math. 72(3), 877–896 (2012)
  * (25) Riaza, R.: Transcritical bifurcation without parameters in memristive circuits.  SIAM J. Appl. Math. 78(1), 395–417 (2018)
  * (26) Riaza, R.: Transcritical bifurcation without parameters in memristive circuits.  SIAM J. APPL. MATH. 78(1), 395–417 (2018)
  * (27) Semenov, V., Korneev, I., Arinushkin, P., Strelkova, G., Vadivasova, T., Anishchenko, V.: Numerical and experimental studies of attractors in memristor-based Chua’s oscillator with a line of equilibria. Noise-induced effects.  Eur. Phys. J. Special Topics 224(8), 1553–1561 (2015)
  * (28) Singh, J., Raj, B.: An accurate and generic window function for nonlinear memristor models.  Journal of Computational Electronics volume 18(2), 640–647 (2019)
  * (29) Tetzlaff, R. (ed.): Memristor and Memristive Systems.  Springer-Verlag New York (2014)
  * (30) Zhou, E., Fang, L., Yang, B.: A general method to describe forgetting effect of memristors.  Phys. Lett. A 383(11), 942–948 (2019)