# Generation of massively entangled bright states of light during harmonic
generation in resonant media

Sili Yi  Max Born Institute, Max-Born Straße 2A, D-12489 Berlin, Germany  Ihar
Babushkin  Max Born Institute, Max-Born Straße 2A, D-12489 Berlin, Germany
Cluster of Excellence PhoenixD (Photonics, Optics, and Engineering -
Innovation Across Disciplines), Welfengarten 1, 30167 Hannover, Germany
Institute for Quantum Optics, Leibniz University of Hannover, Germany  Olga
Smirnova  Max Born Institute, Max-Born Straße 2A, D-12489 Berlin, Germany
Department of Physics, Technical University Berlin, 10623 Berlin, Germany
Technion – Israel Institute of Technology, 3200003 Haifa, Israel  Misha Ivanov
Max Born Institute, Max-Born Straße 2A, D-12489 Berlin, Germany  Department of
Physics, Humboldt University, Newtonstraße 15, D-12489 Berlin, Germany
Blackett Laboratory, Imperial College London, South Kensington Campus, SW7 2AZ
London, United Kingdom  Technion – Israel Institute of Technology, 3200003
Haifa, Israel

###### Abstract

At the fundamental level, full description of light-matter interaction
requires quantum treatment of both matter and light. However, for standard
light sources generating intense laser pulses carrying quadrillions of photons
in a coherent state, classical description of light during intense laser-
matter interaction has been expected to be adequate. Here we show how
nonlinear optical response of matter can be controlled to generate dramatic
deviations from this standard picture, including generation of multiple
harmonics of the incident laser light entangled across many octaves. In
particular, non-trivial quantum states of harmonics are generated as soon as
one of the harmonics induces a transition between different laser-dressed
states of the material system. Such transitions generate an entangled light-
matter wavefunction, which emerges as the key condition for generating quantum
states of harmonics, sufficient even in the absence of a quantum driving field
or material correlations. In turn, entanglement of the material system with a
single harmonic generates and controls entanglement between different
harmonics. Hence, nonlinear media that are near-resonant with at least one of
the harmonics appear to be most attractive for controlled generation of
massively entangled quantum states of light. Our analysis opens remarkable
opportunities at the interface of attosecond physics and quantum optics, with
implications for quantum information science.

##  I Introduction

Quantum nature of light is central to our understanding of the micro-world,
with the quantum description of the photo-electric effect ushering the new era
in physics at the turn of the XX-th century. However, quantum description of
light is typically used in the regime where just a few, or few tens, of
photons interact with quantum matter. In contrast, when an intense incident
light carries quadrillions of photons, counting each one of them individually
is hardly expected to matter.

Recent observation [1] of the Schrödinger cat-type (so-called kitten) states
and related statistics [2] of just such an intense light upon nonlinear-
optical interaction with an atomic gas has upended this conventional wisdom.
Granted, the observed Schrödinger kitten states have only emerged upon
conditioning the observation of the transmitted driving field on detecting the
highly nonlinear optical response to it. Nevertheless, Ref. [2, 1], followed
by [3, 4, 5, 6, 7, 8, 9, 10], have opened an exciting possibility of using
high harmonic emission to generate bright quantum states of light by measuring
the transmitted laser field together with generated harmonics. Quantum
properties of the generated light can be further amplified if the incident
light interacts with a correlated quantum state of a material system, mapping
quantum correlations of a pre-excited massively entangled, e.g. Dicke-like,
state onto a quantum state of the generated light [11, 12, 13].

There is, of course, a different route to generating bright quantum states of
light via nonlinear-optical response: begin with an intense incident light
already in a quantum state [14, 15, 16]. As the quantum properties of the
incident light should affect the quantum properties of the generated
harmonics, the latter should emerge in nontrivial quantum states. The
technological breakthrough in generating the so-called bright squeezed vacuum
states of light [17, 18] with sufficient number of photons has enabled
generation of several harmonics of this light [14].

Here we introduce a completely different avenue. We start with a standard
laser pulse in a coherent state, and a quantum material system in a simple,
uncorrelated ground state. To generate harmonics in a non-trivial quantum
state, we rely only on the excitation of the material system by the generated
harmonic light, a process often ignored in the description of high harmonic
generation. Unsurprisingly, once such excitations are included, the light-
matter state becomes entangled. This entanglement is central to what follows.

First, perhaps unexpectedly, all the components in the nonlinear optical
response correlated to the new excited state of the material system become
entangled with each other. Second, and most importantly, control over the
excitations of the material system gives access to flexible control over the
generated entangled states of light across multiple octaves. Third, such
control over excitations in the material system needs to be neither
complicated nor sophisticated: it can make use of the ubiquitous Stark shifts
induced by the classical laser field, and be as straightforward as modulating
the intensity envelope of the driving classical field.

Our analysis shows that generation of non-trivial quantum states of harmonics
is not an exception: it should occur whenever the material system undergoes
non-adiabatic transitions between its laser-dressed states. It does not need
to rely on pre-exciting strongly correlated material states [12] or on using
non-classical incident light [14, 15, 16]. As for conditioning the
observations on a particular material response [2, 1, 13], this can also be
effectively removed by shaping the transverse modes of the generated harmonics
as discussed below.

We begin our analysis by first identifying the conditions that preclude
generation of nontrivial quantum states of multiple harmonics. This helps us
identify the key conditions for generating highly non-classical states of
harmonics, starting with a purely classical incident light and a simple
material system residing in an uncorrelated ground state. We then discuss one
route to implementing these conditions, which to a theorist’s eye appears
relatively straightforward.

##  II Prerequisites for generating quantum states of high harmonic light

The equation describing the interaction of light with a quantum system is
(atomic units are used throughout):

$\displaystyle i\frac{\partial}{\partial
t}\ket{\Psi}=\hat{\textbf{H}}\Psi=\left[\hat{\textbf{H}}_{A}+\hat{\textbf{x}}\hat{\textbf{F}}+\hat{\textbf{H}}_{F}\right]\ket{\Psi},$
(1)

where $\hat{H}_{A}$ is the Hamiltonian of the material system, $\hat{H}_{F}$
is the Hamiltonian of the quantum field, and
$\hat{\textbf{x}}\hat{\textbf{F}}$ describes their coupling in the dipole
approximation, with $\hat{\textbf{x}}$ representing the negative of the dipole
operator; $\hat{\textbf{F}}$ is the electric field operator acting on all
light modes, here characterized by their frequencies $\omega_{n}=n\omega$,
with $\omega$ the incident frequency and $n=2k+1$ describing its odd
harmonics. The wavefunction starts in the product state
$\ket{\Psi_{0}}=\ket{g}\ket{\alpha_{\omega}}$, with $\ket{g}$ the field-free
ground state of the material system and $\ket{\alpha_{\omega}}$ the initial
coherent state at the fundamental frequency $\omega$.

What happens if the full wavefunction $\ket{\Psi(t)}$ remains, throughout the
whole evolution, in the product state of the material system
$\ket{\psi_{g}(t)}$ and the state of light $\ket{\chi_{{}_{F}}(t)}$,
$\ket{\Psi(t)}=\ket{\psi_{g}(t)}\ket{\chi_{{}_{F}}(t)}$? Here the dressed
quantum state of the material system $\ket{\psi_{g}(t)}$ has evolved from the
field-free ground state and can include all states: ground, bound, and
continuum; its final overlap with the ground state might even be negligible.

Substituting the ansatz
$\ket{\Psi(t)}=\ket{\psi_{g}(t)}\ket{\chi_{{}_{F}}(t)}$ into the Eq.(1) yields
the following equation for the quantum light field,

$\displaystyle i\ket{\dot{\chi}_{{}_{F}}(t)}=$
$\displaystyle\left[\hat{H}_{F}+\textbf{x}_{gg}(t)\hat{\textbf{F}}+\Delta
E_{g}(t)\right]\ket{\chi_{{}_{F}}(t)},$ (2)

where $\textbf{x}_{gg}(t)=\bra{\psi_{g}(t)}\hat{\textbf{x}}\ket{\psi_{g}(t)}$
is the dipole induced in the material system; the time-dependent energy shift
in the material system $\Delta
E_{g}(t)=\bra{\psi_{g}(t)}\hat{H}_{A}-i\frac{\partial}{\partial
t}\ket{\psi_{g}(t)}$ has no influence on the quantum state of light
$|\chi_{F}\rangle$, which starts in the product of vacuum states for each
harmonic $n\omega$. Similar equation obtains for the material system.

For each $n\omega$, Eq.(2) describes a harmonic oscillator $\hat{H}_{F,n}$
starting in its ground state and driven by the time-dependent force
$\textbf{x}_{gg}(t)$. No matter how complex $\textbf{x}_{gg}(t)$ is, the
linear coupling $\textbf{x}_{gg}(t)\hat{\textbf{F}}$ can only shift the
initial vacuum state $\ket{\alpha_{n}=0}$ to a coherent state
$\ket{\alpha_{n}(t)}$ with a higher average photon number. The factorized
nature of the initial state $\ket{\chi_{{}_{F}}}$, which starts as a product
state of individual harmonics, each originally in its vacuum state, is
preserved, as well as the maximally classical nature of the generated
harmonics. What are the possible escape routes from this ”classical corner”,
apart from those listed in the introduction?

Mathematically, one route is to turn the linear Eq.(2) for the quantum state
of a generated harmonic $n$ into nonlinear: the material response
$\textbf{x}_{gg}(t)$ entering Eq.(2) for some harmonic $n$ should be sensitive
to the state $\ket{\chi_{{}_{F},_{n}}}$ of this harmonic. This removes the
curse of limiting the evolution of $\ket{\chi_{{}_{F},_{n}}}$ to a single
coherent state.

Physically, the dependence of the material response on the generated harmonics
implies that they can induce transitions in the material system and change its
nonlinear response. Thus, the quantum states of matter and the state of light
become correlated. The way to incorporate such correlation is to eschew the
factorized ansatz $\ket{\Psi(t)}=\ket{\psi_{g}(t)}\ket{\chi_{{}_{F}}(t)}$,
which brought us into the unwanted ”classical corner” in the first place,
explicitly providing for the possibility of light-matter entanglement. Quantum
correlations and feedback between the generated nonlinear optical response and
the state of the material system are the basis for generating nontrivial
quantum states of all harmonics.

To emphasize the lack of need for quantum light at the input, we shall replace
the field operator at the fundamental frequency $\omega$ with the classical
field $\textbf{F}_{\rm cl}(t)=\bra{\alpha_{1}(t)}\hat{\textbf{F}}_{\rm
1}\ket{\alpha_{1}(t)}$, approximating the full Hamiltonian as

$\displaystyle\hat{\textbf{H}}\simeq\hat{\textbf{H}}_{A}+\hat{\textbf{x}}\textbf{F}_{\rm
cl}(t)+\hat{\textbf{x}}\sum_{n\neq 1}\hat{\textbf{F}}_{n}+\sum_{n\neq
1}\hat{\textbf{H}}_{F,n},$ $\displaystyle\textbf{F}_{\rm
cl}(t)=\bra{\alpha_{1}(t)}\hat{\textbf{F}}_{\rm 1}\ket{\alpha_{1}(t)}.$ (3)

Such approximation is by no means necessary, as the general analysis is just
as straightforward as that developed below, but it allows us to bring to the
fore the key role of light-matter entanglement, which leads to the generation
of quantum light at the output without any quantum light at the input.

Figure 1:  Some schemes for generating harmonics in quantum states. (a) Use of
a cavity. Resonance between two laser-dressed states of a quantum system,
driven by a low-frequency laser field, and a generated harmonic entangles the
material system with the generated light. This also leads to massive
entanglement between all harmonics generated upon the excitation. The
resonance is enhanced by tuning the resonant cavity mode to the resonant
harmonic. (b) Instead of the cavity, one can rely on absorption of some
generated harmonic in a sufficiently dense medium. (c) Transverse shaping of
the resonant cavity mode separates in the far field the highly non-classical
components of the generated harmonics, correlated to the excited state, from
the (less interesting) harmonics correlated to the ground state.

How can we enhance the feedback of the generated harmonic light on the
material system that generates it? Given the relative weakness of generated
harmonics, the most straightforward route is to ensure a one-photon resonance
between some harmonic $n_{0}\omega$ and the transition between the ground
$\ket{g}$ and some excited $\ket{e}$ state of the material system.

Importantly, this one-photon resonance does not have to be present at the zero
fundamental field. On the contrary, we shall see below that it is beneficial
to have such resonance induced by the Stark shifts of the ground and excited
states.

To increase the light-matter coupling further, one can add a cavity resonant
with $n_{0}\omega$, as shown in Fig.1(a). Fig. 1(b) shows a complementary
option of an optically thick medium for some harmonics, which makes sure that
these harmonics are invested in building the correlation between the state of
the material system and the state of the generated light. Fig.1(c) shows that
one can additionally shape the transverse mode of the generated resonant
harmonic, separating in the far field all harmonics correlated to the excited
state of the driven material system from those correlated to its ground state.

