2.5cm2.5cm2.5cm2.5cm

# Invariant circles and phase portraits of cubic vector fields on the sphere

Joji Benny  Department of Mathematics, Indian Institute of Technology Madras,
India<EMAIL_ADDRESS>,  Supriyo Jana  Department of Mathematics, Indian
Institute of Technology Madras, India<EMAIL_ADDRESS>and  Soumen
Sarkar  Department of Mathematics, Indian Institute of Technology Madras,
India<EMAIL_ADDRESS>
###### Abstract.

In this paper, we characterize and study dynamical properties of cubic vector
fields on the sphere
$\mathbb{S}^{2}=\\{(x,y,z)\in\mathbb{R}^{3}~{}|~{}x^{2}+y^{2}+z^{2}=1\\}$. We
start by classifying all degree three polynomial vector fields on
$\mathbb{S}^{2}$ and determine which of them form Kolmogorov systems. Then, we
show that there exist completely integrable cubic vector fields on
$\mathbb{S}^{2}$ and also study the maximum number of various types of
invariant circles for homogeneous cubic vector fields on $\mathbb{S}^{2}$. We
find a tight bound in each case. Further, we also discuss phase portraits of
certain cubic Kolmogorov vector fields on $\mathbb{S}^{2}$.

###### Key words and phrases:

polynomial vector fields, Kolmogorov system, periodic orbit, invariant circle,
invariant great circle, first integral, phase portrait

###### 2020 Mathematics Subject Classification:

34A34, 34C14, 34C40, 34C45, 58J90

##  1\. Introduction

Let $P,Q,R$ be polynomials in $\mathbb{R}[x,y,z]$. Then, the following system
of differential equations

(1.1) $\frac{dx}{dt}=P(x,y,z),\frac{dy}{dt}=Q(x,y,z),\frac{dz}{dt}=R(x,y,z)$

is called a polynomial differential system in $\mathbb{R}^{3}$. The
differential operator

(1.2) $\mathcal{X}=P\frac{\partial}{\partial x}+Q\frac{\partial}{\partial
y}+R\frac{\partial}{\partial z}$

is called the vector field associated with the system (1.1). The degree of the
polynomial vector field in (1.2) is defined to be
$\max\\{\deg(P),\deg(Q),\deg(R)\\}$.

The system (1.1) is called a polynomial Kolmogorov system in $\mathbb{R}^{3}$
when $P=xP^{\prime}$, $Q=yQ^{\prime}$ and $R=zR^{\prime}$ for some
$P^{\prime},Q^{\prime},R^{\prime}\in\mathbb{R}[x,y,z]$. The associated vector
field is called a polynomial Kolmogorov vector field.

An invariant algebraic set for (1.2) is a subset $A\subset\mathbb{R}^{3}$ such
that $A$ is the zero set of some $f(x,y,z)\in\mathbb{R}[x,y,z]$ and
$\mathcal{X}f=Kf$ for some $K\in\mathbb{R}[x,y,z]$. Here, the polynomial $K$
is called the cofactor of $f$. Moreover, $\mathcal{X}$ is also called a vector
field on $A$.

The Darboux theory of integrability and its generalizations [9] posits that if
a system in $\mathbb{R}^{2}$ has sufficiently many invariant algebraic sets,
then it has a rational first integral. This theory can be generalized in
$\mathbb{R}^{n}$ also, see [3]. Determining invariant algebraic curves for
vector fields on surfaces has become a problem of its own; see [8] and the
references therein. Perhaps the simplest invariant algebraic set that a vector
field on a sphere can have is an invariant circle since a circle on
$\mathbb{S}^{2}$ is always an intersection of a plane with $\mathbb{S}^{2}$.
An interesting problem involving invariant circles is determining when an
invariant circle is a limit cycle or at least when it is a periodic orbit.
Note that invariant circles for polynomial vector fields on $\mathbb{S}^{2}$
are discussed in [11] and [16] for quadratic vector fields.

On the other hand, the Kolmogorov system in $\mathbb{R}^{4}_{+}$ has been
studied in [14]. Then, the integrability of a class of Kolmogorov systems in
$\mathbb{R}^{n}$ has been explored in [12]. Importantly, Kolmogorov systems
[6] have found applications in Plasma Physics [7], Economics [1] and other
areas, including Ecology. In this paper, we characterize cubic polynomial
vector fields and study the dynamics of cubic Kolmogorov vector fields on
$\mathbb{S}^{2}$. We remark that the article [4] discusses phase portraits of
certain degree 3 polynomial vector fields on the Poincaré sphere.

This paper is organized in the following way. In Section 2, we discuss the
transformation of the vector fields under the stereographic projection
formulae for projection from the ‘South pole’. Then we recall the definition
of extactic polynomial, first integral and complete integrability. We show
that if $\gamma$ is an orbit of $\mathcal{X}$ on $S^{2}$, then the cone on
$\gamma$ is also invariant for $\mathcal{X}$. We recall the equation of such a
cone when $\gamma$ is a circle.

In Section 3, we give a necessary and sufficient condition when a cubic vector
field in $\mathbb{R}^{3}$ is a vector field on $\mathbb{S}^{2}$, see Theorem
3.2 and further also determine when the invariant great circle is a periodic
orbit. We characterize cubic Kolmogorov vector fields on $\mathbb{S}^{2}$ in
Corollary 3.3. We discuss the existence of an invariant great circle for a
cubic vector field on $\mathbb{S}^{2}$ in Theorem 3.6. We prove that there
exists a large class of completely integrable cubic vector fields on
$\mathbb{S}^{2}$, see Theorem 3.4. We prove that there exist homogeneous and
non-homogeneous cubic vector fields on $\mathbb{S}^{2}$ having invariant great
circle as a periodic orbit, see Proposition 3.7.

In Section 4, the focus of our discussion centers around homogeneous cubic
vector fields on $\mathbb{S}^{2}$. In Theorem 4.1, we find tight bounds on the
number of invariant great circles that a cubic homogeneous vector field on
$\mathbb{S}^{2}$ can have. We also exhibit a condition in Proposition 4.2 when
a cubic homogeneous vector field on $\mathbb{S}^{2}$ has an invariant circle
which is not a great circle.

In Section 5, we study cubic Kolmogorov systems on $\mathbb{S}^{2}$ in detail.
In particular, we draw phase portraits for various possible cubic Kolmogorov
systems and show that a broad class of such systems does not admit a periodic
orbit in Theorem 5.2. We give a sufficient condition when a singular point in
$\mathbb{S}^{2}\setminus\\{z=0\\}$ of a cubic Kolmogorov vector field on
$\mathbb{S}^{2}$ is either center or focus, see Theorem 5.3.

##  2\. Preliminaries

In this section, we discuss the transformation of the vector fields under the
stereographic projection formulae for projection from the ‘South pole’. Note
that these have been obtained in [16]. Then we recall the concept of extactic
polynomial which helps to find invariant hyper-surfaces. We recall the
definition of first integral and complete integrability. We show that if
$\gamma$ is an orbit of $\mathcal{X}$ on $S^{2}$, then the cone on $\gamma$ is
also invariant for $\mathcal{X}$. We recall the equation of such a cone when
$\gamma$ is a circle. Recall that the map

$\Phi\colon\mathbb{S}^{2}-\\{(0,0,-1)\\}\longrightarrow\mathbb{R}^{2}\quad\mbox{defined
by }(x,y,z)\mapsto(u,v)$

where,

$u=\frac{x}{1+z}\quad,\quad v=\frac{y}{1+z}$

is called the stereographic projection from the South pole. For a vector field
$\mathcal{X}=(P,Q,R)$ given by

$P=\sum_{i+j+k\leq n}p_{ijk}x^{i}y^{j}z^{k},\quad Q=\sum_{i+j+k\leq
n}q_{ijk}x^{i}y^{j}z^{k},\quad R=\sum_{i+j+k\leq n}r_{ijk}x^{i}y^{j}z^{k},$

the induced map $\Phi_{*}(P,Q,R)=(\mathcal{P},\mathcal{Q})$ is determined by

$\begin{split}\dot{u}=\mathcal{P}=\tilde{P}(u,v)-u\tilde{R}(u,v),\\\
\dot{v}=\mathcal{Q}=\tilde{Q}(u,v)-v\tilde{R}(u,v),\end{split}$

where,

$\begin{split}\tilde{P}(u,v)=\sum_{i+j+k\leq
n}p_{ijk}(2u)^{i}(2v)^{j}(1-u^{2}-v^{2})^{k}(u^{2}+v^{2}+1)^{n-i-j-k},\\\
\tilde{Q}(u,v)=\sum_{i+j+k\leq
n}q_{ijk}(2u)^{i}(2v)^{j}(1-u^{2}-v^{2})^{k}(u^{2}+v^{2}+1)^{n-i-j-k},\\\
\tilde{R}(u,v)=\sum_{i+j+k\leq
n}r_{ijk}(2u)^{i}(2v)^{j}(1-u^{2}-v^{2})^{k}(u^{2}+v^{2}+1)^{n-i-
j-k}.\end{split}$

One of the best tools in order to search for invariant algebraic hypersurfaces
is the following. Let $W$ be a vector subspace of
$\mathbb{R}[x_{1},...,x_{n}]$ generated by the independent polynomials
$v_{1},...,v_{l}$, i.e., $W=\langle v_{1},...,v_{\ell}\rangle$. The extactic
polynomial of $\mathcal{X}$ associated with $W$ is the polynomial

$\mathcal{E}_{W}(\mathcal{X})=\begin{vmatrix}v_{1}&\dots&v_{\ell}\\\
\mathcal{X}(v_{1})&\dots&\mathcal{X}(v_{\ell})\\\ \vdots&\vdots&\vdots\\\
\mathcal{X}^{{\ell}-1}(v_{1})&\dots&\mathcal{X}^{{\ell}-1}(v_{\ell})\end{vmatrix}$

where $\mathcal{X}^{j}(v_{i})=\mathcal{X}^{j-1}(\mathcal{X}(v_{i}))$ for all
$i,j$. From the properties of determinants, it follows that the definition of
the extactic polynomial is independent of the chosen basis of $W$.

######  Proposition 2.1.

[10, Proposition 1] Let $\mathcal{X}$ be a polynomial vector field in
$\mathbb{R}^{n}$ and $W$ a finite dimensional vector subspace of
$\mathbb{R}[x_{1},x_{2},\ldots,x_{n}]$ with $\dim(W)>1$. If $f=0$ is an
invariant algebraic hypersurface for the vector field $\mathcal{X}$ with $f\in
W$, then $f$ is a factor of $\mathcal{E}_{W}(\mathcal{X})$.

The multiplicity of an invariant algebraic hypersurface $f=0$ with $f\in W$,
is the largest positive integer $k$ such that $f^{k}$ divides the extactic
polynomial $\mathcal{E}_{W}(\mathcal{X})$ when
$\mathcal{E}_{W}(\mathcal{X})\neq 0$, otherwise the multiplicity is infinite.
For more details on this multiplicity, see [2, 15].

######  Definition 2.2.

Let $U$ be an open subset of $\mathbb{R}^{3}$. A non-constant analytic map
$H\colon U\to\mathbb{R}$ is called a first integral of the vector field (1.2)
on $U$ if $H$ is constant on all solution curves of the system (1.1) contained
in $U$ ; i.e., $H(x(t),y(t),z(t))=$ constant for all values of $t$ for which
the solution $(x(t),y(t),z(t))$ is defined and contained in $U$.

Note that $H$ is a first integral of the vector field (1.2) on $U$ if and only
if $\mathcal{X}H=0$ on $U$.

######  Definition 2.3.

The system (1.1) is called completely integrable in an open set $U$ if it has
2 independent first integrals.

If the system (1.1) is completely integrable with 2 independent first
integrals $H_{1}$ and $H_{2}$ then the orbits of the system are contained in
$\\{H_{1}=c_{1}\\}\cap\\{H_{2}=c_{2}\\}$ for some $c_{1},c_{2}\in\mathbb{R}$.
In [13], the complete integrability of vector fields in $\mathbb{R}^{n}$ is
discussed.

Recall that the intersection of a plane $ax+by+cz+d=0$ with $\mathbb{S}^{2}$
is a circle. We can choose $a,b,c,d$ such that $a^{2}+b^{2}+c^{2}=1$ and
$|d|<1.$ Any circle on $\mathbb{S}^{2}$ can be obtained in this way. If the
plane passes through the origin, the intersection is called a great circle.
Note that a vector field on $\mathbb{S}^{2}$ is invariant by $SO(3)$111$SO(3)$
is the group of all rotations about the origin in $\mathbb{R}^{3}$.. Hence, if
a polynomial vector field on $\mathbb{S}^{2}$ has an invariant circle
$\\{ax+by+cz+d=0\\}\cap\mathbb{S}^{2}$, then we can assume the circle to be
$\\{z+d=0\\}\cap\mathbb{S}^{2}$ with $|d|<1$.

We now state a result for homogeneous polynomial vector fields on
$\mathbb{S}^{2}$, which is proved for degree two in [11], and an entirely
similar proof can be given for homogeneous vector fields of any degree.

######  Proposition 2.4.

Let $\gamma=\\{\phi(t)|t\in\mathbb{R}\\}\subset\mathbb{S}^{2}$ be an orbit of
$\mathcal{X}.$ If $\mathcal{X}$ is a homogeneous polynomial vector field on
$\mathbb{S}^{2},$ then $\mathcal{X}$ is tangent to the surface
$S(\gamma)=\\{sp~{}|~{}s\in\mathbb{R},p\in\gamma\\}.$

We want to look at the circles on $\mathbb{S}^{2}$ which are invariant with
respect to the flow of the vector field $\mathcal{X}$ on $S^{2}$. In this
case, Proposition 2.4 will imply that the entire cone on the circle is also
invariant with respect to the flow of $\mathcal{X}$ if the circle is invariant
by the flow of $\mathcal{X}.$ In [11], the authors describe the equation of
the cone of such a circle as

(2.1)
$(a^{2}-d^{2})x^{2}+(b^{2}-d^{2})y^{2}+(c^{2}-d^{2})z^{2}+2abxy+2acxz+2bcyz=0$

if $d\neq 0.$ If $d=0,$ then the cone itself becomes the plane given by
$ax+by+cz=0.$ When $d=0,$ the intersection of $\\{ax+by+cz=0\\}$ with
$\mathbb{S}^{2}$ is a great circle.

##  3\. Cubic vector fields on $\mathbb{S}^{2}$

In this section, we characterize cubic vector fields on the standard sphere in
$\mathbb{R}^{3}$. Then we classify cubic Kolmogorov vector fields on $S^{2}$.
We also study those cubic vector fields which have an invariant great circle.

######  Lemma 3.1.

[5, Lemma 4.1] Let $n\in\mathbb{N}$ and $Q_{1},Q_{2},Q_{3}$ are polynomials in
$\mathbb{R}[x,y,z]$ such that the polynomial
$Q_{1}x^{n}+Q_{2}y^{n}+Q_{3}z^{n}$ is zero. Then $Q_{1}=Ay^{n}+Bz^{n}$,
${Q_{2}=-Ax^{n}+Cz^{n}}$ and $Q_{3}=-Bx^{n}-Cy^{n}$ for some polynomials
$A,B,C$.

######  Theorem 3.2.

Let $\mathcal{X}=(P,Q,R)$ be a cubic polynomial vector field in
$\mathbb{R}^{3}$. Then $\mathcal{X}$ is cubic vector field on $\mathbb{S}^{2}$
if and only if there exist $f,g,h,A,B,C\in\mathbb{R}[x,y,z]$ such that

(3.1) $\begin{split}P&=(1-x^{2}-y^{2}-z^{2})f+Ay+Bz,\\\
Q&=(1-x^{2}-y^{2}-z^{2})g-Ax+Cz,\quad\text{and}\\\
R&=(1-x^{2}-y^{2}-z^{2})h-Bx-Cy\end{split}$

