# Denoising Low-dose Images Using
Deep Learning of Time Series Images

Yang SHAO<EMAIL_ADDRESS>[  Toshie YAGUCHI  [  Toshiaki TANIGAKI  [

###### Abstract

Digital image devices have been widely applied in many fields, including
scientific imaging, recognition of individuals, and remote sensing. As the
application of these imaging technologies to autonomous driving and
measurement, image noise generated when observation cannot be performed with a
sufficient dose has become a major problem. Machine learning denoise
technology is expected to be the solver of this problem, but there are the
following problems. Here we report, artifacts generated by machine learning
denoise in ultra-low dose observation using an in-situ observation video of an
electron microscope as an example. And as a method to solve this problem, we
propose a method to decompose a time series image into a 2D image of the
spatial axis and time to perform machine learning denoise. Our method opens
new avenues accurate and stable reconstruction of continuous high-resolution
images from low-dose imaging in science, industry, and life.

###### keywords:

American Chemical Society, LaTeX

hitachi] Research and Development Group, Hitachi, Ltd. hightech] Hitachi High-
Tech Corporation hitachi] Research and Development Group, Hitachi, Ltd.
IR,NMR,UV

##  1 Introduction

In the realm of scientific research, Electron Microscopy has emerged as an
indispensable tool, offering unprecedented insights into the microscopic
world. Despite its profound contributions, capturing high-quality images under
low-dose conditions remains a formidable challenge. This challenge stems from
the inherent trade-off between image quality and exposure time: while longer
exposure times can yield high-resolution images, they also increase the risk
of sample damage 1. Conversely, shorter exposure times, although safer for the
sample, often result in images plagued by noise and low resolution 2. There
are several attempts to denoise low-dose images by using machine learning 3,
4, 5. However, the denoising for live movies and time series images are still
challenging.

In light of these challenges, this paper proposes a novel approach that
leverages the power of machine learning and three-dimensional (3D) synthesis
to enhance low-dose microscopy images. At the heart of our method is a two-
step process: the generation of synthetic data "(Section 2.2)" followed by
microscope calibration "(Section 2.)". The synthetic data, generated using a
simulator, incorporates various modes of random noise, providing a diverse and
robust training set for our base model. The microscope calibration, performed
prior to actual microscopy observations, involves fine-tuning the base model
with low-noise high-resolution images acquired using intense irradiation.

In addition to this machine learning-based enhancement, we introduce a
computationally efficient 3D image processing technique based on three-
directional synthesis. This technique begins by capturing a three-dimensional
tensor representing the observed scene. The tensor is then sliced along three
orthogonal directions, yielding three continuous image segments. Each segment
undergoes individual noise reduction enhancement processes, resulting in three
separate 3D tensors. The final image is obtained by computing the average of
these tensors, thereby creating a low-noise high-resolution representation of
the original scene.

By combining details from multiple perspectives, our method ensures accurate
and stable reconstruction of continuous high-resolution images from low-dose
observations. Note that this method can be applied to time-resolved electron
microscopy 6, which also requires low-dose imaging. This research opens new
avenues for high-resolution microscopy in challenging lighting conditions. The
potential applications of this work span across various fields, including
materials science, biology, and medical diagnostics.

##  2 Research Method

###  2.1 Pre-training fine-tuning paradigm

In the early days of machine learning, feature engineering and handcrafted
kernel functions were the norm. Techniques like Support Vector Machines (SVMs)
thrived on these manual designs. However, the landscape shifted dramatically
with the rise of deep neural networks. Features no longer demanded human
craftsmanship; instead, the paradigm shifted toward designing neural network
structures tailored to specific problems, employing diverse training
algorithms. Today, machine learning is facing another paradigm shift, and a
new era has emerged - the era of universal basic models. This paradigm
dispenses with the need for bespoke neural network architectures for every
distinct problem. Instead, it harnesses the power of pre-trained base models,
followed by fine-tuning.

The main two phases of pre-training fine-tuning paradigm are as follows. Pre-
training: Constructing a versatile scaffold—a base model—by training it on a
wealth of diverse tasks. This base model learns fundamental features, akin to
a polymath absorbing knowledge across disciplines. It becomes a repository of
general patterns, ready to adapt. Fine-Tuning: Faced with a specific problem,
instead of reinventing the wheel, we fine-tune the base model to adapt it to
the nuances of the task, becoming a specialized tool.

The major advantages and benefits are as follows. Transfer Learning:
Leveraging pre-trained base models accelerates learning. We inherit wisdom
from a multitude of tasks, akin to standing on the shoulders of giants.
Customization: Fine-tuning tailors the base model to specific challenges. It’s
like adjusting a telescope’s focus for celestial clarity. No need to start
from scratch, we build upon existing knowledge. Efficiency: Gone are the days
of designing bespoke architectures and conducting separate training sessions.
Base models provide a head start, and fine-tuning ensures adaptability.

In summary, the pre-training fine-tuning paradigm bridges the gap between the
massive amounts of data required by traditional machine learning methods and
the limited data available for domain-specific needs.

###  2.2 Three-dimensional (3D) synthesis

In recent years, video denoising has gained significant attention due to its
applications in various domains such as surveillance, entertainment, and
medical imaging. Traditional methods often rely on handcrafted features or
filters, which may not fully exploit the rich spatial and temporal information
present in videos. Therefore, deep learning techniques have emerged as
powerful tools for denoising tasks. However, when applying deep neural
networks to video denoising, a common approach involves using three-
dimensional tensors (3D volumes) as input. These 3D tensors represent video
frames along the spatial and temporal dimensions. While effective, this
approach has limitations. The resulting neural network can become quite large,
requiring substantial computational resources for both training and inference.

To address this challenge, we propose a method that leverages sliced three-
dimensional tensors. Here’s how it works:

  * •

Tensor Slicing: Given the original 3D tensor representing a video sequence, we
slice it along three orthogonal directions (X, Y, and T). This results in
three stacks of 2D images, each capturing information from a specific
perspective. The X-slice contains frames along the horizontal spatial
dimension. The Y-slice contains frames along the vertical spatial dimension.
The T-slice contains frames along the temporal dimension.

  * •

Three Separate Models: Next, we build three separate neural network models,
each learned to denoise the 2D images from one of the slices. These models are
smaller and more manageable than a single large network. The X-model processes
the X-slice images. The Y-model processes the Y-slice images. The T-model
processes the T-slice images.

  * •

Fusion and Synthesis: After denoising each slice independently, we combine the
results. Specifically: The denoised X-slice, Y-slice, and T-slice are
reassembled into a 3D tensor. The final denoised video is synthesized from
this reassembled tensor.

By breaking down the problem into three smaller models, our approach reduces
the overall complexity of the denoising process. Training becomes more
efficient, as each model focuses on a specific aspect of denoising. During
inference, the computational speed improves due to the smaller model sizes.
The schematic diagram of three-dimensional slices is shown in Figure 1.

Figure 1: Schematic diagram of three-dimensional slices

###  2.3 Artificial data for base model

We initially used an electromagnetic field generated near the PN junction
through simulation as our ground truth. After achieving satisfactory results,
we further utilized actual electron microscope images of catalysts obtained
from previous captures. To strike a balance between higher observation clarity
and computational time requirements, we temporarily set an observation window
size of 128x128. As computational power increases, we will be able to train
models capable of processing higher-resolution images within the same time
frame. Additionally, we can achieve processing of higher-resolution images
without altering the model by segmenting the image and processing multiple
observation windows, subsequently reassembling the results.

To generate artificial data for training our foundational model, we introduced
varying levels of different types of noise to the ground truth images. As two
typical examples, we added Poisson noise and Gaussian noise. In order to
enhance the foundational model’s ability to handle different types of noise,
we will continue to incorporate additional noise types and explore various
hybrid noise models. An example of actual captured images and images with
artificial noise is shown in Figure 2.

Figure 2: Actual captured images and images with artificial noise

###  2.4 Structure of base model

In the field of image generation, there are currently two prominent
architectural approaches: autoencoder-based and Diffusion-based architectures.
As an initial exploration, our research employed the well-known U-Net
architecture. This architecture is specifically designed for semantic
segmentation tasks. Let’s delve into its components in detail:

  1. 1.

Contracting Path (Encoder):

     * •

The contracting path follows the typical architecture of a convolutional
network.

     * •

It consists of repeated applications of two 3x3 convolutions (unpadded
convolutions), each followed by a rectified linear unit (ReLU).

     * •

Additionally, a 2x2 max pooling operation with a stride of 2 is used for
downsampling.

     * •

At each downsampling step, the number of feature channels is doubled.

  2. 2.

Expansive Path (Decoder):

     * •

The expansive path aims to recover the spatial resolution lost during
downsampling.

     * •

It consists of:

       * –

Upsampling of the feature map.

       * –

A 2x2 convolution ("up-convolution") that halves the number of feature
channels.

       * –

Concatenation with the correspondingly cropped feature map from the
contracting path.

       * –

Two 3x3 convolutions, each followed by a ReLU.

     * •

The cropping step is necessary due to the loss of border pixels in every
convolution.

  3. 3.

Final Layer:

     * •

At the final layer, a 1x1 convolution is used to map each 64-component feature
vector to the desired number of classes.

In total, the U-Net network comprises 23 convolutional layers. While this
architecture served as our initial choice, we plan to explore other more
efficient architectures in our subsequent research.

###  2.5 Training of base model

In our study, we used a dataset comprising 300 sets of 128 high-resolution
images of size 128x128 pixels. These images were artificially corrupted with
noise to serve as the input to our U-Net model as described above. The output
of the U-Net model, which are the denoised images, were then compared with the
corresponding noise-free high-resolution images. The Mean Squared Error (MSE)
between the output and the noise-free images was computed and used as the loss
function for training the U-Net model. The MSE loss function is mathematically
represented as:

$\text{{MSE}}=\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-\hat{Y}_{i})^{2}$

where $Y_{i}$ represents the value of pixels in actual noise-free images,
$\hat{Y}_{i}$ represents the value of pixels in denoised images produced by
the U-Net model, and $n$ is the total number of pixels. In order to improve
the learning effect, we normalized all pixels so that the values of all pixels
fall between 0-1.

To optimize the U-Net model, we employed the Backpropagation algorithm, which
is a widely-used method in training neural networks. Backpropagation computes
the gradient of the loss function with respect to the model parameters and
uses this gradient to update the parameters. The mathematical representation
of the parameter update rule in Backpropagation is:

$\theta=\theta-\eta\cdot\nabla_{\theta}J(\theta)$

where $\theta$ represents the parameters of the model, $\eta$ is the learning
rate, and $\nabla_{\theta}J(\theta)$ is the gradient of the loss function
$J(\theta)$ with respect to the parameters.

For the optimization process, we used the Adam optimizer. Adam, short for
Adaptive Moment Estimation, is an optimization algorithm that computes
adaptive learning rates for each parameter. In contrast to classical
stochastic gradient descent, Adam keeps an exponentially decaying average of
past gradients and past squared gradients, and uses these to scale the
learning rate. We set initial learning rate of Adam as 0.001.

In an effort to balance learning accuracy and efficiency, we employed
PyTorch’s autocast feature in our study. Autocasting is a feature that
automatically chooses the data type - float32 or float16 - for each operation
in order to improve the performance of the model while maintaining the
accuracy. By using autocast, we were able to convert some of the float32 data
to float16 for computation. This allowed us to leverage the benefits of
float16, such as reduced memory usage and increased computational speed,
without significantly compromising the accuracy of our model. This approach is
particularly beneficial when training large models or when memory capacity is
a limiting factor.

Due to the limitation of computing device memory space, we adopted a
relatively small batch size of 4 and conducted learning for 50 epochs.

We trained three different U-Net models for the aforementioned three
directions, x-y, x-t, y-t.

###  2.6 Comparisons

In our study, we employed a comparative approach to evaluate the effectiveness
of our proposed denoising method. To this end, we utilized three well-
established denoising techniques as benchmarks: fastNlMeansDenoising,
GaussianBlur and bilateralFilter.

  1. 1.

fastNlMeansDenoising: This non-local means denoising algorithm works by
replacing the value of a pixel by an average of a selection of other pixels
values. The selection is based on a patch similarity criterion, leading to
effective noise reduction while preserving the structural details of the
image.

  2. 2.

GaussianBlur: Gaussian blur is a widely used effect in graphics software,
typically to reduce image noise and detail. The visual effect of this blurring
technique is a smooth blur resembling that of viewing the image through a
translucent screen.

  3. 3.

bilateralFilter: The bilateral filter is a non-linear, edge-preserving, and
noise-reducing smoothing filter for images. It replaces the intensity of each
pixel with a weighted average of intensity values from nearby pixels. This
weight can be based on a Gaussian distribution, thus preserving sharp edges by
systematically looping through each pixel and adjusting weights to the
adjacent pixels accordingly.

We used the opencv implementation of these methods. They were applied to the
same noisy datasets for a fair comparison. The results of our method were then
juxtaposed with the outcomes of these techniques, providing a comprehensive
evaluation of our proposed method’s performance in image denoising tasks.

##  3 Experiments

We conducted two experiments. One was designed to test the learning
performance of U-Net, and the other was aimed at evaluating the effectiveness
of 3D synthesis.

###  3.1 Comparative Analysis of Denoising Techniques