We shall see that resonant absorption of a single harmonic affects the quantum
state of all other harmonics: single resonant transition is sufficient to
dramatically affect quantum properties of all generated light fields.

##  III General formalism

One-photon resonance with a particular harmonic means that there is also an
$n_{0}$-photon resonance with the fundamental field. Such multiphoton
resonances can be already included in the definition of the laser-driven
ground state $\ket{\psi_{g}(t)}$: we define $\ket{\psi_{g}(t)}$ as the
solution of the time-dependent Schrödinger equation (TDSE) in the classical
incident field $\textbf{F}_{\rm cl}(t)$, for the system starting in the ground
state,

$\displaystyle i\frac{\partial}{\partial
t}\ket{\psi_{g}(t)}=\left[\hat{H}_{A}+\hat{\textbf{x}}\textbf{F}_{\rm
cl}(t)\right]\ket{\psi_{g}(t)},$ $\displaystyle\ket{\psi_{g}(t_{i})}=\ket{g}.$
(4)

We will refer to this state as the ”classically dressed” ground state. One-
photon absorption of the harmonic $n_{0}\omega$ occurs between this state and
the ”classically dressed” excited state $\ket{\psi_{e}(t)}$ defined as the
solution of the standard semi-classical TDSE

$\displaystyle i\frac{\partial}{\partial
t}\ket{\psi_{e}(t)}=\left[\hat{H}_{A}+\hat{\textbf{x}}\textbf{F}_{\rm
cl}(t)\right]\ket{\psi_{e}(t)},$ $\displaystyle\ket{\psi_{e}(t_{i})}=\ket{e}.$
(5)

with the initial condition set to the excited field-free state $\ket{e}$. Note
that these dressed states fully incorporate the incident field, remain
orthogonal to each other, and offer a perfect basis to incorporate effects of
the harmonics on the material system.

We can now write the full light-matter wavefunction as:

$\displaystyle\ket{\Psi(t)}=\ket{\psi_{g}(t)}\ket{\chi_{g}(t)}+\sum_{e}\ket{\psi_{e}(t)}\ket{\chi_{e}(t)},$
(6)

where $\ket{\chi_{g}(t)},\ket{\chi_{e}(t)}$ are the states of light correlated
to the classically dressed states of the material system. Standard semi-
classical treatment of intense light-matter interaction of the quantum
material system with any classical field is limited to the single first term
in the above equation. The second term is responsible for all the new aspects.
To simplify the expressions that follow, we shall omit the summation over all
dressed excited states and only keep the resonant contribution; restoring the
full sum is straightforward.

Substituting Eq.(6) into the full Schrödinger equation yields coupled
equations for the light fields correlated to the two strongly driven states of
the material system,

$\displaystyle i\ket{\dot{\chi}_{g}}=$
$\displaystyle\hat{H}_{F}|\chi_{g}\rangle+\textbf{x}_{gg}(t)\hat{\textbf{F}}\ket{\chi_{g}}+\textbf{x}_{ge}(t)\hat{\textbf{F}}\ket{\chi_{e}},$
(7) $\displaystyle i\ket{\dot{\chi}_{e}}=$
$\displaystyle\hat{H}_{F}\ket{\chi_{e}}+\textbf{x}_{ee}(t)\hat{\textbf{F}}\ket{\chi_{e}}+\textbf{x}_{eg}(t)\hat{\textbf{F}}\ket{\chi_{g}},$
(8)

where $\hat{\textbf{F}}=\sum_{n>1}\hat{\textbf{F}}_{n}$ is the field operator
acting on all harmonics with orders $n\neq 1$,

$\displaystyle\textbf{x}_{gg}(t)=$
$\displaystyle\bra{\psi_{g}(t)}\hat{\textbf{x}}\ket{\psi_{g}(t)},$
$\displaystyle\textbf{x}_{ee}(t)=$
$\displaystyle\bra{\psi_{e}(t)}\hat{\textbf{x}}\ket{\psi_{e}(t)},$ (9)

are the laser-induced dipole in the dressed ground and excited states, and

$\displaystyle\textbf{x}_{ge}(t)=$
$\displaystyle\bra{\psi_{g}(t)}\hat{\textbf{x}}\ket{\psi_{e}(t)},$
$\displaystyle\textbf{x}_{eg}(t)=$
$\displaystyle\bra{\psi_{e}(t)}\hat{\textbf{x}}\ket{\psi_{g}(t)},$ (10)

describe the transition dipole matrix elements between the classically dressed
states. In the limit of weak driving, they coincide with the usual transition
dipoles between the field-free states and keep track of the transition
frequency $\omega_{eg}=E_{e}-E_{g}$,

$\displaystyle\textbf{x}_{eg}(t)=\bra{\psi_{e}(t)}\hat{\textbf{x}}\ket{\psi_{g}(t)}\rightarrow\textbf{x}_{eg}e^{i\omega_{eg}(t-t_{i})},$
(11)

with $t_{i}$ some initial moment of time. In stronger fields, these transition
dipoles can acquire non-trivial time-dependence, e.g. through the time-
dependence of the transition frequency $\omega_{eg}(t)=E_{e}(t)-E_{g}(t)$
which incorporates the Stark shifts of both states. We shall take advantage of
this ubiquitous effect below.

Addition of the couplings between different quantum states of light correlated
to different dressed states of matter is central to turning what is usually a
trivial evolution of the coherent states of all harmonics into an evolution
generating non-trivial quantum states of multiple entangled harmonics.

We now take advantage of the interaction picture, defining
$|\chi_{g,e}(t)\rangle=e^{-i\hat{H}_{F}t}|\phi_{g,e}(t)\rangle$. With this
substitution, we have

$\displaystyle i\ket{\dot{\phi}_{g}(t)}=$
$\displaystyle\textbf{x}_{gg}(t)\hat{\textbf{F}}_{I}(t)\ket{\phi_{g}}+\textbf{x}_{ge}(t)\hat{\textbf{F}}_{I}(t)\ket{\phi_{e}},$
$\displaystyle i\ket{\dot{\phi}_{e}(t)}=$
$\displaystyle\textbf{x}_{ee}(t)\hat{\textbf{F}}_{I}(t)\ket{\phi_{e}}+\textbf{x}_{eg}(t)\hat{\textbf{F}}_{I}(t)\ket{\phi_{g}},$
(12)

where the time dependence of the operator $\hat{\textbf{F}}_{I}(t)$ is given
by the transformation of the creation and annihilation operators for each
frequency $n\omega$, e.g. $\hat{a}_{{}_{I},_{n}}(t)=\hat{a}_{n}e^{-in\omega
t}$. The definition of $\hat{\textbf{F}}_{n}$ in terms of $\hat{a}$,
$\hat{a}^{\dagger}$ is fixed by setting the classical incident field
$\textbf{F}_{\rm cl}(t)=\langle\alpha_{1}(t)|\hat{\textbf{F}}_{\rm
1}|\alpha_{1}(t)\rangle$. If we set $F_{\rm cl}(t)\propto\cos\omega t\propto
e^{-i\omega t}+e^{+i\omega t}$, then in the interaction picture

$\displaystyle\hat{\textbf{F}}_{1,I}(t)=\boldsymbol{\mathcal{F}}_{1}f_{1}(t)\left[\hat{a}^{\dagger}_{1}e^{+i\omega
t}+\hat{a}_{1}e^{-i\omega t}\right],$ (13)

where $\boldsymbol{\mathcal{F}}_{1}$ is the vacuum field amplitude for the
fundamental, and we have added the envelope $f_{1}(t)$ to account for the
finite pulse length, assuming a sufficiently long pulse. We then use the same
convention for the field operators at all other relevant frequencies, with
$\boldsymbol{\mathcal{F}}_{n}$ the vacuum field amplitude for the harmonic
$n$.

The simplest route to move forward with the coupled equations Eq.(12) is to
use the time-dependent perturbation theory. Its zero-order equation describes
the light field correlated to the dressed ground state,

$\displaystyle
i\ket{\dot{\phi}_{g}^{(0)}(t)}=\textbf{x}_{gg}(t)\hat{\textbf{F}}_{I}(t)\ket{\phi_{g}^{(0)}(t)}$
(14)

with the initial condition being the vacuum state for all harmonics. Its
formal solution is
$\ket{\phi_{g}^{(0)}(t)}=\hat{U}^{(0)}_{g}(t,t_{i})\ket{\phi_{g}^{(0)}(t_{i})}$
where $\hat{U}^{(0)}_{g}(t,t_{i})$ is the propagator for the homogeneous
equation Eq.(14) with the driving function $\textbf{x}_{gg}(t)$. It is a
product of propagators for each harmonic,

$\displaystyle\hat{U}^{(0)}_{g}(t,t_{i})=\prod_{n}\hat{U}^{(0)}_{g,n}(t,t_{i}).$
(15)

At any time $t$, the quantum state of generated harmonics remains a product
state. Moreover, for each harmonic $n$ the propagator
$\hat{U}^{(0)}_{g,n}(t,t_{i})$ simply shifts the initial vacuum state
$|\alpha_{n}(t_{i})=0\rangle$ to the coherent state $|\alpha_{n}(t)\rangle$,
$|\alpha_{n}(t)\rangle=\hat{U}^{(0)}_{g,n}(t,t_{i})|\alpha_{n}(t_{i})\rangle$,
where

$\displaystyle\alpha_{n}(t)=-i\boldsymbol{\mathcal{F}}_{n}\int_{t_{i}}^{t}d\tau\textbf{x}_{gg,n}(\tau).$
(16)

Here $\textbf{x}_{gg,n}$ is the complex amplitude of the Fourier component of
the full dipole response $\textbf{x}_{gg}(t)$ at the frequency $n\omega$. For
convenience, we have retained the slow temporal dependence of the amplitudes
$x_{gg,n}=x_{gg,n}(\tau)$ on the envelope of the driving laser pulse, assuming
that the pulse includes sufficiently many oscillations under its envelope to
separate the two time-scales: the sub-cycle scale of the individual
oscillations and the multi-cycle scale of the envelope.

The zero-order solution is then substituted into the equation for the field
correlated to the dressed excited state, yielding

$\displaystyle|\phi_{e}^{(1)}(t)\rangle=$ (17) $\displaystyle=-i\int
dt^{\prime}\hat{U}^{(0)}_{e}(t,t^{\prime})\textbf{x}_{eg}(t^{\prime})\hat{\textbf{F}}_{I}(t^{\prime})\hat{U}^{(0)}_{g}(t^{\prime},t_{i})|\phi_{g}^{(0)}(t_{i})\rangle,$

where the propagator $\hat{U}^{(0)}_{e}(t,\tau)$ solves the homogeneous
equation similar to Eq.(14) but for the laser-induced dipole generated by the
material system in the dressed excited state, $\textbf{x}_{ee}(t)$.

Figure 2:  Generation of entangled harmonics correlated to the dressed excited
state of the material system. (a) Resonant excitation, which can happen at any
time $t^{\prime}$, alters the nonlinear material response and thus alters the
generated product state of harmonics correlated to a specific excitation time
$t^{\prime}$. The overall state $|\phi_{e}(t)\rangle$ is an integral, i.e. a
sum, over all excitation times $t^{\prime}$, and is therefore no longer a
product state. (b) The cartoon diagram describing the generation process, with
the states of the harmonic fields evolving under the propagator $U_{g}^{(0)}$
from $t_{i}$ until $t^{\prime}$ and $U_{e}^{(0)}$ from $t^{\prime}$ to $t$.
Since excitation can happen at any time, quantum mechanics demands summation
over $t^{\prime}$, resulting in the sum over the product states, i.e. in an
entangled state.

The physical meaning of this expression is simple and sketched in Figure 2.
From the initial time $t_{i}$ until some instant $t^{\prime}$ the quantum
fields at all harmonic frequencies are driven by the dipole response
$\textbf{x}_{gg}(t)$ generated by the material system in its dressed ground
state (which includes all dynamics induced by $\textbf{F}_{\rm cl}(t)$,
possibly even the full depletion of the field-free ground state.) At an
instant $t^{\prime}$ all quantum fields at all harmonics are ”lifted” by the
operator $\textbf{x}_{eg}(t^{\prime})\hat{\textbf{F}}_{I}(t^{\prime})$ to the
new electronic ”surface”, where they evolve from $t^{\prime}$ to $t$ driven by
the dipole $\textbf{x}_{ee}(t)$ correlated to the dressed excited state.

Given the near-resonant (or resonant) condition for the harmonic
$n_{0}\omega$, we shall use the resonant (rotating wave) approximation for the
one-photon transition driven by this harmonic. This means that the full
coupling term $\textbf{x}_{eg}(t^{\prime})\hat{\textbf{F}}_{I}(t^{\prime})$ is
approximated by its resonant contribution,

$\displaystyle\textbf{x}_{eg}(t^{\prime})\hat{\textbf{F}}_{I}(t^{\prime})\simeq\tilde{\textbf{x}}_{eg}e^{+i\int_{t_{i}}^{t^{\prime}}d\tau\omega_{eg}(\tau)}\boldsymbol{\mathcal{F}}_{n_{0}}\left[\hat{a}_{n_{0}}e^{-in_{0}\omega
t^{\prime}}\right]=$
$\displaystyle=\tilde{\textbf{x}}_{eg}e^{+i\int_{t_{i}}^{t^{\prime}}d\tau\Delta_{eg}(\tau)}\boldsymbol{\mathcal{F}}_{n_{0}}\hat{a}_{n_{0}},$
(18)

