\section{Introduction}
Quantum computers challenge Moore's law, and their computational power is beyond that of classical computers. Developing powerful features and algorithms for quantum computers are thriving; notable research topics include teleportation \cite{Stefano15}, computation \cite{Cerezo21}, information theory \cite{Horodecki21}, and resource theory \cite{Chitambar19}, and so forth. Within the mentioned topics, entanglement contributes substantially owing to its nonnegligible ability to carry massive data, thereby speeding up algorithms and communication in computation and teleportation. However, entanglement is not always necessary~\cite{Knill98,biham04,Lanyon08}; separability can also increase the capacity in computation and information processing.

Entanglement is a physical phenomenon in which particles in a quantum system come together and share the same information simultaneously. Moreover, when entangled, the states of each particle are closely correlated and cannot be described independently of the state of other particles, which can be embedded in "quantum correlation." Take, for example, Einstein-Podolsky-Rosen (EPR) pairs which are the key components of communication\cite{Gisin07,Mattle96} are well-known entangled states in a bipartite system.

A separable state can fully portray information in composite systems. Separable states are considered to be states of classical correlation because they can be prepared through local operations and classical communication (LOCC), and local measurement operators can precisely understand the information in each joint subsystem. Even though separable states exhibit "classical correlation," some can occasionally accommodate "quantumness."

Detection, certification, and quantification of quantumness are related research topics. Since quantum entanglement possesses a high degree of quantumness, many research topics are desired to determine whether a state is entangled. There are some known criteria \cite{Peres96,Chen03, Rudolph03, Albeverio07} in a bipartite system but still elusive in a multipartite system. Furthermore, quantifying the degree of entanglement through relative entropy\cite{Vedral97} is another approach to understanding.
Besides, some separable states possess quantumness. Take quantum discord for example; in 2001, Ollivier, Harold, and Zurek first developed a concept of a separable state in a system that exhibits quantum correlation \cite{Ollivier01}. In addition, Henderson and Vedral conducted entropy-based quantification \cite{Henderson01}.
Identifying whether a state exhibits quantumness or classicality becomes worth discussing. In \cite{Groisman07}, Groisman, Kenigsberg, and Mor defined classical states and provided a method to gauge the degree of "quantumness" (or "non-classicality") of a given state.
However, this is not an easy task. Determining the separability of a given state is a primary concern before understanding the degree of classicality of a state. Even in a bipartite system, the separability-determined problem was proven NP-hard by Gurvits \cite{Gurvits04, Gharibian10}, showing this problem cannot be tackled by a nondeterministic Turing machine in a polynomial time.

Quantifying nonclassicality or quantumness is difficult, but perhaps possible to approximate it. Viewing the problem as an optimization problem, Chu and Lin\cite{Moody2022} first proposed a projected gradient method to approximate a given state by a linear combination of unit vectors over complex variables.
However, in quantifying classicality, classical-classical states further require orthogonal matrices in each subsystem, which Stiefel manifolds can characterize and facilitates orthonormal-matrix-based optimization problems.

The paper is organized as follows, Section 2 reviews mathematical notations and the foundations of quantum mechanics. It also outlines the structure and properties of the Stiefel manifold. The main problem is addressed using mathematical formation; details are provided at the end of the section. In Section 3, we describe the steps in solving the optimization problem which involves computing the function gradient on the manifold and establishing a descent flow. We ensure the descent flow is convergent and the trajectory stays on Stiefel manifolds through rigorous proofs. In Section 4, numerical results are provided to demonstrate the ability to decompose a quantum state into orthogonal matrices and the stability of our method. Moreover, through numerical results, the proposed method preserves the fundamental properties of quantum states. Finally in this paper, we suggest an optimal selection of the initial guesses to the descent flow.

\section{Preliminaries}
\subsection{Basics of quantum machanics}
In this section, we introduce some fundamental concepts of quantum computation with mathematical notation and address the motivation of our study. Readers can refer to \cite{Nielsen00} for further details. In quantum computation, qubits are relied upon for storing and manipulating information, which is analogous to 'bits' in classical computation. Appropriate mathematical notation can simplify discussions and improve reader understanding. A quantum system is considered as a Hilbert space $\mathcal{H}$ and a state of a qubit is represented as a unit vector $\mathbf{x}\in\mathcal{H}$. Real or complex spaces can be discussed for $\mathcal{H}$, but we consider only real spaces in the present study.\\

\subsubsection{Separable State versus Entangle State}
In a composite system, states can be classified as separable or entangled. Take, a bipartite system $\mathcal{H}_{AB} = \mathcal{H}_{A}\otimes \mathcal{H}_{B}$ for example. A separable state follows the representation,
\begin{equation}\label{eq:sep}
    \rho = \sum_{i}p_i \rho_A^i\otimes \rho_B^i,
\end{equation}
where $\rho_A^i$ and $\rho_B^i$ are density matrices in $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$ respectively, and $p_i$ is some probability distribution. In other words, an entangled state cannot be written in the form \eqref{eq:sep}. The multipartite system follows analogously.
\subsubsection{Classical States}
Some separable states exhibit quantumness. Quantifying the degree of quantumness or non-classicality is a critical also important task in quantum information theory.
In an isolated system $\mathcal{H}_{A}$, a computational basis is said to be classical if it is orthonormal \cite{Groisman07}. Moreover, in a bipartite system, $\mathcal{H}_{AB} = \mathcal{H}_{A}\otimes \mathcal{H}_{B}$, a state $\rho_{AB}$ is said to be classical-classical state if $\rho_{AB} = \displaystyle\sum_{i,j} p_{ij}\mathbf{u}_i\mathbf{u}_i^T\otimes \mathbf{v}_j\mathbf{v}_j^T$,
where $\left\{\mathbf{u}_i\right\}$ and $\left\{\mathbf{v}_i\right\}$ are orthonormal basis in $\mathcal{H}_A$ and $\mathcal{H}_B$ respectively.
\subsection{Distancing}
For the sake of quantifying the "quantumness" or "non-classicality" of a given state, the similarity between the two states must be defined. Different inquiries require different quantifying approaches; we focus on the Frobenius norm to gauge the similarity between states. Suppose $\rho$ and $\sigma$ are two given states in a quantum system. The
Frobenius norm is given as follows:
\begin{equation*}
    D_{F}(\rho, \sigma) = \norm{\rho-\sigma}_F = \sqrt{\sum_{i,j} (\rho_{i,j}-\sigma_{i,j})^2}.
\end{equation*}
Therefore, the "quantumness" or "non-classicality" of $\rho$ can be determined as:
\begin{equation*}
    \mathcal{Q}(\rho) = \min_{\sigma\in \mathcal{D}}D_{F}(\rho, \sigma),
\end{equation*}
where $\mathcal{D}$ is a collection of classical-classical states. A state is said to be classical-classical if $\mathcal{Q}(\rho) = 0$; therefore such a quantum state inherits no "quantumness" from a quantum system.

\subsection{Stiefel Manifold}
The classical-classical formulation requires orthogonal matrices, and the related minimization problem can be discussed on Stiefel manifolds. Stiefel manifold characterizes the geometric details of orthonormal matrices and benefits solving those optimization problems with orthogonal constraints \cite{Edelman1998}. A proper coordinate system on a manifold facilitates designing and applying numerical algorithms by framing abstract objects as concrete objects. When a coordinate system finely characterizes the structure of the Stiefel manifold, many related optimization problems can be resolved by geometric quantities (e.g., gradient, curvature).\\
\subsubsection{Geometric Properties}
The Stiefel manifold $\mathcal{S}_{m,n} = \left\{Y\in\mathbb{R}^{m \times }\rvert Y^T Y = \mathbf{I}_{n\times n}, m\geq n\right\}$ collects the orthonormal matrices in $\mathbb{R}^{m \times n}$. Embedding $\mathcal{S}_{m,n}$ into $mn$-dimensional Euclidean space $\mathbb{R}^{mn}$ is a useful strategy to study it's structure. Furthermore, the tangent and normal space at point $Y$ are given as follows:
\begin{itemize}
    \item (Tangent Space) $\mathcal{T}_{Y}\mathcal{S}_{m,n} = \left\{X\in \mathbb{R}^{m\times n}\rvert X^TY+Y^TX = 0\right\}$,\\
    \item (Normal Space) $\mathcal{N}_{Y}\mathcal{S}_{m,n} = \left\{YS\lvert S\mbox{ is any } n\mbox{-by-}n \mbox{ symmetric matrix} \right\}$.
\end{itemize}

\subsubsection{Riemannian Gradient}
The gradient of a function on a Riemannian manifold must be derived for some optimization methods, such as the conjugate gradient method, the gradient flow method, and so on. The computation of the Riemannian gradient of a function depends on the choice of metrics. Readers may refer to \cite{Lee19} for further details. Herein, the Riemannian gradient of a real-valued function $F$ on the tangent space of $\mathcal{S}_{m,n}$ at point $Y$ is a tangent vector $\mbox{Grad}F\rvert_Y\in \mathcal{T}_{Y}\mathcal{S}_{m,n}$ with the canonical metric $g_c(X_1, X_2) = \inner{X_1, (I-\frac{1}{2}YY^T)X_2}$, where $X_1, X_2\in \mathcal{T}_{Y}\mathcal{S}_{m,n}$, is given as
\begin{equation}\label{eq:gradient}
    \mbox{Grad}F\rvert_Y = 2\mbox{skew} \left(\left(\nabla  F\left|\right._Y\right)Y^T\right)Y,
\end{equation}
where $\nabla F\big\rvert_Y$ is the Fr\'echet derivative of $F$ at $Y$, and $\mbox{skew}(A) = \frac{1}{2}\left(A-A^T\right)$.
\section{Problem Description}
    Quantifying the quantumness in a bipartite system is equivalent to determining whether a given state can be decomposed into orthogonal matrices. The idea motivated us to find the nearest classical-classical state to a given state in a real finite dimensional bipartite system, $\mathcal{H}_A\otimes\mathcal{H}_B$.
    To simplify the discussion, we consider the classical-classical state without the interaction of vectors. In other words, we only discuss the classical-classical state with the following representation,
    \begin{eqnarray}\label{eq:sigma}
        \sigma &=& \sum_{i=1}^N\theta_i \left(\mathbf{x_i}\mathbf{x_i}^T\right)\otimes\left(\mathbf{\mathbf{y_i}}\mathbf{\mathbf{y_i}}\right)^T,\\
        &=& \sum_{i=1}^N\theta_i  \left(\mathbf{\mathbf{x_i}}\otimes \mathbf{\mathbf{y_i}}\right) \left(\mathbf{x_i}\otimes \mathbf{\mathbf{y_i}}\right)^T.\\
    \end{eqnarray}
    where $\theta_i$ are nonnegative real numbers and with unit sum; $\left\{\mathbf{x_i}\right\}\subset \mathcal{H}_A=\mathbb{R}^n$, and  $\left\{\mathbf{y_i}\right\}\subset \mathcal{H}_B=\mathbb{R}^m$ are orthonormal sets.
    The general classical-classical state can also be discussed by extending the $M$ matrix and $N$ matrix in \lemmaref{lemma:MN_rep} to larger matrices that store $M$ and $N$ matrix in the diagonal respectively. Therefore our problem can be rephrased as minimizing the following norm,
    \begin{equation}\label{eq:vector_form}
        \norm{\rho-\sum_{i=1}^N\theta_i \left(\mathbf{x_i}\otimes \mathbf{y_i}\right) \left(\mathbf{x_i}\otimes \mathbf{y_i}\right)^T}_F^2.
    \end{equation}
    Using the Kahati-Rao product of $U \mbox{ and } V$ to simplify the classical-classical formulation:
    \begin{eqnarray}
        \left(U\odot V\right)_i = \left(x_i\otimes y_i\right),
    \end{eqnarray}
    where the column of the Kahati-Rao product stores the tensor product of $x_i$ and $y_i$.
    The degree of quantumness is to minimize the constrained optimization problem,
    \begin{eqnarray}\label{eq:objective_function}
        &&F(U, V, \Sigma) = \min_{\Sigma ,U\in\mathbb{R}^{n\times N}, V\in\mathbb{R}^{m\times N}} \dfrac{1}{2}\norm{\rho-(U\odot V)\Sigma (U\odot V)^T}_F^2\\
        && \mbox{ s.t }  \sum_{i=1}^N \theta_i = 1, \mbox{ and } \theta_i\in(0,1],\\
        && U^T U = \mathbf{I}_{N\times N}, V^T V = \mathbf{I}_{N\times N},
    \end{eqnarray}
    where $U = [\mathbf{x_1}, \mathbf{x_2}, \cdots \mathbf{x_N}]$ and $V = [\mathbf{y_1}, \mathbf{y_2}, \cdots \mathbf{y_N}]$, and $\mathbf{I}_{N\times N}$ is the $N$-by-$N$ identity matrix.