For the first experiment, we utilized electromagnetic field images near the PN
junction generated by a simulator, as well as actual electron microscope
images of catalysts that consist of platinum and carbon particles as our
experimental data. The results of U-Net were compared with those of
fastNlMeans, GaussianBlur, and bilateralFilter. We used the Mean Squared Error
(MSE) of all pixels as the comparison metric. The PN junction data we used to
consist of 5120 pairs of images with artificially added noise and their
corresponding images before noise addition. We used four-fifths of this data
for training and the remaining one-fifth for testing. The real image dataset
captured by the electron microscope consisted of 300 groups of 128 consecutive
images, each of size 128x128. We used five-sixths of this data for training
and the remaining one-sixth for testing.

The results are shown in Table 1. The table presents the Normalized Mean
Squared Error (MSE) of pixels for different methods applied to various
datasets. The MSE is measured in units of $10^{-3}$. For the PN junction in
the x-y plane, the original noised data had an MSE of 175.08. After applying
our U-Net, the MSE was significantly reduced to 8.91, which is the lowest
among all methods. Other methods such as fastNlMeansDenoising, GaussianBlur,
and bilateralFilter also reduced the MSE, but not as effectively as our U-Net.
Similarly, for the PN junction in the x-t plane, the original noised data had
an MSE of 156.12. Our U-Net was able to reduce this to 8.01, again achieving
the best performance compared to other methods. For the catalyst particles in
the x-y plane, the original noised data had a much higher MSE of 268.98. Our
U-Net reduced this to 164.25, outperforming GaussianBlur, which resulted in an
MSE of 176.33. The fastNlMeansDenoising and bilateralFilter methods were not
applied in this case. These results demonstrate the superior performance of
our U-Net in reducing the MSE in different scenarios, outperforming other
commonly used methods. This suggests that our U-Net could be a valuable tool
for noise reduction in these applications.

| Original noised | Our U-Net | fN | GB | bF
---|---|---|---|---|---
PN junction x-y | 175.08 | 8.91 | 10.33 | 9.24 | 10.11
PN junction x-t | 156.12 | 8.01 | 9.48 | 8.56 | 9.91
Catalyst x-y | 268.98 | 164.25 | - | 176.33 | -

Table 1: Normalized MSE of pixels (10e-3)

fN: fastNlMeansDenoising GB: GaussianBlur bF: bilateralFilter

###  3.2 Comparative Analysis of 3D synthesis

Figure 3: The particles in the red box maintain the continuity and stability
of the display in 4 consecutive frames. (Line 1: Noisy in-put. Line 2: Actual
observed image before adding noise. Line 3: Denoising result using only the
x-y direction model. Line 4: Denoising results after 3D synthesis) Figure 4:
Phantom particles appeared in the first frame of the denoising result using
only the x-y direction model. (Line 1: Noisy input. Line 2: Actual observed
image before adding noise. Line 3: Denoising result using only the x-y
direction model. Line 4: Denoising results after 3D synthesis.)

For the second experiment, we utilized a real image dataset composed of 300
groups of 128 consecutive images, each of size 128x128, captured by an
electron microscope. Each group of 128 images was sliced along the t, y, and x
axes, resulting in 38,400 images each for the x-y, x-t, and y-t directions. We
used five-sixths of these images for training and the remaining one-sixth for
testing.

We trained three separate U-Net models to handle slices from the x-y, x-t, and
y-t directions. During testing, we reassembled the 128 slices from the same
group into a three-dimensional tensor. The final synthesized tensor was
obtained by averaging the three-dimensional tensors derived from slices in the
three different directions.

To validate the effectiveness of our 3D synthesis method, we compared the
results of the final synthesized tensor in the x-y direction with the results
obtained by directly applying the x-y U-Net model to the x-y images. This
comparison allowed us to assess the performance and potential advantages of
our 3D synthesis approach.

Although the three-dimensional synthesis method has not yet shown significant
effects in reducing the overall signal-to-noise ratio, we have observed two
significant qualitative effects from the experimental results. Firstly, our
observations demonstrated increased continuity and stability. Specifically,
instances where particles temporarily faded or even disappeared due to noise
interference were significantly reduced. This improvement suggests that the
three-dimensional synthesis method enhances the reliability of our
observations by mitigating the impact of noise. Secondly, we observed a
significant decrease in the occurrence of phantom particles, which are
artifacts that appear due to the aggregation of noise at specific locations.
This reduction in phantom particles further attests to the effectiveness of
the three-dimensional synthesis method in improving the quality of our
observations. Figures 3 and Figure 4 provide examples of these improvements in
continuity and stability, and the reduction of phantom particles,
respectively. For figure 3, in the first and fourth frames of the denoising
result using only the x-y direction model, the particles appear faded compared
to the second and third frames. However, in the noising results after 3D
synthesis, all particles maintain a relatively consistent clarity. For figure
4, in the first frame of the denoising result using only the x-y direction
model, phantom particles appeared. However, this phenomenon did not occur in
any of the noising results after 3D synthesis. These findings underscore the
potential of the three-dimensional synthesis method in enhancing the accuracy
and reliability of observations in our field of study.

##  4 Conclusion

This paper has presented a preliminary exploration into the potential of
enhancing low-dose microscopic images by harnessing the power of machine
learning and three-dimensional (3D) synthesis. Our method, while novel, is
just an initial step in addressing the complex challenge of improving image
quality under low-dose conditions.

Through two carefully designed experiments, we have demonstrated the potential
of our approach. The first experiment, conducted on both simulated PN junction
electromagnetic field images and real-world images of the platinum and carbon
particles, provided encouraging results. Our U-Net-based denoising method
showed promising performance compared to traditional techniques. The second
experiment highlighted the benefits of our three-dimensional synthesis method.
By applying this method to the actual electron microscope images of catalysts
consisting of platinum and carbon particles, we observed a noticeable
improvement in image continuity and stability.

Looking forward, we see several areas for future exploration. We plan to
investigate more sophisticated machine learning models and noise reduction
techniques to further enhance the performance of our method. We also intend to
optimize our 3D synthesis method further to improve computational efficiency
and image quality.

## References

  * Egerton 2013 Egerton, R. Control of radiation damage in the TEM. _Ultramicroscopy_ 2013, _127_ , 100–108
  * Lee et al. 2014 Lee, Z.; Rose, H.; Lehtinen, O.; Biskupek, J.; Kaiser, U. Electron dose dependence of signal-to-noise ratio, atom contrast and resolution in transmission electron microscope images. _Ultramicroscopy_ 2014, _145_ , 3–12
  * Lehtinen et al. 2018 Lehtinen, J.; Munkberg, J.; Hasselgren, J.; Laine, S.; Karras, T.; Aittala, M.; Aila, T. Noise2Noise: Learning image restoration without clean data. _arXiv preprint arXiv:1803.04189_ 2018,
  * Krull et al. 2019 Krull, A.; Buchholz, T.-O.; Jug, F. Noise2void-learning denoising from single noisy images. Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. 2019; pp 2129–2137
  * Li et al. 2023 Li, X.; Li, Y.; Zhou, Y.; Wu, J.; Zhao, Z.; Fan, J.; Deng, F.; Wu, Z.; Xiao, G.; He, J.; others Real-time denoising enables high-sensitivity fluorescence time-lapse imaging beyond the shot-noise limit. _Nature Biotechnology_ 2023, _41_ , 282–292
  * Kuwahara et al. 2022 Kuwahara, M.; Mizuno, L.; Yokoi, R.; Morishita, H.; Ishida, T.; Saitoh, K.; Tanaka, N.; Kuwahara, S.; Agemura, T. Transient electron energy-loss spectroscopy of optically stimulated gold nanoparticles using picosecond pulsed electron beam. _Applied Physics Letters_ 2022, _121_


# Integrability of a Family of Lotka–Volterra Three Species Biological System

Aween Karim1, Azad Amen1,2,3,4 and Waleed Aziz5 1 Department of Mathematics,
Collage of Basic Education, Salahaddin University–Erbil, Kurdistan region–Iraq
2 Department of Mathematics, Basic Education College, Raparin
University–Ranya, Kurdistan region–Iraq  3 Department of Mathematics, Faculty
of Science, Soran University–Soran, Kurdistan region–Iraq  4 Department of
Mathematics, Faculty of Science, Duhok University–Duhok, Kurdistan region–Iraq
5 Department of Mathematics, College of Science, Salahaddin University–Erbil,
Kurdistan region–Iraq<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

The aim of this study is to analyze the integrability problem of
Lotka–Volterra three species biological system. The system which considered in
this work is a biological plausibility or a chemical model. The system has a
complex dynamical behavior because it is chaotic system. We, first show that
the system is a complete integrable when two of the involved parameters in the
system are zero. Second, thorough invariant algebraic surfaces and exponential
factors, the nonintegrability problems have been investigated. Particularly,
we show the non-existence of polynomial, rational, formal series, and Darboux
first integrals when parameters are strictly positive.

###### Key words and phrases:

Integrability, Invariant algebraic surface, Darboux first integral, Formal
first integral.

###### 2010 Mathematics Subject Classification:

Primary 37K10 Secondary 34A05

##  1\. Introduction

During the last fifty years there has been increasing interest in studying the
autonomous differential system, mainly due to their many applications in
natural science. The conservation of ecological and biological systems is a
primary concern for scientists and researchers, and controlling and analyzing
the complex dynamical behavior of these systems is a challenge. The predation
and competition species are the most popular interactions in which nonlinear
polynomials are involved to represent such kind of interactions. The Lotka-
Voltera system is one of the most prominent among existing these models. The
authors in [1] have been demonstrated that three species can possess chaotic
behavior in an ecosystem. Samardzija and Greller [15] derived a three species
Lotka–Volterra biological system, which can be describe by a three dimensinal
system of differential equations

$\displaystyle\dot{x}$ $\displaystyle=x(1-y+cx-axz),$ (1)
$\displaystyle\dot{y}$ $\displaystyle=y(-1+x),$ $\displaystyle\dot{z}$
$\displaystyle=z(-b+ax^{2}).$

In (1), $x$ is the prey population and $y$, $z$ are predator populations, and
by the biological meaning of the system, we only consider the parameters that
are assumed to be non negative real numbers. For a certain value of
parameters, $a=2.9851,b=3,$ and $c=2,$ the authors in [15] showed that the
three-species biological system (1) was chaotic and it had very complicated
dynamical behavior. In 1999, Costello [13] studied chaos synchronization in
the integer-order system (1) for a certain parametric. Elsadany et al. [3]
studied the dynamical behaviors of the system (1). They investigated the
boundedness, existence, and uniqueness of the solutions of the system (1) and
determined the stability and bifurcation of its equilibrium points.

One of the main difficulties in studying these differential systems consists
of controlling the existence and nonexistence of first integrals and complete
integrability. Despite these earlier studies, the integrability of the system
(1) has not yet been investigated. The study of the integrability, the
existence of first integrals, is a main problem in the theory of differential
equations. The existence of first integrals of a system are crucial in
understanding the dynamics of the system. Therefore, it is important to
understand whether a system has first integrals and whether they are
analytical, smooth, etc. We study the existence of first integrals of system
(1), that described either by a formal power series (see [14], chapter one),
rational function or by Darboux function using the Darboux theory of
integrability (see [7], [5] and [19]). Therefore, we are interested in
studying the integrability of system (1) completely and showing that it is
non-integrable for most values of the parameters $a,b$ and $c$.
First, we look at the formal first integrals as it is a classical tool in
studing differential equations. It is alos been used to described solutions
around singularities [2], the existence of first integrals given by formal
series [17], Moussu [16]. The most success in using formal series to study
differential equations has been achieved by Ecalle [8], who used them to prove
the Dulac’s conjecture. Mattei and Moussu [12] are proved that formal
integrability implies analytic integrability in dimension two.

The associated vector field of system (1) is

$\mathcal{X}=x(1-y+cx-axz)\frac{\partial}{\partial
x}+y(-1+x)\frac{\partial}{\partial y}+z(-b+ax^{2})\frac{\partial}{\partial
z}.\\\ $ (2)

Let $U\subseteq\mathbb{C}^{3}$ be an open set. A non-constant function
$H:U\to\mathbb{C}$ is a first integral of the polynomial vector field
$\mathcal{X}$ on $U$ if it stays constant along all solution curves
$(x(t),y(t),z(t))$ of (1). Clearly $H$ is a first integral of $\mathcal{X}$ on
$U$ if and only if $\mathcal{X}H=0$ on $U$. When $U=\mathbb{C}^{3},$ the first
integral $H$ is called a global first integral. When a first integral $H$ is a
rational (polynomial or analytic) function, we say that $H$ is a rational
(polynomial or analytic) first integral. when $H$ is a formal power series in
the variables $x,y$ and $z$, we say that $H$ is a formal first integral.
Finally, a first integral is of Darboux type if it is of the form

$f_{1}^{\lambda_{1}}\dots f_{p}^{\lambda_{p}}E_{1}^{\mu_{1}}\dots
E_{q}^{\mu_{q}},$ (3)

where $f_{1},\dots,f_{p}$ are Darboux polynomials and $E_{1},\dots,E_{q}$ are
exponential factors (see Section 2 for definitions), and
$\lambda_{i},\mu_{j}\in\mathbb{C}$ for $i=1,\dots,p,$ and $j=1,\dots,q.$ The
functions of the form (3) are called Darboux functions, and they are the base
of the Darboux theory of integrability. The Darboux theory of integrability is
an algebraic theory of integrability which based on the an adequate number of
invariant algebraic surfaces and exponential factors associated to polynomial
differential systems. In fact, to every Darboux polynomial there is associated
some invariant algebraic surface, and the exponential factors appear when an
invariant algebraic surface has multiplicity larger than $1$, see [19, 6, 11]
for more details. Several well-known findings on the Darboux theory of
integrability and analytical first integrals can be found in [19, 18, 4]. The
following is our first main results.