where the detuning of resonance $\Delta(t)=\omega_{eg}(t)-n_{0}\omega$
incorporates the Stark shifts of the classically dressed states, and
$\tilde{\mathbf{x}}_{eg}$ is the dipole transition matrix element at the
resonant frequency. This approximation assumes that electronic excitation
between the dressed states is accompanied by annihilating the resonant photon
$n_{0}\omega$. In principle, $\tilde{\mathbf{x}}_{eg}$ can retain a slow
temporal dependence on the pulse envelope, $\tilde{\mathbf{x}}_{eg}(\tau)$,
which incorporates additional distortions of the dressed states (beyond the
simple Stark shift) induced by the driving field $\omega$.

The integrand in Eq.(17) can now be written in a simple form. Indeed, the
initial quantum state of the harmonic light is a product of the vacuum states
for all harmonics, the propagator $\hat{U_{g}}^{(0)}(t^{\prime},t_{i})$ is a
product of linear zero-order propagators for each individual harmonic, and at
$t^{\prime}$ the quantum state of light remains a product state, just as in
the zero order discussed above. Moreover, just as in the zero order, for each
harmonic $n$ the propagator $\hat{U}^{(0)}_{g,n}(t^{\prime},t_{i})$ simply
shifts the vacuum state $|\alpha_{n}(t_{i})=0\rangle$ to the coherent state
$|\alpha_{n}(t^{\prime})\rangle$ with

$\displaystyle\alpha_{n}(t^{\prime})=-i\boldsymbol{\mathcal{F}}_{n}\int_{t_{i}}^{t^{\prime}}d\tau\textbf{x}_{gg,n}(\tau).$
(19)

Next, we apply the transition operator Eq.(18) in Eq.(17). It leaves all
$n\neq n_{0}$ undisturbed, while for $n=n_{0}$ it yields
$\hat{a}_{n_{0}}|\alpha_{n_{0}}(t^{\prime},t_{i})\rangle=\alpha_{n_{0}}(t^{\prime})|\alpha_{n_{0}}(t^{\prime})\rangle$.
Finally, the action of the propagator $\hat{U}^{(0)}_{e,n}(t,t^{\prime})$ on
the coherent state $|\alpha_{n}(t^{\prime})\rangle$ yields another coherent
state at the moment $t$,
$\hat{U}^{(0)}_{e,n}(t,t^{\prime})|\alpha_{n}(t^{\prime})\rangle=|\alpha_{n}(t,t^{\prime})\rangle$,
where

$\displaystyle\alpha_{n}(t,t^{\prime})=$ (20)
$\displaystyle=-i\boldsymbol{\mathcal{F}}_{n}\left[\int_{t_{i}}^{t^{\prime}}d\tau\textbf{x}_{gg,n}(\tau)+\int_{t^{\prime}}^{t}d\tau\textbf{x}_{ee,n}(\tau)\right].$

This yields the final result in the first-order perturbation theory with
respect to the quantum feedback between the generated light and the material
system,

$\displaystyle|\phi_{e}^{(1)}(t)\rangle=$ (21) $\displaystyle=-i\int
dt^{\prime}\tilde{\textbf{x}}_{eg}(t^{\prime})e^{+i\int_{t_{i}}^{t^{\prime}}d\tau\Delta_{eg}(\tau)}\boldsymbol{\mathcal{F}}_{n_{0}}\alpha_{n_{0}}(t^{\prime})\prod_{n>1}|\alpha_{n}(t,t^{\prime})\rangle.$

In general, the above result represents a massively entangled quantum state of
light: indeed, the integral is a sum of the product states
$\prod_{n>1}|\alpha_{n}(t,t^{\prime})\rangle$, with $t^{\prime}$ a continuous
summation index.

Below we discuss several specific scenarios that realize different entangled
states, demonstrating routes to control the resulting outcome. Note that all
generated states emerge from coherent states characterized by their indices
$\alpha_{n}$ Eq.(19). Thus, in what follows, we will characterize the
generated states using the average photon numbers

$\displaystyle
N_{n}(T)=|\boldsymbol{\mathcal{F}}_{n}\int_{t_{i}}^{T}d\tau\textbf{x}_{gg,n}(\tau)|^{2}$
(22)

that the coherent states at frequencies $\omega_{n}=n\omega$ would have
acquired at the end of the laser pulse $T$, had the material system stayed in
the dressed ground state throughout the whole laser pulse. We will use
$N_{3}=N_{3}(T)$ for the third harmonic as the overall reference number. For
other harmonics, the key numbers are then $N_{n}/N_{3}$, which are controlled
by the parameters of the laser pulse and the cavity. For example, if $n=5$ is
resonant with the cavity mode, the choice of the laser frequency, intensity,
and the quality of the cavity can allow one to vary $N_{5}/N_{3}$ in a broad
range. We shall also use $N_{n}(t)$ (for $t<T$) to denote the average number
of photons in a coherent state $|\alpha_{n}\rangle$ that would have been
generated by the material system in its dressed ground state by the time $t$
(given by the same equation Eq.(22) with $t$ in the upper limit of the
integral.)

##  IV Results and Discussion

To appreciate the entangled nature of the generated harmonics state, and the
possibility to control it, we begin by considering the so-called Freeman
resonances [19, 20], where the conventional Stark shift proportional to the
intensity of the classical IR driver brings the states in and out of
resonance, as shown in Fig.3(a). The integral is now accumulated during the
temporal windows near the resonant times $t_{k}$ and can be computed using the
standard stationary phase method, yielding at the end of the laser pulse $T$

$\displaystyle|\phi_{e}^{(1)}(T)\rangle=\sum_{k}C_{k}e^{+i\int_{t_{i}}^{t_{k}}d\tau\Delta_{eg}(\tau)}\prod_{n>1}|\alpha_{n}(T,t_{k})\rangle,$
$\displaystyle C_{k}=\sqrt{\frac{-2\pi
i}{\dot{\Delta}_{eg}(t_{k})}}\tilde{\textbf{x}}_{eg}\boldsymbol{\mathcal{F}}_{n_{0}}\alpha_{n_{0}}(t_{k})\
.$ (23)

Here $\dot{\Delta}_{eg}(t_{k})$ is the time-derivative of the resonance
detuning at the moment of resonance $t_{k}$. The physical meaning of $C_{k}$
is nothing but the amplitude of one-photon excitation by the resonant harmonic
$n_{0}$ during the resonance window $t_{k}$. Here, this amplitude is computed
using first order perturbation theory and is proportional to
$\boldsymbol{\mathcal{F}}_{n_{0}}\alpha_{n_{0}}(t_{k})$, which is simply the
classical electric field amplitude
$\mathbf{F}_{n_{0}}(t_{k})=\boldsymbol{\mathcal{F}}_{n_{0}}\alpha_{n_{0}}(t_{k})$
of the resonant harmonic $n_{0}\omega$ at the moment $t_{k}$. Naturally, the
transition amplitude $C_{k}$ also includes the transition dipole matrix
element $\tilde{\mathbf{x}}_{eg}$. This matrix element may change from its
field-free value because the states $|\psi_{g}\rangle$ and $|\psi_{e}\rangle$
of the material system are dressed by the laser field. Therefore, it depends
on the pulse envelope, and this dependence can be incorporated as
$\tilde{\mathbf{x}}_{eg}=\tilde{\mathbf{x}}_{eg}(t_{k})$. The phase factors
$e^{+i\int_{t_{i}}^{t_{k}}d\tau\Delta_{eg}(\tau)}$ control the relative phases
between the terms in the sum. Overall, the amplitudes $C_{k}$ describe the
non-adiabatic Landau-Zener-Dykhne transitions between the dressed states
driven through the resonance, with the phase factors
$e^{+i\int_{t_{i}}^{t_{k}}d\tau\Delta_{eg}(\tau)}$ responsible for the
Stueckelberg oscillations (named so after Baron Ernst Carl Gerlach
Stueckelberg von Breidenbach zu Breidenstein und Melsbach) in the final
generated state, as the coherent states are not orthogonal. We stress again
that both the number of the terms in the sum and the relative phases between
them are controlled by shaping the intensity envelope of the fundamental laser
pulse.

Figure 3:  Control of resonances and of the generated quantum states of light.
(a) Stark shift follows the intensity envelope of the driving laser field and
brings the ground and excited states in resonance twice during the pulse.
(b,c) Wigner function for one of the harmonics, here $n=3$, upon artificially
removing its entanglement with others. For illustration, we set resonances at
$t_{1}$ and $t_{2}$ to occur when the average number of the harmonic photons
$N_{3}(t)$ generated by the laser-dressed ground state are $N_{3}(t_{1})=25$,
$N_{3}(t_{2})=50$ (b) and $N_{3}(t_{1})=100$, $N_{3}(t_{2})=225$ (c). (d,e)
Intensity of the $X$-quadratures of the entangled state of two harmonics, here
$n=3$ and $n=5$, $|\phi_{e}(X_{3},X_{5})|^{2}$. The parameters are set to
$N_{3}(t_{1})=N_{5}(t_{1})=25$, $N_{5}(t_{2})=N_{3}(t_{2})=50$ in (d) and
$N_{3}(t_{1})=N_{5}(t_{1})=100$, $N_{5}(t_{2})=N_{3}(t_{2})=225$ in (e).

Examples of quantum states of light that can be created this way are shown in
Fig.3, for two resonant windows and hence two terms in the sum. For simplicity
of visual presentation, we assume that only two harmonics, the third (n=3) and
the fifth (n=5), are efficiently generated, so that the entangled wavefunction
contains only two terms in each product state,

$\displaystyle\ket{\phi_{e}^{(1)}(T)}=$
$\displaystyle=C_{1}\ket{\alpha_{3}(T,t_{1})}\ket{\alpha_{5}(T,t_{1})}+C_{2}e^{i\Phi}\ket{\alpha_{3}(T,t_{2})}\ket{\alpha_{5}(T,t_{2})},$
$\displaystyle\Phi=\int_{t_{1}}^{t_{2}}d\tau\Delta_{eg}(\tau).$ (24)

This constraint is only used for simplicity of visualizing the overall quantum
state of light. Note that $\alpha_{n}$ and $\Phi$ are fully controlled by the
parameters of the driving laser pulse and the parameters of the cavity.

To give an impression of the underlying structure of a single frequency mode,
say $n=3$, we first artificially set
$\alpha_{5}(T,t_{1})=\alpha_{5}(T,t_{2})$, which allows us to disentangle the
third harmonic from the fifth and plot its Wigner function. The state is,
clearly, a Schrödinger cat
$\ket{\alpha_{3}(T,t_{1})}+e^{+i\Phi}\ket{\alpha_{3}(T,t_{2})}$. Fig.3 shows
the Wigner functions of this cat state for the case of $\Phi=\pi$ and
$N_{3}(t_{1})=25$, $N_{3}(t_{2})=50$ in panel (b) and $N_{3}(t_{1})=100$,
$N_{3}(t_{2})=225$ in panel (c). Note that both the relative phase $\Phi$ and
the photon numbers are fully controlled by the time delay between the
resonances and the Stark shifts induced by the driving laser field. The
characteristic non-classical features of the Schrödinger cat state are
apparent in Figs. 3(b,c).

We next move to the example where both the third and the fifth harmonics are
generated with similar efficiency. This could be achieved, for example, by
tuning the cavity resonance close to the $n=5$ harmonic and staying within the
perturbative regime with respect to the driving laser field, so that higher
order harmonics are negligible. The 2D plot of the resulting state as a
function of the x-quadratures $X_{3}$ and $X_{5}$ of both harmonics is shown
in Fig.3(d), for $N_{3}(t_{1})=N_{5}(t_{1})=25$ and
$N_{3}(t_{2})=N_{5}(t_{2})=50$, with equal amplitudes $|C_{1}|=|C_{2}|$ and
the relative phase $\Phi=\pi$ between the two terms. Panel (e) shows the same
for $N_{3}(t_{1})=N_{5}(t_{1})=100$ and $N_{3}(t_{2})=N_{5}(t_{2})=225$.

Figure 4:  Examples of the harmonic states correlated to the excited state in
case of resonance at all times. (a,b) Wigner functions for one of the
harmonics, $n=3$, upon removing its entanglement with others, for the average
number of photons $N_{3}(T)=10$ (a) and $N_{3}(T)=100$ (b). (c,d) Intensities
of the entangled wavefunctions in X-quadratures for $N_{3}(T)=N_{5}(T)=10$ (c)
and $N_{3}(t)=N_{5}(t)=100$ (d).