where $f,g,h$ are linear polynomials and $A,B,C$ are quadratic polynomials
without any constant term. Moreover, this vector field has cofactor
$-2(fx+gy+hz)$ for $\mathbb{S}^{2}$.

###### Proof.

We write $P=P^{(3)}+P^{(2)}+P^{(1)}+P^{(0)}$ where $P^{(j)}$ is the degree $j$
homogeneous part of $P$. Similarly, we write
$Q=Q^{(3)}+Q^{(2)}+Q^{(1)}+P^{(0)}$ and $R=R^{(3)}+R^{(2)}+R^{(1)}+R^{(0)}$
also in this fashion.

Suppose $\mathcal{X}=(P,Q,R)$ is a vector field on $\mathbb{S}^{2}$. Then, it
must satisfy

(3.2) $2(Px+Qy+Rz)=K(x^{2}+y^{2}+z^{2}-1).$

for some $K\in\mathbb{R}[x,y,z]$ with $\deg(K)\leq 2$. Assume that $K^{(j)}$
is the $j$ degree homogeneous part of $K$. Then $K^{(0)}$ is $0$ since there
is no constant term on the left side of (3.2). Now, Comparing the degree 4,
degree 3, degree 2 and degree 1 terms in (3.2), we get

$\begin{split}2(P^{(3)}x+Q^{(3)}y+R^{(3)}z)&=K^{(2)}(x^{2}+y^{2}+z^{2}),\\\
2(P^{(2)}x+Q^{(2)}y+R^{(2)}z)&=K^{(1)}(x^{2}+y^{2}+z^{2}),\\\
2(P^{(1)}x+Q^{(1)}y+R^{(1)}z)&=-K^{(2)},\\\
2(P^{(0)}x+Q^{(0)}y+R^{(0)}z)&=-K^{(1)},\end{split}$

respectively. Hence

$P^{(3)}x+Q^{(3)}y+R^{(3)}z=-(P^{(1)}x+Q^{(1)}y+R^{(1)}z)(x^{2}+y^{2}+z^{2}),\quad\text{and}$
$P^{(2)}x+Q^{(2)}y+R^{(2)}z=-(P^{(0)}x+Q^{(0)}y+R^{(0)}z)(x^{2}+y^{2}+z^{2}).$

By Lemma 3.1,

$\begin{split}P^{(3)}&=-P^{(1)}(x^{2}+y^{2}+z^{2})+A_{1}y+B_{1}z,\\\
Q^{(3)}&=-Q^{(1)}(x^{2}+y^{2}+z^{2})-A_{1}x+C_{1}z,\\\
R^{(3)}&=-R^{(1)}(x^{2}+y^{2}+z^{2})-B_{1}x-C_{1}y\end{split}$

where each of $A_{1},B_{1},C_{1}$ is either 0 or a quadratic homogeneous
polynomial in $\mathbb{R}[x,y,z]$, and

$\begin{split}P^{(2)}&=-P^{(0)}(x^{2}+y^{2}+z^{2})+A_{2}y+B_{2}z\\\
Q^{(2)}&=-Q^{(0)}(x^{2}+y^{2}+z^{2})-A_{2}x+C_{2}z\\\
R^{(2)}&=-R^{(0)}(x^{2}+y^{2}+z^{2})-B_{2}x-C_{2}y.\end{split}$

where each of $A_{2},B_{2},C_{2}$ is either 0 or a linear homogeneous
polynomial in $\mathbb{R}[x,y,z]$. Hence

$P=P^{(3)}+P^{(2)}+P^{(1)}+P^{(0)}=(1-x^{2}-y^{2}-z^{2})(P^{(1)}+P^{(0)})+(A_{1}+A_{2})y+(B_{1}+B_{2})z.$

Denoting $f:=P^{(1)}+P^{(0)}$, $A:=A_{1}+A_{2}$ and $B:=B_{1}+B_{2}$; $P$
becomes $(1-x^{2}-y^{2}-z^{2})f+Ay+Bz$. Similarly, we get $Q$ and $R$ as in
the form in (3.1).

Note that if $P,Q,R$ are given by (3.1), then they satisfy (3.2) with
${K=-2(fx+gy+hz)}$. Hence, $\mathcal{X}=(P,Q,R)$ is a vector field on
$\mathbb{S}^{2}$. Thus, the converse part is true. ∎

######  Corollary 3.3.

Suppose $\mathcal{X}=(P,Q,R)$ is the cubic Kolmogorov vector field in
$\mathbb{R}^{3}$. Then $\mathcal{X}$ is a vector field on $\mathbb{S}^{2}$ if
and only if there exist $\alpha,\beta,\gamma,a,b,c\in\mathbb{R}$ such that

(3.3) $\begin{split}P&=x(\alpha(1-x^{2}-y^{2}-z^{2})+ay^{2}+bz^{2}),\\\
Q&=y(\beta(1-x^{2}-y^{2}-z^{2})-ax^{2}+cz^{2}),\quad\text{and}\\\
R&=z(\gamma(1-x^{2}-y^{2}-z^{2})-bx^{2}-cy^{2}).\end{split}$

###### Proof.

Suppose $\mathcal{X}=(P,Q,R)$ is a vector field on $\mathbb{S}^{2}$. Then by
Theorem 3.2, $P,Q,R$ are given by (3.1). Since $\mathcal{X}$ is Kolmogorov,
$x,y,z$ divide $P,Q,R$, respectively. Suppose

$P^{\prime}x=P=(1-x^{2}-y^{2}-z^{2})f+Ay+Bz.$

Then, we get $P^{\prime}(0,0,0)x=f$ by comparing the linear part on both
sides. Assume that ${P^{\prime}(0,0,0)=\alpha}$, we get $f=\alpha x$.
Similarly, for $Q^{\prime}y=Q$ and $R^{\prime}z=R$, one gets $g=\beta y$ and
$h=\gamma z$ for some $\beta,\gamma\in\mathbb{R}$. Therefore, $x$ divides
$Ay+Bz$, $y$ divides $-Ax+Cz$, and $z$ divides $-Bx-Cy$. Suppose

$P_{1}x=Ay+Bz,\quad Q_{1}y=-Ax+Cz\quad\text{and}\quad R_{1}z=-Bx-Cy$

for some polynomials $P_{1},Q_{1},$ and $R_{1}$ with
$\deg(P_{1}),\deg(Q_{1}),\deg(R_{1})\leq 2$. Then, we get that
$P_{1}x^{2}+Q_{1}y^{2}+R_{1}z^{2}=x(Ay+Bz)+y(-Ax+Cz)+z(-Bx-Cy)=0$. So, by
Lemma 3.1, $P_{1}=ay^{2}+bz^{2},Q_{1}=-ax^{2}+cz^{2}$ and
$R_{1}=-bx^{2}-cy^{2}$ for some $a,b,c\in\mathbb{R}$. Thus, we get $P,Q,R$ as
in (3.3).

Suppose $P,Q,R$ are given by (3.3). Then, the associated differential system
of the vector field ${\mathcal{X}=(P,Q,R)}$ is a cubic Kolmogorov vector
field. Also, by Theorem 3.2, $\mathcal{X}$ is a vector field on
$\mathbb{S}^{2}$. Thus, the converse part is true. ∎

######  Theorem 3.4.

There exist completely integrable cubic vector fields on $\mathbb{S}^{2}$.

###### Proof.

By Lemma 3.1, $f^{\prime}x+g^{\prime}y+h^{\prime}z=0$ for some
$f^{\prime},g^{\prime},h^{\prime}\in\mathbb{R}[x,y,z]$ if and only if
$f^{\prime}=\alpha y+\beta z$, $g^{\prime}=-\alpha x+\gamma z$ and
${h^{\prime}=-\beta x-\gamma y}$ for some $\alpha,\beta,\gamma\in\mathbb{R}$.
Suppose $\gamma\neq 0$ and $a\in\mathbb{R}\setminus\\{0\\}$. Then there exist
$b,c\in\mathbb{R}$ such that $\alpha=\frac{c\gamma}{a}$ and
$\beta=-\frac{b\gamma}{a}$. Define $f:=\frac{c\gamma}{a}y-\frac{b\gamma}{a}z$,
$g:=-\frac{c\gamma}{a}x+\gamma z$ and $h:=\frac{b\gamma}{a}x-\gamma y$. Then

$af+bg+ch=0\quad\mbox{and}\quad fx+gy+hz=0.$

Now, we choose $A=\frac{c}{a}C$ and $B=-\frac{b}{a}C$ for any quadratic
polynomial $C$ with $C(0,0,0)=0$. Consider the vector field
$\mathcal{X}=(P,Q,R)$ given by

$\begin{split}P&=(1-x^{2}-y^{2}-z^{2})f+Ay+Bz,\\\
Q&=(1-x^{2}-y^{2}-z^{2})g-Ax+Cz,\quad\text{and}\\\
R&=(1-x^{2}-y^{2}-z^{2})h-Bx-Cy.\end{split}$

By Theorem 3.2, $\mathcal{X}$ is a vector field on $\mathbb{S}^{2}$ with
cofactor $-2(fx+gy+hz)$. Hence, $x^{2}+y^{2}+z^{2}-1$ is a first integral of
$\mathcal{X}$. Again,

$\begin{split}\mathcal{X}(ax+by+cz)=&aP+bQ+cR\\\
=&(1-x^{2}-y^{2}-z^{2})(af+bg+ch)+a(Ay+Bz)+b(-Ax+Cz)+c(-Bx-Cy)\\\ =&(-bA-
cB)x+(aA-cC)y+(aB+bC)z=0.\end{split}$

Hence, $\mathcal{X}$ has two independent first integrals $x^{2}+y^{2}+z^{2}-1$
and $ax+by+cz$, which makes $\mathcal{X}$ completely integrable in
$\mathbb{R}^{3}$. ∎

######  Corollary 3.5.

If one of $a,b,c$ is non-zero, then there exists a cubic vector field on
$\mathbb{S}^{2}$ which has a first integral $ax+by+cz$.

The rest of this section characterizes vector fields on $\mathbb{S}^{2}$ which
have invariant great circles. Suppose $P,Q,R$ are given by (3.1), where

(3.4) $A=\sum\limits_{1\leq i+j+k\leq
2}a_{ijk}x^{i}y^{j}z^{k},~{}~{}B=\sum\limits_{1\leq i+j+k\leq
2}b_{ijk}x^{i}y^{j}z^{k},~{}\mbox{and}~{}C=\sum\limits_{1\leq i+j+k\leq
2}c_{ijk}x^{i}y^{j}z^{k}.$

Then, under the stereographic projection, the system given by (3.1) becomes

$\begin{split}\dot{u}&=\tilde{P}(u,v)-u\tilde{R}(u,v)=\mathcal{P}(u,v),\\\
\dot{v}&=\tilde{Q}(u,v)-v\tilde{R}(u,v)=\mathcal{Q}(u,v),\text{ where
}\end{split}$
$\begin{split}\tilde{P}=&((u^{2}+v^{2}+1)^{2}-(2u)^{2}-(2v)^{2}-(1-u^{2}-v^{2})^{2})\tilde{f}+\tilde{A}(2v)+\tilde{B}(1-u^{2}-v^{2})\\\
=&2\tilde{A}v+\tilde{B}(1-u^{2}-v^{2}).\end{split}$

Similarly, $\tilde{Q}=-2\tilde{A}u+\tilde{C}(1-u^{2}-v^{2})$ and
$\tilde{R}=-2\tilde{B}u-2\tilde{C}v$ where

$\begin{split}\tilde{A}=\sum\limits_{1\leq i+j+k\leq
2}a_{ijk}(2u)^{i}(2v)^{j}(1-u^{2}-v^{2})^{k}(u^{2}+v^{2}+1)^{2-i-j-k},\\\
\tilde{B}=\sum\limits_{1\leq i+j+k\leq
2}b_{ijk}(2u)^{i}(2v)^{j}(1-u^{2}-v^{2})^{k}(u^{2}+v^{2}+1)^{2-i-j-k},\\\
\tilde{C}=\sum\limits_{1\leq i+j+k\leq
2}c_{ijk}(2u)^{i}(2v)^{j}(1-u^{2}-v^{2})^{k}(u^{2}+v^{2}+1)^{2-i-j-k}.\end{split}$

######  Theorem 3.6.

Let $\mathcal{X}=(P,Q,R)$ be a cubic vector field on $\mathbb{S}^{2}$. Assume
$\mathcal{X}$ has an invariant great circle. Without loss of generality, we
can assume that it is $\mathbb{S}^{1}=\\{z=0\\}\cap\mathbb{S}^{2}$. Then, the
following hold.

  1. (a)

The vector field $\mathcal{X}$ can be written as (3.1) with

$B=b_{020}y^{2}+b_{110}xy+b_{010}y+B^{{}^{\prime}}z,C=-b_{110}x^{2}-b_{020}xy-
b_{010}x+C^{{}^{\prime}}z$

where $B^{{}^{\prime}},C^{{}^{\prime}}$ are linear polynomials in
$\mathbb{R}[x,y,z]$.

  2. (b)

The great circle $\mathbb{S}^{1}$ is a periodic orbit of $\mathcal{X}$ if and
only if the hypersurfaces
$\left\\{\sum\limits_{1\leq i+j\leq 2}a_{ij0}u^{i}v^{j}=0\right\\}$ and
$\\{u^{2}+v^{2}=1\\}$ do not intersect in $\mathbb{R}^{2}$.

###### Proof.

(a) If $\mathcal{X}=(P,Q,R)$ given by (3.1) with $A,B,C$ as in (3.4) is a
vector field on $\mathbb{S}^{2}$ then

$xP(x,y,z)+yQ(x,y,z)+zR(x,y,z)=0\text{ for all $(x,y,z)\in\mathbb{S}^{2}$}.$

Hence, after the stereographic projection, we get

$2u\tilde{P}(u,v)+2v\tilde{Q}(u,v)+(1-u^{2}-v^{2})\tilde{R}(u,v)=0\text{ for
all $(u,v)\in\mathbb{R}^{2}$}.$

So, by the above,

(3.5)
$u\dot{u}+v\dot{v}=u\tilde{P}(u,v)+v\tilde{Q}(u,v)-(u^{2}+v^{2})\tilde{R}(u,v)=-\frac{1}{2}(u^{2}+v^{2}+1)\tilde{R}(u,v).$

Under the stereographic projection $\\{z=0\\}\cap\mathbb{S}^{2}$ reduces to
unit circle $u^{2}+v^{2}=1$ in $\mathbb{R}^{2}$. Due to (3.5), one can check
that $u^{2}+v^{2}-1=0$ is invariant of the vector field
$(\mathcal{P},\mathcal{Q})$ if and only if $u^{2}+v^{2}-1$ divides
$\tilde{R}=-2(\tilde{B}u+\tilde{C}v)$. Again, $u^{2}+v^{2}-1$ divides
$\tilde{R}$ if and only if $\tilde{R}|_{u^{2}+v^{2}=1}=0$. Note that

$\begin{split}\tilde{R}|_{u^{2}+v^{2}=1}=&-2(u\tilde{B}|_{u^{2}+v^{2}=1}+v\tilde{C}|_{u^{2}+v^{2}=1})\\\
=&-8\left(\sum\limits_{1\leq i+j\leq
2}b_{ij0}u^{i+1}v^{j}+c_{ij0}u^{i}v^{j+1}\right).\end{split}$

Hence, $\tilde{R}|_{u^{2}+v^{2}=1}=0$ implies that

$b_{100}=b_{200}=c_{010}=c_{020}=0,c_{100}=-b_{010},c_{110}=-b_{020},c_{200}=-b_{110}.$

Hence, the first claim follows.

