\section{Introduction}
\label{intro}
The graphene, the first member of two-dimensional (2D) materials, with honeycomb structure does not exhibit an electronic band gap \cite{neto2009electronic}. This problem has motivated the researchers to explore and design novel two-dimensional semiconductors, such as the transition metal dichalcogenides \cite{radisavljevic2011single}, black phosphorene \cite{guan2014phase}, antimonene \cite{zhang2015atomically} and indium selenide\cite{mudd2013tuning}. Therefore, exploring new 2D materials with proper electronic properties for particular applications is highly demanded.

After monolayer TMDC, 2D transition metal nitride (TMN) have been proposed in the recent years \cite{li2017computational,wang2016strain,ozdemir2018intercalation}. Since TMNs do not contain layered structures, the fabrication of large-area TMNs monolayers has remained challenging. Recently, MoSi$_2$N$_4$ monolayer based on TMN and without 3D counterpart structure has been synthesized\cite{hong2020chemical}. MoSi$_2$N$_4$ is constructed from a MoN$_2$ monolayer sandwiched between two Si-N monolayers. By analyzing MoSi$_2$N$_4$ monolayer, it can be realized that with sandwiching a TMDC-type MZ$_2$ monolayer into InSe-type A$_2$Z$_2$, twelve different structures with MA$_2$Z$_4$ formula have been achieved. They demonstrate different phases such as semiconductor, topological insulators, and Ising superconductors \cite{wang2021}.
Monolayer MoSi$_2$N$_4$ has been successfully fabricated using chemical vapor deposition with large size of up to 15mm$\times$15mm\cite{hong2020chemical}. This compound demonstrates high strength and excellent stability at ambient. MoSi$_2$N$_4$ compound also shows an elastic constant three times larger than monolayer MoS$_2$. In addition, the large electron and hole mobility is estimated to be 270/1200 cm$^2$/Vs that is four to six times higher than MoS$_2$ ones. This high mobility along with outstanding stability in the ambient environment makes this material a promising candidate for future applications.

Other members of MoSi$_2$N$_4$ family with MA$_2$Z$_4$ formula (M = Transition metal, A = Si, Ge and Z = N, P, As) have been extensively investigated \cite{zhong2021strain,wang2021}. These monolayers exhibit a wide range from semiconducting to metallic properties. In addition, some compounds with magnetic transition metal elements also represent magnetic properties\cite{li2021design}. The MoSi$_2$N$_4$ and WSi$_2$N$_4$ monolayers also show high lattice thermal conductivity for thermoelectric applications \cite{mortazavi2021exceptional,yu2021high}.
Monolayer MoSi$_2$N$_4$ and its family also have suitable band gaps up to 1.73eV for potential optical applications in the visible range\cite{bafekry2021mosi2n4,yao2021novel}.

It has been anticipated that MoSi$_2$N$_4$ family compounds have a pair of valley pseudospins similar to monolayer MoS$_2$ \cite{yang2021valley}. These monolayers have been also predicted to exceed monolayer MoS$_2$ due to fascinating valleytronic properties, for instance, multiple-folded valleys in monolayer MoSi$_2$As$_4$.
The lack of inversion symmetry along with strong spin-orbit coupling from transition metal elements in the MA$_2$Z$_4$ family results in two inequivalent valleys (K and K') with sizable spin-splitting and spin-valley coupling \cite{ai2021theoretical}. For example, WSi$_2$N$_4$ demonstrates a large spin-splitting at both conduction and valence bands at K and K'-valleys together with spin-valley coupling.

It has been confirmed theoretically and experimentally that a semiconductor to metal transition occurs in multi-layer MoS$_2$ by applying vertical strain \cite{nayak2014pressure,chi2014pressure,ghobadi2019normal}. Totally, the decrease of interlayer distance may lead to charge redistribution and semiconductor to metal transition. The out-of-plane strain is presented as a powerful tool to tune the electrical properties of bilayer MA$_2$Z$_4$ family \cite{bafekry2021effect}. For instance, the bandgap of bilayer MoSi$_2$N$_4$ decreases with increasing compressive strain and the bandgap closes at the strain of 22$\%$. Such a semiconductor to metal transition also occurs in other MA$_2$Z$_4$ bilayers \cite{zhong2021strain}. The transition pressure is distributed from 2.18 GPa in CrSi$_2$N$_4$ to 32.04 GPa in TiSi$_2$N$_4$. In another work, the biaxial strain is used to tune the band gaps of bilayer MoSi$_2$N$_4$ and WSi$_2$N$_4$ \cite{wu2021semiconductor}. It has been reported that these compounds demonstrate a direct bandgap at compressive strain. Furthermore, it has been shown that strain can also modify the electronic and magnetic properties of VSi$_2$P$_4$ monolayer \cite{guo2020coexistence}. The strain increasing changes the phase of this monolayer from a ferromagnetic metal to a spin-gapless semiconductor, afterward to a ferromagnetic semiconductor, and then come back to spin-gapless semiconductor and finally to a ferromagnetic half-metal.

In this work, the effects of vertical strain on the electrical and spin properties of monolayer MoSi$_2$X$_4$ (X=N, P,  As, and Sb) are investigated. First, the electrical properties such as bandgap, effective mass, and band edge locations are obtained at equilibrium. In the following, vertical strain is applied to all materials, and the variation of band gaps, band edges, and charge densities are explored. Finally, spin-splitting at both equilibrium and strained samples is studied.

\section{Computational details}
In order to investigate the electrical and spin properties of MoSi$_2$X$_4$ (X= N, P, As and Sb) monolayers, density functional calculations are performed using the SIESTA package \cite{soler2002siesta}. The generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) \cite{perdew1981self} functional is employed for the exchange-correlation term. We have adopted fully relativistic pseudopotentials and have taken into account the effect of spin-orbit coupling (SOC). A Monkhorst-Pack k-point grid of $21\times21\times1$ is chosen for the unit-cell. The energy cutoff is set to be 200 Ry and a double-$\zeta$ plus polarization basis-set is used. The total energy is converged to better than $10^{-5}$ eV and the geometries are fully relaxed until the force on each atom is less than 0.02 eV/$\mathrm{\AA}$. A vacuum region of 30 $\mathrm{\AA}$ is added to avoid interactions in the normal direction. To visualize the atomic structures, XCrySDen package has been used \cite{kokalj2003computer}. The out-of-plane strain is defined as $\varepsilon=(t-t_0)/t_0$, where $t_0$ and $t$ are the equilibrium and deformed compound thickness, respectively. The effective masses are calculated by using the following equation \cite{touski2020interplay,ghobadi2020electrical},
\begin{equation}
m^*=\hbar^2/\left(\partial^2E/\partial k^2\right)
\end{equation}
Here, $\hbar$ is the reduced Planck constant, E and k are the energy and wave vector of conduction band minimum and valence band maximum.

\begin{figure}
	\centering
	\includegraphics[width=1.0\linewidth]{Fig1}
	\caption{The schematic of the studied materials with MoSi$_2$X$_4$ formula. The atomic structure of MoSi$_2$N$_4$ is depicted. N atoms are replaced by P, As and Sb in other studied structures.}
	\label{fig:fig1}
\end{figure}

\begin{figure*}
	\centering
	\includegraphics[width=1.0\linewidth]{phonon}
	\caption{Phonon dispersion relations for four compounds.}
	\label{fig:phonon}
\end{figure*}

\begin{figure*}
	\centering
	\includegraphics[width=1.0\linewidth]{con-val}
	\caption{Two-dimensional map of the conduction and valence bands for all compounds. The first Brillouin zone is indicated for MoSi$_2$N$_4$ as a sample. Top and down rows are for the first conduction and valence bands, respectively.}
	\label{fig:con-val}
\end{figure*}

\begin{table*}[t]
	\caption{ The lattice constant (a), the Mo-X (d$_{Mo-X}$) and Si-X (d$_{Si-X}$) bond lengths, the vertical distance between Si and X atoms ($\Delta_{Si-X}$), the thickness and the elastic constants (C$_{11}$ and C$_{12}$) of MoSi$_2$X$_4$ monolayers. \label{tab:tab1}}

	\begin{tabular}{p{1.5cm}p{1.7cm}p{1.7cm}p{1.7cm}p{1.7cm}p{2cm}p{1.7cm}p{1.7cm}}
		\hline
		\hline
		&a($\mathrm{\AA}$) & d$_{Mo-X}$($\mathrm{\AA}$)  & d$_{Si-X}$($\mathrm{\AA}$) & $\Delta_{Si-X}$($\mathrm{\AA}$) & Thickness($\mathrm{\AA}$) & C$_{11}$(N/m) & C$_{12}$(N/m) \\
		\hline
		MoSi$_2$N$_4$ &2.928  &2.111  & 1.763 &0.531  & 7.119 & 472.66 & 142.02    \\

		\hline
		MoSi$_2$P$_4$ &3.486  &2.477  & 2.253 & 1.046 &9.486  & 207.77 & 58.46 \\

		\hline
		MoSi$_2$As$_4$ &3.633  & 2.583 & 2.366 & 1.136 & 10.021 & 178.03 & 54.75\\

		\hline
		MoSi$_2$Sb$_4$ & 3.898 &2.777  &2.589  &1.278  &10.991  & 141.40 &47.72  \\

		\hline
		\hline
	\end{tabular}
\end{table*}

\section{Results and discussion}
\label{Result}
The schematic of the studied compounds is shown in Fig. \ref{fig:fig1}. As it is obvious, a MoN$_2$ monolayer is sandwiched between two SiN monolayers. The unsaturated Si atoms in SiN sub-layer bond with unsaturated N atoms in MoN$_2$. The structural and mechanical properties of four monolayers are presented in Table \ref{tab:tab1}. The lattice constant of MoSi$_2$N$_4$ is 2.928 $\AA$ that is close to experimental result\cite{hong2020chemical}. The lattice constants of MoSi$_2$P$_4$ and MoSi$_2$As$_4$ are also 3.486 and 3.633 $\AA$, respectively, that are compatible with previous studies \cite{guo2020coexistence,ai2021theoretical}. The lattice constant increases for the compounds with heavier X elements. This behavior can be observed in the Mo-X bond length ($d_{Mo-X}$) in MoX$_2$ sub-layer, Si-X bond length, and buckling height of SiN sub-layer.  Thickness of these compounds demonstrates a high value  of 7.119 for MoSi$_2$N$_4$ to 10.991 $\AA$ for MoSi$_2$Sb$_4$. This high thickness confirms that applying vertical strain is feasible for these compounds. The elastic constants C$_{11}$ and C$_{12}$ are also obtained and reported in the table which is compatible with previously reported amounts  \cite{yao2021novel}. These materials are stable and the elastic constants satisfy the Born criteria stability, $0<C_{11}$, $0<C_{22}$, $C_{12}<C_{11}$ and $C_{12}<C_{22}$. In addition, the stability of these compounds is studied by phonon dispersion, see Fig. \ref{fig:phonon}. The positive frequencies confirm the stability of these structures.

The two-dimensional contour plots of the first Brillouin zone for the valence and conduction bands are depicted in Fig. \ref{fig:con-val}. The conduction band minimum (CBM) is located at K-valley for all compounds. In the case of MoSi$_2$N$_4$, K-valley approximately demonstrates isotropic contour in different directions whereas, anisotropic behavior is obvious for the other compounds. The anisotropy of K-valley increases for heavier compounds and MoSi$_2$Sb$_4$ demonstrates the highest anisotropy. At the same time, the valence band maximum (VBM) is located at $\Gamma$-point in the case of MoSi$_2$N$_4$, whereas the energy of K-point is close to the $\Gamma$-point. On the other hand, the VBM of other structures is placed at K-point. MoSi$_2$N$_4$ demonstrates isotropic contours around $\Gamma$-point whereas the other compounds exhibit anisotropic behavior. However, K-point in the valence band shows a high anisotropy with triangular contours for all compounds.