Our next example is shown in Fig.4. Here we consider exact resonance between
the dressed states $|\psi_{e,g}\rangle$ at all times and perform the
integration in Eq.(III) directly. For illustrative purposes we set the
nonlinear response in the dressed excited state to be negligible compared to
that in the dressed ground state, which is typical for a loosely bound excited
state in a not-so-strong field regime, before the onset of recollision-driven
harmonics. We set the resonant transition dipole matrix element
$\mathbf{x}_{eg}$ to a constant. To be able to plot the resulting quantum
states, we also neglect all harmonics apart from $n=3,5$.

As in the previous case of the Freeman resonances, to give an impression of
the underlying structure of a single frequency mode, say $n=3$, we first
artificially suppress the generation of all harmonics except for one and plot
its Wigner function for $N_{3}=|\alpha_{3}(T)|^{2}=10$ and
$N_{3}=|\alpha_{3}(T)|^{2}=100$, see Fig.4(a,b). We see that individual
harmonics become progressively more squeezed as the number of photons in the
harmonic increases. Thus, in contrast to normal expectations and standard
scenarios, the quantum properties of the light state are becoming stronger as
the number of photons grows. This growing squeezing is already interesting,
but the overall quantum state is richer: once we include the generation of
both harmonics $n=3,5$, they become entangled.

The developing quantum correlation is already apparent from the state
intensity $|\phi_{e}(X_{3},X_{5})|^{2}$ shown as a function of the
X-quadrature for both harmonics, which becomes both squeezed and rotated – a
tell-tale sign of a correlated state. To characterize entanglement, we use
both the Schmidt number and the entanglement entropy. Both are obtained from
the Schmidt decomposition of $\phi_{e}(X_{3},X_{5})$ as a sum over the
separable states [21, 22]:
$\phi_{e}(X_{3},X_{5})=\sum_{i}\lambda_{i}\phi_{e,i}(X_{3})\phi_{e,i}(X_{5})$,
where $i=1,2,3,\ldots$, $\phi_{i}(X_{3})$, $\phi_{i}(X_{5})$ represents a set
of single-harmonic wavefunctions, and $\lambda_{i}>0$ are the corresponding
eigenvalues of the Schmidt decomposition, normalized so that
$\sum_{i}\lambda_{i}^{2}=1$. Entanglement entropy $S$ is defined then as
$S=-\sum_{i}\lambda_{i}^{2}\log_{2}\lambda_{i}^{2}$, whereas the Schmidt
number $D$ is defined [22] as $D=1/{\sum{\lambda_{i}^{4}}}$. If the state is
separable, we have $S=0$ and $D=1$, whereas for nonseparable states $S>0$,
$D>1$. Both quantities show how entangled the states are. For the generated
states, these measures are shown in Figure 5, clearly demonstrating strong
photon-number entanglement developed between harmonics. Importantly, the
degree of entanglement and non-classicality of harmonics grows with increasing
the average number of photons per harmonic, a very exciting trend.

Figure 5: Entanglement entropy $S$ and the Schmidt number $D$ as functions of
the average number of photons $N_{3}$ and $N_{5}$. (a) $N_{3}$ is changed,
keeping constant the ratio $N_{5}/N_{3}$. (b) The ratio $N_{5}/N_{3}$ is
modified, keeping $N_{3}$ constant. (c) Both $N_{3}$ and $N_{5}$ are modified
(only entanglement entropy is shown).

##  V Conclusions and Outlook

Our work leads to several important conclusions. We find that non-trivial
quantum states of harmonic light are generated as soon as quantum correlations
develop between the generated light and the generating material system. A
sufficient condition to developing these correlations is rather general and
quite straightforward: (i) the system should undergo non-adiabatic excitations
between laser dressed states and (ii) the different laser-dressed states
should exhibit different non-linear responses. If this is the case,
excitations of the material system into different dressed states would
generate not only an entangled light-matter wavefunction, but also an
entanglement between the different harmonics correlated to the excited state
of the material system. Entanglement between harmonics emerges whenever the
material excitation occurs at different times throughout the laser pulse.

On the one hand, a specific mechanism of resonant excitations between
different laser-dressed states can be quite general and does not need to be
limited to the absorption of a generated harmonic. On the other hand, resonant
transitions in a strongly driven system, with large time-dependent Stark
shifts, are virtually inevitable, making the process considered here rather
typical.

High harmonics correlated to real excitations of matter are often ignored due
to their poor phase matching properties. Indeed, excitation of loosely bound
Rydberg states implies the so-called long trajectories in recollision-driven
[23, 24] high harmonics. The contribution of these long trajectories to the
macroscopic response is generally weak due to the large intensity-dependent
phase associated with these trajectories [25, 26]. However, recent experiments
and theory [27] show that this need not be the case, and that harmonics
generated via excited states can even dominate the far-field spectrum.
Nevertheless, to improve phase matching it might be beneficial to look at
intermediate field regimes before the onset of highly non-perturbative
recollision-type harmonic generation mechanisms, reducing the complexities of
phase-matching large intensity-dependent harmonic phases accumulated during
electron motion in the continuum.

The presence of the cavity is highly beneficial but not strictly necessary:
the essential aspect is significant excitation probability. In this context,
two-dimensional solid state materials, where excitons provide strong
contribution to the optical response, appear an interesting candidate for
generating non-classical light.

An interesting opportunity is offered by alkali atoms in a cavity. This may
look like a suggestion to turn the clock back a few decades, returning to a
well-studied system. However, we suggest to add a new twist: instead of
driving the atomic medium with a resonant field, here we suggest to drive it
with a low-frequency field, tuning the laser frequency into a multiphoton
resonance with the transition between the ground and the first excited states
of an alkali atom.

Control over the strengths and times of resonant excitations becomes a way of
controlling the generated entanglement. We have illustrated this point by
using the so-called Freeman resonances [19], which are both ubiquitous and
unavoidable during intense light-matter interaction. These resonances are
induced by shifting an excited state into a resonance with the ground state
due to the ponderomotive shift $U_{p}=F^{2}/4\omega^{2}$ ($F$ is the field
strength), associated with the wiggling energy of a loosely bound electron
driven by the incident laser field. The resonance is multi-photon for the
driving field and one-photon for one of the generated harmonics. Crucially,
already a single resonant harmonic can trigger generation of non-trivial
quantum states for all other harmonics: all it takes is for the resonant
harmonic to invest its photons in material excitation, benefiting all others.

Even in the simplest case of a continuous resonance occurring throughout the
whole laser pulse, which corresponds to a perturbative regime where the
ponderomotive Stark shift is negligible, all harmonics correlated to the
system left in the excited state are generated in non-trivial quantum states.

What is more, in contrast to normal expectations, the quantum properties of
the generated harmonics become stronger as the number of photons in each
harmonic grows.

We stress that the excited state is by no means obligated to efficiently
generate harmonics on its own. In fact, in the cases considered here, the
nonlinear components of the excited state polarization were negligible
compared to those for the dressed ground state, befitting a loosely bound
state in a moderate-intensity regime (before the onset of recollision-driven
harmonics.)

Since Stark shifts, which control resonances in low-frequency laser fields,
depend on the laser intensity, shaping the intensity envelope of the driving
field shapes the quantum nature of the generated harmonic fields. Here we have
only considered the case of two resonances achieved at the front and at the
back of the pulse. Of course, one can control the number of resonances by
shaping the laser pulse envelope in time.

At the same time, shaping the laser pulse envelope in space will generate
areas in the interaction region where harmonics are generated in different
entangled states. These will propagate into the far field. Thus, spatial
shaping of the generating pulse and/or harmonic modes (e.g. via transverse
modes of the resonant cavity) opens routes for generating highly nontrivial
quantum states of light where the spatial harmonic parameters such as orbital
angular momenta are also entangled.

The use of spatially structured light to shape the far-field macroscopic
response, or the design of cavities to shape the transverse modes of the
generated harmonics are just some examples of the additional opportunities for
controlling the modes of the entangled states of high harmonics. In
particular, employing a nontrivial transverse cavity mode will lead to far-
field separation of harmonics correlated to the excited state and generated in
nontrivial quantum states from harmonics correlated to the system left in the
ground state and generated in a product of coherent states. On the one hand,
this allows one to avoid the conditioning of the generated quantum states on
measuring the material system in the excited state. On the other hand, such
conditioning might also be possible if the excited material system is a
molecule that dissociates: the measurement of light in coincidence with
measuring excited atomic fragments, while extremely challenging, is not
impossible.

Applications of nonclassical, massively entangled states of high harmonics can
be very broad [10]. These include imaging of material properties, where
nonclassical and entangled states of light provide quantum advantage [28],
resulting in, for instance, enhanced absorption probability and resolution, or
selectivity with respect to different quantum pathways in wave-mixing
spectroscopy. Furthermore, entangled combs are known to be a resource for
quantum information processing, such as quantum communication and cryptography
[29]. Whereas typical combs are produced by electro-optic modulations of
cavities and are limited by GHz bandwidth, the approach presented here
promises to expand it to the PHz frequency range.

This discussion demonstrates truly remarkable opportunities arising when
combining the fool toolkit of attosecond physics with the field of quantum
optics, with exciting implications in different fields ranging from ultrafast
and precision spectroscopy to quantum information science.

###### Acknowledgements.

M.I. acknowledges the support of the Horizon 2020 Framework Programme, award
number 899794, and Limati SFB 1777 “Light-matter interaction at interfaces”
project, award number 441234705. O.S. acknowledges support of the ERC grant
’ULISSES’, award number 101054696. I. B. acknowledges support from the
Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy within the
Cluster of Excellence PhoenixD (EXC 2122, Project No. 390833453). M.I.
acknowledges extraordinary hospitality at the Technion, especially after
07.10.2023.