(b) It is well known that the invariant great circle
$\mathbb{S}^{1}\subset\mathbb{S}^{2}$ is a periodic orbit of $\mathcal{X}$ on
$\mathbb{S}^{2}$ if and only if there is no singular point on it. Observe
that, if $\mathbb{S}^{1}$ is invariant then $\tilde{R}|_{u^{2}+v^{2}=1}=0$.
Hence,

$\mathcal{P}|_{u^{2}+v^{2}=1}=2v\tilde{A}|_{u^{2}+v^{2}=1}\text{ and
}\mathcal{Q}|_{u^{2}+v^{2}=1}=-2u\tilde{A}|_{u^{2}+v^{2}=1}.$

Note that $\tilde{A}|_{u^{2}+v^{2}=1}=4\sum\limits_{1\leq i+j\leq
2}a_{ij0}u^{i}v^{j}$. So, $\mathbb{S}^{1}$ is a periodic orbit if and only if
the following two equations do not have a common solution.

$\sum\limits_{1\leq i+j\leq 2}a_{ij0}u^{i}v^{j}=0,~{}u^{2}+v^{2}=1.$

∎

Llibre [16, Theorem 2] showed that a quadratic homogeneous vector field on
$\mathbb{S}^{2}$ has no great circle which is a periodic orbit. This statement
does not hold true for cubic homogeneous vector fields on $\mathbb{S}^{2}$.

######  Proposition 3.7.

  1. (a)

There exists a homogeneous vector field having a great circle as a periodic
orbit.

  2. (b)

There exists a non-homogeneous vector field having a great circle as a
periodic orbit.

###### Proof.

(a) Consider the vector field $\mathcal{X}$ on $\mathbb{S}^{2}$ given by (3.1)
with

$A=x^{2}+y^{2},B=y^{2}+xy,C=-x^{2}-xy\quad\text{and}\quad f=g=h=0.$

By (a) of Theorem 3.6, the intersection $\\{z=0\\}\cap\mathbb{S}^{2}$ is an
invariant great circle for $\mathcal{X}$. Also, by (b) of Theorem 3.6, the
great circle is periodic since $u^{2}+v^{2}=0$ and $u^{2}+v^{2}=1$ do not have
a common solution.

(b) Consider the vector field $\mathcal{X}$ on $\mathbb{S}^{2}$ given by (3.1)
with

$A=x^{2}+y^{2},\quad B=y^{2}+xy,\quad C=-x^{2}-xy$

and at least one of $f,g,h$ is non-zero. By (a) of Theorem 3.6,
$\\{z=0\\}\cap\mathbb{S}^{2}$ is an invariant great circle for $\mathcal{X}$.
Also, by (b) of Theorem 3.6, the great circle is periodic since
$u^{2}+v^{2}=0$ and $u^{2}+v^{2}=1$ do not have a common solution. ∎

##  4\. Cubic homogeneous vector fields on $\mathbb{S}^{2}$

In this section, we study some properties of the vector fields

$\mathcal{X}:=P\frac{\partial}{\partial x}+Q\frac{\partial}{\partial
y}+R\frac{\partial}{\partial z},$

defined on
$\mathbb{S}^{2}=\\{(x,y,z)\in\mathbb{R}~{}|~{}f=x^{2}+y^{2}+z^{2}-1=0\\}$ such
that $P,Q,R$ are homogeneous cubic polynomials in $\mathbb{R}[x,y,z]$. Hence

(4.1) $\mathcal{X}(f)=K_{f}(f)$

where $K_{f}$ is the cofactor. We see that

$\mathcal{X}(f)=2Px+2Qy+2Rz.$

So, $\mathcal{X}$ is homogeneous implies that $\mathcal{X}(f)$ is also
homogeneous, but we see that the right-hand side of (4.1) is not homogeneous
unless $K_{f}=0.$ Therefore for a homogeneous vector field
$\mathcal{X}=(P,Q,R)$ on $\mathbb{S}^{2},$ we have

$Px+Qy+Rz=0.$

Hence by Lemma 3.1, the associated differential system of a homogeneous vector
field $\mathcal{X}=(P,Q,R)$ on $\mathbb{S}^{2}$ will be of the following form.

(4.2) $P=Ay+Bz,\quad Q=-Ax+Cz,\quad\mbox{and}\quad R=-Bx-Cy$

for some polynomials $A,B,C$ in $R[x,y,z]$.

Next, we discuss invariant circles on $\mathbb{S}^{2}$ for homogeneous cubic
vector fields. We see from Proposition 2.4 that a circle on $\mathbb{S}^{2}$
is invariant with respect to the flow of $\mathcal{X}$ if and only if the cone
of the circle is invariant with respect to the flow of $\mathcal{X}.$ This
implies that the circle given by the intersection of $\\{ax+by+cz+d=0\\}$ with
$\mathbb{S}^{2}$ is invariant if and only if

(4.3)
$\begin{split}\mathcal{X}\left((a^{2}-d^{2})x^{2}+(b^{2}-d^{2})y^{2}+(c^{2}-d^{2})z^{2}+2abxy+2acxz+2bcyz\right)=\\\
K\left((a^{2}-d^{2})x^{2}+(b^{2}-d^{2})y^{2}+(c^{2}-d^{2})z^{2}+2abxy+2acxz+2bcyz\right)\end{split}$

for some polynomial $K\in\mathbb{R}[x,y,z]$ where $a^{2}+b^{2}+c^{2}=1$ and
$|d|<1$. We remark that The cones of the invariant circles
$\\{ax+by+cz+d=0\\}\cap\mathbb{S}^{2}$ and
$\\{ax+by+cz-d=0\\}\cap\mathbb{S}^{2}$ are same. Hence, invariant circles
which are not great circles always occur in pairs.

In the case of a great circle $\\{ax+by+cz=0\\}\cap\mathbb{S}^{2},$ the cone
is the plane $\\{ax+by+cz=0\\}$ itself, and then the great circle is invariant
with respect to the flow of $\mathcal{X}$ if and only if

(4.4) $\mathcal{X}(ax+by+cz)=K^{\prime}(ax+by+cz)$

for some polynomial $K^{\prime}\in\mathbb{R}[x,y,z].$

######  Theorem 4.1.

Let $\mathcal{X}=(P,Q,R)$ be a cubic homogeneous vector field on
$\mathbb{S}^{2}$. Assume $\mathcal{X}$ has an invariant great circle. Without
loss of generality, we can assume that it is
$\mathbb{S}^{1}=\\{z=0\\}\cap\mathbb{S}^{2}$. If $\mathcal{X}$ has finitely
many invariant great circles, then the maximum number of invariant great
circles of the form $\\{ax+by+cz=0\\}\cap\mathbb{S}^{2}$ is

  1. (1)

3 if $a=0$ or $b=0$,

  2. (2)

2 if $c=0$,

  3. (3)

2 if $a,b,c$ are non-zero.

Moreover, in each of the above cases, the bound can be reached.

###### Proof.

The vector field $\mathcal{X}$ is given by (4.2). Since $z=0$ is invariant,
$z$ divides $R=-Bx-Cy$. Hence, by Lemma 3.1, $B=Ly+B^{\prime}z$ and
$C=-Lx+C^{\prime}z$ for some $L,B^{\prime},C^{\prime}\in\mathbb{R}[x,y,z]$ of
degree less than or equal to 1. Then

$\begin{split}P&=Ay+(Ly+B^{\prime}z)z=(A+Lz)y+B^{\prime}z^{2}=A^{\prime}y+B^{\prime}z^{2},\\\
Q&=-Ax+(-Lx+C^{\prime}z)z=-(A+Lz)x+C^{\prime}z^{2}=-A^{\prime}x+C^{\prime}z^{2},\\\
R&=-(Ly+B^{\prime}z)x-(-Lx+C^{\prime}z)y=-z(B^{\prime}x+C^{\prime}y).\end{split}$

where $A^{\prime}:=A+Lz$ is a homogeneous polynomial of degree less than or
equal to 2. We write

$\begin{split}A^{\prime}&=a_{1}x^{2}+a_{2}y^{2}+a_{3}z^{2}+a_{4}xy+a_{5}xz+a_{6}yz,\\\
B^{\prime}&=b_{1}x+b_{2}y+b_{3}z,\mbox{ and}\\\
C^{\prime}&=c_{1}x+c_{2}y+c_{3}z.\end{split}$

(1) The extactic polynomial of $\mathcal{X}$ associated with $\langle
y,z\rangle$ is

$\mathcal{E}_{\langle y,z\rangle}(\mathcal{X})=\begin{vmatrix}y&z\\\
Q&R\end{vmatrix}=z(-C^{\prime}(y^{2}+z^{2})+x(A^{\prime}-B^{\prime}y)).$

Suppose $\mathcal{E}_{\langle y,z\rangle}(\mathcal{X})$ has $\ell$ factors of
the form $by+cz$, say $b_{i}y+c_{i}z$ for $i=1,2,3,4$. Then

$z(-C^{\prime}(y^{2}+z^{2})+x(A^{\prime}-B^{\prime}y))=\prod\limits_{i=1}^{4}(b_{i}y+c_{i}z).$

For $x=0$, the above equation becomes
$-zC^{\prime}(0,y,z)(y^{2}+z^{2})=\prod\limits_{i=1}^{4}(b_{i}y+c_{i}z)$. But
this is not possible since $y^{2}+z^{2}$ has no linear factor over
$\mathbb{R}$. So, $\mathcal{E}_{\langle y,z\rangle}(\mathcal{X})$ has at most
3 factors of the form $by+cz$. Hence, for $a=0$ by the Proposition 2.1, the
maximum number of invariant great circles of the form
$\\{by+cz=0\\}\cap\mathbb{S}^{2}$ is 3.

To show that the bound can be reached, consider
$A^{\prime}=\prod\limits_{i=1}^{2}(b_{i}y+c_{i}z)+B^{\prime}y$, $C^{\prime}=0$
for any linear homogeneous polynomial $B^{\prime}$. Then the vector field

$\mathcal{X}=(P,Q,R)=\left(y\prod\limits_{i=1}^{2}(b_{i}y+c_{i}z)+B^{\prime}(y^{2}+z^{2}),-x\prod\limits_{i=1}^{2}(b_{i}y+c_{i}z)-B^{\prime}xy,-B^{\prime}xz\right)$

has invariant hyperplanes $z=0$, $b_{1}y+c_{1}z=0$, and $b_{2}y+c_{2}z=0$.

The proof for $b=0$ is similar to the proof of $a=0$. To show that the bound
can be reached, consider
$A^{\prime}=\prod\limits_{i=1}^{2}(a_{i}x+c_{i}z)-C^{\prime}y$, $B^{\prime}=0$
for any linear homogeneous polynomial $C^{\prime}$. Then the vector field

$\mathcal{X}=(P,Q,R)=\left(y\prod\limits_{i=1}^{2}(a_{i}x+c_{i}z)-C^{\prime}y^{2},-x\prod\limits_{i=1}^{2}(a_{i}x+c_{i}z)+C^{\prime}xy+C^{\prime}z^{2},-C^{\prime}yz\right)$

has invariant hyperplanes $z=0$, $a_{1}x+c_{1}z=0$, and $a_{2}x+c_{2}z=0$.

(2) The extactic polynomial of $\mathcal{X}$ associated with $\langle
x,y\rangle$ is

$\mathcal{E}_{\langle x,y\rangle}(\mathcal{X})=\begin{vmatrix}x&y\\\
P&Q\end{vmatrix}=-A^{\prime}(x^{2}+y^{2})+z^{2}(C^{\prime}x-B^{\prime}y).$

Suppose $\mathcal{E}_{\langle x,y\rangle}(\mathcal{X})$ has $4$ factors of the
form $ax+by$, say $a_{i}x+b_{i}y$ for $i=1,2,..,\ell\leq 4$. Then

(4.5)
$-A^{\prime}(x^{2}+y^{2})+z^{2}(C^{\prime}x-B^{\prime}y)=p\prod\limits_{i=1}^{\ell}(a_{i}x+b_{i}y)$

for some polynomial $p\in\mathbb{R}[x,y,z]$. In particular, for $z=0$, last
equation corresponds that

$-A^{\prime}(x,y,0)(x^{2}+y^{2})=p(x,y,0)\prod\limits_{i=1}^{\ell}(a_{i}x+b_{i}y).$

If $A^{\prime}(x,y,0)\neq 0$ then $\ell\leq 2$. Otherwise,
$A^{\prime}=zA^{\prime\prime}$ and $p=zp^{\prime}$ for some polynomials
$A^{\prime\prime},p^{\prime}\in\mathbb{R}[x,y,z]$. Then (4.5) becomes

$-A^{\prime\prime}(x^{2}+y^{2})+z(C^{\prime}x-B^{\prime}y)=p^{\prime}\prod\limits_{i=1}^{\ell}(a_{i}x+b_{i}y).$

Again, for $z=0$, this equation gives
$-A^{\prime\prime}(x,y,0)(x^{2}+y^{2})=p^{\prime}(x,y,0)\prod\limits_{i=1}^{\ell}(a_{i}x+b_{i}y)$.
If $A^{\prime\prime}(x,y,0)\neq 0$ then $\ell\leq 1$. Otherwise,
$A^{\prime\prime}=\alpha z$ for some $\alpha\in\mathbb{R}$. In that case, we
obtain $\mathcal{E}_{\langle
x,y\rangle}(\mathcal{X})=z^{2}(-\alpha(x^{2}+y^{2})+(C^{\prime}x-B^{\prime}y))$.
So, in this case also $\mathcal{E}_{\langle x,y\rangle}(\mathcal{X})$ has
maximum two factors of the form $ax+by$. So, combining all the cases,
$\mathcal{X}$ has maximum two invariant great circles of the form
$\\{ax+by=0\\}\cap\mathbb{S}^{2}$.

To prove that the bound can be reached, consider
$A^{\prime}=0,B^{\prime}=C^{\prime}=ax+by$. Then the vector field

$\mathcal{X}=(P,Q,R)=((ax+by)z^{2},(ax+by)z^{2},-(ax+by)(x+y)z)$

has invariant hyperplanes $ax+by=0$ and $x-y=0$.

(3) In (4.4), if $\mathcal{X}$ is a homogeneous vector field of degree three,
then $K^{\prime}$ is homogeneous of degree two. Let

$K^{\prime}=k_{1}x^{2}+k_{2}y^{2}+k_{3}z^{2}+k_{4}xy+k_{5}xz+k_{6}yz.$

Then, from (4.4), we get the following by equating the coefficients on its
both sides of each monomial.

$\begin{split}x^{3}:~{}&ba_{1}=-ak_{1}\\\ y^{3}:~{}&aa_{2}=bk_{2}\\\
z^{3}:~{}&ab_{3}+bc_{3}=ck_{3}\\\ x^{2}y:~{}&aa_{1}-ba_{4}=bk_{1}+ak_{4}\\\
x^{2}z:~{}&ba_{5}+cb_{1}=-ck_{1}-ak_{5}\\\
y^{2}z:~{}&aa_{6}-cc_{2}=ck_{2}+bk_{6}\\\
xy^{2}:~{}&aa_{4}-ba_{2}=ak_{2}+bk_{4}\\\
xz^{2}:~{}&-ba_{3}+ab_{1}-cb_{3}+bc_{1}=ak_{3}+ck_{5}\\\
yz^{2}:~{}&aa_{3}+ab_{2}+bc_{2}-cc_{3}=bk_{3}+ck_{6}\\\
xyz:~{}&aa_{5}-ba_{6}-cb_{2}-cc_{1}=ck_{4}+bk_{5}+ak_{6}.\end{split}$

Then, from the first six equations, one can compute $k_{i}$’s. For example
$k_{1}=\frac{ba_{1}}{-a}$, $k_{2}=\frac{aa_{2}}{b}$,
$k_{3}=\frac{ab_{3}+bc_{3}}{c}$, and so on similarly. Substituting these
$k_{i}$’s in the remaining four equations, we get the following four equations
respectively.