The band structures of the four studied materials are drawn in Fig. \ref{fig:fig2}. MoSi$_2$N$_4$ demonstrates an indirect bandgap from $\Gamma$-valley in the valence to K-valley in the conduction band whereas, K-point at the valence band is close to $\Gamma$-point and contributes to the top of the valence band. Three other compounds have a direct bandgap at K-point. In these compounds, $\Gamma$-point at valence band also is close to K-point. On the other hand, K-point displays a considerable spin-splitting whereas, the splitting vanishes for $\Gamma$-point due to the high symmetry. The minimum of the second conduction band with a large spin-splitting is near to the minimum of the first band. In order to further analyze the electronic properties of the
structures, the projected density of states (PDOS) is also depicted in Fig. \ref{fig:fig2}. The states near the VBM and CBM are mainly contributed by the Mo atom. X atoms also contribute to the conduction and valence bands, but their contributions are lower than the Mo atom.

\begin{figure*}
	\centering
	\includegraphics[width=0.9\linewidth]{Fig2}
	\caption{The band structure along with corresponding PDOS of MoSi$_2$X$_4$ structures. The d-orbital of Mo atom and p-orbitals of Si and X atoms have the main contributions to the DOS.}
	\label{fig:fig2}
\end{figure*}

\begin{figure}
	\centering
	\includegraphics[width=1.0\linewidth]{Fig3}
	\caption{The band gap variation as a function of vertical strain.}
	\label{fig:fig3}
\end{figure}

\begin{figure*}
\centering
\includegraphics[width=1.0\linewidth]{Fig-energy}
\caption{The energies of the K-valleys of the first and second conduction bands along with K- and $\Gamma$-points of the valence band versus vertical strain. The up- and down-spin for each band are also indicated. The range of strain is kept in the semiconducting phase and is different for various compounds. }
\label{fig:energy}
\end{figure*}

\begin{figure*}
	\centering
	\includegraphics[width=1.0\linewidth]{charge}
	\caption{Charge distribution for MoSi$_2$P$_4$ as a sample. The charge is plotted for \textGamma- and K-point in the valence band and for the first and second conduction band at K-valley.}
	\label{fig:charge}
\end{figure*}

The electrical and spin properties of four materials are summarized in Table \ref{tab:tab2}. First, the band gaps with and without spin-orbit coupling consideration are compared. In the case of MoSi$_2$N$_4$, the band gaps without and with SOC are approximately the same. These values which are in good agreement with experimental work \cite{hong2020chemical}, are 1.73 and 1.721 eV, respectively. The bandgap of MoSi$_2$P$_4$ without(with) SOC is 0.772(0.648) eV that is near to 0.61 eV reported for with SOC\cite{ai2021theoretical}. The bandgap of MoSi$_2$As$_4$ decreases from 0.664 to 0.525eV with applying spin-orbit coupling, while the values of 0.6eV and 0.41-0.5 eV are reported for without and with spin-orbit coupling consideration, respectively\cite{hong2020chemical,guo2020coexistence,ai2021theoretical,li2020valley}. The SOC has more significant effects on the heavier compounds so that the bandgap reduces by 0.172 eV for MoSi$_2$Sb$_4$ with applying SOC. In addition, the conduction and valence band maximum is presented in the table. In the end, the spin-splitting at K-point in the valence band, the first and second conduction bands are listed in the table. The value of spin-splitting of MoSi$_2$N$_4$ at K-point ($\lambda_{K,V}$) of the valence band is 131 meV which is a little smaller than 140 meV from experimental work. However, 130 meV is reported for $\lambda_{K,V}$ in the theoretical papers \cite{li2020valley}. Furthermore, the calculated value of $\lambda_{K,V}$ for MoSi$_2$As$_4$ is 0.179 eV and in good agreement with previous reported values\cite{guo2020coexistence,li2020valley}. $\lambda_{K,V}$ increases with changing X-atom from N to Sb and reaches 220 meV for MoSi$_2$Sb$_4$. The spin-splitting at the first conduction band ($\lambda_{K,C1}$) is not considerable. On the contrary, the second band demonstrates a significant spin-splitting and for example, is larger than $\lambda_{K,V}$ in MoSi$_2$N$_4$, but the energy of the second band is much higher than the first band in this compound. The spin-splitting at K-valley in the second conduction band ($\lambda_{K,C2}$) decreases from 187 meV in MoSi$_2$N$_4$ to 83 meV in MoSi$_2$P$_4$ and increases by changing X atom from P to Sb. By looking at the conduction band of MoSi$_2$As$_4$ and MoSi$_2$Sb$_4$ at K-point, one can find that the first and third bands do not demonstrate remarkable spin-splitting, while the second and fourth display notable spin-splitting.

\begin{figure*}
	\centering
	\includegraphics[width=1.0\linewidth]{ssc-ssv}
	\caption{Two-dimensional map of the spin-splitting in the first and second conduction bands along with valence band. The first Brillouin zone is highlighted by the dashed line in the first compound. The large spin-splitting at valence and second conduction band is obvious.}
	\label{fig:ssc-ssv}
\end{figure*}

\begin{table*}[t]
	\caption{The band gap calculated with PBE functional (E$^{PBE}_g$) and PBE with spin-orbit coupling included (E$^{SOC}_g$), the energies of conduction (E$_C$-E$_F$) and valence band (E$_V$-E$_F$) edges, the compressive strain ($\varepsilon_{trans}$) and pressure (P$_{trans}$) required for semiconductor to metal transition, the spin-splitting at K point of the valence band ($\lambda_{K_V}$) and first and second bands of conduction ($\lambda_{K_{C1}}$ and $\lambda_{K_{C2}}$). \label{tab:tab2}}
	\begin{tabular}{p{1.5cm}p{1.2cm}p{1.2cm}p{1.2cm}p{1.2cm}p{1.2cm}p{1.2cm}p{1.2cm}p{1.2cm}p{1.2cm}p{1.2cm}p{1.2cm}}
		\hline
		\hline
		& E$^{PBE}_g$ & E$^{SOC}_g$ & VBM & CBM & E$_C$-E$_F$ & E$_V$-E$_F$ &$\varepsilon_{trans}$&$P_{trans}$&  $\lambda_{K_V}$  & $\lambda_{K_{C1}}$  & $\lambda_{K_{C2}}$    \\
		& (eV) & (eV) &  &  & (eV) & (eV)  &($\%$)&(GPa)& (eV) & (eV)  & (eV)  \\
		\hline

		MoSi$_2$N$_4$	&1.73 &1.721 &$\Gamma$ & K& 0.582 &-1.139 &-22&25.3& 0.131 &	0.003&	0.187 \\

		\hline

		MoSi$_2$P$_4$	&0.772 &0.648 &K &K &0.333 &-0.315 &-11&9.1&0.143&	0.004&	0.083 \\

		\hline

		MoSi$_2$As$_4$	&0.664 &0.525 &K &K &0.275 &-0.25 &-11&8.9& 0.179&	0.014&	0.113 \\

		\hline

		MoSi$_2$Sb$_4$	&0.481 &0.309 &K &K &0.179 &-0.131 &-10&7.3& 0.22&	0.018&	0.151 \\

		\hline
		\hline
	\end{tabular}
\end{table*}

\begin{table}[t]
	\caption{The effective masses at K- and $\Gamma$-points of the valence band and K-valley in the conduction band. All reported effective masses are in m$_0$ unit.
		\label{tab:tab3}}
	\begin{tabular}{p{1.5cm}p{1.0cm}p{1.0cm}p{1.0cm}p{1.0cm}p{1.0cm}p{1.0cm}}
		\hline
		\hline
		&$m^{v,*}_{K\rightarrow M}$ & $m^{v,*}_{K\rightarrow \Gamma}$ & $m^{v,*}_{\Gamma\rightarrow K}$ & $m^{v,*}_{\Gamma\rightarrow M}$ & $m^{c,*}_{K\rightarrow M}$ & $m^{c,*}_{K\rightarrow \Gamma}$   \\
		\hline
		MoSi$_2$N$_4$	&0.794&	0.595&	1.375&	1.37 &0.544&	0.485    \\
		\hline
		MoSi$_2$P$_4$	&0.562&	0.388&	1.483&	1.335 &0.697&	0.567    \\
		\hline
		MoSi$_2$As$_4$	&0.701&	0.438&	1.586&	1.407 &1.206&	0.625    \\
		\hline
		MoSi$_2$Sb$_4$	&0.925&	0.447&	2.288&	1.992 &1.843&	0.647    \\
		\hline
		\hline
	\end{tabular}
\end{table}

The effective masses of the conduction and valence bands are reported in Table \ref{tab:tab3}. The effective masses of both $\Gamma$- and K-points of the valence band are presented in the table. K-point demonstrates a lower effective mass with respect to $\Gamma$-one. The anisotropic contour around K-valley leads to two different values at K to M and K to $\Gamma$ paths. The K to M effective mass ($m^{V*}_{K\rightarrow M}$) displays a higher value than K to \textGamma~effective mass ($m^{V*}_{K\rightarrow \Gamma}$). On the other hand, \textGamma-point owns isotropic contour, and effective masses at different paths are approximately equal. In Ref. \onlinecite{hong2020chemical}, the effective mass of the electron is reported to be 0.486 m$_0$, which corresponds to our calculated value 0.485 m$_0$ for K to \textGamma~path. Furthermore, the reported value of the hole effective mass is 0.683 m$_0$ that is in the middle of our calculated effective masses at K-valley. However, \textGamma-point is the VBM and demonstrates a larger effective masses.
Except MoSi$_2$N$_4$, the effective masses for the valence band increase for heavier compounds and MoSi$_2$P$_4$ contains the lowest effective mass. On the other hand, only K-valley contributes to the conduction band. The effective masses at two different paths for this valley are reported in the table. The effective masses in these two paths are in the same range for MoSi$_2$N$_4$ and MoSi$_2$P$_4$, whereas their differences are enhanced for MoSi$_2$As$_4$ and MoSi$_2$Sb$_4$. As mentioned before, the contour around K-valley changes from circular for MoSi$_2$N$_4$ to triangular in MoSi$_2$Sb$_4$. In this regard, MoSi$_2$Sb$_4$ demonstrates high anisotropic effective masses at two different paths. The effective masses of the conduction band increase for the heavier compounds. MoSi$_2$N$_4$ and MoSi$_2$P$_4$ possess the lowest effective mass in the conduction and valence bands, respectively.

The vertical strain is known as a powerful tool to modify electrical and spin properties of two-dimensional compounds \cite{ghobadi2021structural,shamekhi2020band}. This strain can be applied by out-of-plane pressure and 2D materials can be used as a pressure sensor. The vertical strain is applied to the four structures and their electrical and spin properties are studied. First, the variation of the band gaps as a function of vertical strain is depicted in Fig. \ref{fig:fig3}. The band gaps exhibit a maximum value at a specific small tensile strain and then decrease for larger compressive and tensile strains. The band gaps vanish at transition strains ($\varepsilon_{trans}$). The value of compressive $\varepsilon_{trans}$ and its corresponding pressure ($P_{trans}$) are reported in Table \ref{tab:tab2}. The band gaps of MoSi$_2$P$_4$, MoSi$_2$As$_4$ and MoSi$_2$Sb$_4$ vanish around the strain of -10$\%$ or -11$\%$. On the other hand, MoSi$_2$N$_4$ demonstrates a higher bandgap, and $\varepsilon_{trans}$ is two times larger at -22$\%$ strain.
$P_{trans}$ for MoSi$_2$P$_4$, MoSi$_2$As$_4$ and MoSi$_2$Sb$_4$ is in the range of [7.3-9.1] GPa, whereas MoSi$_2$N$_4$ demonstrates much larger $P_{trans}$ of 25.3 GPa. The transition pressures for bilayers are reported from 2.18 GPa in CrSi$_2$N$_4$ to 32.04 GPa in TiSi$_2$N$_4$ \cite{zhong2021strain} that are comparable with MoSi$_2$X$_4$ monolayer. This transition pressure range also proves that applying out-of-strain to monolayer MoSi$_2$X$_4$ is feasible. The band gaps also vanish at tensile strain. MoSi$_2$Sb$_4$ demonstrates higher sensitivity to tensile strain and its band gap vanishes at lower strains. The transition strain at tensile regime increases as X atom changes from Sb to N. One can conclude from the figure that tensile transition strain of MoSi$_2$N$_4$ is much larger than 15$\%$.