## References

  * Lewenstein  _et al._ [2021] M. Lewenstein, M. Ciappina, E. Pisanty, J. Rivera-Dean, P. Stammer, T. Lamprou, and P. Tzallas, Generation of optical Schrödinger cat states in intense laser-matter interactions, Nat. Phys. 17, 1104 (2021).
  * Tsatrafyllis  _et al._ [2017] N. Tsatrafyllis, I. Kominis, I. Gonoskov, and P. Tzallas, High-order harmonics measured by the photon statistics of the infrared driving-field exiting the atomic medium, Nat. Commun. 8, 15170 (2017).
  * Wang  _et al._ [2022] Z. Wang, Z. Bao, Y. Wu, Y. Li, W. Cai, W. Wang, Y. Ma, T. Cai, X. Han, J. Wang, _et al._ , A flying schrödinger’s cat in multipartite entangled states, Sci. Adv. 8, eabn1778 (2022).
  * Rivera-Dean  _et al._ [2022a] J. Rivera-Dean, T. Lamprou, E. Pisanty, P. Stammer, A. F. Ordóñez, A. S. Maxwell, M. F. Ciappina, M. Lewenstein, and P. Tzallas, Strong laser fields and their power to generate controllable high-photon-number coherent-state superpositions, Phys. Rev. A 105, 033714 (2022a).
  * Rivera-Dean  _et al._ [2022b] J. Rivera-Dean, P. Stammer, A. S. Maxwell, T. Lamprou, P. Tzallas, M. Lewenstein, and M. F. Ciappina, Light-matter entanglement after above-threshold ionization processes in atoms, Phys. Rev. A 106, 063705 (2022b).
  * Maxwell  _et al._ [2022] A. S. Maxwell, L. B. Madsen, and M. Lewenstein, Entanglement of orbital angular momentum in non-sequential double ionization, Nat. Commun. 13, 4706 (2022).
  * Stammer [2022] P. Stammer, Theory of entanglement and measurement in high-order harmonic generation, Phys. Rev. A 106, L050402 (2022).
  * Stammer  _et al._ [2022] P. Stammer, J. Rivera-Dean, T. Lamprou, E. Pisanty, M. F. Ciappina, P. Tzallas, and M. Lewenstein, High photon number entangled states and coherent state superposition from the extreme ultraviolet to the far infrared, Phys. Rev. Lett. 128, 123603 (2022).
  * Stammer  _et al._ [2023] P. Stammer, J. Rivera-Dean, A. Maxwell, T. Lamprou, A. Ordóñez, M. F. Ciappina, P. Tzallas, and M. Lewenstein, Quantum electrodynamics of intense laser-matter interactions: A tool for quantum state engineering, PRX Quantum 4, 010201 (2023).
  * Bhattacharya  _et al._ [2023] U. Bhattacharya, T. Lamprou, A. S. Maxwell, A. Ordonez, E. Pisanty, J. Rivera-Dean, P. Stammer, M. F. Ciappina, M. Lewenstein, and P. Tzallas, Strong–laser–field physics, non–classical light states and quantum information science, Rep. Prog. Phys. 86 (2023).
  * Gorlach  _et al._ [2020] A. Gorlach, O. Neufeld, N. Rivera, O. Cohen, and I. Kaminer, The quantum-optical nature of high harmonic generation, Nat. Commun. 11, 4598 (2020).
  * Pizzi  _et al._ [2023] A. Pizzi, A. Gorlach, N. Rivera, A. Nunnenkamp, and I. Kaminer, Light emission from strongly driven many-body systems, Nat. Phys. 19, 551 (2023).
  * Tzallas [2023] P. Tzallas, Quantum correlated atoms in intense laser fields, Nat. Phys. 19, 472 (2023).
  * Spasibko  _et al._ [2017] K. Y. Spasibko, D. A. Kopylov, V. L. Krutyanskiy, T. V. Murzina, G. Leuchs, and M. V. Chekhova, Multiphoton effects enhanced due to ultrafast photon-number fluctuations, Phys. Rev. Lett. 119, 223603 (2017).
  * Even Tzur  _et al._ [2023] M. Even Tzur, M. Birk, A. Gorlach, M. Krüger, I. Kaminer, and O. Cohen, Photon-statistics force in ultrafast electron dynamics, Nat. Photon. 17, 501 (2023).
  * Gorlach  _et al._ [2023] A. Gorlach, M. E. Tzur, M. Birk, M. Krüger, N. Rivera, O. Cohen, and I. Kaminer, High-harmonic generation driven by quantum light, Nat. Phys. 19, 1689–1696 (2023).
  * Iskhakov  _et al._ [2009] T. Iskhakov, M. V. Chekhova, and G. Leuchs, Generation and direct detection of broadband mesoscopic polarization-squeezed vacuum, Phys. Rev. Lett. 102, 183602 (2009).
  * Iskhakov  _et al._ [2012] T. S. Iskhakov, A. Pérez, K. Y. Spasibko, M. Chekhova, and G. Leuchs, Superbunched bright squeezed vacuum state, Opt. Lett. 37, 1919 (2012).
  * Freeman  _et al._ [1987] R. Freeman, P. Bucksbaum, H. Milchberg, S. Darack, D. Schumacher, and M. Geusic, Above-threshold ionization with subpicosecond laser pulses, Phys. Rev. Lett. 59, 1092 (1987).
  * Gibson  _et al._ [1992] G. Gibson, R. Freeman, and T. McIlrath, Verification of the dominant role of resonant enhancement in short-pulse multiphoton ionization, Phys. Rev. Lett. 69, 1904 (1992).
  * Law  _et al._ [2000] C. K. Law, I. A. Walmsley, and J. H. Eberly, Continuous frequency entanglement: Effective finite hilbert space and entropy control, Phys. Rev. Lett. 84, 5304 (2000).
  * Grobe  _et al._ [1994] R. Grobe, K. Rzazewski, and J. Eberly, Measure of electron-electron correlation in atomic physics, J. Phys. B 27, L503 (1994).
  * Schafer  _et al._ [1993] K. Schafer, B. Yang, L. DiMauro, and K. Kulander, Above threshold ionization beyond the high harmonic cutoff, Phys. Rev. Lett. 70, 1599 (1993).
  * Corkum [1993] P. B. Corkum, Plasma perspective on strong field multiphoton ionization, Phys. Rev. Lett. 71, 1994 (1993).
  * Lewenstein  _et al._ [1995] M. Lewenstein, P. Salieres, and A. L’huillier, Phase of the atomic polarization in high-order harmonic generation, Phys. Rev. A 52, 4747 (1995).
  * Gaarde  _et al._ [1999] M. Gaarde, F. Salin, E. Constant, P. Balcou, K. Schafer, K. Kulander, and A. L’Huillier, Spatiotemporal separation of high harmonic radiation into two quantum path components, Phys. Rev. A 59, 1367 (1999).
  * Mayer  _et al._ [2022] N. Mayer, S. Beaulieu, A. Jimenez-Galan, S. Patchkovskii, O. Kornilov, D. Descamps, S. Petit, O. Smirnova, Y. Mairesse, and M. Ivanov, Role of spin-orbit coupling in high-order harmonic generation revealed by supercycle rydberg trajectories, Phys. Rev. Lett. 129, 173202 (2022).
  * Dorfman  _et al._ [2016] K. E. Dorfman, F. Schlawin, and S. Mukamel, Nonlinear optical signals and spectroscopy with quantum light, Rev. Mod. Phys. 88, 045008 (2016).
  * Kues  _et al._ [2019] M. Kues, C. Reimer, J. M. Lukens, W. J. Munro, A. M. Weiner, D. J. Moss, and R. Morandotti, Quantum optical microcombs, Nat. Photon. 13, 170 (2019).


# On discrete-time polynomial dynamical systems on hypergraphs

Shaoxuan Cui1,4, Guofeng Zhang2, Hildeberto Jardón-Kojakhmetov1 and Ming Cao3
1 S. Cui, and H. Jardón-Kojakhmetov are with the Bernoulli Institute for
Mathematics, Computer Science and Artificial Intelligence, University of
Groningen, Groningen, 9747 AG Netherlands {s.cui,
<EMAIL_ADDRESS>G. Zhang is with the Department of Applied
Mathematics, The Hong Kong Polytechnic University, Kowloon 999077, Hong Kong,
China and The Hong Kong Polytechnic University Shenzhen Research Institute,
Shenzhen, Guang Dong 518057, China<EMAIL_ADDRESS>M. Cao is with
the Engineering and Technology institute Groningen, University of Groningen,
Groningen, 9747 AG Netherlands<EMAIL_ADDRESS>S.Cui was supported by China
Scholarship Council.

###### Abstract

This paper studies the stability of discrete-time polynomial dynamical systems
on hypergraphs by utilizing the Perron–Frobenius theorem for nonnegative
tensors with respect to the tensors’ Z-eigenvalues and Z-eigenvectors. First
of all, for a multilinear polynomial system on a uniform hypergraph, we study
the stability of the origin of the corresponding systems. Afterward, we extend
our results to non-homogeneous polynomial systems on non-uniform hypergraphs.
We confirm that the local stability of any discrete-time polynomial system is
in general dominated by pairwise terms. In particular, given the origin is
locally stable, we construct a conservative (but explicit) region of
attraction from the system parameters. Finally, we validate our results via
some numerical examples.

###### Index Terms:

Hypergraphs, Higher-order interactions, Polynomial systems, Z-eigenvalues,
Perron-Frobenius Theorem, Stability

##  I Introduction

Many real complex systems such as those originating from epidemics [1, 2, 3],
biology [4, 5, 6, 7], and engineering [8, 9, 10, 11] are usually modeled as
polynomial systems and studied from a network perspective. However, a
conventional network only captures pairwise information flow, as characterized
by edges of 2-node links of networks, and may lose some higher-order
information [12, 13, 14, 15, 16]. Nowadays, there is abundant evidence that a
hypergraph, a generalization of graphs where each edge contains multiple
($\geq 2$) nodes, is a more powerful modeling tool because a hypergraph can
capture higher-order information [13, 14, 15, 16].

For example, with the introduction of the concept of hypergraphs, an epidemics
model on a hypergraph takes multi-body interactions into account [17, 18, 19]
and is more suitable to describe the process of information diffusion.
Similarly, by introducing higher-order terms into the classical Lotka-Volterra
model, a higher-order Lotka-Volterra model is proposed by [20] and studied
mathematically by [21]. This higher-order Lotka-Volterra model takes indirect
higher-order effects among species into consideration, which better reflects
the reality according to many empirical studies [22, 23]. All such higher-
order systems [18, 19, 20, 21] are, in fact, non-homogeneous polynomial
systems [9].

Although the majority of real systems evolve in continuous time, a discrete-
time system, as an approximation of its continuous-time counterpart, still
attracts much attention. In control engineering, a controller may need to be
implemented in digital hardware [24]. In epidemics, data may be gathered daily
and thus it is easier to deal with the real data when using a discrete-time
model [1, 2, 3]. All these factors motivate us to study a general discrete-
time polynomial system. So far, many related works on discrete-time polynomial
systems [25, 26] rely on the Sum-Of-Squares (SOS) decomposition, which was
first introduced by [27]. The main technique is to find a Lyapunov function in
the form of a sum of squares. Generally, the problem of checking the
nonnegativity of a function is NP-hard [25]. Some other researches rely on the
technique of Kronecker products and Linear Matrix Inequalities [28, 29].
Another research line uses semi tensor product to study the problem [30]. All
these mentioned results may be useful for designing a controller to stabilize
a polynomial system. However, such results are limited when one wants to know
some characteristics of the open-loop system, such as the domain of attraction
of a locally stable equilibrium. Knowing the behavior of an open-loop
autonomous system is often of great significance if we deal with real systems
like epidemics or species populations, where it is very difficult to implement
a controller.

From a modeling perspective, it is well-known that a graph can be captured by
its adjacency matrix. For the analysis of polynomial systems on a graph, the
Perron–Frobenius theorem [31] is a fundamental tool. The Perron–Frobenius
theorem shows that an irreducible nonnegative matrix has a positive eigenvalue
associated with a positive eigenvector. With the development of tensor
algebra, tensor versions of the classical Perron–Frobenius theorem are
provided by, e.g., [32, 33, 34, 35] concerning both H-eigenvalues and
Z-eigenvalues. So, since a hypergraph can be represented by an adjacency
tensor [15, 36], it is reasonable to consider a tensor version of the
Perron–Frobenius theorem as a potential tool. In [37], the author has already
used the tensor’s Z-eigenvalues and tensor decomposition to study the
stability of a discrete-time homogeneous polynomial system. However, we know
that a discrete-time system is often obtained by applying discretization to a
continuous-time system. After discretization, the discrete-time system rarely
remains homogeneous. Thus, the result can be restricted for potential
applications. For continuous-time polynomial systems on hypergraphs, by using
the Perron–Frobenius theorem, the global stability (e.g. of the origin) can be
checked by the tensors’ H-eigenvalues and H-eigenvectors [38]. It is fair to
expect a similar result for a discrete-time system.

_The contributions of the paper are summarized as follows:_ first, regarding a
discrete-time homogeneous polynomial system on a uniform hypergraph, we
achieve a similar result as in [37] with a different proof by using a tensor
version of the Perron–Frobenius theorem (see Section III). Next, concerning
discrete-time non-homogeneous polynomial systems on non-uniform hypergraphs,
under a mild condition, and after knowing an equilibrium is stable, we show
that a domain of attraction even though conservative, can be directly
calculated from the system’s parameters. Lastly, we design some feedback
controllers that can reduce or enlarge the aforementioned domain of
attraction. To validate our analytical results, we also provide some numerical
examples.

_Notation:_ $\mathbb{R}$ ($\mathbb{R}_{++}$) denotes the set of (positive)
real numbers. The superscript in e.g. $\mathbb{R}^{n}$ denotes the dimension
of the space. For a matrix $M\in\mathbb{R}^{n\times r}$ and a vector
$a\in\mathbb{R}^{n}$, $M_{ij}$ and $a_{i}$ denote the element in the $i$th row
and $j$th column and the $i$th entry, respectively. Given a square matrix
$M\in\mathbb{R}^{n\times n}$, $\rho(M)$ denotes the spectral radius of $M$,
which is the largest absolute value of the eigenvalues of $M$. The notation
$|M|$ denotes the matrix whose entry $|M|_{ij}$ is the absolute value of
$M_{ij}$. For any two vectors $a,b\in\mathbb{R}^{n}$, $a>(<)b$ represents that
$a_{i}>(<)b_{i}$, for all $i=1,\ldots,n$; and $a\geq(\leq)b$ means that
$a_{i}\geq(\leq)b_{i}$, for all $i=1,\ldots,n$. These component-wise
comparisons are also used for matrices or tensors with the same dimensions.
The vector $\mathbf{1}$ ($\mathbf{0}$) represents the column vector or matrix
of all ones (zeros) with appropriate dimensions. The previous notations have
straightforward extensions to tensors. In the following section, we introduce
the preliminaries and useful notation regarding tensors.

##  II Preliminaries on tensors and hypergraphs

A tensor $T\in\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{k}}$ is a
multidimensional array, where the order is the positive integer $k$ and each
dimension $n_{i}$, $i=1,\cdots,k$ is a mode of the tensor. A tensor is cubical
if every mode has the same size, that is $T\in\mathbb{R}^{n\times
n\times\cdots\times n}$. We further write a $k$-th order $n$-dimensional
cubical tensor as $T\in\mathbb{R}^{n\times n\times\cdots\times
n}=\mathbb{R}^{[k,n]}$. A cubical tensor $T$ is called supersymmetric if
$T_{j_{1}j_{2}\ldots j_{k}}$ is invariant under any permutation of the
indices.

For the rest of the paper, a tensor always refers to a cubical tensor unless
specified otherwise. We then consider the following notation: for a tensor
$A\in\mathbb{R}^{[m,n]}$, $Ax^{m-1}$ is a vector, whose $i$-th components are

$\displaystyle\left(Ax^{m-1}\right)_{i}$
$\displaystyle=\sum_{i_{2},\ldots,i_{m}=1}^{n}A_{i,i_{2}\cdots
i_{m}}x_{i_{2}}\cdots x_{i_{m}}.$

For a tensor $A\in\mathbb{R}^{[m,n]}$, consider the following equations:

$\begin{split}Ax^{m-1}&=\lambda x,\;\\\ x^{\top}x&=1;\end{split}$ (1)