######  Theorem 1.1.

The unique irreducible Darboux polynomials of system (1) with non-zero
cofactors are $x$ and $y$ and $z$.

This result is the basis of the Darboux theory of integrability, and its proof
is given in Section 4.

######  Theorem 1.2.

The following statements hold for system (1).

  1. 1)

The unique exponential factors of system (1) are $\exp(x+z)$ and $\exp(y)$
with cofactors $cx^{2}-xy-bz+x$ and $y(x-1)$ respectively if $a>0$ and $c>0,$
with an exception if $c=0,$ an extra exponential factor $\exp((x+y+z)^{2})$
can be appear with cofactor $-2(x+y+z)(bz-x-y)$.

  2. 2)

If $a=0,b>0$ and $c>0$, then system (1) has an exponential factor $\exp(z)$
with cofactor $-bz$ with an exception if $c=0,$ an extra exponential factor
$\exp(x+y)$ can be derived with cofactor $x-y.$ Moreover, if $ab=0,$ and
$c>0,$, then it has no exponential factors.

The proof of Theorem 1.2 is given in Section 5. The next two results state the
existence of Darboux first integrals for a set of the values of the parameters
and the absence of formal first integrals for system (1).

######  Theorem 1.3.

If $a=c=0$, then for system (1) the following statements hold.

  1. i)

If $b=0,$ then it is Darboux integrable with the first integrals

$H_{1}=xy\exp(-x-y)\quad\mbox{and}\quad H_{2}=z.$

  2. ii)

If $b>0,$ then it is integrable with the following first integrals

$H_{1}=xy\exp(-x-y)$ and $H_{2}=zh(x,y),$

where $h(x,y)=\exp(\int^{x}\\!{\frac{b}{{\it s\left({\rm
LambertW}\left(-{\frac{xy{{\rm e}^{{\it s}-x-y}}}{{\it
s}}}\right)+1\right)}}}d{\it s}),$ and lambertW computes the principal value
of the Lambert W function.

######  Theorem 1.4.

System (1) with $a>0$ and $c>0$ has no first integrals of Darboux type.

######  Theorem 1.5.

If $a>0$ and $c>0$, then system (1) has no rational first integrals.

Darboux first integrals are not necessarily formal series first integrals, and
vice versa. This demonstrates the independence of these topics.

######  Theorem 1.6.

When $a>0,c>0$, system (1) does not admits any formal first integral.

To prove Theorem 1.6, when $b>0$, we rewrite system (1) as a four dimensional
system in the variables $x,y,z$, and $b$

$\dot{x}=x(1-y+cx-
axz),\quad\dot{y}=y(-1+x),\quad\dot{z}=z(-b+ax^{2}),\quad\dot{b}=0.$ (4)

######  Theorem 1.7.

Suppose that $a>0,$ $c>0,$ and $b>0.$ If $f=f(x,y,z,b)$ is a formal first
integral of system (4), then $f$ is an arbitrary formal power series in the
variable $b$.

From Theorem 1.6, directly the following result for the system (1) is
obtained.

######  Corollary 1.

Suppose that $a>0,$ and $c>0.$ Then, the system (1) has no global analytic
first integrals, no polynomial first integrals and no local analytic first
integrals at the origin.

The proof of Theorems 1.3, 1.6, 1.7, 1.4 and 1.5 are given in Section 6.

Since the Darboux theory of integrability of a polynomial differential system
is based on the existence of Darboux polynomials and their multiplicity, then
the study of the existence or non-existence of Darboux first integrals needs
to seek of the Darboux polynomials. We shall recall the main concepts of the
Darboux theory of integrability.

##  2\. Preliminaries

A polynomial $f(x,y,z)\in\mathbb{C}[x,y,z]$ satisfying the equation

$x(1-y+cx-axz)\frac{\partial f}{\partial x}+y(-1+x)\frac{\partial f}{\partial
y}+z(-b+ax^{2})\frac{\partial f}{\partial z}=Kf,$ (5)

is called a Darboux polynomial of the system (1), where
$k=k(x,y,z)\in\mathbb{C}[x,y,z],$ is a cofactor of $f(x,y,z)$ of degree at
most two. Therefore, without loss of generality, the cofactor can be in the
form

$K=\alpha_{0}+\alpha_{1}x+\alpha_{2}y+\alpha_{3}z+\alpha_{4}x^{2}+\alpha_{5}xy+\alpha_{6}xz+\alpha_{7}y^{2}+\alpha_{8}yz+\alpha_{9}z^{2},$
(6)

where $\alpha_{i}\in\mathbb{C}$ for $i=0,\dots,9$. If $f(x,y,z)$ is a Darboux
polynomial, then the surface $f(x,y,z)=0$ is an invariant algebraic surface.

We recall the following result in [6].

######  Lemma 2.1.

Let $f$ be a polynomial and $f=\prod_{i=1}^{s}f_{j}^{\alpha_{j}}$ be its
decomposition into irreducible factors in $\mathbb{C}[x,y,z].$ Then $f$ is a
Darboux polynomial of system (1) if and only if all the $f_{j}$ are Darboux
polynomials of system (1). Moreover if $K$ and $K_{j}$ are the cofactors of
$f$ and $f_{j},$ then $K=\sum{{}_{j=1}^{s}}\alpha_{j}K_{j}.$

An exponential factor $E$ of system (1) is a function of the form
$E=\exp(\frac{h}{f})$ where $h,f\in\mathbb{C}[x,y,z]$ satisfying $(h,f)=1$ and

$\mathcal{X}E=LE,$ (7)

for some polynomial $L=L(x,y,z)$ of degree at most two. Such function is
called the cofactor of $E$.

######  Proposition 1.

[10, 9]. The following statements hold.

  1. 1)

If $\exp({\frac{g}{h}})$ is an exponential factor for the polynomial
differential system (1) and $h$ is not a constant polynomial, then $h=0$ is an
invariant algebraic surface.

  2. 2)

Eventually $\exp({g})$ can be an exponential factor, coming from the
multiplicity of the infinite invariant plane.

The following is well known result of the Darboux theory of integrability, for
instance see Chapter 3 of [19].

######  Theorem 2.2 (Darboux Theory of Integrability).

Suppose that a polynomial vector field $\mathcal{X}$ defined in
$\mathbb{R}^{n}$ of degree $m$ admits $p$ Darboux polynomials $f_{i}$ with
cofactors $K_{i}$ for $i=1,...,p,$ and $q$ exponential factors
$E_{j}=\exp(g_{j}/h_{j})$ with cofactors $L_{j}$ for $j=1,...,q.$ If there
exist $\lambda_{i},\mu_{j}\in\mathbb{C}$ not all zero such that

$\sum{{}_{i=1}^{p}}\lambda_{i}K_{i}+\sum{{}_{j=1}^{q}}\mu_{j}L_{j}=0,$ (8)

then the following real (multivalued) function is Darboux type

$f_{1}^{\lambda_{1}}...f_{p}^{\lambda_{p}}E_{1}^{\mu_{1}}...E_{q}^{\mu_{q}},$
(9)

and is also a first integral of the vector field $\mathcal{X}$.

The following result is well known, and it can be prove easily by Darboux
theory of integrability.

######  Lemma 2.3.

If system (1) has a rational first integral, then either it has a polynomial
first integral or two Darboux polynomials with the same nonzero cofactor.

##  3\. Auxiliary Results

In this section, we state and prove some auxiliary results that will be used
throughout the paper. These results study the existence of Darboux polynomials
with a non-zero cofactor in two dimensions by restricting the original system
to $y=0$ or $z=0$. It also examines the formal first integral of the system
restricted to $z=0$.

###  3.1. The system restricted to y=0

The system (1) restricted to $y=0$ becomes

$\begin{split}&\dot{x}=x(1+cx-axz),\\\ &\dot{z}=z(-b+ax^{2}).\end{split}$ (10)

Let $\mathcal{Y}=\mathcal{X}|_{y=0}$ be the vector field associated system
(LABEL:system_2),

$\mathcal{Y}=x(1+cx-axz)\frac{\partial}{\partial
x}+z(-b+ax^{2})\frac{\partial}{\partial z}.$

Then the following results can be obtained.

######  Lemma 3.1.

For system (LABEL:system_2), the following statements hold.

  1. i)

If $a>0,$ the unique irreducible Darboux polynomials with non-zero cofactors
are $x$ and $z$.

  2. ii)

If $a=0,$ the unique irreducible Darboux polynomials with non-zero cofactors
are $x,z,$ and $1+cx$.

Proof. It can be easily verified that $x$ and $z$ are Darboux polynomials with
respective cofactors, $1+cx-axz$ and $-b+ax^{2}$. First, suppose $a>0$, and we
show that there is no other Darboux polynomial. Suppose that $f$ is an
irreducible Darboux polynomial of system (LABEL:system_2) with degree greater
than or equal $2$ with a non-zero cofactor
$k_{1}=\beta_{0}+\beta_{1}x+\beta_{2}z+\beta_{3}x^{2}+\beta_{4}xz+\beta_{5}z^{2}$.
Then $f$ must stisfies

$x(1+cx-axz)\frac{\partial f}{\partial x}+z(-b+ax^{2})\frac{\partial
f}{\partial z}=k_{1}f.$ (11)

By restricting (11) to the invariant plane $x=0$ and denoting $f$ by
$\bar{g}$, we obtain

$-bz\frac{d\bar{g}}{dz}=(\beta_{0}+\beta_{2}z+\beta_{5}z^{2})\bar{g}.$ (12)

If $b=0$, we see that $(\beta_{0}+\beta_{2}z+\beta_{5}z^{2})\bar{g}=0$ which
implies $\beta_{0}=\beta_{2}=\beta_{5}=0$. When $b\neq 0$, the solution of
(12) is

$\bar{g}={\it\
\tilde{d_{1}}}\,{z}^{-{\frac{{\it\beta_{0}}}{b}}}{\exp({\,{\frac{-z\left({\it\beta_{5}}\,z+2\,{\it\beta_{2}}\right)}{2b}}}}),\quad\tilde{d_{1}}\in{\mathbb{C}}.$

Given that $\bar{g}$ is a polynomial, then must be $\beta_{2}=\beta_{5}=0$ and
$\beta_{0}=-bm_{0},$ where $m_{0}\in\mathbb{N}\cup\left\\{0\right\\}$. Hence
$k_{1}=-m_{0}b+\beta_{1}x+\beta_{3}x^{2}+\beta_{4}xz$ and
$f=\tilde{d_{1}}z^{m_{0}}+xh(x,z)$, where $h(x,z)$ is a polynomial in the
variables $x$ and $z$. It is obvious that $\tilde{d_{1}}\neq 0$ because $f$ is
irreducible.
Now, restricting (11) to the invariant plane $z=0$ and denoting $f$ by
$\tilde{g}$, we obtain

$x(1+cx)\frac{d\tilde{g}}{dx}=(-m_{0}b+\beta_{1}x+\beta_{3}x^{2})\tilde{g}.$
(13)

Here we consider two distinct cases.

Case 1: Suppose that $c>0$, then

$\tilde{g}={\it\tilde{d}_{2}}\,\left(1+cx\right)^{{\it
m_{0}}\,b+{\frac{{\it\beta_{1}}}{c}}-{\frac{{\it\beta_{3}}}{{c}^{2}}}}{x}^{-{\it
m_{0}}\,b}{\exp({{\frac{{\it\beta_{3}}\,x}{c}}}}),\quad\tilde{d_{2}}\in{\mathbb{C}}.$

Since $\tilde{g}$ is a polynomial, the must be $\beta_{3}=0$ and
$\beta_{1}=c(m_{1}-m_{0}b)$, with $m_{1}\in\mathbb{N}\cup\left\\{0\right\\}$.
Therefore,

$k_{1}=-m_{0}b+c(m_{1}-m_{0}b)x+\beta_{4}xz,$ (14)

and $f={\it\tilde{d}_{2}}\,\left(cx+1\right)^{m_{1}}{x}^{-{\it
m_{0}}\,b}+zT(x,z)$, where $T(x,z)$ is a polynomial, and since $f$ is
irreduced, then must $\tilde{d_{2}}\neq 0$.
We write $f=\sum{{}_{j=0}^{n}}f_{j}{(x,z)}$, where each $f_{j}=f_{j}(x,z)$
denotes a homogeneous polynomial of degree $j$ in $x$ and $z$. Obviously,
$f_{n}\neq 0$. We can deduce from $(\ref{eq4})$ and $(\ref{eq7})$ that the
terms of degree $n+2$ satisfy

$-ax^{2}z\frac{\partial f_{n}}{\partial x}+ax^{2}z\frac{\partial
f_{n}}{\partial z}=\beta_{4}xzf_{n}.$

Solving this linear partial differential equation, we obtain
$f_{n}=G_{n}(x+z)x^{-\frac{\beta_{4}}{a}}$. where $G_{n}(x+z)$ is an arbitrary
function in terms of $x+z$. Since $f_{n}$ must be a homogeneous polynomial, we
get that $\beta_{4}=-am_{2}$ where $m_{2}\in\mathbb{N}\cup\\{0\\}$ and
$G_{n}(x+z)\in{\mathbb{C}[x,z]\backslash\\{0\\}}$. Note that $f_{n}\neq 0$
implies that $G_{n}(x+z)\neq 0$. Furthermore, because $f_{n}$ has degree $n$,
we can write

$f_{n}=\tilde{d_{3}}(x+z)^{n-m_{2}}x^{m_{2}},\quad\tilde{d_{3}}\in\mathbb{C}\backslash\\{0\\},$
(15)

