text:
\section{Introduction}
 During the last decade, the models based on the concepts of critical point symmetries (CPS) related to shape phase transitions provide a very interesting theoretical framework for studies of nuclear structure phenomena. Actually, this interest has increased even more with the insertion of an additional critical point symmetries. Shape phase transitions  have been first considered  in  the framework of the interacting boson model \cite{Iachello1}, which describes collective states of nuclei in terms of collective bosons of angular momentum zero (s-boson) and two (d-boson) in the context of a U(6) overall symmetry, having a dynamical   U(5) (vibrational), SU(3) (prolate deformed rotational or axial rotor) and O(6) ($\gamma$-unstable) as limiting symmetries.  Another important symmetries called E(5) \cite{Iachello2} and X(5) \cite{Iachello3}, which approximate special solutions of the Bohr-Mottelson model \cite{Bohr} with an infinite-well potential and which were offered for the critical points of the shape phase transitions U(5)$\leftrightarrow$O(6) and U(5)$\leftrightarrow$SU(3) respectively, have been realized by Iachello. Later, a $\gamma$ -rigid (with $\gamma=0$ ) version  of the critical symmetry X(5), called  X(3) have been introduced in \cite{Bonatsos1}. Other models considering the extension of X(3) , such as X(3)-$\beta^{2n}$ (n=1,2) \cite{Budaca1} and X(3)-$\beta^6$ \cite{Budaca2} have also been developed not long ago. Besides, several additional attempts have been done to obtain solutions of the Bohr Hamiltonian with a constant mass parameter \cite{Chabab1,Chabab2} as well as within the deformation dependent mass formalism \cite{Chabab3,Bonatsos2}.\\
Recently, a lot of attention has been attracted by the quantum mechanical problems implying a generalized modified commutation relations which includes a minimal length or Generalized	Uncertainty Principle (GUP). Such an important idea was motivated by noncommutative geometry \cite{Witten,Nathan} in the quantum gravity \cite{QGravity1,QGravity2,QGravity3} and the string theory context \cite{StringTh1,StringTh2,StringTh3}. However, the concept of minimal length can be incorporated in the study of  physical systems by considering the deformed canonical commutation relation,
\begin{equation}
\left[X,P\right]=i\hbar\left(1+\alpha^2P^2\right)
\label{E1}
\end{equation}
here $\alpha$  represents the minimal length parameter (a very small positive parameter). This commutation relation leads to the uncertainty relation
\begin{equation}
\Delta X \Delta P\geq \frac{\hbar}{2}\left(1+\alpha\left(\Delta P\right)^2\right)
\label{E2}
\end{equation}
which implies the existence of a minimal length given by $\left(\Delta X\right)_{min}=\hbar \sqrt{\alpha}$.
It should be noted that, since the elaboration of the fundamental principles of the quantum mechanics with GUP in \cite{GUP1,GUP2,GUP3,GUP4}, a much development, in this direction, has been accomplished in order to study   the effect of the minimal length on  quantum systems as well as on  classical ones. Nevertheless, only  few problems are shown to be solved exactly or  approximately. Among them  one can cite the Schrodinger equation for : the harmonic oscillator \cite{Problem1}, the hydrogen atom \cite{Problem2,Problem3,Problem4,Problem5,Problem6}, the inverse square potential \cite{Problem7}, the scattering problem by Yukawa and Coulomb potentials \cite{Problem8} and  square well potential \cite{Problem9}.
In this Letter,  we  study a $\gamma$-rigid version of the Bohr-Mottelson Hamiltonian, within an infinite square well potential  in $\beta$ collective shape variable   as X(3) model, in the presence of a minimal length. Particularly, we investigate the effect of a minimal length on the physical observables such as energy spectrum and eigenfunctions as well as B(E2) electromagnetic transition rates.
	\label{intro}
	\section{Minimal length formalism}
