\section{Introduction}
\noindent

Deterministic chaos has represented a crucial issue in laser physics because the Lorenz system [1963] is formally equivalent to laser equations \cite{haken1975}. The first evidence of chaos in lasers  was given in 1982 in a modulated single mode $CO_2$ laser \cite{PRL1982} and we had to wait 1986 for an evidence of Lorenz type chaos in  a particular and little used laser emitting in far infrared region (the so called class C-laser) \cite{PRL1986}. Commonly used single mode lasers are described by two rate equations and they are intrinsic stable devices (class B-laser). Feedback circuits are frequently used to improve  their stability properties by reducing  residual intensity and frequency fluctuations in order to match specific requests. However, a feedback can have the opposite effect, that is, enhancing the relaxation oscillations around the steady state solution. In other words, a simple linear filtering with the appropriate bandwidth on the laser output intensity can induce chaotic fluctuations on it \cite{PRA1986, PRL1987}. From the dynamical point of view, a feedback increases the dimensionality of the phase space from two to three whence chaos becomes possible.
In the last two decades of the past century, the use of $CO_2$ lasers for demonstrating such a behavior has had advantages, mainly due the convenient time scales associated with the laser intensity and the population inversion. The former, evolves on a time scale regulated by the parameter $k$ (the decay rate of photon number of the laser mode) imposed by the length of the optical cavity and its losses. In our case, $k$ is around $10^7$ $s^{-1}$. The latter, i.e. the decay rate of the population inversion, which is called $\gamma$. For a molecular laser as the $CO_2$ laser, $\gamma$ is of the order $10^3-10^4$ $s^{-1}$. These two parameters imply a resulting time scale given by

\[
\sqrt{k \gamma(p_0 - 1)},
\]

where $p_0$ is the pump strength normalized to the threshold value. Usually, $p_0$ is around 2. To be effective in producing chaos, the feedback loop should act on the above mentioned time scale, in other terms, we have the following condition to be satisfied:

\[
k > \beta > \gamma,
\]

where $\beta$ is the bandwidth of the feedback loop.
A simple three dimensional model accounting for the laser intensity $x$ with decay rate $k$, population inversion $y$ with decay rate $\gamma$ and feedback strength $z$ with decay rate $\beta$, explains qualitatively the observed dynamics but it does not yield an accurate matching with the experiment. The problem is overcome by introducing the so called 4-level model for the $CO_2$ laser. It consists in taking into account two resonant levels with population $N_1$ and $N_2$ and two rotational manifolds with population $M_1$ and $M_2$, respectively. In this refined model the dimensionality of the phase space is increased up to 6 instead of the previous 3. For an introduction to the 6 dimensional model and accurate numerical simulations on it, see for example Freire \textit{et al.} [2015]. To reduce the dimensionality two different approaches can be followed. The first one is to use a reduction based on the Center Manifold Theory (CMT) proposed by Varone \textit{et al.} [1995]. This analytical method implies a reduced four dimensional model with the addition of nonlinear terms whose physical interpretation is difficult to provide. A feasible reduction to three dimensions is imaginable considering that it has been obtained for the five dimensional model of the $CO_2$ laser with cavity losses modulation \cite{ciofini}. The second approach for an equivalent three dimensional model which is more physical and straightforward is to take advantage for the 4 level model only from the correct value of the laser intensity in the stationary regime. This condition implies the use of an artificial value of $\gamma$ which can be two orders of magnitude greater, that is, $10^4$ to $10^5$ $s^{-1}$. As the laser output intensity does not depend on $\beta$, we have to use an effective value of $\beta$ for the feedback variable $z$ up to $10^6$ $s^{-1}$. The advantage consists in keeping the original nonlinearity in the $x$ and $y$ differential equations given by their product $x y$. Considering the above advantages and in view of the extension to different dynamical systems ranging from neuron dynamics in the low frequency region (below 1 Hz) to high frequency domains (fast electronics, opto-electronics, etc.) we adopt the following three dimensional model:

\begin{equation}
\label{eq1}
\begin{aligned}
\dfrac{dx}{dt} & = - k_0 x \left( 1 + k_1 z^2 - y \right), \hfill \\
\dfrac{dy}{dt} & = - \gamma y - 2 \frac{k_0}{\alpha} x y + \gamma p_0, \hfill \\
\dfrac{dz}{dt} & = - \beta \left( z - B_0 + \frac{R}{\alpha} x \right), \hfill \\
\end{aligned}
\end{equation}

where $x$ is the fast variable (laser output intensity), $y$ is the slow variable (population inversion) and $z$ is the feedback variable affecting the fast one in a nonlinear way but regulated in linear way as the result of a low pass filter whose input is the fast variable summed to bias.
The parameter $\alpha$ is a suitable normalization of the fast variable $x$ which however does not alter its form considering that it is homogenous in $x$. If $\alpha = 2k_o / \gamma$, the adopted model is formally equivalent to the original model of the laser with feedback.
The paper is organized as it follows. First we introduce a numerical analysis followed by an experimental part containing an analog implementation of the oscillator. Second we compare the new model with other paradigmatic oscillators. Its potentialities are discussed in conclusions.

\section{Minimal Universal Model}

\subsection{Dimensionless form}

We propose the following change of variables and parameters to recast Eqs. (\ref{eq1}) in a dimensionless form. Let's pose:

\[
y \to p_0 y \mbox{, } t \to \frac{t}{\gamma} \mbox{, } \epsilon_1 = \frac{k_0}{\gamma} \mbox{, } \epsilon_2 = \frac{\beta}{\gamma} \mbox{, } B_1 = \frac{R}{\alpha} = \frac{R \gamma}{2 k_0}.
\]

Thus, the minimal universal model for chaos in laser reads:

\begin{equation}
\label{eq2}
\begin{aligned}
\dfrac{dx}{dt} & = - \epsilon_1 x \left( 1 + k_1 z^2 - p_0 y \right), \hfill \\
\dfrac{dy}{dt} & = - y - x y + 1, \hfill \\
\dfrac{dz}{dt} & =  - \epsilon_2 \left( z - B_0 + B_1 x \right), \hfill \\
\end{aligned}
\end{equation}

Let's notice that with the parameter set used in our experiment and analysis, $\epsilon_1 \gg 1$ and $\epsilon_2 \gg 1$. So, model (\ref{eq2}) is a \textit{slow-fast} dynamical systems involving two \textit{fast} times scales. In the following $B_0$ will play the role of a control parameter.

\subsection{Fixed points}

By using the classical nullclines method, it can be shown that the dynamical system (\ref{eq2}) admits four fixed points only two of which are positive.

\begin{equation}
\label{eq3}
\begin{aligned}
I_1 \left(0, 1, B_0 \right), \hspace{2cm} \\
I_2 \left(x^*, y^* = \frac{1}{1+x}, z^* = B_0 - B_1 x \right), \hfill \\
\end{aligned}
\end{equation}

where the expression of $x^*$ (too large to be explicitly written here since it comes from the solution of a cubic polynomial) depends on the control parameter $B_0$. In this problem all fixed points are supposed to be positive. Thus, starting from the right-hand-side of Eqs. (\ref{eq2}), it can be easily shown that:

\[
0 \leqslant x^* \leqslant p_0 - 1, \mbox{ } \frac{1}{p_0} \leqslant y^*,  \mbox{ } B_0 - B_1 \left( p_0 - 1 \right) \leqslant z^* \leqslant B_0.
\]

\subsection{Jacobian matrix}

The Jacobian matrix of dynamical system (\ref{eq2}) reads:

\begin{equation}
\label{eq4}
J = \begin{pmatrix}
(-1 + p_0 y - k_1 z^2) \epsilon_1 \  &  \  p_0 x \epsilon_1 \  &  \  -2 k_1 x z \epsilon_1 \vspace{6pt} \\
-y \  &  \  - 1 -x \  & \  0 \vspace{6pt} \\
- B_1 \epsilon_2 \  &  \  0 \  &  \  \epsilon_2 \vspace{6pt} \\
\end{pmatrix}
\end{equation}

By replacing the coordinate of the fixed points $I_1$ (\ref{eq3}) in the Jacobian matrix (\ref{eq4}) one obtains the Cayley-Hamilton third degree eigenpolynomial which reads:

\begin{equation}
\label{eq5}
\left( \lambda + 1 \right) \left( \lambda + \epsilon_2 \right) \left[ \lambda + \epsilon_1 \left(1 + B_0^2 k_1  - p_0 \right) \right] = 0
\end{equation}

Thus, provided that:

\begin{equation}
\label{eq6}
p_0 < 1 + B_0^2 k_1,
\end{equation}

the fixed point $I_1$ is a \textit{saddle node}. Moreover, such condition (\ref{eq6}) provides a upper boundary for the control parameter $B_0$:

\begin{equation}
\label{eq7}
B_0 < \sqrt{\frac{p_0 - 1}{k_1}}.
\end{equation}

From the positivity of the fixed points, let's notice that if $B_0 = B_1(p_0 - 1)$, then the second fixed point $I_2$ reads: $I_2 (x^* = p_0 -1, y^* = 1/p_0, z^* = 0)$. In these conditions and according to Eq. (\ref{eq7}), it can be stated that:

\begin{equation}
\label{eq8}
B_1 \left( p_0 - 1 \right) < B_0 < \sqrt{\frac{p_0 - 1}{k_1}}.
\end{equation}

Thus, for $B_0 = B_1(p_0 - 1)$, computation of the eigenvalues of $I_1$ shows that two of them are real and negative and one is real and positive confirming thus the \textit{saddle node} feature of this point while for $I_2$, one is real and negative and the two others are complex conjugate with negative real parts. So, in this case the fixed point $I_2$ is stable and attractive according to Lyapunov. For $B_0 = \sqrt{(p_0 - 1)/k_1}$, we found that both fixed points $I_1$ and $I_2$ are stable and attractive (all the real parts of their eigenvalues are negative). Such a result will enable to explain the limits of the bifurcation diagram presented below (see Fig. 1 \& 2) outside which no attractor can exist. Then, while using the parameters set of our experiment, i.e., for any value of $B_0 \in [B_1 ( p_0 - 1 ), \sqrt{(p_0 - 1)/k_1}]$ it can be shown that $I_2$ is a \textit{saddle-focus} (the first eigenvalue is real and negative while the two others are complex conjugate with positive real parts). It implies that a Hopf bifurcation occurs with this interval \cite{Hopf, Andronov, Marsden, Kuznetsov}. Nevertheless, the fact that the fixed points are solutions of a cubic polynomial precludes from explicit the value of the control parameter $B_0$ for which such a Hopf bifurcation occurs. However, we have found a method (see Appendix) allowing to analytically compute the upper bound of such a parameter. In the following we will use the parameters set in our experiment and analysis:

\[
\epsilon_1 = 200 \mbox{ , } \epsilon_2 = 6 \mbox{ , } k_1 = 12 \mbox{ , } p_0 = 1.208  \mbox{ , } B_1 = 0.555.
\]

\subsection{Bifurcation diagram}

Thus, in order to highlight the effects of the control parameter $B_0$ changes on the topology of the attractor we have  built a bifurcation diagram (see Fig. 1 \& 2) that we have compared to the phase portraits plotted in Fig. 3. First, we observe that for $B_0 \approx 0.12$ a Hopf bifurcation occurs (see Appendix). Then, for $B_0 = 0.123$, a \textit{limit cycle} appears. As $B_0$ increases between $0.123$ and $0.1237$, a ``period doubling cascade'' occurs and so, we observe a \textit{2-periodic limit cycle}. In the interval $0.1237 < B_0 < 0.12425$, the period of the \textit{limit cycle} increases again and becomes equal to the number of branches in the bifurcation diagram (see Fig. 2 and Fig. 3a). For $0.12425 < B_0 < 0.129$, a \textit{stable homoclinic orbit} appears and persists (see Fig. 2d and Fig. 3b).

\begin{figure}[htbp]
\centerline{\includegraphics[width=8cm,height=8cm]{diagbifhopf.eps}}
\caption{Bifurcation diagram $x_{max}$ as function of $B_0$.}
\label{fig1}
\end{figure}

\begin{figure}[htbp]
\centerline{\includegraphics[width=8cm,height=8cm]{diagbifzoom.eps}}
\caption{zoom of the bifurcation diagram $x_{max}$ as function of $B_0$.}
\label{fig2}
\end{figure}

\begin{figure}[htbp]
\centerline{\includegraphics[width=6cm,height=6cm]{phase3.eps}}
\centerline{(a) $B_0 = 0.12413$} \vspace{0.1in}
\centerline{\includegraphics[width=6cm,height=6cm]{phase4.eps}}
\centerline{(b) $B_0 = 0.1243$} \vspace{0.1in}
\centerline{\includegraphics[width=6cm,height=6cm]{phase5.eps}}
\centerline{(c) $B_0 = 0.12463$} \vspace{0.1in}
\caption{Phase portraits of model (\ref{eq2}) for various values $B_0$.}
\label{fig3}
\end{figure}

\begin{figure}[htbp]
\centerline{\includegraphics[width=6cm,height=6cm]{timeseries1.eps}}
\centerline{(a) Time series $x(t)$ and $z(t)$ for $B_0 = 0.12413$} \vspace{0.1in}
\centerline{\includegraphics[width=6cm,height=6cm]{timeseries2best.eps}}
\centerline{(b) Time series $x(t)$ for $B_0 = 0.1243$} \vspace{0.1in}
\centerline{\includegraphics[width=6cm,height=6cm]{timeseries3.eps}}
\centerline{(c) Time series $x(t)$ for $B_0 = 0.12463$} \vspace{0.1in}
\caption{Time series of model (\ref{eq2}) for various values $B_0$.}
\label{fig4}
\end{figure}

In order to confirm such scenario, Lyapunov Characteristic Exponents (LCE) have been computed in each case.

\subsection{Numerical computation of the Lyapunov exponents}