and the cofactor becomes $k_{1}=-m_{0}b+c(m_{1}-m_{0}b)x-am_{2}xz$. Computing
the terms of degree $n+1$ in $(\ref{eq4})$, we obtain

$cx^{2}\frac{\partial f_{n}}{\partial x}-ax^{2}z\frac{\partial
f_{n-1}}{\partial x}+ax^{2}z\frac{\partial f_{n-1}}{\partial
z}=c(m_{1}-m_{0}b)xf_{n}-am_{2}xzf_{n-1}$

Hence

$\displaystyle f_{{n-1}}$ $\displaystyle=$
$\displaystyle{\frac{\tilde{d_{3}}c}{a}}{x}^{{m_{2}}(x+z)^{n-{m_{2}}-1}}((-{m_{0}}\,b+{m_{1}}-n)\ln(-z)+\ln(x)({m_{0}}\,b-{m_{1}}+{m_{2}}))$
$\displaystyle+x^{m_{2}}{G_{n-1}}(x+z),$

where $G_{n-1}(x+z)$ is an arbitrary function in terms of $x+z$. Since
$f_{n-1}$ must be a homogeneous polynomial, we obtain
$G_{n-1}(x+z)\in\mathbb{C}[x,z]\backslash\\{0\\}$ and $-m_{0}b+m_{1}-n=0$, and
$m_{0}b-m_{1}+m_{2}=0$. As a result, $n=m_{2}$, and therefore

$f_{n-1}=x^{n}G_{n-1}(x+z),$

which implies that $G_{n-1}(x+z)=0$ and hence $f_{n-1}=0$. Since
$\beta_{4}=-am_{2}=-a(m_{1}-m_{0}b)$, and computing the terms of degree $n-1$
in (11), that satisfy

$x\frac{\partial f_{n}}{\partial x}-ax^{2}z\frac{\partial f_{n-2}}{\partial
x}+ax^{2}z\frac{\partial f_{n-2}}{\partial z}-bz\frac{\partial f_{n}}{\partial
z}=-a(m_{1}-m_{0}b)xzf_{n-2}-bm_{0}f_{n}.$ (16)

The function

$f_{{n-2}}\left(x,z\right)={\frac{{x}^{-{\it m_{0}}\,b+{\it
m_{1}}}\tilde{d_{3}}{\it
m_{1}}\,\left(\ln\left(x\right)x-x\ln\left(-z\right)-x-z\right)}{a\left(x+z\right)^{2}x}}+{x}^{-{\it
m_{0}}\,b+{\it m_{1}}}{\it G_{n-2}}\left(x+z\right),$

satisfy (16) where $G_{n-2}(x+z)$ is an arbitrary function in terms of $x+z$.
It obvious $\tilde{d_{3}}\neq 0$, otherwise, $f$ becomes reducible. Since
$f_{n-2}$ must be a homogeneous polynomial, we must have $m_{1}=0$. Hence,
$m_{2}=-m_{0}b$, and equation (15) becomes

$f_{n}=\tilde{d_{3}}x^{-m_{0}b},\quad\tilde{d_{3}}\in{\mathbb{C}\backslash\\{0\\}}.$

Since $f$ is a homogeneous polynomial and $m_{0}b\geq 0$, $b\neq 0$, then it
must be $m_{0}=0$, and hence $m_{2}=0$. As a result,
$\beta_{0}=\beta_{1}=\beta_{4}=0$, which implies that the cofactor $k_{1}=0$.

Case 2: When $c=0$. In this case, equation (13) becomes

$x\frac{\partial\tilde{g}}{\partial
x}=(-m_{0}b+\beta_{1}x+\beta_{3}x^{2})\tilde{g}.$ (17)

By solving it, we get $\tilde{g}={\it\tilde{d_{1}}}\,{x}^{-{\it
m0}\,b}{\exp({\frac{1}{2}x\left({\it\beta_{3}}\,x+2\,{\it\beta_{1}}\right)}})$,
where $\tilde{d_{1}}\in{\mathbb{C}\backslash\\{0\\}}$. We must have
$\beta_{1}=\beta_{3}=0$, so

$k_{1}=-m_{0}b+\beta_{4}xz.$ (18)

When $c=0$, from (18) and (11), one can see

$x(1-axz)\frac{\partial f}{\partial x}+z(-b+ax^{2})\frac{\partial f}{\partial
z}=(-m_{0}b+\beta_{4}xz)f.$ (19)

Writing $f=\sum{{}_{j=0}^{n}}f_{j}{(x,z)}$ as a homogeneous polynomial and
repeating the same way used in Case 1, we obtain

$f_{n}=\tilde{d_{2}}(x+z)^{n-m_{2}}x^{m_{2}},\qquad\tilde{d_{2}}\in{\mathbb{C}\backslash\\{0\\}},$
(20)

and

$f_{{n-1}}=x^{m_{2}}{G_{n-1}}\left(x+z\right),$

where $G_{n-1}(x+z)$ is a homogeneous polynomial in the terms of $(x+z)$. Then
we can write

$f_{n-1}=\tilde{d_{3}}(x+z)^{n-m_{2}-1}x^{m_{2}},\qquad\tilde{d_{3}}\in{\mathbb{C}\backslash\\{0\\}}.$

Finally, we calculate the terms of degree $n$ in equation (19), and we find
that

$\displaystyle f_{{n-2}}$ $\displaystyle=$
$\displaystyle{\frac{\tilde{d_{2}}}{a}}{x}^{-1+{m_{2}}}(x+z)^{-2+n-{m_{2}}}(-x({m_{0}}\,b+n)\ln(-z)+x({m_{0}}\,b+n)\ln(x)$
$\displaystyle-(x+z)({m_{0}}\,b+{m_{2}}))+x^{m_{2}}{G_{n-2}}(x+z).$

Since $f_{n-2}$ must be a homogeneous polynomial and also since
$\tilde{d_{2}}\neq 0$, then must be $m_{0}b+n=0$. We have $m_{0}b\geq 0$,
$b\neq 0$ and $n$ is nonnegative, then must be $m_{0}=0$ and so $n=0$. If
$n=0$, from (20), we get $f_{0}=\tilde{d_{2}}(x+z)^{-m_{2}}x^{m_{2}}$. Since
the degree of $f_{0}$ is zero, it must also $m_{2}=0$. Based on these values
of parameters, we get that $\beta_{0}=\beta_{4}=0$, and so $k_{1}=0$. This
concludes the proof of statement (i).

Now considering that system (LABEL:system_2) with $a=0$ has the rational first
integral $H=\frac{zx^{b}}{(1+cx)^{b}}$, Proposition (2.3) implies that all
invariant algebraic curves are in the set
$\\{(x,z)\in{\mathbb{R}^{2}}|Czx^{b}-(1+cx)^{b}=0,\>C\in\mathbb{R}\\}\backslash\\{x=-1/c\\}$.
This concludes the proof of the Lemma.∎

###  3.2. The system restricted to z=0

If system (1) restricted to $z=0$, then it becomes

$\begin{split}\dot{x}&=x(1-y+cx),\\\ \dot{y}&=y(-1+x),\\\ \end{split}$ (21)

The associated vector field of (21) is

$\mathcal{Z}=\mathcal{X}|_{z=0}=x(1-y+cx)\frac{\partial}{\partial
x}+y(-1+x)\frac{\partial}{\partial y}.$

To prove the next lemma, we need the following definition. A polynomial
$F(x,y,z)$ is said to be weight homogeneous of degree $r\in\mathbb{N}$ with
respect to the weight exponent $s=(s_{1},s_{2},s_{3})$ for all
$\mu\in\mathbb{R}\backslash{0}$ if it satisfies

$F(\mu^{s_{1}}x,\mu^{s_{2}}y,\mu^{s_{3}}z)=\mu^{r}F(x,y,z).$

######  Lemma 3.2.

The unique irreducible Darboux polynomials of system (21) with non-zero
cofactors are $x$ and $y$.

Proof. Assume that $f$ is a Darboux polynomial of system (21) with the
cofactor $k_{2}$. Then it must satisfies

$x(1-y+cx)\frac{\partial f}{\partial x}+y(-1+x)\frac{\partial f}{\partial
y}=k_{2}f.$ (22)

It is clear that $k_{2}$ is a polynomial of degree at most one. Without loss
of generality, we set $k_{2}=\alpha_{0}+\alpha_{1}x+\alpha_{2}y$. It is also
obvious that $x$ and $y$ are Darboux polynomial with respective cofactors
$1-y+cx$ and $-1+x$. We now show that the system (21) has no other Darboux
polynomials. Suppose that $f$ is an irreducible Darboux polynomial of system
(21) of degree at least two, and we claim that
$\alpha_{0}=\alpha_{1}=\alpha_{2}=0$.
On the invariant plane $x=0$, the equation (22) reduces to

$-y\frac{\partial\bar{h}}{\partial y}=(\alpha_{0}+\alpha_{2}y)\bar{h},$

where $\bar{h}=f|_{x=0}$, and its solution is
$\bar{h}=d_{1}{exp({-{\it\alpha_{2}}\,y}}){y}^{-{\it\alpha_{0}}},$ where
$d_{1}\in{\mathbb{C}}$. Since $\bar{h}$ must be a polynomial, then
$\alpha_{2}$ must be zero, and $\alpha_{0}=-m_{0}$ where
$m_{0}\in\mathbb{N}\cup\left\\{0\right\\}$. It follows that
$k_{2}=-m_{0}+\alpha_{1}x$ and $f=d_{1}y^{-m_{0}}+xg_{1}(x,y)$, where
$g_{1}(x,y)$ is a polynomial in the variables $x$ and $y$. Note that
$d_{1}\neq 0$, otherwise $f$ would be reducible, which is a contradiction.
Now, if we set $\tilde{h}=f|_{y=0}$, the restriction (22) to $y=0$ is

$x(1+cx)\frac{\partial\tilde{h}}{\partial x}=(-m_{0}+\alpha_{1}x)\tilde{h}.$

Let $c\neq 0$, then
$\tilde{h}=d_{2}\left(cx+1\right)^{{\frac{{\it\alpha_{1}}}{c}}+{\it
m_{0}}}{x}^{-{\it m_{0}}}$, where $d_{2}\in{\mathbb{C}}$. Since $\tilde{h}$
must be a polynomial, we obtain $\alpha_{1}=c(m_{1}-m_{0})$, where
$m_{1}\in\mathbb{N}\cup\left\\{0\right\\}$. So $k_{2}=-m_{0}+c(m_{1}-m_{0})x$
and $f=d_{2}x^{-m_{0}}(cx+1)^{{\frac{{\it\alpha_{1}}}{c}+{\it
m_{0}}}}+yg_{2}(x,y)$, where $g_{2}(x,y)$ is a polynomial in the variables $x$
and $y$. Since $f$ is irreducible, then $d_{2}\neq 0$.
We consider now the change of variables $x=X,y=\mu{Y}$ with
$\mu\in{\mathbb{C}\backslash\\{0\\}}$. Then system (21) becomes

$\begin{split}&\dot{X}=X(1-\mu{Y}+cX),\\\ &\dot{Y}=Y(-1+X),\end{split}$ (23)

where the prime denotes derivative with respect to the variable $T$. We apply
the transformation above and put $F(X,Y)=\mu^{n}f(X,\mu{Y})$, where $n$ is the
highest weight degree in the weight homogeneous components of $f$ in $x$ and
$y$ with the weight $(0,1)$ and $k_{2}=-m_{0}+c(m_{1}-m_{0})X$. Then $F$
satisfies

$\frac{dF}{dT}=\frac{d\mu^{n}f}{dt}=\mu^{n}\frac{df}{dt}=\mu^{n}k_{2}f=k_{2}F.$

Suppose $F=F_{0}+\mu F_{1}+\mu^{2}F_{2}+\dots+\mu^{l}F_{l}$, where $F_{i}$ is
a weight homogeneous polynomial in $X$ and $Y$ with the weight degree $n-i$
for $i=0,1,\dots,l$, and $n\geq l$. Clearly $f=F|_{\mu=1}$. By definition of a
Darboux polynomial, we have

$X(1-\mu{Y}+cX)\sum{{}_{i=0}^{n}}\mu^{i}\frac{\partial F_{i}}{\partial
x}+Y(-1+X)\sum{{}_{i=0}^{n}}\mu^{i}\frac{\partial F_{i}}{\partial
y}=(-m_{0}+c(m_{1}-m_{0})X)\sum{{}_{i=0}^{n}}\mu^{i}F_{i}.$

We claim that $m_{1}=m_{0}=0$. By Equating the terms with $\mu^{i}$, for
$i=0,1$, we obtain

$\begin{split}&L(F_{0})=(-m_{0}+c(m_{1}-m_{0})X)F_{0},\\\
&L(F_{1})=(-m_{0}+c(m_{1}-m_{0})X)F_{1}+XY\frac{\partial F_{0}}{\partial
x},\end{split}$ (24)

where $L$ is a vector field

$L=X(1+cX)\frac{\partial}{\partial X}+Y(-1+X)\frac{\partial}{\partial Y}.$

Solving the first equation in (24) and pulling back its solution to the second
equation, we obtain