where if there is a real number $\lambda$ and a nonzero real vector $x$ that
are solutions of (1), then $\lambda$ is called a Z-eigenvalue of $A$ and $x$
is the Z-eigenvector of $A$ associated with $\lambda$ [39, 32]. Throughout
this paper, the words eigenvalue and eigenvector as well as Z-eigenvalue and
Z-eigenvector are used interchangeably. It is worth mentioning that there are
some other kinds of definitions of the eigenvalues and eigenvectors of a
tensor, e.g. H-eigenvalues [32, 39] and U-eigenvalues [40]. However, in this
paper, we will always refer to Z-eigenvalues and Z-eigenvectors as defined
above.

The tensor $\mathcal{I}_{z}=\left(e_{i_{1}\ldots
i_{m}}\right)\in\mathbb{R}^{[m,n]}$ is the Z-identity tensor [41] defined as a
nonnegative tensor such that

$\mathcal{I}_{z}x^{m-1}=x$

for all $x\in\mathbb{R}^{n}$ with $x^{\top}x=1$. There is no odd-order
Z-identity tensor, only an even-order Z-identity tensor. The choice of even-
order Z-identity tensor is not unique. When $m$ is even, let $B_{i_{1}\cdots
i_{m}}=\delta_{i_{1}i_{2}}\delta_{i_{3}i_{4}}\cdots\delta_{i_{m-1}i_{m}}$,
where $\delta$ is a standard Kronecker delta. Then, $B$ is a Z-identity tensor
[39].

A tensor
$\mathcal{C}=\left(\mathcal{C}_{{i_{1}}\ldots{i_{m}}}\right)\in\mathbb{R}^{[m,n]}$
is called reducible if there exists a nonempty proper index subset
$I\subset\\{1,\ldots,n\\}$ such that

$\mathcal{C}_{i_{1}\cdots i_{m}}=0\quad\forall i_{1}\in I,\quad\forall
i_{2},\ldots,i_{m}\notin I.$

If $\mathcal{C}$ is not reducible, then we call $\mathcal{C}$ irreducible. A
tensor with all non-negative entries is called a non-negative tensor.

We now recall the Perron–Frobenius Theorem for irreducible nonnegative tensors
with respect to the Z-eigenvalue:

######  Lemma 1 (Theorems 4.5 and 4.6 [32])

If $A\in\mathbb{R}_{++}^{[m,n]}$, then there exists a Z-eigenvalue
$\lambda_{0}\geqslant 0$ and a nonnegative Z-eigenvector $x_{0}\neq\mathbf{0}$
of $A$ such that $Ax_{0}^{m-1}=\lambda_{0}x_{0}$. We call $\lambda_{0},x_{0}$
the Perron-Z-eigenvalue and -eigenvector respectively, and refer to
$(\lambda_{0},x_{0})$ as a Perron-Z-eigenpair. If, in addition,
$A\in\mathbb{R}_{++}^{[m,n]}$ is irreducible, then the Perron-Z-eigenpair
$\left(\lambda_{0},x_{0}\right)$ further satisfies

  * 1.

The eigenvalue $\lambda_{0}$ is positive.

  * 2.

The eigenvector $x_{0}$ is positive, that is, all components of $x_{0}$ are
positive.

It is worthwhile to mention that the Perron-Z-eigenpairs
$\left(\lambda_{0},x_{0}\right)$ are generally not unique.

Next, we summarize some definitions regarding hypergraphs.

The concept of a hypergraph is defined in [36]. Here, we use a set of tensors
to denote the weights of a hypergraph. A weighted and directed hypergraph is a
triplet $\mathbf{H}=(\mathcal{V},\mathcal{E},A)$. The set $\mathcal{V}$
denotes a set of vertices and $\mathcal{E}=\\{E_{1},E_{2},\cdots,E_{n}\\}$ is
the set of hyperedges. A hyperedge is defined as an ordered pair
$E=(\mathcal{X},\mathcal{Y})$ of disjoint subsets of vertices; $\mathcal{X}$
is the tail of $E$ and $\mathcal{Y}$ is the head. As a special case, a
weighted and undirected hypergraph is a triplet
$\mathbf{H}=(\mathcal{V},\mathcal{E},A)$, where $\mathcal{E}$ is a finite
collection of non-empty subsets of $\mathcal{V}$ [15]. If all hyperedges of
the hypergraph contain the same number of (tails, heads) nodes, then the
hypergraph is uniform. For details, see [12] for undirected uniform
hypergraphs and see [42] for directed uniform hypergraphs. For modeling
purposes, one directed hyperedge usually denotes the joint influence of a
group of agents on one agent. Thus, it suffices to deal with hyperedges with
one single tail and we assume that each hyperedge has only one tail but
possibly multiple ($\geq 1$) heads. This setting is similar to [42] and has
the advantage that a directed uniform hypergraph can be represented by an
adjacency tensor. Generally, an undirected uniform hypergraph can be
represented by a supersymmetric adjacency tensor. For a non-uniform
hypergraph, we now use the set of tensors $A=\\{A_{2},A_{3},\cdots\\}$ to
collect the weights of all hyperedges, where $A_{2}=[A_{ij}]$ denotes the
weights of all second-order hyperedges, $A_{3}=[A_{ijk}]$ denotes the weights
of all third-order hyperedges, and so on. For instance, $A_{ijkl}$ denotes the
weight of the hyperedge where $i$ is the tail and $j,k,l$ are the heads. For
simplicity, in this paper, we also use the weight (for example, $A_{\bullet}$)
to denote the corresponding hyperedge. If all hyperedges only have one tail
and one head, then the network is a standard directed and weighted graph.

For convenience, throughout this paper, we define a multi-index notation
$I=i_{2},\cdots,i_{k}$, where $k$ is the order of the associated tensor.

##  III Discrete-time polynomial systems on uniform hypergraphs

Here, we consider a discrete-time polynomial system on a uniform hypergraph of
$n$ nodes. The dynamical system is given by

$x^{+}=Ax^{k-1},$ (2)

where $A\in\mathbb{R}^{[m,n]}$, $x\in\mathbb{R}^{n}$ is the state variable
$x(t)$, the time instant $t$ is omitted without ambiguity and the shorthand
$x^{+}$ denotes $x(t+1)$. Component-wise, (2) reads as

$x^{+}_{i}=\sum_{i_{2},\ldots,i_{k}=1}^{n}A_{i,I}x_{i_{2}}\cdots x_{i_{k}}.$
(3)

Now, we are ready to discuss the stability of the origin.

###### Theorem 1

The origin is always an equilibrium of (2). Moreover, if the tensor $A$ is
irreducible, the origin is asymptotically stable with a domain of attraction
$\max_{j}\frac{|x_{j}(0)|}{\delta_{j}}<(\frac{1}{\lambda})^{\frac{1}{k-2}}$,
where $(\lambda,\delta)$ is a Perron-Z-eigenpair of $|A|$.

###### Proof:

The first claim is straightforward. Next, we show the stability result. We
define the Lyapunov function $V=\max_{j}\frac{|x_{j}|}{\delta_{j}}$. Since the
Perron-eigenvector is strictly positive, it holds that $V>0$ for any
$x\neq\mathbf{0}$, and $V=0$ if and only if $x=\mathbf{0}$. Furthermore, for
any $i$, we have

$x_{i}\leq\max_{j}\left(\frac{x_{j}}{\delta_{j}}\right)\delta_{i}=V\delta_{i}.$
(4)

We suppose that at time $t$,
$q=\operatorname*{argmax}_{j}\left(\frac{|x_{j}|}{\delta_{j}}\right)$ and that
at time $t+1$,
$p=\operatorname*{argmax}_{j}\left(\frac{|x^{+}_{j}|}{\delta_{j}}\right)$,

Then, we get

$\begin{split}V^{+}-V&=\frac{1}{\delta_{p}}(|\sum_{i_{2},\ldots,i_{k}=1}^{n}A_{p,I}x_{i_{2}}\cdots
x_{i_{k}}|)-V\\\
&\leq\frac{1}{\delta_{p}}(\sum_{i_{2},\ldots,i_{k}=1}^{n}|A_{p,I}||x_{i_{2}}|\cdots|x_{i_{k}}|)-V\\\
&\leq\frac{1}{\delta_{p}}(\sum_{i_{2},\ldots,i_{k}=1}^{n}|A_{p,I}|V^{k-1}\delta_{i_{2}}\cdots\delta_{i_{k}})-V\\\
&=\frac{1}{\delta_{p}}(V^{k-1}\lambda\delta_{p})-V=V(V^{k-2}\lambda-1),\end{split}$
(5)

where we used (4) to obtain the third line. Once there is a time $t$ such that
$\max_{j}\frac{|x_{j}(t)|}{\delta_{j}}<(\frac{1}{\lambda})^{\frac{1}{k-2}}$,
then $V^{+}<V$ and thus $\max_{j}\frac{|x_{j}(t+1)|}{\delta_{j}}\leq
V<(\frac{1}{\lambda})^{\frac{1}{k-2}}$. By the same analysis, $V(t+2)\leq
V(t+1)$. Thus, by induction, we conclude that the origin is asymptotically
stable with a domain of attraction
$\max_{j}\frac{|x_{j}(0)|}{\delta_{j}}<(\frac{1}{\lambda})^{\frac{1}{k-2}}$. ∎

###### Remark 1

Notice that the domain of attraction in Theorem 1 is derived through a
Lyapunov function and thus may be conservative. In [37], if $A$ is
orthogonally decomposable, a further explicit solution of system (2) is
obtained. The result of Theorem 1 is similar to [37, Proposition 3].

##  IV Discrete-time polynomial systems on non-uniform hypergraphs

To extend the result of the previous section, in this section, we describe the
behavior of discrete-time polynomial systems on non-uniform hypergraphs.

Usually, systems on a non-uniform hypergraph incorporate tensors of different
orders. The dynamical system we consider in this section is given by

$x^{+}=A_{k-1}x^{k-1}+A_{k-2}x^{k-2}+\cdots+A_{1}x,$ (6)

where all the tensors $A_{k-1},\cdots,A_{1}$ are tensors of dimension $n$ and
order $k,k-1,k-2,\cdots$ respectively. To be general, some of the non-leading
order tensor $A_{i}$s with $i\neq k-1$ could be absent and thus be an all-zero
tensor.

###### Theorem 2

The origin is always an equilibrium of (6). Suppose that all the tensors
$|A_{k-1}|,\cdots,|A_{1}|$ are irreducible (some of them could be absent, but
all the remaining tensors must be irreducible) and have a common Perron-Z-
eigenvector $\delta$ and let the corresponding Perron-Z-eigenvalue be
$\lambda(|A_{k-1}|),\cdots,\lambda(|A_{1}|)$. If $A_{1}$ is not an all-zero
matrix and all eigenvalues of $A_{1}$ satisfy $|\tilde{\lambda}(A_{1})|<1$,
then the origin is locally asymptotically stable. If $\lambda(|A_{1}|)<1$,
then the origin is asymptotically stable with a domain of attraction
$\max_{j}\frac{|x_{j}(0)|}{\delta_{j}}<y_{+}$, where $y_{+}$ is the unique
positive solution of the equation $\sum_{i=1}^{k-1}\lambda(|A_{i}|)y^{i-1}=1$,
where $y$ is an unknown scalar.

###### Proof:

The local stability argument when $A_{1}$ is not an all-zero matrix is
straightforward. We see that the Jacobian at the origin is $J=A_{1}$. Thus, if
all eigenvalues of $A_{1}$ satisfy $|\tilde{\lambda}(A_{1})|<1$, the origin is
locally asymptotically stable.

Define the Lyapunov function
$V=\max_{j}\left(\frac{|x_{j}|}{\delta_{j}}\right)$. Consider that
$m=\operatorname*{argmax}_{j}\left(\frac{|x_{j}|}{\delta_{j}}\right)$ and
there may be a swap from time to time. We suppose at time $t$,
$q=\operatorname*{argmax}_{j}\left(\frac{|x_{j}|}{\delta_{j}}\right)$ and at
time $t+1$,
$p=\operatorname*{argmax}_{j}\left(\frac{|x^{+}_{j}|}{\delta_{j}}\right)$,

Then, we get

$\begin{split}&V^{+}-V=\frac{1}{\delta_{p}}\Big{(}\Big{|}\sum_{i_{2},\ldots,i_{k}=1}^{n}(A_{k-1})_{p,I}x_{i_{2}}\cdots
x_{i_{k}}+\\\ &+\sum_{i_{3},\ldots,i_{k}=1}^{n}(A_{k-2})_{p,I}x_{i_{3}}\cdots
x_{i_{k}}+\cdots\Big{|}\Big{)}-V\\\
&\leq\frac{1}{\delta_{p}}(\sum_{i_{2},\ldots,i_{k}=1}^{n}|(A_{k-1})_{p,I}||x_{i_{2}}|\cdots|x_{i_{k}}|+\cdots)-V\\\
&\leq\frac{1}{\delta_{p}}(V^{k-1}\lambda(|A_{k-1}|)\delta_{p}+V^{k-2}\lambda(|A_{k-2}|)\delta_{p}+\cdots)-V\\\
&=V(\sum_{i=1}^{k-1}\lambda(|A_{i}|)V^{i-1}-1);\end{split}$ (7)