The theoretical background of minimal length formalism (MLF)  motivated by a Heisenberg algebra and implying a generalized uncertainty principle (GUP)	has been considered recently in \cite{Problem1,Problem9}. In the framework of this formalism, the generalization of the deformed canonical commutation relation (\ref{E1}) is given by \cite{Problem1,Problem9}
	\begin{equation}
	\left[\hat{X_i},\hat{P_j}\right]=i\hbar\left(\delta_{ij}+\alpha\hat{P}^2\delta_{ij}+\alpha'\hat{P_i}\hat{P_j}\right)
	\label{E3}
	\end{equation}
	where $\alpha'$ is an additional parameter which is of the order of $\alpha$. In this case, the components of the momentum operator commute to one another
	\begin{equation}
	\left[\hat{P_i},\hat{P_j}\right]=0
	\label{E4}
	\end{equation}
	However, the commutator between two position operators is in general different from zero
	\begin{equation}
	\left[\hat{X_i},\hat{X_j}\right]=i\hbar\frac{(2\alpha-\alpha')+(2\alpha+\alpha')\alpha\hat{P}^2}{1+\alpha\hat{P}^2}\left(\hat{P_i}\hat{X_j}-\hat{P_j}\hat{X_i}\right)
	\label{E5}
	\end{equation}
	It is clear that the generalized canonical commutation relation (\ref{E3}) leads to the minimal observable length $(\Delta X_i)_{min}=\hbar\sqrt{3\alpha+\alpha'}$.
	In the same context, we have different representations for the canonical operators $X_i$ and $P_i$. Among these representations, one can cite the momentum space representation \cite{Problem1}:
	\begin{equation}
	\hat{X_i}=i\hbar\left[(1+\alpha p^2)\frac{\partial}{\partial p_i}+\alpha'p_ip_j\frac{\partial}{\partial p_j}\right]+\eta p_i,\ \hat{P_i}=p_i
	\label{E6}
	\end{equation}
	and the position representation given by \cite{Problem3,Problem4}:
	\begin{equation}
	\hat{X_i}=\hat{x_i}+\frac{(2\alpha-\alpha')\left(\hat{p}^2\hat{x_i}+\hat{x_i}\hat{p}^2\right)}{4},\  \hat{P_i}=\hat{p_i}\left(1+\frac{\alpha'}{2}\hat{p}^2\right)
	\label{E7}
	\end{equation}
	where $\hat{x_i}$  and $\hat{p_i}$  are the usual position and momentum operators respectively, which  obey the following
	relations $ \left[\hat{x_i},\hat{p_j}\right]=i\hbar \delta_{ij}$ and $\hat{p}^2=\sum_i \hat{p_i}$. Note that  in the case of $\alpha'=2\alpha$, and for the first order on $\alpha$, the following canonical commutator $ \left[\hat{X_i},\hat{X_j}\right]$ vanishes.  As a consequence, Eq. (\ref{E7}) redruces to
	\begin{equation}
	\hat{X_i}=\hat{x_i},\  \hat{P_i}=\left(1+\alpha\hat{p}^2\right)\hat{p_i}
	\label{E8}
	\end{equation}
	In addition, we can interpret $p_i$ and $P_i$ shown in Eq. (\ref{E8}) according to string theory:
	 $p_i$ is the momentum operator at low energies  and $P_i$ is the momentum operator at high energies. Moreover, $p$ is the magnitude of the $p_i$ vector.
	\section{Bohr-Mottelson model with a minimal length }
	\label{sec:2}
	In the context of the collective geometrical  model of Bohr-Mottelson\cite{Bohr}, the classical expression for the rigid-body kinetic energy associated with the rotation and surface deformations of	a nucleus has the form \cite{Bohr,Sitenko}
	\begin{equation}
	\hat{T} = \frac{1}{2}\sum_{k=1}^{3} {\cal J}_k\, \omega^{\prime2}_k +
	\frac{B_m}{2}\,(\dot{\beta}^2+\beta^2 \dot{\gamma}^2),
	\label{E9}
	\end{equation}
	where $\beta$ and $\gamma$ are the usual collective variables, $B_m$ is the
	mass parameter. Also,
	\begin{equation}
	{\cal J}_k = 4B_m\beta^2 \sin^2\bigl(\gamma - {\textstyle\frac{2}{3}}\pi k\bigr)
	\label{E10}
	\end{equation}
	are the three principal irrotational moments of inertia,
	and $\omega^\prime_k$ ($k=1$, 2, 3) are the components of the angular velocity (angular frequencies)
	on the body-fixed $k$-axes,
	which can be expressed in terms of the time derivatives of the Euler angles
	$\dot{\phi}, \dot{\theta}, \dot{\psi}$ \cite{Sitenko,Zare}
	\begin{eqnarray}
	\omega^\prime_1 &=& -\sin\theta \cos\psi\,\dot{\phi} + \sin\psi\,\dot{\theta},
	\nonumber\\
	\omega^\prime_2 &=& \sin\theta \sin\psi\,\dot{\phi} + \cos\psi\,\dot{\theta},\\
	\omega^\prime_3 &=& \cos\theta\,\dot{\phi} + \dot{\psi}. \nonumber
	\label{E11}
	\end{eqnarray}
	Going further, by  assuming the nucleus to be $\gamma$-rigid (i.e. $\dot{\gamma}=0$), as a non-adiabatic approach proposed by Davydov and Chaban in \cite{Davydov1}, and considering in particular the axially symmetric prolate case of $\gamma=0$, we see that
	the third irrotational moment of inertia ${\cal J}_3$ vanishes, while
	the other two become equal ${\cal J}_1 = {\cal J}_2 = 3B_m\beta^2$, thus the kinetic
	energy of Eq. (\ref{E9}) is simply \cite{Sitenko,Davydov2}
	\begin{align}
	\hat{T} =& \frac{3}{2} B_m\beta^2 (\omega^{\prime2}_1 + \omega^{\prime2}_2)
	+ \frac{B_m}{2}\,\dot{\beta}^2\nonumber \\
	=& \frac{B_m}{2}\Bigl[3\beta^2(\sin^2\theta\,\dot{\phi}^2 + \dot{\theta}^2)
	+ \dot{\beta}^2 \Bigr].
	\label{E12}
	\end{align}
	Since in the case of axial symmetry the nucleus can rotate only about directions perpendicular to the symmetry axis,  the collective motions in the nucleus are characterized by only three degrees of freedom:  $q_1=\phi$, $q_2=\theta$, and $q_3=\beta$.	Having in mind the position space representation (\ref{E8}), the kinetic energy operator, in this case, can be expressed in terms of the Laplacian and bi-Laplacian operators as follows
	\begin{align}
	T=-\frac{\hbar^2}{2B_m}\left(1-2\alpha\hbar^2\frac{1}{\sqrt{g}}\sum_{ij}\frac{\partial}{\partial q_i}\sqrt{g}g_{ij}^{-1}\frac{\partial}{\partial q_j}\right)\nonumber\\ \times\left[\frac{1}{\sqrt{g}}\sum_{ij}\frac{\partial}{\partial q_i}\sqrt{g}g_{ij}^{-1}\frac{\partial}{\partial q_j}\right]
    \label{E13}
	\end{align}
		where the matrix $g_{ij}$ having a diagonal form
		\begin{equation}
			g_{ij} = \left(\begin{array}{ccc} 3\beta^2\sin^2\theta & 0 & 0 \\ 0 &
				3\beta^2 & 0 \\0 & 0 & 1 \end{array}\right)
				\label{E14}
		\end{equation}
		where $g$ is the determinant of the matrix $g_{ij}$ and $g_{ij}^{-1}$ is the inverse matrix of  $g_{ij}$.
	Using the general procedure of quantization (Pauli–Podolsky prescription) in curvilinear coordinates, we obtain, in compact form, the collective Hamiltonian operator, up to the first order of $\alpha$,
	\begin{align}
	\hat{H}=-\frac{\hbar^2}{2B_m}\Delta +\frac{\alpha\hbar^4}{B_m}\Delta^2+ V(\beta)
	\label{E15}
	\end{align}
	with
	\begin{align}
	\Delta = \Biggl[\frac{1}{\beta^2}
	\frac{\partial}{\partial\beta}\beta^2\frac{\partial}{\partial\beta} +
	\frac{1}{3\beta^2} \Delta _\Omega\Biggr]
	\label{E16}
	\end{align}
	where $\Delta_\Omega$ is the angular part of the Laplace operator
	\begin{equation}
	\Delta_\Omega = \frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\sin\theta
	\frac{\partial}{\partial\theta} + \frac{1}{\sin^2\theta}
	\frac{\partial^2}{\partial\phi^2}.
   \label{E17}
	\end{equation}
	The corresponding deformed Schr\"odinger equation to the first order on $\alpha$ reads as
	\begin{equation}
	\Biggl[-\frac{\hbar^2}{2B_m}\Delta +\frac{\alpha\hbar^4}{B_m}\Delta^2+ V(\beta)-E \Biggr] \Psi(\beta,\theta,\phi)=0
   \label{E18}
	\end{equation}
	which is a second order differential equation.
	In addition, we can see, here, that is difficult to obtain analytic solution of this differential equation, because of the bi-Laplacian $\Delta^2\propto p^4$ . However, we can get rid of  the term $\Delta^2$  in equation (\ref{E18}) by introducing an auxiliary wave function $\Phi$ as in \cite{Problem9}, so that
	\begin{equation}
	\Psi(\beta,\theta,\phi)=\left[1-2\alpha(i\hbar)^2\Delta\right]\Phi(\beta,\theta,\phi)
	\label{E19}
	\end{equation}
	Thus, we obtain the following differential equation satisfied by $\Phi$,
	\begin{equation}
	\Biggl[\left(1+4B_m\alpha\left(E-V(\beta)\right)\right)\Delta+\frac{2B_m}{\hbar^2}\left(E-V(\beta)\right)\Biggr] \Phi(\beta,\theta,\phi)=0
	\label{E20}
	\end{equation}
	where $\Delta$ is defined by Eq. (\ref{E16}).
	 The  latter equation can be  solved by using the usual following factorization
	\begin{equation}
	\Phi(\beta,\theta,\phi) = F_{n_{\beta}}(\beta)\, Y_{LM}(\theta,\phi),
	\label{E21}
	\end{equation}
	where $Y_{LM}(\theta,\phi)$ are the spherical harmonics. Then
	the angular part leads to the equation
	\begin{equation}
	\Delta_\Omega Y_{LM}(\theta,\phi) =- L(L+1)Y_{LM}(\theta,\phi),
    \label{E22}
	\end{equation}
	where $L$ is the angular momentum quantum number,
	while  the radial part $F(\beta)$ obeys to:
	\begin{equation}
	\Biggl[\frac{1}{\beta^2}
	\frac{d}{d\beta}\beta^2\frac{d}{d\beta}
	-\frac{L(L+1)}{3\beta^2} + \frac{2B}{\hbar^2}\bar{K}(E,\beta)\Biggr]F_{n_{\beta}}(\beta) = 0.
	\label{E23}
	\end{equation}
	with
	\begin{equation}
	\bar{K}(E,\beta)=\left(\frac{E-V(\beta)}{\left(1+4B_m\alpha\left(E-V(\beta)\right)\right)}\right)
    \label{E24}
	\end{equation}
	and $n_{\beta}$ is the radial quantum number. Eq. (\ref{E23}) is an effective Schrödinger equation including the minimal length. It should be noticed that in the limit $\alpha \rightarrow 0$, Eq. (\ref{E23}) reduces to the ordinary collective Schrödinger equation \cite{Bonatsos1,Budaca1,Budaca2}.\\
	In what concerns the $\beta$ degree of freedom, we will consider here an anharmonic
	behaviour reflected into an infinite square well shape of the potential as in the case of  X(3) symmetry \cite{Bonatsos1}:
	\begin{align}
	V(\beta) = \left\{ \begin{array}{lcl} 0, && \rm{if}\ \beta\leq\beta_{\omega} \\
	\infty, &&  \rm{if}\ \beta > \beta_{\omega} \end{array} \right. ,
	\label{E25}
	\end{align}
	where $\beta_{\omega}$ indicates the width of the well.
	In this case the wave function $F(\beta)$ is a solution of the equation
	\begin{equation}
	\Biggl[\frac{d^2}{d\beta^2} + \frac{2}{\beta}\frac{d}{d\beta}
	+ \Biggl(\bar{k}-\frac{L(L+1)}{3\beta^2}\Biggr)\Biggr]F_{n_{\beta}}(\beta) = 0
	\label{E26}
	\end{equation}
	in the interval $0\leq\beta\leq\beta_{\omega}$,
	where we introduced the reduced energies
	\begin{equation}
	\varepsilon=\bar{k}=\frac{2B_m}{\hbar^2}\cdot\frac{E}{(1+4B_m\alpha E)}
    \label{E27}
	\end{equation}
	 while it vanishes outside.
	Substituting $F_{n_{\beta}}(\beta)=\beta^{-1/2} f_{n_{\beta}}(\beta)$ in Eq.(\ref{E26}), one obtains the Bessel equation
	\begin{equation}
		\Biggl[\frac{d^2}{d\beta^2} +
		\frac{1}{\beta}\frac{d}{d\beta}
		+ \Biggl(\bar{k}^2-\frac{\eta^2}{\beta^2}\Biggr)\Biggr]
		 f_{n_{\beta}}(\beta) = 0,
	\label{E28}
	\end{equation}
	with
	\begin{equation}
		\eta=\left(\frac{L(L+1)}{3}+\frac{1}{4}\right)^{\frac{1}{2}},
	\label{E29}
	\end{equation}
	and the boundary condition being $ f_{n_{\beta}}(\beta_{\omega})=0$.
	The solution of Eq.(\ref{E26}), which is finite at $\beta=0$, is then given by
	\begin{equation}
	F_{n_{\beta}}(\beta)=F_{sL}(\beta) = N_{s,L}\,\beta^{-1/2}
		J_{\eta}(\bar{k}_{s,\eta}\beta),  \ s=n_{\beta}+1
	\label{E30}
	\end{equation}
	with $ \bar{k}_{s,\eta}=\chi_{s,\eta}/\beta_{\omega}$ and $\varepsilon_{s,\eta}=\bar{k}_{s,\eta}^2$,
	where $\chi_{s,\eta}$ is the $s$-th zero of the Bessel function of the first kind
	$J_{\eta}(\bar{k}_{s,\eta}\beta_{\omega})$. $N_{s,L}$ is a normalization constant  to be determined later.
	The corresponding spectrum is then
	\begin{equation}
		E_{s,L} = \frac{\hbar^2}{2B_m}\times\frac{\bar{k}_{s,\eta}^2}{1-2\hbar^2\alpha\bar{k}_{s,\eta}^2},\ \ \bar{k}_{s,\eta}= \frac{\chi_{s,\eta}}{\beta_{\omega}}
		\label{E31}
	\end{equation}
	In the above equation, the term $2\hbar^2\alpha\bar{k}_{s,\eta}^2$ is the correction due
	to the minimal length. Therefore, we conclude that the minimal
	length increases slightly the energy spectrum. In addition, the relative
	correction can be written as
	\begin{equation}
	\frac{\Delta 	E_{s,L}  }{	E_{s,L}^{0} }=\frac{2\hbar^2\alpha\bar{k}_{s,\eta}^2}{1-2\hbar^2\alpha\bar{k}_{s,\eta}^2}
	\label{E32}
	\end{equation}
where  $E_{s,L}^{0}=\lim_{\alpha \rightarrow 0}E_{s,L}$ .\\
	Essentially, the total wave function (\ref{E19}) can be written as
	\begin{align}
		\Psi(\beta,\theta,\phi)&=\left[1-2\alpha(i\hbar)^2\Delta\right]\Phi(\beta,\theta,\phi)\nonumber \\

		                      &=\left(1+2(i\hbar)^2\alpha\bar{k}_{s,\eta}^2\right)F_{n_{\beta}}(\beta) Y_{LM}(\theta,\phi)
    \label{E33}
	\end{align}
	Finally, we have
\begin{align}
\Psi(\beta,\theta,\phi)
=N_{s,L}\left(1+2(i\hbar)^2\alpha\bar{k}_{s,\eta}^2\right)\beta^{-1/2}
J_{\eta}(\bar{k}_{s,\eta}\beta) Y_{LM}(\theta,\phi)
\label{E34}
\end{align}
Using the normalization condition of this function, we easy obtain the factor $N_{s,L}$ :
\begin{equation}
N_{s,L}=\frac{\sqrt{2}}{\beta_{\omega}\,J_{\eta+1}(\chi_{s,\eta})\left(1+2(i\hbar)^2\alpha\bar{k}_{s,\eta}^2\right)}
\label{E35}
\end{equation}
Having the analytical expression of the normalized wave function, one can readily compute the B(E2) transition probabilities. Nevertheless, it should be remarked that  the full normalized wave function does not change by introducing  the concept of minimal length. Therefore, the B(E2) transition probabilities, which are expressed as,
\begin{equation}
	B(E2; sL \to s'L') = \frac{1}{2L+1}
	\left|\langle s'L'||T^{(E2)}||sL\rangle\right|^2
	\label{E36}
\end{equation}
 also remain unchanged by this formalism  and are similar to those obtained in \cite{Bonatsos1} where  $T^{(E2)}_\mu = t\,\beta\,\sqrt{\frac{4\pi}{5}}\,Y_{2\mu}(\theta,\phi)$ is the quadrupole operator for $\gamma= 0$ and $t$ is a scaling factor.
Here, some remarks concerning  X(5)  with a minimal length concept are worth to be mentioned:
\begin{itemize}
\item (1) : notice that, in this case,
the same Eq. (\ref{E28}) occurs, but with $\eta=\left( {L(L+1)\over 3} +{9\over 4}\right)^{1/2} $.
\item (2) : As in the case of X(3) model , the concept of minimal length  has no effect on the B(E2) transition probabilities of X(5).
\end{itemize}
Besides, from the requirement that the wave function be symmetric with respect to the perpendicular plan to the symmetry axis of the nucleus and passing through its center, it follows that only even values of the angular momentum $L$ are allowed. Therefore no $\gamma$ bands appear in the present models as expected, because the $\gamma$ degree of freedom has been initially frozen to $\gamma=0$.

\setlength{\tabcolsep}{5.7pt}
\begin{table}[H]
	\caption{Typical energy levels ( ground state ) of the X(3)-ML and X(5)-ML models, normalized to the $2_{g}^+$ excited state energy  for different
		values of the parameter $\alpha$ with $\hbar=1$.}
	\label{table1}
	\begin{center}
		{\linespread{1.93}
			\footnotesize
			\begin{tabular}{|c|c|c|c|}
				\hline
				L&\multicolumn{3}{|c|}{X(3)-ML}\\
				\hline
				$0^{+}$&   0.000 & 0.000 & 0.000 \\
				$2^{+}$&   1.000 & 1.000 & 1.000 \\
				$4^{+}$&   2.445 & 2.455 & 2.465\\
				$6^{+}$&   4.234 & 4.274 & 4.315\\
				$8^{+}$&   6.348 & 6.448 & 6.551 \\
				$10^{+}$&  8.779 & 8.980 & 9.194\\
				$12^{+}$&  11.520 & 11.880 & 12.270\\
				\noalign{\smallskip}\hline\noalign{\smallskip}
				$\beta_{\omega}$& -- & 60.0& 60.0\\
				$\alpha$& 0 & 0.5& 1\\
				\hline

			\end{tabular}
			}
			{\linespread{1.93}
				\footnotesize
				\begin{tabular}{|c|c|c|c|}
					\hline
					L&\multicolumn{3}{|c|}{X(5)-ML}\\
					\hline
					$0^{+}$&   0.000 & 0.000 & 0.000 \\
					$2^{+}$&   1.000 & 1.000 & 1.000 \\
					$4^{+}$&   2.904 & 2.914 & 2.925\\
					$6^{+}$&   5.430 & 5.477 & 5.526\\
					$8^{+}$&   8.483 & 8.610 & 8.741 \\
					$10^{+}$&  12.03 & 12.290 & 12.570\\
					$12^{+}$&  16.04 & 16.530 & 17.050\\
					\noalign{\smallskip}\hline\noalign{\smallskip}
					$\beta_{\omega}$& -- & 60.0& 60.0\\
					$\alpha$& 0 & 0.5& 1\\
					\hline

				\end{tabular}
			}

				\end{center}
			\end{table}
			\setlength{\tabcolsep}{13.8pt}
			\begin{table}[H]
				\caption{The values of the free parameters used in the calculations.}
				\label{table2}
				\begin{center}
					{\linespread{1.6}
						\footnotesize
						\begin{tabular}{|c|c|c|c|c|c|}
							\hline
							Models&\multicolumn{2}{c|}{X(3)-ML}	&\multicolumn{2}{c|}{X(5)-ML}\\
							\hline
							$Nucleus$&$\alpha$&$\beta_{\omega}$&$\alpha$&$\beta_{\omega}$\\
							\hline
							$ {}^{150}$Nd & 0.961& 29.446& 0.184 & 67.308\\
							$ {}^{176}$Os & 0.421& 42.517& 0.000 & 64.670\\
							$ {}^{178}$Os & 0.444& 38.575& 0.649 & 75.614\\
							$ {}^{180}$Os & 0.999& 21.858& 0.000 & 56.102\\
							$ {}^{156}$Dy & 0.833& 50.763& 0.000 & 95.399\\
							$ {}^{154}$Gd & 0.654& 60.299& 0.233 & 65.648\\
							\hline
						\end{tabular}}
					\end{center}
				\end{table}
				\setlength{\tabcolsep}{4.3pt}
	\section{Model applicability and numerical results}
Because the $\gamma$ degree of freedom has been frozen to $\gamma=0$, the bands in the present models, like in X(3) model, are only classified by the principal quantum number $n_{\beta}$ or $s=n_{\beta}+1$. A few interesting low-lying bands are given as
	\begin{itemize}
		\item i) The energy levels of the ground state band with $s = 1$,
		\item ii)  The $\beta$-vibrational bands with $s > 1$.
	\end{itemize}
In order to avoid any ambiguity of the nomenclature between our models and the existing phenomenological models, namely: X(3) and X(5) ,  we denote X(3)-ML and X(5)-ML  in connection with  X(3) and X(5) respectively.
 The proposed models have two free parameters, namely: the minimal length parameter $\alpha$ and the width of the infinite square well potential $\beta_{\omega}$. Obviously, we do not count the mass parameter $B_m$ since it disappears when calculating the energy ratios. However, according to the general form of the obtained energy spectrum, these parameters  could be dependent from each other and  check a constraint. Indeed, the energy  spectrum  corresponding to our models, where the effect of the minimal length  is considered, is always positive $E_{s,L}\ge 0$ (this is also valid in the ordinary case i.e: without a minimal length scenario). Due to this fact, we can write:
\begin{equation}
1-2\hbar^2\alpha\bar{k}_{s,\eta}^2>0, \ \bar{k}_{s,\eta}= \frac{\chi_{s,\eta}}{\beta_{\omega}}
\label{E37}
\end{equation}
which is a constraint  between $\alpha$ and $\beta_{\omega}$. From practical point of view, it is important to note that the value of $\alpha$ must  be very small compared to the width of the well $\beta_{\omega}$ in order to preserve the mentioned above constraint.
In Fig.~\ref{Fig1}, the energy  of the first $4^+$ and $6^+$  levels, of X(3)-ML and X(5)-ML models, for two values of the width $\beta_{\omega}=5$ and $\beta_{\omega}=40$, are displayed as function of the minimal length parameter $\alpha$ in the interval $[0, 1]$.
In the case of small value of $\beta_{\omega}=5$ , we see that, the energy  ratios of the first $4^+$ and $6^+$  levels, for X(3)-ML, present a singularity nearby  $\alpha=0.3561$ and $\alpha=0.2333$ respectively, because the condition (\ref{E37}), in this case, is not fulfilled. Likewise in the  case of X (5)-ML, but in this time, the singularity occurs around  the following values $\alpha=0.3087$ and $\alpha=0.2149$.
 While in the  case of a large value of $\beta_{\omega}=40$, where the above relationship is very well checked, the energy of the first $4^+$ and $6^+$  levels is very much influenced by $\alpha$ .
In addition, Table.~\ref{table1}  shows  a typical energy levels of ground state  of the X(3)-ML and X(5)-ML models, normalized to the $2_{g}^+$ excited state energy  for $\beta_{\omega}=60$.
From this table, one can see the effect of the minimal length becomes manifest for higher values of the angular momentum. Indeed, such a fact, which results from the uncertainty principle Eq. (\ref{E2}) as expected from string theory, is well illustrated schematically in Fig.~\ref{Fig2} where the evolution of the energy spectrum of the ground state and the $\beta_1$ bands, normalized to the first $2^+$ excited state, is presented. Furthermore, one can see that the effect of the minimal length is more important for the X (3) symmetry than for the X (5) one. Such an effect could be beneficial when trying to reproduce the experimental data for concrete nuclei in comparison, particularly, with the pure X (3) model as it can be seen subsequently. Moreover, from this figure one can see that the ground state band  as well as $\beta_1$ band are very much influenced by $\alpha$ for higher angular momentum. Besides, as is mentioned above, the minimal length effect increases slightly the energy spectrum.
  In Fig.~\ref{Fig3}, we present the variations of the relative correction of our model to the X(3) symmetry given by Eq. (\ref{E32}), as a function of the angular momentum L and  the minimal length  as well as the width $\beta_{\omega}$. The map contour lines are lines with a constant relative correction. The area delimited by two successive contour lines represents the recovery rate of the X(3) symmetry by our model. From Fig.~\ref{Fig3}, one can see that in the vicinity of  $\alpha\rightarrow 0$ and for lower values of the angular momentum L, the recovery area is large . So, in this region our model is identical to the X(3) one. But, as one goes in the same given region to higher values of L, such an area narrows. Also, as $\alpha$ increases, the recovery area starts to contract. So, the gap between our model and the X(3) one increases, as it was  mentioned above in the comment on Table~\ref{table1}. However, this gap between both models is worthwhile for ours insofar as it allows reproducing the experimental data, by our model, with a good precision in comparison with the pure X(3) model as it can be seen from Fig.~\ref{Fig4}, Fig.~\ref{Fig5} and Fig.~\ref{Fig6}.
  In the right panel of Fig.~\ref{Fig3}, given for $\beta_{\omega}=400$, we observe a similar behavior as in the left one for which $\beta_{\omega}=40$ but with a bit more contracted recovery areas and lower values of the relative correction corresponding to the contour lines. This is due to the fact that for a deeper square well, the minimal length becomes smaller in concordance with the constraint (\ref{E37}).
 As a result, the models, developed here, allow to describe properties of nuclei having the  signature $R_{4/2}=E(4_g^+)/E(2_g^+) \geq 2.44$, unlike the model developed in Ref \cite{Budaca1,Budaca2}  which studies a few  properties of nuclei having in this case the ratio $R_{4/2}< 2.44$.
The experimental realization of the models was found to occur in some nuclei $ {}^{150}$Nd,$ {}^{176-180}$Os,$ {}^{156}$Dy and $ {}^{154}$Gd, where the values of the used free parameters  in the calculations are listed in Table ~\ref{table2}.
In Fig.~\ref{Fig4}, Fig.~\ref{Fig5} and Fig.~\ref{Fig6}, we present the numerical results for the energy of the ground state and the $\beta_1$ bands, normalized to the energy of the first excited level, obtained within X(3)-ML and X(5)-ML for these nuclei. From these figures, on can see the X(3)-ML model reproduces well the experimental data in comparison with the pure X(3) one in both bands. While, the X(5)-ML model is generally  identical to the X(5) one and slightly better in the case of ${}^{150}$Nd ,${}^{178}$Os and ${}^{154}$Gd nuclei. Such a difference, in the precision of predictions of both models: X(3)-ML and X(5)-ML, is due to the parameter $\eta$ which enters in the zeros of the Bessel function $\chi_{s,\eta}$ and so defining the energy through the quantity $\bar{k}_{s,\eta}$ Eq. (\ref{E31}).
Indeed, for a given L, the value of $\eta$  is lower in the X(3)-ML than in the X(5)-ML. Hence, the minimal length formalism seems to be more suitable for studying $\gamma$-rigid nuclei in the frame of the X(3) symmetry.

		\begin{figure*}[tbph]
			\centering
			\rotatebox{0}{\includegraphics[height=50mm]{Fig8.eps}}
			\rotatebox{0}{\includegraphics[height=50mm]{Fig9.eps}}
			\rotatebox{0}{\includegraphics[height=50mm]{Fig10.eps}}
			\rotatebox{0}{\includegraphics[height=50mm]{Fig11.eps}}
			\caption{The enegy of the first $4^+$ and $6^+$  levels are ploted as function of  the minimal length parameter $\alpha$. }
			\label{Fig1}
		\end{figure*}
	\begin{figure*}[tbph]
		\centering
		\rotatebox{0}{\includegraphics[height=50mm]{Fig4.eps}}
		\rotatebox{0}{\includegraphics[height=50mm]{Fig5.eps}}
		\rotatebox{0}{\includegraphics[height=50mm]{Fig6.eps}}
		\rotatebox{0}{\includegraphics[height=50mm]{Fig7.eps}}
		\caption{The energy of the ground state  and the $\beta_1$ band, normalized to the energy of the first excited state in the X(3)-ML and X(5)-ML models are ploted as function of  angular momentum $L$ for different values of the minimal length parameter $\alpha$. The ground state is labeled by $R_{L/2}$, while $\beta_1$ band is labeled by $R_{L_{\beta}/2}$. The X(3)  and X(5) predictions are also shown for
			comparison. }
		\label{Fig2}
	\end{figure*}
	 \begin{figure*}[tbph]
	 	\centering
	 	\rotatebox{0}{\includegraphics[height=73mm]{Fig12.eps}}
	 	\rotatebox{0}{\includegraphics[height=73mm]{Fig13.eps}}
	 	\caption{Map contour lines of  the relative correction  (\ref{E32}) for the X(3)-ML model  drawn as a function  of  the angular momentum $L$ and the minimal length parameter $\alpha$ for $\beta_{\omega}=40$ (left) and $\beta_{\omega}=400$ (right).}
	 	\label{Fig3}
	 \end{figure*}
		\begin{figure*}[tbph]
			\centering
			\rotatebox{0}{\includegraphics[height=30mm]{Fig14.eps}}
			\rotatebox{0}{\includegraphics[height=30mm]{Fig15.eps}}
            \rotatebox{0}{\includegraphics[height=30mm]{Fig16.eps}}
            \rotatebox{0}{\includegraphics[height=30mm]{Fig17.eps}}
         	\rotatebox{0}{\includegraphics[height=30mm]{Fig22.eps}}
         	\rotatebox{0}{\includegraphics[height=30mm]{Fig23.eps}}
         	\rotatebox{0}{\includegraphics[height=30mm]{Fig24.eps}}
         	\rotatebox{0}{\includegraphics[height=30mm]{Fig25.eps}}
			\caption{Theoretical results for energy levels of the ground state and  the $\beta_1$-bands of the X(3)-ML and X(5)-ML  models, compared with the available experimental data \cite{dada} for ${}^{150}$Nd and ${}^{176}$Os. The levels of each band are normalized to the $2_1^+$ state. The ground state is labeled by $R_{L/2}$, while $\beta_1$ band is labeled by $R_{L_{\beta}/2}$.  }
			\label{Fig4}
		\end{figure*}
		 \begin{figure*}[tbph]
		 	\centering
		 	\rotatebox{0}{\includegraphics[height=30mm]{Fig26.eps}}
		 	\rotatebox{0}{\includegraphics[height=30mm]{Fig27.eps}}
		 	\rotatebox{0}{\includegraphics[height=30mm]{Fig28.eps}}
		 	\rotatebox{0}{\includegraphics[height=30mm]{Fig29.eps}}
		 	\rotatebox{0}{\includegraphics[height=30mm]{Fig30.eps}}
		 	\rotatebox{0}{\includegraphics[height=30mm]{Fig31.eps}}
		 	\rotatebox{0}{\includegraphics[height=30mm]{Fig32.eps}}
		 	\rotatebox{0}{\includegraphics[height=30mm]{Fig33.eps}}
		 	\caption{Theoretical results for energy levels of the ground state and  the $\beta_1$-bands of the X(3)-ML and X(5)-ML  models, compared with the available experimental data \cite{dada} for ${}^{178}$Os and ${}^{180}$Os. The levels of each band are normalized to the $2_1^+$ state. The ground state is labeled by $R_{L/2}$, while $\beta_1$ band is labeled by $R_{L_{\beta}/2}$.  }
		 	\label{Fig5}
		 \end{figure*}
		  \begin{figure*}[tbph]
		  	\centering
		  	\rotatebox{0}{\includegraphics[height=30mm]{Fig42.eps}}
		  	\rotatebox{0}{\includegraphics[height=30mm]{Fig43.eps}}
		  	\rotatebox{0}{\includegraphics[height=30mm]{Fig44.eps}}
		  	\rotatebox{0}{\includegraphics[height=30mm]{Fig45.eps}}
		  	\rotatebox{0}{\includegraphics[height=30mm]{Fig46.eps}}
		  	\rotatebox{0}{\includegraphics[height=30mm]{Fig47.eps}}
		  	\rotatebox{0}{\includegraphics[height=30mm]{Fig48.eps}}
		  	\rotatebox{0}{\includegraphics[height=30mm]{Fig49.eps}}
		  	\caption{Theoretical results for energy levels of the ground state  and  the $\beta_1$-bands of the X(3)-ML and X(5)-ML  models, compared with the available experimental data \cite{dada}  for ${}^{156}$Dy and ${}^{154}$Gd. The levels of each band are normalized to the $2_1^+$ state.  The ground state is labeled by $R_{L/2}$, while $\beta_1$ band is labeled by $R_{L_{\beta}/2}$.  }
		  	\label{Fig6}
		  \end{figure*}
	\section{Conclusion}
	In this work, we have derived new  solutions of the Bohr-Mottelson Hamiltonian in the $\gamma$-rigid regime  within the minimal length formalism which emerges in many higher dimension theories of quantum physics. The  recall potential of the collective $\beta$-vibrations is assumed to be equal to an infinite square well as in the standard X(3) and X(5) models. So, improved versions of the X(3) and X(5) symmetries being called X(3)-ML and X(5)-ML  are elaborated. Indeed, we have shown, through this work, that  the introduction of the minimal length formalism allows one to enhance the numerical calculation precision of physical observables, particularly the energy spectrum of nuclei in comparison with the X(3) and X(5) models. These later could be easily recovered by taking a null minimal length solutions of our models: X(3)-ML and X(5)-ML.
    \section*{References}
	\bibliographystyle{elsarticle-num}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

An important challenge in  hadron physics research is to understand the
hadron structure and reaction dynamics within QCD.
The electroweak form factors of baryon resonance give important
 informations to step forward the research in this direction.
Those informations are important in testing hadron model and
possibly the lattice QCD calculation.
We have developed a dynamical model for investigating
the pion photoproduction and electroproduction reactions
 in the delta resonance region \cite{sl1,sl2}.
In the dynamical approach\cite{sl1,sl2,sl3,nbl,surya,ky1} we solve the
scattering equations with the interactions defined from the
effective Lagrangian and/or hadron models. The dynamical
approach is different from the approaches based on the
dispersion relations[7-9] or K-matrix approach[10-13]
 in interpreting the data.
In Refs. 1,2, we have not only extracted the $N\Delta$ transition form factors
 from  the data but also we have provided an interpretation of
the  extracted form factors in terms of hadron models.
It was found in Refs. 1 and 2 that the pion cloud effects give
large contribution to the $N\Delta$ transition form factors.
The magnetic dipole $N\Delta$ transition form factor $G_M^{N\Delta}(0)$
is enhanced by about 40\% by meson cloud, which gives explanation on the
long standing discrepancy between the prediction of the constituent
quark model and empirical amplitude analyses. Furthermore,
the long range pion cloud gives soft component of the
$G_M^{N\Delta}(Q^2)$ and very pronounced enhancement of quadrupole
transition form factors $G_E^{N\Delta}(Q^2), G_C^{N\Delta}(Q^2)$
at low $Q^2$.

The axial vector response of the hadron can be studied by the weak
processes. The neutrino-induced pion production reaction
$\nu + N \rightarrow l + \pi + N$ can be used to extract such informations.
The previous investigations of the weak pion production reaction[14-22]
 have been done  by using dispersion relation approach or
K-matrix approach.
We report on our recent progress on the neutrino-induced
pion production reaction in the delta resonance region\cite{sl3}.
The purpose of this work is to develop a dynamical model for
neutrino-induced pion production reaction to extract axial vector
$N\Delta$ form factors by extending our model on pion electroproduction.
In particular we investigate possible role of the dynamical pion
cloud in solving the problem that the $N\Delta$ axial vector form factor
extracted from  neutrino reaction\cite{hhm,lmz} was found to be
 about 30\% larger than the quark model predictions.
In section 2, we will briefly describe our dynamical model.
The results of neutrino-induced pion production and the axial vector
$N\Delta$ form factor is given in section 3.

\section{Dynamical approach}
We start from the Hamiltonian of mesons and baryons fields with the
interaction Hamiltonian $H_I= \sum_{B',M',B}\Gamma^0_{B'M',B}$,
which describes absorption and emission of meson $M$ from baryon $B$ as
\begin{eqnarray}
H & = & H_0 + \Gamma_{B'M',M} + (\mbox{h.c.}).
\end{eqnarray}
Similarly the weak hadron current $J^\mu$ consists of
mesons and  baryons weak currents. To obtain a manageable reaction
theory to describe neutrino-induced pion production reaction, we apply unitary
transformation method\cite{sl1,ksh} up to the second order of the
interaction Hamiltonian.
The idea is that we eliminate 'virtual' interaction $B \rightarrow B'M'$
for $m_B < m_M' + m_B'$ from the Hamiltonian in Eq. (1)
and absorb their effects into many-body potentials.
As a result $N$ and $\pi N$ states decouple with each other up to the
order of our approximation and the effective Hamiltonian consists of
the interactions of the resonance decay/production and the many-body
potentials.
\begin{eqnarray}
H_{eff} = H_0 + v_{\pi N} + \Gamma_{\Delta, \pi N} + \mbox{h.c.}.
\end{eqnarray}
\begin{figure}[h]
\centerline{\epsfig{file=vpin.eps,width=8cm}}
\vspace*{8pt}
\caption{Graphical representation of the interactions
(a) $\Gamma_{\Delta,\pi N}$ and (b) $v_{\pi N}$.}
\end{figure}
\begin{figure}[h]
\vspace*{8pt}
\centerline{\epsfig{file=gmpi.eps,width=8cm}}
\caption{Graphical representation of the effective current
(a) $J^\mu_{\Delta N}$ and (b) $J^\mu_{\pi N}$.}
\end{figure}
Here $v_{\pi N}$ is  non-resonant $\pi N$ potential and
 $\Gamma_{\Delta, \pi N}$ describes  decay of delta into $\pi N$.
The effective hadron current $J^\mu_{eff}$ is obtained by applying the
same unitary trasformation used for $H_{eff}$ on $J^\mu$.
 $J^\mu_{eff}$ consists of the non-resonant pion production current
$J^\mu_{\pi N}$ and $N\Delta$ transition current $J^\mu_{\Delta N}$:
\begin{eqnarray}
J^\mu_{eff} = J^\mu_{\Delta N} + J^\mu_{\pi N}.
\end{eqnarray}
From the effective Hamiltonian and current,
the T-matrix of the pion production reaction can be
easily obtained by solving the coupled channel
Lippman-Schwinger equation within the $\pi N \oplus \Delta$ Fock space.
The resulting matrix element of the current  $\bar{J}^\mu_{\pi N}$,
which includes final state interaction and
 satisfies  the Watson theorem, is given as
\begin{eqnarray}
\bar{J}^\mu_{\pi N} = \bar{J}^\mu_{\pi N}(\mbox{non-res})
 + \frac{\bar{\Gamma}_{\Delta,\pi N}
\bar{J}^\mu_{\Delta N}}{W - m_\Delta - \Sigma}.
\end{eqnarray}
The non-resonant current $\bar{J}^\mu_{\pi N}(\mbox{non-res})$
is calculated only from non-resonant interactions $v_{\pi N}$ and
 $J^\mu_{\pi N}$.
The second term in Eq. (4) is resonant amplitude
with the dressed  $J^\mu + N \rightarrow \Delta$ vertex
 given as
\begin{eqnarray}
\bar{J}^\mu_{\Delta N}(W,q) & = & J^\mu_{\Delta N} +
 \int dk k^2
\frac{\bar{\Gamma}_{\Delta, \pi N}(W,k)J^\mu_{\pi N}(k,q)}
     {W - E_N(k) - E_\pi(k) + i \epsilon}.
\end{eqnarray}
An important feature of the dynamical model is that the
bare vertex  $J^\mu_{\Delta N}$, which
 may be compared with the prediction of the hadron model,
is  modified by the off-shell non-resonant interaction $J^\mu_{\pi N}$
to give the dressed vertex $\bar{J}^\mu_{\Delta N}$.

\section{Weak pion production reaction}

We apply the method described in the previous
section to  $\nu_\mu + N \rightarrow \mu + \pi + N$ reaction.
The low energy effective Lagrangian of the electroweak
standard model is given as
\begin{eqnarray}
 H_{W} & = & \frac{G_F \cos_c}{\sqrt{2}}(V^\mu - A^\mu)l_\mu .
\end{eqnarray}
Here $V^\mu$ and $A^\mu$ are charged currents of hadron and
$l^\mu$ is lepton current. The vector current($V^\mu$) is obtained from
the electromagnetic current in Refs. 1,2 by assuming
CVC and iso-spin rotation. The axial vector hadron current is
obtained from  the chiral Lagrangian.
Following the procedure described in the previous section,
the pion production currents $J^\mu_{eff}$ are calculated
shown in the schematic diagrams of Fig. 2.
We took the constituent quark model relation of
the  coupling constant
of the  axial $N\Delta$ current $A_{\Delta N}^\mu$  is obtained from the
quark model relation $G_A^* = \sqrt{72/25} g_A$.
The $Q^2$ dependence of the axial form factor is assumed as
\begin{eqnarray}
G_A^*(Q^2) = G_A^(0)\frac{1}{(1 + Q^2/m_A^2)^2}R_{SL}(Q^2),
\end{eqnarray}
where the dipole form factor
 is axial vector form factor of nucleon  with $m_A= 1.02GeV$\cite{meissner}
and phenomenological form factor $R_{SL}(Q^2)= (1 + a Q^2)exp(- b Q^2)$
with $a=0.154 GeV^{-2}$
  and $b=1.66 GeV^{-2}$ is obtained by fitting pion electroproduction
reaction cross section at $Q^2=2.8$ and $4(GeV/c)^2$.
All the other coupling constants and the form factors used in
this work are the same value as our model of the pion electroproduction.
Therefore, there is no adjustable parameters  in this calculation.

We first compare the total cross section of $\nu_\mu + p \rightarrow \mu^-
+ \pi^+ + p$ reaction with the data\cite{data79} in Fig. 3.
Our result shown in solid curve agrees reasonably
well with the data. One of the main feature of the dynamical approach is
the bare form factor of $N\Delta$ transition is renormalized by the pion
rescattering effects. The importance of this effect is seen
in Fig. 3.
The full result in solid line is reduced
into dashed line
when we turn off the dynamical pion cloud effect.
The contribution of the pure non-resonant contribution
is shown in dot-dashed line. We next compare the $Q^2$ dependence of the
calculated differential cross section with the data \cite{data90} in Fig. 4.
Our model reproduce very well the data. In the low $Q^2$ region, the
main contribution is the the axial vector current shown in the dashed
line. The BNL data were used in the most recent attempt to extract
$N\Delta$ axial vector form factor.

Finaly we study the effect of dynamical pion cloud effect on the axial
vector form factor. The empirical $N\Delta$ from factor extracted from
the data can be compared with our dressed form factor
$\bar{J}^\mu_{\Delta N}$ shown in the solid line in Fig. 5, while the
contribution of the bare form factor $J^\mu_{\Delta N}$ is shown in the
dotted line. We see sizable contribution of the dynamical pion cloud as
the case of  the electromagnetic form factors.
The differences between solid and
 dotted lines explain the observation that
the quark model prediction of axial $N\Delta$ from
 factor is about 30\% smaller than the empirically extracted form factor.

\begin{figure}
\centerline{\epsfig{file=stot1.ps,width=5cm}}
\vspace*{8pt}
\caption{Total cross section of $p + \nu_\mu \rightarrow \mu^- + \pi^+ + p$
    reaction.}
\end{figure}
\begin{figure}
\centerline{\epsfig{file=bnl-av.ps,width=5cm}}
\vspace*{8pt}
\caption{Differential cross section $d\sigma/dQ^2$ of
 $p + \nu_\mu \rightarrow \mu^- + \pi^+ + p$ reaction}
\end{figure}

\begin{figure}
\centerline{\epsfig{file=ga.ps,width=5cm}}
\vspace*{8pt}
\caption{The $N\Delta$ axial vector form factor.}
\end{figure}

\section{Summary}

In summary, the dynamical model developed in \cite{sl1,sl2,sl3} can
describe extensive data of pion photoproduction, electroproduction
and the neutrino-induced reactions. The predicted $\nu_\mu + N \rightarrow
\pi + N + \mu$ cross sections are in good agreement  with the existing
data. The renormalized axial vector $N\Delta$ form factor contains
large dynamical pion cloud effects and these effects are crucial in
getting agreement with the data. We conclude that the $N\Delta$ transition
axial vector form factor predicted from the quark model is consistent
with the existing data in the delta region. However more extensive and
precise data of neutrino-induced pion production reactions are needed to
further test our model and to pin down the $Q^2$ dependence of the
form factors. The efforts to extend  the current dynamical model
beyond the delta resonance region\cite{sl5} and to improve the effective
Hamiltonian including leading loop corrections\cite{sl4} are in progress.

\section*{Acknowledgments}

This work was supported by the
Japan Society for the Promotion of Science,
Grant-in-Aid for Scientific Research (C) 15540275 and
 the U.S. Department of Energy, Office of
Nuclear Physics Division, under contract no. W-31-109-ENG-38.


text:
\section{Introduction}

	Modern machine learning (ML) approaches based on deep neural networks have recently obtained impressive results in a range of classification tasks, sometimes even outperforming human experts. These successes are made possible by the combination of 1) better learning algorithms, 2) fast, massively parallel computing hardware including graphics processing units, and 3) the availability of large training data sets. However, in many application domains, such large data sets may simply not exist or be extremely expensive to gather. This problem is particularly severe in certain medical applications, where the numbers of patients may be quite small. Typical data sets may contain only hundreds or thousands of samples, while modern ML approaches often require the estimation of many millions of free parameters. Fitting a model with many free parameters to a small set of training samples will likely lead to over-fitting and poor generalization of the learned model. This problem is aggravated if the training data are corrupted by different kinds of noise, which is often unavoidable in biomedical data.

	\begin{figure}[tb]
		\centering
		\includegraphics[width=0.85\linewidth]{spectro1.JPG}
		\caption{Example MRS spectrum from a healthy patient}
		\label{spectrum}
	\end{figure}

	Here, we study the problem of brain tumor detection from magnetic resonance spectroscopy (MRS) data. In clinical practice, MRS is a common tool to identify a brain tumor and distinguish it from other medical conditions. It measures the resonant frequency shift of a chemically bounded hydrogen atom (i.e., a proton), which characterizes different physiological or pathological brain metabolites. There has been increasing interest in MRS for clinical use because of the semiautomatic data acquisition, processing and quantification \cite{kreis2004issues, ranjith2015machine, munteanu2015classification}. An example MRS spectrum from a non-tumor patient is shown in Fig.~\ref{spectrum}. While the interpretation of spectra is traditionally based on the size and location of certain peaks, we here use a novel approach by analysing the pattern of the MR spectrum as a whole in an unbiased fashion with machine learning.

	A common problem with in-vivo MRS data is that they are quite noisy. Noise sources range from heterogeneous magnetic susceptibilities of human tissues over baseline distortions of the spectrum \cite{kreis2004issues} to head movement during the procedure. Hence, the quality of spectra may be inadequate to determine precise metabolite concentrations and artefacts may resemble diagnostic features. As an additional problem, during the tissue selection process, due to the indefinable borders of gliomas, spectra from the tumor-affected hemisphere can be falsely labeled as tumor even though they contain healthy brain tissue.

	Furthermore, depending on the size of the selected region of interest, the number of samples collected from each patient varies substantially. Such a heterogeneous distribution of the individual training samples impedes the generalization of the learned model --- especially in a leave-one-out (LOO) cross validation scheme.

	\begin{figure*}[tb]
		\centering
		\includegraphics[width=0.95\linewidth]{overview.pdf}
		\caption{Overview of the proposed approach. MRS spectra from both classes are obtained. A data distillation network automatically selects representative samples that serve as the basis for creating an augmented data set. The augmented data set it used to train a final network for classification.}
		\label{workflow}
	\end{figure*}

	Scarcity of training data can be a big hurdle when applying DL methods to medical problems. Data augmentation is a common approach to alleviate this problem. It works by synthesizing new training data from the existing data via a variety of methods reviewed below. Here, we propose a new framework using two separate neural networks: a data distillation network to select representative training examples and a final classification network. In a nutshell, our method works by identifying data points that are ``easy'' to classify through the distillation network. Then, these data samples are used to synthesize a large number of new training data samples for training the final classifier. The new training samples are generated by mixing the easy samples with randomly selected data points from the same or other classes. The proposed framework is illustrated in Fig.~\ref{workflow}. Notably, it does not require human supervision to carefully label a small data set as prototype samples \cite{Li2017a, Lee2018}, but learns directly on the noisy labeled data. We show the benefits of this approach by demonstrating that it outperforms state-of-the-art methods and achieves human-expert-level performance. In sum, we make the following contributions:
	\begin{itemize}
		\item We propose a framework for tumor classification based on MRS data that combines deep neural networks with a novel data distillation and augmentation procedure to combat scarcity of the training data and labeling noise.
		\item We quantify the performance of human expert neuroradiologists on \textit{tumor}/\textit{healthy} classification from MRS data and demonstrate that our approach achieves human-expert-level performance.
		\item We show that the network uses prominent features in the data that are commonly used in clinical practice, but also considers features that have not yet received much attention by medical professionals, pointing out new directions for future research.
	\end{itemize}

	The remainder of the paper is organized as follows. In Sect.~\ref{s-relate}, we will briefly review state-of-the-art methods for data augmentation and dealing with noisy labels. In Sect.~\ref{s-data}, we present our data set and the data acquisition and preprocessing. In Sect.~\ref{s-method}, we describe the deep neural network architecture and our new data augmentation technique. Sect.~\ref{s-results} presents and discusses our results, showing that our network can achieve human expert-level performance on this task by using the proposed data augmentation approach. Section~\ref{s-conclusion} concludes the paper.

	\section{Related work}
	\label{s-relate}
	In this section, we provide a brief review of recent research using deep neural networks in medical applications. We focus on the field of oncology and the problems of noisy labels and scarce data.

	\subsection{Deep Neural Networks}
	In recent years, DNN-based methods have gained more and more popularity in the healthcare domain and achieved some impressive results \cite{lin2019super, park2019machine,esteva2017dermatologist,haenssle2018man,hannun2019cardiologist,podnar2019diagnosing,capper2018dna,pereira2016brain,munteanu2015classification, ranjith2015machine}.
	Among different network structures, Convolutional Neural Networks (CNNs) have gained great popularity \cite{simonyan2014very, szegedy2016rethinking, he2016deep}. They are inspired by the information processing mechanism of the visual systems of mammals where individual neurons respond to inputs in a restricted region of the visual field known as their receptive field. In comparison to fully connected neural networks, CNNs have a weight-sharing feature where neurons in different locations have identical receptive fields such that their responses can be calculated via a convolution operation. This design significantly reduces the number of trainable parameters and improves the generalization ability of the network.
	Ng \textit{et~al.} applied a Deep CNN to electrocardiography (ECG) data in heart disease classification and have achieved better performance than human cardiologists \cite{hannun2019cardiologist}. In \cite{haenssle2018man}, a deep CNN was trained on dermascopic melanoma detection and achieved above-dermatologist performance. In the field of oncology, machine learning methods have obtained promising results on problems such as tumor detection, tumor segmentation, tumor progression, etc. \cite{podnar2019diagnosing, lin2019super, capper2018dna, park2019machine, pereira2016brain, munteanu2015classification, ranjith2015machine}. For examples, Pereira \textit{et~al.}, applied a deep CNN for tumor segmentation from MRI data \cite{pereira2016brain}. Podnar \textit{et~al.} used a machine learning predictive model for the diagnosis of brain tumors from routine blood test results\cite{podnar2019diagnosing}. Machine learning methods applied to MRS data, such as in \cite{ranjith2015machine}, obtained good results in tumor grade classification according to the WHO tumor grade standard. However, learning from a larger cohort with multiple medical conditions only from MRS data has not yet been performed.

	\subsection{Learning from Noisy Labels}
	Noisy labels are ubiquitous in the real world. In this study, noisy labeling refers to observed labels that are incorrect, i.e., due to the labeling procedure, the label assigned to the instance does not represent the class membership. Noisy labels are posing a non-trivial problem in deep model learning when an increasing ability to fit noise is accompanied with deeper layers. Given the ubiquity and importance of coping with noisy labeling, many works have been devoted to combat  this problem \cite{Li2017a, Lee2018, han2018co, smyth2019training, VeitLearning}. One promising direction is to learn from a small set of clean labeled data and then use them to update the network \cite{han2018co, Li2017a, VeitLearning}. Another direction is to design models that could learn directly with noisy labels \cite{han2018co, smyth2019training, RolnickDeep}.
	In \cite{Li2017a}, an auxiliary model is trained with a small but clean data set, which was manually labeled by human experts. Then the knowledge obtained by the auxiliary model is guiding the learning of the primary model in the form of one part of the primary training loss being the imitation loss of the primary model to the auxiliary trained model.
	Lee \textit{et~al.} proposed a hybrid system, which requires a small set of representative seed instances with precise labels \cite{Lee2018}. Then, the automated noisy label detection is achieved with a deep CNN.
	Veit \textit{et~al.} proposed a semi-supervised learning framework for multilabel image classification that leverages small sets of clean labels in conjunction with large amounts of noisy labels \cite{ VeitLearning}. Small sets of clean labels facilitate the learning of the mapping between noisy and clean labels, which not only reflects the noisy patterns but also the labeling structure.
	Han \textit{et~al.} proposed a co-teaching framework where two DNNs were trained simultaneously \cite{han2018co} with the whole data set. The networks train each other using small-loss training instances, since they are likely to be clean annotations. The intuition is that since two neural networks are different and are equipped with different learning capabilities, during learning they may capture diverse features from the training samples and through the small-loss instances they are filtering out ``clean'' samples for each other.
	Smyth \textit{et~al.} proposed a DNN-based framework, Training-ValueNet, which evaluates the contribution of one sample to the whole learning process and then discards those that negatively contribute to the learning \cite{smyth2019training}.

	\subsection{Data augmentation}

	A number of techniques have been explored to alleviate the problem of small training data sets. Data augmentation is a very effective way to expand the existing training set with more and diverse data in order to improve the generalization ability and incorporate invariance. Usually, data augmentation methods are domain- and dataset-specific. The fundamental rule of data augmentation is that the meaning of the target samples should be maintained regardless of the augmentation methods applied. The trained model should be reliable enough to predict the same class even when the samples are perturbed. One common class of data augmentation methods especially applicable to image data is based on different data transformations such as cropping, rotating, flipping, shearing, etc. \cite{Krizhevsky2017, simonyan2014very}. Another class of methods is referred to as adversarial training where models are trained with generated adversarial samples \cite{Volpi2018, tran2017bayesian}. In \cite{Volpi2018}, the authors were concerned with the problem of generalizing learning from only one single source distribution to the unseen data domain. They augment the training set with generated adversarial samples. Tran \textit{et~al.} proposed a joint learning scheme where a Bayesian data generator is trained with existing training samples and continuously generates new training samples for further classification \cite{tran2017bayesian}. In \cite{perez2017effectiveness}, images in different styles are generated through a CycleGAN model and then used for further image classification.

	In another line of thinking, one augmentation strategy is blending two or more training samples to generate new ones, though still in its early stage \cite{summers2019improved}. Inoue \textit{et~al.} propose a data augmentation method by mixing randomly selected images from the training set \cite{inoue2018data}. Jaderberg \textit{et~al.} presented a framework for recognizing natural scene text \cite{Jaderberg2014}. In this work, a larger text corpus is generated with font rendering, creating and coloring with a background image-layer, a foreground image-layer, and an optional shadow image-layer. A natural data blending process is applied, where a random crop of an image from the training dataset is blended with each layer of the synthesized image. The three image layers are also blended together randomly to give a single output image. Summers \textit{et~al.} investigated various example-mixing methods in generating new samples and found that all mixing-based data augmentation methods resulted in an improvement of baseline performance \cite{summers2019improved} . In this work, the algorithm learned that mixing several samples of certain classes in a nonlinear way results in an improvement of the generalization ability of the learned model. However, data blending requires more delicate considerations compared to traditional data augmentation methods with various image transformations. Questions such as blending what together, how much of each component should be used, etc., need to be carefully addressed.

	\section{Dataset}
	\label{s-data}

	\paragraph{Data acquisition}
	1H-MR-spectroscopy data from 435 patients recorded in the Institute for Neuroradiology of the University Hospital, Frankfurt during the time interval from 01/2009 to 3/2019 were reviewed retrospectively. The spectroscopy was performed on a clinical 3T MR Scanner (Skyra, Siemens Medical Solutions, Erlangen, Germany) using a phased array head coil with 20 arrays with TE = 30 ms; TR~=~1500 ms; flip angle 90$^{\circ}$; scan time of 6:11 min. These patients were suffering from either glial or glioneuronal first diagnosed tumors (the \textit{tumor} group) or other non-neoplastic lesions e.g. demyelination, gliosis, focal cortical dysplasia, enlarged Virchow-Robin spaces or similar (the non-tumor/\textit{healthy} group). The tumor group included all spectra from the tumor-affected hemisphere. The non-tumor group consisted of all spectra from both hemispheres of the patients. As a result, 7442 spectra (3388 non-tumor and 4054 tumor) were selected for further analysis. The obtained MRS examples are saved as column vectors ($288\times1$), shown in Fig.~\ref{spectrum}, where the $y$-axis shows signal intensities of different metabolites, and the $x$-axis represents the chemical shift positions in ppm indicating the various metabolites.

	\section{Methods}
	\label{s-method}
	In this section, we formulate our problem of classifying MRS data collected from patients with and without brain tumors into \textit{tumor} and \textit{healthy} classes with deep neural networks. We first outline the challenges we face in this work and propose solutions. Then, we describe in detail the network structure and the corresponding parameters. To evaluate how well our proposed method works, we construct a performance comparison in a realistic clinical setting with eight neuroradiologists.
	\subsection{Challenges}
	In our particular problem, we face several challenges regarding the dataset, i.e., noisy labeling, data shortage and imbalanced classes.
	\paragraph{Noisy labeling}
	Infiltrative growth is an important feature of gliomas, which distinguishes them from expansively growing tumors, such as metastases. The real borders of gliomas are indefinable, which can strongly confound the selection and labeling of the voxels from multivoxel spectroscopy. One source of labeling noise is introduced when spectra from the tumor-affected hemisphere are falsely labeled as tumor-containing voxel although they contain healthy brain tissue.

	\paragraph{Data Shortage and Class Imbalance}
	A large amount of training data is one of the most essential factors in training DL models successfully. However, as mentioned before, the amount of such MRS data is limited by the number of patients with the medical conditions of interest. Furthermore, as the size of the selected region of interest varies substantially for each patient, the number of samples collected from each patient also varies. Such imbalance can negatively affect the training of a classifier.

	\subsection{Proposed Solutions}
	To tackle the problems mentioned above, we propose the following solutions.
	\paragraph{Data Distillation}\label{disillation}
	To deal with the noisy labeling problem, we propose to distillate the data before they are fed into the deep neural network. In this way, the deep neural network  will start by learning a subset of most representative samples from each class, i.e., ``clean-labeled'' samples \cite{arpit2017closer}. Specifically, we propose a distillation network to collect samples that this network is very certain about during initial training epochs. The intuition is that these samples are highly representative and most likely to be correctly labeled.
	The detailed data distillation procedure is defined in Algorithm~\ref{distillation}. In the algorithm, $\mathcal{D}=\{{x_1, \dots, x_N}\}$, is the training set with a total number of $N$ training samples. We first train a distillation network, which has the same structure and configuration as the final model. In principle, any type of network can be used for the distillation purpose. The softmax output from the DNN of the $j$-th sample $x_j$ is a two-element vector of classification probabilities over the two possible classes: \textit{healthy} and \textit{tumor}. With a threshold $\theta$, set empirically to $\theta=0.99$, we collect the certain samples where the maximum probability is greater or equal than $\theta$ at the end of each training epoch.

	{\centering
		\begin{minipage}{.55\linewidth}
			\begin{algorithm}[H]
				\KwIn{The training set $\mathcal{D}=\{{x_1, \dots, x_N}\}$}
				\KwOut{$\mathcal{C}$  \quad // collection of certain samples}
				Initialize a DNN\;
				$\mathcal{C} \gets \{\}$  \;
				\For{$E\gets 1$ \KwTo $max_E$}
				{
					$\mathcal{C}_{E} \gets \{\}$\;
					\For {$j\gets 1$ \KwTo $N$}{
						\If{$\max(P(j)) \geq \theta$}{
							$\mathcal{C}_{E}\gets \mathcal{C}_{E} \cup \{x_j\}$\;
						}
					}
					$\mathcal{C} \gets \mathcal{C} \cup \mathcal{C}_{E}$\;
				}
				\caption{Data Distillation}
				\label{distillation}
			\end{algorithm}
		\end{minipage}
		\par
	}

	Here, we add $x_j$ to the certain sample set $\mathcal{C} \subseteq \mathcal{D}$, when the maximum value in $P(j)$ is greater or equal to the threshold $\theta$. As mentioned that the ``clean'' samples are learned first at the initial training process, $max_E$ is the maximum training epoch that is still considered as in the initial training. The set $\mathcal{C}_E$, where $E \in \{{1, 2, \dots, max_E}$\}, is the collection of certain samples from epoch $E$. $\mathcal{C}$ is defined as the union of all $\mathcal{C}_E$: $ \mathcal{C} = \mathcal{C}_1 \cup \ldots \cup \mathcal{C}_E.$

	\paragraph{Data Augmentation}
	With an increasing number of distillation epochs, more and more samples are collected. However, with the maximum number of $\textit{E}$, the collected samples are still only a small part of the whole training set. To increase the number of training samples and improve the performance of our primary learning model, we propose a data augmentation method. Let $\mathcal{A}$ denote the augmented set of samples.

	A new augmented sample is created according to
	\begin{equation}\label{mix_factor}
	x^{\mathcal{A}}_i = (1 - \alpha)\cdot x^{\mathcal{C}}_j + \alpha \cdot x^{\mathcal{C}}_k,
	\end{equation}
	where $\alpha \in [0, 1]$ is the mixing weight with a value of 0.5 as default, $x^{\mathcal{A}}_i$ is the newly generated sample, $x^{\mathcal{C}}_j \in \mathcal{C}$ is the target sample which will be augmented, and $x^{\mathcal{C}}_k \in \mathcal{C}$ is the sample that is used to be mixed into the target sample. The total number of augmented samples in set $\mathcal{A}$ divided by that of the set $\mathcal{C}$ is termed the augmentation factor. We propose three different augmentation strategies: augment with the same class (aug-with-same), augment with the opposite class (aug-with-other) and augment with both classes (aug-with-both). Based on the choice of the augmentation strategy, $x^{\mathcal{C}}_k$ could be randomly selected from either class groups or both. During the augmentation, the label of $x^{\mathcal{A}}_i$ is defined as the label of $x^{\mathcal{C}}_j$. To deal with the class imbalance, we apply the method of oversampling the minority class described in \cite{Chawla2002}.

	\subsection{Deep Neural Network Structure}
	In our implementation, we apply the residual neural network proposed by He \textit{et~al.} as the backbone \cite{he2016deep}. Residual neural network feature skip-connections, which connect the input of one layer and the pre-activation of another layer skipping multiple layers in between. This structure is usually termed a residual block.  One block usually consists of multiple computational layers such as convolutional or dense layers with batch normalization \cite{ioffe2015batch}, drop-out \cite{srivastava2014dropout}, and a non-linear activation transformation \cite{maas2013rectifier}. The input to the residual block is split into two branches: the main branch with convolution or dense matrix multiplication, batch-normalization, drop-out and the other branch usually with the identity transformation or max-pooling. The combination of the outputs of these two branches is passed through a non-linear activation function as the input of the next block.
	\begin{table}[t]
		\caption{Proposed network structure. The \textbf{Config} column shows the configuration in convolutional and dense layers (filter size $32\times1$ and the number of filters, or the number of units in the dense layer). The number of filters is increased every other block by a factor of 2. Every other block subsamples its input by a factor of 2, indicated by the value of \textbf{Stride}. Here, the batch size at the first dimension is omitted in the output shape column. GAP: global average pooling.}
		\label{structure}
		\centering
		\begin{tabular}{lccc}
			\toprule
			\textbf{Name}     & \textbf{Config}      & \textbf{Stride} & \textbf{Output size}\\
			\midrule
			Conv & $\left[ \begin{array}{cc} 32\times1, & 16  \end{array}\right]$ & 1 &   [batch size, 288, 1, 16] \\
			\midrule
			ResBlock 1 & $\left[ \begin{array}{cc} 32\times1, & 16  \\ 32\times1, & 16  \end{array}\right]$   & 1 &   [batch size, 144, 1, 16]  \\
			\midrule
			ResBlock 2 & $\left[ \begin{array}{cc} 32\times1, & 16  \\ 32\times1, & 16 \end{array}\right]$   & 1 &   [batch size, 144, 1, 16]  \\
			\midrule
			ResBlock 3     & $\left[ \begin{array}{cc} 32\times1, & 32 \\ 32\times1, & 32 \end{array}\right]$ & 2 &  [batch size, 72, 1, 32]     \\
			\midrule
			ResBlock 4     & $\left[ \begin{array}{cc} 32\times1, & 32 \\ 32\times1, & 32 \end{array}\right]$  & 1 &   [batch size, 72, 1, 32]    \\
			\midrule
			ResBlock 5     & $\left[ \begin{array}{cc} 32\times1, & 64 \\ 32\times1, & 64 \end{array}\right]$  & 2 &   [batch size, 36, 1, 64]    \\
			\midrule
			ResBlock 6     & $\left[ \begin{array}{cc} 32\times1, & 64 \\ 32\times1, & 64 \end{array}\right]$  & 1 &   [batch size, 36, 1, 64]    \\
			\midrule
			ResBlock 7     & $\left[ \begin{array}{cc} 32\times1, & 128 \\ 32\times1, & 128 \end{array}\right]$  & 2 &   [batch size, 18, 1, 128]    \\
			\midrule
			ResBlock 8     & $\left[ \begin{array}{cc} 32\times1, & 128 \\ 32\times1, & 128 \end{array}\right]$  &1 &   [batch size, 18, 1, 128]    \\
			\midrule
			GAP          &      &    & [batch size, 128] \\
			\midrule
			Dense        & 2   &     &   [batch size, 2]  \\
			\bottomrule
		\end{tabular}
	\end{table}
	We implement a deep residual neural network with 8 residual blocks following the classic structure from \cite{he2016deep}, including 17 convolutional layers and skip connections. It is inspired by the network architecture in \cite{hannun2019cardiologist}. Each residual block consists of two convolutional layers with batch normalization, drop out and ReLU non-linear activation functions. The convolutional layers have a filter width of $32 \times 1$. Experimenting with different kernel sizes, 32 gives good performance. The number of filters increases by a factor of 2 in every other block starting from 16. There is a sub-sample layer of factor 2 in every other block occurring at the same time when increasing the number of filters. We apply a dropout rate of 0.55 in all blocks. A global average pooling (GAP) layer follows the last convolutional layer to provide further visualization, which is termed a class activation map (CAM) \cite{Zhou2016}. The GAP layer is followed by a soft-max layer, which outputs a probability distribution over the two possible classes. The detailed parameters of the network structure are shown in Table~\ref{structure}.

	\subsection{Visualization through Class Activation Maps}
	Modern DL techniques are often viewed as black-box methods, where the decision making process is difficult to understand for humans. It raises worrying questions and hinders the practical deployment of such techniques. Much effort has been devoted to develop explainable and interpretable DL approaches \cite{Kindermans2018, Yosinski2015a, SebastianOn, Zhou2016}.

	In our work, we apply a GAP layer to reduce the risk of over-fitting and provide further visualization of the network decision making processes. The GAP squashes the output of each feature map with the shape $h\times w\times d$ from the previous layer into one single value with the shape of $1\times 1\times d$ reducing the number of features by $h\times w$ fold. The output of the GAP layer is fed directly to the final classification layer. Intuitively, the GAP operation converts feature maps into weights that represent the ``importance'' of all feature maps, namely the CAMs. An added value of this method is that we can easily trace back the ``importance'' to the input space and visualize how much of each part of the input contributes to the final classification decision.

	\subsection{Quantifying Performance of Human Experts}
	To illustrate how well our proposed method works in comparison to routine clinical diagnostic, a classification task on the same test set is conducted for both the network and human neuroradiologists. Eight experts with different levels of experience in the 1H-MR spectroscopy (from resident to specialist of neuroradiology), were given 844 randomly selected spectra (around 105 per person). They were asked to classify each spectrum as originating from the tumor or from non-tumor tissue reviewing only the spectral lines. They were blinded to any additional information such as T2-weighted images or similar. The overall performance of neuroradiologists is regarded as a collective effort. Inter-rater reliability is not applicable here, since every radiologist received different subsets of the data to classify.

	\section{Results}
	\label{s-results}
	To evaluate performance, we use the receiver operating characteristic (ROC) curve, which is a gold standard to evaluate the discriminative ability of a classifier. It is constructed by varying the classification threshold and calculating the false positive rate (FPR), i.e., 1 - specificity (SPE) with $\text{FPR} = \frac{\text{FP}}{\text{FP} + \text{TN}}$ and the true positive rate (TPR), i.e., sensitivity (SEN) with $\text{TPR} = \frac{\text{TP}}{\text{TP} + \text{FN}}$. The area under the curve (AUC) is a scalar value between zero and one which characterizes the goodness of the classifier. We also compare our results with three baseline methods: a support vector machine (SVM), a random forest (RF), and a plain deep residual CNN (ResCNN) with the same network structure but without distillation and data augmentation. The SVM classifier is implemented with the $sklearn.svm.SVC$ function with default configurations. The RF is implemented with the $sklearn.RandomForestClassifier$ function with 500 trees in the forest and entropy as the measure of the quality of a split. The AUCs of the SVM and the RF are obtained without data distillation and data augmentation as a baseline comparison.

	\begin{table*}[t]
		\caption{Performance measures with default configurations. The performance of the neuroradiologists is computed on one randomly selected cross validation set. The notion ``partial'' in \textbf{Dataset} indicates the test set that was evaluated by the neuroradiologists, whereas ``whole'' refers to averaging across all ten cross validation sets. Results are given as mean $\pm$ standard deviation.}
		\centering
		\begin{tabular}{cccccc}
			\toprule
			&  \textbf{Dataset }& \textbf{SEN}  & \textbf{SPE} & \textbf{AUC} & \textbf{Patient-wise} \\
			&   &   &  &  &  \textbf{accuracy} \\
			\midrule
			Neuroradiologists& partial &  $0.54$  & $\textbf{0.88}$ & $0.68$ & $0.69$\\
			SVM & partial &  $0.66 \pm 0.00$  & $0.58 \pm 0.00$ & $0.62 \pm 0.00$ & $0.57 \pm 0.00$\\
			RF & partial &  $0.63 \pm 0.01$  & $0.60 \pm 0.01$ & $0.61 \pm 0.01$ & $0.63 \pm 0.01$\\
			Plain ResCNN & partial &  $\textbf{0.74} \pm\textbf{ 0.03}$  & $0.58 \pm 0.04$ & $0.68 \pm 0.02$ & $0.69 \pm 0.15$\\
			Proposed model &partial &  $0.73 \pm 0.20$  & $0.76 \pm 0.30$ & $\textbf{0.77} \pm \textbf{0.08}$ & $\textbf{0.76} \pm \textbf{0.15}$\\
			\midrule
			SVM &whole &  $0.69 \pm 0.08$  & $0.61 \pm 0.17$ & $0.63 \pm 0.05$ & $0.67 \pm 0.07$\\
			RF &whole &  $0.70 \pm 0.08$  & $0.60 \pm 0.15$ & $0.64 \pm 0.05$ & $0.70 \pm 0.09$\\
			Plain ResCNN & whole &  $0.70 \pm 0.03$  & $0.60 \pm 0.03$ & $0.69 \pm 0.04$ & $0.67 \pm 0.50$\\
			Proposed model & whole &  $\textbf{0.71} \pm \textbf{0.05}$  & $\textbf{0.69} \pm \textbf{0.07}$ & $\textbf{0.73} \pm \textbf{0.07}$ & $\textbf{0.73} \pm\textbf{0.30}$\\
			\bottomrule
		\end{tabular}
		\label{without}
	\end{table*}

	\subsection{Training procedure}
	The ability of the classifier to generalize to new previously unseen patients is of great clinical importance. Therefore, we apply a 10-fold leave-subjects-out cross validation scheme. To be specific, we divide the patient list into 10 sub-lists each with around 40 patients. In each cross validation set, we withhold the data from the patients of one sub-list, while we train and validate on the data from the other sub-lists. The patient-wise accuracy is computed in each leave-out test set. For each patient, the classification probability of all voxels are averaged to get the probabilities of each class. Then, the patient-wise diagnosis is obtained as the class that has the highest probability. The patient-wise accuracy is defined by the number of correct patient-wise diagnoses divided by the total number of patients in that set.
	We randomly select one cross validation set which consists of 844 spectra from 40 patients for the final test against human neuroradiologists.
	The network is trained with randomly initialized weights using the Adam optimizer with default parameters $\beta_1=0.9$ and $\beta_2=0.999$ and a mini-batch size of 32. If not specified, the default parameters for data distillation and data augmentation are as follows: the primary learning model is trained with the certain samples collected from the first epoch up to the fifth epoch ($E=5$) from the data distillation network; the training data is augmented with both classes five-fold with the mixing weight $\alpha=$0.5. With our network configuration, at the epoch $E=5$, the average number of collected certain samples is around 2000, providing a total of around 10,000 samples with five-fold augmentation.

	To get an average performance of the effect of the proposed distillation process, we train the whole framework twice, i.e., a distillation network, which collects certain samples and a primary classifier with proposed data augmentation on all 10 cross-validation sets, with different random seeds. The results are averaged across the 10 cross-validation sets as well as the two runs. The overall performance is reported in Table.~\ref{without}. It shows that our proposed method slightly outperforms the human neuroradiologists.

	\begin{figure*}[tb]
		\centering
		\includegraphics[width=0.95\linewidth]{clustering-landscape.pdf}
		\caption{K-means clustering on the whole data set. \textbf{A}. Cross-tab relation of the clustering results. \textbf{B-H}. Mean spectra of each cluster (dark blue) with standard deviation (cyan). }
		\label{clustering}
	\end{figure*}

	\begin{figure*}[t]
		\centering
		\includegraphics[width=0.98\linewidth]{distill-hist.pdf}
		\caption{Distillation effect visualization. A: 2-D t-SNE visualization of all training samples  (light blue and light pink) and distilled samples from both classes (blue and red). B: Euclidean distance histogram of all training samples and distilled samples from both classes. Dashed lines are the medians in each class with and without distillation. }
		\label{hist}
	\end{figure*}

	\subsection{K-Means Clustering of the Data}
	To get an overview of the data we use in this task, we performed k-means clustering on the whole data set $\mathcal{D}$, which has 7442 spectra (3388 healthy and 4054 tumor). The euclidean distance is used as the criterion to cluster the data. The number of clusters is determined by the elbow point in the inertia curve \cite{kodinariya2013review}, where the within cluster distance does not decrease significantly with an increasing number of clusters (a number of seven is chosen in this study). We did find a large overlap between two classes as we expected.

	The clustering results are shown in Fig.~\ref{clustering}. The cross-tab relation, which is a frequency count of one variable (\textit{healthy} or \textit{tumor}) in each cluster is shown in Fig.~\ref{clustering}. A. For example, cluster 1 contains 18.5\% of the healthy spectra and 8.7\% of the tumor spectra We can see that 1) there are samples from healthy and tumor group in every cluster, 2) there are roughly equal amounts of healthy and tumor samples in clusters 2, 3, 5, 6 and 7, 3) the majority of samples in cluster 1 are showing typical features of \textit{healthy} (no Lip1 or Lip2 concentration \cite{Fan2006Magnetic}) and those of cluster 4 are mainly typical \textit{tumor} (high Lip peaks, an elevated Cho peak, high Glx region, etc. \cite{RaeRE}), and 4) the majority of the spectra are neither typical healthy nor tumor, rather somewhere in between. The positions of typical metabolites are demonstrated in Fig.~\ref{CAM}-A. The mean spectra of those clusters illustrate commonly applied clinical assessment criteria:  in healthy tissues, there is a dominant peak at NAA and almost no mobile lipids to be detected since they are mostly confined to the membrane \cite{RaeRE}. In tumor tissues, there are elevated Cho and Lip peaks \cite{Fan2006Magnetic}. A median to high Cho peak with easily visible Cr peaks can contribute to the identification of a tumor \cite{Fan2006Magnetic}. The clustering results support our argument that the labeling process is noisy, so the spectra from both classes are largely mixed with each other.

	\subsection{Distillation Analysis}
	To better understand the effect of the distillation process, we investigate what samples are collected through the first distillation network. To simplify the visualization, we randomly select one training set and its corresponding distilled set after five epochs of training. First, we present a visualization of the whole training set through a 2-D projection via t-SNE \cite{maaten2008visualizing}, shown in Fig.~\ref{hist}A. Light blue and light pink dots represent original training samples from class 0 and class 1, respectively. Then, we highlight the samples collected by the distillation process with darker colors (blue and red for class 0 and class 1, respectively). One can observe that the overall distributions of the whole training samples from both classes largely overlap with each other. However, the distillation process drives the distributions of both classes further apart. For example, the collected \textit{healthy} data samples lie mostly on the left side of the 2-D projection and the collected \textit{tumor} samples reside away from the center of the projetion of the \textit{healthy} class. Second, we quantify the distribution of Euclidean distances of samples from both classes with respect to the center of the opposite class in the original data space. The histograms of the Euclidean distances are shown in Fig.~\ref{hist}B. Indeed, the distillation process collects samples that are further away from the other class's center.

	\subsection{Data Augmentation}
	In this section, we discuss different effects on learning resulting from different options including the mixing weight $\alpha$, the augmentation factor, the index of the last source epoch from which we collect the certain samples and the three augmentation strategies (aug-with-same, aug-with-other and aug-with-both).

	We measure the AUC of the ROC curve with different parameter options for different augmentation strategies. The results are averaged across all 10 cross-validation sets with two different initial distillation networks.

	Adding noise to augment data is a common practice in image data enrichment. Here, we also report results for the case when Gaussian noise is added to augment the data (noise augmentation). We explore different hyper-parameters such as noise amplitude and augmentation factor, and report the performance under the parameters that yielded the best result during the exploration.

	\begin{figure}[!tbp]
		\centering
		\begin{minipage}[b]{0.48\textwidth}
			\includegraphics[width=\textwidth]{auc_as_mix.pdf}
			\caption{The AUC as a function of the mixing weight $\alpha$ (from 0.05 to 0.95) for three different augmentation methods. The default augmenting factor is 5 and the index of the last source epoch is $E=5$. Error bars represent one standard deviation across all cross-validation sets. RF: random forest, SVM: support vector machine.}
			\label{auc_mix}
		\end{minipage}
		\hfill
		\begin{minipage}[b]{0.48\textwidth}
			\includegraphics[width=\textwidth]{auc_as_epoch.pdf}
			\caption{The AUC as a function of the source epoch index $E$ whose certain samples we used for data distillation for the three different augmentation strategies. The default augmenting factor is 5 and the mixing weight is $\alpha = 0.5$. Error bars represent one standard deviation across all cross-validation sets.}
			\label{auc_epoch}
		\end{minipage}
	\end{figure}

	In Fig.~\ref{auc_mix},	we show the AUC as a function of the mixing weight $\alpha$ for different augmentation scenarios. We note a number of observations. First, the plain ResCNN without any data augmentation already outperforms the simple SVM and RF classifiers (horizontal lines). Second, the common data augmentation technique of adding noise to the data points actually deteriorates performance compared to the plain ResCNN. Third, of the three proposed data distillation and augmentation approaches, the ``other'' method performs comparably to the plain ResCNN, sometimes giving slightly better or slightly worse results depending on the mixing weight. Fourth, best performance is achieved with the ``both'' and ``same'' methods, which show improved performance compared to the plain ResCNN over a range of intermediate mixing weights.

	In Fig.~\ref{auc_epoch}, we show the AUC as a function of the index of the source epoch $E$, when the collection of distilled samples stops (see Algorithm ~\ref{distillation}). Here, $E = 1, \dots, max_E$, where $max_E$ is set to be 10. The mixing weight $\alpha$ was 0.5 and the augmentation factor was 5. In this task,
	The AUC shows no clear preference among augmenting with ``same''and ``both'' cases, which are always superior to the plain ResCNN method. However, the average performance in the ``other'' case drops (with large standard deviation) as $E$ increases.

	\begin{figure}[!tbp]
		\centering
		\begin{minipage}[t]{0.48\textwidth}
			\includegraphics[width=\textwidth]{auc_as_factor.pdf}
			\caption{The AUC as a function of the augmentation factor in three different augmentation strategies during the test. The default source epoch index is 5 and the mixing weight is 0.5.  Error bars represent one standard deviation across all cross-validation sets.}
			\label{auc_fold}
		\end{minipage}
		\hfill
		\begin{minipage}[t]{0.48\textwidth}
			\includegraphics[width=\textwidth]{comparison.pdf}
			\caption{Comparison of the model and neuroradiologists on one randomly selected cross validation set. The individual and collective performance of neuroradiologists are shown as red dots and red star, respectively. The average ROC curve of the plain ResCNN model and our proposed model with default augmentation parameters (aug-with-both method, augmentation factor is five, and the mixing weight is $ \alpha = 0.5$) are depicted in dashed and solid blue lines, respectively.}
			\label{comparison}
		\end{minipage}
	\end{figure}

	The AUC as a function of the augmentation factor is shown in Fig.~\ref{auc_fold}. The index of the source epoch $E$ is five, and the mixing weight is $\alpha = 0.5$. It shows that with an increasing augmenting factor up to 3 the result shows a slight increase of the performance. Then, when it goes beyond five it shows a slight decrease in aug-with-\textit{same} and aug-with-\textit{other} augmentation scenarios.

	\subsection{Human vs. Machine}
	To assess how well our proposed method works in a more realistic clinical setting, we compared it to human neuroradiologists on one randomly selected test set. The result is shown in Fig.~\ref{comparison}.
	The performance of each individual neuroradiologist is denoted as a red dot, the collective performance is shown as a red star. The model without data augmentation has an AUC of 0.71 (dashed blue line). Our method achieves an AUC of 0.77 (solid blue), which encompasses most of the neuroradiologists in the ROC plot. It shows that our proposed method slightly outperforms the group of neuroradiologists as a whole (AUC: 0.77 vs. 0.68; sensitivity 0.73 vs. 0.54; accuracy: 0.76 vs. 0.69).

	\subsection{Feature Visualization}
	As described in section \ref{s-method}, we apply a GAP layer after the convolutional layers to prevent over-fitting and benefit from the possibility of visualizing class activation maps. These show how the network is making the final decision by assigning different weights, which can be interpreted as ``importance'', to different regions in the input data.

	\begin{figure}[th]
		\centering
		\includegraphics[width=0.85\linewidth]{CAMs_selected.pdf}
		\caption{Class activation maps of examples from both classes during test. At the top is an example where various metabolite peaks have been marked. Examples from class \textit{healthy} and \textit{tumor} are color coded in green and purple, respectively. Solid lines are original examples and the dashed line below is the corresponding CAM. Cr1, Cr2: creatine, Ins: myo-inositol. Cho: choline. NAA: n-acetylaspartic acid. Glx: glutamine, Lip: lipid}
		\label{CAM}
	\end{figure}

	In Fig.~\ref{CAM}, we show some examples of CAMs with original MRS samples from both classes. The results show that the CAMs vary with regard to specific samples. To interpret these CAMs, one must not only focus on the highest peak but rather the overall shape together with the original spectra. Note that the ``importance" does not reflect whether the signal intensity of the corresponding metabolite is high or low. The co-occurrence of high ``importance" regions provides insights in the CAM interpretation. We can see that the network considers various common metabolites during the classification. Interestingly, for the \textit{healthy} class the network also pays more attention to the plateau left of the dominant NAA peak. In the \textit{tumor} spectra, this part of the spectrum appears as a rising slope, and represents the oncometabolite 2-hydroxyglutatat \cite{Andronesi,Yen} as well as tumor associated metabolites like glutamine \cite{Li}. The Ins peak together with Cr2 and Lip regions are highly interesting. A high Ins peak with above-baseline Lip peaks highly suggests tumor presence and a low Ins concentration with almost no free lipids suggests the \textit{healthy} class \cite{Lipid}. In cases where a high ``importance" is assigned to the Cho region, the \textit{tumor} spectra show a high Cho peak flanked by other tumor-associated metabolite peaks (glycine, myo-inositol) \cite{naaVScholin,Hattingen2009}. On the other hand, the \textit{healthy} group shows a similar or smaller Cho peak as the Cr1 peak.

	\section{Conclusion}
	\label{s-conclusion}
	In this paper, we present a DNN-based framework, which achieves above human-level performance on a realistic clinical task of classifying tumor and non-tumor tissues based on MRS data. We construct an effective data cleaning and augmentation framework consisting of two steps: 1) a data distillation network to clean noisy labeled data, 2) a data augmentation process, which enlarges the data set acquired in the first step for training a primary neural network for the final classification. Due to its generality, this data augmentation method could be used in various other research domains. By exploring various configurations of the proposed data augmentation method, we further demonstrate that data augementation by mixing samples from both classes is more stable and yields better results. A deep residual neural network is used as the primary learning model and a global average pooling (GAP) layer at the end of all convolutional layers provides us with a visualization of how much each part of the input contributes to the final classification decision. Our proposed framework outperforms neuroradiologists on sensitivity and patient-wise diagnosis accuracy with an area under the ROC curves of 0.77. With an improved capability of coping with noisy labeling and the scarcity of the training data, we believe that the framework proposed in this work could improve clinical practice, ultimately leading to more effective and accurate diagnosis of brain tumors in patients.

	\section*{Acknowledgment}
	This work is supported by the China Scholarship Council (No. [2016]3100), the LOEWE Center for Personalized Translational Epilepsy Research (CePTER), and the Johanna Quandt Foundation. Special thanks to Charles Wilmot for inspiring discussions. Furthermore a particular appreciation goes to Marija Radović for her ideas on automating the data export.

	\bibliographystyle{unsrt}