(4.6) $\begin{split}E_{1}:~{}&b^{2}a_{1}+a^{2}a_{2}-aba_{4}=0,\\\
E_{2}:~{}&bc^{3}a_{1}+a^{2}bca_{3}-abc^{2}a_{5}-ac(a^{2}+c^{2})b_{1}+a^{2}(a^{2}+c^{2})b_{3}-a^{2}bcc_{1}+a^{3}bc_{3}=0,\\\
E_{3}:~{}&ac^{3}a_{2}+ab^{2}ca_{3}-abc^{2}a_{6}+ab^{2}cb_{2}-ab^{3}b_{3}+bc(b^{2}+c^{2})c_{2}-b^{2}(b^{2}+c^{2})c_{3}=0,\\\
E_{4}:~{}&b^{2}c(a^{2}+2b^{2})a_{1}-a^{4}ca_{2}-ab^{3}ca_{4}-ab^{2}(a^{2}+b^{2})a_{5}+a^{2}b(a^{2}+b^{2})a_{6}-ab^{3}cb_{1}\\\
&+a^{2}b^{2}cb_{2}+a^{2}b^{2}cc_{1}-a^{3}bcc_{2}=0.\end{split}$

If the vector field $\mathcal{X}$ has $\\{ax+by+cz=0\\}\cap\mathbb{S}^{2}$ as
an invariant great circle, then (4.6) must be satisfied. In order for there to
be additional great circles $\\{px+qy+rz=0\\}\cap\mathbb{S}^{2},$ and
$\\{sx+ty+uz=0\\}\cap\mathbb{S}^{2},$ the vector field $\mathcal{X}$ must
satisfy (4.6) with $(p,q,r)$ and $(s,t,u)$ in place of $(a,b,c)$ in addition
to (4.6) itself. Thus we have a system of twelve equations, and solving them
(with the help of SAGEMATH222 Sagemath official website:
https://www.sagemath.org/), we get

$a_{1}=a_{2}=a_{4}=a_{5}=a_{6}=b_{1}=b_{3}=c_{2}=c_{3}=0,\text{and
}a_{3}=-b_{2}=c_{1}.$

For these coefficients, the vector field becomes $\mathcal{X}=(0,0,0)$. Hence,
the maximum number of invariant great circles might be 2.

To prove that the bound can be reached, consider the following example.
Suppose we have the polynomials
$A^{\prime}=x^{2}+y^{2}+2xy+xz+yz,B^{\prime}=-y+z$, $C^{\prime}=x-z$. Then

$\begin{split}\mathcal{X}=&(P,Q,R)\\\
=&(A^{\prime}y+B^{\prime}z^{2},-A^{\prime}x+C^{\prime}z^{2},-B^{\prime}xz-C^{\prime}xz)\\\
=&(x^{2}y+y^{3}+2xy^{2}+xyz+y^{2}z-yz^{2}+z^{3},-x^{3}-xy^{2}-2x^{2}y-x^{2}z-xyz+xz^{2}-z^{3},-xz^{2}+yz^{2})\end{split}$

is a vector field on $\mathbb{S}^{2}$. One can check that $x+y+z=0$ and
$x+y-z=0$ are invariant hyperplanes for $\mathcal{X}$ with corresponding
cofactors $-x^{2}+y^{2}$ and $-x^{2}+y^{2}-2xz+2yz$ respectively. Hence,
$\\{x+y+z=0\\}\cap\mathbb{S}^{2}$ and $\\{x+y-z=0\\}\cap\mathbb{S}^{2}$ are
invariant great circles of $\mathcal{X}$. ∎

Now, we shall look at invariant circles, which are not great circles.

######  Proposition 4.2.

Assume that a cubic homogeneous vector field $\mathcal{X}=(P,Q,R)$ on
$\mathbb{S}^{2}$ has an invariant circle which is not a great circle. Without
loss of generality, assume the invariant circle to be
$\\{z+d=0\\}\cap\mathbb{S}^{2}$ with $0<d<1$. Then
$R=(px+qy)(-d^{2}(x^{2}+y^{2}+z^{2})+z^{2})$ for some $p,q\in\mathbb{R}$.

###### Proof.

Recall that $\mathcal{X}=(P,Q,R)$ is a cubic homogeneous vector field on
$\mathbb{S}^{2}$ if and only if $\mathcal{X}$ is given by (4.2), i.e.,
$P=Ay+Bz,Q=-Ax+Cz,R=-Bx-Cy$ for some polynomials $A,B,C\in\mathbb{R}[x,y,z]$.

The equation of the cone corresponding to the circle
$\mathcal{C}:=\\{z+d=0\\}\cap\mathbb{S}^{2}$ is

$g_{\mathcal{C}}:=-d^{2}(x^{2}+y^{2}+z^{2})+z^{2}.$

Therefore, the condition for the circle $\mathcal{C}$ to be invariant is

(4.7) $\mathcal{X}g_{\mathcal{C}}=K_{\mathcal{C}}g_{\mathcal{C}},$

where the polynomial $K_{\mathcal{C}}$ is the cofactor. Note that
$K_{\mathcal{C}}$ is a homogeneous polynomial. Now,

$\begin{split}&\mathcal{X}g_{\mathcal{C}}=K_{\mathcal{C}}g_{\mathcal{C}}\\\
\implies&-2d^{2}(Px+Qy+Rz)+2Rz=K_{\mathcal{C}}g_{\mathcal{C}}\\\
\implies&2Rz=K_{\mathcal{C}}g_{\mathcal{C}}\quad\text{since
$Px+Qy+Rz=0$}.\end{split}$

Since $d\neq 0$, $z$ does not divide $g_{\mathcal{C}}$. So, the last
implication tells us that $z$ divides $K_{\mathcal{C}}$, assume that
$K_{\mathcal{C}}=2z(px+qy+rz)$ for some $p,q,r\in\mathbb{R}$. Hence,
$R=(px+qy+rz)(-d^{2}(x^{2}+y^{2}+z^{2})+z^{2})$. Since $R=-Bx-Cy$, we get
$r(1-d^{2})=0$. So, either $r=0$ or $d=\pm 1$. Since $d<1$, $r$ must be zero.
Then $R=(px+qy)(-d^{2}(x^{2}+y^{2}+z^{2})+z^{2})$. ∎

######  Remark 4.3.

Proposition 4.2 can also be proved for degree $n(>3)$ homogeneous vector
fields on $\mathbb{S}^{2}$. In that case,
$R=(Mx+Ny)(-d^{2}(x^{2}+y^{2}+z^{2})+z^{2})$ for some
$M,N\in\mathbb{R}[x,y,z]$ such that $\deg M=(n-3)$ or $\deg N=(n-3)$.

##  5\. Cubic Kolmogorov vector fields on $\mathbb{S}^{2}$

We recall that the form of a cubic Kolmogorov vector field on $\mathbb{S}^{2}$
has been examined in Corollary 3.3. This section aims to study the dynamical
properties of these cubic Kolmogorov vector fields on $\mathbb{S}^{2}.$ For
convenience, we shall rewrite the general form of a cubic Kolmogorov vector
field, $\mathcal{X}=(P,Q,R)$ defined on $\mathbb{S}^{2},$

(5.1) $\begin{split}P&=x(\alpha(1-x^{2}-y^{2}-z^{2})+Ay^{2}+Bz^{2}),\\\
Q&=y(\beta(1-x^{2}-y^{2}-z^{2})-Ax^{2}+Cz^{2}),~{}\mbox{and}\\\
R&=z(\gamma(1-x^{2}-y^{2}-z^{2})-Bx^{2}-Cy^{2}),\end{split}$

where $\alpha,\beta,\gamma$ and $A,B,C$ are constants. In what follows, we
shall assume that $A,B,C\neq 0$ unless it is specified. Notice from (5.1) that
the hyperplanes $\\{x=0\\}$, $\\{y=0\\}$, and $\\{z=0\\}$ are invariant with
respect to the flow of $\mathcal{X}.$ This gives that on $\mathbb{S}^{2},$ the
great circles determined by the intersection of these hyperplanes with
$\mathbb{S}^{2},$ that is, $\\{x=0\\}\cap\mathbb{S}^{2}$,
$\\{y=0\\}\cap\mathbb{S}^{2}$, and $\\{z=0\\}\cap\mathbb{S}^{2}$ are invariant
with respect to the flow of $\mathcal{X}.$ From this, it implies that the
points where these great circles intersect are singular points. Hence $(\pm
1,0,0),$ $(0,\pm 1,0),$ and $(0,0,\pm 1)$ are singular points of
$\mathcal{X}.$ We want to determine if $\mathcal{X}$ has any other singular
points on $\mathbb{S}^{2}.$ Note that from (5.1), additional singular points
can be obtained by solving the following equations.

$\begin{split}Ay^{2}+Bz^{2}=0,\\\ -Ax^{2}+Cz^{2}=0,\\\
~{}~{}-Bx^{2}-Cy^{2}=0.\end{split}$

These equations, together with the fact that we are on the sphere,
$x^{2}+y^{2}+z^{2}=1,$ give

(5.2) $x=\pm\sqrt{\frac{-C}{B-A-C}},\quad
y=\pm\sqrt{\frac{B}{B-A-C}},\quad\mbox{and}\quad
z=\pm\sqrt{\frac{-A}{B-A-C}}.$

Observe that $(x,y,z)$ is a point in $\mathbb{R}^{3}$ when $A,C>0,$ $B<0$ or
when $A,C<0,$ $B>0.$ If neither of these conditions are true, $\mathcal{X}$
has only $(\pm 1,0,0),$ $(0,\pm 1,0),$ and $(0,0,\pm 1)$ as singularities.

Our objective is to determine the phase portraits of cubic Kolmogorov vector
fields on $\mathbb{S}^{2}.$ From (5.1), we notice that the antipodal
transformation $(x,y,z)\mapsto(-x,-y,-z)$ on $\mathbb{R}^{3}$ causes
$(P,Q,R)\mapsto(-P,-Q,-R)$. This implies that trajectories in the southern
hemisphere mirror trajectories in the northern hemisphere with the direction
of flow reversed. Thus, in order to determine the phase portrait of these
vector fields on $\mathbb{S}^{2}$, it is sufficient to determine the phase
portrait on one of the hemispheres (including its boundary). In this section,
we shall use the stereographic projection $\Phi$ of Section 2 from the South
pole, to determine the phase portrait on the northern hemisphere. The northern
hemisphere along with its boundary is mapped to the closed unit disk
$D^{2}:=\\{(u,v)\in\mathbb{R}^{2}~{}|~{}u^{2}+v^{2}\leq 1\\}$ by $\Phi.$ We
find the converted vector field in $\mathbb{R}^{2}$ of the cubic Kolmogorov
vector field on $\mathbb{S}^{2}$ given by (5.1). For that, we have the
following:

$\begin{split}\tilde{P}&=2u(4Av^{2}+B(1-u^{2}-v^{2})^{2}),\\\
\tilde{Q}&=2v(-4Au^{2}+C(1-u^{2}-v^{2})^{2}),~{}\mbox{and}\\\
\tilde{R}&=-4(1-u^{2}-v^{2})(Bu^{2}+Cv^{2}).\end{split}$

This gives

(5.3)
$\begin{split}&\mathcal{P}=\dot{u}=2u(4Av^{2}+B(1-u^{2}-v^{2})^{2})-u\tilde{R},~{}\mbox{and}\\\
&\mathcal{Q}=\dot{v}=2v(-4Au^{2}+C(1-u^{2}-v^{2})^{2})-v\tilde{R}.\end{split}$

We shall now study the behaviour of this system near the singular points
$(1,0,0),$ $(0,1,0),$ and $(0,0,1).$ Note that the behaviour of the system
near the singular points $(-1,0,0),$ $(0,-1,0),$ and $(0,0,-1)$ is identical
(with reverse arrows) to the behaviour around their antipodal counterparts by
the discussion above.

In the case(s) where these are the only singular points of $\mathcal{X},$ we
shall be able to draw phase portraits. We shall compute the Jacobian, $J$ at
singular points of the vector field
$\Phi_{*}(\mathcal{X})=(\mathcal{P},\mathcal{Q})$ in $\mathbb{R}^{2}.$ We have

(5.4) $J=\begin{pmatrix}\mathcal{P}_{u}&\mathcal{P}_{v}\\\
\mathcal{Q}_{u}&\mathcal{Q}_{v}\end{pmatrix}.$

From (5.3), we have

$\begin{split}&\mathcal{P}_{u}=2(4Av^{2}+B(1-u^{2}-v^{2})^{2})-8Bu^{2}(1-u^{2}-v^{2})-\tilde{R}-u\tilde{R}_{u},\\\
&\mathcal{P}_{v}=2u(8Av-4Bv(1-u^{2}-v^{2}))-u\tilde{R}_{v},\\\
&\mathcal{Q}_{u}=2v(-8Au-4Cu(1-u^{2}-v^{2}))-v\tilde{R}_{u},\\\
&\mathcal{Q}_{v}=2(-4Au^{2}+C(1-u^{2}-v^{2})^{2})-8Cv^{2}(1-u^{2}-v^{2})-\tilde{R}-v\tilde{R}_{v}.\end{split}$

Also,

$\begin{split}&\tilde{R}_{u}=8u(Bu^{2}+Cv^{2})-8Bu(1-u^{2}-v^{2}),\\\
&\tilde{R}_{v}=8v(Bu^{2}+Cv^{2})-8Cv(1-u^{2}-v^{2}).\end{split}$

Now we compute the Jacobian for the singular points $(1,0,0),$ $(0,1,0),$ and
$(0,0,1).$

  1. (1)

For $(1,0,0)$: The stereographic projection $\Phi$ maps this point to $(1,0)$
on the plane. We compute

$\mathcal{P}_{u}(1,0)=-8B,\quad\mathcal{P}_{v}(1,0)=0,\quad\mathcal{Q}_{u}(1,0)=0,\quad\mathcal{Q}_{v}(1,0)=-8A.$

Therefore,

$J_{(1,0)}=8\begin{pmatrix}-B&0\\\ 0&-A\end{pmatrix}.$

  2. (2)

For $(0,1,0)$: The stereographic projection $\Phi$ maps this point to $(0,1)$
on the plane. We compute

$\mathcal{P}_{u}(0,1)=8A,\quad\mathcal{P}_{v}(0,1)=0,\quad\mathcal{Q}_{u}(0,1)=0,\quad\mathcal{Q}_{v}(0,1)=-8C.$

Therefore,

$J_{(0,1)}=8\begin{pmatrix}A&0\\\ 0&-C\end{pmatrix}.$

  3. (3)

For $(0,0,1)$: The stereographic projection $\Phi$ maps this point to $(0,0)$
on the plane. We compute

$\mathcal{P}_{u}(0,0)=2B,\quad\mathcal{P}_{v}(0,0)=0,\quad\mathcal{Q}_{u}(0,0)=0,\quad\mathcal{Q}_{v}(0,0)=2C.$

Therefore,

$J_{(0,0)}=2\begin{pmatrix}B&0\\\ 0&C\end{pmatrix}.$

Assume that $A,B,C\neq 0.$ Let

(5.5) $(a)~{}~{}A,C>0,B<0,\quad\mbox{and}\quad(b)~{}~{}A,C<0,B>0.$

Thus, if neither (a) nor (b) is true, then $(\pm 1,0,0),$ $(0,\pm 1,0),$ and
$(0,0,\pm 1)$ are the only singular points of $\mathcal{X}$ on
$\mathbb{S}^{2}$, and we draw the phase portraits for each of these cases in
Figure 1.

i When $A,B,C>0$.

ii When $A<0,$ $B,C>0$.

iii When $A,B>0,$ $C<0$.

iv When $A,B,C<0$.

v When $A>0,$ $B,C<0$.

vi When $A,B<0,$ $C>0$.

Figure 1. Phase portraits when neither (a) nor (b) of (5.5) is true.

######  Remark 5.1.

When one of $A,B,C$ is zero, then the singular points are no longer isolated.
For example, we plot the phase portrait when $A=0$ and $B,C>0$ in Figure 2. In
this case, all points on the boundary of the disk (i.e., the unit circle) are
singular.

Figure 2.

From the above discussion, we see that when neither (a) nor (b) is true in
(5.5), any Kolmogorov vector field of the form in (5.1) has phase portrait
that is topologically equivalent to either (I) or (III) in Figure 1. So we
have the following theorem.

######  Theorem 5.2.

A cubic Kolmogorov vector field on $\mathbb{S}^{2}$ of the form given in (5.1)
admits no periodic orbit if $A,B,C$ in (5.1) do not satisfy the following.

$(1)~{}A,C>0,B<0,\quad\mbox{and}\quad(2)~{}A,C<0,B>0.$

In the remaining, we discuss what happens when condition (a) in (5.5) is true.
Note that the same analysis will also apply for (b), albeit with a change in
the sense of the orbits in the phase portrait. That is when $A,C>0$ and $B<0.$
We know that in this case, there exist extra singularities given in (5.2).
Applying the stereographic projection to these singular points, we have,

$u=\frac{\pm\sqrt{-C}}{\sqrt{B-A-C}\pm\sqrt{-A}}\quad,\quad
v=\frac{\pm\sqrt{B}}{\sqrt{B-A-C}\pm\sqrt{-A}}.$

Let us fix a singular point

$u_{0}=\frac{\sqrt{-C}}{\sqrt{B-A-C}+\sqrt{-A}}\quad,\quad
v_{0}=\frac{\sqrt{B}}{\sqrt{B-A-C}+\sqrt{-A}}.$

We note that the analysis below can be carried out entirely similarly for the
other singular points. We need to evaluate the Jacobian in (5.4) at
$(u_{0},v_{0}).$ After the calculation, we get the following.

$\begin{split}\mathcal{P}_{u}(u_{0},v_{0})&=\frac{8AB}{D^{2}}+2B\left(\frac{D^{2}+C-B}{D^{2}}\right)^{2},\\\
\mathcal{P}_{v}(u_{0},v_{0})&=\frac{16A\sqrt{-BC}}{D^{2}}+8\frac{(C-B)\sqrt{-BC}}{D^{2}}\left(\frac{D^{2}+C-B}{D^{2}}\right),\\\
\mathcal{Q}_{u}(u_{0},v_{0})&=\frac{-16A\sqrt{-BC}}{D^{2}}-8\frac{(C-B)\sqrt{-BC}}{D^{2}}\left(\frac{D^{2}+C-B}{D^{2}}\right),\\\
\mathcal{Q}_{v}(u_{0},v_{0})&=\frac{8AC}{D^{2}}+2C\left(\frac{D^{2}+C-B}{D^{2}}\right)^{2},\end{split}$

where $D:=\sqrt{B-A-C}+\sqrt{-A}.$ Notice that with our choice of the
coefficients $A,B$ and $C$, the term $D$ is imaginary and so $D^{2}<0.$ We
take

$F=\frac{8A}{D^{2}}+2\left(\frac{D^{2}+C-B}{D^{2}}\right)^{2},$

then the characteristic polynomial $c(\lambda)$ for the Jacobian matrix
$J_{(u_{0},v_{0})}$ is

(5.6)
$c(\lambda)=\lambda^{2}-(B+C)F\lambda+BCF^{2}+\mathcal{P}_{v}^{2}(u_{0},v_{0}).$

The roots of this polynomial are by definition the eigenvalues of the
Jacobian, and from these eigenvalues, we shall determine the nature of the
singular point. The discriminant, $\Delta$, of $c(\lambda)$ is

$\begin{split}\Delta:=&(B+C)^{2}F^{2}-4BCF^{2}-4\mathcal{P}_{v}^{2}(u_{0},v_{0})\\\
=&(C-B)^{2}F^{2}-4\mathcal{P}_{v}^{2}(u_{0},v_{0})\\\
=&((C-B)F+2\mathcal{P}_{v}(u_{0},v_{0}))((C-B)F-2\mathcal{P}_{v}(u_{0},v_{0})).\end{split}$

Now,

$\begin{split}&(C-B)F+2\mathcal{P}_{v}(u_{0},v_{0})\\\
=&\frac{8A(C-B)}{D^{2}}+2(C-B)\left(\frac{D^{2}+C-B}{D^{2}}\right)^{2}+\frac{32A\sqrt{-BC}}{D^{2}}+\frac{16(C-B)\sqrt{-BC}}{D^{2}}\left(\frac{D^{2}+C-B}{D^{2}}\right)\\\
=&\frac{8A(C-B+4\sqrt{-BC})}{D^{2}}+\frac{2(C-B)(D^{2}+C-B+8\sqrt{-BC})}{D^{2}}\left(\frac{D^{2}+C-B}{D^{2}}\right).\end{split}$

For large enough $A,$ the expression $(C-B)F+2\mathcal{P}_{v}(u_{0},v_{0})$ is
negative. In fact

$\begin{split}&\lim_{A\rightarrow\infty}\\{(C-B)F+2\mathcal{P}_{v}(u_{0},v_{0})\\}\\\
=&\lim_{A\rightarrow\infty}\left(\frac{8A[(\sqrt{C}+\sqrt{-B})^{2}+2\sqrt{-BC}]}{D^{2}}+\frac{2(C-B)[D^{2}+(\sqrt{C}+\sqrt{-B})^{2}+6\sqrt{-BC}]}{D^{2}}\left(\frac{D^{2}+C-B}{D^{2}}\right)\right)\\\
=&\lim_{A\rightarrow\infty}\left(\frac{8A[(C-B+4\sqrt{-BC}]}{B-2A-C+2\sqrt{-A}\sqrt{B-A-C}}+\frac{2(C-B)[-2A+2\sqrt{-A}\sqrt{B-A-C}+8\sqrt{-BC}]}{B-2A-C+2\sqrt{-A}\sqrt{B-A-C}}\left(\frac{-2A+2\sqrt{-A}\sqrt{B-A-C}}{B-2A-C+2\sqrt{-A}\sqrt{B-A-C}}\right)\right)\\\
=&-8\sqrt{-BC}.\end{split}$

In a similar manner

$\begin{split}&(C-B)F-2\mathcal{P}_{v}(u_{0},v_{0})\\\
=&\frac{8A(C-B)}{D^{2}}+2(C-B)\left(\frac{D^{2}+C-B}{D^{2}}\right)^{2}-\frac{32A\sqrt{-BC}}{D^{2}}-\frac{16(C-B)\sqrt{-BC}}{D^{2}}\left(\frac{D^{2}+C-B}{D^{2}}\right)\\\
=&\frac{8A(C-B-4\sqrt{-BC})}{D^{2}}+\frac{2(C-B)(D^{2}+C-B-8\sqrt{-BC})}{D^{2}}\left(\frac{D^{2}+C-B}{D^{2}}\right).\end{split}$

For large enough $A,$ the expression $(C-B)F-2\mathcal{P}_{v}$ is positive. In
fact

$\lim_{A\rightarrow\infty}\\{(C-B)F-2\mathcal{P}_{v}\\}=8\sqrt{-BC}.$

From the above, we conclude that for $A$ large enough, the discriminant
$\Delta$ of $c(\lambda)$ is negative and hence the Jacobian matrix at
$(u_{0},v_{0})$ has a pair of complex conjugate eigenvalues, say $\mu\pm
i\nu$. Thus we get the following result using [17, Theorem 4, Page-143] and
[17, Corollary to Theorem 5, Page-145].

######  Theorem 5.3.

For $A$ large enough, the singular point $(u_{0},v_{0})$ is either a center or
a focus.

######  Remark 5.4.

If the real part of the eigenvalue $\mu>0,$ then the singular point
$(u_{0},v_{0})$ is an unstable focus and if $\mu<0,$ then the singular point
$(u_{0},v_{0})$ is a stable focus, see [17, Theorem 4, Page-143].

######  Example 5.5.

We compute the discriminant and also determine the nature of the extra
singularities when $A=5,$ $B=-1$ and $C=2.$ In this case, the computation
gives that $(B+C)F=-0.0001$ and $\Delta=-123.54.$

From (5.6), we see that the eigenvalues are given by

$\frac{(B+C)F\pm\sqrt{\Delta}}{2}.$

Therefore, the eigenvalues are complex conjugate with a negative real part.
This implies that the singular point is a stable focus. A similar calculation
will also show that the remaining singularities belonging to the interior of
$D^{2}$ are also stable foci. In Figure 3, we plot the phase portrait of this
case.

Figure 3. When $A=5$ and $B=-1$, and $C=2$.

Acknowledgments. The first author was supported by a Senior Research
Fellowship of the University Grants Commission of India for the duration of
this work. The second author is supported by the Prime Minister’s Research
Fellowship, Government of India. The third author thanks ‘ICSR office IIT
Madras’ for SEED research grant.

## References

  * [1] L. Apedaille, H. Freedman, S. Schilizzi, and M. Solomonovich.  Equilibria and dynamics in an economic predator-prey model of agriculture.  Mathematical and Computer Modelling, 19(11):1–15, 1994.
  * [2] C. Christopher, J. Llibre, and J. V. Pereira.  Multiplicity of invariant algebraic curves in polynomial vector fields.  Pacific Journal of Mathematics, 229(1):63–117, 2007.
  * [3] F. Dumortier, J. Llibre, and J. C. Artés.  Qualitative theory of planar differential systems, volume 2.  Springer, 2006.
  * [4] A. Gasull and R. Prohens.  Quadratic and cubic systems with degenerate infinity.  J. Math. Anal. Appl., 198(1):25–34, 1996.
  * [5] S. Jana and S. Sarkar.  Quadratic, homogeneous and Kolmogorov vector fields on $S^{1}\times S^{2}$ and $S^{2}\times S^{1}$.  arXiv preprint arXiv:2307.09439, 2023.
  * [6] A. Kolmogorov.  Sulla teoria di Volterra della lotta per lesistenza.  Gi. Inst. Ital. Attuari, 7:74–80, 1936.
  * [7] M. Leconte, P. Masson, and L. Qi.  Limit cycle oscillations, response time, and the time-dependent solution to the Lotka-Volterra predator-prey model.  Physics of Plasmas, 29(2):022302, 02 2022.
  * [8] J. Llibre.  Open problems on the algebraic limit cycles of planar polynomial vector fields.  Bul. Acad. Ştiinţe Repub. Mold. Mat., (1):19–26, 2008.
  * [9] J. Llibre and Y. Bolaños.  Rational first integrals for polynomial vector fields on algebraic hypersurfaces of $\mathbb{R}^{n+1}$.  Internat. J. Bifur. Chaos Appl. Sci. Engrg., 22(11):1250270, 11, 2012.
  * [10] J. Llibre and J. C. Medrado.  On the invariant hyperplanes for $d$-dimensional polynomial vector fields.  Journal of Physics A: Mathematical and Theoretical, 40(29):8385, 2007.
  * [11] J. Llibre and C. Pessoa.  Homogeneous polynomial vector fields of degree 2 on the 2-dimensional sphere.  Extracta Math., 21(2):167–190, 2006.
  * [12] J. Llibre, R. Ramírez, and V. Ramírez.  Integrability of a class of $N$-dimensional Lotka–Volterra and Kolmogorov systems.  Journal of differential equations, 269(3):2503–2531, 2020.
  * [13] J. Llibre, R. Ramírez, and V. Ramírez.  Complete integrability of vector fields in $\mathbb{R}^{N}$.  Journal of Geometry and Physics, 168:104308, 2021.
  * [14] J. Llibre and D. Xiao.  Dynamics, integrability and topology for some classes of kolmogorov hamiltonian systems in $\mathbb{R}_{+}^{4}$.  Journal of Differential Equations, 262(3):2231–2253, 2017.
  * [15] J. Llibre and X. Zhang.  Darboux theory of integrability in $\mathbb{C}^{n}$ taking into account the multiplicity.  Journal of Differential Equations, 246(2):541–551, 2009.
  * [16] J. Llibre and Y. Zhao.  On the polynomial vector fields on $\mathbb{S}^{2}$.  Proc. Roy. Soc. Edinburgh Sect. A, 141(5):1055–1069, 2011.
  * [17] L. Perko.  Differential equations and dynamical systems, volume 7 of Texts in Applied Mathematics.  Springer-Verlag, New York, third edition, 2001.


# Laser-painted cavity-mediated interactions in a quantum gas

Mariano Bonifacio  Max-Planck-Institut für Physik komplexer Systeme, 01187
Dresden, Germany Francesco Piazza<EMAIL_ADDRESS>Max-Planck-Institut für
Physik komplexer Systeme, 01187 Dresden, Germany Theoretical Physics III,
Center for Electronic Correlations and Magnetism, Institute of Physics,
University of Augsburg, 86135 Augsburg, Germany Tobias Donner  Institute for
Quantum Electronics, Eidgenössische Technische Hochschule Zürich, Otto-Stern-
Weg 1, CH-8093 Zurich, Switzerland

###### Abstract

Experimental platforms based on ultracold atomic gases have significantly
advanced the quantum simulation of complex systems, yet the exploration of
phenomena driven by long-range interactions remains a formidable challenge.
Currently available methods utilizing dipolar quantum gases or multi-mode
cavities allow to implement long-range interactions with a $1/r^{3}$ character
or with a spatial profile fixed by the mode-structure of the vacuum
electromagnetic field surrounding the atoms, respectively. Here we propose an
experimental scheme employing laser-painted cavity-mediated interactions,
which enables the realization of atom-atom interactions that are fully tunable
in range, shape, and sign. Our approach combines the versatility of cavity
quantum electrodynamics with the precision of laser manipulation, thus
providing a highly flexible platform for simulating and understanding long-
range interactions in quantum many-body systems. Our analytical predictions
are supported by numerical simulations describing the full dynamics of atoms,
laser, and cavity. The latter demonstrate that there is a wide and
experimentally accessible parameter regime where our protocol optimally works.
The methodology not only paves the way for exploring new territories in
quantum simulation but also enhances the understanding of fundamental physics,
potentially leading to the discovery of novel quantum states and phases.

Quantum simulators hold the potential to unravel intricate quantum states and
dynamics of many-body systems that are challenging to analyze using
traditional methods and impractical for classical numerical simulations [1].
Specifically ultracold atoms have been established as a powerful experimental
platform for simulating models ranging from condensed matter [2, 3, 4] over
high-energy physics [5] to cosmology [6]. This success is based on the almost
perfect control over internal and external degrees of freedom of the
particles, paired with their fully tunable local interactions provided by
elastic atomic collisions.

A large number of phenomena in nature however arises due to long-range
interactions between its constituents or due to competing long- and short-
range interactions, and are thus out of reach for the traditional tool set of
simulations using quantum gases [7]. For example the formation of crystalline
or glassy phases [8, 9, 10, 11, 12, 13, 14, 15], the existence of exotic many-
body phenomena such as topological or frustrated magnetism [16, 17, 18, 19],
supersolidity [20, 21, 22, 23, 24] and roton excitations[25, 26], or quantum
droplets [27, 28, 29] usually require the presence of long-range interactions.
Not only interactions with a simple power-law decay $1/r^{\alpha}$ as a
function of distance $r$ between the particles are of interest, with for
example $\alpha=1$ for Coulomb interactions and $\alpha=3$ for dipolar
interactions, but specifically also more complex interactions such as the
Ruderman-Kittel-Kasuya-Yosida (RKKY)-type interaction [30, 31, 32, 33] which
oscillates in space and is a fundamental model in condensed matter physics
also applied to describe giant magnetoresistance.

Recently, tunable spin-exchange interactions have been implemented in ion
chains and in coupled thermal atomic ensembles [34, 35]. However, approaches
to engineer long-range interactions in quantum gases usually involve atoms
with large magnetic dipole moments [36, 37], polar molecules [38], and Rydberg
atoms [39], which all lead to a $1/r^{3}$ scaling of the interactions due to
the underlying dipole-dipole coupling. An alternative approach to engineer
long-range interactions in quantum gases makes use of hybrid systems combining
quantum gases with high-finesse optical cavities, where photons from an
external drive effectively mediate interactions between the involved atoms
[40]. Microscopically, a photon from the drive field impinging on the atomic
cloud is scattered into the cavity mode at a first atom and then back into the
drive at a second atom, mediating the effective interactions. For a single-
mode cavity, where the photon is delocalized over the entire cavity mode, this
results in global-range interactions which have been extensively exploited to
realize lattice supersolids [41, 42], quantum optics spin models [43], or
extended Hubbard-models [44, 45]. Coupling simultaneously to two cavity modes
allowed the realization of a supersolid breaking a continuous translational
symmetry [21], and extending the coupling to many degenerate modes enabled the
realization of finite-range [46] and sign-changing interactions [47], and the
simulation of phonons in a lattice model [48]. However, while the cavity-
mediated interaction range becomes indeed tunable to some degree, their shape
is still dictated by the involved cavity mode patterns and thus fixed by a
given choice of cavity geometry.

In this article, we introduce a novel experimental scheme based on laser-
painted cavity-mediated interactions, enabling the realization of atom-atom
interactions fully tunable in range, shape, and sign. We detail the basic
setup and establish a comprehensive framework for the effective atom-atom
interaction, encompassing both cavity-mediated long-range and collisional
short-range interactions. Through numerical simulations in 1D, we illustrate
the practicality, resilience, and boundaries of our proposed methodology. Our
approach includes an integrated tomographic feature for the in-situ analysis
of the atomic density profile through the cavity’s output. The proposed scheme
offers a gateway to explore the diverse realms of long-range interaction
physics, not only making quantum simulations of purely theoretical models
feasible but also opening doors to uncharted phenomenological exploration.

## Results

### Experimental scheme

The experimental concept is visualized in Figure 1. A Bose-Einstein condensate
(BEC) is dispersively coupled with single-atom vacuum Rabi rate $g_{0}$ to a
single optical cavity mode with resonance frequency $\omega_{c}$, wavelength
$\lambda_{c}$, and field decay rate $\kappa$. The oblate BEC with lateral size
$L$ is confined in a plane orthogonal to the cavity axis such that it resides
at an anti-node of the cavity mode and the coupling to the cavity is
approximately constant across the entire cloud. Also schemes using a different
geometry are feasible as we demonstrate with our numerical simulations. The
BEC is illuminated by a focused laser field with waist $w\ll L$ and local Rabi
frequency $\Omega(\vec{r},t)$. This laser spot is swept as a function of time
across the BEC following a path $\vec{r}_{0}(t)$ along horizontal lines
separated by a distance $l$.

Concentrating for the moment only on the region illuminated by the laser beam
at a certain moment in time, effective interactions between two atoms are
generated by a cavity-mediated dipole-dipole coupling [40]. Since the
frequency $\omega_{L}$ of the laser is chosen far red detuned by
$\Delta_{a}=\omega_{a}-\omega_{L}\gg\Gamma$ with respect to the atomic
resonance frequency $\omega_{a}$, where $\Gamma$ is the decay rate of the
atomic transition, the atoms are not electronically excited, but an
oscillating electric dipole moment is induced in each atom. A free-space
coupling between two induced dipoles had been suggested as a mechanism to
generate long-range interactions [49, 50]. These direct dipole-dipole
interactions decay with $1/r^{3}$, and in addition exhibit the same
$(\Gamma/\Delta_{a})^{2}$ scaling as spontaneous scattering [51]. As a result,
interesting quantum phases induced by this interaction are endangered by
heating [52, 53]. The presence of a strongly coupled cavity mode however
boosts the interaction between the induced dipoles by the finesse of the
cavity [54]: The field radiated by an atomic dipole will predominantly drive
the cavity mode if the detuning $\Delta_{c}=\omega_{c}-\omega_{L}$ is
sufficiently small. The cavity field then will interfere with the driving
laser field at the position of the other atomic dipole, mediating an
interaction between the two. Only atoms simultaneously coupled to the cavity
and illuminated by the laser spot will participate in this coupling, since the
cavity field alone is too small to significantly influence the dipole moment
of an atom without amplification by the interference term.

Figure 1: Schematic of the cavity-QED experimental setup for the
implementation of finite-range light-mediated interactions between atoms. A
Bose-Einstein condensate (grey, not to scale) is strongly coupled to an
optical cavity mode (light red mode confined by blue mirrors) with resonance
frequency $\omega_{c}$, and tightly confined in a plane perpendicular to the
axis of the cavity. The atomic cloud is driven by a focused laser field with
mode function $\Omega(\vec{r},t)$ and frequency $\omega_{L}$ that scans across
the atomic cloud periodically following a path of horizontal lines. Photons
leaking from the cavity mode at rate $\kappa$ provide real-time information on
the atomic system. The inset shows the microscopic process of light-mediated
interactions: Only atoms within the spot illuminated by the transverse laser
field interact with each other via the exchange of photons. The interaction
potential is determined by the shape of the laser spot. Time-averaging of the
process results in a effective atomic interaction in the entire cloud.

This process induces for $\Delta_{a},\Delta_{c}<0$ attractive “local”
interactions between atoms at the position of the laser spot, and repulsive
ones for $\Delta_{c}>0$. The range and shape of these interactions is set by
the field distribution $\Omega(\vec{r},t)$ of the laser mode, which can be
freely engineered in the experiment using beam shaping techniques. Even far
distant atoms can be made to interact if two regions of the BEC are
simultaneously illuminated.

In order to generate a translational invariant interaction between the atoms
in the entire BEC, the focused beam is rapidly scanned across the cloud at
rate $\omega_{\mathrm{scan}}$. Analogous to the creation of time-averaged
optical potentials [55], this process generates a time-averaged interaction
between the atoms. Two time scales are important when considering this scheme:
Atomic density patterns emerge and decay on a time scale given by the recoil
frequency $\omega_{R}=2\pi^{2}/M\lambda_{c}^{2}$, such that the beam has to be
scanned across the cloud much faster than this time scale, with $M$ being the
atomic mass. On the other hand, the cavity decay rate $\kappa$ determines if
the field scattered into the cavity by the atoms illuminated at one position
has decayed before the spot has moved to the next position. This is required
such that the interaction is only induced where the scanning laser spot
resides. If the beam scans across the cloud before the field decays, the
system would behave the same as a globally illuminated cloud, i.e. giving rise
to infinite-range interactions. Assuming typical values of $\omega_{R}\approx$
kHz and $\kappa\approx$ MHz [40], leaves a broad frequency window for the
scanning frequency $\omega_{\mathrm{scan}}$.

### Effective atom-atom interaction

We consider the situation sketched in Figure 1 with $N$ atoms coupled to a
single mode cavity and subject to a scanning laser field that is far red
detuned from atomic resonance, $\Delta_{a}<0$. With $\Delta_{a}$ being the
largest energy scale, the excited state of the atoms can be adiabatically
eliminated and the dynamics of the system is well described in mean-field
approximation by the Gross-Pitaevskii equation for the condensate wavefunction
$\psi(\vec{r},t)$, normalized to $N$ [40, 56, 57]

$\begin{split}i\frac{\partial{\psi}(\vec{r},t)}{\partial
t}=\Big{[}-\frac{\hbar\nabla^{2}}{2M}+\frac{g_{\text{aa}}}{\hbar}|{\psi}(\vec{r},t)|^{2}+\frac{\Omega^{2}(\vec{r},t)}{\Delta_{a}}+&\\\
+\frac{g_{0}^{2}}{\Delta_{a}}|\alpha(t)|^{2}+\frac{g_{0}\Omega(\vec{r},t)}{\Delta_{a}}(\alpha(t)+\alpha^{*}(t))\Big{]}{\psi}(\vec{r},t)&\end{split}$
(1)

which is coupled to the equation for the evolution of the coherent field
amplitude $\alpha(t)$ of the cavity

$\begin{split}i\frac{\partial\alpha(t)}{\partial
t}=&(-\Delta_{c}-i\kappa)\alpha(t)+\\\ &+\frac{1}{\Delta_{a}}\int
g_{0}[g_{0}\alpha(t)+\Omega(\vec{r},t)]|{\psi}(\vec{r},t)|^{2}d\vec{r}\,,\end{split}$
(2)

where $M$ is the atomic mass, $g_{\text{aa}}$ is the effective contact
interaction strength between atoms associated with the s-wave scattering and
${\Omega(\vec{r},t)=\tilde{\Omega}(\vec{r}-\vec{r}_{0}(t))}$ is the mode
function of the laser beam.

We assume the characteristic time scale $1/\kappa$ of the cavity field to be
fast with respect to the atomic motion and adiabatically eliminate the cavity
field (see Methods). In addition, we consider a regime of parameters in which
the movement of the laser is faster than the atoms’ dynamics, this is when the
angular frequency of the laser sweep is large compared to both the recoil
frequency and the contact interaction
${\omega_{\mathrm{scan}}\gg\omega_{R},\hbar^{-1}g_{\text{aa}}|\psi|^{2}}$. In
this situation one can take the average over one period of the laser
displacement $T=2\pi/\omega_{\mathrm{scan}}$, which results in the equation of
motion

$\begin{split}i\frac{\partial{\psi}(\vec{r},t)}{\partial
t}=\Big{[}-\frac{\hbar\nabla^{2}}{2M}+\frac{g_{\text{aa}}}{\hbar}|{\psi}(\vec{r},t)|^{2}+\frac{\overline{\Omega^{2}(\vec{r},t)}}{\Delta_{a}}+&\\\
+\int
V(\vec{r},\vec{r}\,^{\prime})|{\psi}(\vec{r}\,^{\prime},t)|^{2}d\vec{r}\,^{\prime}\Big{]}{\psi}(\vec{r},t)&\end{split}$
(3)

with effective interaction

$V(\vec{r},\vec{r}\,^{\prime})=\frac{2\Delta_{c}g_{0}^{2}}{\Delta_{a}^{2}(\kappa^{2}+\Delta_{c}^{2})}\overline{\Omega(\vec{r},t)\Omega(\vec{r}\,^{\prime},t)}\,,$
(4)

where

$\overline{\Omega(\vec{r},t)\Omega(\vec{r}\,^{\prime},t)}=\frac{1}{T}\int_{0}^{T}dt\tilde{\Omega}(\vec{r}-\vec{r}_{0}(t))\tilde{\Omega}(\vec{r}\,^{\prime}-\vec{r}_{0}(t))\,.$
(5)

The term $\frac{\overline{\Omega^{2}(\vec{r},t)}}{\Delta_{a}}$ in Equation (3)
is a constant position-independent potential for the atoms and does not
contribute to the dynamics. For simplicity, we assume a centrosymmetric laser
spot, $\tilde{\Omega}(-(x,y))=\tilde{\Omega}((x,y))$, for which the above
expression takes the form of a convolution
${\frac{1}{L^{2}}(\tilde{\Omega}\ast\tilde{\Omega})(\vec{r}-\vec{r}\,^{\prime})}$,
see Methods. As a result, the atom-atom interaction in momentum space is
determined by the square of the Fourier transform of the shape of the laser
spot,

$\hat{V}(\vec{k})=\frac{2\Delta_{c}g_{0}^{2}}{\Delta_{a}^{2}(\kappa^{2}+\Delta_{c}^{2})}\frac{1}{L^{2}}\hat{\tilde{\Omega}}(\vec{k})^{2}\,.$
(6)

With the current experimental state of the art, this interaction potential can
be essentially designed at will: Rapid scanning of a laser beam is routinely
achieved using galvo-mirrors, polygon scanners, or acousto-optical deflectors,
while arbitrary and even time dependent beam shapes can be produced using
digital mirror devices (DMDs). DMDs are specifically suited if placed in the
Fourier plane of the optics projecting the laser spot onto the atoms, such
that the interaction potential, Equation (6), can be implemented
straightforward in momentum space.

Thanks again to the hierarchy of timescales
$\kappa\gg\omega_{\mathrm{scan}}\gg\omega_{R},Ng_{\text{aa}}/(\hbar L^{2})$,
the light field leaking from the cavity allows to efficiently probe the atomic
density profile in real-time. Indeed, since the laser field $\Omega$ changes
much faster than the atomic wavefunction $\psi$, the former performs a full
scan over the cloud before the latter has evolved. This tomographic signal is
directly accessible as a time-modulation of the complex-valued cavity field.

### Example in 1 dimension

Figure 2: (a) One dimensional realization of the protocol. An elongated cloud
is coupled to a standing wave cavity and illuminated by a narrow laser beam
which periodically scans along the $z$-direction across the cloud. (b-c)
Comparison between the atomic density at time $t=500\omega_{R}^{-1}$ obtained
by the real-time dynamics for scanning frequency of the laser
$\omega_{\mathrm{scan}}=62.8\omega_{R}$ (green) and
$\omega_{\mathrm{scan}}=628\omega_{R}$ (red) and the ground state obtained by
the imaginary-time evolution with the effective interaction (dashed black) for
parameters $L=10\lambda_{c}$, $w=1\lambda_{c}$, $N=10^{5}$,
$g_{\text{aa}}=10^{-4}\hbar\omega_{R}\lambda_{c}$, $\tau=10\omega_{R}^{-1}$,
and $V_{0}=-1.55\times 10^{-4}\omega_{R}$. This value of the effective cavity-
mediated interaction strength can be achieved with different combinations of
the microscopic parameters, as long as the separation of time scales is
fulfilled. In particular, we chose $g_{0}=600\omega_{R}$,
$\kappa=250\omega_{R}$, $\Delta_{c}=-2500\omega_{R}$, $\Delta_{a}=-2.5\times
10^{7}\omega_{R}$ and $\Omega_{0}=5.2\times 10^{4}\omega_{R}$. Notice that
rescaling $g_{0}$, $\Omega_{0}$ and $\sqrt{\Delta_{a}}$ by the same factor
leaves the equations of motion invariant. (d-e) Kinetic energy and density-
wave order parameter as a function of time (color solid lines) and the average
of these quantities in the time interval
$[450\omega_{R}^{-1},500\omega_{R}^{-1}]$ (dashed lines). The solid black line
indicates the value obtained for the ground state by imaginary-time evolution.
The time $t_{c}$ indicates the moment at which the critical intensity of the
laser is reached for the transition from the homogeneous to the modulated
phase.

In order to demonstrate the protocol and to test its regime of validity, we
focus on one of the simplest possible situations, which further is
experimentally realizable in a straight-forward way. We consider a quasi-1D
cigar-shaped cloud aligned with the axis of the cavity. The derivation
provided in the previous section remains valid, but one needs to replace
$g_{0}\rightarrow g(\vec{r})=g_{0}\cos(2\pi z/\lambda_{c})$ to take into
account the spatial modulation of the cavity mode. We employ a Gaussian laser
profile $\tilde{\Omega}(z)=\Omega_{0}e^{-z^{2}/w^{2}}$, for which the
effective interaction (4) takes the form

$V(z,z^{\prime})=V_{0}\cos(2\pi z/\lambda_{c})\cos(2\pi
z^{\prime}/\lambda_{c})e^{-|z-z^{\prime}|^{2}/2w^{2}}$ (7)

with

$V_{0}=\frac{\sqrt{2\pi}\Delta_{c}g_{0}^{2}\Omega_{0}^{2}}{\Delta_{a}^{2}(\kappa^{2}+\Delta_{c}^{2})}\frac{w}{L}\,.$

For negative values of $\Delta_{c}$ this interaction is attractive and
competes with the repulsive contact interaction and the kinetic energy. For
small values of $V_{0}$, the ground state of the system thus remains
homogeneous. However, by increasing $V_{0}$ it can undergo a phase transition
to a phase with periodically modulated density. We compute the ground state
resulting from the effective interaction (7) by evolving Eq. (3) in imaginary
time. In order to verify that our laser-painting protocol does reproduce the
physics of the effective atom-only model (3), we compare the latter with the
full dynamics obtained from the evolution of Eqs. (1) and (2) in real time. As
initial conditions, we choose a homogeneous atomic density and zero cavity
field. The laser intensity is ramped up according to
$\tilde{\Omega}(z,t)=\Omega_{0}\tanh(t/\tau)e^{-z^{2}/w^{2}}$, and periodic
boundary conditions are employed.

We analyze the transition between the homogeneous and the modulated phase with
help of the kinetic energy per atom

$E_{\text{kin}}=\frac{1}{N}\frac{\hbar^{2}}{2M}\int_{-L/2}^{L/2}\,dz|\partial_{z}\psi(z)|^{2}$
(8)

and the density-wave order parameter capturing the overlap of the atomic
density with a perfect $\lambda_{c}$-periodic modulation favoured by the long-
range interaction,

$O_{\text{d-w}}=\frac{1}{N}\int_{-L/2}^{L/2}\,dz\cos(2\pi
z/\lambda_{c})|\psi(z)|^{2}.$ (9)

Figure 2 characterizes the protocol, comparing two different regimes for the
scanning rate $\omega_{\mathrm{scan}}$ of the laser spot. For intermediate
rates (green lines), $\omega_{\mathrm{scan}}=62.8\omega_{R}$, the protocol
works as intended. Initially, the system is in the homogeneous state, but some
time after the laser strength has been ramped up across the critical value,
the atomic system transitions to the density-wave phase, showing a modulation
in agreement with one found in the ground state of the effective model (3)
(black-dashed lines in panel (b)). Residual deviations are due to the finite-
time ramp and are present also within the effective model by comparing the
long-time state with the ground state. This is demonstrated in panels (d) and
(e), showing the time evolution of kinetic energy and density-wave order
parameter as green solid lines. Their long-time average has a small deviation
from the ground-state value (solid black line), but coincides with the long-
time average computed using the effective model (3) (green dashed line).

For large scanning rates instead, $\omega_{\mathrm{scan}}=628\omega_{R}$, the
cavity is unable to follow adiabatically the dynamics of the laser (red
lines). The combined dynamics of cavity and laser now has beatings at a
frequency much slower than the respective characteristic timescales, which can
in turn efficiently excite the atoms. These beatings show up in the
oscillations of the laser-cavity-interference potential felt by the atoms
(last term in Equation (1)). This leads to the heating manifested as monotonic
growth of the kinetic energy (red line in Fig. 2(d)). Correspondingly, the
atomic density gets modulated at multiple wavelengths and thus the density-
wave order parameter (9) decreases and oscillates around a value that vanishes
in the thermodynamic limit, as more and more wavelengths (momenta) get
excited.

Figure 3: Kinetic energy per particle (blue) and order parameter (black)
averaged in the time interval $[450\omega_{R}^{-1},500\omega_{R}^{-1}]$ for
different scan rates of the laser, considering the same parameters as in Fig.
2. The dashed lines indicate the value obtained for the ground state of the
static effective model (7). The scan-rate parameter range in which our laser-
painted interaction protocol works is indicated by the green background. Scan
rates which lead to heating (increased kinetic energy) are indicated by the
orange background, while the parameter range in which the interaction becomes
of global range is indicated by the blue background.

In order to provide the full characterization of the protocol with respect to
the robustness of the hierarchy of timescales, we show in Figure 3 the kinetic
energy and order parameter averaged over the time interval
$[450\omega_{R}^{-1},500\omega_{R}^{-1}]$ as a function of the laser scanning
rate $\omega_{\mathrm{scan}}$. For slow scan rates
($\omega_{\mathrm{scan}}\ll\omega_{R}$), we observe the formation of a single
peak in the density produced by the negative potential contribution
$\Omega(\vec{r},t)^{2}/\Delta_{a}$ in equation (1), which decays once the
laser moves on. Since the system is perturbed in all its length during the
ramping of the laser intensity, a strong increment of the kinetic energy that
persists in time is induced. For very fast scan rates
($\omega_{\mathrm{scan}}\gg\Delta_{c}$), however, the laser acts like it is
spatially constant, resulting in a modulation of the atomic density that
cannot be distinguished from global cavity-mediated interactions. For scan
rates comparable to one of the natural time scales of the system
($\omega_{R},\kappa$,$\Delta_{c}$), energy is efficiently transferred to the
atoms, either due to Floquet heating induced by the local laser beam moving
periodically at a rate $\omega_{\rm scan}\sim\omega_{R}$, or due to the
excitation of the cavity by the laser whenever $\omega_{\rm scan}$ is
comparable to a cavity timescale. As discussed above (see Fig.2), in the
latter case heating is caused by beating oscillations in the cavity-laser
interference potential. A plateau of very small kinetic energy and large order
parameter, approaching the respective ground-state values, appears in a broad
range of intermediate scan rates, demonstrating the robustness of our
protocol. For the parameters employed in the figure, this plateau lies roughly
between $15\omega_{R}$ and $500\omega_{R}$. We note also that, when the cavity
loss rate is by far the largest scale, cavity cooling/heating is inefficient
[58].

### Cavity-based tomography of domains

The light leaking from the cavity mirrors can be used to obtain tomographic
real-time information about the density of the atoms, which we demonstrate
using the example from the previous section. Applying the cavity field
elimination (see Methods) to the 1D situation, we obtain

$\alpha(t)=\frac{\int dzg_{0}\cos(2\pi
z/\lambda_{c})\Omega(z,t)|{\psi}(z)|^{2}}{\Delta_{a}(\Delta_{c}+i\kappa)}\,,$
(10)

where the atomic density is considered static as it evolves with the slowest
time scale. In the density-wave phase, the atoms spontaneously break a
$\mathbb{Z}_{2}$ symmetry [40] with density maxima or minima at integer values
of $z/\lambda_{c}$. If the range of the interaction is sufficiently smaller
than the size of the system ($w\ll L$), domains with different spontaneous
breaking of the symmetry can be formed as the critical point is crossed within
a finite time. According to the Kibble-Zurek mechanism [59, 60], we thus
expect in our case the formation of domains with different position of the
density-modulation maxima when we ramp the laser strength across the
homogeneous-to-density-wave critical point. This is in stark contrast to a
situation with global-range interactions which does not allow for any domain
formation.

\begin{overpic}[width=390.25534pt]{Fig_Domains_Map.png} \end{overpic}

\begin{overpic}[width=390.25534pt]{Fig_Domains_Ground_State.pdf} \end{overpic}

\begin{overpic}[width=390.25534pt]{Fig_Domains_Tomography.pdf}
\put(18.5,42.0){(c)} \end{overpic}

Figure 4: (a) Atomic density as a function of time starting from the
homogeneous density for parameters $L=30\lambda_{c}$, $w=1\lambda_{c}$,
$N=10^{5}$, $g_{\text{aa}}=3\times 10^{-4}\hbar\omega_{R}\lambda_{c}$,
$V_{0}=-4.9\times 10^{-4}\omega_{R}$, $\tau=10\omega_{R}^{-1}$, and
$\omega_{\mathrm{scan}}=20.9\omega_{R}$. The blue scale is used at integer
positions and the red one at semi-integer positions. (b) Atomic density
profile at time $t=100\omega_{R}^{-1}$. (c) Amplitude and phase of the cavity
mode $\alpha$ as a function of time during one period of the laser sweep.

The dynamics of domains can be characterized by analysing the cavity field
$\alpha(t)$ during one period of the laser sweep. Figure 4 shows results
computed with the same parameters as in the previous section, but for a larger
system, $L=30\lambda_{c}$, scanning frequency
$\omega_{\mathrm{scan}}=20.9\omega_{R}$ and slightly different interaction
strengths: $g_{\text{aa}}=3\times 10^{-4}\hbar\omega_{R}\lambda_{c}$ and
$V_{0}=-4.9\times 10^{-4}\omega_{R}$. As shown in panels (a) and (b), the
atomic wave function develops two long-lived domains.

How the information about domains is stored in the cavity field dynamics can
be understood analytically. By approximating the atomic-density modulation as
$|\psi(z)|^{2}=(2N/L)\cos(\pi z/\lambda_{c})^{2}$ or
$|\psi(z)|^{2}=(2N/L)\sin(\pi z/\lambda_{c})^{2}$ for the phase with maxima or
minima at integer values of $z/\lambda_{c}$, respectively, we obtain

$\alpha=\pm\frac{g_{0}\Omega_{0}\sqrt{\pi}wN}{\Delta_{a}(\Delta_{c}+i\kappa)2L}(1+\mathcal{O}(e^{-\pi^{2}w^{2}/\lambda_{c}^{2}}))$
(11)

where $\pm$ corresponds to the 2 different cases. The time phase of the cavity
field is thus shifted by $\pi$ as the laser moves from one domain to the
other. This is confirmed in Fig. 4 (c) where the amplitude takes the same
value in both domains but the phase is shifted by $\pi$. Naturally, the
amplitude vanishes when the laser is illuminating the region between domains.
For this simple case of a Gaussian laser beam, the spatial resolution of the
cavity tomography is set by the laser spot itself, i.e. by $w$. While cavity
microscopes rely on overlapping multiple cavity modes for increasing the
resolution [61], our cavity tomograph relies on focusing the laser beam
instead. In both devices, the resolution is at the same time the range of the
interaction, and thus sets also the minimum size of the structures to be
probed.

## Discussion

We theoretically demonstrated that the combination of a single-mode Fabry-
Pérot cavity and a focused, scanning laser allow to implement effective
interactions between atoms which are in principle arbitrarily tuneable in
shape, strength, and sign. Our protocol is well within reach of current
experimental capabilities and also features a built-in tomographic probing
tool based on the cavity output signal. It is important to stress at this
point that the whole protocol works equally also for fermionic atoms [62, 63,
64].

As an experimentally realizable proof of principle, we characterized the
protocol for the case of a Gaussian laser beam sweeping across a quasi-1D BEC
trapped along the axis of a standing-wave cavity. This realizes a periodically
sign-changing interaction, with a Gaussian envelope setting an interaction
range given by the laser spot. We observed a transition between a homogeneous
and a density-modulated BEC, the crossing of which generates different domains
that can be probed in a space- and time-resolved fashion by the built-in
tomographic tool.

Since the interaction potential can be freely painted pixel-by-pixel using
digital beam shaping devices for the laser beam, the range of possibilities
for future investigations is hard to delimit. Elusive ground states of matter
that impose very specific requirements on the interaction potential can be
stabilized and studied, including for instance strongly-correlated solitons
and droplets induced by light [65, 66], quantum cluster quasi-crystals [15,
67], requiring very specific combinations of lengthscales to appear in the
potential, or quantum paramagnets [68, 69], relying on strongly frustrated
interactions in concomitance with specific symmetries. Our scheme thus allows
the experimental implementation of so-far-purely-theoretical interaction
potentials of paradigmatic nature, but might also lead to the discovery of
yet-unknown phenomenology.

## Methods

### Elimination of the cavity field

If the characteristic time scale $1/\kappa$ of the cavity field is fast with
respect to the atomic motion, the cavity field can be adiabatically eliminated
which results in

$\alpha(t)\approx\frac{\int
g_{0}\Omega(\vec{r},t)|{\psi}(\vec{r},t)|^{2}d\vec{r}}{\Delta_{a}(\Delta_{c}+i\kappa)}\,.$
(12)

Inserting this expression into Equation (1) we find an effective description
of the atomic subsystem

$\begin{split}&i\frac{\partial{\psi}(\vec{r},t)}{\partial
t}=\Big{[}-\frac{\hbar\nabla^{2}}{2M}+\frac{g_{\text{aa}}}{\hbar}|{\psi}(\vec{r},t)|^{2}+\frac{\Omega^{2}(\vec{r},t)}{\Delta_{a}}+\\\
&+\frac{2\Delta_{c}g_{0}\Omega(\vec{r},t)}{\Delta_{a}^{2}(\kappa^{2}+\Delta_{c}^{2})}\int
g_{0}\Omega(\vec{r}\,^{\prime},t)|{\psi}(\vec{r}\,^{\prime},t)|^{2}d\vec{r}\,^{\prime}\Big{]}{\psi}(\vec{r},t)\,,\end{split}$
(13)

where we used $|\Delta_{a}\Delta_{c}|\gg g_{0}^{2}N$, and neglected
$\frac{g_{0}^{2}}{\Delta_{a}}|\alpha|^{2}$ given the low occupation of the
cavity mode. From here it is apparent that two atoms at positions $\vec{r}$
and $\vec{r}\,^{\prime}$ only interact if they are simultaneously illuminated
by the laser.

### Averaging over the scan

The laser has velocity $v$ and moves periodically through a path of horizontal
lines separated by a distance $l$. In this way the period is $T=L^{2}/lv$ and
the average (5) can be parameterised as

$\frac{1}{T}\sum_{n=1}^{L/l}\int_{t_{n-1}}^{t_{n}}dt\tilde{\Omega}(\vec{r}-\vec{r}_{0}^{\,n}(t))\tilde{\Omega}(\vec{r}\,^{\prime}-\vec{r}_{0}^{\,n}(t))\,,$
(14)

where the path for each line is given by
${\vec{r}_{0}^{\,n}(t)=(x_{n}(t),y_{n})=(-\frac{L}{2}+v(t-t_{n-1}),-\frac{L}{2}+nl)}$
and the times are $t_{n}=nL/v$. With the previous expression for the path, the
integral in time can be written as an integral in space with the variable $x$.
Furthermore, if the separation between lines is small compared to the
characteristic size of the laser spot $w$, the sum over $n$ can be turned into
an integral in $y$ resulting in

$\frac{1}{L^{2}}\int_{-L/2}^{L/2}dy\int_{-L/2}^{L/2}dx\,\tilde{\Omega}(\vec{r}-(x,y))\tilde{\Omega}(\vec{r}\,^{\prime}-(x,y))\,.$
(15)

Given that the laser has finite range, the integrand vanishes when the point
$(x,y)$ is far away from $\vec{r}$ and $\vec{r}\,^{\prime}$. Consequently, the
integrals in space can be extended to $\pm\infty$ without changing the result.
We can additionally make the change of coordinates
$(x,y)\rightarrow\vec{r}\,^{\prime}+(x,y)$ and obtain the expression

$\frac{1}{L^{2}}\int_{-\infty}^{\infty}dy\int_{-\infty}^{\infty}dx\,\tilde{\Omega}(\vec{r}-\vec{r}\,^{\prime}-(x,y))\tilde{\Omega}(-(x,y)).$
(16)

## Acknowledgement

We are thankful for insightful discussions with Alexander Baumgärtner, David
Baur, Simon Hertlein, Gabriele Natale, and Justyna Stefaniak. T.D.
acknowledges funding from the Swiss National Science Foundation SNF (Project
No. IZBRZ2 186312), and funding from the Swiss State Secretariat for
Education, Research and Innovation (SERI).

## References

  * Feynman [1982] R. P. Feynman, Simulating physics with computers, International Journal of Theoretical Physics 21, 467 (1982).
  * Bloch  _et al._ [2012] I. Bloch, J. Dalibard, and S. Nascimbène, Quantum simulations with ultracold quantum gases, Nature Physics 8, 267 (2012).
  * Gross and Bloch [2017] C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, Science 357, 995 (2017).
  * Schäfer  _et al._ [2020] F. Schäfer, T. Fukuhara, S. Sugawa, Y. Takasu, and Y. Takahashi, Tools for quantum simulation with ultracold atoms in optical lattices, Nature Reviews Physics 2, 411 (2020).
  * Zohar  _et al._ [2016] E. Zohar, J. I. Cirac, and B. Reznik, Quantum simulations of lattice gauge theories using ultracold atoms in optical lattices, Reports on Progress in Physics 79, 014401 (2016).
  * Viermann  _et al._ [2022] C. Viermann, M. Sparn, N. Liebster, M. Hans, E. Kath, Á. Parra-López, M. Tolosa-Simeón, N. Sánchez-Kuntz, T. Haas, H. Strobel, S. Floerchinger, and M. K. Oberthaler, Quantum field simulator for dynamics in curved spacetime, Nature 611, 260 (2022).
  * Defenu  _et al._ [2023] N. Defenu, T. Donner, T. Macrì, G. Pagano, S. Ruffo, and A. Trombettoni, Long-range interacting quantum systems, Rev. Mod. Phys. 95, 035002 (2023).
  * Dalmonte  _et al._ [2010] M. Dalmonte, G. Pupillo, and P. Zoller, One-Dimensional Quantum Liquids with Power-Law Interactions: The Luttinger Staircase, Phys. Rev. Lett. 105, 140401 (2010).
  * Pupillo  _et al._ [2010] G. Pupillo, A. Micheli, M. Boninsegni, I. Lesanovsky, and P. Zoller, Strongly Correlated Gases of Rydberg-Dressed Atoms: Quantum and Classical Dynamics, Phys. Rev. Lett. 104, 223002 (2010).
  * Cinti  _et al._ [2014] F. Cinti, T. Macrì, W. Lechner, G. Pupillo, and T. Pohl, Defect-induced supersolidity with soft-core bosons, Nature communications 5, 3235 (2014).
  * Mendoza-Coto  _et al._ [2015] A. Mendoza-Coto, D. A. Stariolo, and L. Nicolao, Nature of Long-Range Order in Stripe-Forming Systems with Long-Range Repulsive Interactions, Phys. Rev. Lett. 114, 116101 (2015).
  * Rossotti  _et al._ [2017] S. Rossotti, M. Teruzzi, D. Pini, D. E. Galli, and G. Bertaina, Quantum Critical Behavior of One-Dimensional Soft Bosons in the Continuum, Phys. Rev. Lett. 119, 215301 (2017).
  * Angelone  _et al._ [2016] A. Angelone, F. Mezzacapo, and G. Pupillo, Superglass Phase of Interaction-Blockaded Gases on a Triangular Lattice, Phys. Rev. Lett. 116, 135303 (2016).
  * Cinti and Macrì [2019] F. Cinti and T. Macrì, Thermal and Quantum Fluctuation Effects in Quasiperiodic Systems in External Potentials, Condensed Matter 4, 93 (2019).
  * Pupillo  _et al._ [2020] G. Pupillo, P. c. v. Ziherl, and F. Cinti, Quantum cluster quasicrystals, Phys. Rev. B 101, 134522 (2020).
  * Micheli  _et al._ [2006] A. Micheli, G. K. Brennen, and P. Zoller, A toolbox for lattice-spin models with polar molecules, Nature Physics 2, 341 (2006).
  * van Bijnen and Pohl [2015] R. M. W. van Bijnen and T. Pohl, Quantum Magnetism and Topological Ordering via Rydberg Dressing near Förster Resonances, Physical Review Letters 114, 243002 (2015).
  * Samajdar  _et al._ [2021] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Quantum phases of Rydberg atoms on a kagome lattice, Proceedings of the National Academy of Sciences of the United States of America 118, 1 (2021).
  * Yao  _et al._ [2013] N. Y. Yao, A. V. Gorshkov, C. R. Laumann, A. M. Läuchli, J. Ye, and M. D. Lukin, Realizing Fractional Chern Insulators in Dipolar Spin Systems, Physical Review Letters 110, 185302 (2013).
  * Scarola and Das Sarma [2005] V. W. Scarola and S. Das Sarma, Quantum Phases of the Extended Bose-Hubbard Hamiltonian: Possibility of a Supersolid State of Cold Atoms in Optical Lattices, Physical Review Letters 95, 033003 (2005).
  * Léonard  _et al._ [2017] J. Léonard, A. Morales, P. Zupancic, T. Esslinger, and T. Donner, Supersolid formation in a quantum gas breaking a continuous translational symmetry, Nature 543, 87 (2017).
  * Tanzi  _et al._ [2019] L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, Observation of a Dipolar Quantum Gas with Metastable Supersolid Properties, Physical Review Letters , 130405 (2019).
  * Böttcher  _et al._ [2019] F. Böttcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, Transient Supersolid Properties in an Array of Dipolar Quantum Droplets, Physical Review X 9, 011051 (2019).
  * Chomaz  _et al._ [2019] L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, and F. Ferlaino, Long-Lived and Transient Supersolid Behaviors in Dipolar Quantum Gases, Physical Review X 9, 021012 (2019).
  * Santos  _et al._ [2003] L. Santos, G. V. Shlyapnikov, and M. Lewenstein, Roton-Maxon Spectrum and Stability of Trapped Dipolar Bose-Einstein Condensates, Physical Review Letters 90, 250403 (2003).
  * Mottl  _et al._ [2012] R. Mottl, F. Brennecke, K. Baumann, R. Landig, T. Donner, and T. Esslinger, Roton-Type Mode Softening in a Quantum Gas with Cavity-Mediated Long-Range Interactions, Science 336, 1570 (2012).
  * Chomaz  _et al._ [2016] L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Quantum-Fluctuation-Driven Crossover from a Dilute Bose-Einstein Condensate to a Macrodroplet in a Dipolar Quantum Fluid, Physical Review X 6, 041039 (2016).
  * Ferrier-Barbut  _et al._ [2016] I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel, and T. Pfau, Observation of Quantum Droplets in a Strongly Dipolar Bose Gas, Physical Review Letters 116, 215301 (2016).
  * Cabrera  _et al._ [2018] C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, Quantum liquid droplets in a mixture of Bose-Einstein condensates, Science 359, 301 (2018).
  * Kasuya [1956] T. Kasuya, A Theory of Metallic Ferro- and Antiferromagnetism on Zener’s Model, Progress of Theoretical Physics 16, 45 (1956).
  * Yosida [1957] K. Yosida, Magnetic Properties of Cu-Mn Alloys, Physical Review 106, 893 (1957).
  * Ruderman and Kittel [1954] M. A. Ruderman and C. Kittel, Indirect Exchange Coupling of Nuclear Magnetic Moments by Conduction Electrons, Physical Review 96, 99 (1954).
  * Argüello-Luengo  _et al._ [2022] J. Argüello-Luengo, A. González-Tudela, and D. González-Cuadra, Tuning Long-Range Fermion-Mediated Interactions in Cold-Atom Quantum Simulators, Physical Review Letters 129, 083401 (2022).
  * Manovitz  _et al._ [2020] T. Manovitz, Y. Shapira, N. Akerman, A. Stern, and R. Ozeri, Quantum Simulations with Complex Geometries and Synthetic Gauge Fields in a Trapped Ion Chain, PRX Quantum 1, 020303 (2020).
  * Periwal  _et al._ [2021] A. Periwal, E. S. Cooper, P. Kunkel, J. F. Wienand, E. J. Davis, and M. Schleier-Smith, Programmable interactions and emergent geometry in an array of atom clouds., Nature 600, 630 (2021).
  * Chomaz  _et al._ [2023] L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. Laburthe-Tolra, B. L. Lev, and T. Pfau, Dipolar physics: a review of experiments with magnetic quantum gases, Reports on Progress in Physics 86, 026401 (2023).
  * Böttcher  _et al._ [2021] F. Böttcher, J. N. Schmidt, J. Hertkorn, K. S. Ng, S. D. Graham, M. Guo, T. Langen, and T. Pfau, New states of matter with fine-tuned interactions: Quantum droplets and dipolar supersolids, Reports on Progress in Physics 84 (2021).
  * Baranov  _et al._ [2012] M. A. Baranov, M. Dalmonte, G. Pupillo, and P. Zoller, Condensed matter theory of dipolar quantum gases, Chemical Reviews 112, 5012 (2012).
  * Adams  _et al._ [2020] C. S. Adams, J. D. Pritchard, and J. P. Shaffer, Rydberg atom quantum technologies, Journal of Physics B: Atomic, Molecular and Optical Physics 53, 012002 (2020).
  * Mivehvar  _et al._ [2021] F. Mivehvar, F. Piazza, T. Donner, and H. Ritsch, Cavity QED with quantum gases: new paradigms in many-body physics, Advances in Physics 70, 1 (2021).
  * Baumann  _et al._ [2010] K. Baumann, C. Guerlin, F. Brennecke, and T. Esslinger, Dicke quantum phase transition with a superfluid gas in an optical cavity, Nature 464, 1301 (2010).
  * Klinder  _et al._ [2015a] J. Klinder, H. Keßler, M. Wolke, L. Mathey, and A. Hemmerich, Dynamical phase transition in the open Dicke model, Proceedings of the National Academy of Sciences of the United States of America 112, 3290 (2015a).
  * Ferri  _et al._ [2021] F. Ferri, R. Rosa-Medina, F. Finger, N. Dogra, M. Soriente, O. Zilberberg, T. Donner, and T. Esslinger, Emerging dissipative phases in a superradiant quantum gas with tunable decay, Physical Review X 11, 41046 (2021).
  * Klinder  _et al._ [2015b] J. Klinder, H. Keßler, M. R. Bakhtiari, M. Thorwart, and A. Hemmerich, Observation of a Superradiant Mott Insulator in the Dicke-Hubbard Model, Physical Review Letters 115, 230403 (2015b).
  * Landig  _et al._ [2016] R. Landig, L. Hruby, N. Dogra, M. Landini, R. Mottl, T. Donner, and T. Esslinger, Quantum phases from competing short- and long-range interactions in an optical lattice, Nature 532, 476 (2016).
  * Vaidya  _et al._ [2018] V. D. Vaidya, Y. Guo, R. M. Kroeze, K. E. Ballantine, A. J. Kollár, J. Keeling, and B. L. Lev, Tunable-Range, Photon-Mediated Atomic Interactions in Multimode Cavity QED, Physical Review X 8, 011002 (2018).
  * Guo  _et al._ [2019] Y. Guo, R. M. Kroeze, V. D. Vaidya, J. Keeling, and B. L. Lev, Sign-Changing Photon-Mediated Atom Interactions in Multimode Cavity Quantum Electrodynamics, Physical Review Letters 122, 193601 (2019).
  * Guo  _et al._ [2021] Y. Guo, R. M. Kroeze, B. P. Marsh, S. Gopalakrishnan, J. Keeling, and B. L. Lev, An optical lattice with sound, Nature 599, 211 (2021).
  * O’Dell  _et al._ [2003] D. H. J. O’Dell, S. Giovanazzi, and G. Kurizki, Rotons in Gaseous Bose-Einstein Condensates Irradiated by a Laser, Physical Review Letters 90, 110402 (2003).
  * Giovanazzi  _et al._ [2002] S. Giovanazzi, D. O’Dell, and G. Kurizki, Density Modulations of Bose-Einstein Condensates via Laser-Induced Interactions, Physical Review Letters 88, 4 (2002).
  * Grimm  _et al._ [2000] R. Grimm, M. Weidemüller, and Y. B. Ovchinnikov, Optical Dipole Traps for Neutral Atoms, in  _Advances in Atomic, Molecular, and Optical Physics_ , Vol. 42 (2000) pp. 95–170.
  * Ostermann  _et al._ [2016] S. Ostermann, F. Piazza, and H. Ritsch, Spontaneous Crystallization of Light and Ultracold Atoms, Physical Review X 6, 021026 (2016).
  * Dimitrova  _et al._ [2017] I. Dimitrova, W. Lunden, J. Amato-Grill, N. Jepsen, Y. Yu, M. Messer, T. Rigaldo, G. Puentes, D. Weld, and W. Ketterle, Observation of two-beam collective scattering phenomena in a Bose-Einstein condensate, Physical Review A 96, 051603 (2017).
  * Münstermann  _et al._ [2000] P. Münstermann, T. Fischer, P. Maunz, P. W. H. Pinkse, and G. Rempe, Observation of Cavity-Mediated Long-Range Light Forces between Strongly Coupled Atoms, Phys. Rev. Lett. 84, 4068 (2000).
  * Schnelle  _et al._ [2008] S. K. Schnelle, E. D. van Ooijen, M. J. Davis, N. R. Heckenberg, and H. Rubinsztein-Dunlop, Versatile two-dimensional potentials for ultra-cold atoms, Opt. Express 16, 1405 (2008).
  * Piazza  _et al._ [2013] F. Piazza, P. Strack, and W. Zwerger, Bose–Einstein condensation versus Dicke–Hepp–Lieb transition in an optical cavity, Annals of Physics 339, 135 (2013).
  * Lang  _et al._ [2017] J. Lang, F. Piazza, and W. Zwerger, Collective excitations and supersolid behavior of bosonic atoms inside two crossed optical cavities, New Journal of Physics 19, 123027 (2017).
  * Piazza and Strack [2014] F. Piazza and P. Strack, Quantum kinetics of ultracold fermions coupled to an optical resonator, Physical Review A 90, 043823 (2014).
  * Kibble [1976] T. W. Kibble, Topology of cosmic domains and strings, Journal of Physics A: Mathematical and General 9, 1387 (1976).
  * Zurek [1985] W. H. Zurek, Cosmological experiments in superfluid helium?, Nature 317, 505 (1985).
  * Kroeze  _et al._ [2023] R. M. Kroeze, B. P. Marsh, K.-Y. Lin, J. Keeling, and B. L. Lev, High cooperativity using a confocal-cavity–qed microscope, PRX Quantum 4, 020326 (2023).
  * Zhang  _et al._ [2021] X. Zhang, Y. Chen, Z. Wu, J. Wang, J. Fan, S. Deng, and H. Wu, Observation of a superradiant quantum phase transition in an intracavity degenerate Fermi gas, Science 373, 1359 (2021).
  * Helson  _et al._ [2023] V. Helson, T. Zwettler, F. Mivehvar, E. Colella, K. Roux, H. Konishi, H. Ritsch, and J.-P. Brantut, Density-wave ordering in a unitary Fermi gas with photon-mediated interactions, Nature , 1 (2023).
  * Mivehvar  _et al._ [2017] F. Mivehvar, H. Ritsch, and F. Piazza, Superradiant topological peierls insulator inside an optical cavity, Physical review letters 118, 073602 (2017).
  * Fraser and Piazza [2019] K. A. Fraser and F. Piazza, Topological soliton-polaritons in 1D systems of light and fermionic matter, Communications Physics 2, 48 (2019).
  * Karpov and Piazza [2022] P. Karpov and F. Piazza, Light-Induced Quantum Droplet Phases of Lattice Bosons in Multimode Cavities, Phys. Rev. Lett. 128, 103201 (2022).
  * Mivehvar  _et al._ [2019a] F. Mivehvar, H. Ritsch, and F. Piazza, Emergent quasicrystalline symmetry in light-induced quantum phase transitions, Physical Review Letters 123, 210604 (2019a).
  * Mivehvar  _et al._ [2019b] F. Mivehvar, H. Ritsch, and F. Piazza, Cavity-quantum-electrodynamical toolbox for quantum magnetism, Physical Review Letters 122, 113603 (2019b).
  * Chiocchetta  _et al._ [2021] A. Chiocchetta, D. Kiese, C. P. Zelle, F. Piazza, and S. Diehl, Cavity-induced quantum spin liquids, Nature Communications 12, 5901 (2021).