\begin{figure*}
	\centering
	\includegraphics[width=0.8\linewidth]{orbitalbandstructure}
	\caption{Projected band structures for the studied materials. The d-orbital of Mo along with the p-orbital of Si and X atoms have the main contributions around the band gaps.  }
	\label{fig:orbitalbandstructure}
\end{figure*}

The energies of the K-valleys for the first and second conduction bands, and K- and $\Gamma$-points of the valence band versus vertical strain is plotted in Fig. \ref{fig:energy}. The energies of the up- and down-spins are separately shown in the figure. The strain ranges of the figures are different and the strain ranges are selected in the semiconducting phase. At large compressive strains, the energy of $\Gamma$-point of the valence band is higher than K-point and K-valley determines the VBM. By increasing strain from compressive to the tensile regime, the energy of $\Gamma$-point decreases, and the energy of K-point increases. One can find the large spin-splitting at K-point, whereas $\Gamma$-point does not demonstrate any splitting. On the other hand, two bands contribute to the conduction band. We have called the bands with low and large spin-splitting as $\alpha$-band and $\beta$-band, respectively. The $\alpha$-band contributes as CBM at compressive and small tensile strains. The larger tensile strains bring down the $\beta$-band and change the location of these two bands. Therefore, CBM is contributed by $\beta$-band at a large tensile strain regime.

The charge distribution around K- and $\Gamma$-points of the valence band and K-valleys of the first and second conduction bands in the case of MoSi$_2$P$_4$ is shown in Fig. \ref{fig:charge}. A single layer of MoSi$_2$P$_4$ can be divided into three sub-layers. The charge is significantly localized around the Mo atom in the K-point of valence and conduction bands. $\Gamma$-point of the valence band has bonding character between P and Si atoms of two sub-layers. In addition, the charge is distributed around the most external P atoms. The out-of-plane strain does not significantly affect charge distribution at this point. In the K-point of the valence band, the charge is highly distributed in the internal MoN$_2$ sub-layer and a small charge is around outer P atoms. At compressive out-of-plane strain, the charge around outer P atoms is expanded. Furthermore, the first and second conduction bands at K-valley are compared in the figure. At equilibrium, the charge is localized in the internal sub-layer for the first band, whereas the charge is distributed in the whole thickness of the compound for the second band. The compressive strain does not affect the charge distribution, whereas tensile strain changes the location and charge density of these two bands.

\begin{figure*}
	\centering
	\includegraphics[width=1.0\linewidth]{SS}
	\caption{The spin-splitting of the conduction and valence bands at the K-point. The valence band along with the first and second conduction bands are studied. }
	\label{fig:sssb}
\end{figure*}

Two-dimensional maps of the spin-splitting in the first Brillouin zone for all compounds have been shown in Fig. \ref{fig:ssc-ssv}. Both first and second conduction bands along with the valence band are indicated in the figure. The positive and negative spin-splitting around K and K'-valleys are obvious in the figures, respectively. The first conduction bands of all compounds except MoSi$_2$N$_4$ display a low spin-splitting. In the case of MoSi$_2$N$_4$, a large spin-splitting is around M- and \textGamma-valley, but these two valleys do not contribute to the CBM. One can observe that the spin-splitting of the first conduction band is approximately zero on the boundary of the first Brillouin zone for MoSi$_2$P$_4$ and MoSi$_2$As$_4$. On the contrary, the second conduction band exhibits a high spin-splitting which is located at K- and K'- valleys. MoSi$_2$N$_4$ and MoSi$_2$Sb$_4$ demonstrate a higher spin-splitting compared to the others. ُThe valence band shows a tidy spin-splitting so that the splitting is maximum around K- and K'-valleys and it vanishes elsewhere. The spin-splitting increases for heavier compounds and MoSi$_2$Sb$_4$ show the highest value. We have seen from Table \ref{tab:tab2} that the valence band of MoSi$_2$Sb$_4$ has the highest spin-splitting as high as 220 meV. The spin-splitting of the valence band is isotropic around the K-valley of MoSi$_2$P$_4$ and other materials specially MoSi$_2$Sb$_4$ demonstrate anisotropic spin-splitting. In this regard, the contour around the K-valley of MoSi$_2$P$_4$ is circular, whereas MoSi$_2$Sb$_4$ has a triangular shape.

To clarify the origin of spin-splitting in the valence and conduction bands, the projected band structure of four structures is plotted in Fig. \ref{fig:orbitalbandstructure}. As it is obvious, the first and second conduction bands at K-valley are mainly affected by the d-orbitals of molybdenum atoms. The p-orbitals of X and Si atoms have a low contribution on the first and second bands, respectively. The large spin-splitting of the second band comes from the d-orbitals, whereas the first conduction band with a high d-orbital contribution shows a small splitting. Maybe, this is related to the odd and even bands because the third one also displays a negligible spin-splitting in MoSi$_2$As$_4$ and MoSi$_2$Sb$_4$. In addition, the valence band at K- and \textGamma-points is highly formed by the d-orbitals of Mo atoms. The p-orbitals of X atoms also have a considerable influence at the K-point, whereas p-orbitals of Si atoms contribute to the \textGamma-point. Therefore, the d-orbitals of Mo atoms, the p-orbitals of X and Si atoms have the most impact on the \textGamma-point of the valence band, respectively. The large spin-splitting at valence band is originated from the d-orbitals, whereas the p-orbitals of X atoms enhance the spin-splitting for heavier compounds.

The spin-splitting of the K-point of the valence band and the first and second bands of the conduction band versus vertical strain is shown in Fig. \ref{fig:sssb}. $\lambda_{K,V}$ increases slightly when strain increases from compressive to tensile regime. The increase of the spin-splitting of the valence band with respect to strain has also observed in the group-III monochalcogenides \cite{ariapour2020strain}. One can find from the band structures and Table \ref{tab:tab2}, that the first conduction band has a negligible spin-splitting, while the second band demonstrates a considerable spin-splitting. $\lambda_{K,C1}$ increases when strain increases from compressive to tensile, whereas $\lambda_{K,C2}$ decreases. The slope of $\lambda_{K,C2}$ variation is higher in lighter compounds. As one can observe the curves of $\lambda_{K,C1}$ and $\lambda_{K,C2}$ change with each other at the tensile strain. As we have seen in Fig. \ref{fig:energy}, the $\beta$-band moves downward and becomes lower than $\alpha$-band for tensile strain. So, in the tensile regime, $\beta$-band with large spin-splitting is the first conduction band. This band changing is also observed in the figure.

\section{conclusion}
\label{Sec:Conclusions}
The electrical and spin properties of monolayer MoSi$_2$X$_4$ (X = N, P, As, and Sb) under vertical strain are investigated. MoSi$_2$N$_4$ demonstrates an indirect bandgap from $\Gamma$ at valence to K at conduction band, whereas K-valley at valence band is close to $\Gamma$-point and contributes to the valence band. On the other hand, three other compounds have direct bandgap at K-point and $\Gamma$-point at valence band is also close to K-point. Although the first conduction band demonstrates a small spin-splitting, the second conduction band has a large spin-splitting and its minimum is near to the first band.
The PDOS demonstrates that the conduction and valence bands are mostly affected by the d-orbitals of Mo atoms. X atoms have a lower contribution to the valence band. The projected band structure also demonstrates that the conduction and valence bands are mainly constructed from the d-orbitals of Mo atoms. The p-orbitals of X and Si atoms also contribute to the valence and second conduction bands, respectively.
The spin-splitting at K-point ($\lambda_{K,V}$) of valence band increases with changing X-atom from N to Sb and reaches 220 meV for MoSi$_2$Sb$_4$.
The spin-splitting at K-valley in the second conduction band ($\lambda_{K,C2}$) increases from 83 meV in MoSi$_2$P$_4$ to 151 meV for MoSi$_2$Sb$_4$. In the following, the effect of vertical strain has been investigated. The band gaps exhibit a maximum value at small tensile strains, then they decrease for compressive and larger tensile strains.
The band gaps close at compressive transition strains ($\varepsilon_{trans}$) of around -10$\%$ for all compounds except MoSi$_2$N$_4$ that its bandgap vanishes at the larger strain of about -22$\%$.
Transition pressures are also distributed from 9.1 to 7.3 GPa for MoSi$_2$P$_4$, MoSi$_2$As$_4$ and MoSi$_2$Sb$_4$, whereas MoSi$_2$N$_4$ displays a larger $P_{trans}$ of 25.3 GPa. In addition, the charge is mostly localized around Mo atoms in both valence and conduction bands. The spin-splitting around K-valley at the valence band is isotropic and circular in the case of MoSi$_2$P$_4$, whereas the anisotropy increases for heavier compounds and MoSi$_2$Sb$_4$ exhibits a triangular shape contour. The spin-splitting of K-valley of the valence band increases when strain increases from compressive to tensile regime. $\lambda_{K,C1}$ increases when strain increases from compressive to tensile, and at the same time, $\lambda_{K,C2}$ decreases. The slope of $\lambda_{K,C2}$ variation is higher for the lighter compounds. $\lambda_{K,C1}$ and $\lambda_{K,C2}$ change with each other at the tensile strain. In addition to the electrical properties, the vertical strain can also control the spin properties.

\bibliographystyle{apsrev}


\section{Introduction}\label{introduction}

There has been an increasing interest in the theory and properties of trapped
Bose-Einstein condensates at the limit of an infinite number of particles \cite{CASTIN,LIEB1,LIEB2,ERDOS1,ERDOS2,VAR1,VAR2,CEDERBAUM1,
MIX1,MIX2,MIX3,CEDERBAUM2}.
Here, one may divide the research questions into two, inter-connected groups.
The first group of research questions deals with rigorous results,
mainly proving when many-body and mean-field,
Gross-Pitaevskii theories coincide at this limit,
whereas the second group of questions deals with characterizing correlations
in a trapped Bose-Einstein condensate
based on the difference between many-body and mean-field
properties at the infinite-particle-number limit.
In \cite{LIEB1},
it has been shown that the ground-state energy per particle and density per particle
computed at the many-body level of theory coincide with the respective mean-field results.
The infinite-particle-number limit is defined
such that the interaction parameter, i.e., the product of the number of bosons times the interaction strength
(proportional to the scattering length), is held constant.
Similarly, in \cite{LIEB2} it has been shown that the reduced one-body and any reduced finite-$n$-body density matrix \cite{RDMB1,RDMB2} per particle are 100\% condensed,
and that the leading natural orbital boils down to the Gross-Pitaevskii single-particle function.
Analogous results connecting time-dependent many-body and mean-field theories are given in \cite{ERDOS1,ERDOS2}
and developments for mixtures in \cite{MIX1,MIX2,MIX3}.
In \cite{CEDERBAUM2},
the many-boson wave-function at the infinite-particle-number limit has been constructed explicitly.

The difference between many-body and mean-field theories at the limit of an infinite number of particles,
which as stated above coincide at the level of the energy, densities, and reduced density matrices per particle,
starts to show up in variances of many-particle observables \cite{VAR1,VAR2}.
Of course, the wave-functions
themselves differ and their overlap is smaller than one \cite{CEDERBAUM1}.
In evaluating the variances of many-particle observables
two-body operators emerge whose combination with
the elements of the reduced two-body density matrix can pick up even
the tiniest depletion, which always exist due to the inter-particle interaction \cite{VAR1,VAR2}.
Here, the quantitative difference between the many-body and mean-field variances
is a useful tool to benchmark many-body numerical approaches \cite{RMP},
whereas the qualitative differences serve to define and characterize the nature
of correlations in 100\% condensed bosons at the infinite-particle-number limit.