The algorithm developed by Marco Sandri [1996] for Mathematica$^{\mbox{\scriptsize{\textregistered}}}$ has been used to perform the numerical calculation of the Lyapunov characteristics exponents (LCE) of dynamical system (\ref{eq2}) in each case. LCEs values have been computed within each considered interval ($B_0 \in [0.123, 0.1234]$ and $[0.1235, 0.125]$). As an example, for $B_0 = 0.123, 0.124$ and $0.1246$, Sandri's algorithm has provided respectively the following LCEs $(0, -0.56, -6.63)$, $(+0.27, 0, -7.47)$ and $(+0.2, 0, -7.32)$. Then, following the works of Klein and Baier [1991], a classification of (autonomous) continuous-time attractors of dynamical system (\ref{eq2}) on the basis of their Lyapunov spectrum is presented in Tab. 1. LCEs values have been also computed with the Lyapunov Exponents Toolbox (LET) developed by Steve Siu for MatLab$^{\mbox{\scriptsize{\textregistered}}}$ and involving the two algorithms proposed by Wolf \textit{et al.} [1985] and Eckmann and Ruelle [1985] (see https://fr.mathworks.com/matlabcentral/fileexchange/233-let). Results obtained by both algorithms are consistent.

\begin{table}[h]
\centering
\caption{Lyapunov characteristics exponents of dynamical system (\ref{eq2}) for various values of $B_0$.}
{\begin{tabular}{c c c}\\[-2pt]
{\hspace{11mm} $m$} & LCE spectrum & Dynamics of the attractor  \\[6pt]
\hline\\[-2pt]
{$0.1230 < B_0 < 0.1234$} & ($ 0, -, - $) & Periodic Motion \\[1pt]
{$0.1235 < B_0 < 0.124$} & ($ 0, -, - $) & $n$-Periodic Motion  \\[2pt]
{$0.1241 < B_0 < 0.125$} & ($ +, 0, - $) & Homolinic Chaos  \\[2pt]
\end{tabular}}
\label{tab1}
\end{table}

\section{Experimental part}

The  set up used in our experiment is illustrated in Fig. 7. Briefly, it consists of three integrators $I_1$, $I_2$ and $I_3$ (LT1114 by Analog Devices), whose outputs are the signals $x$, $y$ and $z$ and $y = \dot{x}$ contained in Eqs. \ref{eq1}. The other integrator $I_4$, is employed for an inverting amplifier with unitary gain. The two nonlinearities are implemented by means of three analog multipliers $M_1$, $M_2$ and $M_3$ (MLT04, by Analog Devices). The first one yields the product $x y$, while the other two multipliers implement the product $x z^2$. The semplification of the proposed scheme is evident when compared with the one obtained using a Field  Programmable Analog  Array (FPAA) circuit\cite{fortunafrasca2005}.

\begin{figure}[htbp]
\centerline{\includegraphics[width=6cm,height=6cm]{Fig5aBEST.eps}}
\vspace{0.1in}
\centerline{(a) Phase portrait for $V(B_0) = -2.550V$}
\vspace{0.1in}
\centerline{\includegraphics[width=6cm,height=6cm]{Fig5bBEST.eps}}
\vspace{0.1in}
\centerline{(b) Phase portrait for $V(B_0) = -2.596V$}
\vspace{0.1in}
\centerline{\includegraphics[width=6cm,height=6cm]{Fig5cBEST.eps}}
\vspace{0.1in}
\centerline{(c) Phase portrait for $V(B_0) = -2.628V$}
\vspace{0.1in}
\caption{Oscilloscope snapshots of phase portraits for $V(p_0) = -2.983 V$ and various values of $V(B_0)$.}
\label{fig5}
\end{figure}

\begin{figure}[htbp]
\centerline{\includegraphics[width=6cm,height=6cm]{Fig6a.eps}}
\vspace{0.1in}
\centerline{(a) Time series $x(t)$ and $z(t)$ for $V(B_0) = -2.550V$}
\vspace{0.1in}
\centerline{\includegraphics[width=6cm,height=6cm]{Fig6b.eps}}
\vspace{0.1in}
\centerline{(b) Time series $x(t)$ for $V(B_0) = -2.596V$}
\vspace{0.1in}
\centerline{\includegraphics[width=6cm,height=6cm]{Fig6c.eps}}
\vspace{0.1in}
\centerline{(c) Times series $x(t)$ for $V(B_0) = -2.628V$}
\vspace{0.1in}
\caption{Oscilloscope snapshots of time series for $V(p_0) = -2.983 V$ and various values of $V(B_0)$.}
\label{fig6}
\end{figure}

\begin{figure}[htbp]
\centerline{\includegraphics[width=9 cm,height=11cm]{circuitdiag.eps}}
\caption{Circuit diagram of the ``laser with feedback''.}
\label{fig7}
\end{figure}

Considering the limits imposed by analog simulations, the  desired dynamics of the oscillator is obtained by fine adjustments of  two bias voltages $V(p_0)$ and $V(B_0)$ accounting for the  parameters, $p_0$ and $B_0$  in Eqs. \ref{eq1}. This condition  has been achieved by means of two potentiometers $P_1$ and $P_2$ connected to a fixed negative voltage source $-Vs$. The relaxation rates of the three variables (reciprocal of the integration times of the three integrators), according to Eqs. \ref{eq1}  and after a temporal rescaling of three orders of magnitude, are selected as follows. For the $x$ integrator we set  $k_0 = 1/(R_1 C_1)= 2.463 \times 10^4s^{-1}$ where $R_1 = 33 k \Omega$ and $C_1 = 1.67nF$, ($R_3 = R_1 = 33 k \Omega$). For the $y$  integrator we set  $\gamma = 1/ (R_4 \times C_2) = 1.00 \times 10^2s^{-1}$, where $R_4 = 10 k \Omega$, $C_2 = 1 \mu F$. For the $z$ integrator we select  $\beta = 1/ (R_7 \times C_3 ) = 1.00 \times 10^3s^{-1}$, with $R_7 = 10 k \Omega$ and $C_3 = 100nF$. For the other parameter values we have:$ k_1 = R_1/(R_2 \times 2.5) = 13.2$ where $R_2 = 1 k \Omega$; $p_0 = R_4/R_6 = 1.196$ where $R_6=8.36 k \Omega$;  $B_0 = R_7/R_9 = 0.115$ where $R_9 = 86.6 k \Omega$; $R/\alpha = R_7/R_8 = 0.222$, where $R_8 = 45 k \Omega$.

The attractors in the $x-z$ phase space for different values of the $B_0$ parameter value are reported in Figs. 5.  In panel a) the local dynamics emerging after an Hopf bifurcation is shown. An increase of the control parameter $B_0$ leads to a condition of homoclinic chaos due to a homoclinic orbit around which a chaotic regime characterized by pulses nearly of the same high but erratically separated in time due to the rejection mechanism around the local chaos shown in a). A successive increment of the control parameter leads to the stabilization of the pulsed regime which are typical in a relaxation oscillator. In the presented attractors (see b and c) we observe a small distortion for high amplitudes of the $x$ signal due to the frequency limitation of the analog multipliers. This limitation does not induce severe limitations to the global dynamics as shown in the temporal behavior of the $x$ signal, when we consider applications in the low frequency range below 1 kHz. For applications in high frequency regimes the nonlinearities must be implemented by using CMOS devices (range 1-10 MHz).

\section{Discussion}