$F_{0}\left(X,Y\right)={\it
G_{0}}\left(YX\left(Xc+1\right)^{{\frac{-1-c}{c}}}\right)\left(Xc+1\right)^{{\it
m_{1}}}{X}^{-{\it m_{0}}},$ $\begin{split}F_{1}(X,Y)=\frac{(Xc+1)^{\it
m_{1}}X^{-\it
m_{0}}}{\Gamma({\frac{2c-1}{c}})\Gamma({\frac{c-1}{c}})}\Bigg{(}-D(G_{0})\left(YX(Xc+1)^{({\frac{c-1}{c}})}\right)X^{3}\Gamma({\frac{c-1}{c}})\\\
\Gamma({\frac{2c-1}{c}})Y^{2}(Xc+1)^{({\frac{-2c-2}{c}})}\mbox{hypergeom}([1,{\frac{2c-1}{c}}],[2],-Xc)c\\\
+M_{1}(Xc+1)^{({\frac{-c-1}{c}})}+\Gamma({\frac{c-1}{c}})\Big{(}{\frac{1}{2}}M_{2}D(G_{0})\left(YX(Xc+1)^{({\frac{-c-1}{c}})}\right)\\\
Y^{2}X(Xc+1)^{({\frac{-2c-2}{c}})}+G_{1}\left(YX(Xc+1)^{{\frac{-c-1}{c}}}\right)\Gamma({\frac{2c-1}{c}})\Big{)}\Bigg{)},\end{split}$

where

$\begin{split}M_{1}=G_{0}\big{(}YX\left(Xc+1\right)^{\frac{-c-1}{c}}\big{)}Y\Gamma({\frac{2c-1}{c}})\Bigg{(}c\Big{(}\Psi({\frac{2c-1}{c}})m_{0}+m_{0}\ln(X)+m_{0}\ln(c)\\\
+\mbox{hypergeom}\left([1,1,{\frac{2c-1}{c}}],[2,2],-Xc\right)\left(m_{0}-m_{1}\right)Xc+m_{0}(\Gamma-1)\Big{)}X\Gamma\left({\frac{2c-1}{c}}\right)\\\
-{\frac{1}{2}}m_{0}\mbox{hypergeom}\left([1,1,{\frac{3c-1}{c}}],[2,3],-Xc\right)\Bigg{)}\Gamma\left({\frac{2c-1}{c}}\right)X^{2}c^{2}\\\
-\Gamma\left({\frac{c-1}{c}}\right)\Bigg{(}Xc(m_{0}-m_{1})\Psi({\frac{c-1}{c}})+Xc\left(m_{0}-m_{1}\right)\ln(X)\\\
+Xc(m_{0}-m_{1})\ln(c)+c\Gamma\left(m_{0}-m_{1}\right)X-m_{0}\Bigg{)},\end{split}$

and

$\begin{split}M_{2}=\Bigg{(}\Big{(}2X(c-1)\Psi({\frac{2c-1}{c}})+2X(c-1)\ln(X)+2X(c-1)\ln(c)-2\\\
+2\Gamma(c-1)X\Big{)}\Gamma({\frac{2c-1}{c}})+\Big{(}(-2\Psi({\frac{3c-1}{c}})-2\ln(X)-2\ln(c)\\\
+(-2\mbox{hypergeom}([1,1,{\frac{3c-1}{c}}],[2,2],-Xc)c\\\
+2\mbox{hypergeom}([1,1,{\frac{3c-1}{c}}],[2,2],-Xc))X-2\Gamma+2)\Gamma({\frac{3c-1}{c}})\\\
+\Gamma({\frac{4c-1}{c}})Xc\>\mbox{hypergeom}([1,1,{\frac{4c-1}{c}}],[2,3],-Xc)\Big{)}cX\Bigg{)}.\\\
\end{split}$

where $G_{0},G_{1}\in{\mathbb{C}[X,Y]}$ and have weight degrees $n$ and $n-1$,
respectively. Since $F_{1}$ is a polynomial and $c\neq 0$, we must have
$m_{0}=0$ and $m_{1}-m_{0}=0$. This gives, $k_{2}=0$.
When $c=0$, the equation (22) takes the form

$x(1-y)\frac{\partial f}{\partial x}+y(-1+x)\frac{\partial f}{\partial
y}=-m_{0}f.$ (25)

Writing $f=\sum{{}_{j=0}^{n}f_{j}}$, where each $f_{j}=f_{j}(x,y)$ represents
a homogeneous polynomial of degree $j$ in $x$ and $y$, in which $f_{n}\neq 0$.
By computing the terms of degree $n+1$ in (25), we obtain

$-xy\frac{\partial f_{n}}{\partial x}+xy\frac{\partial f_{n}}{\partial y}=0,$

and whose solution is that is

$f_{n}(x,y)=T_{1}(x+y),$

where $T_{1}(x+y)$ is a homogeneous polynomial of degree $n$. The terms of
degree $n$ in (25) satisfies

$x\frac{\partial f_{n}}{\partial x}-xy\frac{\partial f_{n-1}}{\partial
x}-y\frac{\partial f_{n}}{\partial y}+xy\frac{\partial f_{n-1}}{\partial
y}=-m_{0}f_{n}.$

Solving the partial differential differential above, we get

$\begin{split}f_{{n-1}}\left(x,y\right)=&-\frac{1}{x+y}\left(\ln\left(x\right)+\ln\left(-y\right)\right)\left(x+y\right)\mbox{D}\left(T_{1}\right)\left(x+y\right)\\\
&-\ln\left(-y\right)m_{{0}}T_{1}\left(x+y\right)+\left(x+y\right)T_{2}\left(x+y\right)\\\
&+\ln\left(x\right)m_{{0}}T_{1}\left(x+y\right).\\\ \end{split}$ (26)

Since $f_{{n-1}}$ must be a polynomial, we obtain $m_{0}=0$, which lead to
$\alpha_{0}=0$, and so $k_{2}=0$. The proof now is complete.∎

######  Proposition 2.

The system (21) has no formal first integrals when $c>0$.

Proof. We assume that $f(x,y)$ is a formal first integral of system (21) with
$c\neq 0$. Then it satisfies

$x(1-y+cx)\frac{\partial f}{\partial x}+y(-1+x)\frac{\partial f}{\partial
y}=0.$ (27)

We write $f$ as $f=\sum{{}_{k\geq 0}}f_{k}(x)y^{k}$, where each $f_{k}(x)$ is
a formal power series in the variable $x$. Denoting the restriction of $f$ to
$y=0$ by $f_{0}=f_{0}(x)$ in equation (27). Then

$x(1+cx)\frac{df_{0}}{dx}=0,$

and its solution is $f_{0}=d_{0}$, where $d_{0}$ is a constant. Then we can
write $f=f_{0}+yg(x,y)=d_{0}+yg(x,y)$ where $g=\sum{{}_{k\geq
0}}f_{k+1}(x)y^{k}$. Now, the function $g$ must satisfies the equation

$x(1-y+cx)\frac{\partial g}{\partial x}+y(-1+x)\frac{\partial g}{\partial
y}=-(-1+x)g.$ (28)

It is suffices to show that

$f_{k+1}(x)=0\quad\mbox{for}\quad k\geq 0.$ (29)

Since the restriction of $g$ to $y=0$ is $f_{1}=f_{1}(x)$, then from (28), we
obtain

$x(1+cx)\frac{df_{1}}{dx}=-(-1+x)f_{1},$

and its solution is

$f_{1}=d_{1}x(1+cx)^{-1-\frac{1}{c}},$ (30)

where $d_{1}$ is a constant. By calculating the coefficient of $y$ in (28), we
see

$x(1+cx)\frac{df_{2}}{dx}-x\frac{df_{1}}{dx}=-2(-1+x)f_{2},$

or equivalently,

$x(1+cx)\frac{df_{2}}{dx}+d_{1}x(-1+x)(1+cx)^{-2-\frac{1}{c}}+2(-1+x)f_{2}=0.$
(31)

Multiplying (31) by $(1+cx)^{2+\frac{1}{c}}$, gives

$x(1+cx)^{3+\frac{1}{c}}\frac{df_{2}}{dx}+d_{1}x(-1+x)+2(-1+x)(1+cx)^{2+\frac{1}{c}}f_{2}=0.$
(32)

Evaluating (32) on $x=-\frac{1}{c}$, we obtain that
$-\frac{d_{1}(-1-\frac{1}{c})}{c}=0$, and since $c>0$, this implies $d_{1}=0$.
Then, from (30), we obtain $f_{1}=0$. We now assume that (29) is true for
$k=0,\dots,l-2$, and we will prove it is also true for $k=l-1$. By the
induction hypothesis, we have

$f_{1}=\dots=f_{l-2}=0\quad\mbox{and}\quad g=\sum{{}_{k\geq
0}}f_{k+1}(x)y^{k}=y^{l-1}\sum{{}_{k\geq 0}}f_{k+l}(x)y^{k}.$

By introducing $g$ in (28) and determining the coefficient of $y^{l-1}$, we
see that $f_{l}$ satisfies

$x(1+cx)\frac{df_{l}}{dx}+(-1+x)(l-1)f_{l}=-2(-1+x)f_{l},$

which yields

$f_{l}=d_{l}x^{l+1}(1+cx)^{-\frac{(1+c)(1+l)}{c}},$ (33)

where $d_{l}$ is a constant. By calculating the coefficient of $y^{l}$ in
(28), we obtain

$x(1+cx)\frac{df_{l+1}}{dx}-x\frac{df_{l}}{dx}+l(-1+x)f_{l+1}=-2(-1+x)f_{l+1},$

or equivalently,

$\begin{split}d_{{l}}{x}^{l+1}\left(-1+x\right)\left(l+1\right)\left(cx+1\right)^{-{\frac{\left(1+c\right)\left(l+1\right)}{c}}-1}+\left(c{x}^{2}+x\right){\frac{\rm
d}{{\rm d}x}}f_{{l+1}}\left(x\right)\\\
+f_{{l+1}}\left(x\right)\left(-1+x\right)\left(l+2\right)=0.\\\ \end{split}$
(34)

Multiplying (34) by $(1+cx)^{1+\frac{(1+c)(l+1)}{c}}$, implies

$\begin{split}\left(\left(c{x}^{2}+x\right){\frac{\rm d}{{\rm
d}x}}f_{{l+1}}\left(x\right)+f_{{l+1}}\left(x\right)\left(-1+x\right)\left(l+2\right)\right)\left(cx+1\right)^{{\frac{lc+l+2\,c+1}{c}}}\\\
+{x}^{l+1}d_{{l}}\left(-1+x\right)\left(l+1\right)=0.\\\ \end{split}$ (35)

Evaluating (35) on $x=-\frac{1}{c}$, we obtain that
$-\frac{d_{l}(l+1)(-\frac{1}{c})^{l}(-1-\frac{1}{c})}{c}=0$, and since $c>0$,
this implies $d_{l}=0$. Then, from (33), we obtain $f_{l}=0$, for $k=l-1$. So
we have by induction that $f_{k}=0$ for all $k\geq 1$, which yields $g=0$.
Consequently, $f=d_{0}$, and hence the system (21) has no formal first
integrals.∎

##  4\. Darboux polynomials with non-zero cofactors

It is obvious that $x$,$y$, and $z$ are Darboux polynomials of (1). We show
that system (1) has no any Darboux polynomial of degrees greater than one. The
proof of Theorem 1.1 will follows by the following Lemmas.

Proof of Theorem 1.1. We consider two cases.

Case 1. When $a=0$, the system (1) reduces to

$\begin{split}&\dot{x}=x(1-y+cx),\\\ &\dot{y}=y(-1+x),\\\ &\dot{z}=-bz,\\\
\end{split}$ (36)

then any cofactor of (36) must be the form

$k=\alpha_{0}+\alpha_{1}x+\alpha_{2}y+\alpha_{3}z.$ (37)

######  Lemma 4.1.

Let $f$ be an irreducible Darboux polynomial of degree greater than one with
non-zero cofactor (37). Then $\alpha_{0}=\alpha_{1}=\alpha_{2}=\alpha_{3}=0$.

Proof. Assume that $f$ is a Darboux polynomial of degree $n\geq 2$ in system
(36) with a non-zero cofactor $k$. Then $f$ satisfies

$x(1-y+cx)\frac{\partial f}{\partial x}+y(-1+x)\frac{\partial f}{\partial
y}-bz\frac{\partial f}{\partial z}=kf.$ (38)

Suppose that $\bar{f}$ is a restriction of $f$ to $x=0$, then from equation
(38), $\bar{f}$ satisfies

$-y\frac{\partial\bar{f}}{\partial y}-bz\frac{\partial\bar{f}}{\partial
z}=(\alpha_{0}+\alpha_{2}y+\alpha_{3}z)\bar{f},\\\ $ (39)

and $\bar{f}$ satisfies

$\bar{f}={y}^{-{\it\alpha_{0}}}{\it\
F1}\left(z{y}^{-b}\right){\exp({-{\it\alpha_{2}}\,y-{\frac{{\it\alpha_{3}}}{b}}z}}).\\\
$ (40)

Since $\bar{f}$ must be a polynomial, we get $\alpha_{2}=\alpha_{3}=0$. We now
expand $f$ in powers of the variable $z$ so that
$f=\sum{{}_{j=0}^{n}}f_{j}{z}^{j}$, where each $f_{j}=f_{j}(x,y)$ is a
polynomial in the variables $x$ and $y$. The Darbox polynomial of system (36)
with the restriction $z=0$ is given by $f=f_{0}$, and according to the
irreducibility condition for $f$, $f_{0}\neq 0$. Moreover, system (36)
restricted to $z=0$ yields system (21), and so Lemma 3.2 guarantees that
$\alpha_{0}=\alpha_{1}=0$. Therefore, $k=0$, which means that the system (36)
has no Darboux polynomials of greater than one with non-zero cofactor.∎

Case 2. $a>0$. To prove this case, we shall need the following Lemmas.

