\section{Introduction}

Large scale mappings of Faraday rotation measure (RM) of extragalactic linearly polarized radio sources have been achieved extensively in the decades (Taylor et al 2009; Oppermann et al. 2012), with which various analyses have been obtained to investigate structures of galactic as well as intergalactic magnetic fields (e.g., review by Akahori et al. 2018).

Local magnetic fields in the Solar vicinity have been also studied using these RM data as well as polarization observations of the Galactic radio emission (Mao et al. 2012; {Wolleben et al.} {2010}; Stil et al. 2011; Sun et al. 2015; Sofue and Nakanishi 2017; {Liu et al.} {2017}; {Van Eck et al.} {2017} Alves et al. 2018).

Synchrotron radio emission is a tool to measure the total strength of magnetic field on the assumption that the magnetic energy-density (pressure) is in equipartition with the thermal and cosmic ray energy densities (e.g., Sofue et al. 1986). This method requires information about the depth of emitting region in order to calculate the synchrotron emissivity per volume, as the intensity is an integration of the emissivity along the line of sight (LOS).

Rotation measure is an integration of the parallel component of magnetic field multiplied by thermal electron density along the line-of-sight (LOS). It is related not only to thermal (free-free) radio emission, but also to \HI column density through thermal electron fraction in the neutral interstellar medium (ISM).

Determination of the LOS depth is, therefore, a key to measure the magnetic strength from synchrotron emission and the parallel magnetic component from RM. The depth is also required to estimate the volume densities of \HI and thermal electrons from observed \HI and thermal radio intensities. Emission measure and \HI column density are useful to estimate the LOS depth, given a relation between the thermal and \HI gas densities is appropriately settled.

In this paper, correlation analyses are obtained among various radio astronomical observables (RM, \HI column density, thermal and synchrotron radio brightness) in order to determine physical quantities of the ISM such as the magnetic strength, gas densities, and  LOS depth. One of the major goal of the present hybrid analysis method will be to obtain whole-sky maps of the total strength and parallel component of the magnetic field in the local Galactic disk.

\begin{figure}
\begin{center}
\includegraphics[width=70mm]{fig1.eps}
\end{center}
\caption{All-sky maps of $RM$, $|RM|$ ($\uRM$)  (Taylor et al. 2009), $\Tsyn$ and $\Tff$ (mK) at 23 GHz (Gold et al 2011), and  $\NHI$ ($\uNHI$) (Kalberla et al. 2005). }
\label{mapData}
\end{figure}

\section{Observables and correlations}

\subsection{Data}

The observational data for the Faraday rotation were taken from all-sky RM survey by Taylor et al. (2009), \HI data from the Leiden-Argentine-Bonn (LAB) survey by Kalberla et al. (2005), synchrotron and free-free emissions at 23 GHz from the 7-years result of Wilkinson Microwave Anisotropy Probe (WMAP) project by Gold et al (2011).  Figure \ref{mapData} shows the employed map data for $RM$,  $|RM|$, \HI column density $\NHI$, thermal (free-free) radio brightness temperature $\Tff$ and synchrotron radio brightness $\Tsyn$, both at 23 GHz.

Figure \ref{All-Lat} shows plots of the same data in figure \ref{mapData} against the latitude $b$ and cosec $|b|$. The global similarity of the latitudinal variations indicates that the four observed quantities are deeply coupled with each other. It is remarkable that all the plots beautifully obey the cosec $|b|$ relation shown by the full lines. This fact indicates that these radio observed quantities are tightly coupled with the line-of-sight depth (LOS) through galactic disk composed of a plane parallel layer in the first approximation.

A more detailed inspection of the figures reveals that, besides the global common cosec $|b|$ property, there exist systematic differences in the latitudinal variations among the quantities. The rotation measure, $RM$, shows milder increase toward the galactic plane than the other quantities. The saturation of $|RM|$ near the galactic plane suggests that the magnetic field directions are reversing there. On the other hand, synchrotron intensity has much sharper peak at the plane, and shows similar variation to \HI intensity. Another remarkable property is the sharper increase of the thermal emission toward the galactic plane than \HI. This manifests stronger dependence of the thermal emission on the ISM density through the emission measure $EM \p \ne^2 L$ than that of the \HI column $\NHI \p \nHI L$.

\begin{figure*}
\begin{center}
\includegraphics[width=13cm]{fig2.eps}
\end{center}
\caption{$|RM| \ (\uRM)$, $\Tsyn$ (mK), $0.1\x \Tff$ (mK) and $10^{-25}\x \NHI (\uHI)$  plotted against latitude $b$. Points are reduced by 1 per 10. Note the beautiful dependence of the plotted quantities on cosec $|b|$ relation shown by the full lines, indicating that they are tightly coupled with the line-of-sight depth of a plane parallel layer.
 Right panel shows the same, but against cosec $|b|$. Open circles are averages of the plots in every $2\deg$ bin of absolute latitude. Bars show standard deviations (sd) of the original values, which are given only in upper sides in order to avoid logarithm of negative values for large sd.
}
\label{All-Lat}
\end{figure*}

The radio observables are related to the ISM quantities as follows. Faraday rotation measure $RM$ is related to the thermal electron density $\ne$ and line-of-sight (LOS) component of the magnetic field strength $\Bpara$ through
\be
\(RM\/\uRM \)\sim 0.81\(\la\ne\ra\/\cc\)\(\la \Bpara \ra\/\muG\) \( L\/\pc\).
\label{RM}
\ee
The emission measure $EM$ is rewritten by the volume density $\ne$ of thermal electrons as
\be
\(EM\/\uEM \) \sim \( \avne \/ \cc \)^2 \(L \/ \pc \).
\ee
The \textsc{HI} column density $\NHI$ is given by the \HI volume density $\nHI$ as
\be
\(\NHI\/{\rm cm^{-2}}\)
\sim 3.086\x 10^{18} \(\avnHI\/\unHI\) \(L\/\pc\).
\ee
Here, $\la ~ \ra$ denotes LOS average, and is defined and described later.
The synchrotron radio brightness $\Sigma_\nu$, as observed by the brightness temperature $\Tsyn$, is related to the volume emissivity $\ep$, frequency $\nu$, and $L$ as
\be
\Sigma_\nu={2k\Tsyn\/\lambda^2} \sim {L\/4\pi}{ d\ep\/d\nu}\sim {L\/4 \pi}{ \ep \/ \nu},
\label{Sigma}
\ee
which may be rewritten in a practical way as
\be
\(\Tsyn \/{\rm K}\) \sim 1.472\times 10^{11}\(\nu \/{\rm GHz}\)^{-3}
\(\ep\/{\rm erg\ cm^{-3}s^{-1}}\)  \(\L\/{\rm pc}\),
\label{Tsynbyep}
\ee
where $\lambda=c/\nu$ is the wavelength and $k$ is the Boltzmann constant.

We assume that the Galactic disk is composed of four horizontal layers (disks) of \textsc{HI} gas, thermal electrons, magnetic fields and cosmic rays, which have the same half thickness (scale height) $\zdisk$. This means that the LOS depth $L$ of the four quantities are equal. This assumption may not be good enough for the synchrotron emission that may originate from a thicker magnetic halo. However, it may be considered that the contribution of magnetic halo to $RM$ and synchrotron emission is much smaller than that of the disk because of weaker magnetic strength and electron density by an order of magnitude. The here used depth $L$ is an effective depth, and is related to the geometrical depth through a volume filling factor, as will be described later in detail.

Based on these considerations, we assume that the ISM quantities are smooth functions of the effective LOS depth $L$ for the first approximation, and an average of any quantity $f$ over $L$ satisfies the following relation,
\be
\la f \ra ={\int_0^L f dx\/\int_0^L dx}={\int_0^L f dx \/L}.
\label{f}
\ee
We also assume for any quantities $f$ and $g$
\be
\la f \ra \sim  \la f^2 \ra ^{1/2},
\label{ff}
\ee
and
\be
\la fg \ra \sim \la f \ra \la g \ra.
\label{fg}
\ee

\subsection{Correlation among Radio Observables}

\subsubsection{Free-Free to \textsc{HI} tight relation}

Figure \ref{ISM-ISM}(a) shows a plot of $\Tff$ against the square of $\NHI$. The straight line indicates a $\Tff \p \NHI^2$, and plots on the log-log space well obeys this proportionality. This relation indicates that the thermal electron density $\avne$ is approximately proportional to $\avnhi$, if $L$ is not strongly variable from point to point, which is indeed the case except for the high $\NHI$ region close to the galactic plane. This correlation will be used to estimate the electron density from \textsc{HI} column.

\begin{figure*}
\begin{center}
\includegraphics[width=16cm]{fig3.eps}
\end{center}
\caption{Correlation of $\Tff$ and $\NHI$ with a line of power index 2, $\Tsyn$ and $\NHI$ with a line of index 8/7.  Vertical scaling of the lines are arbitrary.
}
\label{ISM-ISM}

\begin{center}
\includegraphics[width=16cm]{fig4.eps}
\end{center}
\caption{Correlation of $|RM|$ with (a) HI, (b) FF and (c) Synchrotron emissions. }
\label{RM-ISM}
\end{figure*}

\subsubsection{Synchrotron to ISM relation}

Figure \ref{ISM-ISM}(b) is a plot of $\Tsyn$ against $\NHI$. High intensity region is approximately represented  by a power law of index 7/4, $\Tsyn \propto \NHI^{7/4}$, as expected from frozen-in magnetic field into the ISM and energy-density equipartition between the magnetic field, cosmic rays, and ISM  (see Appendix). Plot of $\Tsyn$ against $\Tff$ in (c) shows a similar relation, where $\Tsyn \propto \Tff^{7/8}$ is expected from the equipartition, because $\NHI \propto \Tff^{1/2}$. In both plots, the synchrotron emission tends to exceed the energy equipartition lines at low intensity regions (high latitudes). This yields larger uncertainty of the estimated magnetic strength at high latitudes.

\subsubsection{RM to ISM relation}

Figure \ref{RM-ISM} shows plots of the absolute RM values against (a) $\NHI$, (b) $\Tsyn$ and (c) $\Tff$. It is impressible that the plots are more scattered than those in figure \ref{ISM-ISM}. This is because the rotation measure is an integrated function of the magnetic field strength along the LOS including the reversal of field direction. This scattered characteristics of RM is useful to derive the spatial variation of the LOS field direction and strength $\Bpara$.

Although $|RM|$ is scattered against the other ISM observables in the whole-sky data, it may better be correlated in a narrower restricted region. Figure \ref{All_to_HI} shows an example of plots of $|RM|$, $\Tsyn$ and $\Tff$ against $\NHI$ in a small area in the 1st quadrant of the Galaxy at $30\deg\le l \le 50\deg$ and $b\ge 0\deg$. Shown by green circles are latitudes corresponding to individual $\NHI$ data points, indicating the tight dependence of $\NHI$ on the latitude through LOS depths. This figure demonstrates how $|RM|$ is tightly correlated to $\NHI$ in a restricted area.

\begin{figure}
\begin{center}
\includegraphics[width=7cm]{fig5.eps}
\end{center}
\caption{Correlations of $|RM|$, $\Tff$ and $\Tsyn$ with $\NHI$ in a restricted area on the sky ($30\deg\le l \le 50\deg$, $b\ge 0\deg$). Shown by green circles are latitudes  with the values indicated by the right axis}.
\label{All_to_HI}
\end{figure}

Figure \ref{RM-HI-Lat} shows plots of $RM$ against galactic latitude and \textsc{HI} column in the 1st quadrant of the Galaxy. Grey dots show all points, and blue and red dots represent those in northern and southern two small regions in the same quadrant at $30\deg\le l \le 60\deg$ and $b \ge 10\deg$ and at $30\deg\le l \le 60\deg$ and $b \le -10\deg$.

Absolute RM value increases toward the galactic plane, the sign of RM changes from negative to positive as the latitude increases, which indicates sudden reversal of the magnetic field direction. The bottom panel in the figure represents the same phenomenon in terms of the column density of \textsc{HI} gas.

Linear relation of RM with \textsc{HI} column is found in low $|RM|$ and $\NHI$ regions. However, the linearity is lost toward the galactic plane at $|b|\le \sim 10\deg$ with increasing $\NHI$. This represents decrease in the LOS component of magnetic strength $|\Bpara|$ toward the plane, which indicates rapid change of the field direction near the plane.

\begin{figure}
\begin{center}
\includegraphics[width=7cm]{fig6.eps}
\end{center}
\caption{RM against galactic latitude (top) and HI column density (bottom) in the 1st quadrant of the Galaxy. Grey circles show all data points, blue dots are for a small region at $30\deg\le l \le 60\deg$ and $b \ge 10\deg$ and red for $30\deg \le l \le 60\deg$ and $b \le -10\deg$  in the same quadrant.  }
\label{RM-HI-Lat}
\end{figure}

\section{Gas densities, Line-of-sight depth, and disk thickness}