Numerical analysis of the proposed model and its electronic implementation confirm its potentialities that promote it to the class of  as other well-known paradigmatic models as the Lorenz's \cite{lorenz1963}, Chua's \cite{matsumoto1984,chua1986} Chen's \cite{guanrong2010,chen2002},  Roessler's \cite{roessler1976} and other simple three dimensional systems reported by Sprott \cite{sprott}. The Lorenz model is formally equivalent to the laser equations (the so called class C-laser), as demonstrated by H. Haken \cite{haken1975} so the link with laser dynamics is direct. The intrinsic symmetries imply a chaotic trajectory visiting two saddle foci to be contrasted with the laser with feedback where the competing steady states are a saddle focus and a saddle node which contributes to the creation of the homoclinic connection. The comparison with the Chua's circuit which relies on a ``locally active resistor'' with static nonlinear characteristic is of another kind. This element is the Chua's diode and it can be implemented in different ways. The chaotic attractors from the Chua circuit perfectly reproduce Lorenz chaos even though the correspondence with laser equations is difficult to draw. The chaotic Chen system is similar but not equivalent to the Lorenz one as recently pointed out by Chen \cite{chen2020} introducing of the concept ``generalized Lorenz systems family''. In both cases the nonlinearities are two of the quadratic type in models which have seven terms on the right hand side. The Roessler circuit is simpler when compared with Lorenz systems because it possesses only one quadratic term and seven terms on the right hand side. It does not seem that both Chua and Roessler systems are related to laser dynamics. The above cited systems played a crucial role in chaos synchronization demonstrated by Pecora  and Carroll \cite{pecoracarroll1990} using Lorenz and Roessler circuits.
The increased complexity in the laser with feedback is due to feedback process which implies the additional cubic term $xz^2$ to the two quadratic nonlinearities $xy$. This is the first time that the cubic term is treated in its simplest form approximating $\sin(z)^2$ with $z^2$. It is important to stress that adiabatic elimination of the fast variable corresponding to laser polarization in the Lorenz system does not lead to rate model of the laser equations the so called class B laser, the large class of the available lasers, including semiconductor lasers and solid state lasers. As far as laser semiconductor lasers are concerned, it is important to note that their dynamics is well described by the Lang and Kobayashi  (LK) model \cite{Lang} accounting for the effects of delayed optical feedback acting on the timescale of the intrinsic semiconductor laser. The LK equations  describe the complex  dynamics of the complex electric field $E$ and the inversion (number of electron-hole pairs) $N$ inside the laser. The fast chaotic dynamics from these laser has been largely used in applications  for secure communication systems (\cite{Fischer, roy, donati, Ohtsubo}).

\section{Conclusions}

In conclusion, we retain that the laser with feedback, whose historical origins dates back in the same period of other chaotic oscillators will be appropriate for describing instances of local chaos, reached after subharmonic bifurcations, global bifurcations as homoclinic chaos and relaxation type oscillation behavior or regular spiking behavior when a control parameter is changed. Another valuable advantage over the nonlinear oscillators described above is related to the fact that  elimination of the feedback variable leaves unchanged its potentialities and chaos can be reached by modulation of the parameter $k$ provided the general condition on the timescales fulfilled. In this framework, interesting perspectives are also foreseen for competition population dynamics ruled by  Volterra-Lotka models \cite{volterra1926, volterra1931, lotka1910, lotka1920} and  chaotic epidemiological models \cite{SEIR}, in both cases the underlying nonlinearities are of the direct product of two variables.

\section{Appendix}

This appendix presents a method allowing to provide a upper bound for the Hopf bifurcation parameter for any three-dimensional autonomous dynamical system for which the fixed points coordinates cannot be easily expressed analytically as it is the case for the dynamical system (\ref{eq2}) for which the coordinates of the fixed point $I_2$ are the roots of a cubic polynomial. Let's suppose that the three eigenvalues of the Jacobian matrix $J$ of this dynamical system evaluated at the fixed point ($I_2$ in our case) are real and complex conjugate $\lambda_1$, $\lambda_{2, 3} = \alpha \pm i \omega$. The Cayley-Hamilton eigenpolynomial reads:

\begin{equation}
\label{eq9}
\lambda^3 -\sigma_1 \lambda^2 + \sigma_2 \lambda - \sigma_3 = 0
\end{equation}

where $\sigma_1 = Tr \left( J \right)$, $\sigma_2 = \sum_{i = 1}^3 M_{ii} \left( J \right)$ is the sum of all first-order diagonal minors of $J$ and $\sigma_3 = Det \left( J \right)$. Thus, we have:

\begin{equation}
\label{eq10}
\begin{aligned}
\sigma_1 & = Tr \left( J \right) = \lambda_1 + \lambda_2 + \lambda_3 = \lambda_1 + 2 \alpha, \hfill \\
\sigma_2 & = \sum_{i = 1}^3 M_{ii}  \left( J \right) = \lambda_1 \lambda_2 + \lambda_1 \lambda_3 + \lambda_2 \lambda_3 = 2 \alpha \lambda_1 + \beta, \hfill \\
\sigma_3 & = Det \left( J \right) = \lambda_1 \lambda_2 \lambda_3 = \lambda_1 \beta, \hfill \\
\end{aligned}
\end{equation}

where $\beta = \alpha^2 + \omega^2$. In order to analyze the stability of fixed points according to a control parameter value ($B_0$ here), i.e. the occurrence of Hopf bifurcation, we propose to use the Routh-Hurwitz' theorem \cite{Routh1877,Hurwitz1893} which states that if $D_1 = \sigma_2$ and $D_2 = \sigma_3 - \sigma_2 \sigma_1$ are both positive then \textit{eigenpolynomial} equation (\ref{eq9}) would have eigenvalues with negative real parts. From Eqs. (\ref{eq10}) it can be stated that:

\begin{equation}
\label{eq11}
\alpha = \frac{\sigma_1 \sigma_2 - \sigma_3}{\lambda_1^2 + \sigma_2}
\end{equation}

Thus, $\alpha = 0$ provided that $\kappa = \sigma_1 \sigma_2 - \sigma_3 = 0$. For dynamical system (\ref{eq2}), we obtain:

\begin{equation}
\label{eq12}
\kappa = \left( 1 + x \right) \left[ p_0 x y \epsilon_1 + \left( 1 + x + \epsilon_2 \right) \epsilon_2 \right] -2 B_1 k_1 x z \epsilon_1 \epsilon_2^2
\end{equation}

then, by replacing $x$, $y$ and $z$ the coordinates (\ref{eq3}) of the fixed point $I_2$, we have:

\begin{equation}
\label{eq13}
B_0 = \frac{ 1 + \epsilon_2}{2 B_1 k_1 x^* \epsilon_1 \epsilon_2} + \frac{p_0 \epsilon_1 + \epsilon_2  (2 + \epsilon_2)}{2 B_1 k_1 \epsilon_1 \epsilon_2^2} + \frac{(1 + 2 B_1^2 k_1 \epsilon_1 \epsilon_2)x^*}{2 B_1 k_1 \epsilon_1 \epsilon_2}
\end{equation}

Positivity of fixed points has led to $x^* \leqslant p_0 -1$ which implies that $\max\left( x^* \right) = p_0 - 1$. Thus, by posing $x^* = p_0 - 1$ in (\ref{eq13}) and while using the parameters sets of our experiment, i.e., $\epsilon_1 = 200$, $\epsilon_2 = 6$, $k_1 = 12$, $p_0 = 1.208$ and $B_1 = 0.555$, we find:

\[
B_0^{Hopf} \leqslant 0.12057
\]

The numerical computation of the Hopf bifurcation parameter value has been found equal to $0.12036$ which is below and very near the upper bound analytically obtained.


\section{INTRODUCTION}
The fission process emits an enormous amount of energy per event, of the order of 200 MeV per nucleus for actinides. The majority of the energy released in the fission of actinide nuclei appears in the form of the kinetic energy of the two fission fragments. According to the conservation of momentum, the partitioning of energy between the two fragments is dependent on their masses; the lighter fragment garners more kinetic energy than its heavier counterpart. The sum of the two fragments’ kinetic energy is the total kinetic energy (TKE). The TKE is an important fundamental feature of fission, and studies of TKE can help to elucidate some of the physics behind the large scale collective motion of the nucleus. \newline
\indent	Two dynamics dictate the magnitude of the TKE release. The first is the large Coulomb repulsion between the proto-fragments at the scission point,  and the second is the conversion of the collective energy of the nucleus to kinetic energy of the fragments as the nucleus transitions from the saddle point to the scission point. Based on these features, we can write TKE as a function of the fragments charge and radius, assuming R = $r_o\;^{1/3}$.

\let\saveeqnno\theequation
\let\savefrac\frac
\def\displaystyle\savefrac{\displaystyle\savefrac}
\begin{eqnarray}
\let\frac\displaystyle\savefrac
\gdef\theequation{1}
\let\theHequation\theequation
\label{dfg-b55e055b7286}
\begin{array}{@{}l}\begin{array}{l}TKE\;=\frac{Z_1Z_2e^{2}}{r_0(A_1^{1/3}+A_2^{1/3})}\\\end{array}\end{array}
\end{eqnarray}
\global\let\theequation\saveeqnno
\addtocounter{equation}{-1}\ignorespaces
Equation 1 is derived from the Coulomb forces acting on the fragments. While this equation is reasonable at making first order approximations of TKE, it slightly overestimates the impact of Coulomb repulsion because the equation assumes the proto-fragments are perfect spheres at the scission point. In reality, they are highly deformed. Such deformation means the proto-fragments are more distant from each other at the scission point, and therefore experience less Coulomb repulsion. A more widely used equation is
\let\saveeqnno\theequation
\let\savefrac\frac
\def\displaystyle\savefrac{\displaystyle\savefrac}
\begin{eqnarray}
\let\frac\displaystyle\savefrac
\gdef\theequation{2}
\let\theHequation\theequation
\label{dfg-1a22b593dbe7}
\begin{array}{@{}l}TKE\;=\;\frac{(0.1189\;\pm\;0.0011)Z^{2}}{A^{1/3}}+7.3\;\pm\;1.5\;MeV\end{array}
\end{eqnarray}
\global\let\theequation\saveeqnno
\addtocounter{equation}{-1}\ignorespaces

derived by Viola, et al. \cite{1} where the second term represents the conversion of energy stored in deformation to kinetic energy.

\subsection{Previous studies of of $^{237} Np (n,f)$}
While the dependency of TKE on thermal neutron-induced fission has been extensively studied in many actinide systems, the same cannot be said with high incident neutron energies. There have been a limited number of studies on the $^{237}$ Np (n,f)\ reaction \cite{2,3} in which TKE release and mass distributions were reported. These studies are summarized in Table 1. Until this work, there have been no measurements made on the $^{237}$ Np (n,f)reaction for En > 5.5 MeV. As such the measurements presented in this paper substantially extend our knowledge of the TKE release for the $^{237} $Np (n,f) reaction.

\begin{table}[!htbp]
\caption{{Previous work on the 237Np (n,f) reaction.} }
\label{table-wrap-93af179abb8747a0978484aea4b9e73d}
\begin{ruledtabular}
\ignorespaces
\centering
\begin{tabular}{ll}

\textbf{E}\ensuremath{_{\textbf{n}}} \textbf{ range (MeV)} &
  \\
 0.3 - 5.5 &
  \unskip~\cite{3}\\
 0.8 and 5.5 &
  \unskip~\cite{2}\\

\end{tabular}\par
\end{ruledtabular}
\end{table}
\section{EXPERIMENTAL}
This experiment was carried out using the 15R beamline at Los Alamos National Laboratory-Los Alamos Neutron Science Center at the Weapons Neutron Research (LANL-LANSCE-WNR) facility over a 7 day period in December 2019. The LANSCE-WNR provides a high fluence of “white spectrum” neutrons generated by spallation of a thick tungsten target by 800 MeV protons, producing ~ $10^5 – 10^6$ n/s up to En = 100 MeV \cite{4, 5}. The proton beam consists of  625 $\mu$s macropulses. Each macropulse is composed of 250 ps-wide micropulses that are spaced 1.8 $\mu$s apart. The macropulses had a repetition rate of 100 Hz. The pulsed nature of the proton beam enables the measurement of the time-of-flight (energy) of the incident neutrons. \newline
\indent The general arrangement of the experiment is shown in Figure 1. The neutron flight path i.e. the distance from the spallation source to the $^{237}$Np target, was measured to be 13.85 m.  The neutron beam was collimated to a 1 cm diameter. An evacuated, thin-walled aluminum scattering chamber held the $^{237}$Np target and the Si PIN diode fission detectors. At the end of the 7 day irradiation period, a total of 8,618 coincident fission events were detected.

\begin{figure*}
 \includegraphics[angle=-90,width=.75 \textwidth, trim=6cm 1.5cm 3cm 1.3cm, clip=true]{Figure1}
 \centering
\vspace{-	2 cm}
\caption{Schematic diagram of the detection set-up at the LANSCE facility}
\label{fig:1}
\end{figure*}

\subsection{Fission Detectors}
The scattering chamber was aligned with the neutron beam using a combination of neutron-exposed photo plates and an alignment laser. Inside the scattering chamber, two arrays of nine Si PIN photodiodes (Hamamatsu S3590-09), arranged in a 3x3 array, were used to detect the fission products. The arrays were arranged on opposite sides of the beam at angles of $65^o$ and $115^o$ from the incident beam. Each individual PIN diode has an area of 1 $cm^2$, and the detectors are placed 3.5 cm from the actinide target. Coincident fission fragments in both the left and right array were required to process a signal. The pulse height defect of the detectors for fission fragments was determined by employing the Schmitt method \cite{6} with a ${^252}$Cf source. The fragment energy loss in the actinide deposition and the carbon foil target backing was determined by Northcliff- Schilling energy loss tables on an event by event basis \cite{7}.\newline
\indent The time of flight (TOF) of each interacting neutron was measured using a timing pulse from a Si PIN diode and the accelerator RF signal. The neutron energies were calculated based on the relativistic relationship between the distance traveled by the neutrons, the position of the photo fission peak in the fission time of flight spectrum, and the observed time difference between the neutron timing signal and the accelerator RF signal \cite{8}. The position of the photofission peak was measured to have an uncertainty of 0.24 ns, with a relative uncertainty of 0.52\%. The 1.8 $\mu$s proton pulse rate allows for a small fraction of “wrap around” neutrons to reach the scattering chamber among the neutrons of the previous pulse. Accordingly, these neutrons would either be of very high energy, or very low energy. For the 13.85 m flight-path used in this experiment, the low-energy “wrap around” neutrons would have En $<$ 0.4 MeV. The near-negligible cross section for the $^{237}$Np(n, f) reaction at this energy means that these low-energy “wrap around” neutrons were ignored. Similarly, the high-energy “wrap around” neutrons would be of sufficient energy as to appear with the slowest neutrons from the previous beam burst. Wrap around events of this nature set the lower limit for the neutron energy resolution. Since no events below the fission threshold of En $\leq$ 0.7 MeV \cite{9} were detected, the high energy “wrap around” neutrons were also ignored. The neutron energies were binned logarithmically to give bins of equal associated uncertainty for the neutron energy.
\subsection{The 2E Method}
The 2E method is used to measure the fission product masses and their TKE by detecting the fission products and determining their energies after neutron emission. (We use the convention that the term “fission fragments” refers to pre-neutron emission while the term “fission product refers to post neutron emission.) Conservation of momentum and mass for the coincident fission products gives:
let\saveeqnno\theequation
\let\savefrac\frac
\def\displaystyle\savefrac{\displaystyle\savefrac}
\begin{eqnarray}
\let\frac\displaystyle\savefrac
\gdef\theequation{3}
\let\theHequation\theequation
\label{dfg-afdea13e84b0}
\begin{array}{@{}l}M_H^\ast V_H^\ast\;=\;M_L^\ast V_L^\ast\end{array}
\end{eqnarray}
\global\let\theequation\saveeqnno
\addtocounter{equation}{-1}\ignorespaces