where $\sum_{i=1}^{k-1}\lambda(A_{i})V^{i-1}-1$ yields the equation
$\sum_{i=1}^{k-1}\lambda(A_{i})y^{i-1}=1$. We now show that there is a unique
positive solution for the equation
$f(y)=\sum_{i=1}^{k-1}\lambda(|A_{i}|)y^{i-1}-1=0$. We notice that
$f(0)=\lambda(A_{1})-1<0$ and $\lim_{y\rightarrow\infty}f(y)>0$, since all
$\lambda(|A_{i}|)$ are positive. Note that $f(x)$ is a continuous function.
Then, by the intermediate value theorem, there must be at least one positive
solution. Furthermore, we see that
$f^{{}^{\prime}}(y)=\sum_{i=2}^{k-1}(i-1)\lambda(|A_{i}|)y^{i-2}>0$ for any
positive $y$. This ensures the uniqueness. We further see that if there is a
unique positive solution $y_{+}$ for the equation
$\sum_{i=1}^{k-1}\lambda(A_{i})y^{i-1}=1$, then for the inequality
$\sum_{i=1}^{k-1}\lambda(A_{i})y^{i-1}-1<0$, the solution set is
$\\{x|x<y_{+}\\}$. By a similar induction argument as in Theorem 1, one
completes the proof. ∎

###### Remark 2

Let $C$ be an irreducible non-negative matrix and $D$ an arbitrary matrix.
According to Wielandt’s theorem [43], if $|D|\leq C$, then any eigenvalue
$\tilde{\lambda}(D)$ satisfies $|\tilde{\lambda}(D)|\leq\rho(C)$. Now, let
$C=|A_{1}|$ and $D=A_{1}$. Wielandt’s theorem then shows that
$|\tilde{\lambda}(A_{1})|\leq\lambda(|A_{1}|)<1$, since in the matrix case the
Perron-Z-eigenvalue is the spectral radius. Thus, the condition
$\tilde{\lambda}(A_{1})\leq\lambda(|A_{1}|)<1$ tells us that the pairwise
interaction must be stable. Theorem 2 further provides a conservative region
of attraction. From the definition of irreducibility, the condition that all
the tensors $|A_{k-1}|,\cdots,|A_{1}|$ are irreducible is equivalent to the
condition that all the tensors $A_{k-1},\cdots,A_{1}$ are irreducible. As a
special case, $A_{1}$ could be absent, then $\lambda(|A_{1}|)<1$ holds
automatically.

A clear disadvantage of Theorem 2 is that the condition that all tensors are
to be associated with a common Perron-Z-eigenvector is too strict, and cannot
be always satisfied. We proceed to relax such a restriction.

###### Theorem 3

Consider system (6). Suppose that all the tensors $A_{k-1},\cdots,A_{2},A_{1}$
are irreducible (some of them could be absent, but all the remaining tensors
must be irreducible). Then, if $\max_{p}(\sum_{l}|A_{1}|_{p,l})<1$, the origin
is asymptotically stable with a domain of attraction
$\max_{j}{|x_{j}(0)|}<\min_{p}y_{p+}$, where $y_{p+}$ is the unique positive
scalar solution of
$\sum_{i=1}^{k-1}(\sum_{i_{2},\cdots,i_{k}=1}^{n}|(A_{i})_{p,I}|y^{i-2})=1$
corresponding to the index $p$.

###### Proof:

We know that $|A_{i}|$ is an irreducible nonnegative tensor. Thus,
$|A_{i}|\delta^{i}>\mathbf{0}$ for any positive vector $\delta>\mathbf{0}$.
Since $\delta$ is an arbitrary positive vector, we let
$\delta=(\epsilon,\epsilon,\cdots,\epsilon)^{\top}$ for an arbitary
$\epsilon>0$. Then, $|A_{i}|\delta^{i}>\mathbf{0}$ still holds.

Similar to the proof of Theorem 2, we define
$V=\max_{j}\left(\frac{|x_{j}|}{\delta_{j}}\right)=\frac{\max_{j}|x_{j}|}{\epsilon}$.
We suppose that at time $t$,
$q=\operatorname*{argmax}_{j}\left(\frac{|x_{j}|}{\delta_{j}}\right)$ and at
time $t+1$,
$p=\operatorname*{argmax}_{j}\left(\frac{|x^{+}_{j}|}{\delta_{j}}\right)$.
Then, we get

$\begin{split}V^{+}-V&=\frac{1}{\epsilon}\Big{(}\Big{|}\sum_{i_{2},\cdots,i_{k}}(A_{k-1})_{p,I}x_{i_{2}}x_{i_{3}}\cdots
x_{i_{k}}\\\ &+\sum_{i_{3},\cdots,i_{k}}(A_{k-2})_{p,I}x_{i_{3}}\cdots
x_{i_{k}}+\cdots\Big{|}\Big{)}-V\\\
&\leq\frac{1}{\epsilon}\left(\sum_{i_{2},\cdots,i_{k}}|(A_{k-1})_{p,I}||x_{i_{2}}||x_{i_{3}}|\cdots|x_{i_{k}}|\right.\\\
&\left.+\sum_{i_{3},\cdots,i_{k}}|(A_{k-2})_{p,I}||x_{i_{3}}|\cdots|x_{i_{k}}|+\cdots\right)-V\\\
&\leq\frac{1}{\epsilon}\left(\sum_{i_{2},\cdots,i_{k}}|(A_{k-1})_{p,I}|V^{k-1}\epsilon^{k-1}+\cdots\right.\\\
&\left.+\sum_{i_{2}}|(A_{2})_{p,I}|V\epsilon\right)-V\\\
&=V\left(\sum_{i=1}^{k-1}\left(\sum_{I}|(A_{i})_{p,I}|(V\epsilon)^{i-2}\right)-1\right).\end{split}$
(8)

Using similar arguments as in the proof of Theorem 2, we know that there is a
unique positive solution $y_{p+}$ for the equation
$\sum_{i=1}^{k-1}(\sum_{i_{2},\cdots,i_{k}=1}^{n}|(A_{i})_{p,I}|y^{i-2})=1$ if
$\max_{p}(\sum_{l}|A_{1}|_{p,l})<1$. Making $V^{+}-V<0$ requires (i):
$V\epsilon=\max_{j}|x_{j}|<y_{p+}$ and
(ii):$\sum_{i2=1}^{n}|(A_{k-2})_{p,I}|-1<0$. Notice that the index $p$ may
change from time to time. However, once $\max_{j}{|x_{j}(0)|}<y_{p+},\forall
p=1,\cdots,n$ and $\sum_{l}|A_{1}|_{p,l}<1,\forall p=1,\cdots,n$, then the
conditions (i) and (ii) are satisfied. This completes the proof. ∎

###### Remark 3

Theorem 3 is very useful for potential application. It just requires that each
tensor $A_{i}$ is an irreducible tensor. Moreover, $y_{p+}$ is only related to
the absolute value of the entries of each tensor $A_{i}$. If we know all
tensors, it is easy to compute the domain of attraction. Furthermore, from the
Gershgorin circle theorem, the condition $\max_{p}(\sum_{l}|A_{1}|_{p,l})<1$
guarantees that $|\tilde{\lambda}(A_{1})|\leq\lambda(|A_{1}|)<1$ and thus the
origin must be locally stable.

Next, let us consider the following quadratic system:

$x^{+}=A_{2}x^{2}+A_{1}x.$ (9)

Using Theorem 3 we have a more concrete result:

###### Corollary 1

The origin is always an equilibrium of (9). Suppose that both tensors
$A_{2},A_{1}$ are irreducible. Then, if $\max_{p}(\sum_{l}|A_{1}|_{p,l})<1$,
the origin is asymptotically stable with a domain of attraction
$\max_{j}{|x_{j}(0)|}<\min_{p}\frac{1-\sum_{l}|A_{1}|_{p,l}}{\sum_{l,m}|A_{2}|_{p,lm}}$.

Similarly, we can look at the cubic system:

$x^{+}=A_{3}x^{3}+A_{2}x^{2}+A_{1}x.$ (10)

###### Corollary 2

The origin is always an equilibrium of (10). Suppose that all the tensors
$A_{3},A_{2},A_{1}$ are irreducible. Then, if
$\max_{p}(\sum_{l}|A_{1}|_{p,l})<1$, the origin is asymptotically stable
within a domain of attraction
$\max_{j}{|x_{j}(0)|}<\min_{p}\frac{-C_{3}+\sqrt{C_{2}^{2}-4C_{3}(C_{1}-1)}}{2C_{3}}$,
where
$C_{3}=\sum_{l,m,q}|A_{3}|_{p,lmq},\,C_{2}=\sum_{l,m}|A_{2}|_{p,lm},\,C_{1}=\sum_{l}|A_{1}|_{p,l}$.

We don’t go further to higher orders since there do not exist closed-form
solutions for general fifth or higher-order polynomial equations. Instead, we
rather deal with polynomial systems on non-uniform hypergraphs with constant
terms. That is, the dynamics we now consider are given by

$x^{+}=A_{k-1}x^{k-1}+A_{k-2}x^{k-2}+\cdots+A_{1}x+b,$ (11)

where $b$ is a constant vector.

Suppose (11) has a fixed point $x=a$, which leads to
$x^{+}-x=0=A_{k-1}a^{k-1}+A_{k-2}a^{k-2}+\cdots+A_{1}a+b-a$. We calculate the
error dynamic with the shift $y=x-a$ in the following form:

$\begin{split}y^{+}&=x^{+}-a=A_{k-1}(y+a)^{k-1}+\cdots+A_{1}(y+a)+b-a\\\
&={A}_{k-1}y^{k-1}+\cdots+{A}_{1}y+\cdots\\\
&+\underbrace{A_{k-1}a^{k-1}+A_{k-2}a^{k-2}+\cdots+A_{1}a+b}_{=a}-a\\\
&=\tilde{A}_{k-1}y^{k-1}+\cdots+\tilde{A}_{1}y.\end{split}$ (12)

We can see that the constant term vanishes. Therefore, we can factorize the
remaining terms and thus we get a system in the form of
$y^{+}=\tilde{A}_{k-1}y^{k-1}+\cdots+\tilde{A}_{1}y$, where the tensors
$\tilde{A}_{i}$ depend on $a$ and the $A_{i}$’s.

###### Example 1

Consider a simple system in the form of (11): $x^{+}=x^{2}-2x+2$. We see that
$x^{+}-x=x^{2}-3x+2$. Thus, $x=1$ and $x=2$ are both fixed points of
$x^{+}=x^{2}-2x+2$. Let $y=x-1$, this yields $y^{+}=y^{2}$, which is in the
form of (12). Let $z=x-2$, this yields that $z^{+}=z^{2}+2z$, which is also in
the form of (12).

###### Remark 4

Once we rewrite the system (11) as (12), all results regarding the origin of
the system (12) apply to system (11) at $x=a$.

##  V Feedback control strategies

In this section, we propose some feedback control strategies for systems (6)
to manipulate the conservative domain of attraction in Theorems 2. We consider
the close loop system with a feedback controller
$g(u)=s\tilde{\mathcal{I}}x^{l-1}$ with $l$ an even number:

$x^{+}=A_{k-1}x^{k-1}+A_{k-2}x^{k-2}+\cdots+A_{1}x+g(u),$ (13)

where $\tilde{\mathcal{I}}$ is a tensor of diemnsion $n$ and order $l>2$ we
intentionally designed and $s$ is a real scalar.

We design
$\tilde{\mathcal{I}}_{i_{1},i_{2},\cdots,i_{n}}=({\mathcal{I}}_{z})_{i_{1},i_{2},\cdots,i_{n}}$
if $({\mathcal{I}}_{z})_{i_{1},i_{2},\cdots,i_{n}}=0$; otherwise,
$\tilde{\mathcal{I}}_{i_{1},i_{2},\cdots,i_{n}}=\operatorname*{\mathsf{sgn}}((A_{l})_{i_{1},i_{2},\cdots,i_{n}})|({\mathcal{I}}_{z})_{i_{1},i_{2},\cdots,i_{n}}|$
if $({\mathcal{I}}_{z})_{i_{1},i_{2},\cdots,i_{n}}\neq 0$. We assume that
$\max_{i_{1},i_{2},\cdots,i_{n}}|A_{l}|_{i_{1},i_{2},\cdots,i_{n}}>|s|$ if
$s<0$. By this construction, we have that
$(|A_{l}+s\tilde{\mathcal{I}}|)_{i_{1},i_{2},\cdots,i_{n}}=(|A_{l}|)_{i_{1},i_{2},\cdots,i_{n}}+s({\mathcal{I}}_{z})_{i_{1},i_{2},\cdots,i_{n}}$
such that $|A_{l}+s\tilde{\mathcal{I}}|=|A|+s{\mathcal{I}}_{z}$. Then,
$(|A_{l}|+s{\mathcal{I}}_{z})\delta^{k-1}=(\lambda(|A_{l}|)+s)\delta$, where
$(\lambda,\delta)$ is a Perron-Z-eigenpair of $|A_{l}|$. This means
$\lambda(|A_{l}|)$ in Theorem 2 is substituted into $\lambda(|A_{l}|)+s$.
Consider the equation $f(y)=\sum_{i=1}^{k-1}\lambda(|A_{i}|)y^{i-1}-1=0$. If
we enlarge (decrease) $\lambda(|A_{l}|)$ with a $s>0$ ($s<0$), since
$f(0)=\lambda(|A_{1}|)-1$ is unchanged but $f^{{}^{\prime}}(y)$ increases
(decreases) as $\lambda(|A_{l}|)$ increases (decreases), thus $y^{+}$
decreases (increases). In this way, we can manipulate the conservative domain
of attraction. For the even-order multilinear system (2), since the system (2)
is a special case of (6), the control strategy is still applicable when we
choose $g(u)=s\tilde{\mathcal{I}}x^{k-1}$.