Given the four radio observables ($RM$, $\Tff$, $\Tsyn$, and $\NHI$), four ISM parameters (averaged electron density $\avne\sim \xe \nHI$,  effective LOS depth $L$, total magnetic intensity $B$, and LOS component of magnetic field $\Bpara$) can be estimated as follows. Figure \ref{ISMhybrid} illustrates the flow of the analysis, which we call the ISM hybrid.   The effective depth $L$ will be related to geometrical scale height of the disk through volume filling factor.

\begin{figure}
\begin{center}
\includegraphics[width=5.5cm]{fig7.eps}
\end{center}
\caption{ISM hybrid for determination of physical quantities in the local disk from multiple  radio observables.  }
\label{ISMhybrid}
\end{figure}

\subsection{Thermal electron density}

Thermal electron density is assumed to be proportional to the \textsc{HI} gas density as
\be
\ne = \xe \nHI,
\ee
where, $\xe$ is the thermal electron fraction in the neutral ISM,
\be
\xe={\avne\/\avne+\avnhi}\sim {\avne \/\avnhi},
\ee
which is assumed to be $\xe\sim 0.1$ after Foster et al. (2013), who obtained $\xe\sim 0.08$.
 Local \textsc{HI} gas is considered to be in cold phase from recent measurement of spin temperature (Sofue 2017, 2018). Even if warm \textsc{HI} is contaminated, its density is an order of magnitude lower, so that the column density is not much affected by warm HI, unless the scale height of warm \textsc{HI} is an order of magnitude greater than that of cold HI.

Electron density is obtained by dividing the emission measure by column density of electrons, which is related to \textsc{HI} column, or $\ne \sim (\ne^2 L)/(\ne L) \sim (\ne^2L)/(\xe \nHI L)\sim EM/(\xe \NHI)$. Thus, we have
\begin{eqnarray}
\({\avne\/\cc}\)
\sim 3.09\x 10^{-3}\({EM\/\uEM}\) \(\xe\NHI\/\uNHI\)^{-1}\\ \nonumber
\sim 6.82\times 10^2 \xe^{-1} \(\NHI\/\uNHI\)^{-1} \(\Tff\/{\rm K}\)_{23{\rm GHz}}
\label{avne}
\end{eqnarray}
at $\nu=23$ GHz, where the following relations were used.
The emission measure is expressed by $\Tff$, electron temperature ($\Te \sim 10^4$ K), and observing frequency ($\nu=23$ GHz) as
\begin{eqnarray}
\(EM \/\rm pc \ cm^{-6}\)
= {3.05\x 10^2} \(\Tff\/K\) \( \Te \/10^4\K \)^{0.35}\( \nu \/\GHz \)^{2.1} \\ \nonumber
 \sim 2.21\x 10^5\({\Tff\/\K}\)_{23 \GHz}, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\label{EMinTb}
\end{eqnarray}
for $\Te=10^4$ K and $\nu=23$ GHz, where optical depth is given by (Oster 1961),
\be
\tau=3.28\x 10^{-7} \( \Te \/10^4\K \)^{-1.35}\( \nu \/\GHz \)^{-2.1}\(EM \/{\rm pc \ cm^{-6}}\),
\ee
which is related to $\Tff$ for optically thin case as
\be
\Tff =(1-e^{-\tau})\Te \simeq \tau \Te.
\ee

Figure \ref{mapDen}  shows a calculated all-sky map of $\nHI$, which is equal to $\ne/\xe$, using equation (\ref{avne}). The derived \textsc{HI} density at $|b|>\sim 10\deg$ has nearly a constant value around $\sim 1-2\ \cc$, except for clumpy regions and the GC.

\begin{figure}
\begin{center}
\includegraphics[width=85mm]{fig10.eps}
\end{center}
\caption{Logarithmic display of the volume density map of \textsc{HI} and thermal electrons. Shown is $\avnHI\ [\cc]$, which is assumed to be proportional to thermal electron density as $\avne \sim \xe \avnHI\ [\cc]$ with $\xe \sim 0.1$.  }
\label{mapDen}
\end{figure}

\subsection{LOS depth, scale height, and volume filling factor}

Recalling that $\NHI\sim L \nHI$ and $EM\sim L\ne^2 \sim L\xe^2\nHI^2$, the LOS depth $L$ is given by
\begin{eqnarray}
\({L\/\pc}\)
\sim 1.05 \x 10^5 \( \xe \NHI \/ \uNHI \)^2\(EM\/\uEM\)^{-1} \nonumber \\
\sim  0.475\({\xe\NHI \/ \uNHI}\)^2\({\Tff\/\K} \)^{-1}_{23 \GHz}.
\label{L_in_Tff}
\label{L}
\end{eqnarray}
Figure \ref{L-NHI} shows the derived $L$ plotted against $\NHI$ for $\xe=0.1$. The plot roughly obeys the linear relation, $\NHI \propto L$, indicated by the straight line, while points are largely scattered.

\begin{figure}
\begin{center}
 \includegraphics[width=7cm]{fig8.eps}
\end{center}
\caption{Effective LOS depth $L$ of HI gas for $\xe=0.1$ plotted against $\NHI$. The straight line indicates a linear relation, $\NHI \propto L$, with arbitrary vertical scaling.  }
\label{L-NHI}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=7cm]{fig9.eps}
\end{center}
\caption{"Effective" half thickness (scale height) of the \textsc{HI} disk defined by $\zhalf=L\ \sin \ |b|$ plotted against latitude, which is related to geometrical half thickness $\zdisk$ by $\zhalf \sim \eta^{1/3} \zdisk$ with $\eta$ being the volume filling factor of th gas. Circles are averages of the neighboring latitude points with standard errors by bars. Points without bar indicate errors greater than the averaged values.  }
\label{z-b}
\end{figure}

The here defined $L$ is an "effective (physical)" LOS depth, and is related to the effective half thickness (scale height) of the disk, $\zhalf$, by
\be
L=\zhalf \cosec \ |b|.
\label{L}
\label{Ltozdisk}
\ee
The effective half thickness is further related to the 'geometrical' half thickness
\be
\zhalf=\eta^{1/3}\zdisk,
\label{L}
\label{Ltozdisk}
\ee
where $\eta$ is the volume filling factor of the ISM. The factor $\eta$ will be determined using these relations referring to independent measurement of the \textsc{HI} disk scale height.

Figure \ref{z-b} shows calculated $\zhalf$ against latitude. Original data points are shown by gray dots, and averaged values in latitudinal interval of $\Delta b= \pm 5\deg$ are shown by big dots with standard errors. Points without error bar indicate those, whose errors are greater than the averaged values. The averaged effective half thickness tends to a constant of $\zhalf\sim 60\pm 8$ pc at $40\deg \le |b| \le 80\deg$.

From the current measurements of \HI disk half thickness
(160 pc, Lockman  (1984);
150 pc, Wouterloot et al. (1990);
200 pc, Levine et al. (2006);
173 pc, Kalberla et al. (2007);
200 pc, Nakanishi and Sofue (2016);
217 pc, Marasco et al. (2017)),
we adopt a simple average of the authors' values, $\zdisk=183\pm 26$ pc. In order for the present determination of $\zhalf$ to satisfy equation (\ref{L}), we obtain $\eta \sim (\zhalf/\zdisk)^3=0.035\pm 0.007$. This value agrees with the recent determination for cold HI gas by Fukui et al. (2018).

\section{Magnetic Fields}

\subsection{Parallel component}

The LOS (or parallel) component of the magnetic field can be obtained by dividing RM by the column density of thermal electrons as
\be
\({\avB \/\muG}\)
\sim 3.81\x 10^{-3}\({RM\/{\rm rad\ m^{-2}}}\) \(\xe\NHI\/\uNHI\)^{-1}.
\label{Bpara}
\ee
It is stressed that this formula yields $\avB$ directly from the observables $RM$ and $\NHI$ without employing the LOS depth $L$. So, $\avB$ is the most accurate quantity determined in this paper.

Equation (\ref{Bpara}) is particularly simple and useful to estimate the parallel component of magnetic strength, because it includes only two observables, $RM$ and $\NHI$, where the effective LOS depth $L$ has been canceled out, leaving $\xe$ as one parameter to be assumed. This relation is now applied for mapping of the $\Bpara$ value on the sky assuming $\xe=0.1$.

First, the sky is binned into $1\deg$ degree grids in longitude and latitude, and calculate averages of RM and $\NHI$ values within $\pm 1\deg$ about each grid point for $|b|\le 50\deg$ region, and within $\pm 2\deg$ for $|b|>50\deg$. At each point on the grids, $\avB$ is calculated with the aid of equation (\ref{Bpara}). By this procedure a $361\x 181$ meshed map of $\avB$ is obtained on the sky with a resolution of $1-2\deg$.

Figure \ref{mapBpara} shows the thus obtained all-sky map of $\avB$.
  Positive value in red color indicates a magnetic field away from the observer, and negative with blue indicates a field approaching the observer.

Let us remember that the RM map was strongly affected by the peaked line-of-sight depth near the galactic plane, causing large positive and negative values near the plane. This caused steep latitudinal gradient of RM due to the field reversal from north to south, resulting in RM singularity along the galactic plane.

On the other hand, the $\avB$ map is not affected by the LOS depth, so that it exhibits the field strength and direction only, so that the RM singularity along the galactic plane does not appear. The map reveals a widely extended arched region with positive magnetic strength of $\avB\sim +5 \muG$ in the north from $(l,b)\sim (40\deg,5\deg)$ to $(210\deg,0\deg)$. This arch seems to be continued by a negative strength arch with $\avB\sim -5\muG$ in the south from $(l,b)\sim (50\deg,-5\deg)$ to $(160\deg,-30\deg)$.

It may be possible to connect the positive and negative $\Bpara$ arches to draw a giant loop, or a shell, from $l\sim 40\deg$ to $220\deg$ with the field direction being reversed from north to south. Alternatively, the positive arch may be traced through the empty sky around the south pole in the present data (Taylor et al. 2009), where the improved map shows positive RM (Oppermann et al. 2012). If this is the case, the RM arches may trace a sinusoidal belt from the southern hemisphere in the 1st and 2nd quadrants to northern in the 3rd and 4th quadrants, drawing an $\infty$ shaped belt on the sky, with the necks in the galactic plane at $l\sim 30\deg$ and $240\deg$.
The arched magnetic region along the Aquila Rift from $(l,b)\sim (30\deg,-10\deg)$ to $(300\deg,+30\deg)$ with $\avB\sim +2$ to $-4 \muG$ could be a part of the $\infty$ belt.

 It is also interesting to note that both the northern and southern polar regions show positive $\Bpara$ with $\avB\sim +1 \ \muG$, indicating that the vertical (zenith) field directions are pointing away from the Sun.

\begin{figure*}
\begin{center}
\includegraphics[width=13cm]{fig11.eps}
\end{center}
\caption{All-sky maps of  $\avBpara$ with contours at interval of $1$ $\muG$.  Positive value with red color indicates a field away from the observer, and negative with blue, approaching. }
\label{mapBpara}
\end{figure*}

\begin{figure*}
\begin{center}
\includegraphics[width=13cm]{fig12.eps}
\end{center}
\caption{All-sky maps of  $\avBtot/Q$ with contours at interval of 2 $\muG$. The $\eta$-corrected strength is obtained by $\BQ\sim 0.72\times$ map values.}
\label{mapBtot}
\end{figure*}

\subsection{Total  intensity}

The total magnetic intensity $\Btot$ is calculated by assuming that the magnetic and cosmic ray energy densities are in equipartition as (see Appendix)
\be
B^2/8\pi \sim \Ncr \Ecr,
\label{Bequip}
\ee
where $\Ncr$ is the cosmic-ray electron number density and $\Ecr$ is representative energy of radio emitting cosmic rays. The magnetic strength is then related to the frequency $\nu$ and volume emissivity $\ep$ as
\be
\(\avBt\/ \muG \)
\sim 3.03 \x 10^8   Q \({\nu\/\GHz}\)^{-1/7}{\ep}^{2/7},
\label{B}
\ee
where, $Q$ is an equipartition factor, which depends on various assumed conditions and source models. There have been decades of discussion about $Q$ since Burbidge (1956), which includes dependence on such parameters as the spectral index, cut-off frequencies, proton-to-electron density ratio, volume filling factor, field orientation, and/or degree of alignment (e.g., Beck and Krause 2005).
The emissivity $\ep$ is related to $L$ and $\Tsyn$ through equation (\ref{Tsynbyep}),
\be
\(\ep\/{\rm erg\ cm^{-3}s^{-1}}\) \sim
6.796\times 10^{-12}\(\Tsyn \/{\rm K}\) \(\nu \/{\rm GHz}\)^3 \(\L\/{\rm pc}\)^{-1}.
\label{ep}
\ee
The emissivity depends on the filling factor through $\ep \propto L^{-1}\propto \eta^{-1/3}$. So, we here introduce an $\eta$-corrected magnetic strength,
\be
\BQ\sim \eta^{2/21}\avBtot/Q.
\ee
For $\eta\sim 0.035$ as measured in the previous section, we obtain $\BQ \sim 0.73 \avBtot/Q$.