\section{Optimization}
The gradient descent method is widely considered to solve optimization problems. Herein, we establish
 a dynamic system that is inspired by the gradient descent method, to solve the aforementioned optimization problem\eqref{eq:objective_function}.
    First, we compute the gradient of the objective function on variable $U$ and $V$, $\nabla_{U, V} F = \left(\dfrac{\partial F}{\partial U}, \dfrac{\partial F}{\partial V}\right)$, and project it onto Stiefel manifolds to establish descent flows. Subsequently, the descent flows to  $\theta_i$ are obtained by taking the function gradient of $\theta_i$ to satisfy the constraint. Finally, the minimizer of \eqref{eq:objective_function} is obtained by solving the dynamic system.\\
    \subsection{Gradients of Objective Function}
    The gradient of a function on a Stiefel manifold is formulated by computing each variable's derivative and then transforming them to the Riemannian gradient on the manifold. \\
    The following lemma presents an explicit representation of the partial derivatives of each variable.
    \begin{lemma}\label{lemma:MN_rep}
    Let $M$, and $N$ satisfy the following representation,
        \begin{equation*}
        \left\{
            \begin{array}{cc}
                \vectorize(X\odot V) = M \vectorize(X),\\
                \\
                \vectorize(U\odot Y) = N \vectorize(Y).\\
            \end{array}
        \right.
    \end{equation*}
        The gradient $\nabla F$ of the objective function \eqref{eq:objective_function} is given as,
        \begin{equation}
            \left\{
                \begin{array}{ccl}
                    \dfrac{\partial F}{\partial U} & = &  -2\left(\tilde{M}-\hat{M}\right),\\
                    \\
                    \dfrac{\partial F}{\partial V} & = &  -2\left(\tilde{N}-\hat{N}\right),\\
                    \\
                    \dfrac{\partial F}{\partial \theta_i} & = &\theta_i-\inner{(\mathbf{x}_i\otimes \mathbf{y}_i)(\mathbf{x}_i \otimes \mathbf{y}_i)^T,\rho},
                \end{array}
                \hspace{10pt} i\in \left\{1\cdots N\right\},
            \right.
        \end{equation}
     where
    \begin{equation*}
        \left\{
            \begin{array}{ccc}
                \tilde{M} & = & \reshape \left(M^T \vectorize(\rho(U\odot V)\Sigma), [n, N]\right),\\
                \\
                \hat{M} & = & \reshape \left(M^T \vectorize\left((U\odot V)\Sigma^2\right), [n, N]\right),\\
                \\
                \tilde{N} & = & \reshape \left(N^T \vectorize\left(\rho (U\odot V\right)\Sigma), [m, N]\right),\\
                \\
                \hat{N} & = & \reshape \left(N^T \vectorize\left((U\odot V)\Sigma^2\right), [m, N]\right),\\
            \end{array}
        \right.
    \end{equation*}

    \end{lemma}
    Note that "$\vectorize$" converts a matrix into a column vector by stacking the columns of a given matrix on top of one another, and "
    $\reshape$" transforms a matrix into another matrix of a given shape in which elements are taken and restored column-wisely from the original matrix. \\
    \begin{proof}

        First, consider the Fr\'echet derivative of \eqref{eq:objective_function} to $U$ and  $ V$. \\
        Let
        \begin{eqnarray*}
            F(U,V, \Sigma) = \dfrac{1}{2}\norm{\rho-(U\odot V)\Sigma (U\odot V)^T}_F^2 = \dfrac{1}{2}\norm{\beta(U,V, \Sigma)}_F^2,
        \end{eqnarray*}
        where $\beta (U,V, \Sigma) = \rho-(U\odot V)\Sigma(U\odot V)^T$.\\

        By the chain rule and the product rule, the partial Fr\'echet derivatives of $F$ at $U$ and $V$ acting on matrices $X\in\mathbb{R}^{n\times N}$ and $Y\in\mathbb{R}^{m\times N}$ respectively are given as follows,
        \begin{equation*}
            \left\{
                \begin{array}{ccl}
                    \dfrac{\partial F}{\partial U}\bullet(X)

                    & = &  \left<\beta(U,V), -(U\odot V)\Sigma(X\odot V)^T-(X\odot V)\Sigma(U\odot V)^T \right>,
                    \\
                    \\
                    \dfrac{\partial F}{\partial V}\bullet(Y)

                    & = &  \left<\beta(U,V), -(U\odot V)\Sigma(U\odot Y)^T-(U\odot Y)\Sigma(U\odot V)^T\right>,
                \end{array}
            \right.
        \end{equation*}

        Through direct calculation, we have
        \begin{eqnarray*}
            \dfrac{\partial F}{\partial U}\bullet(X) &=& -2\left<\rho,(X\odot V)\Sigma(U\odot V)^T\right>\\
                &+& 2\left< (X\odot V)\Sigma(U\odot V)^T, (U\odot V)\Sigma(U\odot V)^T\right>,\\
                &=& \left<-2(\tilde{M}-\hat{M}), X\right>,
        \end{eqnarray*}
        and
        \begin{eqnarray*}
            \dfrac{\partial F}{\partial V}\bullet(X) &=& -2\left<\rho, (U\odot Y)\Sigma(U\odot V)^T\right>\\
                &+& 2\left<(U\odot Y)\Sigma(U\odot V)^T, (U\odot V)\Sigma(U\odot V)^T\right>,\\
                &=& \left<-2(\tilde{N}-\hat{N}), Y\right>,
        \end{eqnarray*}
        Regarding the partial derivative $\dfrac{\partial F}{\partial \Sigma}$, because $\Sigma$ is a diagonal matrix, we only consider $\dfrac{\partial F}{\partial \theta_i}$. \\
        Recalling \eqref{eq:objective_function}, the partial derivative of $F$ with respect to $\theta_i$ is given as,
        \begin{equation*}
            \dfrac{\partial F}{\partial \theta_i} = -\inner{(\mathbf{x_i}\otimes \mathbf{y_i})(\mathbf{x_i} \otimes \mathbf{y_i})^T,\rho-\sum_{j=1}^N\theta_j \left( \left(\mathbf{x_j}\otimes \mathbf{y_j}\right) \left(\mathbf{x_j}\otimes \mathbf{y_j}\right)^T\right)}.
        \end{equation*}
        Because $(\mathbf{x_i}\otimes \mathbf{y_i})$ are mutually orthogonal, $\dfrac{\partial F}{\partial \theta_i}$ becomes
        \begin{equation*}
            \dfrac{\partial F}{\partial \theta_i} = \theta_i
            -\inner{(\mathbf{x_i}\otimes \mathbf{y_i})(\mathbf{x_i} \otimes \mathbf{y_i})^T,\rho},
        \end{equation*}
        which completes the derivation.
    \end{proof}
    As \eqref{eq:objective_function} is restricted on the Stiefel manifolds where $U\in\mathcal{S}_{n,N}$ and $ V\in\mathcal{S}_{m,N}$. Therefore, we compute the Riemannian gradients $\left(\mbox{Grad}_U F, \mbox{Grad}_V F\right)$ of our objective function. We project the partial derivative of $U$ and $V$ to a matrix on the tangent space of each Stiefel manifold to specify the relevant Riemannian gradient respectively. The following lemma provides further details.
        \begin{lemma}
            The Riemannian gradients $ \left(\mbox{Grad}_U F, \mbox{Grad}_V F\right) = \left(\mbox{Proj}_{\mathcal{S}_{n, N}}\dfrac{\partial F}{\partial U}, \mbox{Proj}_{\mathcal{S}_{m, N}}\dfrac{\partial F}{\partial V}\right)$ on $\mathcal{T}_U\mathcal{S}_{n,N}$ and $\mathcal{T}_V\mathcal{S}_{m,N}$ respectively follow the expressions,

            \begin{equation}\label{eq:proj}
                \left\{
                    \begin{array}{ccl}
                        \mbox{Proj}_{\mathcal{S}_{n, N}}\dfrac{\partial F}{\partial U} & = & 2\mbox{ skew }\left(\left(\dfrac{\partial F}{\partial U}\right)U^T\right)U,\\
                        \\
                        \mbox{Proj}_{\mathcal{S}_{m, N}}\dfrac{\partial F}{\partial V} & = & 2\mbox{ skew }\left(\left(\dfrac{\partial F}{\partial V}\right)V^T\right)V.\\
                    \end{array}
                \right.
            \end{equation}
        \end{lemma}

    Taking the negative Riemannian gradient in \eqref{eq:proj}, we define the dynamic system,
    \begin{proposition}
    Consider the dynamic system
        \begin{equation}\label{eq:dy_sys}
            \left\{
                \begin{array}{ccl}
                    \dfrac{dU}{dt} = -2\mbox{ skew }\left(\left(\dfrac{\partial F}{\partial U}\right)U^T\right)U,\\
                    \\
                    \dfrac{dV}{dt} = -2\mbox{ skew }\left(\left(\dfrac{\partial F}{\partial V}\right)V^T\right)V,\\
                \end{array}
            \right.
        \end{equation}
        where $t$ represents a dimensionless parameter at a given time.\\
        Let $U(0)\in\mathcal{S}_{n, N}$ and $V(0)\in\mathcal{S}_{m, N}$ be the initial data of the differential system \eqref{eq:dy_sys}, Then the trajectories of $U(t), V(t)$ stay on the manifold $\mathcal{S}_{n, N}$ and $\mathcal{S}_{m, N}$ respectively.
    \end{proposition}
    \begin{proof}
        We only prove for $U$, since $V$ can be discussed similarly.\\
        Consider the equation
        \begin{equation*}
            g(t) = U^T(t)U(t),
        \end{equation*}
        we want to prove $g(t) = g(0) = \mathbf{I}_N$ to ensure the trajectories of $U(t)$ will stay on the manifold.
        Differentiate both sides of the equation,
        \begin{equation}\label{eq:prop_U}
            g'(t) = \dot{U}^T(t)U(t) + U^T \dot{U}(t),
        \end{equation}
        where $\dot{U}(t) =\dfrac{d}{dt}U(t)$.
        Thus substituting $\dot{U}(t) = -2\mbox{ skew }\left(\left(\dfrac{\partial F}{\partial U}\right)U^T\right)U$ to \eqref{eq:prop_U}, we have $g'(t) = 0$ which shows $g(t) = g(0)$. Because the initial data starts on the Stiefel manifold, in which $g(0) = U^T(0) U(0) =\mathbf{I}_N$ and together with $g'(t)=0 $, we conclude $g(t) = g(0) = \mathbf{I}_N$.\\
    \end{proof}
    Third, as a dynamic system on Stiefel manifolds is defined by the proposition above,
    \begin{equation}\label{eq:gradient_flow}
        \left\{
            \begin{array}{ccl}
                \dfrac{dU}{dt} & = &  -4\left(\mbox{ skew }\left(\left(-\tilde{M}+\hat{M}\right)U^T\right)U\right),\\
                \\
                \dfrac{dV}{dt} & = & -4\left(\mbox{ skew }\left(\left(-\tilde{N}+\hat{N}\right)U^T\right)V\right),\\

            \end{array}
        \right.
    \end{equation}

    The stationary points $U$ and $V$ of the flow \eqref{eq:gradient_flow}, which are the minimizers of the objective function are discussed by the following lemma.
    \begin{lemma}
        Let $U\mbox{ and } V$ be the stationary points. Then
        \begin{equation}
            \left\{
                \begin{array}{cll}
                    U^T(\tilde{M}-\hat{M})-(\tilde{M}-\hat{M})^TU &=& 0,\\
                    \\
                    V^T(\tilde{N}-\hat{N})-(\tilde{N}-\hat{N})^TV &=& 0 .
                \end{array}
            \right.
        \end{equation}
    \end{lemma}
    \begin{proof}
        We first consider the flow of $U$, and the result of $V$ can be derived analogously.
        The stationary point of $U$ satisfies
        \begin{eqnarray}\label{eq:nothing}
            0 = \dfrac{dU}{dt} =  -4\mbox{ skew }\left(\left(-\tilde{M}+\hat{M}\right)U^T\right)U,
        \end{eqnarray}
        for a fixed $V$.\\
        Next, by direct calculation,
        \begin{eqnarray*}
            0 &=& \left(-\tilde{M}+\hat{M}\right)U^TU-U\left(-\tilde{M}+\hat{M}\right)^TU,\\
            &=& \left(-\tilde{M}+\hat{M}\right)-U\left(-\tilde{M}+\hat{M}\right)^TU\\
            &=& U^T\left(-\tilde{M}+\hat{M}\right)-\left(-\tilde{M}+\hat{M}\right)^TU
        \end{eqnarray*}
        and thus completes the proof.
    \end{proof}
    \begin{lemma}[Theorem 2.2 \cite{Absil2005}]\label{lemma:conv_lemma}
        Let $\phi$ be a real analytic function.
        Assume that there exists a $\delta >0$ and a real $\tau$ such that for $t> \tau$, $x(t)$ satisfies the angle condition
        \begin{equation*}
            \dfrac{d\phi (x(t))}{dt} = \left<\nabla \phi(x(t)), \dot{x}(t)\right> \leq -\delta \norm{\nabla \phi(x(t))}\norm{\dot{x}(t)}
        \end{equation*}
        and a weak decrease condition
        \begin{equation*}
           \dfrac{d}{dt}\phi(x(t)) = 0\Rightarrow \dot{x}(t) = 0.
        \end{equation*}
        Then, either $\lim_{t\rightarrow\infty}\norm{x(t)} = \infty$, or there exists $x^*$ such that $\lim_{t\rightarrow\infty} x(t) = x^*$.
    \end{lemma}
    Note that an inner product on the product space $\mathbb{R}^{n\times N}\times \mathbb{R}^{m\times N}$ can be defined through the Frobenius inner product via the relationship,
    \begin{equation*}
            \inner{(a_1,b_1), (a_2, b_2)} = \inner{a_1, a_2}+\inner{b_1,b_2}.
    \end{equation*}
    where $a_i\in \mathbb{R}^{n\times N}$ and $b_j\in \mathbb{R}^{m\times N}$ for $i,j = 1,2$.
    Based on Lemma \ref{lemma:conv_lemma}, we now show the existence of the stationary points $U^*$ and $V^*$ of the dynamics system \eqref{eq:gradient_flow}.
    \begin{theorem}\label{thm:stationary}
        The dynamic system \eqref{eq:gradient_flow} guarantees the existence of stationary points $U^*$ and $V^*$ such that $\lim_{t\rightarrow\infty} U(t) = U^*$ and $\lim_{t\rightarrow\infty} V(t) = V^*$.
    \end{theorem}
    \begin{proof}
        We set $\left(\dfrac{\partial F}{\partial U}, \dfrac{\partial F}{\partial V}\right) = (f_1, f_2)$.\\

        Observe that $\inner{(f_1, f_2), (\dot{U}(t),\dot{V}(t))} = \inner{f_1, \dot{U}(t)}+\inner{f_2, \dot{V}(t)}$.\\
        For $\inner{f_1, \dot{U}(t)}$,
        \begin{eqnarray*}
            \inner{f_1, \dot{U}(t)}  &=& -2\inner{f_1, \mbox{ skew}\left( f_1U^T\right)U}, \\
            &=& -\left(\inner{f_1, f_1} -\inner{f_1, U f_1^T U }\right).
        \end{eqnarray*}
        And we consider
        \begin{eqnarray*}
            \norm{\dot{U}(t)}_F^2 &=&\inner{f_1, f_1}-2\inner{U f_1^T U, f_1}+\inner{f_1^T U, f_1^T U}\\
            & = & \inner{f_1, f_1}-\inner{U f_1^T U, f_1}+\inner{f_1^T U, f_1^T U}-\inner{U f_1^T U, f_1},\\
            & = & \inner{f_1, f_1}-\inner{U f_1^T U, f_1}+\inner{f_1^T U, f_1^T U}-\inner{ f_1^T U, U^T f_1},\\
            & = & -\inner{f_1, \dot{U}(t)}+\norm{f_1^T U}_F^2-\inner{ f_1^T U, U^T f_1}.
        \end{eqnarray*}
        By Cauchy-Schwartz inequality,
        \begin{eqnarray*}
            \left| \inner{ f_1^T U, U^T f_1}\right|\leq \norm{f_1^T U}_F \norm{U^T f_1}_F = \norm{f_1^T U}_F^2.
        \end{eqnarray*}
        We then have $\norm{\dot{U}(t)}_F^2 \leq -\inner{f_1, \dot{U}(t)}$.\\
        Similarly,
        \begin{equation*}
            \norm{\dot{V}(t)}_F^2\leq -\inner{f_2, \dot{V}(t)}.
        \end{equation*}
        Thus
        \begin{eqnarray*}
            \inner{\left(\dfrac{\partial F}{\partial U}, \dfrac{\partial F}{\partial V}\right), \left(\dot{U}(t), \dot{V}(t)\right)} \leq -\left(\norm{\dot{U}(t)}_F^2+\norm{\dot{V}(t)}_F^2\right),
        \end{eqnarray*}
        By definition,
        \begin{eqnarray*}
            \norm{\left(\dfrac{\partial F}{\partial U}, \dfrac{\partial F}{\partial V}\right)}_F &=& \sqrt{\norm{f_1}_F^2+\norm{f_2}_F^2},\\
            \norm{\left(\dot{U}(t), \dot{V}(t)\right)}_F &=& \sqrt{\norm{\dot{U}(t)}_F^2+\norm{\dot{V}(t)}_F^2}.\\
        \end{eqnarray*}
        By arithmetic-geometric mean inequality,
        \begin{eqnarray*}
            \norm{\left(\dfrac{\partial F}{\partial U}, \dfrac{\partial F}{\partial V}\right)}_F\norm{\left(\dot{U}(t),
            \dot{V}(t)\right)}_F&\leq&\dfrac{1}{2}\left(\norm{f_1}_F^2+\norm{f_2}_F^2+\norm{\dot{U}(t)}_F^2+\norm{\dot{V}(t)}_F^2\right).
        \end{eqnarray*}
        Since $\dot{U}(t)$ and $\dot{V}(t)$ are on the tangent spaces $\mathcal{T}_U\mathcal{S}_{n,N}\mbox{, and } \mathcal{T}_V\mathcal{S}_{m,N}$ respectively, then we have,
        \begin{eqnarray*}
            \norm{\dot{U}(t)}_F = \norm{f_1}_F\cos \theta_1,\\
            \norm{\dot{V}(t)}_F = \norm{f_2}_F\cos \theta_2,
        \end{eqnarray*}
        for some $\theta_1, \theta_2\in\left[0, \dfrac{\pi}{2}\right)$.\\
        Henceforth, by choosing $\delta = \max\left\{\dfrac{1}{\cos^2 \theta_1}, \dfrac{1}{\cos^2 \theta_2}\right\}$, and then we have
        \begin{eqnarray*}
            \norm{\left(\dfrac{\partial F}{\partial U}, \dfrac{\partial F}{\partial V}\right)}_F\norm{\left(\dot{U}(t),
            \dot{V}(t)\right)}_F&\leq&\dfrac{\delta+1}{2}\left(\norm{\dot{U}(t)}_F^2+\norm{\dot{V}(t)}_F^2\right).
        \end{eqnarray*}
        \begin{eqnarray*}
            \inner{\left(\dfrac{\partial F}{\partial U}, \dfrac{\partial F}{\partial V}\right), \left(\dot{U}(t), \dot{V}(t)\right)} &\leq& -\left(\norm{\dot{U}(t)}_F^2+\norm{\dot{V}(t)}_F^2\right),\\
            &\leq& -\dfrac{2}{1+\delta} \norm{\left(\dfrac{\partial F}{\partial U}, \dfrac{\partial F}{\partial V}\right)}_F\norm{\left(\dot{U}(t), \dot{V}(t)\right)}_F.
        \end{eqnarray*}
        To prove the weak decrease condition,
        \begin{eqnarray*}
            &&\inner{\left(\dfrac{\partial F}{\partial U}, \dfrac{\partial F}{\partial V}\right), \left(\dot{U}(t), \dot{V}(t)\right)} = 0, \\
            \Rightarrow &&-\norm{\dot{U}(t)}_F^2-\norm{\dot{V}(t)}_F^2 = 0,\\
            \Rightarrow &&\dot{U}(t) = \dot{V}(t) = 0.
        \end{eqnarray*}
        Finally, because $\norm{U}_F = \norm{V}_F = N\leq \infty$, we guarantee the existence of stationary points $U^*$ and $V^*$ from Lemma \ref{lemma:conv_lemma}.
    \end{proof}
    \subsection{Sum-to-One Property}
     The fundamental property of quantum mechanics requires eigenvalues in a density matrix must be positive and sum to one. Therefore the flow to $\theta_i$ must admit this rule and is generated by the idea in \cite{Moody2022}. \\
    Consider the constraint to $\theta_i$,
    \begin{equation}\label{eq:sum_to_1}
        \displaystyle\sum_{i=1}^N\theta_i(t) = 1\Rightarrow \displaystyle\sum_{i=1}^N\dfrac{d\theta_i(t)}{dt} = 0 \mbox{ for any } t\geq 0.
    \end{equation}
    Denote
    $e_i = \theta_i-\inner{(\mathbf{x_i}\otimes \mathbf{y_i})(\mathbf{x_i} \otimes \mathbf{y_i})^T,\rho}$ and restrict the initial eigenvalues $\left\{\theta_i(0)\right\}_{i=1}^N$ such that they are positive and satisfy $\displaystyle\sum_{i=1}^N\theta_i(0) = 1$.  The flow of $\theta_i$ is reformulated as
    \begin{equation}\label{eq:eigenflow}
        \dfrac{d\theta_i}{dt} = -\left(e_i-\dfrac{1}{N}\displaystyle\sum_{j=1}^N e_j\right),
    \end{equation}
    Finally, we have the flow,
    \begin{equation}\label{eq:flow}
        \left\{
            \begin{array}{ccl}
                \dfrac{d\theta_i}{dt} & = &-\left(e_i-\dfrac{1}{N}\displaystyle\sum_{j=1}^N e_j\right),  \\
                \dfrac{dU}{dt} & = &  -4\mbox{ skew }\left(\left(-\tilde{M}+\hat{M}\right)U^T\right)U,\\
                \\
                \dfrac{dV}{dt} & = & -4\mbox{ skew }\left(\left(-\tilde{N}+\hat{N}\right)V^T\right)V.
            \end{array}
        \right.
    \end{equation}
    We then must demonstrate the objective function \eqref{eq:objective_function} is decreasing through the trajectory of \eqref{eq:flow} with given initial data $\left(U(0), V(0), \Sigma(0)\right)$.
    \begin{lemma}
        The objective function \eqref{eq:objective_function} is decreasing along the trajectory of \eqref{eq:flow}.
    \end{lemma}
    \begin{proof}
        First, consider
        \begin{eqnarray*}
            \dfrac{dF}{dt} &=& \dfrac{\partial F}{\partial U}\dfrac{dU}{dt}+\dfrac{\partial F}{\partial V}\dfrac{dV}{dt}+\sum_{i=1}^N \dfrac{\partial F}{\partial \theta_i}\dfrac{d\theta_i}{dt}\\
            & \leq & -\dfrac{2}{1+\delta} \norm{\left(\dfrac{\partial F}{\partial U}_F, \dfrac{\partial F}{\partial V}\right)}\norm{\left(\dot{U}(t), \dot{V}(t)\right)}_F+\sum_{i=1}^N \dfrac{\partial F}{\partial \theta_i}\dfrac{d\theta_i}{dt},\\
        \end{eqnarray*}
        which follows from \textbf{Theorem~\ref{thm:stationary}} directly.\\
        Next, by Cauchy-Schwarz inequality,
        \begin{eqnarray*}
            \sum_{i=1}^N \dfrac{\partial F}{\partial \theta_i}\dfrac{d\theta_i}{dt} = -\left(\sum_{i=1}^N e_i^2-\dfrac{1}{N}\sum_{i=1}^N e_i\sum_{j=1}^N e_j\right) \leq 0.
        \end{eqnarray*}
        Finally, we conclude that the objective function $F$ decreases along the flow \eqref{eq:flow}.
    \end{proof}

\section{Numerical Implementation}
    Three examples were designed and are used to illustrate \textbf{the descent property of the objective function}, \textbf{the rank adjustment property}, and \textbf{the maintenance of the sum-to-one property}. In these examples, we outline the ability to \textbf{decompose a classical-classical state} and \textbf{approximate the nearest classical-classical state to a given state}. These examples were accomplished using Python 3.8 on a MacBook Pro laptop with an M1 core 3.2-GHz processor and 16 GB of RAM. The \textbf{scipy.integrate.solve_ivp} package was applied to solve the problem of our dynamic system. We generally use \textbf{RK45} as our solver; yet if the system is stiff, \textbf{Radau} is an alternative choice. The aforementioned package allows advance user-defined arguments; for example, \textbf{Events},  \textbf{AbsTol} = $10^{-12}$, and \textbf{RelTol} = $10^{-12}$ help us to determine when the integration should be terminated.\\
    Consider a general bipartite system in our examples, where  $\left\{\mathbf{x}_i\right\}_{i=1}^N\in\mathbb{R}^{n}$ and $\left\{\mathbf{y}_i\right\}_{i=1}^N\in\mathbb{R}^{m}$. $m,n$ are arbitrary positive integers denoting the size of the system, and $N$ with $m,n\geq N$ is a positive integer denoting the number of qubits. \\
    In our examples, we consider general cases and remove the requirement that $n, m$ be the power of 2.\\

    \subsection{Example 1}
    This example examines the descent property of the objective function and the ability to decompose a given classical-classical state. A rank-3 synthetic classical-classical state $\rho\in \mathbb{R}^{16}\otimes \mathbb{R}^{8}$ is prepared. The state is composed of two random orthogonal matrices $\mathbf{U}\in\mathcal{S}_{16, 3}$ and  $\mathbf{V}\in\mathcal{S}_{8, 3}$  and $\mathbf{\Sigma} = \mbox{ diag}(\theta_1, \theta_2, \theta_3),$ which stores three random positive eigenvalues and satisfies the sum-to-one property:\\
    \begin{equation*}
        \rho = \left(\mathbf{U}\odot\mathbf{V}\right)\mathbf{\Sigma}\left(\mathbf{U}\odot\mathbf{V}\right)^T.
    \end{equation*}
    \begin{figure}
        \centering
       \begin{subfigure}[htbp]{.6\textwidth}
            \centering\includegraphics[width=\textwidth]{example 1/example_1_reduction_log.png}

      \end{subfigure}
      \caption{Descent of objective values}
      \label{fig:descent_obj_val}
    \end{figure}
    \begin{figure}
        \centering
      \begin{subfigure}[htbp]{.6\textwidth}
            \centering\includegraphics[width=\textwidth]{example 1/example_1_zoomin_logit.png}

      \end{subfigure}
      \caption{Evolution of initially guessed eigenvalues}
      \label{fig:eigen_flow}
    \end{figure}

    We assume that we have no information on the rank of the synthetic state, and we determine it by our method. The initial data of $U\mbox{ and } V$ are generated by applying singular value decomposition to two random matrices of size $16\times 8$ and $8 \times 8$ respectively, and eight random positive numbers are the initial eigenvalues which satisfy $\sum_{i=1}^8\theta_i = 1$. The objective values and the evolution of the eigenvalues are recorded in each iteration. The objective function decreases along the flow \eqref{eq:flow} in \figref{fig:descent_obj_val}.
    During the iteration, discard the eigenvalue when it evolves to zero. The iteration is terminated when either the condition \textbf{RelTol} or \textbf{AbsTol} is achieved or the integration successfully reaches the end of the time interval.
    Finally, we observe the objective value nearly descends to zero, which indicates that our algorithm can sufficiently decompose a given classical-classical state. \\
    The evolution of eigenvalues is plotted in \figref{fig:eigen_flow}; each curve reports an evolution process. For \figref{fig:eigen_flow}, eight eigenvalues were prepared initially for the algorithm, but five were discarded during the evolution, leaving only three eigenvalues when terminated. The small window in \figref{fig:eigen_flow} shows the scenarios in which eigenvalues reach zero.\\
    Based on our results, we conclude that our method can descend the objective function and correctly decompose a given classical-classical state.
    \subsection{Example 2}\label{example2}
    The second example clarifies the algorithm's preservation of the sum-to-one property. We arrange two states; both are generated randomly and are symmetric positive definite matrices with traces equal to 1. One state is of full rank, and the other is of rank 3.
    We denote $\rho_e$ and $\rho\in \mathbb{R}^{40\times 40}$ as full rank and rank 3 matrices, respectively. We approximate these matrices by $U\in \mathcal{S}_{8, 5} $ and$\ V\in \mathcal{S}_{5, 5}$ , and $\left\{\theta_i\right\}_{i=1}^5$.

    \begin{figure}
       \begin{subfigure}[htbp]{.44\textwidth}
            \centering\includegraphics[width=\textwidth]{example 2/example_2_full_error_log.png}
            \caption{Descent property of objective values in each trial of the full-rank matrix}
            \label{fig:descent_obj_trial_full}
      \end{subfigure}
      \hfill
      \begin{subfigure}[htbp]{.44\textwidth}
            \centering\includegraphics[width=\textwidth]{example 2/example_2_full_eigen_sum_to_one.png}
            \caption{Sum-to-one preservation of the full-rank matrix}
            \label{fig:eigen_sum_full}
      \end{subfigure}
      \caption{}
    \end{figure}

    \begin{figure}
       \begin{subfigure}[htbp]{.44\textwidth}
            \centering\includegraphics[width=\textwidth]{example 2/example_2_non_full_error_log.png}
            \caption{Descent property of objective values in each trial of the non-full-rank matrix}
            \label{fig:descent_obj_trial_non_full}
      \end{subfigure}
      \hfill
      \begin{subfigure}[htbp]{.44\textwidth}
            \centering\includegraphics[width=\textwidth]{example 2/example_2_non_full_eigen_sum_to_one.png}
            \caption{Sum-to-one preservation of the non-full-rank matrix}
            \label{fig:eigen_sum_non_full}
      \end{subfigure}
      \caption{}
    \end{figure}

    The figures of objective values are plotted in the log-log scale.\\
    First, the descent property and the consistency can be verified in \figref{fig:descent_obj_trial_full} and \figref{fig:descent_obj_trial_non_full}. The validation of sum-to-one preservation is presented in \figref{fig:eigen_sum_full} and \figref{fig:eigen_sum_non_full}.
    The results in \figref{fig:descent_obj_trial_full} and \figref{fig:descent_obj_trial_non_full} indicate that the objective function exhibits descending property. Furthermore, we reveal the consistency of our method by demonstrating that the optimal objective values approach nearly the same value with a small standard deviation. This evidence indicates that the minimizer of \eqref{eq:objective_function} can be stably detected through our method.\\
    On the other hand, \figref{fig:eigen_sum_full} and \figref{fig:eigen_sum_non_full} plot the sum of the eigenvalues of each iterative trial of the rank 3 matrix. These figures present oscillations where these outcomes are generated by floating point precision and numerical solvers during numerical integration. However, the oscillations are on the scale of $10^{-12}$; we can still reasonably claim that the sum-to-one property is maintained.\\

    \subsection{Example 3}
    The third example suggests the size and rank of the initial data in the dynamic system. \\
    Consider an arbitrary state $\rho$ of a bipartite system, namely, a positive definite matrix with its trace equal to 1. Let $\rho\in \mathbb{R}^{N\times N}$, and let $N$ be a non-prime number. Define a set $F$ that collects the factors of $N$, $F = \left\{(n, m)|2\leq n,m\in \mathbb{N,} \mbox{ and } nm = N\right\}$ and $r$, the rank of matrices; this set is required to satisfy $r\leq \min\left\{n, m\right\}$. Therefore, initial data for the dynamic system can be chosen from the following collection:
    \begin{eqnarray*}
        \left\{U\in \mathcal{S}_{n, r}, V\in \mathcal{S}_{m, r}|(n, m)\in F, r\leq n\right\}.
    \end{eqnarray*}

    \begin{figure}
       \begin{subfigure}[htbp]{\textwidth}
            \centering\includegraphics[width=\textwidth]{example 3/example_3_rank_adj_fwd.png}
            \caption{Objective Values with Different Matrix Sizes and Ranks}
            \label{fig:matrix_size_fwd}
      \end{subfigure}
    \hfill
       \begin{subfigure}[htbp]{\textwidth}
            \centering\includegraphics[width=\textwidth]{example 3/example_3_rank_adj_bkd.png}
            \caption{Objective Values with Different Matrix Sizes and Ranks}
            \label{fig:matrix_size_bkd}
      \end{subfigure}
    \end{figure}

    We select the initial guesses from the collection and plot the results in \figref{fig:matrix_size_fwd} and \figref{fig:matrix_size_bkd}. The figures indicate that the best approximation is obtained when $U\in \mathcal{S}_{15,4}$ and $V\in\mathcal{S}_{4,4}$. \\

    From this example, if a full-rank matrix is in one of the initial data, then it will produce the minimal objective value because each case in \figref{fig:matrix_size_fwd} and \figref{fig:matrix_size_bkd} indicates that the optimal objective value decreases when the rank of the initial data increases.
\section{Conclusion}
Identifying the nearest classical-classical state offers a means of understanding the degree of quantum correlation, even if quantifying the exact value is an NP-hard problem. In this study, we consider an approximation method to gauge the quantity. The approximation is resolved by solving a gradient-driven descent flow on Stiefel manifolds. Thanks to numerous existing ODE solvers, the aforementioned flow can be solved efficiently. At the end of this article, we discuss the convergence and descent properties theoretically and further provide three examples to illustrate the descent property, the preservation of a sum-to-one property, and the consistency of the objective value.\\

\begin{small}
\textbf{Acknowledgement}\\
This research was supported in part
by the National Center for Theoretical Sciences of Taiwan and
by the Ministry of Science and Technology of Taiwan
under grant 110-2636-M-006-006 and grant 111-2115-M-006-019. We also acknowledge Science College of National Cheng Kung University (NCKU Science) and Ministry of Science and Technology, Taiwan for a fellowship to support BZL’s PhD study
\end{small}
\bibliographystyle{elsarticle-num}


\section{Introduction}\label{sec1}

In physics equations of motion are derived by using the Lagrangian formalism,
which requires knowledge of a function called Lagrangian. By specifying a
Lagrangian for a given system, its equation of motion is derived by substituting
this Lagrangian into the Euler-Lagrange (E--L) equation. There are three different
families of Lagrangians, namely, standard, non-standard, and null.

Standard Lagrangians represent the difference between the kinetic and potential
energy [1], and they play a central role in Classical Mechanics [2-4] as well as
other areas of physics [5]. Non-standard Lagrangians (NSLs) are different in form
from standard Lagrangians but also give the same equation of motion as the NSLs.
The main difference is that NSLs lack terms that clearly discern energy-like forms [6].
There are also null Lagrangians and they identically satisfy the E-L equation and
can be expressed as the total derivative of any scalar function [7].  Since null
Lagrangians have no effect on equations of motion, they are not discussed any
further in this paper.

Different methods to derive standard [8-15] and nonstandard [12-26] Lagrangians
were developed. The methods were used to obtain both types of Lagrangians for
many different physical systems. In the first attempts to obtain Lagrangians for
selected biological systems, the Lagrangians were found by guessing their forms
[27-29]. Then, they were derived using the Method of Jacobi Last Multiplier [30,31].
More recently, the first standard Lagrangians were derived for five population
dynamics models and their implications for the models were discussed [32].

This paper aims to extend this work by deriving NSLs for selected population
dynamics models considered in [32], which are: the Lotka-Volterra [33,34],
Verhulst [35], Gompertz [36], Host-Parasite [37], and SIR [38] models. The
models are described by second-order nonlinear and damped ODEs, for which a
new method is required to find NSLs. In this paper, a novel method is developed
and then used to derive NSLs for all considered population models. The derived
NSLs are compared to those previously obtained [29,30,32] and used to gain new
insights into the role they play in the models, and in the symmetries underlying
them. For the derived NSLs, their corresponding null Lagrangians (NLs) and gauge
functions are also calculated, and these NLs in the population dynamics
is discussed.

The paper is organized as follows. Section 2 presents our method to construct
non-standard Lagrangians and its limitations; in Section 3, the population
dynamics models are described and our method to derive their non-standard
Lagrangians for these models are presented; then, non-standard Lagrangians
for the population dynamics models are given in Section 4; null Lagrangians
and their gauge functions are derived and discussed in Section 5; the obtained
results are compared to those previously obtained and discussed in Section 6;
and our Conclusions are given in Section 7.

\section{Lagrangian formalism}\label{sec2}

\subsection{Action and non-standard Lagrangians}\label{sec2.1}

The functional $\mathcal{S}[x(t)]$ is called action and is defined by
an integral over a scalar function $L$ that depends on the differentiable
function $x(t)$ that describes the time evolution of any dynamical system,
and on the time derivative $\dot x(t) = dx/dt$.  The function $L(\dot{x},
x)$ is called the Lagrangian function or simply {\it Lagrangian}, and it
is a map from the tangent bundle $T\mathcal{Q}$ to the real line
$\mathbf{R}$, or $\mathcal{L}:T\mathcal{Q} \rightarrow \mathbf{R}$,
with $\mathcal{Q}$ being a configuration manifold [3].   In general,
the Lagrangian may also depend explicitly on $t$, which can be written
as $L(\dot{x},x,t)$ and requires $\mathcal{L}:T \mathcal{Q}\times
\mathbf{R}\rightarrow \mathbf{R}$.

According to the principle of least action, or Hamilton's principle [2-6],
the action $\mathcal{S}[x(t)]$ must obey the following requirement
$\delta\mathcal{S}=0$, which guarantees that the action is stationary, or
has either a minimum, maximum, or saddle point.  The necessary condition
that $\delta\mathcal{S} = 0$, is known as the E--L equation, whose operator
$\hat {EL}$ acts on the Lagrangian, $\hat {EL} [ \mathcal{L} (\dot x, x, t) ]
= 0$ and gives a second-order ODE that becomes an equation of motion for
a dynamical system with $\mathcal{L} (\dot x, x, t)$. The process of deriving
the equation of motion is called the Lagrangian formalism [2-4], and it has
been extensively used in physics to derive its fundamental classical and
quantum equations [5].

The Lagrangian formalism is valid for both standard and non-standard
Lagrangians but only the latter are considered in this paper.  There are
different forms of NSLs [12-22] but only two of them are considered in
this paper.  The first considered form of NSL is

\begin{equation}
L_{ns1} (\dot x,x,t)=\frac{1}{A(t)\dot x+B(t)x+C(t)}\ ,
\label{eq1a}
\end{equation}
and it is applied to the equation of motion given by
\begin{equation}
\ddot x+a(t)\dot x^2+b(t)\dot x+c(t)x= 0\ ,
\label{eq1b}
\end{equation}
where the coefficients $a(t)$, $b(t)$ and $c(t)$ are at least twice differentiable
functions of the independent variable $t$, and the functions $A(t)$, $B(t)$
and $C(t)$ are expressed in terms of these coefficients by substituting
$L_{ns1} (\dot x,x,t)$ into the E--L equation.

The second form of the considered NSL is
\begin{equation}
L_{ns2}(\dot x,x)=\frac{1}{f(x)\dot x+g(x)x+h(x)}\ ,
\label{eq2a}
\end{equation}
and it is used for the following equation of motion
\begin{equation}
\ddot x+\alpha (x)\dot x^2+\beta(x)\dot x+\gamma(x)x=0\ ,
\label{eq2b}
\end{equation}
where the coefficients $\alpha(x)$, $\beta(x)$ and $\gamma(x)$
are at least twice differentiable functions of the dependent variable
$x(t)$, and the functions $f(x)$, $g(x)$ and $h(x)$ are determined
by substituting $L_{ns2} (\dot x,x)$ into the E--L equation, which
allows expressing these functions in terms of the coefficients.

The NSL $L_{ns1} (\dot x,x,t)$ was extensively studied and applied
to several dynamical systems [12,13,20].  However, studies of
$L_{ns2}(\dot x,x)$ were limited to only some special cases [20].
Therefore, one of the objectives of this paper is to perform a detailed
analysis of $L_{ns2}(\dot x,x)$ and its applicability to the population
dynamics models; so far, only other forms of NSLs been applied to
the models [27-31].  Moreover, in a recent work [32], the standard
Lagrangians were constructed for the same sample of the population

\subsection{Limits on construction of non-standard Lagrangians}\label{sec2.2}

To obtain a NSL means to determine the unknown functions in
Eqs (\ref{eq1a}) and (\ref{eq2a}), which requires that an equation
of motion is given.  As shown in [32], the equations of motion
for the population dynamics models can be written in the following form
\begin{equation}
\ddot x+\alpha (x)\dot x^2+\beta(x)\dot x+\gamma(x)x= F(x)\ ,
\label{eq3}
\end{equation}
which generalizes Eq. (\ref{eq2b}) by taking into account a driving
force $F(x)$.  For different models, the coefficients $\alpha(x)$,
$\beta(x)$ and $\gamma(x)$ are given, and $F(x)$ is also known.
Since Eq. (\ref{eq3}) does not depend explicitly on time, the
Lagrangian that can be used for this equation is $L_{ns2}(\dot x,x)$
given by Eq. (\ref{eq2a}).  From now on, we shall take $L_{ns2}
(\dot x,x)\equiv L_{ns}(\dot x,x)$ and use the latter throughout
this paper.

It is seen that Eq. (\ref{eq3}) has both linear and quadratic damping
terms, that can become nonlinear depending on the form of
$\gamma(x)x$, and that it is driven.  In a previous study [32], the
standard Lagrangian for Eq. (\ref{eq3}) was found in case $ F(x)$
= const.  To construct NSL for Eq. (\ref{eq3}), we take $L_{ns}(\dot x,x)$
given by Eq. (\ref{eq2a}) and determine the functions $f(x)$, $g(x)$
and $h(x)$ in terms of the known coefficients $\alpha(x)$, $\beta(x)$
and $\gamma(x)$.  Our approach generalizes the previous work [18],
and it allows us to investigate the applicability of $L_{ns}(\dot x,x)$
to Eq. (\ref{eq3}) and its limitations.\\ \\ \\
{\bf Proposition 1:}  Let a general non-standard Lagrangian be given by
\begin{equation}
L_{ns}(\dot x,x)=\frac{1}{f(x)\dot x + G(x)}\ ,
\label{eq4}
\end{equation}
where $G(x) = g(x)x+h(x)$.  Then, $L_{ns}(\dot x,x)$ becomes the
NSL for Eq. (\ref{eq3}) if, and only if,
\begin{equation}
f (x) = e^{\int^x \alpha (\xi) d \xi} \equiv e^{I_{\alpha} (x)}\ ,
\label{eq5}
\end{equation}
and
\begin{equation}
G (x) = 3 \left [ \frac{\gamma (x) x -  F(x)}{\beta (x)} \right ]
e^{I_{\alpha} (x)}\ ,
\label{eq6}
\end{equation}
and either one of the following conditions
\begin{equation}
\beta (x) = \frac{9}{2} e^{- I_{\alpha} (x)} \frac{d}{dx} \left [
\left ( \frac{\gamma (x) x -  F(x)}{\beta (x)} \right )
e^{I_{\alpha} (x)} \right ]\ ,
\label{eq7}
\end{equation}
and
\begin{equation}
\gamma (x) x - F(x) = \frac{2}{9} \beta (x) e^{- I_{\alpha} (x)}
\int^x \beta (\xi) e^{I_{\alpha} (\xi)} d\xi\ ,
\label{eq8}
\end{equation}
is satisfied.\\ \\
{\bf Proof:} Substituting Eq. (\ref{eq3}) into the E--L equation,
we obtain
\begin{equation}
\ddot x + \left [ \frac{f^{\prime (x)}}{f(x)} \right ] \dot x^2 +
\left [ \frac{3}{2f (x)} \frac{G^{\prime} (x)}{f(x)} \right ] \dot x
+ \left [ \frac{1}{2 x} \frac{G^{\prime} (x)}{f(x)} \frac{G (x)} {f (x)}
\right ] x = 0\ .
\label{eq9}
\end{equation}
By comparing this equation to Eq. (\ref{eq3}), we find the following
relationships
\begin{equation}
\alpha (x) = \frac{f^{\prime (x)}}{f(x)}\ ,
\label{eq10a}
\end{equation}
\begin{equation}
\beta (x) = \frac{3}{2f (x)} \frac{G^{\prime} (x)}{f(x)}\ ,
\label{eq10b}
\end{equation}
and
\begin{equation}
\gamma (x) = \frac{1}{2 x} \frac{G^{\prime} (x)}{f(x)}
\frac{G (x)}{f(x)} + \frac{F (x)}{x}\ .
\label{eq10c}
\end{equation}
Using Eq. (\ref{eq10a}), we obtain $f(x)$ given by Eq. (\ref{eq5}).
However, by combining Eqs (\ref{eq10b}) and (\ref{eq10c}), and
using Eq. (\ref{eq10a}), we find $G(x)$ given by Eq. (\ref{eq6}).
Then, the conditions expressed by Eqs (\ref{eq7}) and (\ref{eq8})
are easy to derive from Eqs (\ref{eq10a}), (\ref{eq10b}) and
(\ref{eq10c}).  Moreover, substituting Eq. (\ref{eq8}) into Eq.
(\ref{eq7}) shows that both conditions are equivalent, thus, it is
sufficient to use only one of them. This concludes the proof.\\ \\
{\bf Corollary 1:} The coefficients $\beta (x)$ and $\gamma (x)$
mutually depend on each other through the function $G(x)$ and,
in addition, they also depend on $\alpha (x)$ through the function
$f(x)$.\\ \\
{\bf Corollary 2:} The conditions given by Eqs (\ref{eq7}) and
(\ref{eq8}) are equivalent, which means that if one of them is not
satisfied then $L_{ns}(\dot x,x)$, with $f(x)$ and $G(x)$ (see
Eqs \ref{eq5} and \ref{eq6}), cannot be considered to be the NSL
for Eq. (\ref{eq3}).\\ \\
{\bf Corollary 3:} Only $\alpha (x)$ is required to uniquely
determine the function $f(x)$.

\subsection{From nonlinear to linear equations}

The main result of Proposition 1 is that any NSL of the form of Eq.
(\ref{eq4}) can only be constructed if one of the conditions given by
Eqs (\ref{eq7}) and (\ref{eq8}) is satisfied.  We now explore other
consequences of the condition given by Eq. (\ref{eq8}) in the
following proposition.\\
{\bf Proposition 2:}  Let the equation of motion be
\begin{equation}
\ddot x+\alpha (x)\dot x^2+\beta(x)\dot x+H(x) = 0\ ,
\label{eq11}
\end{equation}
where $H (x) = \gamma (x) x  - F (x)$, and let
\begin{equation}
z (x) = \int^x \beta (\xi) e^{I_{\alpha} (\xi)} d\xi\ ,
\label{eq12}
\end{equation}
be a new dependent variable that can be expressed in terms of
a new variable $\eta$ that is related to the original variable $x$
by $d \eta = \beta (x) dx$.  Then, Eq. (\ref{eq11}) becomes
\begin{equation}
z^{\prime \prime} + z^{\prime} + z = 0\ ,
\label{eq13}
\end{equation}
with $z^{\prime} = dz / d\eta$ and $z^{\prime \prime} = d^2 z /
d\eta^2$, if, and only if,
\begin{equation}
H (x) = \frac{2}{9} \beta (x) e^{- I_{\alpha} (x)}
\int^x \beta (\xi) e^{I_{\alpha} (\xi)} d\xi\ .
\label{eq14}
\end{equation}
{\bf Proof:} Using Eq. (\ref{eq12}), we find
\begin{equation}
z^{\prime} = \dot x e^{I_{\alpha} (x)}\ ,
\label{eq15a}
\end{equation}
and
\begin{equation}
z^{\prime \prime} = \frac{1}{\beta(x)} \left [ \ddot x + \alpha (x)
\dot {x}^2 \right ] e^{I_{\alpha} (x)}\ .
\label{eq15b}
\end{equation}
Substituting Eqs (\ref{eq15a}), (\ref{eq15b}) and (\ref{eq12}) into
Eq. (\ref{eq13}), we obtain the original equation of motion (see Eq.
\ref{eq11}) with the function $H(x)$ given by Eq. (\ref{eq14}).
This concludes the proof.\\ \\
{\bf Corollary 4:} Equation (\ref{eq13}) is the Sturm-Liouville equation
[43] whose solutions are well-known and given as $z (\eta) = c_1 [ \sin
(\sqrt{3} \eta / 2) + c_2 \cos (\sqrt{3} \eta / 2) ] e^{-\eta /2}$.\\ \\
{\bf Corollary 5:} The non-standard Lagrangian given by Eq. (\ref{eq4})
can be constructed without any limits for those nonlinear ODEs that can
be converted into linear ODEs, which requires that the condition given
by Eq. (\ref{eq8}) is satisfied.\\

\subsection{Non-standard Lagrangians without limits}

The results of Proposition 2 and Corollaries 4 and 5 demonstrate that
the existence of NSL of the form of Eq. (\ref{eq4}) is determined by
the conditions given by Eqs (\ref{eq7}) and (\ref{eq8}), which must
be satisfied in order for the NSL to exist.  According to Corollaries 1-3,
the only two terms in the equation of motion given by Eq. (\ref{eq3})
that can be uniquely derived without any limits by the NSL are the terms
with $\ddot x (t)$ and $\dot x^2 (t)$.  The coefficients in terms with
$\dot x(t)$ and $x (t)$ are mutually related, thus, they are dependent
on each other as shown by either Eq. (\ref{eq7}) or Eq. (\ref{eq8}).

Because of these limitations on the construction of NSLs for Eq. (\ref{eq3}),
let us now propose another method by extending the previous work [32].
The basic idea of this work is to write Eq. (\ref{eq3}) in the following form
\begin{equation}
\ddot x+\alpha (x)\dot x^2 = {\mathbb F}_{dis} (\dot x, x)\ ,
\label{eq16}
\end{equation}
where
\begin{equation}
{\mathbb F}_{dis} (\dot x, x) = F(x) - \beta(x)\dot x - \gamma(x)x\ ,
\label{e17}
\end{equation}
becomes a dissipative force because of its dependence on $\dot x (t)$.
In this case, the E--L equation [2,3,32] can be written as
\begin{equation}
{d\over dt}\biggr({\partial L\over\partial\dot x}\biggr)-{\partial L
\over\partial x}={\mathbb F}_{dis} (\dot x,x)e^{2I_{\alpha}(x)}\ ,
\label{eq18}
\end{equation}
where the force on the RHS of this equation is known at the Rayleigh
force [2,3]

Based on the results presented in Section 2.2, the NSL for Eq. (\ref{eq11})
is
\begin{equation}
L_{ns}(\dot x,x)=\frac{1}{\dot x e^{I_{\alpha} (x)} + C_o}\ ,
\label{eq19}
\end{equation}
where the constant $C_o$ replaces the function $G (x)$ in Eq. (\ref{eq4});
the constant can have any real value, and it is not required to determine
its value.  This non-standard Lagrangian has no restrictions or limitations
and it exists for any differentiable $\alpha (x)$ regardless of the forms of
$\beta (x)$ and $\gamma (x)$.  Therefore, $L_{ns}(\dot x,x)$ will be
used to find NSLs for the population dynamics models (see Section 3).
Moreover, the validity of one of the conditions given by Eqs (\ref{eq7})
and (\ref{eq8}) must be verified for all the models to determine whether
any of these models allows for the NSL given by Eq. (\ref{eq4}).

\section{Population dynamics models and methods}

\subsection{Selected models}

We consider the following population dynamics models: the Lotka-Volterra
[33,34], Verhulst [35], Gompertz [36] and Host-Parasite [37] models that
describe two interacting (preys and predators) species.  However, the SIR
model [38] describes the spread of a disease in a given population.

\begin{table}
\caption{Biological models and their equations of motion.\label{tab1}}
\centering
\begin{tabular}{llccc}
\hline\hline
{\bf Population models}	&{\bf Equations of Motion}\\
\hline
Lotka-Volterra Model&$\dot{w_1}=w_1(a\:+\:bw_2\:)$\\
                    &$\dot{w_ 2}=w_2(A+Bw_1)$\\ \hline
Verhulst Model      &$\dot{w}_1=w_1(A+\:Bw_1\:+f_1 w_2)$\\
                    &$\dot{w}_2=w_2(\:a+\:\:bw_2\:\:+f_2 w_1)$\\ \hline
Gompertz Model      &$\dot{w}_1=w_1{(}A\log{\bigr(}{w_1\over m_1}{\bigr)}+B{w_2}{)}$\\
                    &$\dot{w}_2=w_2{(}\:a\log{\bigr(}{w_2\over m_2}{\bigr)}+\:b{w_1})$\\ \hline
Host-Parasite Model &$\dot{w}_1=w_1(a\:-\:\:bw_2)$\\
			    &$\dot{w}_2=w_2(A-B{w_2\over w_1})$ \\ \hline
SIR Model           &$\dot w_1=-b w_1w_2$\\
			      &$\dot w_2=\:\:\:b w_1w_2-aw_2$\\
\hline
\end{tabular}
\end{table}
\bigskip
The variables $w_1 (t)$ and $w_2 (t)$ in the first four models of Table 1
represent prey and predators of the interacting species, respectively.
Moreover, the time derivatives $\dot w_1 (t)$ and $\dot w_2 (t)$ describe
changes of these species in time.  The interaction between the species in
each model is given by the coefficients $a$, $A$, $b$, $B$ $f_1$, $f_2$,
$m_1$ and $m_2$, which are real and constant.  However, in the SIR
model, the variables $w_1(t)$ and $w_2(t)$ describe susceptible and
infectious populations, with $a$ and $b$ being the recovery and infection
rates, respectively.

Among the models given in Table 1, the Lotka-Volterra, Verhulst and
Gompertz models are symmetric, and the remaining two models are
asymmetric, where being symmetric means that the dependent variables
can be replaced by each other, and the same can be done with the
coefficients. Obviously, this cannot be done for the asymmetric models.

\subsection{Methods to construct non-standard Lagrangians}

According to Table 1, each model is described by a coupled first-order
nonlinear ODEs that can be cast into a second-order nonlinear ODE of
the following form
\begin{equation}
\ddot w_{i}+\alpha_{i}(w_{i})\dot w_{i}^2+\beta_{i}(w_{i})\dot w_{i}
+ \gamma_{i}(w_{i})w_{i}=F_{i}(w_i)\ ,
\label{eq20}
\end{equation}

where $i=1$ and $2$.  This equation is of the same form as that given
by Eq. (\ref{eq3}).  The non-standard Lagrangian for this equation is
\begin{equation}
L_{ns}(\dot w_i,w_i)=\frac{1}{f(w_i)\dot w_i + G(w_i)}\ ,
\label{eq21}
\end{equation}
with the functions $f(w_i)$ and $G(w_i)$ being given by Eqs (\ref{eq5})
and (\ref{eq6}), respectively.  Thus, for each considered model and for
each variable in this model, we determine these functions and obtain
$L_{ns}(\dot w_i,w_i)$.

Then, we verify the validity of one of the conditions given by Eqs (\ref{eq7})
and (\ref{eq8}).  Since these conditions are equivalent, we take only one of
them and select the condition on $\gamma (w_i)$ (see Eq. \ref{eq8}).  The
explicit form of this condition used in our calculations is
\begin{equation}
\gamma (w_i) w_i - F_i (w_i) = \frac{2}{9} \beta (w_i) e^{- I_{\alpha} (w_i)}
\int^{w_i} \beta (\xi) e^{I_{\alpha} (\xi)} d\xi\ ,
\label{eq22}
\end{equation}
This condition is used to verify whether the derived NSL (see Eq. \ref{eq4})
can be constructed for Eq (\ref{eq3}) or not (see Proposition 1).  Moreover,
the validity of this condition determines uniquely that the equation of motion
can be converted into a linear second-order ODE, whose solutions are easy to
find (see Proposition 2).

As shown in Section 2.3, another method to construct NSLs is to cast Eq.
(\ref{eq20}) into the form
\begin{equation}
\ddot w_{i}+\alpha_{i}(w_{i})\dot w_{i}^2 ={\mathbb F}_{dis, i}(\dot w_{i},
w_{i})\ ,
\label{eq23}
\end{equation}
where
\begin{equation}
{\mathbb F}_{dis,i}(\dot w_{i}, w_{i}) = F_i (w_i) - \beta_{i}(w_{i})\dot w_{i}
- \gamma_{i}(w_{i})w_{i}\ .
\label{eq24}
\end{equation}
Then, according to Eq. (\ref{eq14}), the NSL for Eq. (\ref{eq23}) is given by
\begin{equation}
L_{ns,i}(\dot w_i,w_i)=\frac{1}{\dot w_i e^{I_{\alpha} (w_i)} + C_o}\ .
\label{eq25}
\end{equation}
This non-standard Lagrangian is not constrained by any conditions and
it can be derived for all considered population dynamics models.  If this
NSL is substituted into the following EL equation
\begin{equation}
{d\over dt}\biggr({\partial L_{ns,i}\over\partial\dot w_{i}}\biggr)-{\partial
L_{ns,i}\over\partial w_{i}}={\mathbb F}_{dis,i} (\dot w_{i}, w_{i})
e^{2I_{\alpha}(w_i)}\ ,
\label{eq26}
\end{equation}
then, the equation of motion for the considered model is obtained.
In the following, this method is used to construct the NSLs for all population
dynamics models considered in this paper.

\section{Models and their non-standard Lagrangians}

\subsection{\bf Lotka-Volterra Model}

The model is symmetric and it is represented mathematically by a system of
coupled nonlinear first-order ODEs given in Table 1.  We cast the first-order
ODEs into the second-order ODEs of the form given by Eq. (\ref{eq23}), and
obtain
\begin{equation}
\ddot w_1 - {1\over w_1} \dot w_1^2 = {\mathbb F}_{dis, 1}(\dot w_{1},
w_{1})\ ,
\label{eq27a}
\end{equation}
where
\begin{equation}
{\mathbb F}_{dis, 1}(\dot w_{1}, w_{1}) =(Bw_1+A) \dot w_1 +
a (Bw_1+A) w_1\ ,
\label{eq27b}
\end{equation}
and
\begin{equation}
\ddot w_2 - {1\over w_2} \dot w_2^2 = {\mathbb F}_{dis, 2}(\dot w_{2},
w_{2})\ ,
\label{eq27c}
\end{equation}
with
\begin{equation}
{\mathbb F}_{dis, 2}(\dot w_{2}, w_{2}) = (b w_2+a) \dot w_2 +
A (bw_2+a) w_2\ .
\label{eq27d}
\end{equation}
Using Eq. (\ref{eq5}), the factors $e^{I_{\alpha} (w_i)}$ for both models
can be calculated, and the obtained results are substituted into Eq. (\ref{eq25})
to give
\begin{equation}
L_{ns,1}(\dot w_1,w_1)=\frac{1}{\dot w_1 w_1^{-1} + C_o}\ .
\label{eq28a}
\end{equation}
and
\begin{equation}
L_{ns,2}(\dot w_2,w_2)=\frac{1}{\dot w_2 w_2^{-1} + C_o}\ ,
\label{eq28b}
\end{equation}
which are the NSLs for the Lotka-Volterra model.  It is easy to verify that
by substituting them into the E--L equation given by Eq. (\ref{eq26}),
the evolution equations describing the Lotka-Volterra model are obtained
(see Eqs \ref{eq27a} and \ref{eq27c}).

\subsection{\bf Verhulst Model}

The system of coupled nonlinear ODEs given in Table 1 shows
that the model is symmetric. The second-order equations for
the dynamical variables of this model are:
\begin{equation}
\ddot w_1 - (1 + b) {1\over w_1} \dot w_1^2 = {\mathbb F}_{dis, 1}(\dot w_{1},
w_{1})\ ,
\label{eq29a}
\end{equation}
where
\[
{\mathbb F}_{dis, 1}(\dot w_{1}, w_{1}) = -\dot w_1(2b-1)Bw_1 -
f_2w_1^2+(2Ab-a)+(f_2-b)B
\]
\begin{equation}
\hskip0.75in + [(f_2-b)Bw_1^2 + (Af_2-2Ab-a)w_1 + A(a-Ab)] w_1\ ,
\label{eq29b}
\end{equation}
and
\begin{equation}
\ddot w_2 - (1 + B) {1\over w_2} \dot w_2^2 = {\mathbb F}_{dis, 2}(\dot w_{2},
w_{2})\ ,
\label{eq29c}
\end{equation}
with
\[
{\mathbb F}_{dis, 2}(\dot w_{2}, w_{2}) = -\dot w_2(2B-1)bw_2 -
f_1w_2^2 + (2aB-A) + (f_1-B)b
\]
\begin{equation}
\hskip0.75in  + [(f_1-B) bw_2^2 + (af_1 - 2aB - A)w_2 + a(A - aB) ] w_2\ .
\label{eq29d}
\end{equation}
The non-standard Lagrangians for the evolution equations describing this model
are:
\begin{equation}
L_{ns,1}(\dot w_1,w_1)=\frac{1}{\dot w_1 w_1 ^ {-(1+b)} + C_0}\ .
\label{eq30a}
\end{equation}
and
\begin{equation}
L_{ns,2}(\dot w_2,w_2)=\frac{1}{\dot w_2 w_2 ^ {-(1+B)} + C_o}\ ,
\label{eq30b}
\end{equation}

The derived NSLs give the original evolution equations for the model
(see Eqs \ref{eq29a} and \ref{eq29c}) after they are substituted into
the E--L equation given by Eq. (\ref{eq26}).

\subsection{\bf Gompertz Model}
The mathematical representation of this model given by the coupled and nonlinear
ODEs in Table 1 show that the model is symmetric.

The equation describing the time evolution of each model variable is given as
\begin{equation}
\ddot w_1 - {1\over w_1} \dot w_1^2 = {\mathbb F}_{dis, 1}(\dot w_{1},
w_{1})\ ,
\label{eq31a}
\end{equation}
where
\[
{\mathbb F}_{dis, 1}(\dot w_{1}, w_{1}) = Am_1 + bw_1 + g_1(\dot w_1,w_1)
\dot w_1 - g_1(\dot w_1,w_1) Aw_1
\]
\begin{equation}
\hskip0.75in + A\log({w_1\over m_1})\ w^2_1\ ,
\label{eq31b}
\end{equation}
and
\begin{equation}
\ddot w_2 - {1\over w_2} \dot w_2^2 = {\mathbb F}_{dis, 2}(\dot w_{2},
w_{2})\ ,
\label{eq31c}
\end{equation}
with
\[
{\mathbb F}_{dis, 2}(\dot w_{2}, w_{2}) = am_2 + Bw_2 + g_2
(\dot w_2,w_2) \dot w_2-g_2(\dot w_2,w_2)aw_2 +
\]
\begin{equation}
\hskip0.75in a\log({w_2\over m_2})\ w^2_2 \ .
\label{eq31d}
\end{equation}
\begin{equation}
L_{ns,1}(\dot w_1,w_1)=\frac{1}{\dot w_1 w_1^{-1} + C_o}\ .
\label{eq32a}
\end{equation}
and
\begin{equation}
L_{ns,2}(\dot w_2,w_2)=\frac{1}{\dot w_2 w_2^{-1} + C_o}\ ,
\label{eq32b}
\end{equation}

To obtain the original evolution equations given by Eqs \ref{eq31a}
and \ref{eq31c}) for this model, it is necessary to substitute the
derived NSLs into the E--L equation (see Eq. \ref{eq26}).

\subsection{\bf Host-Parasite Model}

This model describes the interaction between a host and its parasite.
The model takes into account the non-linear effects of the host population
size on the growth rate of the parasite population [22].  The system
of coupled nonlinear ODEs (see Table 1) is asymmetric in the dependent
variables $w_1$ and $w_2$.  The time evolution equations for these
variables are:
\begin{equation}
\ddot w_1 - {1\over w_1}\Bigr(1+{B\over bw_1}\Bigr) \dot w_1^2 =
{\mathbb F}_{dis, 1}(\dot w_{1},
w_{1})\ ,
\label{eq33a}
\end{equation}
where
\begin{equation}
{\mathbb F}_{dis, 1}(\dot w_{1}, w_{1}) = B{a^2\over b} +
\Bigr(A-{2aB\over bw_1}\Bigr)\dot w_1 + aAw_1  \ ,
\label{eq33b}
\end{equation}
and
\begin{equation}
\ddot w_2 - {2\over w_2} \dot w_2^2 = {\mathbb F}_{dis, 2}(\dot w_{2},
w_{2})\ ,
\label{eq33c}
\end{equation}
with
\begin{equation}
{\mathbb F}_{dis, 2}(\dot w_{2}, w_{2}) =
(bw_2-a-A)\dot w_2 + A(bw_2-a) w_2 \ .
\label{eq33d}
\end{equation}

\begin{equation}
L_{ns,1}(\dot w_1,w_1)=\frac{1}{\dot w_1  w_1^{-1}e^{B \over w_1} + C_o}\ .
\label{eq34a}
\end{equation}
and
\begin{equation}
L_{ns,2}(\dot w_2,w_2)=\frac{1}{\dot w_2 w_2^{-2} + C_o}\ ,
\label{eq34b}
\end{equation}

Then, the original evolution equations for this model (see Eqs \ref{eq33a}
and \ref{eq33c}) is obtained by substituting the derived NSLs into the E--L
equation (see Eq. \ref{eq26}).

\subsection{\bf SIR Model}
The equation describing the time evolution of each model variable is given as
\begin{equation}
\ddot w_1 - {1\over w_1} \dot w_1^2 = {\mathbb F}_{dis, 1}(\dot w_{1},
w_{1})\ ,
\label{eq35a}
\end{equation}
where
\begin{equation}
{\mathbb F}_{dis, 1}(\dot w_{1}, w_{1}) = (bw_1 - a )\dot w_1 \ ,
\label{eq35b}
\end{equation}
and
\begin{equation}
\ddot w_2 - {1\over w_2} \dot w_2^2  = {\mathbb F}_{dis, 2}(\dot w_{2},
w_{2})\ ,
\label{eq35c}
\end{equation}
with
\begin{equation}
{\mathbb F}_{dis, 2}(\dot w_{2}, w_{2}) = -bw_2 \dot w_2 + abw_2^2\ .
\label{eq35d}
\end{equation}
Using Eq. (\ref{eq5}), the factors $e^{I_{\alpha} (w_i)}$ for both models
can be calculated, and the obtained results are substituted into Eq. (\ref{eq25})
to give
\begin{equation}
L_{ns,1}(\dot w_1,w_1)=\frac{1}{\dot w_1 w_1^{-1} + C_o}\ .
\label{eq36a}
\end{equation}
and
\begin{equation}
L_{ns,2}(\dot w_2,w_2)=\frac{1}{\dot w_2 w_2^{-1} + C_o}\ ,
\label{eq36b}
\end{equation}

The SIR model is asymmetric and the derived NSLs give the original
evolution equations for this model (see Eqs \ref{eq35a} and \ref{eq35c})
after the NSLs are substituted into the E--L equation given by Eq.
(\ref{eq26}).

\section{Null Lagrangians for the population models}

In our previous work [32], we showed how to construct standard Lagrangians
for the population dynamics models.  Moreover, in this paper, we constructed
the NSLs for the same models. However, there is another family of Lagrangians
called null Lagrangians (NLs), which make the E--L equation identically zero
and are given as the total derivative of a scalar function [39]; the latter is
called here a gauge function.  The NLs were extensively studied in mathematics
(e.g., [39-44]) and recently in physics [45-49], but to the best of our knowledge,
NSLs have not yet been introduced to biology and, specifically, to its population
dynamics.  Therefore, in the following, we present the first applications of NLs to
biology and its population dynamics.

Recent studies of null Lagrangians demonstrated that there is a different
condition that is obeyed by NLs and that this condition plays the same role
for NLs as the E--L equation plays for standard and non-standard Lagrangians
[48,49].  The condition can be written as
\begin{equation}\
     \frac{d L_{null,i}(\dot w_i,w_i)}{dt} =
     \frac{\partial L_{null,i}}{\partial t} + \dot w_i
     \frac{\partial L_{null,i}}{\partial w_i} + \ddot w_i
     \frac{\partial L_{null,i}}{\partial \dot w_i} = 0\ ,
\label{eq37}
\end{equation}
and it shows that the substitution of any NL into Eq. (\ref{eq37}) results in
an equation of motion.  However, the resulting equations of motion may be
limited because their coefficients are required to obey relationships that are
different for different equations of motion [48,49].  The previous work also
demonstrated that an inverse of any null Lagrangian generates a non-standard
Lagrangian, whose substitution into the E--L equation gives a new equation
of motion.  However, the reverse is not always true, which means that not
all NSLs have their corresponding NLs because of the so-called null condition
that must be satisfied [48,49].

For any non-standard Lagrangian of the form
\begin{equation}
L_{ns,i}(\dot{w_i},w_i,t) = \frac{1}{B_i (w_i,t) \dot{w_i} + C_i (w_i,t) w_i}\ ,
\label{eq38}
\end{equation}
the null condition [48,49] is
\begin{equation}
\left ( \frac{\partial B_i (w_i,t)}{\partial t} \right ) = \frac{\partial [w_i C_i (w_i,t)]}
{\partial w_i}\ .
\label{eq39}
\end{equation}
Comparison of Eq. (\ref{eq38}) to Eq. (\ref{eq25}) shows that $B_i (w_i) =
exp[ I_{\alpha}(w_i)]$ and $w_i C_i (w_i,t) = C_o$, which means that
the denominator of Eq. (\ref{eq25}) satisfies the null condition and, therefore,
the null Lagrangians for all the considered population dynamics systems
are of the form
\begin{equation}
L_{null,i}(\dot w_i,w_i)=\dot w_i e^{I_{\alpha} (w_i)} + C_o\ .
\label{eq40}
\end{equation}
It is easy to verify that substitution of $L_{null,i}(\dot w_i,w_i)$
into the E--L equation gives identically zero.

Now, to derive the equation of motion given by Eq. (\ref{eq26}),
the null condition of Eq. (\ref{eq37}) must be modified to
account for the dissipative force.  Then, the null condition is
\begin{equation}\
     \frac{d L_{null,i}(\dot w_i,w_i)}{dt} = {\mathbb F}_{dis,i}
(\dot w_{i}, w_{i}) e^{2I_{\alpha}(w_i)}\ .
\label{eq42}
\end{equation}
Substitution of $L_{null, i}(\dot w_i,w_i)$ into this equation gives
the required equation of motion (see Eq. \ref{eq26}).  In the
following, we use Eq. (\ref{eq40}) to find the NLs for all the
population dynamics models considered in this paper.

Another important characteristic of NLs is the fact that they
can always be expressed as the total derivative of a scalar
function [39], which has been called a gauge function [45-49].
Thus, we may write
\begin{equation}
L_{null,i}(\dot w_i,w_i) = \frac{d \Phi_i(w_i, t)}{dt}=
\frac{\partial \Phi_i}{\partial t} + \dot w_i \frac{\partial
\Phi_i}{\partial w_i}\ ,
\label{eq43}
\end{equation}
and use it to determine the gauge functions $\Phi_i(w_i, t)$
for all derived null Lagrangian $L_{null,i} (\dot w_i,w_i)$.

The results of our derivations of the null Lagrangians (see Eq.
\ref{eq40}) and their corresponding gauge functions (see Eq.
\ref{eq43}) are presented in Table 2. To the best of our knowledge,
these are the first null Lagrangians and gauge functions ever
presented for any biological systems, especially for the
population dynamics models.\\
\begin{table}
\begin{center}
\caption{Null Lagrangians and their gauge functions for the population models}\label{tab2}%
\begin{tabular}{@{}lll@{}}
\hline
{\bf Population models}	&{\bf Null Lagrangians } &{\bf Gauge functions}\\
\hline
Lotka-Volterra  &	$L_{null,1} = \dot{w_1} w_1^{-1} + C_0$
                &   $\Phi_{1} = \ln \vert w_1 \vert + C_0 t$ \\
			&	$L_{null,2} = \dot{w_2} w_2^{-1} + C_0$
                &   $\Phi_{2} = \ln \vert w_2 \vert + C_0 t$\\ \hline

Verhulst  	    &	$L_{null,1} = \dot w_1 w_1 ^ {-(1+b)} + C_0$
                &   $\Phi_{1} = - b^{-1} w_1^{-b} + C_0 t$ \\
			&	$L_{null,2} = \dot w_2 w_2 ^ {-(1+B)} + C_0$
                &   $\Phi_{2} = - B^{-1} w_1^{-B} + C_0 t$\\ \hline

Gompertz  	    & 	$L_{null,1} = \dot{w_1} w_1^{-1} + C_0$
                &   $\Phi_{1} = \ln \vert w_1 \vert + C_0 t$ \\
			&	$L_{null,2} = \dot{w_2} w_2^{-1} + C_0$
                &   $\Phi_{2} = \ln \vert w_2 \vert + C_0 t$\\ \hline

Host-Parasite	&	$L_{null,1} = \dot w_1  w_1^{-1}e^{B \over w1} + C_0$
                &   $\Phi_{1} = - Ei ({B \over w_1}) + C_0 t$ \\
			&	$L_{null,2} = \dot{w_2} w_2^{-2} + C_0$
                &   $\Phi_{2} = - w_2^{-1} + C_0 t$\\ \hline

SIR  		    &	$L_{null,1} = \dot{w_1} w_1^{-1} + C_0$
                &   $\Phi_{1} = \ln \vert w_1 \vert + C_0 t$ \\
	       	&	$L_{null,2} = \dot{w_2} w_2^{-1} + C_0$
                &   $\Phi_{2} = \ln \vert w_2 \vert + C_0 $t\\ \hline
\end{tabular}
\end{center}
\end{table}

\section{Discussion and comparison to the previous work}

Different methods to construct standard and nonstandard Lagrangians
have been developed and applied to various dynamical systems.  In the
constructed standard Lagrangians, the presence of the kinetic and
potential energy-like terms is evident [1-15,32].  However, the
previously constructed NSLs have different forms, in which such
energy-like terms cannot be easily recognized [12-31]. Some
previously obtained NSLs were of the forms given by either Eq.
(\ref{eq1a}) or Eq. (\ref{eq2a}) [12-15,21,26], but other NSLs
had significantly different forms [17-19,22-25]. It is important
to note that the NSLs obtained using the Jacobi Last Multiplier
method developed by Nucci and Leach (e.g., [17-19]), and those
obtained by El-Nabulsi (e.g., [22-24]) using a different method,
have been applied to many various dynamical systems in applied
mathematics, physics, and astronomy.

In the first attempts to obtain Lagrangians for selected biological
systems, the Lagrangians were found by guessing their forms for
selected population dynamics models [27-29]. Later, those Lagrangians
were formally derived by Nucci and Tamizhmani [30], using the
method of Jacobi Last Multiplier; see also [31] for other applications
of this method to biological systems.  Interesting recent work was
done by Carinena and Fernandez-Nunez [50], who considered systems
of first-order equations and derived the NSLs that are linear, or
more generally affine, in velocities using the method of Jacobi
Last Multiplier. Among the applications of their results to
different dynamical systems, they also included several
population dynamics models [50].

In our previous work [32], we derived the first standard Lagrangians
for five population dynamics models and discussed physical and
biological implications of these Lagrangians for the models.
In this paper, we derived the NSLs by using a method that
significantly modified the one previously developed [12,20].
The results obtained in this paper are now compared to those
previously obtained by Nucci and Tamizhmani [30], as they
are the most relevant to our work. The comparison shows that
there are some advantages of using the method of Jacobi Last
Multiplier, such as the method does not require introducing
forcing functions and gives directly the same Lagrangians as
those found earlier in [27-29]. However, the advantage of
the method developed in this paper is that the derived NSLs
are directly related to their null Lagrangians (NLs) and gauge
functions (GFs).

Having obtained the NSLs for the considered population dynamics
models, we use them to find their corresponding NLs and then
to derive the gauge functions (see Table 2); to the best of
our knowledge, these NLs and GFs are the first obtained for
biological systems, and specifically, for the population
dynamics models.  An interesting result is that the NLs
and their gauge functions are identical for the Lotka-Volterra,
Gompertz, and SIR models, which are caused by the same forms of
NSLs for these three models.  It must also be noted that the
gauge functions are given by the logarithmic functions, which
is a novel form among the previously obtained gauge functions
(e.g., [45-49]). In the approach presented in this paper,
the main differences between the models are shown by the
introduced dissipative functions.

Now, the NSLs, NLs and GFs for the Verhulst and Host-Parasite
models have significantly different forms when compared to the
other three models (see Table 2).  Moreover, there are also
differences in their NSLs and NLs. A new result is the form
of the gauge function $\Phi_1$ for the Host-Parasite model as
this function is given by the exponential function $Ei(B/w_1)$,
which also appeared in the standard Lagrangian found for this
system in [32]; this implies that there exist some relationships
between standard, non-standard and null Lagrangians, and that the
gauge function plays an important role in such relationships;
further exploration of such relationships will be done in a
separate paper.

\section{Conclusions}
A new method to derive non-standard Lagrangians for second-order
nonlinear differential equations with damping is developed and
applied to the Lotka-Volterra, Verhulst, Gompertz, Host-Parasite,
and SIR population dynamics models.  For the considered models,
the method shows some limitations, which are explored and it is
demonstrated that these limitations do not exist for the models
whose equations of motion can be converted into linear ones.

The obtained non-standard Lagrangians are different than those
previously obtained for the same models [30], and the main
difference is that the previously used Jacobi Last Multiplier
method does not require introducing dissipative forcing
functions, which were defined in this paper. However, the
advantage of the method developed in this paper is that it
allows us to use the derived non-standard Lagrangians to
obtain first null Lagrangians and their gauge functions for
the population dynamics models.

By following the recent work [48,49], the presented results
also demonstrate how the derived null Lagrangians and gauge
functions can be used to obtain the equations of motion for
the considered population dynamics models. Our approach to
solving the inverse calculus of variation problem and deriving
non-standard and null Lagrangians is applied to the models of
population dynamics. However, the presented results show that the method can be easily
extended to other biological or physical dynamical systems
whose equations of motion are known.\\

\noindent
{\bf Acknowledgement:} We thank Rupam Das and Lesley Vestal for
discussions of null Lagrangians and gauge functions and their
applications to different dynamical systems.\\