##  VI Numerical examples and further discussions

Consider system (9) with $A_{1}=\left[\begin{array}[]{cc}0.1&0.1\\\
0.1&0.1\end{array}\right]$ and $A_{2}$ being a cubical tensor of order $3$
dimension $2$, whose entries are all one except
$(\tilde{A}_{2})_{112}=(\tilde{A}_{2})_{212}=(\tilde{A}_{2})_{121}=(\tilde{A}_{2})_{221}=0.5$.
Component-wise, the dynamic reads as

$x_{i}^{+}=0.1x_{1}+0.1x_{2}+x_{1}^{2}+x_{1}x_{2}+x_{2}^{2},\quad i=1,2.$ (14)

We can then calculate, from Corollary 1, that the domain of attraction is
$\max_{j}{|x_{j}(0)|}<\min_{p}\frac{1-\sum_{l}|A_{1}|_{p,l}}{\sum_{l,m}|A_{2}|_{p,lm}}=\frac{4}{15}$.
Then, we use Matlab to scatter the domain of attraction for the origin. If
from an initial condition, the solution of (9) finally converges to the
origin, we mark the point of the initial condition as an ’o’. Otherwise, it is
out of the domain of attraction for the origin, we mark it as an ’x’. The
result is shown in figure 1. We can see that the conservative region of
attraction according to the Theorem 3 is indeed a region of attraction of the
origin.

Figure 1: Points with ’o’ are within the domain of attraction for the origin
and points with ’x’ are out of the domain of attraction. The square area
(without the boundary) denotes the conservative region of attraction according
to the Theorem 3 and Corollary 1.

Now, let $\tilde{A}_{2}$ be a cubical tensor with the order $3$ and dimension
$2$, whose entries are all one except
$(\tilde{A}_{2})_{112}=(\tilde{A}_{2})_{212}=1.5$ and
$(\tilde{A}_{2})_{121}=(\tilde{A}_{2})_{221}=-0.5$. We further see that the
system (9) with $A_{1}$ and $\tilde{A}_{2}$ is equivalent to the system (9)
with $A_{1}$ and $A_{2}$, because they all yield the same component-wise
representation (14). However, from Corollary 1, for the system with
$\tilde{A}_{2}$, the conservative region of attraction is
$\max_{j}{|x_{j}(0)|}<\min_{p}\frac{1-\sum_{l}|A_{1}|_{p,l}}{\sum_{l,m}|A_{2}|_{p,lm}}=\frac{0.8}{4}=0.2<\frac{4}{15}$
and is smaller than the conservative region of attraction calculated by
$A_{2}$. We see that although one deals with the same system, the choice of
tensors will also influence the conservative region of attraction. To make the
region as large as possible, one needs to choose appropriate higher order
($\geq 3$) tensors such that
$\sum_{i_{2},\cdots,i_{k}=1}^{n}|(A_{i})_{p,I}|,i\geq 3$ is as small as
possible.

##  VII Conclusions

In this paper, we investigate general discrete-time non-homogeneous polynomial
dynamical systems on non-uniform hypergraphs. In particular, we give a simple
way to calculate the conservative region of attraction of the origin directly
with system parameters. Concerning the equilibrium other than the origin, one
just needs to calculate the error dynamic and then the main result can be
still applied to that error dynamic. Furthermore, we develop a feedback
control strategy that can manipulate the conservative region of attraction.
Finally, the main results are illustrated via numerical examples.

## References

  * [1] P. E. Paré, J. Liu, C. L. Beck, B. E. Kirwan, and T. Başar, “Analysis, estimation, and validation of discrete-time epidemic processes,” _IEEE Transactions on Control Systems Technology_ , vol. 28, no. 1, pp. 79–93, 2018.
  * [2] F. Liu, C. Shaoxuan, X. Li, and M. Buss, “On the stability of the endemic equilibrium of a discrete-time networked epidemic model,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 2576–2581, 2020.
  * [3] S. Cui, F. Liu, H. Jardón-Kojakhmetov, and M. Cao, “Discrete-time layered-network epidemics model with time-varying transition rates and multiple resources,” _arXiv preprint arXiv:2206.07425_ , 2022.
  * [4] P. Duarte, R. L. Fernandes, and W. M. Oliva, “Dynamics of the attractor in the lotka–volterra equations,” _journal of differential equations_ , vol. 149, no. 1, pp. 143–189, 1998.
  * [5] A. Slavík, “Lotka–volterra competition model on graphs,” _SIAM Journal on Applied Dynamical Systems_ , vol. 19, no. 2, pp. 725–762, 2020.
  * [6] B. Goh, “Global stability in two species interactions,” _Journal of Mathematical Biology_ , vol. 3, no. 3, pp. 313–318, 1976.
  * [7] ——, “Stability in models of mutualism,” _The American Naturalist_ , vol. 113, no. 2, pp. 261–275, 1979.
  * [8] P. Donnell and M. Banaji, “Local and global stability of equilibria for a class of chemical reaction networks,” _SIAM Journal on Applied Dynamical Systems_ , vol. 12, no. 2, pp. 899–920, 2013.
  * [9] G. Craciun, “Polynomial dynamical systems, reaction networks, and toric differential inclusions,” _SIAM Journal on Applied Algebra and Geometry_ , vol. 3, no. 1, pp. 87–106, 2019.
  * [10] D. Angeli, “A tutorial on chemical reaction networks dynamics,” in _2009 European Control Conference (ECC)_. IEEE, 2009, pp. 649–657.
  * [11] W. Ji and S. Deng, “Autonomous discovery of unknown reaction pathways from data by chemical reaction neural network,” _The Journal of Physical Chemistry A_ , vol. 125, no. 4, pp. 1082–1092, 2021.
  * [12] C. Chen, A. Surana, A. M. Bloch, and I. Rajapakse, “Controllability of hypergraphs,” _IEEE Transactions on Network Science and Engineering_ , vol. 8, no. 2, pp. 1646–1657, 2021.
  * [13] C. Chen and I. Rajapakse, “Tensor entropy for uniform hypergraphs,” _IEEE Transactions on Network Science and Engineering_ , vol. 7, no. 4, pp. 2889–2900, 2020.
  * [14] M. M. Wolf, A. M. Klinvex, and D. M. Dunlavy, “Advantages to modeling relational data using hypergraphs versus graphs,” in _2016 IEEE High Performance Extreme Computing Conference (HPEC)_. IEEE, 2016, pp. 1–7.
  * [15] C. Bick, E. Gross, H. A. Harrington, and M. T. Schaub, “What are higher-order networks?” _SIAM Review_ , vol. 65, no. 3, pp. 686–731, 2023.
  * [16] Y. Wang, Y. Wei, G. Zhang, and S. Y. Chang, “Algebraic riccati tensor equations with applications in multilinear control systems,” _arXiv preprint arXiv:2402.13491_ , 2024.
  * [17] I. Iacopini, G. Petri, A. Barrat, and V. Latora, “Simplicial models of social contagion,” _Nature communications_ , vol. 10, no. 1, p. 2485, 2019.
  * [18] P. Cisneros-Velarde and F. Bullo, “Multigroup sis epidemics with simplicial and higher order interactions,” _IEEE Transactions on Control of Network Systems_ , vol. 9, no. 2, pp. 695–705, 2021.
  * [19] S. Cui, F. Liu, H. Jardón-Kojakhmetov, and M. Cao, “General sis diffusion process with indirect spreading pathways on a hypergraph,” _arXiv preprint arXiv:2306.00619_ , 2023.
  * [20] A. D. Letten and D. B. Stouffer, “The mechanistic basis for higher-order interactions and non-additivity in competitive communities,” _Ecology letters_ , vol. 22, no. 3, pp. 423–436, 2019.
  * [21] S. Cui, Q. Zhao, H. Jardon-Kojakhmetov, and M. Cao, “Species coexistence and extinction resulting from higher-order lotka-volterra two-faction competition,” in _2023 62nd IEEE Conference on Decision and Control (CDC)_. IEEE, 2023, pp. 467–472.
  * [22] M. M. Mayfield and D. B. Stouffer, “Higher-order interactions capture unexplained complexity in diverse communities,” _Nature ecology & evolution_, vol. 1, no. 3, pp. 1–7, 2017.
  * [23] P. A. Abrams, “Arguments in favor of higher order interactions,” _The American Naturalist_ , vol. 121, no. 6, pp. 887–891, 1983.
  * [24] K. J. Åström and B. Wittenmark, _Computer-controlled systems: theory and design_. Courier Corporation, 2013.
  * [25] M. S. M. Saat, S. K. Nguang, and A. Nasiri, _Analysis and synthesis of polynomial discrete-time systems: an SOS approach_. Butterworth-Heinemann, 2017.
  * [26] K. Tanaka, H. Yoshida, H. Ohtake, and H. O. Wang, “A sum-of-squares approach to modeling and control of nonlinear dynamical systems with polynomial fuzzy systems,” _IEEE Transactions on Fuzzy systems_ , vol. 17, no. 4, pp. 911–922, 2008.
  * [27] P. A. Parrilo, _Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization_. California Institute of Technology, 2000.
  * [28] R. Mtar, M. Belhaouane, H. B. Ayadi, and N. B. Braiek, “An lmi criterion for the global stability analysis of nonlinear polynomial systems,” _Nonlinear Dynamics and Systems Theory_ , vol. 9, no. 2, pp. 171–183, 2009\.
  * [29] M. M. Belhaouane, M. F. Ghariani, H. Belkhiria Ayadi, N. Benhadj Braiek _et al._ , “Improved results on robust stability analysis and stabilization for a class of uncertain nonlinear systems,” _Mathematical Problems in Engineering_ , vol. 2010, 2010.
  * [30] D. Cheng and H. Qi, “Global stability and stabilization of polynomial systems,” in _2007 46th IEEE Conference on Decision and Control_. IEEE, 2007, pp. 1746–1751.
  * [31] F. Bullo, _Lectures on Network Systems_ , 1.6 ed. Kindle Direct Publishing, 2022. [Online]. Available: http://motion.me.ucsb.edu/book-lns
  * [32] K. Chang, L. Qi, and T. Zhang, “A survey on the spectral theory of nonnegative tensors,” _Numerical Linear Algebra with Applications_ , vol. 20, no. 6, pp. 891–912, 2013.
  * [33] K.-C. Chang, K. Pearson, and T. Zhang, “Perron-frobenius theorem for nonnegative tensors,” _Communications in Mathematical Sciences_ , vol. 6, no. 2, pp. 507–520, 2008.
  * [34] Y. Yang and Q. Yang, “Further results for perron–frobenius theorem for nonnegative tensors,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 31, no. 5, pp. 2517–2530, 2010.
  * [35] Q. Yang and Y. Yang, “Further results for perron–frobenius theorem for nonnegative tensors ii,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 32, no. 4, pp. 1236–1250, 2011.
  * [36] G. Gallo, G. Longo, S. Pallottino, and S. Nguyen, “Directed hypergraphs and applications,” _Discrete applied mathematics_ , vol. 42, no. 2-3, pp. 177–201, 1993.
  * [37] C. Chen, “On the stability of multilinear dynamical systems,” _arXiv preprint arXiv:2105.01041_ , 2021.
  * [38] S. Cui, G. Zhang, H. Jardón-Kojakhmetov, and M. Cao, “On metzler positive systems on hypergraphs,” _arXiv preprint arXiv:2401.03652_ , 2024.
  * [39] L. Qi, “Eigenvalues of a real supersymmetric tensor,” _Journal of Symbolic Computation_ , vol. 40, no. 6, pp. 1302–1324, 2005.
  * [40] M. Zhang, G. Ni, and G. Zhang, “Iterative methods for computing u-eigenvalues of non-symmetric complex tensors with application in quantum entanglement,” _Computational Optimization and Applications_ , vol. 75, pp. 779–798, 2020\.
  * [41] C. Mo and Y. Wei, “On nonnegative solution of multi-linear system with strong mz-tensors,” _Numer Math Theor Meth Appl_ , vol. 14, pp. 176–193, 2021.
  * [42] J. Xie and L. Qi, “Spectral directed hypergraph theory via tensors,” _Linear and Multilinear Algebra_ , vol. 64, no. 4, pp. 780–794, 2016.
  * [43] I. S. Gradshteyn and I. M. Ryzhik, _Table of integrals, series, and products_. Academic press, 2014.