Figure \ref{mapBtot} shows an all-sky map of the calculated total magnetic intensity $\avBtot/Q$. Except for discrete radio sources including radio spurs and GC, the map shows a smooth local magnetic intensity within $\sim 200$ pc.
 Local magnetic strengths were calculated in intermediate latitude regions at $+30\deg \le b \le +70\deg$ and $-70\deg \le b \le  -30\deg$ to obtain $\Btot/Q=10.7\pm 3.1\ \muG$ and $10.3 \pm 2.7 \muG$, respectively. Combining the two regions, we obtain $\Btot/Q=10.5\pm 3.0 \muG$.
By correcting for the volume filling factor ($\eta^{2/21}=0.73$), we obtain a representative local field strength of $\BQ=7.6\pm 2.1 \muG$ for $Q=1$.

Given $\avBtot$ and $\avBpara$ maps, the perpendicular component of the magnetic field is easily calculated by $\avBperp\sim \sqrt{\avBtot^2-\avBpara^2}$. However, the accuracy of the above estimated $\avBtot$ would be too poor to obtain a meaningful map of the perpendicular component.

\section{Discussion}

\subsection{Summary}

The latitudinal plots in figure \ref{All-Lat} indicate that the four observables, $RM$, $\Tsyn$, $\Tff$, and $\NHI$, are tightly correlated with each other through their cosec $|b|$ variations. This indicates that the distributions of the sources and their physical parameters are also tightly correlated with each other. Based on this fact, the sources of these emissions and Faraday rotation are assumed to be distributed in a single local disk in the Galaxy.

On this assumption, some useful relations were derived for calculating the local ISM quantities such as magnetic strength $\avBtot$, and LOS component of magnetic field $\avBpara$, thermal electron density $\avne$, \textsc{HI} density $\avnHI$, and LOS depth $L$, or the scale thickness $\zhalf$ and $\zdisk$ with the volume filling factor $\eta$. It was emphasized that determination of $L$ plays an essential role in the present hybrid method to calculate the physical quantities, while only $\avBpara$ can be directly calculated from $RM$ and $\NHI$ without being affected by $L$.

Applying the method to archival radio data, all-sky maps of $\avBpara$ and $\avBtot$ were obtained, which revealed a detailed magnetic structure in the local interstellar space within the Galactic disk near the Sun. The $\avBpara$ map showed that the magnetic direction varies sinusoidally along a giant arch-shaped belt on the sky, changing its LOS direction from north to south and vise versa every two galactic quadrants. Maximum parallel component of $\sim \pm 5-6 \muG$ was observed on the belt at intermediate latitudes. The $\avBtot$ map showed that the total magnetic strength is smoothly distributed on the sky, and the averaged value was obtained to be $\avBtot/Q \sim 10.5 \muG$ in the intermediate latitude region. Assuming an equipartition factor of $Q=1$, we obtained an $\eta$-corrected field strength of $\BQ\sim 7.6\muG$ for the measured volume-filling factor of $\eta\sim 0.035$.

\subsection{Dependence on the thermal electron fraction}

The proportionality of the densities of thermal electrons and \textsc{HI} gas is confirmed through the tight correlation between $\NHI$ and $\Tff\ (EM)$ by figures \ref{All-Lat} and \ref{ISM-ISM}. On this basis, we assumed a constant thermal electron fraction of $\xe\sim 0.1$ close to the current measurement on the order of $\sim 0.08$ (Foster et al. 2013). However, $\xe$ affects the result through equation (\ref{L}), where $L \p \xe^2$ and it propagates to the other quantities as $\avBpara \p \xe^{-1}$, $\avBtot \p \xe^{-4/7}$, $\avne \p \xe^{-1}$, and $\avnHI \p \xe^{-2}$. Namely, the ISM quantities are generally proportional inversely to $\xe$, with strongest effect on $\avnHI$ and weakest on $\avBtot$.

\subsection{Uncertainty from energy equipartition}

The most uncertain point in the present analysis is the estimation of magnetic strength from the energy equipartition of cosmic-ray electrons and magnetic field. Since the equipartition factor $Q$ is still open to discussion, the obtained total magnetic intensities should be taken only as a reference to see the relative distribution of the strength on the sky.

ALso, the single disk assumption for synchrotron and thermal components may break at high latitudes. As in figure \ref{ISM-ISM}, the plots of $\Tsyn$ against $\NHI$ and $\Tff$ bend at high latitudes, showing an order of magnitude excess at high latitudes over smooth extension from the disk component. The synchrotron excess over that expected from frozen-in assumption is about $\delta \Tsyn \sim 0.1$ K.

In figure \ref{syn-ism-bgrm} we plot $\Tsyn-\delta \Tsyn$ mimicking {\it halo-subtracted} synchrotron emission, which is well fitted by a power law expected from low and intermediate latitude regions. This fact suggests that the energy-equipartition holds inside the disk, whereas a non-thermal halo at $\sim 0.1$ mK level at 23 GHz is extending outside the gas disk.

\begin{figure}
\begin{center}
\includegraphics[width=6cm]{fig13.eps}
\end{center}
\caption{Correlation of  $\Tsyn-0.1$ mK with $\NHI$ and $\Tff$, which are approximately represented by a power law of index 7/4 and 7/8, respectively. Vertical scaling of the lines are arbitrary. Original data are shown in violet. }
\label{syn-ism-bgrm}
\end{figure}

\subsection{Effect of inhomogeneity}

From the tight cosec $|b|$ relation of the used radio observables, for which no extinction problem exists, we assumed a uniform layered disk of ISM. In more realistic conditions, however, the layer may be more or less not uniform, and the assumption made in equations (\ref{f}-\ref{fg}) may not hold, or must be modified.

However, it is emphasized that 'clumpy' inhomogeneity does not affect the 'effective' LOS depth $L$ by definition, because $L$ already includes the volume filling factor. Hence, the determined values of the ISM, which are averages of values 'inside' the clumps (or within $L$), are not affected by the inhomogeneity.

The assumption made for equations \ref{f} to \ref{fg} will not hold exactly in a disk with globally varying density with the height. For example, if the functions $f$ and $g$ are represented by a Gaussian function of the height from the galactic plane, we have $\la f \ra \sim  1.1 \la f^2 \ra ^{1/2}$ and $\la fg \ra \sim 1.2 \la f \ra \la g \ra$. For a cosh$^{-2}$ function, as for self-gravitating disk, the factors are 1.2 and 1.4, respectively. These factors propagate onto the results, yielding uncertainty by a factor of $\sim 1.1-1.2$ for quantities having linear dependence on the distance, and by $1.2-1.4$ to those with non-linear dependence such as $RM$ and $EM$. However, the finally determined $B$ and $\ne$ or $\nHI$ are more linearly dependent on the distance, and hence their uncertainties may be about a factor of $\sim 1.1-1.2$ at most.

\subsection{Local bubble}

\def\rbub{r_{\rm bub}}
\def\tan{{\rm tan}}

A large-scale irregularity of the ISM has been reported as a local bubble ({Bochkarev} {1992}; Lallement et al. 2003; Liu et al. 2017; Alves et al. 2018), which makes a cavity around the Sun of radius $\rbub \sim 100-200$ pc widely open to the galactic halo. Then, a difficulty is encountered to explain the tight cosec $|b|$ relation in figure \ref{All-Lat}. In order for the cosec relation to hold up to $b\sim \pm 80\deg$ at least, the scale height of the galactic disk must be greater than $\sim \rbub \tan \ 80\deg \sim $ 600 pc to 1.2 kpc, which is obviously not the case. If the disk scale height is $\sim 200$ pc as measured in HI, the open cavity should result in huge empty sky in radio around the galactic poles, which also appears not the case.

In figure \ref{bubble} we compare the cosec $|b|$ relation observed for $\NHI$ with cavity models mimicking the local bubble. Line A indicates a Gaussian disk of scale height 200 pc without bubble; B represents a case with a spherical bubble of radius 100 pc in the Gaussian disk of scale height 200 pc, and C and D for cylindrical cavity of radius 100 and 200 pc, respectively. Model C may be compared with the result by Lallement et al. (2003), which appears significantly displaced from the cosec $|b|$ relation. Such is found not only in HI, but also in thermal and synchrotron emissions, and Faraday RM (figure \ref{All-Lat}). Therefore, the relation between the local bubble and the cosec $|b|$ disk in radio remains as a question.

\begin{figure}
\begin{center}
\includegraphics[width=7cm]{fig14.eps}
\end{center}
\caption{  Cosec $|b|$ relation of the \textsc{HI} disk from figure \ref{All-Lat} (averaged values every $2\deg$ by circles) compared with (A) a model of Gaussian disk of scale height 200 pc without bubble; (B) same disk with a spherical bubble of radius 100 pc; (C) and (D) disk with vertical cylindrical cavity of radius 100 and 200 pc, respectively. Values are normalized to the observation at $b\sim 5\deg$. }
\label{bubble}
\end{figure}

\subsection{Other observables}

Molecular gas has not been taken into account in this study, because the nearest molecular clouds within LOS depths concerned in this paper are rather few (Knude and Hog 1998). Comparison with a local bubble surrounded by dusty clouds (e.g., Lallement et al. 2003), besides the cosec $|b|$ problem, would be an interesting subject, although the present analysis gives only averaged values along the LOS within $L$, and hence cannot be directly compared with the 3D study.

Polarization data in radio and infrared observations were not used, although they are obviously useful to improve the present hybrid analysis. Inclusion of these observables is beyond the scope of this paper, for which more sophisticated analyses would be required.

\subsection{Local magnetic topology}

Despite of the various uncertainties as above, we emphasize that the projected topology of $\Bpara$ mapped in figure \ref{mapBpara} is rather certain. Although the {\Large $\propto$}-shaped variation of $RM$ might sound a bit strange, we could speculate possible topology of the magnetic lines of force in the local space.

The field direction reverses about the galactic plane from north to south in a wide range from $l\sim  30\deg$ to $\sim 210\deg$. The $\Bpara$ value attains its maximum and minimum at both intermediate latitudes around $(l,b)\sim (130\deg,\pm 30\sim 50\deg)$. Such $RM$ behavior on the sky could be explained by a reversed topology of local field as illustrated in figure \ref{magtopo}.

\begin{figure}
\begin{center}
\includegraphics[width=7cm]{fig15.eps}
\end{center}
\caption{ Possible topology of local magnetic lines of force to explain $\Bpara$ map in figure \ref{mapBpara}.  }
\label{magtopo}
\end{figure}

\section*{Acknowledgments}
We thank the authors of the LAB \textsc{HI} survey (Dr. Kalberla et al. ), all-sky rotation measure map (Dr. Taylor et al. ), and the WMAP 7 years maps (Dr. Gold et al.) for the archival data. The data analyses were performed on a computer system at the Astronomical Data Center of the National Astronomical Observatories of Japan.


\section{Introduction}
Nowadays it is generally accepted that a huge number of real processes arising in physics, biology, chemistry, etc. can be described by nonlinear PDEs. And the most powerful methods for costruction of exact solutions for a wide ranges of nonlinear PDEs are symmetry-based methods, and these methods originated from the Lie method, so the basic part of the theory, is infinitesimal method of Sophus Lie~\cite{Lie}, that is connection between continuous transformation groups and algebras of their infinitesimal generators. This method leads to techniques in the group-invariant solutions and conservation laws of differential  equations~\cite{O1,Ib3,Ov1}.

In fact the method that we proposed, the method of preliminary group classification, is a conclusion of Lie infinitesimal  method, and is defined  and related to the theory of group classification of differential equations. this method is proposed in~\cite{Akhatov} and is developed for deferential equation in~\cite{Ib2,Bihlo}.

The main idea of preliminary group classification is based on extension of the kernel of admitted Lie groups that
are obtained by the transformations from the corresponding equivalence Lie group. The problem of finding inequivalent cases of such extension of symmetry can reduce to the classification of inequivalent subgroups of the equivalence Lie group(In particular, if a Lie group is finite-parameter, then one can use an optimal systems of its subgroups). we use equivalence transformations and the theory of classification of finite-dimensional Lie algebras. In this paper we study point symmetry and equivalence classification of $\mathrm{HESI}$ equation, leading to a number of new interesting nonlinear invariant models associated to non-trivial invariance algebras.

A complete list of these models are given for a finite-dimensional equivalence algebra derived for $\mathrm{HESI}$ equation. To obtain these goals we perform algorithms that is explained in references ~\cite{Ov1,Bila,Ib4,Cherniha}, and we use  similar works in ~\cite{Karn,Ib2,mah1,Dr1,Dr2,Ib1,Song}.