Qualitative differences between the many-body and mean-field variances per particle
depend on both the system and observable under investigation
and emerge because the 100\% condensed bosons are interacting.
In the ground-state of a one-dimensional double-well potential,
the mean-field position variance per particle increases monotonously with the interaction parameter
whereas, once about a single particle is excited outside the condensed mode,
the many-body position variance per particle starts to decrease \cite{VAR1}.
In the analogous time-dependent setup of a bosonic Josephson junction,
the mean-field variance is oscillating and bound by the size of the junction,
whereas the many-body variance increases to `sizes' several times larger \cite{VAR2}.
In two spatial dimensions additional features come out.
The position variance per particle in a thin annulus can exhibit a different dimensionality \cite{VAR3}
and both the position and momentum variances can exhibit opposite anisotropies \cite{VAR4}
when computed at the many-body and mean-field levels of theory
in an out-of-equilibrium quench dynamics.
In two spatial dimensions the many-particle variance
of the $\hat L_Z$ component of the angular-momentum operator becomes available,
and used to analyze the lack of conservation of symmetries in the mean-field dynamics \cite{VAR5}.

In the present work we analyze the many-particle position, momentum, and angular-momentum
variances of a three-dimensional anisotropic trapped Bose-Einstein condensate at
the limit of an infinite number of particles,
focusing on three-dimensional scenarios that do not have (one-dimensional and) two-dimensional analogs.
Mainly, the available permutations between the three Cartesian components of a many-particle operator,
such as the position and momentum operators,
allow one for various different
anisotropies of the respective mean-field and many-body variances than in two spatial dimensions \cite{VAR4}.
Furthermore,
anisotropy of the angular-momentum variance can only be investigated when
there is more than one component,
and this occurs with the three Cartesian components of the angular-momentum operator
in three spatial dimensions.

The structure of the paper is as follows.
In Sec.~\ref{theory} theory and definitions are developed.
In Sec.~\ref{applications} we present two applications where a common methodological line of investigation
is that the variances at the many-body level of theory can be computed analytically.
In Subsec.~\ref{position}, the anisotropy of the position and momentum variances in the out-of-equilibrium breathing dynamics of a Bose-Einstein condensate in a three-dimensional anisotropic harmonic
potential are analyzed,
and in Subsec.~\ref{angular},
a solvable model is devised which allows one to analyze
the anisotropy of the angular-momentum variances in the ground state
of interacting bosons in a three-dimensional anisotropic harmonic potential.
All quantities are computed at the infinite-particle-number limit.
Finally, we summarize in Sec.~\ref{summary}.

\section{Theory}\label{theory}

The variances per particle of a many-particle
observable $\hat O = \sum_j^N \hat o_j$
computed at the many-body (MB) and mean-field, Gross-Pitaevskii (GP) levels of theory
are connected at the
limit of an infinite number of particles
by the following relation \cite{VAR1}:
\begin{equation}\label{VAR_INF}
{\rm Var|_{MB}}(\hat O) = {\rm Var|_{GP}}(\hat O) + {\rm Var|_{correlations}}(\hat O).
\end{equation}
Here, ${\rm Var|_{MB}}(\hat O)=\lim_{N \to \infty}
\left[\frac{1}{N}(\langle\Psi|\hat O^2|\Psi\rangle-\langle\Psi|\hat O|\Psi\rangle^2)\right]$,
where $\Psi$ is the solution of the many-particle Schr\"odinger equation,
and
${\rm Var|_{GP}}(\hat O) = \langle\psi_{GP}|\hat o^2|\psi_{GP}\rangle-\langle\psi_{GP}|\hat o|\psi_{GP}\rangle^2$,
where $\psi_{GP}$ is the solution of the corresponding Gross-Pitaevskii equation.
Recall that the infinite-particle-number limit is defined
such that the interaction parameter, i.e., the product of the number of bosons times the interaction strength, is kept fixed.
The correlations term, ${\rm Var|_{correlations}}(\hat O)$, quantifying the difference between the mean-field and many-body variances, depends on the elements of the reduced two-body density matrix
where at least one of the indexes corresponds to a natural orbital higher than the condensed mode
\cite{VAR1}.
For non-interacting bosons the correlations term obviously vanishes.
As stated above,
one is interested in qualitative differences between ${\rm Var|_{GP}}(\hat O)$ and ${\rm Var|_{MB}}(\hat O)$
and their origin.

Consider a Bose-Einstein condensate for which the many-body variances of the three Cartesian components
of, say, the position operator are different and satisfy, without loss of generality,
the anisotropy
\begin{subequations}\label{VAR_ANI_3D}
\begin{equation}\label{VAR_ANI_3D_MB}
{\rm Var|_{MB}}(\hat X) > {\rm Var|_{MB}}(\hat Y) > {\rm Var|_{MB}}(\hat Z).
\end{equation}
We define the following classification with respect to the possible different anisotropies of the
respective mean-field position variances:
\begin{eqnarray}\label{VAR_ANI_3D_GP}
\{1\}: \quad & & \Big\{ {\rm Var|_{GP}}(\hat X) > {\rm Var|_{GP}}(\hat Y) > {\rm Var|_{GP}}(\hat Z) \Big\}, \nonumber \\
\{1,2\}: \quad & & \Big\{ {\rm Var|_{GP}}(\hat Y) > {\rm Var|_{GP}}(\hat X) > {\rm Var|_{GP}}(\hat Z); \nonumber \\
         \quad & & \ \ {\rm Var|_{GP}}(\hat X) > {\rm Var|_{GP}}(\hat Z) > {\rm Var|_{GP}}(\hat Y); \nonumber \\
         \quad & & \ \ {\rm Var|_{GP}}(\hat Z) > {\rm Var|_{GP}}(\hat Y) > {\rm Var|_{GP}}(\hat X) \Big\}, \nonumber \\
\{1,2,3\}: \quad & & \Big\{ {\rm Var|_{GP}}(\hat Y) > {\rm Var|_{GP}}(\hat Z) > {\rm Var|_{GP}}(\hat X); \nonumber \\
         \quad & & \ \ {\rm Var|_{GP}}(\hat Z) > {\rm Var|_{GP}}(\hat X) > {\rm Var|_{GP}}(\hat Y) \Big\}. \
\end{eqnarray}
\end{subequations}
Naturally, the classification (\ref{VAR_ANI_3D_GP}) follows the classes of the $S_3$ permutation group denoted by
$\{1\}$, $\{1,2\}$, and $\{1,2,3\}$.
If the mean-field variances exhibit anisotropy other than the anisotropy of the respective many-body variances,
i.e., the ordering of the former does not belong to the class $\{1\}$,
we may interpret that the mean-field and many-body morphologies
of the Bose-Einstein condensate with respect to the operators under investigation are distinct.
In the present work we investigate manifestations of definition (\ref{VAR_ANI_3D})
utilizing the many-particle position ($\hat X, \hat Y, \hat Z$), momentum ($\hat P_X, \hat P_Y, \hat P_Z$),
and angular-momentum ($\hat L_X, \hat L_Y, \hat L_Z$) operators
for classifying the morphology of 100\% condensed trapped bosons at the limit of an infinite number of particles.

\section{Applications}\label{applications}

\subsection{Position and momentum variances in an out-of-equilibrium dynamics of a
three-dimensional trapped Bose-Einstein condensate}\label{position}

Consider $N$ structureless bosons trapped in a three-dimensional anisotropic harmonic potential
and interacting by a general two-body interaction $\hat W({\bf r}-{\bf r}')$.
The frequencies of the trap satisfy, without loss of generality, $\omega_x < \omega_y < \omega_z$.
We work with dimensionless quantities, $\hbar=m=1$.
Using Jacobi coordinates,
${\bf Q}_k = \frac{1}{\sqrt{k(k+1)}} \sum_{j=1}^k ({\bf r}_{k+1} - {\bf r}_j), \ k=1,\ldots,N-1$,
${\bf Q}_N = \frac{1}{\sqrt{N}} \sum_{j=1}^N {\bf r}_j$,
where ${\bf r}_1,\ldots,{\bf r}_N$ are the coordinates in the laboratory frame,
the Hamiltonian can be written as:
\begin{eqnarray}\label{HAM_3D}
\hat H({\bf Q}_1,\ldots,{\bf Q}_N) = -\frac{1}{2} \frac{\partial^2}{\partial {\bf Q}^2_N} +
\frac{1}{2} \left( \omega_x^2 Q^2_{N,x} + \omega_y^2 Q^2_{N,y} + \omega_z^2 Q^2_{N,z} \right) +
\hat H_{rel}({\bf Q}_1,\ldots,{\bf Q}_{N-1}).
\end{eqnarray}
The `relative' Hamiltonian $\hat H_{rel}$ collects all terms depending on the relative coordinates ${\bf Q}_1,\ldots,{\bf Q}_{N-1}$,
and ${\bf Q}_k=(Q_{k,x},Q_{k,y},Q_{k,z})$.
Suppose now that the bosons are prepared in the ground state of the non-interacting system.
The ground-state is separable in the Jacoby coordinates and reads
$\Phi({\bf Q}_1,\ldots,{\bf Q}_N) = \left(\frac{\omega_x}{\pi}\right)^{\frac{N}{4}}\left(\frac{\omega_y}{\pi}\right)^{\frac{N}{4}}
\left(\frac{\omega_z}{\pi}\right)^{\frac{N}{4}}
e^{-\frac{1}{2}\sum_{k=1}^N \left(\omega_x Q_{k,x}^2 + \omega_y Q_{k,y}^2 + \omega_z Q_{k,z}^2\right)}$,
where the relation $\sum_{j=1}^N {\bf r}_j^2 = \sum_{k=1}^N {\bf Q}_k^2$ connecting the laboratory
and Jacoby coordinates is used.
The solution of the time-dependent many-boson Schr\"odinger equation,
$\hat H \Psi(t) = i \frac{\partial \Psi(t)}{\partial t}$,
where $\Phi$ is the initial condition,
reads $\Psi(t) = e^{-i \hat H t} \Phi$.
Consequently,
because of the center-of-mass separability of the Hamiltonian $\hat H$
and of $\Phi$,
the position and momentum variances per particle
of the time-dependent state $\Psi(t)$,
for a general inter-particle interaction $\hat W({\bf r}-{\bf r}')$,
are those of the static, non-interacting system:
\begin{eqnarray}\label{POS_MOM_VAR}
& &
{\rm Var|_{MB}}(\hat X) = \frac{1}{2\omega_x}, \quad
{\rm Var|_{MB}}(\hat Y) = \frac{1}{2\omega_y}, \quad
{\rm Var|_{MB}}(\hat Z) = \frac{1}{2\omega_z}, \quad \nonumber \\
& &
{\rm Var|_{MB}}(\hat P_X) = \frac{\omega_x}{2}, \quad
{\rm Var|_{MB}}(\hat P_Y) = \frac{\omega_y}{2}, \quad
{\rm Var|_{MB}}(\hat P_Z) = \frac{\omega_z}{2}. \quad \
\end{eqnarray}
In other words,
the anisotropies
${\rm Var|_{MB}}(\hat X) > {\rm Var|_{MB}}(\hat Y) > {\rm Var|_{MB}}(\hat Z)$
of the position operator
and
${\rm Var|_{MB}}(\hat P_Z) > {\rm Var|_{MB}}(\hat P_Y) > {\rm Var|_{MB}}(\hat P_X)$
of the momentum operator,
when computed at the many-body level of theory,
hold for all times during the out-of-equilibrium dynamics,
see the constant-value (dashed) curves in Figs.~\ref{f1} and \ref{f2}.

What happens at the Gross-Pitaevskii level of theory?
Can the mean-field variances have different orderings than the many-body variances,
i.e., belong to other anisotropy classes based on the $S_3$ permutation group,
see (\ref{VAR_ANI_3D_GP}),
than to $\{1\}$?
If yes, then why and how?
The Gross-Pitaevskii or non-linear Schr\"odinger equation
is given by\break\hfill
$\left[-\frac{1}{2}\frac{\partial^2}{\partial {\bf r}^2} +
\frac{1}{2} \left( \omega_x^2 x^2 + \omega_y^2 y^2 + \omega_z^2 z^2 \right)
+ g|\psi_{GP}({\bf r},t)|^2 \right] \psi_{GP}({\bf r},t) =
i\frac{\partial \psi_{GP}({\bf r},t)}{\partial t}$,
where $g=4\pi N a_s$ is the coupling constant and $a_s$
the $s$-wave scattering length of the
above two-body interaction $\hat W({\bf r}-{\bf r}')$.
The initial condition, as above, is the ground state of the non-interacting system,
$\psi_{GP}({\bf r},0) = \left(\frac{\omega_x}{\pi}\right)^{\frac{1}{4}}\left(\frac{\omega_y}{\pi}\right)^{\frac{1}{4}}
\left(\frac{\omega_z}{\pi}\right)^{\frac{1}{4}}
e^{-\frac{1}{2}\left(\omega_x x^2 + \omega_y y^2 + \omega_z z^2\right)}$.
The Gross-Pitaevskii equation does not maintain the center-of-mass separability
of the initial condition because of its non-linear term,
which, therefore, can lead to variations of the position and momentum variances
when computed at the mean-field level of theory.

\begin{figure}[!]
\begin{center}
\vglue -2.0 truecm
\hglue -1.0 truecm
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig1a.eps}
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig1b.eps}
\vglue -0.25 truecm
\hglue -1.0 truecm
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig1c.eps}
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig1d.eps}
\vglue -0.25 truecm
\hglue -1.0 truecm
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig1e.eps}
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig1f.eps}
\vglue -0.25 truecm
\hglue -1.0 truecm
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig1g.eps}
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig1h.eps}
\end{center}
\vglue -0.275 truecm
\caption{Many-particle position ($\hat X$, $\hat Y$, and $\hat Z$; in red, green, and blue)
variance per particle as a function of time
computed at the limit of an infinite number of particles
within many-body (dashed lines) and mean-field (solid lines) levels of theory in an interaction-quench scenario.
The harmonic trap is 10\% anisotropic in panels (a), (c), (e), (g)
and 20\% anisotropic in panels (b), (d), (f), (h).
The coupling constant $g$ is indicated in each panel.
Different anisotropy classes of the position variance emerge with time.
See the text for more details.
The quantities shown are dimensionless.}
\label{f1}
\end{figure}