######  Lemma 4.2.

Let $f$ be an irreducible Darboux polynomial of greater than one with non-zero
cofactor as in (6). Then $\alpha_{5}=\alpha_{7}=\alpha_{8}=\alpha_{9}=0$, and
$\alpha_{4}=N_{4}a,\alpha_{6}=-N_{6}a$, where
$N_{4},N_{6}\in{\mathbb{N}\cup\\{0\\}}$.

Proof. Let $f(x,y,z)=\sum{{}_{i=0}^{n}}f_{i}{(x,y,z)}$ be an irreducible
Darboux polynomial of system (1) where $n\geq 2$, where each $f_{i}$ is a
homogeneous polynomial of degree $i$ for $i=0,1,\dots,n$. Clearly $f_{n}\neq
0$. Then $f$ satisfies

$\displaystyle x(1-y+cx-axz)\sum{{}_{i=0}^{n}}\frac{\partial f_{i}}{\partial
x}+y(-1+x)\sum{{}_{i=0}^{n}}\frac{\partial f_{i}}{\partial y}$
$\displaystyle+z(-b+ax^{2})\sum{{}_{i=0}^{n}}\frac{\partial f_{i}}{\partial
z}=K\sum{{}_{i=0}^{n}}f_{i}.$

Computing the homogeneous terms of degree $n+2$, we get

$-ax^{2}z\frac{\partial f_{n}}{\partial x}+ax^{2}z\frac{\partial
f_{n}}{\partial
z}=(\alpha_{4}x^{2}+\alpha_{5}xy+\alpha_{6}xz+\alpha_{7}y^{2}+\alpha_{8}yz+\alpha_{9}z^{2})f_{n}.$

The solution this partial differential equation above is

$\displaystyle
f_{n}=T_{n}(y,x+z)(-z)^{{\frac{{\it\alpha_{4}}}{a}}+{\frac{{\it\alpha_{5}}\,y}{a\left(x+z\right)}}+{\frac{{\it\alpha_{7}}\,{y}^{2}}{a\left(x+z\right)^{2}}}}{x}^{-{\frac{{\it\alpha_{5}}\,y}{a\left(x+z\right)}}-{\frac{{\it\alpha_{6}}}{a}}-{\frac{{\it\alpha_{7}}\,{y}^{2}}{a\left(x+z\right)^{2}}}+{\frac{{\it\alpha_{9}}}{a}}}$
$\displaystyle{\exp({{\frac{{\it\alpha_{7}}\,{y}^{2}}{a\left(x+z\right)x}}+{\frac{{\it\alpha_{8}}\,y}{ax}}+{\frac{\left(x+z\right){\it\alpha_{9}}}{ax}}}}),$

where $T_{n}$ is an arbitrary function of $y$ and $x+z$. Since $f_{n}$ is a
homogeneous polynomial of degree $n$, we must have
$\alpha_{5}=\alpha_{7}=\alpha_{8}=\alpha_{9}=0,$ and
$\alpha_{4}=N_{4}a,\alpha_{6}=-N_{6}a$ where
$N_{4},N_{6}\in{\mathbb{N}\cup\\{0\\}}$.∎
Therefore,

$K=\alpha_{0}+\alpha_{1}x+\alpha_{2}y+\alpha_{3}z+N_{4}ax^{2}-N_{6}axz.$ (41)

######  Lemma 4.3.

Let $f$ be an irreducible Darboux polynomial of degree greater than one with
non-zero cofactor $K$ as in (41). Then
$\alpha_{0}=\alpha_{1}=\alpha_{3}=N_{4}=N_{6}=0.$

Proof. Let $f\in{\mathbb{C}[x,y,z]}$ be an irreducible Darboux polynomial of
system (1) of degree $n\geq 2$. Writing $f$ in powers of the variable $y$,
that is, $f=\sum{{}_{i=0}^{n}}f_{i}(x,z)y^{i}$, where $f_{i}=f_{i}(x,z)$ is a
polynomial in variables $x$ and $z$ for each $i$. Given that $f$ is
irreducible, we conclude that $f_{0}\neq 0$, otherwise, $f$ would be divisible
by $y$, which is a contradiction. Moreover, $f_{0}$ is a Darboux polynomial of
system (1) restricted to $y=0$. Then, $f_{0}$ satisfies

$x(1+cx-axz)\frac{\partial f_{0}}{\partial x}+z(-b+ax^{2})\frac{\partial
f_{0}}{\partial
z}=(\alpha_{0}+\alpha_{1}x+\alpha_{3}z+N_{4}ax^{2}-N_{6}axz)f_{0},$ (42)

which is (11) with
$\beta_{0}=\alpha_{0},\beta_{1}=\alpha_{1},\beta_{2}=\alpha_{3},\beta_{3}=N_{4}a,\beta_{4}=-N_{6}a,\beta_{5}=0$.
Proceeding exactly the same way as in Lemma 3.1, we get that
$\alpha_{0}=\alpha_{1}=\alpha_{3}=N_{4}=N_{6}=0$.∎

By Lemma 4.3, we obtain that $K=\alpha_{2}y$.

######  Lemma 4.4.

Let $f$ be an irreducible Darboux polynomial of greater than one with non-zero
cofactor $K=\alpha_{2}y$. Then $\alpha_{2}=0$.

Proof. Let $f$ be an irreducible Darboux polynomial of system (1) of degree
$n\geq 2$. Without loss of generality, let
$f=\sum{{}_{i=0}^{n}}f_{i}{(y,z)x^{i}}$, where each $f_{i}=f_{i}(y,z)$ are
polynomials in the variables $y,z$. It is clear that $f_{0}\neq 0$, otherwise,
$f$ would be divisible by $x$, which is a contradiction. Since $f_{0}$ is
independent of $x$, it is a Darboux polynomial of system (1) restricted to
$x=0$ that satisfies

$-y\frac{\partial f_{0}}{\partial y}-bz\frac{\partial f_{0}}{\partial
z}=(\alpha_{2}y)f_{0}.$ (43)

We write $f_{0}$ as the sum of its homogeneous parts as
$f_{0}=\sum{{}_{i=0}^{m}}f_{0,i}(y,z)$, where each $f_{0,i}=f_{0,i}(y,z)$ is a
homogeneous polynomial in its variables of degree $i$ and $0\leq m\leq n$.
Since $f_{0}\neq 0$, then $f_{0,m}\neq 0$. Computing the terms of degree $m+1$
in (43), we have $(\alpha_{2}y)f_{0,m}=0$, and since $f_{0,m}\neq 0$, so must
be $\alpha_{2}=0$. This completes the proof of the Theorem 1.1.∎

##  5\. Exponential Factors

Proof of Theorem 1.2. Let $E=\exp(\frac{h}{f})$ be an exponential factor of
the system (1) with cofactor $L$ where $h,f\in{\mathbb{C}[x,y,z]}$ with
$(h,f)=1$. Then, from the definition of exponential factor and in view of
Proposition 1, we have either $f$ is a constant, in this case we can take
$f=1$, or $f$ is a Darboux polynomial of system (1), in this case from Theorem
1.1, $E$ can be of the form $E=\exp(\frac{h}{x^{s_{1}}y^{s_{2}}z^{s_{3}}})$
with $s_{1},s_{2}$, and $s_{3}$ are non-negative integers and such that
$h\in{\mathbb{C}[x,y,z]}$ is coprime with $x,y,z$. Obviously, when
$s_{1}=s_{2}=s_{3}=0$, we are in the previous case.
We first prove that the system (1) has no exponential factors of the form
$\exp(\frac{h}{x^{s_{1}}y^{s_{2}}z^{s_{3}}})$. For this purpose, we assume
that $\exp(\frac{h}{x^{s_{1}}y^{s_{2}}z^{s_{3}}})$ is an exponential factor of
the system and that at least one of the $s_{1},s_{2}$, or $s_{3}$ is positive,
and then we get a contradiction. Clearly, by (7), $h$ satisfies,

$\dot{x}\frac{\partial h}{\partial x}+\dot{y}\frac{\partial h}{\partial
y}+\dot{z}\frac{\partial h}{\partial z}-h(\frac{\dot{x}s_{1}}{\partial
x}+\frac{\dot{y}s_{2}}{\partial y}+\frac{\dot{z}s_{3}}{\partial
z})=Lx^{s_{1}}y^{s_{2}}z^{s_{3}},$ (44)

We distinguish the following cases.

Case 1: When $s_{1}>{0}$. Evaluating (44) on $x=0$ and setting
$\bar{h}=h|_{x=0}$, then $\bar{h}$ satisfies

$-y\frac{\partial\bar{h}}{\partial y}-bz\frac{\partial\bar{h}}{\partial
z}=\bar{h}\big{(}{s_{1}}(1-y)-{s_{2}}-b{s_{3}}\big{)}.$

and clearly $\bar{h}\neq 0$ because $(h,x)=1$. The right hand side of the
equation above has one degree more than the left hand side, and this is a
contradiction.

Case 2: When $s_{1}=0$ and $s_{2}>{0}$. On $y=0$ and setting
$\tilde{h}=h|_{y=0}$, then from (44) it must satisfies

$x(1+cx-axz)\frac{\partial\tilde{h}}{\partial
x}+z(-b+ax^{2})\frac{\partial\tilde{h}}{\partial
z}=\tilde{h}\big{(}{s_{2}}(-1+x)+{s_{3}}(-b+ax^{2})\big{)}.$ (45)

Clearly, $\tilde{h}\neq 0$, since $(h,y)=1$. Moreover, we note that
$\tilde{h}$ is a Darboux polynomial of the vector field $\mathcal{Y}$. Then,
from Lemma 3.1, we can consider the following subcases below.

Subcase 2.1. If $a>0$. Then system (LABEL:system_2) has Darboux polynomials
$x$ and $z$, so we can write $\tilde{h}=Gx^{m_{1}}z^{m_{2}}$ with $G$ is a
constant and $m_{1}$, $m_{2}$ are non-negative integers such that
$m_{1}+m_{2}>0$. We substitute $\tilde{h}$ in equation (45), and we deduce

$-am_{1}xz+am_{2}x^{2}+cm_{1}x-bm_{1}x-bm_{2}+m_{1}=as_{3}x^{2}-bs_{3}+s_{2}x-s_{2}.$

Then clearly $am_{1}=0$, and then $m_{1}=0$. This implies that
$\tilde{h}(x,z)=\tilde{h}(z)=Gz^{m_{2}}$. By substituting $\tilde{h}(z)$ in
(45), and simplifying it, we can derive

$(ax^{2}-b)m_{2}=as_{3}x^{2}-bs_{3}+s_{2}x-s_{2}.$ (46)

It is also easy to clear that $s_{2}=0$, which is a contradiction.

Subcase 2.2. If $a=0$. In this case, $\tilde{h}$ must be of the form
$\tilde{h}=G(1+cx)^{bv}-x^{bv}z^{v}$, where $G$ is a constant and $v$ is a
non-negative integer. The cofactor of $\tilde{h}$ is $vbc\,x$, and from (45),
we obtain

${s_{2}}(-1+x)+{s_{3}}(-b)=vbc\,x,$

which implies $s_{2}=-bs_{3}\leq 0$, and this a contradiction.

Case 3: When $s_{1}=s_{2}=0$ and $s_{3}>{0}$. From equation (44) with $z=0$,
we get

$x(1-y+cx)\frac{\partial\hat{h}}{\partial
x}+y(-1+x)\frac{\partial\hat{h}}{\partial y}=\hat{h}{s_{3}}(-b+ax^{2}),$ (47)

where $\hat{h}=h|_{z=0}$ and $\hat{h}\neq{0}$ because $h$ is coprime with $z$.
We emphasize that $\hat{h}$ is a Darboux polynomial of the vector field
$\mathcal{Z}$ with a non-zero cofactor. Therefore, by Lemma 3.2, we get that
$\hat{h}=Gx^{m_{1}}y^{m_{2}}$ where $G$ as a constant and $m_{1}$ and $m_{2}$
as non-negative integers such that $m_{1}+m_{2}>0$. Then from (47), we see
that

$(c\,m_{1}+m_{2})x-m_{1}y+m_{1}-m_{2}=s_{3}(-b+ax^{2}).$ (48)

We consider two subcases depending on $a$ and $b$. Obviously, at least one of
$a>0$ or $b>0$ is required for the existence of the Darboux polynomial $z$.

Subcase 3.1. If $a>0$. Computing the coefficient of $x^{2}$ in (48), we get
$as_{3}=0$. Hence, $s_{3}=0$, which is a contradiction.

Subcase 3.2. If $a={0}$ and $b>0$. Then equation (48) implies

$(c\,m_{1}+m_{2})x-m_{1}y+m_{1}-m_{2}=-bs_{3}.$ (49)

Computing the coefficients of $y$ and $x$ in (49), we obtain $m_{1}=m_{2}=0$.
Consequently, $s_{3}=0$, which is a contradiction.

Case 4: $s_{1}=s_{2}=s_{3}={0}$. In this case, equation (44) becomes

$x(1-y+cx-axz)\frac{\partial h}{\partial x}+y(-1+x)\frac{\partial h}{\partial
y}+z(-b+ax^{2})\frac{\partial h}{\partial z}=L,$ (50)

where
$L=l_{0}+l_{1}x+l_{2}y+l_{3}z+l_{4}x^{2}+l_{5}xy+l_{6}xz+l_{7}y^{2}+l_{8}yz+l_{9}z^{2}$,
with $l_{i}\in{\mathbb{C}}$ for $i=0,\dots,9$. Here, we consider different
subcases depending on $a$ and $c$.