\let\saveeqnno\theequation
\let\savefrac\frac
\def\displaystyle\savefrac{\displaystyle\savefrac}
\begin{eqnarray}
\let\frac\displaystyle\savefrac
\gdef\theequation{4}
\let\theHequation\theequation
\label{dfg-ba6bb0ee3049}
\begin{array}{@{}l}M_L^\ast\;+\;M_H^\ast\;=\;M_{CN}\end{array}
\end{eqnarray}
\global\let\theequation\saveeqnno
\addtocounter{equation}{-1}\ignorespaces

\let\saveeqnno\theequation
\let\savefrac\frac
\def\displaystyle\savefrac{\displaystyle\savefrac}
\begin{eqnarray}
\let\frac\displaystyle\savefrac
\gdef\theequation{5}
\let\theHequation\theequation
\label{dfg-e48bfbe4ac93}
\begin{array}{@{}l}M_H^\ast E_H^\ast\;=\;M_L^\ast E_L^\ast\end{array}
\end{eqnarray}
\global\let\theequation\saveeqnno
\addtocounter{equation}{-1}\ignorespaces

\let\saveeqnno\theequation
\let\savefrac\frac
\def\displaystyle\savefrac{\displaystyle\savefrac}
\begin{eqnarray}
\let\frac\displaystyle\savefrac
\gdef\theequation{6}
\let\theHequation\theequation
\label{dfg-8bb973e37bcf}
\begin{array}{@{}l}M_L^\ast\;=\;M_{CN}\;\frac{E_H^\ast}{E_H^\ast+E_L^\ast}\end{array}
\end{eqnarray}
\global\let\theequation\saveeqnno
\addtocounter{equation}{-1}\ignorespaces
where H represents the heavy fragment, L represents the light fragment, CN represents the compound nucleus, and ‘indicates a quantity that pertains to pre-neutron emission. Our detection apparatus, due to the sub-femtosecond timescale of the compound nucleus’ descent to scission, measures post neutron evaporation fission products \cite{10}. Determination of mass yields requires pre-neutron evaporation masses. As such, we must have a sensible approximation for prompt neutron multiplicity as a function of excitation energy and post neutron fragment mass. We therefore make the following corrections for neutron emission.

\let\saveeqnno\theequation
\let\savefrac\frac
\def\displaystyle\savefrac{\displaystyle\savefrac}
\begin{eqnarray}
\let\frac\displaystyle\savefrac
\gdef\theequation{7}
\let\theHequation\theequation
\label{disp-formula-group-60bdec3559a94c37b6252425a37d1dad}
\begin{array}{@{}l}E_{L,H}^\ast\;=\;E_{L,H}\;\left(1+\frac{\nu_{post}(m_{L,H}^\ast)}{m_{L,H}}\right)\;\end{array}
\end{eqnarray}
\global\let\theequation\saveeqnno
\addtocounter{equation}{-1}\ignorespaces

\let\saveeqnno\theequation
\let\savefrac\frac
\def\displaystyle\savefrac{\displaystyle\savefrac}
\begin{eqnarray}
\let\frac\displaystyle\savefrac
\gdef\theequation{8}
\let\theHequation\theequation
\label{disp-formula-group-065cf2a2eb3a4b5abedc5c2afbd3a7bb}
\begin{array}{@{}l}m_{L,H\;}=\;m_{L,H}^\ast\;-\;\nu_{post\;L,H}\;\left(E_{n,}m_{L,H}^\ast\right)\end{array}
\end{eqnarray}
\global\let\theequation\saveeqnno
\addtocounter{equation}{-1}\ignorespaces

where ${v_{post} (A)}$ is the post fission prompt neutron multiplicity calculated from the General Description of Fission Observables (GEF) model \cite{9, 11}. In the laboratory system post neutron emission fission product energies $(E_{L,H})$ are calculated by correcting the measured fission product energies for detector pulse height defects and fission product energy losses in the target material and backing \cite{8}. To deduce the post neutron emission mass distributions an iterative correction procedure \cite{12} is applied. This iterative procedure assumes zero charged particle emission and ceases when the differences in fragment mass between iterations is less than 0.1 amu.

\subsection{Corrections for PHD and Fragment energy loss}
We measure, on an event by event basis, the pulse heights of the coincident fission fragments. To accurately transform these pulse heights into energies, we must account for the pulse height defect (PHD). A PHD occurs when a fission fragment strikes the PIN diode with such high energy that recombination occurs in the tracks of the incident particle along with any losses of energy in the detector dead layer. To make corrections for this effect we apply the Schmitt method \cite{6}. The Schmitt method exploits the known fission spectra of $^{252}$Cf to determine a set of constants the describe the PHD for a given detector. For Si surface barrier detectors, the magnitude of the PHD generally ranges from 3-13 MeV \cite{13,14,15} The energy loss of the fission products in the target deposit and the carbon foil were calculated using the Northcliffe-Schilling tables \cite{7}.

\begin{table*}[!htbp]
\caption{{TKE and $_{post}$TKE in the $^{237}$Np (n,f) reaction as a function of incident neutron energy. The neutron bin limits are given in the first column while the second column is the geometric mean of the neutron energies. The last column is the number of events N in the bin.} }
\label{table-wrap-67ef239583c14a9fa1a524d36be95d2f}
\begin{ruledtabular}
\ignorespaces
\centering
\begin{tabular}{llllll}