\begin{figure}[!]
\begin{center}
\vglue -2.0 truecm
\hglue -1.0 truecm
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig2a.eps}
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig2b.eps}
\vglue -0.25 truecm
\hglue -1.0 truecm
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig2c.eps}
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig2d.eps}
\vglue -0.25 truecm
\hglue -1.0 truecm
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig2e.eps}
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig2f.eps}
\vglue -0.25 truecm
\hglue -1.0 truecm
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig2g.eps}
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig2h.eps}
\end{center}
\vglue -0.275 truecm
\caption{Many-particle momentum ($\hat P_X$, $\hat P_Y$, and $\hat P_Z$; in red, green, and blue)
variance per particle as a function of time
computed at the infinite-particle-number limit
within many-body (dashed lines) and mean-field (solid lines) levels of theory in an interaction-quench scenario.
The harmonic trap is 10\% anisotropic in panels (a), (c), (e), (g)
and 20\% anisotropic in panels (b), (d), (f), (h).
The coupling constant $g$ is indicated in each panel.
Different anisotropy classes of the momentum variance emerge with time.
See the text for more details.
The quantities shown are dimensionless.}
\label{f2}
\end{figure}

Figs.~\ref{f1} and \ref{f2} display the Gross-Pitaevskii dynamics of
the position and momentum variances per particle, respectively,
for four coupling constants,
$g=0.18, 9.0, 18.0$, and $27.0$.
To integrate the three-dimensional Gross-Pitaevskii equation we
use a box of size $[-10,10) \times [-10,10) \times [-10,10)$, a Fourier-discrete-variable-representation with
$128^3$ grid points and periodic boundary conditions,
and the numerical implementation embedded in \cite{GP_in_MCTDHB_Package}.
The dynamics are computed for the four
coupling constants $g$ and depicted
by the oscillating (solid) curves in Figs.~\ref{f1} and \ref{f2}.
The left columns are for a 10\% anisotropy of the harmonic trap,
i.e., $\omega_x = 0.9$,
$\omega_y = 1.0$, and
$\omega_z = 1.1$,
and the right columns are for a 20\% anisotropy of the harmonic trap,
namely, $\omega_x = 0.8$,
$\omega_y = 1.0$, and
$\omega_z = 1.2$.
We remark that the expectation values per particle of the position ($\hat X, \hat Y, \hat Z$)
and momentum ($\hat P_X, \hat P_Y, \hat P_Z$) operators
computed at the mean-field and many-body levels of theory
coincide at the limit of an infinite number of particles
and are all equal to zero in the present scenario.

For the smallest coupling constant, $g=0.18$,
we see that the mean-field variances
oscillate with very small amplitudes around the respective constant values of the many-body variances.
This means that the mean-field anisotropy of the position variance,
${\rm Var|_{GP}}(\hat X) > {\rm Var|_{GP}}(\hat Y) > {\rm Var|_{GP}}(\hat Z)$,
and its many-body anisotropy,
${\rm Var|_{MB}}(\hat X) > {\rm Var|_{MB}}(\hat Y) > {\rm Var|_{MB}}(\hat Z)$,
are alike.
A similar situation is found for the momentum variance,
namely,
that the mean-field momentum anisotropy,
${\rm Var|_{GP}}(\hat P_Z) > {\rm Var|_{GP}}(\hat P_Y) > {\rm Var|_{GP}}(\hat P_X)$,
and the many-body anisotropy,
${\rm Var|_{MB}}(\hat P_Z) > {\rm Var|_{MB}}(\hat P_Y) > {\rm Var|_{MB}}(\hat P_X)$,
are the same.
Consequently,
we may conclude that
for small coupling constants the anisotropy class of the position operator is $\{1\}$
and, likewise, the anisotropy class of the momentum operator is $\{1\}$,
see (\ref{VAR_ANI_3D}).

The situation becomes more interesting for the larger coupling constants,
$g=9.0, 18.0$, and $27.0$.
We begin with the position variances, Fig.~\ref{f1}.
The variances are found to oscillate prominently,
with much larger amplitudes than for $g=0.18$,
and, subsequently, to cross each other.
There are three ingredients that enable and govern this crossing dynamics.
The first, is that the amplitudes of oscillations
of ${\rm Var|_{GP}}(\hat X)$, ${\rm Var|_{GP}}(\hat Y)$, and ${\rm Var|_{GP}}(\hat Z)$
are slightly different already at short times,
with the former being the larger
and the latter being the smaller
(more prominent for 20\% than for 10\% trap anisotropy).
The second, is that the respective frequencies of oscillations
are also slightly different at short times,
with the former being the smaller
and the latter being the larger.
Both features correlate with the
ordering of the frequencies of the trap,
$\omega_x < \omega_y < \omega_z$.
The third ingredient is that the three Cartesian components are coupled to
each other during the dynamics,
what impacts the oscillatory pattern at intermediate and later times
(more prominent for 10\% than for 20\% trap anisotropy, see Fig.~\ref{f1}).

Combing the above,
we find for the 10\% trap anisotropy that
around $t=4.0$
${\rm Var|_{GP}}(\hat Z) > {\rm Var|_{GP}}(\hat Y) > {\rm Var|_{GP}}(\hat X)$ takes place,
around $t=6.0$
${\rm Var|_{GP}}(\hat X) > {\rm Var|_{GP}}(\hat Y) > {\rm Var|_{GP}}(\hat Z)$ holds (again),
and
around $t=7.0$
${\rm Var|_{GP}}(\hat Z) > {\rm Var|_{GP}}(\hat Y) > {\rm Var|_{GP}}(\hat X)$ occurs.
In other words,
the anisotropy class of the position variance starts as $\{1\}$ for $t=0$,
changes to $\{1,2\}$ around $t=4$,
is back to $\{1\}$ around $t=6$,
and becomes $\{1,2,3\}$ around $t=7$.
Furthermore,
this pattern is found to be robust for different, increasing coupling constants, see Fig.~\ref{f1}.
For 20\% trap anisotropy we find a different crossing patten of the position variances.
The anisotropy class begins as $\{1\}$ for $t=0$,
changes at around $t=3.0$ to $\{1,2\}$,
and immediately after,
at around $t=3.75$, it is $\{1,2,3\}$.
Now, around $t=7.0$ there is a broad regime
of anisotropy class $\{1,2\}$.
Another difference of a geometrical origin between the
dynamics in the 20\% and 10\% trap anisotropies
can be seen for $g=9$, see Fig.~\ref{f1}c,d.
Here, the coupling constant is sufficiently large
to lead to crossing of all position variances for the 10\% anisotropy trap,
and, consequently,
to the position anisotropy class $\{1,2,3\}$
(at around $t=7$).
On the other hand,
for the 20\% anisotropy trap
the coupling constant is just short
of allowing all position variances to cross each other
and, clearly, the anisotropy class $\{1,2,3\}$ cannot occur
(as it happens at around $t=3.75$ for the further larger coupling constants, $g=18$ and $27$).
All in all,
we have demonstrated in a rather common (out-of-equilibrium quench) scenario
the emergence of anisotropy classes other than $\{1\}$,
i.e., $\{1,2\}$ and $\{1,2,3\}$,
for the position operator
of a Bose-Einstein condensate
at the infinite-particle-number limit.

The results for the momentum variances per particle, see Fig.~\ref{f2},
follow similar and corresponding trends as those for the position operator,
albeit the crossings of the respective momentum curves
take place during slightly narrower time windows than for the position operator
for the parameters used.
Thus, for the 10\% trap anisotropy we have
${\rm Var|_{GP}}(\hat P_X) > {\rm Var|_{GP}}(\hat P_Y) > {\rm Var|_{GP}}(\hat P_Z)$ around $t=4.25$,
${\rm Var|_{GP}}(\hat P_Z) > {\rm Var|_{GP}}(\hat P_Y) > {\rm Var|_{GP}}(\hat P_X)$ in the vicinity of $t=6.0$,
and
${\rm Var|_{GP}}(\hat P_Y) > {\rm Var|_{GP}}(\hat P_X) > {\rm Var|_{GP}}(\hat P_Z)$ around $t=7.25$.
Therefore,
the anisotropy class of the momentum variance starts at $\{1\}$ for $t=0$,
turns to $\{1,2\}$ around $t=4.25$,
returns to $\{1\}$ for a wider time window around $t=6.0$,
and changes to $\{1,2,3\}$ around $t=7.25$.
For the 20\% trap anisotropy we find,
starting from the anisotropy class $\{1\}$ for $t=0$,
the class $\{1,2,3\}$ at around $t=4.0$,
the class $\{1,2\}$ at around $t=4.75$,
again the class $\{1\}$ at around $t=6$,
and once more the class $\{1,2,3\}$ at around $t=9.5$.
As for the position variance,
we find the pattern to be robust for different, increasing coupling constants, see Fig.~\ref{f2}.
Furthermore,
the above-discussed difference of a geometrical origin between the
position-variance dynamics in the 20\% and 10\% trap anisotropies for $g=9$
emerges also for the momentum variance, see Fig.~\ref{f2}c,d.
Here, the coupling constant is sufficiently large
to lead to crossing of all momentum variances for the 10\% trap anisotropy,
but not for the 20\% trap anisotropy.
As a result,
the former system exhibits also the anisotropy class $\{1,2,3\}$ for the momentum variance,
whereas the latter only the $\{1,2\}$ anisotropy class.
Summarizing,
we have demonstrated in a simple scenario, of an out-of-equilibrium breathing dynamics,
the emergence of anisotropy classes other than $\{1\}$,
namely, $\{1,2\}$ and $\{1,2,3\}$,
for the many-particle position as well as many-particle momentum operators
of a trapped Bose-Einstein condensate condensate at the limit of an infinite number of particles.