\textbf{For the local solution:} The existence of  $\textit{C}^{\infty}$ local solutions of $\mathrm{HESI}$ equation In $\mathbb{R}^3$ is studied in~\cite{Tian}, and the solution is in the following form:
\begin{equation*}
u(x,y,z)=\frac{1}{2}(\tau_1x^2+\tau_2y^2+\tau_3z^2)+\varepsilon^5\omega(\varepsilon^{-2}(x,y,z)),
\end{equation*}
where $\varepsilon$ and $\tau_i$  are arbitrary constants, and $\omega$ is a given smooth function. So at the end of this work, with the symmetry group of the equation, this solutions  transforms to another solution of $\mathrm{HESI}$ equation.

\textbf{About $\mathrm{HESI}$ equation:} Based on refrences~\cite{fros,Tian}, $k$-Hessian equations are a family of PDEs in $n$-dimensional space equations that can be written as ${\cal S}_k[u]=f$, where $1\leqslant k \leqslant n$, ${\cal S}_k[u]=\sigma_k(\lambda({\cal D}^2u))$, and $\lambda({\cal D}^2u)=(\lambda_1,\cdots,\lambda_n)$, are the eigenvalues of the \textit{Hessian matrix} ${\cal D}^2u$ $\big((\partial_i \partial_ju)_{1\leqslant i,j \leqslant n}\big)$ and  $\sigma_k(\lambda)=\sum_{i_1<\cdots<i_k}\lambda_{i_1}\cdots\lambda_{i_k}$, is a $k$th elementry symmetric polynomial.

The  $k$-Hessian equations include the Laplace equation, when $k=1$, And the Monge-Ampere equation, when $k=n$.

Here we study $2-$hessian equation in three dimensions, and f is an arbitrary function of $x,y,z$ ($\mathrm{HESI}$ equation) :
\begin{equation}
{\cal S}_2[u] := u_{xx}u_{yy}+u_{xx}u_{yy}+u_{yy}u_{zz}-u_{xy}^2-u_{yz}^2-u_{xz}^2,
\end{equation}
this equation is a fully nonlinear elliptic partial diferential equation, that is related to intrinsic curvature for three-dimensional manifolds.

In fact, the $2-$hessian equation is unfamiliar outside Riemannian geometry and elliptic regularity theory, that is closely related to the scalar curvature operator, which provides an intrinstic curvature for a three-dimensional manifold. Geometric PDEs have been used widely in image analysis ~\cite{sapiro}. In particular, the Monge-Ampere equation in the context of optimal transportation has been used in three dimensional volume-based image registration ~\cite{haker}.

The $2-$hessian operator  also appears in conformal mapping problems. Conformal surface mapping has been used for two-dimmensional image registration ~\cite{angenent,Gu}, but does not generalize directly to three dimensions. Quasi-conformal maps have been used in three dimensions ~\cite{wang,zeng}. However, these methods are still being developed.
\section{Principal Lie Algebra}
The symmetry approach to the classification of admissible partial differential equations depends heavily on a useful way of describing transformation groups that keep invariant the form of a given partial differential equation. This is done via the well-known infinitesimal method developed by Sophus Lie~\cite{O1,O2,Ov1}.Given a partial differential equation, the problem of finding its maximal (in some sense) Lie invariance group reduces to solving \textit{determining equations} that is an over-determined system of linear partial differential equation.
We consider the $2-$hessian equation as the form:
 \begin{equation}
\mathrm{HESI}:\; {\cal S}_2[u]=f(x,y,z),\label{Hesi}
 \end{equation}
where $u=u(x,y,z)$ is dependent variable and $x,y,z$ are independent variables and $f$ is arbitrary function. Considering the total space $ E=X\times U$ with local coordinate $(x,y,z,u)$ which $x,y,z\in X$ and $u\in U$. The solution space of equation (\ref{Hesi}) is a subvariety $S_{\Delta}\subset J^2(\mathbb{R}^3,\mathbb{R})$ of the second order of jet bundle of $3-$dimensional submanifolds of $E$. The $1-$parameter Lie group of infinitesimal transformations on $E$ is as follows:
\begin{equation}
\begin{aligned}
\tilde{x}&=x+t\xi(x,y,z,u)+O(t^2), & \tilde{y}&=y+t\zeta(x,y,z,u)+O(t^2), \\
\tilde{z}&=z+t\eta(x,y,z,u)+O(t^2), & \tilde{u}&=u+t\phi(x,y,z,u)+O(t^2),
\end{aligned}
\end{equation}
where $t$ is the group parameters and $\xi,\zeta,\phi$ and $\eta$ are the infinitesimals of the transformations for the independent and dependent variables, resp. So the corresponding infinitesimal generators have the following form generally
\begin{equation}
V=\xi(x,y,z,u)\partial_x +\zeta\partial_y +\eta\partial_z +\phi\partial_u.
\end{equation}
So based on Theorem 2.31 of~\cite{O2}, $V$ is a invariant point transformation if ${\bf pr}^{(2)}V[\mathrm{HESI}]=0$. Where ${\bf pr}^{(2)}V$ is the second order prolongation of the vector field $V$, that means:
\begin{equation}
 {\bf pr}^{(2)}V = V+{\phi}^x(x,y,z,u^{(2)})\,{\partial}_{u_x}+

 \cdots+{\phi}^{xz}(x,y,z,u^{(2)})\,{\partial}_{u_{xz}},
\end{equation}
 in which $u^{(2)}=(u,u_x,u_y,u_z,u_{xx},u_{xy},u_{xz},u_{yy},u_{yz},u_{zz})$ and
 \begin{align}
 {\phi}^J(x,y,z,u^{(2)})=\textbf{D}_J\bigg(\phi-\sum_{i=1}^{3}{\xi}^i\frac{\partial u}{\partial x_i}\bigg)+\sum_{i=1}^{3}{\xi}^i\frac{\partial {u_J}}{\partial x_i}
 \end{align}
 where $({\xi}^1,{\xi}^2,{\xi}^3)=(\xi,\zeta,\eta)$ and $(x_1,x_2,x_3)=(x,y,z)$, further $J=(j_1,\cdots,j_k)$ is a $k$-th order multi-index, and $j_i$s adopt $x$ or $y$ or $z$, for each $1\leqslant i \leqslant k$, then $\textbf{D}_J$ denotes the total derivatives for the multi-index $J$, and the $J$-th total derivative is as
 $\textbf{D}_J=\textbf{D}_{j_1}\cdots\textbf{D}_{j_k}$, where
\begin{equation*}
\textbf{D}_i={\partial}_{x_i}+\sum_{J}\frac{\partial {u_J}}{\partial x_i}\,{\partial}_{u_J}, \qquad (x_1,x_2,x_3)=(x,y,z).
\end{equation*}
So ${\bf pr}^{(2)}V$ acts on Eq.\eqref{Hesi} and with replacing $u_{yy}$ with equivalent expression of $\mathrm{HESI}$ equation we have the following system as determining equation:
\begin{eqnarray}
\begin{aligned}
&{\phi}_{xx}= {\phi}_{xy}= {\phi}_{xz}={\phi}_{xu}= {\phi}_{yy}= {\phi}_{yz}={\phi}_{yu}= {\phi}_{zz}= {\phi}_{zu}=0,\\
&{\phi}_{uu}= {\xi}_{u}= {\xi}_{yy}={\xi}_{yz}= {\xi}_{zz}= {\zeta}_{u}={\zeta}_{zz}= {\eta}_{u}= {\eta}_{zz}=0,\\
&{\xi}_{x}={\eta}_{z}, \;\; {\zeta}_{x}=-{\xi}_{y}, \;\; {\zeta}_{y}={\eta}_{z},\;\;{\eta}_{x}=-{\xi}_{z}, \;\; {\eta}_{y}=-{\zeta}_{z}, \\
&f_x{\xi}+f_y{\zeta}+f_z{\eta}+2f(2{\eta}_z-{\phi}_u)=0.
\end{aligned}
\end{eqnarray}
where $f$ is arbitrary function.with solving above relations we have:
\begin{eqnarray} \label{con}
\begin{aligned}
&\xi=c_6x+c_7y+c_8z+c_9, \qquad \zeta=c_{10}z+c_6y-c_7x+c_{11},\\
&\eta=-c_{10}y+c_6z-c_8x+c_{12}, \qquad \phi=c_1x+c_2u+c_3y+c_4z+c_5, \\
&f_x{\xi}+f_y{\zeta}+f_z{\eta}+2f(2{\eta}_z-f{\phi}_u)=0.
\end{aligned}
\end{eqnarray}
which $c_i,i=1,\cdots,12$ are arbitrary costants.

So if $f(x,y,z)=0$ the last equation of relations \eqref{con} will be removed and based on the first four equations in relations \eqref{con}, we have $12$-dimensional symmetry group, but if $f(x,y,z)\neq0$ we substitute the first four equations of \eqref{con} in the last one, and obtain the following \textit{condition}:
\begin{eqnarray} \label{con2}
\begin{aligned}
(c_6x+c_7y&+c_8z+c_9)f_x+( c_{10}z+c_6y-c_7x+c_{11})f_y \\
&+(-c_{10}y+c_6z-c_8x+c_{12})f_z+(-2c_2+4c_6)f=0.
\end{aligned}
\end{eqnarray}
There isn't $c_1,c_3,c_4,c_5$ in condition \eqref{con2}, so these coeficients are free,that means equation \eqref{Hesi} have $4-$dimentional symmetry group at minimum. So we conclude the following theorem from above relations:
\begin{theorem}:
The $\mathrm{HESI}$ equation(Eq.\eqref{Hesi}) admits symmetry group of dimension 4 to 12, for different choises of given function $f(x,y,z)$. These equations have the common following vectors as infinitesimal generators:
\begin{align}\label{principle}
V_1&=\partial_u , & V_2&=x\partial_u , & V_3&=y\partial_u , & V_4&=z\partial_u.
\end{align}
\end{theorem}

Then the Lie algebra $\mathfrak{g}$ generated with the vectors \eqref{principle} is called the principal Lie algebra for Eq.\eqref{Hesi}. Now we want to specify the coefficient $f$ such that Eq.\eqref{Hesi} admits an extension of the principal algebra $\mathfrak{g}$. therefore, we do not solve the determining equation, instead we obtain a partial group classification of Eq.\eqref{Hesi} via so-called method of  \textit{preliminary group classification}.

This method was suggested in ~\cite{Akhatov} and applied when an equivalence group is generated by a finite-dimensional Lie algebra $\mathfrak{g}_{\mathscr{E}}$. The essential part of the method is the classification of all nonsimilar subalgebras of $\mathfrak{g}_{\mathscr{E}}$. Actually the classification is based on finite-dimentional equivalence algebra $\mathfrak{g}_{\mathscr{E}}$.
\section{Equivalence Transformations}

with a nondegenerate change of the variables $x,y,z$ an equation of the form $\mathrm{HESI}$ equation convert to an equation of the same form, but with different $f(x,y,z)$. The set of all equivalence transformatioms forms an equivalence group $E$. We shal find a subgroup $E_c$ of it with infinitesimal method.

We suppose an operator of the group $E_c$ is in the form:
 \begin{equation}
 Y=\xi(x,y,z,u)\partial_x +\zeta\partial_y +\eta\partial_z +\phi\partial_u +\psi(x,y,z,u,f)\partial_f.\label{Y}
 \end{equation}
So from the invariance conditions of Eq.\eqref{Hesi} written as the following system:
\begin{align}
{\cal S}_2[u]=f(x,y,z),\qquad f_u=0.
\end{align}