\textbf{E}\ensuremath{_{\textbf{n}}} \textbf{ range (MeV)} &
  \textbf{E}\ensuremath{_{\textbf{n}}} \textbf{(MeV)} &
  \textbf{\ensuremath{_{post}}} \textbf{TKE(MeV)} &
  \textbf{\ensuremath{_{pre}}}\textbf{TKE} \textbf{(MeV)} &
  \textbf{s}\textbf{\ensuremath{^{2}}}\textbf{\ensuremath{_{post}}}\textbf{(MeV}\textbf{\ensuremath{^{2}}} \textbf{)} &
  \textbf{Events}\\
 2.60-3.31 &
   2.96 &
   170.8 \ensuremath{\pm} 0.4 &
   173.1 \ensuremath{\pm} 0.4 &
   92.78 \ensuremath{\pm} 0.33 &
   887\\
 3.32-4.38 &
   3.83 &
   171.6 \ensuremath{\pm} 0.4 &
   173.8 \ensuremath{\pm} 0.4 &
   91.36 \ensuremath{\pm} 0.33 &
   868\\
 4.39-6.09 &
   5.19 &
   171.0 \ensuremath{\pm} 0.4 &
   173.6 \ensuremath{\pm} 0.4 &
   89.40 \ensuremath{\pm} 0.33 &
   865\\
 6.10-8.06 &
   7.06 &
   171.2 \ensuremath{\pm} 0.4 &
   173.3 \ensuremath{\pm} 0.4 &
   109.83 \ensuremath{\pm} 0.35 &
   861\\
 8.07-11.66 &
   9.7 &
   170.2 \ensuremath{\pm} 0.4 &
   172.2 \ensuremath{\pm} 0.4 &
   107.95 \ensuremath{\pm} 0.34 &
   867\\
 11.67-18.56 &
   14.8 &
   168.0 \ensuremath{\pm} 0.4 &
   170.9 \ensuremath{\pm} 0.4 &
   109.62 \ensuremath{\pm} 0.42 &
   861\\
 18.57-29.66 &
   23.9 &
   167.1 \ensuremath{\pm} 0.4 &
   169.5 \ensuremath{\pm} 0.4 &
   107.33 \ensuremath{\pm} 0.35 &
   863\\
 29.67-45.72 &
   37.1 &
   164.9 \ensuremath{\pm} 0.4 &
   167.5 \ensuremath{\pm} 0.4 &
   118.85 \ensuremath{\pm} 0.39 &
   862\\
 45.73-68.78 &
   56.5 &
   166.0 \ensuremath{\pm} 0.4 &
   168.6 \ensuremath{\pm} 0.4 &
   113.85 \ensuremath{\pm} 0.36 &
   870\\
 68.79-100 &
   83.8 &
   164.8 \ensuremath{\pm} 0.5 &
   168.0 \ensuremath{\pm} 0.4 &
   126.34 \ensuremath{\pm} 0.44 &
   814\\

\end{tabular}\par
\end{ruledtabular}
\end{table*}

\begin{figure}
 \includegraphics[width=.5\textwidth]{Figure2}
\caption{Plots for the post TKE distributions sorted into the energy bins described in Table 2. The plotted curves (line) represent the results of fitting the distributions with Gaussian functions. (Color online).}
\label{fig:2}
\end{figure}

\section{Results}
In Table 2, we report the pre- and post-neutron TKE and its variance as a function of incident neutron energy for the $^{237}$Np (n,f) reaction. This data set significantly increases the range of energies studied compared to previous measurements that only explored this reaction to a maximum value of 5.5 MeV in neutron energy \cite{2,3}. In  Figure 2, we present the post neutron evaporation TKE distributions sorted by neutron energy bin. The data were binned as to have a similar number of events per neutron energy bin. The thin red lines in each plot represent the results of fitting the data with Gaussian distributions.

\subsection{TKE vs $E_n$}
 In Figure 3, we show the $_{post}$TKE response up to En = 6 MeV from this work and that documented in the literature \cite{2, 3}. In doing so, we restrict our attention to systems where only first chance fission is expected. The experiment found the TKE released in the low En energy range to be ~1.7 and 0.9 MeV higher than that reported by Hambsch \cite{3} and Naqvi \cite{2} respectively. It is important to note that Hambsch normalized their data to a $^{235}$U pre-neutron TKE of 179.5 MeV, a value ~0.5 MeV lower than the compilation of Madland \cite{16}. Similarly, Naqvi compared their $^{237}$Np TKE data to a $^{235}$U post-neutron TKE value of 167.45 MeV. This value is 1.7 MeV lower than what is reported by Madland. As a result, both the Naqvi and Hambsch data likely underestimate the TKE release in $^{237}$Np. \newline
 \indent The $_{post}$TKE release for the fast neutron-induced fission of $^{237}$Np over the entire measured energy range is given in Figure 4. In Figure 4, we show our data, as well as comparisons to other $^{237}$Np (n,f) reaction data and the GEF model. The overall trend of decreasing TKE with increasing excitation is consistent with other fast neutron TKE studies of actinides \cite{17,18,19,20,21,22,23,24,25,26}. While the TKE release for the lowest energy neutrons is approximately linear, a first order log$_{10}(En)$ polynomial fit is needed to describe the fast neutron energy region (En$ >$ 1 MeV).  The fit to our data (Figure 4) is given by the equation TKE(MeV) = (174.38$\pm$  0.72) – (5.11 $\pm$ 0..5821) log$_{10}$ En. \newline
 \indent An interesting feature apparent in Figure 4 is a “bump” in TKE in the En   = 46-69 MeV energy bin. This peculiar bump has been observed in other actinides within this neutron energy range, particularly with $^{232}$Th, $^{235}$U and $^{239}$Pu \cite{17, 18, 25}. This anomaly is not predicted in models. The observation of such a “bump” across several actinides suggests either a systematic flaw in our measurement or a gap in our understanding of nuclear physics in this energy range.

 \subsection{$_{post}$TKE comparison to the GEF model}
When comparing the GEF model postTKE values to those measured in this work, GEF systematically overestimates postTKE by ~1.2$\%$ for all data points in the region En $\leq$ 50 MeV. This is equivalent to approximately 2 MeV per energy bin (Figure 4). This overestimation is likely the result of the GEF model misattributing a small fraction
of the energy of the fissioning $^{237}$Np nucleus to the kinetic energy of the fission fragments instead of particle evaporation. A study of the high energy fission of $^{235}$U(n,f) found that the  GEF framework routinely underestimated the total neutron evaporation by 5-10$\%$  for En = 0- 50 MeV \cite{27}. A similar underestimation is seen for the $^{239}$Pu (n,f) reaction, where the GEF prompt neutron multiplicity underestimates the known values for A $>$ 147 amu, and/or En $>$ 12 MeV \cite{28,29,30,31,32}. \newline
\indent The GEF model underestimates the $_{post}$TKE for En $>$ 50 MeV by ~ 0.7$\%$ or 1.1 MeV. This transition from the GEF model overestimation to sudden underestimation is due to the “bump” in the experimental data and further highlights the importance of being able to understand this anomaly.

\begin{figure}
 \includegraphics[width=.5\textwidth]{Figure3}
\caption{Comparison of the average postTKE measured for the 237Np (n,f) reaction as a function of incident neutron energy measured in this work to previous studies and a fit to our data (dashed line)(Color-online)}
\label{fig:3}
\end{figure}

\begin{figure}
 \includegraphics[angle=-90,width=.5\textwidth]{Figure4}
 \vspace{-1cm}
\caption{Post neutron TKE released for the 237Np (n,f) reaction from En =  2.6 -100 MeV. TKE error bars are statistical (Color-online)}
\label{fig:4}
\end{figure}