\subsection{Angular-momentum variance in the ground state of a
three-dimensional trapped Bose-Einstein condensate}\label{angular}

The possibility to learn on the relations governing correlations
and variance anisotropy between the different components
of the angular-momentum operator opens up only in three spatial dimensions.
Here, in the context of the present work, the challenge is to find a many-particle model
where angular-momentum properties can be treated analytically
at the many-body level of theory and
in the limit of an infinite number of particles.
Such a model is the three-dimensional anisotropic harmonic-interaction model,
and the results presented below build on and clearly extends
the investigation of the two-dimensional anisotropic harmonic-interaction model
reported in \cite{VAR6}.
The harmonic-interaction model has been used quite extensively including to model Bose-Einstein condensates
\cite{H0,H1,H2,H3,H4,H5,H6,H7,H8,H9,H10,H11,H12,H13,H14}.
Finally, and as a bonus, we mention that the
three-dimensional anisotropic harmonic-interaction model
can also be solved analytically at the mean-field level of theory, which is useful for the analysis.

In the laboratory frame
the three-dimensional anisotropic harmonic-interaction model reads:
$\hat H({\bf r}_1,\ldots,{\bf r}_N) =
\sum_{j=1}^N \left[-\frac{1}{2}\frac{\partial^2}{\partial {\bf r}_j^2} + \frac{1}{2}
\left(\omega_x^2 x_j^2 + \omega_y^2 y_j^2 + \omega_z^2 z_j^2 \right)\right]
+ \lambda_0 \sum_{1 \le j < k}^N ({\bf r}_j-{\bf r}_k)^2$, i.e.,
it is obtained from the Hamiltonian (\ref{HAM_3D}) when the two-body interaction is $\hat W({\bf r}-{\bf r}') = \lambda_0 ({\bf r}-{\bf r}')^2$.
Then, the `relative' Hamiltonian is given explicitly by
\begin{eqnarray}\label{HAM_3D_HIM}
& & \hat H_{rel}({\bf Q}_1,\ldots,{\bf Q}_{N-1}) = \\
& & \quad = \sum_{k=1}^{N-1} \left[
-\frac{1}{2} \frac{\partial^2}{\partial {\bf Q}^2_k} +
\frac{1}{2} \left( \sqrt{\omega_x^2+2N\lambda_0} Q^2_{k,x} + \sqrt{\omega_y^2+2N\lambda_0} Q^2_{k,y} + \sqrt{\omega_z^2+2N\lambda_0} Q^2_{k,z} \right) \right]. \nonumber \
\end{eqnarray}
The many-body ground state of $\hat H$ is readily obtained and given by
\begin{eqnarray}\label{HIM_3D_GS}
& & \Psi({\bf Q}_1,\ldots,{\bf Q}_N) = \nonumber \\
& & = \left(\frac{\omega_x}{\pi}\right)^{\frac{1}{4}}
\left(\frac{\omega_y}{\pi}\right)^{\frac{1}{4}}
\left(\frac{\omega_z}{\pi}\right)^{\frac{1}{4}}
\left(\frac{\sqrt{\omega_x^2+2N\lambda_0}}{\pi}\right)^{\frac{N-1}{4}}
\left(\frac{\sqrt{\omega_y^2+2N\lambda_0}}{\pi}\right)^{\frac{N-1}{4}}
\left(\frac{\sqrt{\omega_z^2+2N\lambda_0}}{\pi}\right)^{\frac{N-1}{4}}
\times \nonumber \\
& &
e^{-\frac{1}{2} \left(\sqrt{\omega_x^2+2N\lambda_0} \sum_{k=1}^{N-1} Q_{k,x}^2 + \omega_x Q_{N,x}^2 \right)}
e^{-\frac{1}{2} \left( \sqrt{\omega_y^2+2N\lambda_0} \sum_{k=1}^{N-1} Q_{k,y}^2 + \omega_y Q_{N,y}^2 \right)}
e^{-\frac{1}{2} \left( \sqrt{\omega_z^2+2N\lambda_0} \sum_{k=1}^{N-1} Q_{k,z}^2 + \omega_z Q_{N,z}^2 \right)}.
\nonumber \\ \
\end{eqnarray}
As states above,
it is also possible to solve analytically
the three-dimensional anisotropic harmonic-interaction model
at the mean-field level of theory by generalizing \cite{VAR6,H3}.
The final result for
the mean-field solution of the ground state reads
\begin{eqnarray}\label{HIM_3D_GP}
& & \psi_{GP}({\bf r}) =
\left(\frac{\sqrt{\omega_x^2+2\Lambda}}{\pi}\right)^{\frac{1}{4}}
\left(\frac{\sqrt{\omega_y^2+2\Lambda}}{\pi}\right)^{\frac{1}{4}}
\left(\frac{\sqrt{\omega_z^2+2\Lambda}}{\pi}\right)^{\frac{1}{4}} \times \nonumber \\
& & \times
e^{-\frac{1}{2}\sqrt{\omega_x^2+2\Lambda} \, x^2}
e^{-\frac{1}{2}\sqrt{\omega_y^2+2\Lambda} \, y^2}
e^{-\frac{1}{2}\sqrt{\omega_z^2+2\Lambda} \, z^2},
\end{eqnarray}
where $\Lambda=\lambda_0(N-1)$ is the interaction parameter.
For reference,
$\psi_{GP}({\bf r})$ solves the Gross-Pitaevskii equation
$\left[-\frac{1}{2}\frac{\partial^2}{\partial {\bf r}^2} +
\frac{1}{2} \left( \omega_x^2 x^2 + \omega_y^2 y^2 + \omega_z^2 z^2 \right) +
\Lambda \int d{\bf r}' |\psi_{GP}({\bf r}')|^2 ({\bf r}-{\bf r}')^2 \right] \psi_{GP}({\bf r}) =
\mu \psi_{GP}({\bf r})$,
where $\mu$ is the chemical potential.
Note that both many-body and mean-field solutions can be written as products of
the respective solutions in one dimension along the $x$, $y$, and $z$ directions.

Before we arrive at the angular-momentum variances and for our needs, see below,
we make a stopover and compute the position and momentum variances per particle in the model.
At the many-body level we obviously have the result (\ref{POS_MOM_VAR}),
since for the interacting ground-state the center-of-mass is separable and, hence,
the position and momentum variances are independent of the two-body interaction.
At the mean-field level we readily find from (\ref{HIM_3D_GP}) the result
\begin{eqnarray}\label{POS_MOM_VAR_GP}
& &
{\rm Var|_{GP}}(\hat X) = \frac{1}{2\sqrt{\omega_x^2+2\Lambda}}, \quad
{\rm Var|_{GP}}(\hat Y) = \frac{1}{2\sqrt{\omega_y^2+2\Lambda}}, \quad
{\rm Var|_{GP}}(\hat Z) = \frac{1}{2\sqrt{\omega_z^2+2\Lambda}}, \quad \nonumber \\
& &
{\rm Var|_{GP}}(\hat P_X) = \frac{\sqrt{\omega_x^2+2\Lambda}}{2}, \quad
{\rm Var|_{GP}}(\hat P_Y) = \frac{\sqrt{\omega_y^2+2\Lambda}}{2}, \quad
{\rm Var|_{GP}}(\hat P_Z) = \frac{\sqrt{\omega_z^2+2\Lambda}}{2}. \quad \
\end{eqnarray}
The mean-field variances (\ref{POS_MOM_VAR_GP})
depend on the interaction parameter $\Lambda$,
unlike the respective many-body variances (\ref{POS_MOM_VAR}).
It turns out that this property would be instrumental when analyzing
the anisotropy of the angular-momentum variance below.
We briefly comment on the anisotropies of the position and momentum variances
in the model.
Comparing the mean-field (\ref{POS_MOM_VAR_GP}) and
many-body (\ref{POS_MOM_VAR}) variances per particle
we find that the former belong to the anisotropy class $\{1\}$
independently of the interaction parameter $\Lambda$
both for the position and momentum operators.
For the mean-field variance of the ground state at the infinite-particle-number limit
to belong to an anisotropy class other than $\{1\}$,
one would have to go beyond the simple single-well geometry,
see the anisotropy of the position variance in a double-well potential
in two spatial dimensions \cite{VAR7}.

We can now move to the expressions for the angular-momentum variances
at the limit of an infinite number of particles,
by generalizing results obtained in two spatial dimensions \cite{VAR6} two three spatial dimensions.
The calculation at the mean-field level using (\ref{HIM_3D_GP}) readily gives
\begin{eqnarray}\label{VAR_Lx_Ly_Lz_GP}
& &
{\rm Var|_{GP}}(\hat L_X) = \frac{1}{4} \frac{\left(\sqrt{\omega_y^2+2\Lambda}-\sqrt{\omega_z^2+2\Lambda}\right)^2}{\sqrt{\omega_y^2+2\Lambda} \sqrt{\omega_z^2+2\Lambda}}
= \frac{1}{4}\frac{\left(\frac{\sqrt{\omega_y^2+2\Lambda}}{\sqrt{\omega_z^2+2\Lambda}}-1\right)^2}{\frac{\sqrt{\omega_y^2+2\Lambda}}{\sqrt{\omega_z^2+2\Lambda}}}, \nonumber \\
& &
{\rm Var|_{GP}}(\hat L_Y) = \frac{1}{4} \frac{\left(\sqrt{\omega_z^2+2\Lambda}-\sqrt{\omega_x^2+2\Lambda}\right)^2}{\sqrt{\omega_z^2+2\Lambda} \sqrt{\omega_x^2+2\Lambda}}
= \frac{1}{4}\frac{\left(\frac{\sqrt{\omega_z^2+2\Lambda}}{\sqrt{\omega_x^2+2\Lambda}}-1\right)^2}{\frac{\sqrt{\omega_z^2+2\Lambda}}{\sqrt{\omega_x^2+2\Lambda}}}, \nonumber \\
& &
{\rm Var|_{GP}}(\hat L_Z) = \frac{1}{4} \frac{\left(\sqrt{\omega_x^2+2\Lambda}-\sqrt{\omega_y^2+2\Lambda}\right)^2}{\sqrt{\omega_x^2+2\Lambda} \sqrt{\omega_y^2+2\Lambda}}
= \frac{1}{4}\frac{\left(\frac{\sqrt{\omega_x^2+2\Lambda}}{\sqrt{\omega_y^2+2\Lambda}}-1\right)^2}{\frac{\sqrt{\omega_x^2+2\Lambda}}{\sqrt{\omega_y^2+2\Lambda}}}. \
\end{eqnarray}
In the absence of interaction these
expressions boil down, respectively, to
$\frac{1}{4}\frac{\left(\frac{\omega_y}{\omega_z}-1\right)^2}{\frac{\omega_y}{\omega_z}}$,
$\frac{1}{4}\frac{\left(\frac{\omega_z}{\omega_x}-1\right)^2}{\frac{\omega_z}{\omega_x}}$,
and
$\frac{1}{4}\frac{\left(\frac{\omega_x}{\omega_y}-1\right)^2}{\frac{\omega_x}{\omega_y}}$,
the angular-momentum variances of a single particle in a three-dimensional anisotropic harmonic potential.
We see that for non-interacting particles and at the mean-field level
the angular-momentum variances per particle depend on the ratios of frequencies,
not on their absolute magnitudes.
In the first case these are the bare frequencies of the harmonic trap
whereas in the second case these are the interaction-dressed frequencies (\ref{HIM_3D_GP})
resulting from the non-linear term.

The computation of the many-body variances is lengthier and using \cite{VAR6} the final expressions for
the correlations terms (\ref{VAR_INF}) are
\begin{eqnarray}\label{VAR_Lx_Ly_Lz_CORR}
& & {\rm Var|_{correlations}}(\hat L_X) =
\frac{1}{4} \frac{\left(\sqrt{\omega_y^2+2\Lambda}-\sqrt{\omega_z^2+2\Lambda}\right)^2}{\sqrt{\omega_y^2+2\Lambda} \sqrt{\omega_z^2+2\Lambda}}
\left(\sqrt{1+\frac{2\Lambda}{\omega_y^2}} - 1\right)\left(\sqrt{1+\frac{2\Lambda}{\omega_z^2}} - 1\right), \nonumber \\
& & {\rm Var|_{correlations}}(\hat L_Y) =
\frac{1}{4} \frac{\left(\sqrt{\omega_z^2+2\Lambda}-\sqrt{\omega_x^2+2\Lambda}\right)^2}{\sqrt{\omega_z^2+2\Lambda} \sqrt{\omega_x^2+2\Lambda}}
\left(\sqrt{1+\frac{2\Lambda}{\omega_z^2}} - 1\right)\left(\sqrt{1+\frac{2\Lambda}{\omega_x^2}} - 1\right), \nonumber \\
& & {\rm Var|_{correlations}}(\hat L_Z) =
\frac{1}{4} \frac{\left(\sqrt{\omega_x^2+2\Lambda}-\sqrt{\omega_y^2+2\Lambda}\right)^2}{\sqrt{\omega_x^2+2\Lambda} \sqrt{\omega_y^2+2\Lambda}}
\left(\sqrt{1+\frac{2\Lambda}{\omega_x^2}} - 1\right)\left(\sqrt{1+\frac{2\Lambda}{\omega_y^2}} - 1\right). \nonumber \\
\end{eqnarray}
Hence, adding (\ref{VAR_Lx_Ly_Lz_GP}) and (\ref{VAR_Lx_Ly_Lz_CORR}) we readily have from (\ref{VAR_INF})
the many-body variances per particle at the infinite-particle-number limit,
${\rm Var|_{MB}}(\hat L_X)$,
${\rm Var|_{MB}}(\hat L_Y)$,
and
${\rm Var|_{MB}}(\hat L_Z)$.
We remark that the expectation values per particle of the angular-momentum operator ($\hat L_X, \hat L_Y, \hat L_Z$),
as well as the respective expectation values of the position and momentum operators,
computed at the mean-field and many-body levels of theory
coincide at the limit of an infinite number of particles
and are all equal to zero in the ground state.

We investigate and discuss an example.
Let the frequencies of the three-dimensional anisotropic harmonic trap be
$\omega_x = 0.7$,
$\omega_y = 5.0$,
and
$\omega_z = 10.5$.
Their ratios from large to small are:
$\frac{\omega_z}{\omega_x} = 15$,
$\frac{\omega_y}{\omega_x} = 7\frac{1}{7}$,
and
$\frac{\omega_z}{\omega_y} = 2\frac{1}{10}$.
Then,
the values of the angular-momentum variances per particle
at zero interaction parameter, $\Lambda=0$, are given from large to small by
${\rm Var|_{GP}}(\hat L_Y) = {\rm Var|_{MB}}(\hat L_Y) = \frac{7^2}{15} \approx 3.267$,
${\rm Var|_{GP}}(\hat L_Z) = {\rm Var|_{MB}}(\hat L_Z) = \frac{43^2}{1400} \approx 1.321$,
and
${\rm Var|_{GP}}(\hat L_X) = {\rm Var|_{MB}}(\hat L_X) = \frac{11^2}{840} \approx 0.144$.
Indeed, as the ratio of frequencies with respect to two axes is bigger,
the corresponding angular-momentum variance per particle with respect to the third axis is larger,
and vise versa.

What happens as the interaction sets in?
Fig.~\ref{f3}a depicts the many-body and mean-field
angular-momentum variances as a function of the interaction parameter parameter $\Lambda$.
We examine positive values of $\Lambda$ which correspond to the attractive sector of the
harmonic-interaction model, see, e.g., \cite{VAR6,H3,H7}.
Let us analyze the observations.
With increasing interaction parameter the density narrows, along the $x$, $y$, and $z$ directions.
This is clear because the interaction between particles is attractive,
and is manifested by the monotonously decreasing values of the
position variances per particle (\ref{POS_MOM_VAR_GP}).
Furthermore,
the density becomes less anisotropic,
because the ratios of the dressed frequencies
$\frac{\sqrt{\omega_z^2+2\Lambda}}{\sqrt{\omega_x^2+2\Lambda}}$,
$\frac{\sqrt{\omega_y^2+2\Lambda}}{\sqrt{\omega_x^2+2\Lambda}}$,
and
$\frac{\sqrt{\omega_z^2+2\Lambda}}{\sqrt{\omega_y^2+2\Lambda}}$
monotonously decreases with increasing $\Lambda$.
Consequently, the angular-momentum variances per particle decrease with the interaction parameter as well,
see Fig.~\ref{f3}a.
The mean-field angular-momentum variances (\ref{VAR_Lx_Ly_Lz_GP}) are monotonously decreasing
because of the just-described decreasing ratios of the dressed frequencies.
The many-body angular-momentum variances are decreasing,
at least for the values of interaction parameters studied here,
because the positive-value correlations terms (\ref{VAR_Lx_Ly_Lz_CORR})
grow slower
than the mean-field angular-momentum variances decrease with $\Lambda$.

All in all, the anisotropy of the angular-momentum variance can now be determined.
We find the anisotropy
${\rm Var|_{MB}}(\hat L_Y) > {\rm Var|_{MB}}(\hat L_Z) > {\rm Var|_{MB}}(\hat L_X)$
to hold for all interaction parameters at the many-body level of theory.
At the mean-field level of theory we find the same anisotropy, namely,
${\rm Var|_{GP}}(\hat L_Y) > {\rm Var|_{GP}}(\hat L_Z) > {\rm Var|_{GP}}(\hat L_X)$,
to hold for small interaction parameters.
But then, at just about $\Lambda=5.0$ the mean-field anisotropy changes to
${\rm Var|_{GP}}(\hat L_Y) > {\rm Var|_{GP}}(\hat L_X) > {\rm Var|_{GP}}(\hat L_Z)$,
and this anisotropy continues for larger interaction parameters.
Hence, we have found that
the anisotropy of the angular momentum operator
in the ground state of the
three-dimensional anisotropic harmonic-interaction model
at the infinite-particle-number limit
changes as a function of the interaction parameter
from the anisotropy class $\{1\}$ to $\{1,2\}$,
see Fig.~\ref{f3}a.

\begin{figure}[!]
\begin{center}
\hglue -1.0 truecm
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig3a.eps}
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig3b.eps}
\hglue -1.0 truecm
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig3c.eps}
\includegraphics[width=0.345\columnwidth,angle=-90]{Fig3d.eps}
\end{center}
\vglue 0.75 truecm
\caption{Many-particle angular-momentum ($\hat L_X$, $\hat L_Y$, and $\hat L_Z$; in red, green, and blue)
variance per particle as a function of the interaction parameter $\Lambda$
computed at the limit of an infinite number of particles
within many-body (dashed lines) and mean-field (solid lines) levels of theory for the ground state
of the three-dimensional anisotropic harmonic-interaction model.
The frequencies of the trap are $\omega_x = 0.7$, $\omega_y = 5.0$, and $\omega_z = 10.5$.
Results at several translations ${\bf r}_0$ of the center of the trap are shown in the panels:
(a) ${\bf r}_0=(0,0,0)$;
(b) ${\bf r}_0=(0.25,0,0)$;
(c) ${\bf r}_0=(0,0.25,0)$;
(d) ${\bf r}_0=(0.25,0.25,0)$.
Different anisotropy classes of the angular-momentum variance emerge with the interaction parameter.
See the text for more details.
The quantities shown are dimensionless.}
\label{f3}
\end{figure}

Still, can the above-found picture of angular-momentum variance anisotropy be made richer?
The answer is positive and requires one to dive deeper into to properties of
angular-momentum variances under translations.
To this end,
we employ and extend prior work in two spatial dimensions \cite{VAR6}.
Suppose now that the harmonic trap is located not in the origin but translated to a general point ${\bf r}_0=(x_0,y_0,z_0)$.
The expectation values per particle of the momentum and angular-momentum operators are still zero,
whereas the expectation value of the position operator is, of course, ${\bf r}_0$.
Whereas the position and momentum variances are invariant to translations,
the angular-momentum variances are not, which open up another degree of freedom to investigate the anisotropy class
of the angular-momentum variances in three spatial dimensions.
Using the transformation properties of the angular-momentum operator $(\hat L_X, \hat L_Y, \hat L_Z)$ under translations,
see appendix \ref{APP} and \cite{VAR6},
the final expressions
for the translated angular-momentum variances per particle
read explicitly
\begin{eqnarray}\label{VAR_Lx_Ly_Lz_MB_TRANS_FINAL}
& &
{\rm Var|_{MB}}(\hat L_X;{\bf r}_0) =
\frac{1}{4} \frac{\left(\sqrt{\omega_y^2+2\Lambda}-\sqrt{\omega_z^2+2\Lambda}\right)^2}{\sqrt{\omega_y^2+2\Lambda} \sqrt{\omega_z^2+2\Lambda}}
\left[1 +
\left(\sqrt{1+\frac{2\Lambda}{\omega_y^2}} - 1\right)\left(\sqrt{1+\frac{2\Lambda}{\omega_z^2}} - 1\right)\right] +
\nonumber \\
& & \quad + y_0^2 \frac{\omega_z}{2} + z_0^2 \frac{\omega_y}{2},
\nonumber \\
& &
{\rm Var|_{MB}}(\hat L_Y;{\bf r}_0) =
\frac{1}{4} \frac{\left(\sqrt{\omega_z^2+2\Lambda}-\sqrt{\omega_x^2+2\Lambda}\right)^2}{\sqrt{\omega_z^2+2\Lambda} \sqrt{\omega_x^2+2\Lambda}}
\left[1 +
\left(\sqrt{1+\frac{2\Lambda}{\omega_z^2}} - 1\right)\left(\sqrt{1+\frac{2\Lambda}{\omega_x^2}} - 1\right)\right] +
\nonumber \\
& & \quad + z_0^2 \frac{\omega_x}{2} + x_0^2 \frac{\omega_z}{2},
\nonumber \\
& &
{\rm Var|_{MB}}(\hat L_Z;{\bf r}_0) =
\frac{1}{4} \frac{\left(\sqrt{\omega_x^2+2\Lambda}-\sqrt{\omega_y^2+2\Lambda}\right)^2}{\sqrt{\omega_x^2+2\Lambda} \sqrt{\omega_y^2+2\Lambda}}
\left[1 +
\left(\sqrt{1+\frac{2\Lambda}{\omega_x^2}} - 1\right)\left(\sqrt{1+\frac{2\Lambda}{\omega_y^2}} - 1\right)\right] +
\nonumber \\
& & \quad + x_0^2 \frac{\omega_y}{2} + y_0^2 \frac{\omega_x}{2}
\end{eqnarray}
at the many-body level of theory and
\begin{eqnarray}\label{VAR_Lx_Ly_Lz_GP_TRANS_FINAL}
& &
{\rm Var|_{GP}}(\hat L_X;{\bf r}_0) = \frac{1}{4} \frac{\left(\sqrt{\omega_y^2+2\Lambda}-\sqrt{\omega_z^2+2\Lambda}\right)^2}{\sqrt{\omega_y^2+2\Lambda} \sqrt{\omega_z^2+2\Lambda}} +
y_0^2 \frac{\sqrt{\omega_z^2+2\Lambda}}{2} + z_0^2 \frac{\sqrt{\omega_y^2+2\Lambda}}{2}, \nonumber \\
& &
{\rm Var|_{GP}}(\hat L_Y;{\bf r}_0) = \frac{1}{4} \frac{\left(\sqrt{\omega_z^2+2\Lambda}-\sqrt{\omega_x^2+2\Lambda}\right)^2}{\sqrt{\omega_z^2+2\Lambda} \sqrt{\omega_x^2+2\Lambda}} +
z_0^2 \frac{\sqrt{\omega_x^2+2\Lambda}}{2} + x_0^2 \frac{\sqrt{\omega_z^2+2\Lambda}}{2}, \nonumber \\
& &
{\rm Var|_{GP}}(\hat L_Z;{\bf r}_0) = \frac{1}{4} \frac{\left(\sqrt{\omega_x^2+2\Lambda}-\sqrt{\omega_y^2+2\Lambda}\right)^2}{\sqrt{\omega_x^2+2\Lambda} \sqrt{\omega_y^2+2\Lambda}} +
x_0^2 \frac{\sqrt{\omega_y^2+2\Lambda}}{2} + y_0^2 \frac{\sqrt{\omega_x^2+2\Lambda}}{2} \
\end{eqnarray}
at the mean-field level of theory.
Let us examine expressions (\ref{VAR_Lx_Ly_Lz_MB_TRANS_FINAL}) and (\ref{VAR_Lx_Ly_Lz_GP_TRANS_FINAL})
more closely.
The terms added to the translated angular-momentum variances at the many-body level
depend on the corresponding components of the translation vector but not on the interaction parameter,
whereas the added terms at the mean-field level of theory depend on
and increase with $\Lambda$, see appendix \ref{APP} for more details.
The combined effects can be seen in Fig.~\ref{f3},
compare panels Fig.~\ref{f3}b,c,d with panel Fig.~\ref{f3}a.

For the translation by ${\bf r}_0=(0.25,0,0)$,
the angular-momentum variances ${\rm Var|_{MB}}(\hat L_X;{\bf r}_0)$ and ${\rm Var|_{GP}}(\hat L_X;{\bf r}_0)$
are invariant quantities,
${\rm Var|_{MB}}(\hat L_Y;{\bf r}_0)$ and ${\rm Var|_{MB}}(\hat L_Z;{\bf r}_0)$
are shifted by interaction-independent values,
and
${\rm Var|_{GP}}(\hat L_Y;{\bf r}_0)$ and ${\rm Var|_{GP}}(\hat L_Z;{\bf r}_0)$
increase by interaction-dependent values,
see (\ref{VAR_Lx_Ly_Lz_MB_TRANS_FINAL}) and (\ref{VAR_Lx_Ly_Lz_GP_TRANS_FINAL}),
and compare Fig.~\ref{f3}b and Fig.~\ref{f3}a.
The combined effect is that now
both
${\rm Var|_{MB}}(\hat L_Y;{\bf r}_0) > {\rm Var|_{MB}}(\hat L_Z;{\bf r}_0) > {\rm Var|_{MB}}(\hat L_X;{\bf r}_0)$
and
${\rm Var|_{GP}}(\hat L_Y;{\bf r}_0) > {\rm Var|_{GP}}(\hat L_Z;{\bf r}_0) > {\rm Var|_{GP}}(\hat L_X;{\bf r}_0)$
hold
for all interaction parameters $\Lambda$ in the range studied.
Consequently,
the anisotropy class of the angular-momentum variance per particle for ${\bf r}_0=(0.25,0,0)$
is $\{1\}$ only.

Next,
for ${\bf r}_0=(0,0.25,0)$,
the angular-momentum variances ${\rm Var|_{MB}}(\hat L_Y;{\bf r}_0)$ and\break\hfill
${\rm Var|_{GP}}(\hat L_Y;{\bf r}_0)$
are invariant quantities,
${\rm Var|_{MB}}(\hat L_X;{\bf r}_0)$ and ${\rm Var|_{MB}}(\hat L_Z;{\bf r}_0)$
are shifted by interaction-independent values,
and
${\rm Var|_{GP}}(\hat L_X;{\bf r}_0)$ and ${\rm Var|_{GP}}(\hat L_Z;{\bf r}_0)$
grow by interaction-dependent values,
contrast Fig.~\ref{f3}c and Fig.~\ref{f3}a.
The combined effect is, of course, different than with ${\bf r}_0=(0.25,0,0)$,
and we discuss its main features,
focusing on the regime of interaction parameters larger than about $\Lambda=6$,
for which the two many-body curves ${\rm Var|_{MB}}(\hat L_X;{\bf r}_0)$ and ${\rm Var|_{MB}}(\hat L_Y;{\bf r}_0)$
cross, see Fig.~\ref{f3}c.
The many-body angular-momentum variances satisfy
${\rm Var|_{MB}}(\hat L_X;{\bf r}_0) > {\rm Var|_{MB}}(\hat L_Y;{\bf r}_0) > {\rm Var|_{MB}}(\hat L_Z;{\bf r}_0)$
for all studied interaction parameters.
So, the effect of this translation is to alter the
order of the many-body variances,
i.e., to change the many-body anisotropy.
Now,
at the mean-field level,
we find
${\rm Var|_{GP}}(\hat L_X;{\bf r}_0) > {\rm Var|_{GP}}(\hat L_Y;{\bf r}_0) > {\rm Var|_{GP}}(\hat L_Z;{\bf r}_0)$
up to about $\Lambda=12$
and
${\rm Var|_{GP}}(\hat L_X;{\bf r}_0) > {\rm Var|_{GP}}(\hat L_Z;{\bf r}_0) > {\rm Var|_{GP}}(\hat L_Y;{\bf r}_0)$
for the interaction parameters larger than about $\Lambda=12$.
All in all,
for the translation by
${\bf r}_0=(0.25,0,0)$ the above-described relations
correspond, respectively, to
the anisotropy classes $\{1\}$ and $\{1,2\}$
of the angular-momentum variances per particle.

Finally,
for the translation by ${\bf r}_0=(0.25,0.25,0)$
non of the angular-momentum variances is invariant,
see (\ref{VAR_Lx_Ly_Lz_MB_TRANS_FINAL}) and (\ref{VAR_Lx_Ly_Lz_GP_TRANS_FINAL}).
The many-body angular-momentum variances
${\rm Var|_{MB}}(\hat L_X;{\bf r}_0)$,
${\rm Var|_{MB}}(\hat L_Y;{\bf r}_0)$,
and
${\rm Var|_{MB}}(\hat L_Z;{\bf r}_0)$
are shifted by interaction-independent values
and
the mean-field angular-momentum quantities
${\rm Var|_{GP}}(\hat L_X;{\bf r}_0)$,
${\rm Var|_{GP}}(\hat L_Y;{\bf r}_0)$,
and
${\rm Var|_{GP}}(\hat L_Z;{\bf r}_0)$
increase by interaction-dependent values,
see Fig.~\ref{f3}d.
We examine the overall effect,
concentrating on the regime of interaction parameters larger than about $\Lambda=2.5$,
where the two many-body curves ${\rm Var|_{MB}}(\hat L_X;{\bf r}_0)$ and ${\rm Var|_{MB}}(\hat L_Z;{\bf r}_0)$
cross, see Fig.~\ref{f3}d.
The many-body angular-momentum variances satisfy
${\rm Var|_{MB}}(\hat L_Y;{\bf r}_0) > {\rm Var|_{MB}}(\hat L_X;{\bf r}_0) > {\rm Var|_{MB}}(\hat L_Z;{\bf r}_0)$
for all interaction parameters in the range studied.
Once again,
the effect of the translation is to change the
order of the many-body variances,
thereby altering the many-body anisotropy,
compare panels Fig.~\ref{f3}a,c,d.
At the mean-field level
one finds
${\rm Var|_{GP}}(\hat L_Y;{\bf r}_0) > {\rm Var|_{GP}}(\hat L_X;{\bf r}_0) > {\rm Var|_{GP}}(\hat L_Z;{\bf r}_0)$
up to about $\Lambda=15$,
then
${\rm Var|_{GP}}(\hat L_Y;{\bf r}_0) > {\rm Var|_{GP}}(\hat L_Z;{\bf r}_0) > {\rm Var|_{GP}}(\hat L_X;{\bf r}_0)$
till about $\Lambda=22.5$,
and
${\rm Var|_{GP}}(\hat L_Z;{\bf r}_0) > {\rm Var|_{GP}}(\hat L_Y;{\bf r}_0) > {\rm Var|_{GP}}(\hat L_X;{\bf r}_0)$
for all the interaction parameters larger than about $\Lambda=22.5$ studied,
see Fig.~\ref{f3}d.
Therefore,
for ${\bf r}_0=(0.25,0.25,0)$ the above-discussed findings
imply that all anisotropy classes can be attained
by the angular-momentum variances per particle
in the three-dimensional harmonic-interaction model at the infinite-particle-number limit.
These are, respectively,
$\{1\}$, $\{1,2\}$, $\{1,2,3\}$,
which rounds off the present work.

\section{Summary}\label{summary}

The present work deals with a connection between anisotropy of and correlations
in a three-dimensional trapped Bose-Einstein condensate.
The merit of treating the limit of an infinite number of bosons is appealing,
since the system is known to be 100\% condensed in this limit,
and some of its properties,
notably the density per particle,
are identical at the many-body and mean-field levels of theory.

We have analyzed the variances per particle of the three Cartesian components of the
position ($\hat X, \hat Y, \hat Z$),
momentum ($\hat P_X, \hat P_Y, \hat P_Z$),
and angular-momentum ($\hat L_X, \hat L_Y, \hat L_Z$) operators
at the many-body and mean-field levels of theory.
In generally, for small interaction parameters
the differences between the many-body and mean-field quantities
are quantitative
whereas for larger interaction parameters qualitative differences emerges.
We define the anisotropy class of the variance according
to the different orderings from large to small,
or permutations,
of the three respective many-body and mean-field quantities.
The anisotropy class $\{1\}$ implies the same ordering,
the class $\{1,2\}$ implies that two of the components are permuted,
and the anisotropy class $\{1,2,3\}$ means that the three components
are permuted.

Two relatively transparent applications are presented,
the first is the breathing of an anisotropic three-dimensional trapped Bose-Einstein condensate,
and the second is the ground state of the anisotropic three-dimensional harmonic-interaction model.
The former exhibits different anisotropy classes of the position and momentum variances per particle,
because at the many-body level
the variances are constant,
while at the mean-field level
the variance of each component oscillates with a different amplitude and frequency,
and consequently different anisotropy classes occur in time.
The latter application shows different anisotropy classes of the angular-momentum variance per particle,
owing to the intricate transformation properties of the angular-momentum variance
when the wave-packet is translated from the origin.

To sum up, the anisotropy or morphology of a three-dimensional trapped Bose-Einstein condensate
can look quite different when examined through the `glasses'
of many-body and mean-field theories,
even for 100\% condensed bosons at the limit of an infinite number of particles.
It would be interesting to conduct the investigation presented here
and analyze the results
in more complicated numerical many-particle scenarios.
Last but not least,
it is possible to envision classifying the morphology of a Bose-Einstein
condensate beyond that emanating from the variances
of the position, momentum, and angular-momentum operators.

\section*{Acknowledgements}

This research was supported by the Israel Science Foundation (Grants No. 600/15 and 1516/19).
We thank Anal Bhowmik and Lorenz Cederbaum for discussions.
Computation time on the BwForCluster,
the High Performance Computing system Hive of the Faculty of Natural Sciences at University of Haifa,
and the Hawk at the High Performance Computing Center
Stuttgart (HLRS) is gratefully acknowledged.