Note that $f$ and $u$ are considered as differential variables; $u$ on the space $(x,y,z)$ and $f$ on the space $(x,y,z,u)$. The coordinates $\xi,\zeta,\eta,\phi$ of operator \eqref{Y} are funtions of $x,y,z,u$, while the coordinate $\psi$ is function of $x,y,z,u,f$. as usual way we should solve the following system that obtained of the invariance conditions:
\begin{align}
{\bf pr}^{(2)}Y( {\cal S}_2[u]=f(x,y,z)),\qquad {\bf pr}^{(2)}Y(f_u)=0.
\end{align}
Where ${\bf pr}^{(2)}Y$ is the second order prolongation of the vector field $Y$.

 \textbf{But}, to obtain the operator $Y$ of the group $E_c$ we use of \textit{N. Bila's method} in ref.~\cite{Bila}. The base of our procedure is theorem (1) of ~\cite{Bila}, then this theorem and it's results can be summarized as the following three-steps procedure:

 \textbf{step 1:} Find the determining equations of the extended classical symmetries related to the Eq.\eqref{Hesi}. For the meaning of the extended classical symmetries, a vector
 \begin{equation}
 V=\xi(x,y,z,u,f)\partial_x +\zeta\partial_y +\eta\partial_z +\phi\partial_u +\psi\partial_f.\label{v}
 \end{equation}
 is said the extended classical symmetry operator assosiated with $\mathrm{HESI}$ Equation and the determining equations of the extended classical symmetries related to the $\mathrm{HESI}$ Equation is the following equation:
 \begin{equation}
 {\bf pr}^{(2)}V[\mathrm{HESI}]=0,\label{sys}
 \end{equation}
 where $\xi$, $\zeta$, $\eta$, $\phi$ and $\psi$ are functions of $x$, $y$, $z$, $u$ and $f$, and ${\bf pr}^{(2)}V$ is
 \begin{align}
 V+\sum_{J}{\phi}^J(x,y,z,u^{(2)},f^{(2)})\,{\partial}_{u_J}+\sum_{J}{\psi}^J(x,y,z,u^{(2)},f^{(2)})\,{\partial}_{f_J}
 \end{align}
where $u^{(2)}=(u, u_x,\cdots,u_{zz})$ and $f^{(2)}=(f, f_x,\cdots, f_{zz})$:
 and the coefficients obtain from:
 \begin{align}
 {\phi}^J(x,y,z,u^{(2)},f^{(2)})=\textbf{D}_J\bigg(\phi-\sum_{i=1}^{3}{\xi}^i\frac{\partial u}{\partial x_i}\bigg)+\sum_{i=1}^{3}{\xi}^i\frac{\partial {u_J}}{\partial x_i}\nonumber\\
 {\psi}^J(x,y,z,u^{(2)},f^{(2)})=\textbf{D}_J\bigg(\psi-\sum_{i=1}^{3}{\xi}^i\frac{\partial f}{\partial x_i}\bigg)+\sum_{i=1}^{3}{\xi}^i\frac{\partial {f_J}}{\partial x_i}
 \end{align}
 where $({\xi}^1,{\xi}^2,{\xi}^3)=(\xi,\zeta,\eta)$ and $(x_1,x_2,x_3)=(x,y,z)$, further $J=(j_1,\cdots,j_k)$ is a $k$-th order multi-index, and $j_i$s adopt $x$, $y$ or $z$, for $1\leqslant i \leqslant k$, then $\textbf{D}_J$ denotes the total derivatives for the multi-index $J$, and the $J$-th total derivative is as
$\textbf{D}_J =\textbf{D}_{j_1}\cdots\textbf{D}_{j_k}$, that total derivative operator with respect to $i$ is as following
 \begin{equation*}
\textbf{D}_i=\partial_{x_i}+\sum_{J}\frac{\partial {u_J}}{\partial x_i}\,{\partial}_{u_J}+\sum_{J}\frac{\partial {u_f}}{\partial x_i}\,{\partial}_{f_J},
 \end{equation*}
\textbf{Note:} ${\bf pr}^{(2)}V$ is determined by taking into account that $u$ and $f$ are both dependent variables, exactly as one would proceed in finding the classical Lie symmetries for a system without arbitrary functions.

So with solving equation \eqref{sys} we gain:
\begin{equation}\label{sys2}
\begin{aligned}
& \xi=c_6x+c_9y+c_7z+c_8, & & \zeta=c_{10}z+c_6y-c_9x+c_{11}, \\
& \psi=2f(-2c_6+c_3), & & \eta=-c_{10}y+c_6z-c_7x+c_{12}, \\
& \phi=c_1x+c_3u+c_2y+c_5z+c_4,
\end{aligned}
\end{equation}
which $c_i,i=1,\cdots,12$ are arbitrary costants.

\textbf{step 2:} Augment the system of step 1 with the following conditions:
\begin{equation}
\frac{\partial {\xi}}{\partial u}=0,\qquad \frac{\partial {\zeta}}{\partial u}=0,\qquad\frac{\partial {\eta}}{\partial u}=0,\qquad \frac{\partial {\psi}}{\partial u}=0.
 \end{equation}
As we seen in relations \eqref{sys2} above conditions are satiesfied.

\textbf{step 3:} Augment the system of steps 1 and 2 with the following conditions:
\begin{equation}
\frac{\partial {\xi}}{\partial f}=0,\qquad \frac{\partial {\zeta}}{\partial f}=0,\qquad\frac{\partial {\eta}}{\partial f}=0,\qquad \frac{\partial {\psi}}{\partial f}=0.
 \end{equation}
As we seen in relations \eqref{sys2} above conditions are satisfied too.

Ultimately, The class of equations \eqref{Hesi} has a finite continuous group of equivalence transformations generated by the following infinitesimal operators:
\begin{equation} \label{vectors2}
\begin{aligned}
Y_1 &=\partial_x , & Y_2 &=\partial_y , & Y_3 &=\partial_z ,\\
Y_4 &=\partial_u , & Y_5 &=x\partial_u , & Y_6 &=y\partial_u ,\\
Y_7 &=z\partial_u , & Y_8 &=z\partial_x -x\partial_z , & Y_9 &=y\partial_x -x\partial_y ,\\
Y_{10} &=z\partial_y -y\partial_z , & Y_{11} &=u\partial_u +2f\partial_f, & Y_{12} &=x\partial_x +y\partial_y +z\partial_z -4f\partial_f.
\end{aligned}
\end{equation}
Moreover, in the group of equivalence transformations are included also discrete transformations, i.e., reflections
$(x,y,z,u,f)\mapsto-(x,y,z,u,f)$.
\section{Sketch of the method of preliminary group classification}
In many applications of group analysis, most of extensions of the principal Lie algebra admitted by an equation are obtained from the equivalence algebra $\mathfrak{g}_{\mathscr{E}}$. We call these extension ${\mathscr{E}}$-\textit{extension of the principal Lie algebra}. The classification of all nonequivalent equations admitting ${\mathscr{E}}$-extension of the principal Lie algebra is called \textit{a preliminary qroup classification}. What we obtain also is not necessarily the largest equivalence group but, it can be any subgroup of the qroup of all equivalence transformations.

The application of this method is effective and simple when it is based on a finite-dimensional equivalence algebra $\mathfrak{g}_{\mathscr{E}}$. So we take finite dimensional algebra $\mathfrak{g}_{12}$ spanned on the basis \eqref{vectors2} and use it for preliminary group classification.

The function $f$ of Eq.\eqref{Hesi} depends on the variables $x,y,z$, so we don't construct any prolongations of operators \eqref{Y}. But we take projections on the space $(x,y,z,f)$.

The nonzero projections of \eqref{vectors2} are:
\begin{align}
Z_1 &={\bf pr}(Y_1)=\partial_x , & Z_2 &={\bf pr}(Y_2)=\partial_y ,\nonumber\\
Z_3 &={\bf pr}(Y_3)=\partial_z , & Z_4 &={\bf pr}(Y_8)=z\partial_x -x\partial_z , \nonumber\\
Z_5 &={\bf pr}(Y_9)=y\partial_x -x\partial_y , & Z_6 &={\bf pr}(Y_{10})z\partial_y -y\partial_z ,\nonumber\\
Z_7 &={\bf pr}(Y_{11})=2f\partial_f, & Z_8 &={\bf pr}(Y_{12})=x\partial_x +y\partial_y +z\partial_z -4f\partial_f.\label{vectors3}
\end{align}
It's clear that there aren't the minimal infinitesimal generators \eqref{principle}, among above vectors.

The Lie algebra generated with the vectors in \eqref{vectors3} is denoted by $\mathfrak{g}_{8}$.

The essence of the preliminary method is based on the following two proposition:

\begin{proposition}:\label{proposition1}
Let $\mathfrak{g}_{m}$ be a $m-$dimensional subalgebra of $\mathfrak{g}_{8}$. Suppose $Z^{(i)}$, $i=1,\cdots,m$ be a basis of $\mathfrak{g}_{m}$ and $Y^{(i)}$ is the elements of the algebra $\mathfrak{g}_{12}$, such that $Z^{(i)}={\bf pr}(Y^{(i)})$, that means, if
\begin{equation}\label{1}
 Z^{(i)}=\sum_{\alpha=1}^{8}e^\alpha_iZ_{\alpha},
 \end{equation}
 then with respect to \eqref{vectors2} and \eqref{vectors3}:
 \begin{align}\label{2}
 Y^{(i)}=e^{1}_iY_{1}+e^{2}_iY_{2}+e^{3}_iY_{3}+e^{4}_iY_{8}+e^{5}_iY_{9}+e^{6}_iY_{10}+
 e^{7}_iY_{11}+e^{8}_iY_{12}.
\end{align}
If function $f=f(x,y,z)$ be invariant with respect to the algebra $\mathfrak{g}_{m}$, then the $\mathrm{HESI}$ equation admits the operators
 \begin{equation}\label{3}
 X^{(i)}=\mbox{projection of}\;\;Y^{(i)} \mbox{on}\;\; (x,y,z,u).
\end{equation}
\end{proposition}
\begin{proposition}:\label{proposition2}
Let equations
\begin{eqnarray}
 {\cal S}_2[u]=f(x,y,z), \label{eq1} \\
{\cal S}_2[u]=f^{\prime}(x,y,z),\label{eq2}
\end{eqnarray}
 be constructed according to proposition \eqref{proposition1} with subalgebras $\mathfrak{g}_{m}$ and $\mathfrak{g}_{m^{\prime}}$, respectively. If $\mathfrak{g}_{m}$ and $\mathfrak{g}_{m^{\prime}}$ are similar subalgebras in $\mathfrak{g}_{12}$ then equations \eqref{eq1} and \eqref{eq2} are equivalent with respect to the equivalence group $G_{12}$ generated by $\mathfrak{g}_{12}$.
\end{proposition}

According to above propositions, continuation of the preliminary group classification of Eq.\eqref{Hesi} with respect to the finite-dimensional algebra $\mathfrak{g}_{12}$, is reduced to the algebraic problem of constructing of nonsimilar subalgebras of $\mathfrak{g}_{8}$, or optimal systems of subalgebras.

\textbf{note:} In this paper we just solve the problem of preliminary group classification with respect to one-dimensional subalgebras.
\section{Adjoint group for algebra $\mathfrak{g}_{8}$ }
We determine a list or optimal system, of conjuacy inequivalent subalgebras with the property that any other subalgebra is equivalent to a unique member of the list under some element of the adjoint representation, i.e. $\bar{\mathfrak{h}}={\bf Ad }(g)\mathfrak{h}$ for some g of a considered Lie group, see~\cite{O1,O2,Ov1}.

The adjoint action is given by the Lie series
 \begin{equation}
 {\bf Ad }(\exp(\varepsilon Y_{i}))Y_{j}=Y_{j}-\varepsilon [Y_{i},Y_{j}]+\frac{\varepsilon^2}{2}[Y_{i},[Y_{i},Y_{j}]]-\cdots,
 \end{equation}
The commutator and adjoint representations of $\mathfrak{g}_{8}$ are listed in tables \ref{tab1} and \ref{tab2}.

\begin{table}[ht]
\caption{Commutators table for $\mathfrak{g}_{8}$: $[Z_{i},Z_{j}]$} \label{tab1}
\begin{center}
\begin{tabular}{c|cccccccc}
 & $Z_{1}$ & $Z_{2}$ & $Z_{3}$ & $Z_{4}$ & $Z_{5}$ & $Z_{6}$ & $Z_{7}$ & $Z_{8}$\\
\hline
$Z_{1}$ & 0 & 0 & 0 & $-Z_{3}$ & $-Z_{2}$ & 0 & 0 & $Z_{1}$\\
$Z_{2}$ & 0 & 0 & 0 & 0 & $Z_{1}$ & $-Z_{3}$ & 0 & $Z_{2}$\\
$Z_{3}$ & 0 & 0 & 0 & $Z_{1}$ & 0 & $Z_{2}$ & 0& $Z_{3}$\\
$Z_{4}$ & $Z_{3}$ & 0 & $-Z_{1}$ & 0 & $-Z_{6}$ & $Z_{5}$ & 0 & 0\\
 $Z_{5}$ & $Z_{2}$ & $-Z_{1}$ & 0 & $Z_{6}$ & 0 & $-Z_{4}$ & 0 & 0\\
$Z_{6}$ & 0 & $Z_{3}$ & $-Z_{2}$ & $-Z_{5}$ & $Z_{4}$ & 0 & 0 & 0\\
$Z_{7}$ & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
 $Z_{8}$ & $-Z_{1}$ & $-Z_{2}$ & $-Z_{3}$ & 0 & 0 & 0 & 0 & 0\\
\end{tabular}
\end{center}
\end{table}
\begin{table}[ht]
\begin{center}
\caption{Adjoint table for $\mathfrak{g}_{8}$: ${\bf Ad }(\exp(\varepsilon_i Y_{i}))Y_{j}$} \label{tab2}
\scalebox{0.8}{
\begin{tabular}{c|cccccccc}
  & $Z_{1}$ & $Z_{2}$ & $Z_{3}$ & $Z_{4}$ & $Z_{5}$ & $Z_{6}$ & $Z_{7}$ & $Z_{8}$\\
\hline
$Z_{1}$ & $Z_{1}$ & $Z_{2}$ & $Z_{3}$ & $\varepsilon_1Z_{3}+Z_{4}$ & $\varepsilon_1Z_{2}+Z_{5}$ & $Z_{6}$ & $Z_{7}$ & $-\varepsilon_1Z_{1}+Z_{8}$\\
\hline
$Z_{2}$ & $Z_{1}$ & $Z_{2}$ & $Z_{3}$ & $Z_{4}$ & $-\varepsilon_2Z_{1}+Z_{5}$ & $\varepsilon_2Z_{3} +Z_{6}$ & $Z_{7}$ & $-\varepsilon_2Z_{2}+Z_{8}$\\
\hline
$Z_{3}$ & $Z_{1}$ & $Z_{2}$ & $Z_{3}$ & $-\varepsilon_3Z_{1}+Z_{4}$ & $Z_{5}$ & $-\varepsilon_3Z_{2} +Z_{6}$ & $Z_{7}$ & $-\varepsilon_3Z_{3}+Z_{8}$\\
\hline
$Z_{4}$ & $\displaystyle  \cos(\varepsilon_4)Z_{1} \atop \displaystyle -\sin(\varepsilon_4)Z_{3}$ & $Z_{2}$ & $\displaystyle  \sin(\varepsilon_4)Z_{1} \atop \displaystyle   +\cos(\varepsilon_4)Z_{3}$ & $Z_{4}$ & $\displaystyle \cos(\varepsilon_4)Z_{5} \atop \displaystyle +\sin(\varepsilon_4)Z_{6}$ & $\displaystyle -\sin(\varepsilon_4)Z_{5} \atop \displaystyle  +\cos(\varepsilon_4)Z_{6}$ & $Z_{7}$ & $Z_{8}$ \\
\hline
 $Z_{5}$ & $\displaystyle  \cos(\varepsilon_5)Z_{1} \atop \displaystyle  -\sin(\varepsilon_5)Z_{2}$ & $\displaystyle  \sin(\varepsilon_5)Z_{1} \atop \displaystyle  +\cos(\varepsilon_5)Z_{2}$ & $Z_{3}$ & $\displaystyle  \cos(\varepsilon_5)Z_{4} \atop \displaystyle -\sin(\varepsilon_5)Z_{6}$ & $Z_{5}$ & $\displaystyle  \sin(\varepsilon_5)Z_{4} \atop \displaystyle  +\cos(\varepsilon_5)Z_{6}$ & $Z_{7}$ & $Z_{8}$ \\
 \hline
$Z_{6}$ & $Z_{1}$ & $\displaystyle  \cos(\varepsilon_6)Z_{2} \atop \displaystyle  -\sin(\varepsilon_6)Z_{3}$ & $\displaystyle  \sin(\varepsilon_6)Z_{2} \atop \displaystyle   +\cos(\varepsilon_6)Z_{3}$ & $\displaystyle  \cos(\varepsilon_6)Z_{4} \atop \displaystyle  +\sin(\varepsilon_6)Z_{5}$ & $\displaystyle  -\sin(\varepsilon_6)Z_{4} \atop \displaystyle   +\cos(\varepsilon_6)Z_{5}$ & $Z_{6}$ & $Z_{7}$ & $Z_{8}$ \\
\hline
$Z_{7}$ & $Z_{1}$ & $Z_{2}$ & $Z_{3}$ & $Z_{4}$ & $Z_{5}$ & $Z_{6}$ & $Z_{7}$ & $Z_{8}$\\
\hline
$ Z_{8}$ & $e^{\varepsilon_8}Z_{1}$ & $e^{\varepsilon_8}Z_{2}$ & $e^{\varepsilon_8}Z_{3}$ & $Z_{4}$ & $Z_{5}$ & $Z_{6}$ & $Z_{7}$ & $Z_{8}$
\end{tabular}}
\end{center}
\end{table}
\section{Construction of the optimal system of one-dimensional subalgebras of $\mathfrak{g}_{8}$}
\begin{theorem}:\label{opti}
An optimal system of one-dimensional Lie algebras of $\mathfrak{g}_{8}$ in $\mathrm{HESI}$ equation are as follows:
\begin{equation}\label{op}
\begin{aligned}
& A^1=Z_{7}, & A^2&=\pm Z_{1}+Z_{7},\\
& A^3=\gamma_1Z_{6}+Z_{7}, & A^4&=\pm Z_{1}+\gamma_2Z_{6}+Z_{7},\\
& A^5=\alpha_1Z_{4}+Z_{7}, & A^6&=\pm Z_{2}+\alpha_2Z_{4}+Z_{7},\\
&A^7=\alpha_3Z_{4}+\gamma_3Z_{6}+Z_{7}, & A^8&=\pm Z_{1}+\alpha_4Z_{4}+\gamma_4Z_{6}+Z_{7},\\
& A^9=\alpha_5Z_{4}+\beta_1Z_{5}+Z_{7}, & A^{10}&=\pm Z_{3}+\alpha_6Z_{4}+\beta_2Z_{5}+Z_{7},\\
& A^{11}=\alpha_7Z_{4}+\beta_3Z_{5}+\gamma_5Z_{7}+Z_{8}, & A^{12}&=\pm Z_{2}+\alpha_8Z_{4}+\beta_4Z_{5}+\gamma_6Z_{7}+Z_{8},
\end{aligned}
\end{equation}
Where $\alpha_i$, $i=1,...8$ and $\beta_i$, $i=1,\cdots,4$ and $\gamma_i$, $i=1,\cdots,6$ are arbitrary constants.
\end{theorem}
\medskip \noindent {\it Proof:}
We will start with $Z=\sum_{i=1}^{8}a_iZ_{i}$, suppose Z is a nonzero vector field of $\mathfrak{g}_{8}$, we want simplify as many of the coefficients $a_i$, $i=1,\cdots,8$ as possible through proper Adjoint applications on Z. We proceed this simplifications through following cases:

\textbf{note:} The coefficients $a_7$ and $a_8$ don't change at all.

\textbf{Case 1:} At first, we assume that $a_8=0$, so with scalling on Z, we can suppose that $a_7=1$, then we have $Z=\sum_{i=1}^{6}a_iZ_{i}+Z_7$. therefore for different values of $a_5=0$, when it is either zero or nonzero, we have cases 1.1 and 1.2.

\textbf{Case 1.1:} If $a_8=a_5=0$, so $Z=a_1Z_{1}+a_2Z_{2}+a_3Z_{3}+a_4Z_{4}+a_6Z_{6}+Z_{7}$. Then for different values of $a_4$, when it is either zero or nonzero, we have cases 1.1.a and 1.1.b.

\textbf{Case 1.1.a:} If $a_8=a_5=a_4=0$, so $Z=a_1Z_{1}+a_2Z_{2}+a_3Z_{3}+a_6Z_{6}+Z_{7}$. Then for different values of $a_6$, when it is either zero or nonzero, we have cases 1.1.a1 and 1.1.a2.

\textbf{Case 1.1.a1:} If $a_8=a_5=a_4=a_6=0$, so $Z=a_1Z_{1}+a_2Z_{2}+a_3Z_{3}+Z_{7}$. Then for different values of $a_3$, when it is either zero or nonzero, the coefficient can be vanished; when $a_3\neq 0$, with effecting ${\bf Ad }(\exp(\cot^{-1}(a_1/{a_3})Z_4))$ on Z. Then we have $Z=a_1Z_{1}+a_2Z_{2}+Z_{7}$.

 Now if $a_2=0$ or $a_2\neq0$; by effecting ${\bf Ad }(\exp(\cot^{-1}(a_1/{a_2})Z_5))$ on Z, we can make the coefficient of $Z_2$ vanished.Then we have $Z=a_1Z_{1}+Z_{7}$.

 So if $a_1=0$, then $Z=Z_7$, so we have $A^1$.

 And if $a_1\neq 0$, with ${\bf Ad }(\exp(\ln(\pm1/{a_1})Z_8))$ change the coefficient of $Z_1$ equal $\pm 1$, so $Z=\pm Z_1+Z_7$, therefore we have $A^2$.

 \textbf{Case 1.1.a2:} If $a_8=a_5=a_4=0$, but $a_6\neq0$, so $Z=a_1Z_{1}+a_2Z_{2}+a_3Z_{3}+a_6Z_{6}+Z_{7}$. Then for different values of $a_3$, when it is either zero or nonzero, the coefficient can be vanished; when $a_3\neq 0$, with applying ${\bf Ad }(\exp(-a_3/{a_6})Z_2)$ on Z. So we have $Z=a_1Z_{1}+a_2Z_{2}+a_6Z_{6}+Z_{7}$.

 Similarly, the coefficient $a_2$ is either zero or we make it vanished with effecting ${\bf Ad }(\exp(a_2/{a_6})Z_3)$ on Z. Then we have $Z=a_1Z_{1}+a_6Z_{6}+Z_{7}$.

 So $a_1=0$ or $a_2\neq0$, if $a_1=0$, so $Z=a_6Z_{6}+Z_7$, and we have $A^3$.

 And if $a_1\neq 0$, with ${\bf Ad }(\exp(\ln(\pm1/{a_1})Z_8))$ change the coefficient of $Z_1$ equal $\pm 1$, so $Z=\pm Z_1++a_6Z_{6}+Z_7$, therefore we have $A^4$.

 \textbf{Case 1.1.b:} If $a_8=a_5=0$, but $a_4\neq0$, so $Z=a_1Z_{1}+a_2Z_{2}+a_3Z_{3}+a_4Z_{4}+a_6Z_{6}+Z_{7}$. Then for different values of $a_3$, when it is either zero or nonzero, the coefficient can be vanished; when $a_3\neq 0$, with effecting ${\bf Ad }(\exp(-a_1/{a_4})Z_1))$ on Z. Then we have $Z=a_1Z_{1}+a_2Z_{2}+a_4Z_{4}+a_6Z_{6}+Z_{7}$.