Subcase 4.1. When $a>0$ and $c>0$. We write $h$ in the form
$h=\sum{{}_{j=0}^{n}}h_{j}(x,y,z)$, where each $h_{j}$ is a homogeneous
polynomial of degree $j$ in the variables $x,y,z$. Firstly, assume that $n\geq
3$. The terms of degree $n+2$ in (50) are

$-ax^{2}z\frac{\partial h_{n}}{\partial x}+ax^{2}z\frac{\partial
h_{n}}{\partial z}=0,$

whose solution is $h_{n}=W_{n}(y,x+z)$, where $W_{n}$ is an arbitrary $C^{1}$
function. Using the fact that $h_{n}$ is a homogeneous polynomial of degree
$n$, we can write $h_{n}=\sum{{}_{j=0}^{n}}a_{j}y^{n-j}(x+z)^{j}$, where
$a_{j}\in{\mathbb{C}}$. The terms of degree $n+1$ in (50) satisfy

$x(-y+cx)\frac{\partial h_{n}}{\partial x}-ax^{2}z\frac{\partial
h_{n-1}}{\partial x}+xy\frac{\partial h_{n}}{\partial y}+ax^{2}z\frac{\partial
h_{n-1}}{\partial z}=0.$

or equivalently,

$\displaystyle-ax^{2}z\frac{\partial h_{n-1}}{\partial
x}+ax^{2}z\frac{\partial h_{n-1}}{\partial
z}+(cx^{2}-xy)\sum{{}_{j=0}^{n}}ja_{j}y^{n-j}(x+z)^{j-1}$
$\displaystyle+xy\sum{{}_{j=0}^{n}}(n-j)a_{j}y^{n-j-1}(x+z)^{j}=0,$

The function

$\displaystyle
h_{n-1}=\frac{1}{2a(x+z)}(-4\operatorname{arctanh}(\frac{x-z}{x+z})yB+A((2xc+2zc-4y)\operatorname{arctanh}(\frac{x-z}{x+z})$
$\displaystyle-\ln(-(x+z)xz)c(x+z)))+W_{n-1}(y,x+z),$

satisfy the partial differential equations above where
$A=\sum{{}_{j=0}^{n}}ja_{j}y^{n-j}(x+z)^{j-1}$,
$B=\sum{{}_{j=0}^{n}}a_{j}(-n+j)y^{n-j-1}(x+z)^{j}$, and $W_{n-1}$ is an
arbitrary function of $y$ and $x+z$. Since $h_{n-1}$ must be a homogeneous
polynomial of degree $n-1$, we must have $A=0$ and $B=0$, which implies that
$a_{j}=0$ for $j=0,1,\dots,n$, and thus
$h_{n}=\sum{{}_{j=0}^{n}}a_{j}y^{n-j}(x+z)^{j}=0$. This contradiction with the
fact that $h_{n}$ is a polynomial of degree $n\geq{3}$. Then we must have
$n\leq{2}$. Therefore

$h=h_{0}+h_{1}x+h_{2}y+h_{3}z+h_{4}x^{2}+h_{5}xy+h_{6}xz+h_{7}y^{2}+h_{8}yz+h_{9}z^{2}.$

Substituting $h$ in (50), we get that

$\displaystyle
l_{0}=0,l_{1}=h_{1},l_{2}=-h_{2},l_{3}=-bh_{3},l_{4}=ch_{1}+2h_{4},l_{5}=-h_{1}+h_{2},l_{6}=-bh_{6}+h_{6},$
$\displaystyle l_{7}=-2h_{7},l_{8}=-bh_{8}-h_{8},l_{9}=-2bh_{9},$

where
$h_{0}=h_{0},h_{1}=h_{1},h_{2}=h_{2},h_{3}=h_{1},h_{4}=h_{5}=h_{6}=h_{7}=h_{8}=h_{9}=0$.
Then $h=h_{0}+(x+z)h_{1}+yh_{2}$ with
$L=(cx^{2}+(-y+1)x-bz)h_{1}+y(x-1)h_{2}$.

Subcase 4.2. When $a>0$ and $c=0$. From equation (50), we have

$x(1-y-axz)\frac{\partial h}{\partial x}+y(-1+x)\frac{\partial h}{\partial
y}+z(-b+ax^{2})\frac{\partial h}{\partial z}=L.$ (51)

By proceeding in a similar way as above, we get

$h_{n-1}=\frac{y}{a(x+z)}(A+B)(\ln(x)-\ln(-z))+W_{n-1}(y,x+z),$

where
$A=\sum{{}_{j=0}^{n}}ja_{j}y^{n-j}(x+z)^{j-1},B=\sum{{}_{j=0}^{n}}a_{j}(-n+j)y^{n-j-1}(x+z)^{j}$,
and $W_{n-1}$ is an arbitrary polynomial in $y$ and $x+z$. Since $h_{n-1}$
just admits polynomial solutions, we must have $A+B=0$, and
$a_{j}=\frac{n!}{j!(n-j)!}A_{0}$ for $j=0,1,\dots,n$. Hence,
$h_{n}=\sum{{}_{j=0}^{n}}\frac{n!}{j!(n-j)!}A_{0}y^{n-j}(x+z)^{j}$ and
$h_{n-1}=W_{n-1}(y,x+z)$. Since $h_{n-1}$ must be a homogeneous polynomial of
degree $n-1$, we write $h_{n-1}=\sum{{}_{j=0}^{n}}b_{j}y^{n-j-1}(x+z)^{j}$,
where $b_{j}\in{\mathbb{C}}$. Now, computing the homogeneous part of degree
$n$ in (51) yields

$-ax^{2}z\frac{\partial h_{n-2}}{\partial x}-xy\frac{\partial
h_{n-1}}{\partial x}+x\frac{\partial h_{n}}{\partial x}-y\frac{\partial
h_{n}}{\partial y}+xy\frac{\partial h_{n-1}}{\partial y}-bz\frac{\partial
h_{n}}{\partial z}+ax^{2}z\frac{\partial h_{n-2}}{\partial z}=0,$

which is

$\begin{split}-ax^{2}z\frac{\partial h_{n-2}}{\partial
x}-xy\frac{\partial\sum{{}_{j=0}^{n}}b_{j}y^{n-j-1}(x+z)^{j}}{\partial
x}+x\frac{\partial\sum{{}_{j=0}^{n}}\frac{n!}{j!(n-j)!}A_{0}y^{n-j}(x+z)^{j}}{\partial
x}\\\
-y\frac{\partial\sum{{}_{j=0}^{n}}\frac{n!}{j!(n-j)!}A_{0}y^{n-j}(x+z)^{j}}{\partial
y}+xy\frac{\partial\sum{{}_{j=0}^{n}}b_{j}y^{n-j-1}(x+z)^{j}}{\partial y}\\\
-bz\frac{\partial\sum{{}_{j=0}^{n}}\frac{n!}{j!(n-j)!}A_{0}y^{n-j}(x+z)^{j}}{\partial
z}+ax^{2}z\frac{\partial h_{n-2}}{\partial z}=0.\\\ \end{split}$

Solving it, we get

$\begin{split}h_{n-2}=\frac{1}{ax(x+z)^{2}}\big{(}-2\operatorname{arctanh}(\frac{x-z}{x+z})xy(x+z)C-2\operatorname{arctanh}(\frac{x-z}{x+z})xy(x+z)D\\\
+\left(-x\left(x-y+z\right)\ln\left(-z\right)+x\left(x-y+z\right)\ln\left(x\right)+\left(x+z\right)\left(bx+bz+y\right)\right)\\\
nA_{{0}}\left(x+z+y\right)^{n-1}\big{)}+{\it W_{n-2}}\left(y,x+z\right),\\\
\end{split}$

where $C=\sum{{}_{j=0}^{n}}(-n+j+1)b_{j}y^{n-j-2}(x+z)^{j}$,
$D=\sum{{}_{j=0}^{n}}jb_{j}y^{n-j-1}(x+z)^{j-1}$, and $W_{n-2}$ is an
arbitrary polynomial in $y$ and $x+z$. Since $h_{n-2}$ must be a homogeneous
polynomial of degree $n-2$ and $n\geq 3$, then must be $A_{0}=0$, so
$h_{n}=0$, which is a contradiction with the fact that $h_{n}$ is a polynomial
of degree $n\geq{3}$. Then we must have $n\leq{2}$. Thus,
$h=h_{0}+h_{1}x+h_{2}y+h_{3}z+h_{4}x^{2}+h_{5}xy+h_{6}xz+h_{7}y^{2}+h_{8}yz+h_{9}z^{2}$.
By substituting $h$ in (51), we can obtain

$h=h_{0}+(x+y+z)^{2}h_{1}+(x+z)h_{2}+yh_{3},$

and

$L=-2(x+y+z)(bz-x+y)h_{1}+\big{(}(h_{3}-h_{2})y+h_{2}\big{)}x-bzh_{2}-yh_{3}.$

This completes the proof of the statement a.

Subcase 4.3. When $a={0}$ and $c>0$. In this case, we have

$x(1-y+cx)\frac{\partial h}{\partial x}+y(-1+x)\frac{\partial h}{\partial
y}+z(-b)\frac{\partial h}{\partial z}=l_{0}+l_{1}x+l_{2}y+l_{3}z,$ (52)

where $l_{i}\in{\mathbb{R}}$ for $i=0,\dots,3$. First, we assume that $b>0$.
We decompose $h$ as a sum of polynomials in the variable $x$ as
$h=\sum{{}_{j=0}^{n}}h_{j}(y,z)x^{j}$, where each $h_{j}\in{\mathbb{C}[y,z]}$.
Assume that $n\geq 1$, then the coefficient of $x^{n+1}$ in (52), is

$c\,n\,h_{n}+y\frac{\partial h_{n}}{\partial y}=0,$

and its solution is $h_{n}=W_{n}(z)y^{-cn}$. Since $h_{n}$ must be a
polynomial and $n\geq{1}$, $c>0$, then must be $h_{n}=0$. As a result,
$h_{n}=0$ if $n\geq{1}.$ For $n=0,$ we have $h=h_{0}(y,z),$ and the
coefficients of $x$ from (52) are

$\begin{split}y\frac{\partial h_{0}}{\partial y}&=l_{1},\\\ -y\frac{\partial
h_{0}}{\partial y}-bz\frac{\partial h_{0}}{\partial
z}&=l_{0}+l_{2}y+l_{3}z.\\\ \end{split}$ (53)

The solution of the first equation in (53), is

$h_{0}=l_{1}\ln(y)+W_{0}(z).$

The function $h_{0}$ is a polynomial if $l_{1}=0$, and putting it in the
second equation in (53), we obtain

$W_{0}(z)=(-\frac{l_{2}}{b}y-\frac{l_{0}}{b})\ln(z)-\frac{l_{3}}{b}z.$

Since $W_{0}(z)$ must be a polynomial, we get $l_{0}=l_{2}=0$ and
$W_{0}=-\frac{l_{3}}{b}z$. Then $h=-\frac{l_{3}}{b}z$ with $L=l_{3}z$.

If $b=0$, then equation (52) becomes

$x(1-y+cx)\frac{\partial h}{\partial x}+y(-1+x)\frac{\partial h}{\partial
y}=l_{0}+l_{1}x+l_{2}y,$ (54)

so by eliminating the variable $z$ and repeating the previous steps, we get

$c\,n\,h_{n}+y\frac{\partial h_{n}}{\partial y}=0,\quad\mbox{so that,}\quad
h_{n}=d_{1}y^{-cn},\;d_{1}\in{\mathbb{C}}.$

Hence, $h_{n}=0$ if $n\geq{1}$. For $n=0$, we have $h=h_{0}(y)$, and by
comparing the coefficients of $x$ from (54), we obtain

$\begin{split}y\frac{\partial h_{0}}{\partial y}&=l_{1},\\\ -y\frac{\partial
h_{0}}{\partial y}&=l_{0}+l_{2}y.\\\ \end{split}$ (55)

Solving the first equation, gives $h_{0}=l_{1}{\ln(y)}$, and it is a
polynomial if $l_{1}=0$. This implies that $h=h_{0}=0$. Eventually, the system
(1) has no exponential factors for $a=b=0$, and $c>0$.

Subcase 4.4. When $a={0}$ and $c=0$. The equation (50) takes the form

$x(1-y)\frac{\partial h}{\partial x}+y(-1+x)\frac{\partial h}{\partial
y}+z(-b)\frac{\partial h}{\partial z}=l_{0}+l_{1}x+l_{2}y+l_{3}z.$ (56)

If $b=0$, then the system has an exponential factor, as shown in the part (a)
of Theorem 1.3. So from now on, we assume that $b\neq 0$. We proceed in a
similar way to the proof of Subcase 4.1. We write $h$ in the form
$h=\sum{{}_{j=0}^{n}}h_{j}(x,y,z)$, where each $h_{j}$ is a homogeneous
polynomial of degree $j$ in its variables. Assume $n\geq 2$. As before, we
will show that $h$ has degree $n\leq 1$. The terms of degree $n+1$ in (56) are

$-xy\frac{\partial h_{n}}{\partial x}+xy\frac{\partial h_{n}}{\partial y}=0,$

and whose solution is, $h_{n}=W_{n}(x+y,z)$, where $W_{n}$ is an arbitrary
polynomial in $x+y$ and $z$. Since $h_{n}$ is homogeneous of degree $n$, we
can write $h_{n}=\sum{{}_{j=0}^{n}}a_{j}z^{n-j}(x+y)^{j}$, where
$a_{j}\in{\mathbb{C}}$. By computing the terms of degree $n$, in equation
(56), we can obtain