\subsection{TKE distribution Variances ($\sigma$)}
In Figure 2, the TKE data are fit with Gaussian distributions. The variances of the Gaussian fits, given in Table 2, are sensitive markers of multi-chance fission. At energies where multi-chance fission is possible, a new fission pathway opens, and the presence of this new channel is reflected in the broadening of the Gaussian distribution. This phenomenon is well documented, particularly in $^{235}$U (n,f) studies \cite{17,33}. In Figure 5, we limit our attention to measured variances below 20 MeV to highlight the dependence of the measured variances of the TKE distributions on the number of fission channels present. As expected, we see a significant increase in the measured variance at the onset of second chance fission, ~7 MeV.

\begin{figure}
 \includegraphics[width=.5\textwidth, trim=2cm 1.5cm 1cm 1.3cm, clip=true]{Figure5}
\vspace{-2cm}
\caption{Top: cross section for $^{237}$Np (n,f) from ENDF/B-VIII, ENDF/B-VI. Bottom: measured $TKE  \sigma^2$ vs En for this study. (Color online)}
\label{fig:5}
\end{figure}

\subsection{Mass Yield Distributions}
In actinide fission, the fission fragment mass distribution is bimodal at low excitation energies. As E$_n$ is increased the valley between the asymmetric fission peaks begins to fill. In Figure 6 we plot the mass distributions measured in this work, normalized to 200$\%$,  and the predictions of the GEF model. The resolution of our measurement is $\pm$ 5 amu. It is evident in  Figure 6 that as End increases the two asymmetric peaks merge reflecting the emergence of symmetric fission. The rate of symmetric ingrowth, and its relative proportion to the total number of events, is in excellent agreement with that predicted by GEF \cite{9}. \newline
\indent It is well documented in the literature \cite{34} that the constancy of the heavy fragment mass peak is due to the influence of the spherical N=82 and deformed N=88 neutron shells and Z = 50 proton shell. In this system the mass of the heavy fragment is anchored at ~ 139.4 amu. This is in agreement with \cite{3}.

\begin{figure}
 \includegraphics[angle = 90, width=.55\textwidth]{Figure6}
\caption{Comparison of the mass distributions measured in this work (red squares), normalized to 200$\%$,  and the predictions of the GEF$_{post}$  model (blue line)\cite{9}  (Color-online)}
\label{fig:6}
\end{figure}

 \subsection{Fission Channel Symmetry and TKE}
 It has been demonstrated in many actinide systems that the general trend of decreasing TKE with increasing excitation is a consequence of the increasing probability of symmetric fission as En increases. For low energy neutron-induced fission, asymmetric fission dominates due to nuclear shell effects. As En increases, symmetric fission begins to emerge, and eventually, symmetric fission dominates. Figure 7 shows the fragment mass and TKE distributions for $^{237}$Np at various energies. At the lowest neutron energies, the fissioning system is almost purely asymmetric. As En increases, more fission fragments with A = 120 amu are detected indicating a symmetric fission event.

 \begin{figure}
\includegraphics[angle=-90,width=.5 \textwidth]{Figure7}
\caption{$^{237}$Np (n,f) product  mass TKE distributions at various En.  (Color online)}
\label{fig:7}
\end{figure}

The symmetry of the fissioning system was determined via the calculated masses of the fission products. For an event to fall into the asymmetric bin, the following condition had to be meant for one of the fission fragments,
\let\saveeqnno\theequation
\let\savefrac\frac
\def\displaystyle\savefrac{\displaystyle\savefrac}
\begin{eqnarray}
\let\frac\displaystyle\savefrac
\gdef\theequation{9}
\let\theHequation\theequation
\label{disp-formula-group-67de663e2f3b460bb818280cfe82ae45}
\begin{array}{@{}l}A_{range}\;=\;\left(\frac{A_{CN}\;-\;v_{tot,En}}2\;+\;20\right)\pm\;5\;amu\end{array}
\end{eqnarray}
\global\let\theequation\saveeqnno
\addtocounter{equation}{-1}\ignorespaces
where $_{vtot,En}$ is the total neutron multiplicity as calculated by GEF \cite{9} at the corresponding En. The second term (+ 20) is the calculated difference between the A$_H$  peak in asymmetric fission and the single peak corresponding to symmetric fission. If both fission products fall within a finite mass range ($S_{range}$) the event passed into the symmetric bin, where $S_{range}$  is defined as,

\ensuremath{_{}}is defined as,
\let\saveeqnno\theequation
\let\savefrac\frac
\def\displaystyle\savefrac{\displaystyle\savefrac}
\begin{eqnarray}
\let\frac\displaystyle\savefrac
\gdef\theequation{10}
\let\theHequation\theequation
\label{disp-formula-group-80ec323a08ab4b359c054f9a64b27603}
\begin{array}{@{}l}\begin{array}{l}\begin{array}{l}S_{range}\;=\;\frac{A_{CN}\;-\;v_{tot,En}}2\;\pm\;5\;amu\\\\\end{array}\\\end{array}\end{array}
\end{eqnarray}
\global\let\theequation\saveeqnno
\addtocounter{equation}{-1}\ignorespaces

\begin{figure}
 \includegraphics[width=.5 \textwidth, trim=5cm 1.5cm 3cm 4.3cm, clip=true]{Figure8}
 \vspace{-3cm}
\caption{(a)  TKE vs. Em with fission events separated into symmetric and asymmetric events. (b)Dependence of asymmetric TKE events on neutron energy.(Color online)}
\label{fig:8}
\end{figure}

In Figure 8, we show the TKE$_{sym}$ and TKE$_{asym}$ events and plot them against En.  In Figure 8 (inset) we  see that TKE$_{asym}$ drops off rather sharply above En = 10 MeV. Across the entire En range measured, TKE$_{asym}$ decreased from 172.5 MeV to 165.6 MeV. However , the TKEs release remained virtually constant. This suggests that the overall trend of decreasing TKE is not exclusively a consequence of the onset of symmetric fission but is instead also impacted by a decreasing TKEasym as well. Hambsch theorized that this decreasing asymmetric TKE trend observed in the$^{237}$Np(n,f) reaction could involve a more drastic deformation of fission fragments, particularly for 130 $<$ A $<$ 145 amu range, as En increases \cite{3}.

\section{CONCLUSION}
Previous TKE studies of $^{237}$Np (n,f) are limited and were restricted to  incident neutron energies less than 6 MeV. The goal of this work was to build upon the existing TKE data available for $^{237}$Np by extending measurements into the fast neutron energy range, up to 100 MeV. \newline
\indent A vapor deposited $^{237}$NpF$_4$ target was irradiated at LANSCE-WNR to measure the TKE dependence for the $^{237}$Np (n,f) reaction over  En 2.6 -100 MeV. The TKE values  reported in this work were compared to that of \cite{2,3} and the prediction of the GEF model \cite{9}. We also report fission fragment mass distributions and gained additional detail about the fissioning system by separating the fission path by symmetry regimes. We attribute the observed trend of decreasing TKE to the emergence of symmetric fission, as well as the rapid drop off of TKE$_{asym}$ as excitation increases. The peculiar “bump” in TKE at En = 45-70 MeV, while consistent with other actinide studies within this framework, is not yet understood and needs further exploration.