Therefore, for different values of $a_6$, when it is either zero or nonzero, we have cases 1.1.b1 and 1.1.b2.

\textbf{Case 1.1.b1:} If $a_8=a_5=a_3=a_6=0$, but $a_4\neq0$, so $Z=a_1Z_{1}+a_2Z_{2}+a_4Z_{4}+Z_{7}$. Then either $a_1$ is zero or nonzero, but
the coefficient can be vanished; when $a_1\neq 0$, with effecting ${\bf Ad }(\exp(a_1/{a_4})Z_3))$ on Z. Then we have $Z=a_2Z_{2}+a_4Z_{4}+Z_{7}$.

 Ultimately, $a_2=0$ or $a_2\neq0$; If $a_2=0$, so we have $Z=a_4Z_{4}+Z_{7}$.Then we have $A^5$.

 And if $a_2\neq 0$, with ${\bf Ad }(\exp(\ln(\pm1/{a_2})Z_8))$ change the coefficient of $Z_2$ equal $\pm 1$, so $Z=\pm Z_2+a_4Z_{4}+Z_7$, therefore we have $A^6$.

 \textbf{Case 1.1.b2:} If $a_8=a_5=a_3=0$, but $a_6\neq0$ and $a_4\neq0$, so $Z=a_1Z_{1}+a_2Z_{2}+a_4Z_{4}+a_6Z_{6}+Z_{7}$. Then for different values of $a_2$, when it is either zero or nonzero, the coefficient can be vanished; when $a_2\neq 0$, with applying ${\bf Ad }(\exp(a_2/{a_6})Z_3)$ on Z. So we have $Z=a_1Z_{1}+a_4Z_{4}+a_6Z_{6}+Z_{7}$.

 Similarly, the coefficient $a_1$ is either zero or nonzero, if $a_1=0$, so $Z=a_4Z_{4}+a_6Z_{6}+Z_7$, and we have $A^7$.

 And if $a_1\neq 0$, with ${\bf Ad }(\exp(\ln(\pm1/{a_1})Z_8))$ change the coefficient of $Z_1$ equal $\pm 1$, so $Z=\pm Z_1+a_4Z_{4}++a_6Z_{6}+Z_7$, therefore we have $A^8$.

 \textbf{Case 1.2:} If $a_8=0$ but $a_5\neq0$, so $Z=a_1Z_{1}+a_2Z_{2}+a_3Z_{3}+a_4Z_{4}+a_5Z_{5}+a_6Z_{6}+Z_{7}$. So we have different values of $a_2$, when it is either zero or nonzero, the coefficient can be vanished; when $a_2\neq 0$, with applying ${\bf Ad }(\exp(-a_2/a_5)Z_1)$ on Z. Then again we have two cases $a_1=0$ or $a_1\neq0$; If $a_1\neq0$, with applying ${\bf Ad }(\exp(a_1/a_5)Z_2)$ on Z, make it vanished. so we have $Z=a_3Z_{3}+a_4Z_{4}+a_5Z_{5}+a_6Z_{6}+Z_{7}$.

 Again we have two cases $a_6=0$ or $a_6\neq0$; If $a_6\neq0$, with applying ${\bf Ad }(\exp(\cot^{-1}(a_4/{a_6}))Z_5)$ on Z, make it vanished. so we have $Z=a_3Z_{3}+a_4Z_{4}+a_5Z_{5}+Z_{7}$.

 Ultimately, If $a_3=0$, then $Z=a_4Z_{4}+a_5Z_{5}+Z_{7}$.Then we have $A^9$.

 And if $a_3\neq0$, then with ${\bf Ad }(\exp(\ln(\pm1/{a_3})Z_8))$ change the coefficient of $Z_3$ equal $\pm 1$, so $Z=\pm Z_3+a_4Z_{4}+a_5Z_{5}+Z_{7}$.Then we have $A^{10}$.

 \textbf{Case 2:} If $a_8\neq0$, so with scalling on Z, we can suppose that $a_8=1$, then we have $Z=\sum_{i=1}^{7}a_iZ_{i}+Z_8$. So we have different values of $a_1$, when it is either zero or nonzero, the coefficient can be vanished; when $a_1\neq 0$, with applying ${\bf Ad }(\exp(a_1Z_1))$ on Z. So we reduce Z to $Z=a_2Z_{2}+a_3Z_{3}+a_4Z_{4}+a_5Z_{5}+a_6Z_{6}+a_7Z_{7}+Z_{8}$.

 Now if $a_3=0$ or $a_3\neq0$; by effecting ${\bf Ad }(\exp(\cot^{-1}(a_2/{a_3})Z_6)$ on Z, make the coefficient of $Z_3$ vanished. so we have $Z=a_2Z_{2}+a_4Z_{4}+a_5Z_{5}+a_6Z_{6}+a_7Z_{7}+Z_{8}$.

 Again we have two cases $a_6=0$ or $a_6\neq0$; If $a_6\neq0$, with applying ${\bf Ad }(\exp(-\cot^{-1}(a_5/{a_6}))Z_4)$ on Z, make it vanished. So we have $Z=a_2Z_{2}+a_4Z_{4}+a_5Z_{5}+a_7Z_{7}+Z_{8}$.

 Ultimately, If $a_2=0$, then $Z=a_4Z_{4}+a_5Z_{5}+a_7Z_{7}+Z_{8}$. Then we have $A^{11}$.

 And if $a_2\neq0$, then with ${\bf Ad }(\exp(\ln(\pm1/{a_2})Z_8))$ change the coefficient of $Z_2$ equal $\pm 1$, so $Z=\pm Z_2+a_4Z_{4}+a_5Z_{5}+a_7Z_{7}+Z_{8}$. Then we have $A^{12}$.

 There is not any more possible cases, and the proof is complete. \hfill\ $\Box$
\section{Equations admitting an extension by one of the principal Lie algebra}
Now based on propositions \eqref{proposition1} and \eqref{proposition2}, and with the optimal system \eqref{op}, we obtain all nonequivalent equations of the form equation \eqref{Hesi}, that admitting extension of principal Lie algebra $\mathfrak{g}$ by one operator $V_5$, that means every equation of the form equation \eqref{Hesi} admits symmetry group of dimension $4$ with infinitesimal generators \eqref{principle}, also together with a fifth operator $V_5$. for every case, when this extension occurs, we indicate the corresponding coefficients $f$ and additional operator $V_5$.

The algorithm will be clarified with these examples:
 %
\paragraph{First example:} Consider the operator $A^3=\gamma_1Z_{6}+Z_{7}$, which $\gamma_1\neq 0$, from \eqref{op}, so \begin{equation}
\begin{aligned}
 A^3=\gamma_1z\partial_y -\gamma_1y\partial_z +2f\partial_f.
 \end{aligned}
 \end{equation}

 Invariants are found from the following equation (see ~\cite{O2}):
 \begin{equation}
 \begin{aligned}
 \frac{dy}{\gamma_1z}=-\frac{d z}{\gamma_1y}=\frac{d f}{2f},
 \end{aligned}
 \end{equation}
 and are the following functions:
 \begin{equation}
 \begin{aligned}
 I_1=x,\qquad I_2=y^2+z^2,\qquad I_3=f\exp(-\frac{2}{\gamma_1}{\tan}^{-1}(\frac{y}{z})).
 \end{aligned}
 \end{equation}
It follows
 \begin{equation}
 \begin{aligned}
f={\exp(\frac{2}{\gamma_1}{\tan}^{-1}(\frac{y}{z}))}{H(x,y^2+z^2)},
 \end{aligned}
 \end{equation}
 where $H$ is arbitrary function.

 By applying the formulas \eqref{1}, \eqref{2} and \eqref{3} on the operator $A^3$ we obtain the additional operator $V_5=\gamma_1z\partial_y -\gamma_1y\partial_z +u\partial_u$. Thus, the equation
 \begin{equation}
 \mathrm{HESI}: \; {\cal S}_2[u]={\exp(\frac{2}{\gamma_1}{\tan}^{-1}(\frac{y}{z}))}{H(x,y^2+z^2)},
 \end{equation}
admits the five-dimensional algebra ${\mathfrak{g}}_5$ , that is generated with the following vectors
 \begin{align}
V_1&=\partial_u , & V_2&=x\partial_u , &V_3&=y\partial_u , \nonumber\\
V_4&=z\partial_u , & V_5&=\gamma_1z\partial_y -\gamma_1y\partial_z +u\partial_u.
\end{align}
\paragraph{Second example:} Consider the operator $A^1=Z_{7}$ from \eqref{op}, so $A^1=2f\partial_f$.

Invariants are the following functions:
 \begin{equation}\label{inv}
 \begin{aligned}
 I_1=x,\qquad I_2=y,\qquad I_3=z.
 \end{aligned}
 \end{equation}
 So, there are no invariant functions $f=f(x,y,z)$ because the necessary condition for existence of invariant solutions based on ref.~\cite{Ov1} (section 19.3) is not satisfied; that means invariants \eqref{inv} can't be solved with respect to $f$.

We continue calculations on some operators of \eqref{op}, and show results in table \ref{tab3}, that is the preliminary group classification of equation \eqref{Hesi}, which admit an extension ${\mathfrak{g}}_5$ of the principal Lie algebra ${\mathfrak{g}}$.

The results of classification:
\begin{longtable}{c|cc}
\caption{The equation ${\cal E}:\;{\cal S}_2[u]=f$ has $V_5$ as its additional operator v.r.t $A^s$}  \label{tab3}\\
           & $f$ & $V_5$ \\ \hline
$A^2$ & $ {\exp(\pm 2x)}{H(y,z)}$ &  $ \pm\partial_x +u\partial_u $\\
\hline
${A^3}_{({\gamma_1}\neq 0)}$ & $\displaystyle \exp(\frac{2}{\gamma_1}{\tan}^{-1}(\frac{y}{z})) H(x,y^2+z^2)$ &  $ \gamma_1z\partial_y -\gamma_1y\partial_z +u\partial_u $\\
\hline
${A^4}_{({\gamma_2}\neq 0)}$ &  $\displaystyle  \exp(\frac{2}{\gamma_2}{\tan}^{-1}(\frac{y}{z})) \atop \displaystyle  H(y^2+z^2,x\mp \frac{1}{\gamma_2}{\tan}^{-1}(\frac{y}{z}))$ &  $ \pm\partial_x +\gamma_2z\partial_y -\gamma_2y\partial_z +u\partial_u $\\
\hline
${A^5}_{({\alpha_1}\neq 0)}$ &  $ \exp(\frac{2}{\alpha_1}{\tan}^{-1}(\frac{x}{z}))H(y,x^2+z^2)$ &  $  \alpha_1z\partial_x -\alpha_1x\partial_z +u\partial_u $\\
\hline
${A^6}_{({\alpha_2}=0)}$ &  $ {\exp(\pm 2y)}{H(x,z)}$ & $ \pm\partial_y +u\partial_u $\\
\hline
${A^6}_{({\alpha_2}\neq 0)}$ & $\displaystyle  \exp(\frac{2}{\alpha_2}{\tan}^{-1}(\frac{x}{z})) \atop \displaystyle  H(x^2+z^2,y\mp \frac{1}{\alpha_2}{\tan}^{-1}(\frac{x}{z}))$ &  $ \pm\partial_y +\alpha_2z\partial_x -\alpha_2x\partial_z +u\partial_u $\\
\hline
${A^7}_{({\alpha_3}\neq 0,{\gamma_3}\neq 0)}$ &  $\displaystyle  \exp(\frac{2}{\sqrt{{\alpha_3}^2+{\gamma_3}^2}}{\tan}^{-1}{(\frac{\alpha_3x+\gamma_3y}{z\sqrt{\alpha_3^2+\gamma_3^2}})}) \atop \displaystyle  H(y-\frac{\gamma_3}{\alpha_3}x,(
(1-\frac{{\gamma_3}^2}{{\alpha_3}^2})x^2+\frac{2\gamma_3}{\alpha_3}xy+z^2))$ &  $\displaystyle  \alpha_3z\partial_x +\gamma_3z\partial_y \atop \displaystyle   -(\alpha_3x+\gamma_3y)\partial_z +u\partial_u $\\
\hline
${A^9}_{({\alpha_5}=0,{\beta_1}\neq 0)}$ &  $ \exp(\frac{2}{\beta_1}{\tan}^{-1}(\frac{x}{y}))H(z,x^2+y^2)$ &  $ \beta_1y\partial_x -\beta_1x\partial_y +u\partial_u$\\
\hline
${A^9}_{(\alpha_5\neq 0,\beta_1\neq 0)}$ &  $\displaystyle  \exp(\frac{-2}{\sqrt{\beta_1^2+\alpha_5^2}}\tan^{-1}(
\frac{\beta_1y+\alpha_5z}{x\sqrt{\alpha_5^2+\beta_1^2}})) \atop \displaystyle  H(z-\frac{\alpha_5}{\beta_1}y, (x^2+(1-\frac{\alpha_5^2}{\beta_1^2})y^2+\frac{2\alpha_5}{\beta_1}yz))$ &  $\displaystyle  (\alpha_5z+\beta_1y)\partial_x \atop \displaystyle   -\beta_1x\partial_y -\alpha_5x\partial_z +u\partial_u $\\
\hline
 ${A^{10}}_{(\alpha_6=\beta_2=0)}$ &  $ {\exp(\pm 2z)}{H(x,y)}$ &  $ \pm\partial_z +u\partial_u $\\
 \hline
${A^{10}}_{(\alpha_6\neq 0,\beta_2=0)}$ &  $\displaystyle  \exp(\frac{2}{\beta_2}{\tan}^{-1}(\frac{x}{y})) \atop \displaystyle  H(x^2+y^2,z\mp \frac{1}{\beta_2}{\tan}^{-1}(\frac{x}{y}))$ &  $ \pm\partial_z +\beta_2y\partial_x -\beta_2x\partial_y +u\partial_u $\\
\hline
${A^{11}}_{(\alpha_7=\beta_3=\gamma_5=0)}$ &  $ x^{-4}H(\frac{y}{x},\frac{z}{x})$ &  $ x\partial_x +y\partial_y +z\partial_z $\\
\hline
${A^{11}}_{(\alpha_7=\beta_3=0,\gamma_5\neq 0)}$ &  $ x^{2c-4}H(\frac{y}{x},\frac{z}{x})$ & $ x\partial_x +y\partial_y +z\partial_z $\\
\hline
${A^{12}}_{(\alpha_8=\beta_4=\gamma_6=0)}$ &  $ x^{-4}H(\frac{y\pm 1}{x},\frac{z}{x})$ &  $ x\partial_x +(y\pm 1)\partial_y +z\partial_z $\\
\hline
${A^{12}}_{(\alpha_8=\beta_4=0,\gamma_6\neq 0)}$ &  $ x^{2c-4}H(\frac{y\pm 1}{x},\frac{z}{x})$ &  $ x\partial_x +(y\pm 1)\partial_y +z\partial_z $
\end{longtable}
\section{Some Local Solutions }
Based on the following theorem of~\cite{Tian}, The $2-$hessian equation in $\mathbb{R}^3$,
\begin{equation}\label{2hesi}
{\cal S}_2[u]=f(x,y,z,u,{\cal D}u) \qquad  \mbox{on }  \Omega\subset\mathbb{R}^3,
\end{equation}
which $f\in\textit{C}^{\infty}(\Omega\times\mathbb{R}\times\mathbb{R}^3)$, ${\cal D}u=(\partial_1u,\cdots,\partial_nu)$, has $\textit{C}^{\infty}$ local solutions in $\mathbb{R}^3$, and the solution is in the following form:
\begin{equation}\label{solution}
u(x,y,z)=\frac{1}{2}(\tau_1x^2+\tau_2y^2+\tau_3z^2)+\varepsilon^5\omega(\varepsilon^{-2}(x,y,z)),
\end{equation}
where $\varepsilon$ and $\tau_i$  are arbitrary constants, and $\omega$ is a given smooth function.
\begin{theorem}
Assume that  $f\in\textit{C}^{\infty}(\Omega\times\mathbb{R}\times\mathbb{R}^3)$, then for any $Z_0=(x_0,u_0,p_0)\in\Omega\times\mathbb{R}\times\mathbb{R}^3$, we have that

 $(1)$ If $f(Z_0)=0$, then  \eqref{2hesi} admits a $1$-convex $\textit{C}^{\infty}$ local solution which is not convex.

   $(2)$ If $f\geq0$ near $Z_0$, then  \eqref{2hesi} admits a $2$-convex $\textit{C}^{\infty}$ local solution which is not convex.

   If $f(Z_0)>0$,  \eqref{2hesi} admits a convex $\textit{C}^{\infty}$ local solution.

   $(3)$ If $f(Z_0)<0$, then  \eqref{2hesi} admits a $1$-convex $\textit{C}^{\infty}$ local solution which is not $2$-convex. \\
Moreover, the equation \eqref{2hesi} is uniformly elliptic with respect to the above local solutions.
\end{theorem}

In this part we obtain the one-parameter groups generated by some operators, and since these groups are symmetry groups of $\mathrm{HESI}$ equation, then the above solution  transforms to another solution of $\mathrm{HESI}$ equation.

\begin{itemize}
\item\textbf{Case 1:} The operator $V_1=\partial_u$  that produces one-parameter group $G_1=(x,y,z, t+u)$, $t\in\mathbb{R}$, for equation ${\cal S}_2[u]=f(x,y,z)$, transforms solution\eqref{solution} to the following solution
\begin{align*}
u(x)=\frac{1}{2}(\tau_1x^2+\tau_2y^2+\tau_3z^2)-t+\varepsilon^5\omega(\varepsilon^{-2}(x,y,z)),
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 2:} The operator $V_2=x\partial_u$  that produces one-parameter group $G_2=(x,y,z, tx+u)$, $t\in\mathbb{R}$, for equation ${\cal S}_2[u]=f(x,y,z)$, transforms solution\eqref{solution} to the following solution
\begin{align*}
u(x)=\frac{1}{2}(\tau_1x^2+\tau_2y^2+\tau_3z^2)-tx+\varepsilon^5\omega(\varepsilon^{-2}(x,y,z)),
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
In such a manner are the following cases:
\item\textbf{Case 3:} The operator $V_3=y\partial_u$, $G_3=(x,y,z, ty+u)$, $t\in\mathbb{R}$, for equation ${\cal S}_2[u]=f(x,y,z)$; So
\begin{align*}
u(x)=\frac{1}{2}(\tau_1x^2+\tau_2y^2+\tau_3z^2)-ty+\varepsilon^5\omega(\varepsilon^{-2}(x,y,z)),
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 4:} The operator $V_4=z\partial_u$, $G_4=(x,y,z, tz+u)$, $t\in\mathbb{R}$, for equation ${\cal S}_2[u]=f(x,y,z)$; So
\begin{align*}
u(x)=\frac{1}{2}(\tau_1x^2+\tau_2y^2+\tau_3z^2)-tz+\varepsilon^5\omega(\varepsilon^{-2}(x,y,z)),
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 5:} The operator $V_5=\pm\partial_x +u\partial_u$, $G_5=(x\pm t,y,z, tu+u)$, $t\in\mathbb{R}$,\\for equation ${\cal S}_2[u]=exp(\pm2x)H(y,z)$; So
\begin{align*}
u(x)&=\frac{1}{2(t+1)}\big[\tau_1(x\pm{t})^2+\tau_2y^2+\tau_3z^2+2\varepsilon^5\omega(\varepsilon^{-2}(x\pm t,y,z))\big],
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 6:} The operator $V_5=\gamma_1z\partial_y -\gamma_1y\partial_z +u\partial_u$, \\$G_5=(x,z\sin(\gamma_1t)+y\cos(\gamma_1t),z\cos(\gamma_1t)-y\sin(\gamma_1t), tu+u)$, $t\in\mathbb{R}$,\\
 for equation ${\cal S}_2[u]=\displaystyle \exp(\frac{2}{\gamma_1}{\tan}^{-1}(\frac{y}{z})) H(x,y^2+z^2)$; So