$h_{n-1}=-2\,b\,B\,arctanh(\frac{x-y}{x+y})-A\,\ln(-xy)+W_{n-1}(x+y,z),$

where $A=\sum{{}_{j=0}^{n}}ja_{j}(x+y)^{j-1}z^{n-j}$,
$B=\sum{{}_{j=0}^{n}}a_{j}(n-j)(x+y)^{j-1}z^{n-j}$ and $W_{n-1}$ is an
arbitrary function of $x+y$ and $z$. Since $h_{n-1}$ must be homogeneous
polynomial of degree $n-1$ we must have $A=0$ and $B=0$ which implies that
$a_{j}=0$ for $j=0,1,\dots,n$ and thus $h_{n}=0$, which is a contradiction.
Then must be $n\leq{1}$. Therefore $h=h_{0}+h_{1}x+h_{2}y+h_{3}z$.
Substituting $h$ in (56) we get that $h=h_{2}(x+y)+h_{3}z$, and
$L=h_{2}(x-y)-b\,h_{3}z$. This completes the proof of the statement b.∎

##  6\. First Integrals

The proof of Theorem 1.3 follows by direct computations, and it is omitted.∎

Proof of Theorem 1.4. According to Theorem 2.2, the system (1) has a Darboux
first integral if and only if there exist $\lambda_{i},\mu_{j}\in{\mathbb{C}}$
not all zero, such that equation (8) is satisfied. It follows from Theorem 1.1
that the system (1) have three Darboux polynomials with cofactors of
$K_{1}=1-y+cx-axz,K_{2}=-1+x$, and $K_{3}=-b+ax^{2}$. Now, by Theorem 1.2,
when $a>0$ and $c>0$, there are two exponential factors with cofactors
$L_{1}=cx^{2}-xy-bz+x$ and $L_{2}=y(x-1)$. So equation (8) is equivalent to

$\lambda_{1}(1-y+cx-
axz)+\lambda_{2}(-1+x)+\lambda_{3}(-b+ax^{2})+\mu_{1}(cx^{2}-xy-
bz+x)+\mu_{2}(y(x-1))=0.$

Solving it, we get $\lambda_{1}=\lambda_{2}=\lambda_{3}=\mu_{1}=\mu_{2}=0$. In
short there are no first integrals of Darboux type in this case.∎

Proof of Theorem 1.5. The proof of Theorem 1.5 can be obtain easily from
Theorem 1.1 and corollary 1 as well as Lemmas 2.1 and 2.3. ∎

Proof of Theorem 1.6. If $b=0$, then system (1) has the form

$\begin{split}&\dot{x}=x(1-y+cx-axz),\\\ &\dot{y}=y(-1+x),\\\
&\dot{z}=ax^{2}z,\\\ \end{split}$ (57)

Let $f$ be a formal first integral of system (57). Without loss of generality,
we can assume that $f$ has no constant terms. We write $f=\sum{{}_{j\geq
0}}f_{j}(x,y)z^{j}$, where every $f_{j}(x,y)$ is a formal power series in the
variables $x$ and $y$. We consider two cases.

Case 1: If $f$ is not divisible by $z$. In this case, we have that
$f_{0}=f_{0}(x,y)$ is a formal first integral of system (57) restricted to
$z=0$. Note that system (57), restricted to $z=0$, becomes system (21).
Consequently, $f_{0}$ is also a formal first integral of system (21). However,
we proved in Proposition 2 that system (21) has no formal first integral.
Therefore, $f_{0}$ is not a formal first integral of system (21), and we have
a contradiction.

Case 2: If $f$ is divisible by $z$. In this case, we write $f=z^{l}g$ where
$l\geq 1$, and $g$ is not divisible by $z$. Furthermore, $g$ is a formal power
series that satisfy,

$x(1-y+cx-axz)\frac{\partial g}{\partial x}+y(-1+x)\frac{\partial g}{\partial
y}+ax^{2}z\frac{\partial g}{\partial z}=a\,l\,x^{2}g.$

Now we introduce the variable $g=\exp(x)T$. Then, since $g$ is a formal series
in the variables $x,y$, and $z$, we have that $T$ is a formal power series in
the same variables, and it satisfies

$x(1-y+cx-axz)\frac{\partial T}{\partial x}+y(-1+x)\frac{\partial T}{\partial
y}+ax^{2}z\frac{\partial T}{\partial z}=(-ax^{2}z+(al+c)x^{2}-xy+x)T,$ (58)

after dividing by $\exp(x)$. We write $T$ as $T=\sum{{}_{j\geq 0}}T_{j}z^{j}$,
where every $T_{j}=T_{j}(x,y)$ is a formal power series in the variables $x$
and $y$. Furthermore, since $g$ is not divisible by $z$, we have that $T$ is
not divisible by $z$, so $T_{0}=T_{0}(x,y)\neq 0$. It follows that at least
one of $\tilde{T}_{0}=\tilde{T}_{0}(x)\neq 0$ or
$\bar{T}_{0}=\bar{T}_{0}(y)\neq 0$ must hold, where $\tilde{T}_{0}(x)$ is the
restriction of $T_{0}$ to $y=0$, and $\bar{T}_{0}(y)$ is the restriction of
$T_{0}$ to $x=0$. Otherwise, $T$ would be divisible by $z$, a contradiction.
Without loss of generality, we can assume that $\tilde{T}_{0}\neq 0$.
Moreover, if we restrict equation (58) to $y=z=0$ and simplify it, we obtain
that

$(1+cx)\frac{d\tilde{T}_{0}}{dx}=-(1+(al+c)x)\tilde{T}_{0}.$ (59)

Now we consider two cases.
Subcase 2.1 $\tilde{T}_{0}$ is not divisible by $(1+cx)$. In this case, since
$\tilde{T}_{0}\neq 0$ and $al\neq 0$ from (59), we get a contradiction.
Subcase 2.2 $\tilde{T}_{0}$ is divisible by $(1+cx)$. We write
$\tilde{T}_{0}(x)=(1+cx)^{m}h(x)$ with $m\geq 1$ and $h(x)$ being a formal
power series in the variable $x$ that is not divisible by $(1+cx)$. Then from
(59) we obtain $h$ must satisfy, after dividing by $(1+cx)^{m}$,

$(1+cx)\frac{dh}{dx}=-(1+(al+c)x-cm)h.$ (60)

Since $cm>0$ and $al>0$, we have from (60) that $h$ must be divisible by
$(1+cx)$, a contradiction. So system (1) has no formal first integral if
$b=0$. ∎

Now we consider system (1) with $a\neq 0$,$c\neq 0,$ and $b\neq 0$. We provide
the proof of Theorem 1.7 in order to prove Theorem 1.6 when $b>0$.

Proof of Theorem 1.7. Since $b$ is one of the parameters in the system, we can
rewrite it as

$\begin{split}&\dot{x}=x(1-y+cx-axz),\\\ &\dot{y}=y(-1+x),\\\
&\dot{z}=z(-b+ax^{2}),\\\ &\dot{b}=0.\\\ \end{split}$ (61)

In other words, we add a new variable, $b$, that was a parameter in system
(1). Note that a non constant function $f=f(b)$ is a first integral of system
(61), but it is not a first integral of the system (1).
We assume that $f=f(x,y,z,b)$ is a formal power series first integral of
system (61). Expanding $f$ in powers of the variable $b$, we get
$f=\sum{{}_{k\geq 0}}f_{k}(x,y,z)b^{k}$, where each $f_{k}$ is a formal series
in its variables. We can write $f=f_{0}+bg$ where $g=\sum{{}_{k\geq
0}}f_{k+1}b^{k}=\sum{{}_{k\geq 0}}f_{k+1}(x,y,z)b^{k}$ is a formal series in
variables $x,y$, and $z$. Since $f(x,y,z,0)$ is a formal first integral of the
system (1) with $b=0$, and since $a\neq 0$, and $c\neq 0$, we are in the
assumptions of Theorem 1.6 for $b>0$, applying it, we get that
$f(x,y,z,0)=f_{0}=d_{0}$. We claim that

$f_{k+1}=d_{k+1},\quad k\geq 0,$ (62)

where $d_{k+1}$ are constants for $k\geq 0$. Then, since $f$ is a first
integral, it satisfies ${\mathcal{X}}f=0$. So, the function $g$ satisfies the
equation

$x(1-y+cx-axz)\frac{\partial g}{\partial x}+y(-1+x)\frac{\partial g}{\partial
y}+(b+ax^{2}z)\frac{\partial g}{\partial z}=0.$ (63)

Moreover, $f_{1}=g(x,y,z,0)$ satisfies (63) restricted to $b=0$, that is,

$x(1-y+cx-axz)\frac{\partial f_{1}}{\partial x}+y(-1+x)\frac{\partial
f_{1}}{\partial y}+(ax^{2}z)\frac{\partial f_{1}}{\partial z}=0.$

Hence, $f_{1}$ is a formal power series first integral of system (1) with
$b=0$, and by assumption, additionally, we have that $a\neq 0$, and $c\neq 0$.
So, from Theorem 1.6 for $b>0$, we obtain that $f_{1}=d_{1}$, a constant. This
proves (62) for $k=0$. Now we assume (62) is true for $k=0,\dots,j-1$ with
$j\geq 1$, and we will prove it for $k=j$. Since $g=bd_{1}+b\sum{{}_{k\geq
1}}f_{k+1}b^{k}$, then $f=d_{0}+bd_{1}+b\sum{{}_{k\geq 1}}f_{k+1}b^{k}$, and
by the induction hypothesis, we have

$f=\sum{{}_{k=0}^{j-1}}d_{k}b^{k}+b^{j+1}\sum{{}_{k\geq j}}f_{k+1}b^{k-j},$

so $\sum{{}_{k\geq j}}f_{k+1}b^{k-j}$ is a formal first integral of system
(4). Consequently, $f_{j+1}$ is a formal first integral of system (1) with
$b=0$, and from part (a), with $a\neq 0$,$c\neq 0$, we get $f_{j+1}=d_{j+1}$,
and this proves the claim for $k=j$. Then, from (62), we get that
$f=\sum{{}_{k\geq 0}}d_{k+1}b^{k}$, which finishes the proof of the theorem.∎

## References

  * [1] A. Arneodo, P. Coullet, C. Tresser, Occurence of strange attractors in three-dimensional Volterra equations, Physics Letters A. 79 259–263 (1980).
  * [2] A. Seidenberg, Reduction of singularities of the differential equation Ady= Bdx, American Journal of Mathematics. 90 248–269 (1968).
  * [3] AA. Elsadany, AE. Matouk, AG. Abdelwahab, HS. Abdallah, Dynamical analysis, linear feedback control and synchronization of a generalized Lotka-Volterra system, International Journal of Dynamics and Control. 6 328–338 (2018).
  * [4] AI. Amen, Integrability Analysis of the Smallest 3D Biochemical Reaction Model, MATCH Commun. Math. Comput. Chem. 90 333–356 (2023).
  * [5] C. Christopher, J. Llibre, Integrability via invariant algebraic curves for planar polynomial differential systems, Annals of Differential Equations. 16 5–19 (2000).
  * [6] C. Christopher, J. Llibre, JV. Pereira, Multiplicity of invariant algebraic curves and Darboux integrability, Pacific J. Math. 229 63–117 (2007).
  * [7] G. Darboux, Mémoire sur les équations différentielles algébriques du premier ordre et du premier degré, Bulletin des Sciences Mathématiques et Astronomiques. 2 151–200 (1878).
  * [8] J. Écalle, Introduction aux fonctions analysables et preuve constructive de la conjecture de Dulac, Hermann. (1992).
  * [9] J. Llibre, X. Zhang, Darboux theory of integrability for polynomial vector fields in Rn taking into account the multiplicity at infinity, Bulletin des sciences mathematiques. 133 765–778 (2009).
  * [10] J. Llibre, X. Zhang, Darboux theory of integrability in Cn taking into account the multiplicity, Journal of Differential Equations. 246 541–551 (2009).
  * [11] J. Llibre, X. Zhang, On the Darboux integrability of polynomial differential systems, Qual. Theory Dyn. Syst. 11 129–144 (2012).
  * [12] JF. Mattei, R. Moussu, Holonomie et intégrales premieres, Annales scientifiques de l’École Normale Supérieure. 13 469–523 (1980).
  * [13] JS. Costello, Synchronization of chaos in a generalized Lotka–Volterra attractor, Nonlinear J. 1 11–17 (1999).
  * [14] KH. Bachmann, Henrici, P., Applied and Computational Complex Analysis, Bd. I, 682 S., New York-London-Sydney-Toronto. John Wiley & sons. 1974.£ 13, 50., Zeitschrift Angewandte Mathematik und Mechanik. 57 352–352 (1977).
  * [15] N. Samardzija, LD. Greller, Explosive route to chaos through a fractal torus in a generalized Lotka-Volterra model, Bulletin of Mathematical Biology. 50 465–491 (1988).
  * [16] R. Moussu, Symétrie et forme normale des centres et foyers dégénérés, Ergodic Theory and Dynamical Systems. 2 241–251 (1982).
  * [17] VV. Nemytskii, Qualitative theory of differential equations, Princeton University Press. 2083 (2015).
  * [18] W. Aziz, A. Amen, C. Pantazi, Integrability and linearizability of a family of three-dimensional quadratic systems, Dynamical Systems. 36 317–331 (2021).
  * [19] X. Zhang, Integrability of dynamical systems: algebra and analysis, Springer. 47 (2017).