\begin{align*}
u(x)&=\frac{1}{2(t+1)}\big[\tau_1x^2+\tau_2(z\sin(\gamma_1t)+y\cos(\gamma_1t))^2+\tau_3(z\cos(\gamma_1t)-y\sin(\gamma_1t))^2\\
&+2\varepsilon^5\omega(\varepsilon^{-2}(x,z\sin(\gamma_1t)+y\cos(\gamma_1t),z\cos(\gamma_1t)-y\sin(\gamma_1t)))\big],
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 7:} The operator $V_5=\pm\partial_x +\gamma_2z\partial_y -\gamma_2y\partial_z +u\partial_u$,\\
$G_5=(x\pm t,z\sin(\gamma_2t)+y\cos(\gamma_2t),z\cos(\gamma_2t)-y\sin(\gamma_2t), tu+u)$, $t\in\mathbb{R}$,\\ for equation ${\cal S}_2[u]=\displaystyle  \exp(\frac{2}{\gamma_2}{\tan}^{-1}(\frac{y}{z}))\displaystyle  H(y^2+z^2,x\mp \frac{1}{\gamma_2}{\tan}^{-1}(\frac{y}{z}))$; So
\begin{align*}
u(x)&=\frac{1}{2(t+1)}\big[\tau_1(x\pm t)^2+\tau_2(z\sin(\gamma_2t)+y\cos(\gamma_2t))^2+\tau_3(z\cos(\gamma_2t)-y\sin(\gamma_2t))^2\\
&+2\varepsilon^5\omega(\varepsilon^{-2}(x\pm t,z\sin(\gamma_2t)+y\cos(\gamma_2t),z\cos(\gamma_2t)-y\sin(\gamma_2t)))\big],
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 8:} The operator $V_5=\alpha_1z\partial_x -\alpha_1x\partial_z +u\partial_u $,\\ $G_5=(z\sin(\alpha_1t)+x\cos(\alpha_1t),y,z\cos(\alpha_1t)-x\sin(\alpha_1t), tu+u)$, $t\in\mathbb{R}$,\\
 for equation ${\cal S}_2[u]=\exp(\frac{2}{\alpha_1}{\tan}^{-1}(\frac{x}{z}))H(y,x^2+z^2)$; So
\begin{align*}
u(x)&=\frac{1}{2(t+1)}\big[\tau_1(z\sin(\alpha_1t)+x\cos(\alpha_1t))^2+\tau_2y^2+\tau_3(z\cos(\alpha_1t)-x\sin(\alpha_1t))^2\\
&+2\varepsilon^5\omega(\varepsilon^{-2}(z\sin(\alpha_1t)+x\cos(\alpha_1t),y,z\cos(\alpha_1t)-x\sin(\alpha_1t)))\big],
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 9:} The operator $V_5=\pm\partial_y +u\partial_u$, $G_5=(x,y\pm t,z, tu+u)$, $t\in\mathbb{R}$,\\ for equation ${\cal S}_2[u]=exp(\pm2y)H(x,z)$; So
\begin{align*}
u(x)=\frac{1}{2(t+1)}\big[\tau_1x^2+\tau_2(y\pm t)^2+\tau_3z^2+2\varepsilon^5\omega(\varepsilon^{-2}(x,y\pm t,z))\big],
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 10:} The operator $V_5= \pm\partial_y +\alpha_2z\partial_x -\alpha_2x\partial_z +u\partial_u$,\\ $G_5=(z\sin(\alpha_2t)+x\cos(\alpha_2t),y\pm t,z\cos(\alpha_2t)-x\sin(\alpha_2t), tu+u)$, $t\in\mathbb{R}$,\\
 for equation ${\cal S}_2[u]=\displaystyle  \exp(\frac{2}{\alpha_2}{\tan}^{-1}(\frac{x}{z}))\displaystyle  H(x^2+z^2,y\mp \frac{1}{\alpha_2}{\tan}^{-1}(\frac{x}{z}))$; So
\begin{align*}
u(x)&=\frac{1}{2(t+1)}\big[\tau_1(z\sin(\alpha_2t)+x\cos(\alpha_2t))^2+\tau_2(y\pm t)^2+\tau_3(z\cos(\alpha_2t)-x\sin(\alpha_2t))^2\\
&+2\varepsilon^5\omega(\varepsilon^{-2}(z\sin(\alpha_2t)+x\cos(\alpha_2t),y\pm t,z\cos(\alpha_2t)-x\sin(\alpha_2t)))\big],
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 11:} The operator $V_5=\beta_1y\partial_x -\beta_1x\partial_y +u\partial_u $,\\
$G_5=(y\sin(\beta_1t)+x\cos(\beta_1t),y\cos(\beta_1t)-x\sin(\beta_1t),z, tu+u)$, $t\in\mathbb{R}$,\\
 for equation ${\cal S}_2[u]= \exp(\frac{2}{\beta_1}{\tan}^{-1}(\frac{x}{y}))H(z,x^2+y^2)$; So
\begin{align*}
u(x)&=\frac{1}{2(t+1)}\Big[\tau_1(y\sin(\beta_1t)+x\cos(\beta_1t))^2+\tau_2(y\cos(\beta_1t)-x\sin(\beta_1t))^2+\tau_3z^2\\
&+2\varepsilon^5\omega(\varepsilon^{-2}(y\sin(\beta_1t)+x\cos(\beta_1t),y\cos(\beta_1t)-x\sin(\beta_1t),z))\Big],
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 12:} The operator $V_5=\pm\partial_z +u\partial_u$,  $G_5=(x,y,z\pm t, tu+u)$, $t\in\mathbb{R}$,\\ for equation ${\cal S}_2[u]= {\exp(\pm 2z)}{H(x,y)}$; So
\begin{align*}
u(x)=\frac{1}{2(t+1)}\Big[\tau_1x^2+\tau_2y^2+\tau_3(z\pm t)^2+2\varepsilon^5\omega(\varepsilon^{-2}(x,y,z\pm t))\Big],
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 13:} The operator $V_5=\pm\partial_z +\beta_2y\partial_x -\beta_2x\partial_y +u\partial_u$,\\  $G_5=(y\sin(\beta_2t)+x\cos(\beta_2t),y\cos(\beta_2t)-x\sin(\beta_2t),z\pm t, tu+u)$, $t\in\mathbb{R}$,\\
 for equation ${\cal S}_2[u]=\displaystyle  \exp(\frac{2}{\beta_2}{\tan}^{-1}(\frac{x}{y}))\displaystyle  H(x^2+y^2,z\mp \frac{1}{\beta_2}{\tan}^{-1}(\frac{x}{y}))$; So
\begin{align*}
u(x)&=\frac{1}{2(t+1)}\big[\tau_1(y\sin(\beta_2t)+x\cos(\beta_2t))^2+\tau_2(y\cos(\beta_2t)-x\sin(\beta_2t))^2+\tau_3(z\pm t)^2\\
&+2\varepsilon^5\omega(\varepsilon^{-2}(y\sin(\beta_2t)+x\cos(\beta_2t),y\cos(\beta_2t)-x\sin(\beta_2t),z\pm t))\big],
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 14:} The operator $V_5=x\partial_x+y\partial_y+z\partial_z$,\\
$G_5=(e^tx,e^ty,e^tz,u)$, $t\in\mathbb{R}$, for equation ${\cal S}_2[u]= x^{2c-4}H(\frac{y}{x},\frac{z}{x})$; So
\begin{align*}
u(x)=\frac{1}{2}(\tau_1e^{2t}x^2+\tau_2e^{2t}y^2+\tau_3e^{2t}z^2)+\varepsilon^5\omega(\varepsilon^{-2}(e^tx,e^ty,e^tz)),
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\item\textbf{Case 15:} The operator $V_5=x\partial_x+(y\pm 1)\partial_y+z\partial_z$, \\
$G_5=(e^tx,e^ty\pm e^t\mp1,e^tz,u)$, $t\in\mathbb{R}$, for equation ${\cal S}_2[u]= x^{2c-4}H(\frac{y\pm 1}{x},\frac{z}{x})$; So
\begin{align*}
u(x)&=\frac{1}{2}(\tau_1e^{2t}x^2+\tau_2(e^ty\pm e^t\mp1)^2+\tau_3e^{2t}z^2)+\varepsilon^5\omega(\varepsilon^{-2}(e^tx,e^ty\pm e^t\mp1,e^tz)),
\end{align*}
where $\varepsilon,\tau_i\in\mathbb{R}$.
\end{itemize}
\section{Conclusion}
In this paper,  we performed preliminary group classification on equation, by studying the class of $3$-dimensional nonlinear $2-$hessian equations ${\cal S}_2[u]=f(x,y,z)$, and investigating the algebraic structure of the symmetry groups for the equation. Then, we obtained an optimal system of one-dimensional Lie subalgebras of this equation, with the aid of propositions \eqref{proposition1} and \eqref{proposition2}. The result of these work is a wide class of equations which summarized in table \ref{tab3}. And at the end of this work, some exact solutions of $2-$hessian equation are presented.
 Of course, the results in table \ref{tab3} can be continued for remainder vectors, and it is possible to obtain the corresponding reduced equations for all the cases in the classification in up comming works.