\section{Introduction}

Nuclear matter in the inner crust of neutron stars consists of a lattice of nuclear clusters immersed in superfluid neutron gas\cite{Chamel-Haensel2008,Haensel-book,Pethick-Ravenhall95}.
The coexistence of nuclear clusters and superfluid neutrons brings about rich many-body phenomena.
A typical example might be the pinning and unpinning of superfluid vortices, which is considered to be responsible for the glitch in
the rotation frequency of a neutron star\cite{Anderson-Itoh1975,Alpar1977,Pines-Alpar92}. Another interest is found in collective excitations of inner crust matter.
The lattice vibration, the collective displacement motion of nuclear clusters, may be  microscopic
entity of the seismic oscillation observed as the quasi-periodic oscillation in the giant flares\cite{Duncan1998,Samuelson2007,Andersson09}.
Recently the superfluid phonon\cite{Anderson58,Bogoliubov59,Galitskii58}, called also the Anderson-Bogoliubov
phonon,
attracts new attentions\cite{Aguilera09,Pethick10,Cirigliano11,Page-Reddy2012,Kobyakov13,Kobyakov14,Chamel13,Martin14}.
It may influence the heat capacity (due to this additional degrees of freedom)\cite{Khan05,Martin14},
the seismic oscillation  (through the coupling between the superfluid and the lattice phonons)\cite{Pethick10,Chamel13},
and the thermal conductivity (as a new heat carrier)\cite{Aguilera09,Pethick10,Cirigliano11,Page-Reddy2012}.
Possibility of forming a exotic crystalline structure
is also suggested\cite{Kobyakov13,Kobyakov14}.
The coupling between the nuclear cluster and the superfluid phonon is regarded as a key issue in these studies.

Theoretical approaches to the collective excitation of inner crust matter have been undertaken in
two different directions. One is macroscopic approaches focusing on low-frequency limit\cite{Aguilera09,Pethick10,Cirigliano11,Page-Reddy2012,Kobyakov13,Kobyakov14,Chamel13,Magierski04,Martin16},
where the collective degrees of freedom of the superfluid phonon is explicitly introduced
 in the thermodynamic limit or in the hydrodynamic framework
while the nucleon degrees of freedom are treated implicitly.
Another is microscopic approaches based on the many-body quantum mechanics for interacting
nucleons.  In the latter direction
the linear response methods (the quasiparticle random phase approximation QRPA) based
on nuclear density functional models or the Hartree-Fock-Bogoliubov (HFB) models  are often adopted\cite{Khan05,Grasso08,Baroni2010,Martin14,Inakura17}.
Properties of the superfluid phonon mode are studied in detail within the QRPA
for pure and uniform neutron superfluid\cite{Martin14}. However, if one intends to describe
 the inhomogeneous system consisting of
nuclear cluster and neutron superfluid,  one resorts to QRPA calculations
assuming a spherical Wigner-Seitz cell enclosing a single cluster\cite{Khan05,Grasso08,Baroni2010,Inakura17}.
 In a preceding publication (Part I hereafter)\cite{Inakura17}
we have studied along this line the dipole excitation
in the cell since  the displacement motion of cluster has the multipole $L=1$ and
the dipole channel is responsible for the coupling
between the superfluid phonon and the lattice phonon. This study predicts that
 the superfluid phonon with the dipole mulitpolarity
is influenced strongly by the presence of cluster in such a way that the amplitude of superfluid phonon
is significantly hindered inside the nuclear cluster. This result suggests that the
coupling between the superfluid phonon and the displacement motion of cluster
(i.e. the lattice phonon) might be weak, implying possible enhancement of the
thermal conductivity of the inner crust.

In the present work we analyze the quadrupole excitation in the spherical
Wigner-Seitz cell of inner crust with use of the same model as that in Part I.
In a pioneering work by Khan et al., where a HFB plus QRPA model
is applied to quadrupole excitations, a low-lying
collective mode called the supergiant resonance is predicted\cite{Khan05}. Its relation
to the superfluid phonon is suggested, but only in a qualitative argument.
In another study using a Hartree-Fock plus RPA, in which the nucleon pair correlation is
neglected and the superfluid phonon is not expected,
low-lying collective modes with a character of the multipole shape vibration of cluster are predicted\cite{Baroni2010}.
Because of this situation, a consistent picture is not obtained yet
and questions and problems remain.
What is the relation among
the superfluid phonon, the supergiant resonance in Khan et al.\cite{Khan05} and
the shape vibration of cluster?
How is the superfluid phonon affected by the presence of cluster? Is there difference between
the dipole and quadrupole modes of the superfluid phonon?
Does the cluster shape vibration exist or how is it affected when it is immersed in  superfluid neutron gas?
The purposes of the present study is to answer these
 questions.
In short, we intend to  clarify the nature of the low-lying excitation modes of superfluid inner crust matter
with focus on quadrupole excitations around a spherical nuclear cluster.

In the following we call the superfluid phonon mode of neutron
superfluid  the Anderson-Bogoliubov phonon mode  (abbreviated as AB phonon) to keep the same nomenclature as
in Part I\cite{Inakura17}.
The paper is organized as follows:
In Sec. II, we briefly explain the model we use to describe the ground state and the excitation modes of inner crust matter in the
Wigner-Seitz approximation. In Sec. III,
we show numerical results of the QRPA calculation performed for the quadrupole excitations in the spherical Wigner-Seitz
cell. We will show that there emerge two low-lying collective excitations, which are related to
both the AB phonon mode of superfluid neutron gas
and the quadrupole shape vibration of the nuclear cluster.  We will investigate in detail basic properties of these
low-lying collective modes to clarify the interplay between the AB phonon mode and the shape vibration of cluster.
It will be turned out that the quadrupole shape collectivity of the cluster exhibits a peculiar shell effect which arises
not from neutrons bound to clusters, but from unbound neutrons permeating the nuclear cluster and neutron gas.
Sec. IV is devoted to the conclusion.

\section{Model}

We employ the HFB model and the QRPA theory.
The calculation framework and the model parameters are
the same as those in Part I\cite{Inakura17}.
Here we recapitulate the adopted HFB+QRPA formalism briefly.
For other details, we refer the readers to Refs.~\cite{Inakura17,Matsuo01,Serizawa09,Matsuo10,Nakatsukasa2016}.

To calculate the properties of the inner crust, we introduce the spherical Wigner-Seitz approximation
with a Wigner-Seitz radius of $R_\mathrm{WS}=20$~fm.
Following the standard prescriptions adopted in the HF/HFB
calculations\cite{Negele73,Barranco1998,Sandulescu04a,Sandulescu04b,Baldo2005,Baldo2006,Grill11,Pastore11,Pastore12},
we impose the von-Neumann-Dirichlet boundary condition\cite{Negele73}.
Nucleon configuration of the Wigner-Seitz cell is specified with the neutron chemical potential (Fermi energy)
$\lambda_n$  and a fixed integer number $Z$ of protons in the cell.

The ground state configuration
 of inner crust matter for given $Z$ and $\lambda_n$ is obtained
by solving the HFB equation in the coordinate representation.
We adopt the zero-range effective Skyrme SLy4 interaction~\cite{SLy4} for the Hartree-Fock
potential (in the particle-hole channel).
A version of density-dependent delta interaction~\cite{Matsuo06,Matsuo07} is used to obtain the
 pairing potential (in the pairing channel)
\begin{eqnarray}
\Delta_\tau(\mbox{\boldmath$r$}) = V_\mathrm{pair} \left[ 1 - \eta \left(\frac{\rho_\tau(\mbox{\boldmath$r$})}{\rho_c}\right)^\alpha \right] \tilde{\rho}_\tau(\mbox{\boldmath$r$})
\label{pairing.Delta}
\end{eqnarray}
and the pair density
\begin{eqnarray}
\tilde{\rho}_\tau(\mbox{\boldmath$r$})=\langle \Psi_0 | \psi     (\mbox{\boldmath$r$}\uparrow)   \psi     (\mbox{\boldmath$r$}\downarrow) | \Psi_0 \rangle.
\end{eqnarray}
Here $\tau=n,\, p$ denotes isospin, $\uparrow$ and $\downarrow$ are up and down components of spin $\sigma$, and $ \rho_c=0.08$ fm$^{-3}$.
We set the cut-off $j_\mathrm{max}= (75/2)\hbar$ for the single-particle partial waves and we include all the quasiparticle states under the cut-off quasiparticle energy $E_\mathrm{cut}=$ 60 MeV. To characterize the neutron superfluid
outside the nuclear cluster
we evaluate its average neutron density and average pair potential
 as
\begin{eqnarray}
\rho_\mathrm{ext}   = \frac{\int^{R_2}_{R_1} dr \, r^2 \rho_n(r)  }{\int^{R_2}_{R_1} dr \, r^2} \,,\quad
\Delta_\mathrm{ext} = \frac{\int^{R_2}_{R_1} dr \, r^2 \Delta_n(r)}{\int^{R_2}_{R_1} dr \, r^2}
\label{matter.Delta.rho}
\end{eqnarray}
with $R_1 =$ 12 fm and $R_2 =$ 18 fm.

The QRPA calculation is performed on top of the HFB ground state $\Psi_0$.
In the present study we focus on the quadrupole excitations.
We utilize the linear response formalism~\cite{Matsuo01,Serizawa09,Matsuo10} for the QRPA.
We describe responses of the system with respect to the quadrupole operator $Q$,
and the pair addition and removal operators, $P_\mathrm{add}$ and $P_\mathrm{rm}$, defined by
\begin{eqnarray}
&& Q              = \sum_{\sigma} \int \!d\mbox{\boldmath$r$}\, r^2Y_{2M}(\hat{\mbox{\boldmath$r$}})\, \psi^\dag(\mbox{\boldmath$r$}\sigma)\psi(\mbox{\boldmath$r$}\sigma)  \,, \nonumber \\
&& P_\mathrm{add} =         \int \!d\mbox{\boldmath$r$}\,  Y_{2M}(\hat{\mbox{\boldmath$r$}})\, \psi^\dag(\mbox{\boldmath$r$}\downarrow)\psi^\dag(\mbox{\boldmath$r$}\uparrow) \,, \nonumber \\
&& P_\mathrm{rm}  =         \int \!d\mbox{\boldmath$r$}\,  Y_{2M}(\hat{\mbox{\boldmath$r$}})\, \psi     (\mbox{\boldmath$r$}\uparrow)  \psi     (\mbox{\boldmath$r$}\downarrow) \,.
\label{eq:operators}
\end{eqnarray}
We solve the QRPA linear response equations for fluctuations of the nucleon density $\rho(\mbox{\boldmath$r$})$, the
nucleon pair density $\tilde{\rho}(\mbox{\boldmath$r$})$  and its complex conjugate $\tilde{\rho}^*(\mbox{\boldmath$r$})$.
The spectral representation is adopted for the density response function and
all the quasiparticle states adopted in the HFB calculation are included.
We calculate the strength function
\begin{eqnarray}
S(O;E) = \sum_{Mi} \delta(E - E_i) \left|\langle \Psi_i^{2M}| \hat{O} | \Psi_0\rangle\right|^2 \,,
\end{eqnarray}
for the quadrupole excited states $ \Psi_i^{2M} $ with respect to
the operators $\hat{O}=Q, P_\mathrm{add}$ and $P_\mathrm{rm}$
of neutrons and protons separately.
With a small imaginary
constant $\epsilon$ in the energy argument, the
delta function peaks in
the strength functions are smeared with the Lorentzian function having
 the FWHM of $2\epsilon$.
We evaluate the strength
$B(O)=\sum_M\left|\langle \Psi_i^{2M}| \hat{O} | \Psi_0\rangle\right|^2$
of each excited state by integrating the strength function in an energy interval $E \in [E_i - 10 \epsilon, E_i + 10 \epsilon]$ around its peak energy $E_i$. We employ $\epsilon=10$ keV.
Three transition densities from the HFB ground state $\Psi_0$ to the $i$th QRPA excited state $\Psi_i^{2M}$,
\begin{eqnarray}
&& \delta       \rho _\mathrm{ph}(\mbox{\boldmath$r$}) =   \langle \Psi_0 | \sum_\sigma\psi^\dag(\mbox{\boldmath$r$}\sigma)     \psi     (\mbox{\boldmath$r$}\sigma)     | \Psi_i^{2M} \rangle = Y_{2M}(\hat{\mbox{\boldmath$r$}})\delta \rho _\mathrm{ph}(r)  \,,\nonumber \\
&& \delta\tilde{\rho}_\mathrm{pp}(\mbox{\boldmath$r$}) = \langle \Psi_0 | \psi     (\mbox{\boldmath$r$}\uparrow)   \psi     (\mbox{\boldmath$r$}\downarrow) |
\Psi_i^{2M}  \rangle = Y_{2M}(\hat{\mbox{\boldmath$r$}})\delta\tilde{\rho}_\mathrm{pp}(r) \,,\nonumber \\
&& \delta\tilde{\rho}_\mathrm{hh}(\mbox{\boldmath$r$}) =  \langle \Psi_0 | \psi^\dag(\mbox{\boldmath$r$}\downarrow) \psi^\dag(\mbox{\boldmath$r$}\uparrow)   |
\Psi_i^{2M}  \rangle  = Y_{2M}(\hat{\mbox{\boldmath$r$}})\delta\tilde{\rho}_\mathrm{hh}(r)\,,
\label{eq:densities}
\end{eqnarray}
are obtained from the corresponding fluctuating densities at the peak energy $E_i$.
Note that the calculated spectrum is discretized because of the boundary condition.

As the residual interaction which enters in the QRPA calculation,
we adopt the same effective pairing interaction for the pairing channel.
Concerning the residual interaction in the particle-hole channel,
we adopt the Landau-Migdal approximation~\cite{Serizawa09,Matsuo10,Khan04,Khan09,Paar07} with a renormalization
scheme often employed in this approximation. Namely we replace the
self-consistent particle-hole interaction $v_{\mathrm{ph}}$ by the
Landau-Migdal interaction $f \times v_{\mathrm{LM}}$ derived from the Skyrme
interaction and renormalized with a factor $f$.
The factor $f$ is determined to bring a peak corresponding to the displacement motion
to zero energy in the drip-line nucleus for a given $Z$.
The same value of $f$ is used for the Wigner-Seitz cells with the same $Z$.

\section{Results and discussion}

We have performed the HFB and QRPA calculations systematically for various configurations
specified by the neutron chemical potential $\lambda_n$ and
the proton number $Z$, with changing $Z$ and $\lambda_n$, in order to study basic properties of excitation modes of
 inner crust matter.
The adopted proton numbers are $Z=$ 20, 28, 40, and 50, chosen to cover the range predicted
in the previous HFB or HF calculations for the equilibrium \cite{Negele73,Baldo2005,Baldo2006,Grill11}.
For the neutron chemical potential, we vary it in the range
$\lambda_n = 1-6$ MeV, which corresponds to the density of neutron matter
$\rho_n \approx 5 \times 10^{-4} - 1 \times 10^{-2}$ fm$^{-3}$.
The ground states obtained by the HFB calculation are same as in Part I~\cite{Inakura17}.

\subsection{Low-lying quadrupole excitations}

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.490\textwidth,keepaspectratio]{rpa.NM.Z28.DDDIbare.E0400.neutron.proton.Padd.Premove.upto30MeV.fixed_f.eps}
\caption{(a) Strength function $S(Q_n;E)$ for neutron quadrupole operator,
(b) $S(Q_p;E)$ for proton quadrupole operator, (c) $S(P_\mathrm{add};E)$ for neutron pair-addition operator
and (d) $S(P_\mathrm{rm};E)$ for neutron pair-removal operator, calculated
for $Z=28$ system with $\lambda_n =$ 4.0 MeV,
plotted as a function of the excitation energy $E$.
Arrows indicate the threshold energy $2\Delta_\mathrm{ext}$.
The smearing parameter is $\epsilon=$ 10 keV. Dashed curve
in (b) is the proton quadrupole strength function  $S(Q_p;E)$ obtained with
$\epsilon=$ 500 keV (with the scale on the right axis).
}
\label{neutron.proton.Padd.Prm}
\end{center}
\end{figure}

We first discuss the case of $Z=28$ and $\lambda_n =$ 4.0 MeV  chosen as representative.

Figures~\ref{neutron.proton.Padd.Prm}(a)-(d) show the strength functions
$S(O;E)$ for the neutron quadrupole moment, $O=Q_n$, the proton quadrupole moment  $Q_p$,
the neutron pair-addition operator $P_\mathrm{add}$ and the neutron pair-removal operator $P_\mathrm{rm}$,
respectively,
for the $Z=28$ system with $\lambda_n = 4.0$ MeV.  Two distinct peaks at $E=$ 1.57 MeV and
at $E=$ 2.06 MeV are seen in the low-lying region.
These peaks are located below $2\Delta_\mathrm{ext}$ (twice of the neutron pair gap, indicated with arrows).
Above this threshold energy, the QRPA excited states
are densely distributed (and would become continuous if $R_\mathrm{WS} \rightarrow \infty$).
It is noted here that the peak at $E=$ 1.57 MeV
is significant in the strength functions
for the neutron operators $Q_n$, $P_\mathrm{add}$  and $P_\mathrm{rm}$, but not in the
strength function for the proton operator $Q_p$. The second peak at $E=$ 2.06 MeV, however,
has the largest strength in the proton quadrupole moment $Q_p$, but the strengths for
the neutron operators are much smaller than those of the first peak at $E=$ 1.57 MeV.
Apparently there exist two kinds of low-lying collective
modes which have
different characters.

Let us examine the peak at $E=$ 1.57 MeV in detail.
This state has the strengths $B(Q_n)= 1.49 \times 10^6$ fm$^4$ ,
$B(P_\mathrm{add})= 183$ and $B(P_\mathrm{rm})= 63$ for the neutron operators, which
are the largest among the QRPA excited states and exhaust
67 \%, 21 \%, and 42 \%, respectively,
of the total strengths (non-energy-weighted sums integrated up to 50 MeV).
Note that the neutron quadrupole strength $B(Q_n)=1.49 \times 10^6$ fm$^4$ is larger by a factor of  $10^{2-3}$
than a typical value of the low-lying collective $2_1^+$ state in isolated nuclei,
$B(E2)/e^2 \sim 10 B_\mathrm{W}$ ($B_\mathrm{W}$ being the Weisskopf unit\cite{Ring-Schuck}).
Having this extraordinary quadrupole strength, we consider that this state may correspond to what is called
the supergiant quadrupole resonance in Khan et al.\cite{Khan05}.

Figure~\ref{drho}(a) shows the neutron and proton particle-hole transition densities,
$\delta\rho^\nu_\mathrm{ph}(r)$ and $\delta\rho^\pi_\mathrm{ph}(r)$, and the neutron pair-addition and -removal
transition densities, $\delta\tilde{\rho}^\nu_\mathrm{pp}(r)$ and $\delta\tilde{\rho}^\nu_\mathrm{hh}(r)$,
of the 1.57 MeV state. It is seen that the neutron particle-hole transition density
$\delta\rho^\nu_\mathrm{ph}(r)$ as well as the neutron pair-addition/removal
transition densities, $\delta\tilde{\rho}^\nu_\mathrm{pp}(r)$ and $\delta\tilde{\rho}^\nu_\mathrm{hh}(r)$,
have finite amplitudes in the whole region of neutron gas $r \m@thcombine>\sim 6$ fm. This is the origin of the extremely large
neutron strengths. The three transition densities far outside the nucleus $r \m@thcombine>\sim 9$ fm exhibits
a signature of the AB phonon mode: They share an approximate
sinusoidal shape in the neutron gas region while the pair transition densities $\delta\tilde{\rho}^\nu_\mathrm{pp}(r)$ and $\delta\tilde{\rho}^\nu_\mathrm{hh}(r)$ have an opposite phase to one another\cite{Inakura17}.
The neutron and proton particle-hole
transition densities  $\delta\rho^\nu_\mathrm{ph}(r)$ and $\delta\rho^\pi_\mathrm{ph}(r)$, however,
show enhanced amplitude in the surface region of the cluster, 4 fm  $\m@thcombine<\sim r \m@thcombine<\sim$ 8 fm,
with the phase coherence. This behaviour in the surface region is similar to the one known in the
low-lying quadrupole shape oscillation of isolated nuclei.

The 2.06 MeV state is recognized in the proton quadrupole
strength function, but it is hard to be seen in the neutron strength functions.
The proton strength $B(Q_p)=$ 508 fm$^4$ of the 2.06 MeV state is the same order of magnitude as
the typical $B(E2)$ values of the low-lying quadrupole vibration in isolated nuclei.
(Note $B(E2)=e^2B(Q_p)$.) The neutron quadrupole strength $B(Q_n)=3.46 \times 10^5$ fm$^4$
is two orders of magnitude larger than those in isolated nuclei, but significantly smaller than
that of the 1.57 MeV state.
The transition densities of the 2.06 MeV state, shown in Fig.~\ref{drho}(b), displays a similar
behavior as those of the 1.57 MeV state, but some difference is clearly seen: The neutron
and proton particle-hole amplitudes around the
cluster surface is much larger than those of the 1.57 MeV state and it suggests a closer relation to the
quadrupole shape vibration of the cluster. The neutron amplitudes outside the cluster
displays a node at $r \approx$ 12fm.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.4850\textwidth,keepaspectratio]{transition.density.NM.Z28.E0400.E2.fixed_f.2.eps}
\caption{Transition densities of the low-lying quadrupole collective modes whose excitation energies are
(a) $E=$ 1.57 MeV, (b) 2.06 MeV, and (c) 3.03 MeV, obtained for $Z=$ 28 system with $\lambda_n=$ 4.0 MeV.
Red, gray, green, and blue curves are $\delta\rho^\nu_\mathrm{ph}(r)$, $\delta\rho^\pi_\mathrm{ph}(r)$,
$\delta\tilde{\rho}^\nu_\mathrm{pp}(r)$, and $\delta\tilde{\rho}^\nu_\mathrm{hh}(r)$, respectively.
The arrows indicate the surface position of the cluster\cite{Inakura17}. See text for details.
}
\label{drho}
\end{center}
\end{figure}

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.4730\textwidth,keepaspectratio]{rpa.drho.NM.Z28.DDDIbare.E0400.E2.Padd.residual.fixed_f.2.eps}
\caption{
Uppers: (a) Strength function $S(Q_n;E)$, (b) $S(Q_p;E)$, and (c) $S(P_\mathrm{add};E)$
calculated with the full residual interaction (red curve), with only the residual pairing interaction (green dashed), and
with only the residual particle-hole interpretation (blue long-dashed).
Unperturbed results (magenta long-dashed dotted) are also plotted.
The system is  $Z=28$ and $\lambda_n = 4.0$ MeV.
Lowers: Transition densities of the low-lying collective mode obtained in calculations
(d) with only residual pairing interaction and (e) with only the residual particle-hole interaction.
Red, gray, green, and blue curves are  $\delta\rho^\nu_\mathrm{ph(r)}$, $\delta\rho^\pi_\mathrm{ph}(r)$,
$\delta\tilde{\rho}^\nu_\mathrm{pp}(r)$, and $\delta\tilde{\rho}^\nu_\mathrm{hh}(r)$, respectively.
See also the captions of Figs.~\ref{neutron.proton.Padd.Prm} and \ref{drho}.
}
\label{residual}
\end{center}
\end{figure}

\subsection{Coexistence of AB phonon  and quadrupole cluster vibration}

In order to investigate the characteristic features of these two collective excitations,
we look into roles of the residual interaction.
For this purpose we shall drop off some parts of the residual interaction which enters in the QRPA calculation.
Figure~\ref{residual} shows the strength functions $S(Q_n;E)$, $S(Q_p;E)$ and $S(P_\mathrm{add};E)$,
1) which are obtained using
only the residual pairing interaction (i.e. the particle-hole part of the residual interaction is neglected)
plotted with green curves, 2) those (blue curves) obtained with only the residual particle-hole interaction
(the residual pairing interaction is neglected instead), and 3) unperturbed results (magenta curves) in which
all the
residual interactions are neglected.
The obtained strengths functions are shown in the panels (a), (b) and (c).
The calculation using only the residual pairing interaction produces
a collective state which appears as a peak at 2.18 MeV
with large neutron strengths. This state has apparent characters of the
AB phonon mode as are seen
in its transition densities, shown in Fig.~\ref{residual}(d):
the oscillatory
behavior with long wave length, $\delta\tilde{\rho}^\nu_\mathrm{pp}(r) \sim  A j_2(qr) $ and
$\delta\tilde{\rho}^\nu_\mathrm{hh}(r) \sim - A j_2(qr) $,
is clearly seen, especially outside the nuclear cluster $r \m@thcombine>\sim 8$ fm.
When only the particle-hole interaction included (and the residual pairing interaction is neglected),
the quadrupole shape vibration of the cluster emerges seen as
a peak at 2.23 MeV in the proton strength function. In fact,
transition densities of this state [Fig.~\ref{residual}(e)] exhibit clearly features of
 isoscalar quadrupole shape vibration of the cluster:
The particle-hole transition densities $\delta\rho^\nu_\mathrm{ph}(r)$ and $\delta\rho^\pi_\mathrm{ph}(r)$
of neutrons and protons both have large amplitudes with the same phase
around the surface region of the cluster, 4 fm $\m@thcombine<\sim r \m@thcombine<\sim$ 8 fm,
while all the transition densities diminish outside the cluster $r \m@thcombine>\sim 10$ fm.

Having these observations, we can interpret the two excitation modes
at $E=$1.57 MeV and 2.06 MeV in the full QRPA calculation
 as mixtures or coupled modes of the AB phonon mode and the quadrupole shape vibration of the cluster. The transition
densities of the 1.57 MeV state, Fig.~\ref{drho}(a),  appear  an in-phase superposition of those
in Figs.~\ref{residual}(d) and (e) while the component of the AB phonon mode is dominant.
The dominance of the AB phonon is seen also in the neutron strengths [Fig.~\ref{residual}(a)(c)],
where the peak heights (red curves) are
comparable to those of the pure AB phonon mode (green curves) obtained with only the pairing interaction.
Therefore the 1.57 MeV state may be characterized as a AB-phonon-dominant mode.
Concerning the 2.06 MeV state, on the other hand,
the transition densities in Fig.~\ref{drho}(b) appear an approximate out-of-phase superposition of
those in Fig.~\ref{residual}(d) and (e) while the component of the cluster shape vibration is dominant.
Therefore it is characterized as a cluster-vibration-dominant mode.
It is not surprising that the AB phonon mode of superfluid neutron gas and the shape vibration of the cluster
couple with one another as the neutron gas and the nuclear cluster share the
same nucleon degrees of freedom. In addition the low-lying shape vibration of nuclei is known to be
influenced by the pair correlation, and hence the same is expected for the cluster shape vibration.
 A remarkable point here is that the
characters of the AB phonon mode and the cluster shape vibration remain well  in spite of the mixing.
This suggests that the AB phonon mode and the cluster shape vibration coexist and
the coupling between the two fundamental excitation modes is rather weak
in the case of the quadrupole excitations. This is consistent with our finding in the dipole case\cite{Inakura17},
where weak coupling between the dipole AB phonon mode and the displacement motion of cluster is suggested.

We briefly mention other types of collective excitation modes which are identified in
the present QRPA calculation.  The third largest peak in the neutron quadrupole strength function
[Fig.~\ref{neutron.proton.Padd.Prm} (a)] has an excitation energy of 3.03 MeV.
The transition densities of this state, shown in Fig.~\ref{drho}(c),
suggest that this state may be interpreted as the AB phonon mode of the second harmonics with one node
at $r \approx 12$ fm with small mixing of the cluster shape vibration. A bunch of states around
$E \approx$ 11 to 16 MeV having significant proton quadrupole strength is a high-lying
collective excitation of the cluster, which corresponds to the isoscalar giant quadrupole resonance in isolated nuclei.
A distribution of the proton strengths around 25 MeV is an analog of
the isovector giant quadrupole resonance.
Detailed analysis of these collective excitations are out of scope of the present paper, and is left for future study.

We remark here relation between the low-lying collective excitation modes found in the present calculation
and those discussed in the preceding works\cite{Khan05,Baroni2010}.
Khan et al.\cite{Khan05}  adopted the HFB plus QRPA model, which is similar to the present one.
They discuss collective excitations in terms of the strength function for the isoscalar
quadrupole operator $Q_{IS}=Q_n + Q_p$, and found a low-lying significant peak, named
the supergiant resonance. We find in the present calculation that the isoscalar quadrupole
strength function $S(Q_{IS};E)$ looks similar to the neutron strength function $S(Q_{n};E)$
when plotted in the same scale. (This is because the proton strength is smaller by a factor of $10^{-3}$
than the neutron strength, see Figs.~\ref{drho}(a) and (b).) Therefore the supergiant resonance of
Khan et al.\cite{Khan05} may correspond to the AB-phonon-dominant mode in the present calculation.
It is likely that the cluster-vibration-dominant mode is not visible in the analysis of Khan et al. due to
their use of relatively large smoothing constant (150keV), which may hide a small peak nearby the large peak.
Baroni et al. performed the HF plus RPA calculation, which neglects the nucleon pair correlation\cite{Baroni2010}.
The HF plus RPA model predicts low-lying quadrupole shape vibration mode as well as the high-lying mode
which corresponds to isoscalar
giant quadrupole resonance of the cluster. The low-lying shape vibration mode in Ref.~\cite{Baroni2010}, however,
exhibits some differences from the cluster-vibration-dominant mode in the present model: The former
emerges as a bump of strength distribution while the latter is a single discrete state. These differences are attributed to the
neutron pair correlation and the presence of the pair gap. Note also that due to the lack of the
pair correlation the AB phonon mode is not
expected to exist in Baroni et al. as we have seen in
the reduced QRPA calculation including only the residual particle-hole interaction (blue curve) shown in Fig.~\ref{residual}.

\subsection{Systematics}

\begin{figure}[tbh]
\begin{center}
\includegraphics[width=0.4850\textwidth,keepaspectratio]{rpa.NM.Z28.DDDIbare.E.neutron.proton.Padd.E2.upto6MeV.fixed_f.eps}
\caption{
Strength functions (a) $S(Q_n;E)$, (b) $S(Q_p;E)$, and (c) $S(P_\mathrm{add};E)$
for $Z=28$ system with $\lambda_n =$ 1.0 MeV (red curves),
2.0 MeV (brown), 3.0 MeV (yellow green), 4.0 MeV (blue green), 5.0 MeV (blue), and 6.0 MeV (purple).
Colored arrows indicate the threshold energy $2\Delta_\mathrm{ext}$ for corresponding $\lambda_n$.
See text for details.
}
\label{ABmode.quadrupole.lambda}
\end{center}
\end{figure}

We shall discuss systematic behaviour of the two kinds of low-lying collective modes by varying $\lambda_n$ and $Z$.

\begin{figure}[tbh]
\begin{center}
\includegraphics[width=0.4850\textwidth,keepaspectratio]{transition.density.NM.Z28.E.E2.ABdominant.fixed_f.eps}
\caption{
Transition densities
$\delta\rho^\nu_\mathrm{ph}(r)$, $\delta\rho^\pi_\mathrm{ph}(r)$,
$\delta\tilde{\rho}^\nu_\mathrm{pp}(r)$, and $\delta\tilde{\rho}^\nu_\mathrm{hh}(r)$ of the
AB-phonon-dominant-mode for $Z=28$ systems with  $\lambda_n=$
(a) 1.0 MeV, (b) 2.0 MeV, (c) 3.0 MeV, (d) 5.0 MeV and (e) 6.0 MeV.
See text for details.
}
\label{ABdominant}
\end{center}
\end{figure}

We first examine dependence on the neutron chemical potential $\lambda_n$.
In Fig.~\ref{ABmode.quadrupole.lambda} we plot the strength functions
$S(Q_n;E)$, $S(Q_p;E)$, and $S(P_\mathrm{add};E)$ for $Z=28$ systems with $\lambda_n=$ 1.0, 2.0, $\cdots$, 6.0 MeV.
It is evident in the neutron strength functions, Figs.~\ref{ABmode.quadrupole.lambda}(a)(c), that
a low-lying collective state having large neutron strengths
appears around 1 -- 2 MeV of excitation energy in all the cases
of $\lambda_n$.
The transition densities of these states, shown in Fig.~\ref{ABdominant}, indicate that they have the
character of the AB phonon mode in the neutron gas region $r \m@thcombine>\sim $ 8 fm. Although admixture of
the cluster shape vibration is seen in the transition densities,
the proton quadrupole strength of these states is not significant (see
Fig.~\ref{ABmode.quadrupole.lambda}(b)).  Therefore all these states are interpreted as the
AB-phonon-dominant mode.  The excitation energy and the neutron strengths increase monotonically with
increase of the neutron chemical potential $\lambda_n$ .

The proton quadrupole strength, shown in Fig.~\ref{ABmode.quadrupole.lambda}(b),
exhibits different behaviours. The low-lying collective state with enhanced proton strength are found
only for $\lambda_n=$ 3.0 MeV and 4.0 MeV.
These are the cluster-vibration-dominant mode discussed in the subsection just above.
There exists no cluster-vibration-dominant mode in the low-lying region
in the cases of $\lambda_n=$ 2.0 MeV, 5.0 MeV and 6.0 MeV.
(In case of  $\lambda_n=$ 1.0 MeV, we find a state which has partial character of the cluster vibration
at $E=$ 1.18 MeV, which lies above the continuum threshold $2 \Delta_\mathrm{ext}$.)

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.49990\textwidth,keepaspectratio]{quadrupole.Z28.lambda.Esurf.BE2.eps}
\caption{Excitation energy $E_{AB}$ of the lowest AB-phonon-dominant mode (open circle) and
the energy $E_{cl}$ of the cluster-vibration-dominant mode (filled diamond) for $Z=$ 28 systems,
plotted as a function of the neutron chemical potential $\lambda_n$.
Blue dashed curve is the threshold energy $2\Delta_\mathrm{ext}$.
Red dots connected with line denote the proton quadrupole strengths $B(Q_p)$ of the cluster-vibration-dominant mode.
The results for near-drip-line nuclei $^{80,82,\cdots,88}$Ni are also plotted.
The thin black line is the hydrodynamic estimate $\hbar\omega_{\rm hyd}$ of the quadrupole
AB phonon mode in uniform superfluid neutron gas in the same Wigner-Seitz cell.
}
\label{ABmode.Z28}
\end{center}
\end{figure}

Figure~\ref{ABmode.Z28} shows more detailed information on the excitation energy $E_\mathrm{AB}$ of
the AB-phonon-dominant mode
and the energy $E_\mathrm{cl}$ of
the cluster-vibration-dominant mode, plotted as a function of $\lambda_n$.
The proton quadrupole strength $B(Q_p)$ of the cluster-vibration-dominant mode is also plotted.
For comparison we also plot the excitation energy and $B(Q_p)$ of the $2_1^+$ state of the near-drip-line nuclei $^{80-88}$Ni
which have negative $\lambda_n$.

The excitation energy $E_\mathrm{AB}$ of the AB-phonon-dominant mode has smooth dependence on $\lambda_n$.
It agrees well with an  hydrodynamic estimate of the excitation energy of the quadrupole AB mode, $\hbar\omega_\mathrm{hyd}$,
which is obtained assuming  that the phonon amplitude is expressed by the spherical Bessel function, $\phi(r) \propto j_2(qr)$
with boundary condition $d\phi(r)/dr =0$ at $r=R_\mathrm{WS}=20$ fm.\footnote{
In this estimate, we assume $\phi(r) \propto j_2(qr)$ together with the boundary condition $\left. d\phi(r)/dr \right|_{r=R_\mathrm{box}}=0$.
Combination of the phonon dispersion relation $\omega = c_\mathrm{ph} q$ and
the hydrodynamic estimate of the phonon velocity $c_\mathrm{ph}=\sqrt{(\partial\lambda_n/\partial\rho)\rho/m}$
gives an estimated excitation energy $\hbar\omega_\mathrm{hyd}$.
Here, we use $\rho(\lambda_n)$ evaluated for uniform neutron matter with the same Skyrme functional SLy4.
}
This indicates that the presence of the nuclear cluster
affects only to a small extent the excitation
energy of the AB-phonon-dominant mode. It is also consistent with our observation
in Figs.~\ref{neutron.proton.Padd.Prm} and \ref{residual} that the coupling between the quadrupole
shape vibration of cluster and the quadrupole AB phonon is weak.
For $0< \lambda_n <$ 1 MeV, we do not plot the results since
the collective modes do not appear at low excitation energies
below $2\Delta_\mathrm{ext}$.
Note however that the AB phonon mode would appear below
$2\Delta_\mathrm{ext}$ if we adopted a large Wigner-Seitz cell.

The cluster-vibration-dominant mode displays strong, non-monotonic and oscillatory $\lambda_n$-dependence.
As $\lambda_n$ increases slightly above 1.3 MeV the excitation energy $E_\mathrm{cl}$
drops off sharply. The $B(Q_p)$ value increases dramatically. This behaviour reminds us
 of the softening of low-lying quadrupole
vibration in isolated nuclei, and the same appears to happen for the nuclear cluster in the
Wigner-Seitz cell.
The cluster-vibration-dominant mode is not found
in the interval of $\lambda_n \approx 1.6 - 2.0$ MeV, and it suggests instability of
the spherical nuclear cluster
 toward the quadrupole deformation. (In this case
 the QRPA energy of the quadrupole shape vibration is expected to be  ``imaginary".)
In the interval of $\lambda_n \approx  2.1 \sim 4.8$ MeV, the cluster-vibration-dominant mode emerges at finite
excitation energy.
The softening is seen above $\lambda_n \m@thcombine>\sim 4.5$ MeV and the instability occurs again in the interval
$\lambda_n \approx 4.8 \sim 6$ MeV.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.4380\textwidth,keepaspectratio]{quadrupole.Z20.Z40.Z50.lambda.Esurf.BE2.eps}
\caption{The same as Fig.~\ref{ABmode.Z28}, but for systems with $Z=$ 20, 40, and 50.
}
\label{ABmode.Zall}
\end{center}
\end{figure}

Figure~\ref{ABmode.Zall} shows the same plot as Fig.~\ref{ABmode.Z28}, but for systems
with $Z=$ 20, 40 and 50. The energy of the AB-phonon-dominant mode agrees well with the
hydrodynamical estimate while the strong and oscillatory $\lambda_n$-dependence of the
cluster-vibration-dominant mode is present in all the cases.

\subsection{Resonance shell effect on cluster shape vibration}

\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.480\textwidth,keepaspectratio]{spe.Z28.neutron.proton.paper.eps}
\caption{(a) Neutron  single-particle energies for $Z=$ 28 systems with $\lambda_n=$ 1.0, 2.0, $\cdots$, 6.0 MeV.
Thick back bars connected with dashed lines denote the resonant single-particle states with quantum
numbers  $1h_{11/2}$, $1h_{9/2}$,  $1i_{13/2}$, and $1i_{11/2}$. Bound high-$j$ orbits are also connected with
thin dashed lines for reference.
Red horizontal bars connected with dashed line denote the neutron Fermi energy (chemical potential) $\lambda_n$, and blue
bars the neutron Hartree-Fock potential  $U_\mathrm{ext}$ in the region of neutron gas.
 (b) The same as (a), but for protons.
Numbers 28, 50 and 82 inserted in levels are the spherical magic numbers.
}
\label{spe}
\end{center}
\end{figure}

The strong and oscillatory $\lambda_n$-dependence suggests possible shell effect in the cluster-vibration-dominant mode.
However, the conventional idea of the shells characterized by the magic numbers $N=$ 50, 82, 126, etc. is not
appropriate here as  neutron single-particle orbits around the
neutron Fermi energy $\lambda_n$ are unbound states having a quasi-continuum spectrum.

In order to clarify the suggested shell effect, we look into the single-particle eigenstates in the Hartree-Fock potential.
Figure~\ref{spe} shows the neutron and proton single-particle energies for the $Z=$ 28 systems with
$\lambda_n= 1.0, 2.0, \cdots, 6.0$ MeV. The thick blue bar indicates the depth $U_{\rm ext}$ of the Hartree-Fock
potential in the region of neutron gas, namely in the external region outside the nuclear cluster. The single-particle
states with energy below $U_{\rm ext}$ are localized states bound to the nuclear cluster while those above
$U_{\rm ext}$ are unbound states, which also have discrete energies because of the Wigner-Seitz approximation with
finite cell radius.
Concerning the bound single-particle states, the single-particle energy gaps at
the neutron number $N=50$ and at the proton
numbers $Z=$ 28, 50 and 82 are clearly visible.   For unbound single-particle states, on the other hand,
the single-particle energy spacings are much smaller than those of bound states, and the spacings approach zero
if we take larger Wigner-Seitz radius.
We have examined wave functions of all the unbound
single-particle states of neutrons. Among the discretized unbound states, we find a few characteristic states
whose wave function has large amplitude inside the nuclear cluster ($r \m@thcombine<\sim 6$ fm), at least larger than amplitude
in the exterior region.
Such states correspond to resonant single-particle states and they are marked with thick lines in Fig.~\ref{spe}(a).
The resonant single-particle states are high-$l$ orbits with no radial node inside the cluster,
which can be labeled as $1h_{11/2}$, $1h_{9/2}$,
$1i_{13/2}$ etc.
Note that these high-$l$ orbits have  large
centrifugal barriers which allow
resonances. The
barrier height for $l=5$ partial waves ($h_{11/2,9/2}$) is about 10 MeV in the case of $\lambda_n=5.0$ MeV.

We found that these resonant high-$l$ orbits are the origin
of the shell effect seen in the cluster shape vibration.

Thick red bars in Fig.~\ref{spe}(a) indicate the neutron Fermi energy (chemical potential) $\lambda_n$.
As $\lambda_n$ increases, the relative position of $\lambda_n$ moves upward
 with respect to the spectrum of unbound
neutron single-particle states. Then, the neutron Fermi energy $\lambda_n$ crosses
with the resonant $1h_{11/2}$ neutron state
 at $\lambda_n \approx 1.8$ MeV while the other resonant
$1h_{9/2}$ and $1i_{13/2}$ orbits crosses with $\lambda_n$
at $\lambda_n \approx 5.0$ and 5.5 MeV, respectively.
It is seen that the softening and the instability of the cluster shape vibration mode occurs
when these crossings occur.

With these observations, we infer the mechanism of the shell effect seen in the
quadrupole shape vibration of cluster. A resonant high-$l$ orbit has wave function
concentrated  in the region of the nuclear cluster, and hence it
may play a similar role of a single-$j$ shell in nuclear structure problems.
Furthermore we can assume that non-resonant unbound single-particle orbits do not play
major roles as they are single-particle states contributing mainly to neutron gas, but not to
the nuclear cluster.
In this situation  enhancement of the quadrupole collectivity occurs when  the
Fermi energy is located close to the energy of the single-$j$ shell (the
resonant high-$l$ state in the present case)  and the shell
is approximately half-filled.  This shell effect  caused by the resonant single-particle states
is very different from the conventional
one, in which
the shells refer to bound orbits in nuclei and the associated magic numbers
20, 28, 50, 82, etc. It also differs from the
continuum shell effect\cite{BulgacMagierski2001,MagierskiHeenen2002} which is
present in the level density of  the non-resonant unbound states.
 To distinguish it from the existing concepts of the shell effect,  we
shall call it {\it resonance shell effect}. Note that Baroni et al. \cite{Baroni2010} points out
importance of neutron resonant high-$l$ single-particle states as one of
major contributors to the low-lying and high-lying quadrupole modes. Their discussion
is related to the present resonance shell effect,  but the roles of the Fermi energy
and the pair correlation need to be noticed to explain
the oscillatory  behaviour of  the quadrupole
collectivity as a function of the neutron Fermi energy.

The strong and oscillatory $\lambda_n$-dependence of the
cluster-vibration-dominant mode in the $Z=$ 20, 40 and 50 systems, shown in
Fig.~\ref{ABmode.Zall}, is also nicely explained with the resonance shell effect.
The softening of the shape vibration around $\lambda_n \approx$ 3.3 MeV and 4.8 MeV
in the $Z=20$ system can be related to $1g_{7/2}$
and $1h_{11/2}$ resonant states of neutrons, respectively.
For the softening around $\lambda_n \approx$ 1.4 MeV, 2.6 MeV  and 5.4 MeV
in the $Z=40$ system,  neutron $1h_{9/2}$, $1i_{13/2}$, and $1i_{11/2}$ resonant states, respectively,
are the keys. The instability around $\lambda_n \approx$ 3.7 MeV in
 the $Z=50$ system is due to neutron $1i_{11/2}$  resonant state.

In the present study, we have analyzed the systems with $Z=$ 20, 28, 40 and 50
in which the protons have closed-shell configurations. In realistic calculations of
inner crust\cite{Negele73,Baldo2005,Baldo2006,Grill11}, the predicted proton number of
nuclear cluster is not limited to these
magic numbers, but rather the protons often have open shell configurations
where the proton pair correlation  also takes part in. In such cases, the
quadrupole collectivity of the cluster shape vibration is even more enhanced
due to the contribution of open shell protons
than in the present case of proton closed shell, and there will be more possibilities
of the quadrupole deformation of nuclear clusters in situations where
the resonance shell effect is expected. It is also remarked that
the occurrence of quadrupole-deformed
nuclear cluster caused by the resonance shell effect is expected even
at low densities (as well as high densities) of superfluid neutron gas.
This is in contrast to the deformation mechanism originating from the
Coulomb lattice energy, which is responsible for
the manifestation of various pasta structures\cite{Pethick-Ravenhall95,
RavenhallPethick1983, HashimotoSeki1984, Oyamatsu1993}. It is interesting to
explore deformed clusters at low neutron-gas densities, $0 < \rho_{n,{\rm ext}} \m@thcombine<\sim \rho_0/10$,
by means of
quantum mechanical approach such as the HFB and the
Bogoliubov-de-Genne-Kohn-Sham density functional models. Note that the previous HF calculations
analyzing  deformed clusters\cite{MagierskiHeenen2002,GegeleinMuther2007,NewtonStone2009, PaisStone2012, FattoyevHorowitz2017}
focus on the pasta structures in the bottom layers of inner crust where the neutron-gas density or
the average baryon density is rather high $\rho \m@thcombine>\sim \rho_0/10$.

\section{Conclusions}

We have investigated collective excitations in the inner crust
of neutron stars in the framework of the nuclear density functional theory.
We consider configurations where spherical clusters are immersed in superfluid neutron gas, and
describe low-lying quadrupole collective excitations around the cluster by means of
the HFB plus QRPA model formulated in a spherical Wigner-Seitz cell.
We adopt the Skyrme SLy4 functional and the density-dependent delta pairing interaction.
Numerical analysis is performed by varying the neutron chemical potential (Fermi energy)
$\lambda_n =1 - 6$ MeV changing the density of superfluid neutrons and
the proton number of the cluster.

Our calculation demonstrates that there coexist the Anderson-Bogoliubov (AB) phonon mode of neutron superfluid and the
quadrupole shape vibration of the cluster, and these two fundamental modes mix with each other only weakly.
The AB-phonon-dominant-mode, i.e.  the AB phonon mode with small mixing of the
quadrupole shape vibration of cluster, emerges systematically, and its excitation energy exhibits
smooth and monotonic dependence on $\lambda_n$, which agrees fairly well with the hydrodynamical estimate
of the AB phonon in uniform neutron superfluid.
The results point to the weak coupling
between the quadrupole AB phonon mode of neutron gas and the quadrupole shape vibration of clusters.
This is  consistent with our previous study\cite{Inakura17} suggesting that
the dipole AB phonon and the displacement motion of clusters is weak.

Another remarkable finding is that the excitation energy of the cluster shape vibration varies in the range of
$0-3$ MeV with strong and oscillatory $\lambda_n$-dependence, in contrast to the AB phonon mode.
This behavior is traced to a new kind
of shell effect, named the resonance shell effect, which is associated with
unbound but resonant single-particle states of neutrons orbiting around the cluster with large orbital angular momentum.
The resonance shell effect causes enhancement of collectivity, softening and  even instability of the quadrupole shape  degrees
of freedom of clusters. This suggests that occurrence of quadrupole-deformed nuclear cluster is possible in any layers of inner crust
if the resonance shell effect induces the instability.

\section*{Acknowledgments}

This work was supported by JSPS KAKENHI Grant Number JP17K05436, and by Grant-in-Aid for Scientific
Research on Innovative Areas, No. 24105008, by The Ministry of Education,
Culture, Sports, Science and Technology, Japan. \\


\section{Supplmental material}

\renewcommand{\theequation}{S\arabic{equation}}
\renewcommand{\thefigure}{S\arabic{figure}}
\setcounter{equation}{0}
\setcounter{figure}{0}

\subsection{Numerical method to solve the Schr\"odinger equation}
The state vector at time $t$ is given by
$|\Psi(t)\rangle = \hat{U}(t) |\Psi(0)\rangle$ with
the unitary time evolution operator
\begin{equation}
\hat{U}(t) \equiv e^{-i \hat{H}t} .
\end{equation}
When the whole set of eigenvectors $\{|a\rangle\}$ of
$\hat{H}$ is available,
it is given by $|\Psi(t)\rangle = \sum_a c_a e^{-i E_a t}
| a \rangle$ with $c_a = \langle a|\Psi(0)\rangle$.
This is the exact method.

\begin{figure}[h]
\includegraphics*[width=0.5\columnwidth]{sfig1.eps}
\caption{Illustration of the XXZ Hamiltonian with the nearest and next
nearest interactions under the open boundary condition.}
\label{sfig1}
\end{figure}

When the dimension of the Hilbert space is too large, the exact method is
impractical and one needs to rely on an approximate method.
We will explain an efficient method to simulate the time evolution.
The XXZ Hamiltonian of $N$ spins under the open boundary condition
is given by~(see Fig.~\ref{sfig1})
\begin{equation}
\hat{H} = \sum_{l=1}^{N-1} \hat{h}_{l,l+1} + \sum_{l=1}^{N-2} \hat{h}_{l,l+2}
\end{equation}
with the short-hand notation
\begin{equation}
\begin{aligned}
\hat{h}_{l,m} \equiv h(\hat{\sigma}_l,\hat{\sigma}_m) &= -\frac{J}{2} \left( \hat{\sigma}_l^x\hat{\sigma}_m^x +
                                      \hat{\sigma}_l^y\hat{\sigma}_m^y +
                        \Delta \hat{\sigma}_l^z \hat{\sigma}_m^z \right) \\
  & = -J \left( \hat{\sigma}_l^+ \hat{\sigma}_m^- +
                \hat{\sigma}_l^- \hat{\sigma}_m^+
      + \frac{\Delta}{2} \hat{\sigma}_l^z \hat{\sigma}_m^z \right) ,
\end{aligned}
\end{equation}
where $\hat\sigma_l^{\pm} = (\hat\sigma_l^x \pm i \hat\sigma_l^y)/2$.

Note that the Hamiltonian is slightly different from that considered in the
paper.
Generalizations to the Hamiltonian with bond-dependent $J$ and $\Delta$
are straightforward.

The time evolution operator satisfies $\hat{U}(\epsilon M) =
\hat{U}(\epsilon)^{M}$. Thus, the state vector at $t = \epsilon M$
is obtained by multiplying $\hat{U}(\epsilon)$ to the initial
state vector $M$ times repeatedly.
For infinitesimal $\epsilon$, $\hat{U}(\epsilon)$ can be approximated
efficiently.  First, we decompose the Hamiltonian as the sum
$\hat{H} = \hat{H}_A + \hat{H}_B + \hat{H}_C + \hat{H}_D$ with
\begin{equation}\label{decomp}
\begin{split}
\hat{H}_A & =  \sum_{l} \hat{h}_{2l-1,2l} \\
\hat{H}_B & =  \sum_{l} \hat{h}_{2l,2l+1} \\
\hat{H}_C & =  \sum_{l} ( \hat{h}_{4l-3,4l-1} + \hat{h}_{4l-2,4l})\\
\hat{H}_D & =  \sum_{l} ( \hat{h}_{4l-1,4l+1} + \hat{h}_{4l,4l+2} ) .
\end{split}
\end{equation}
In Fig.~\ref{sfig1}, the interaction bonds contained in
$\hat{H}_A$, $\hat{H}_B$, $\hat{H}_C$, and $\hat{H}_D$
are represented by the solid, dashed, dotted, and dash-dotted lines,
respectively.
Then, using the generalized Lie-Suzuki-Trotter formula of the second order
$e^{\epsilon(\hat{X}+\hat{Y})} = e^{\frac{\epsilon}{2}\hat{X}} e^{\epsilon
\hat{Y}} e^{\frac{\epsilon}{2}\hat{X}} +
O(\epsilon^3)$~[30], we obtain
\begin{equation}
\hat{U}(\epsilon) =  e^{-i\frac{\epsilon}{2}\hat{H}_A}
                     e^{-i\frac{\epsilon}{2}\hat{H}_B}
                     e^{-i\frac{\epsilon}{2}\hat{H}_C}
                     e^{-{i\epsilon}\hat{H}_D}
                     e^{-i\frac{\epsilon}{2}\hat{H}_C}
                     e^{-i\frac{\epsilon}{2}\hat{H}_B}
                     e^{-i\frac{\epsilon}{2}\hat{H}_A} + O(\epsilon^3)
\end{equation}

The important feature of the decomposition \eqref{decomp} is that all
$\hat{h}_{l,m}$'s contained in a partial Hamiltonian commute with each
others. It yields
\begin{equation}
e^{-i\gamma \hat{H}_A} = \prod_{l} e^{-i \gamma \hat{h}_{2l-1,2l}}
\end{equation}
and similar relations for the others. Consequently, $\hat{U}(\epsilon)$ can
be approximated as the ordered product of $\hat{U}_{l,m}(\gamma) = e^{-i \gamma
\hat{h}_{l,m}}$ with $\gamma=\epsilon$ or $\epsilon/2$ with an error $O(\epsilon^3)$.
If one chooses
$\{ |-\rangle_l\otimes |-\rangle_m,
|-\rangle_l\otimes |+\rangle_m,
|+\rangle_l\otimes |-\rangle_m,
|+\rangle_l\otimes |+\rangle_m \}$ as the basis states for the spins at
sites $l$ and $m$, $\hat{U}_{l,m}(\gamma) = \hat{u}_{l,m}(\gamma) \otimes
\mathsf{I}_{N-2}$ with the $4\times 4$ unitary matrix
\begin{equation}
\hat{u}_{l,m}(\gamma) =
\begin{pmatrix}
e^{i \phi } & 0 & 0 & 0 \\
0 & e^{-i\phi} \cos\gamma J& - i e^{-i\phi} \sin \gamma J & 0 \\
0 & -i e^{-i\phi} \sin\gamma J&  e^{-i\phi} \cos \gamma J & 0  \\
0 & 0 & 0 & e^{i \phi}
\end{pmatrix}
\end{equation}
with $\phi = \gamma J \Delta /2$.
The overall error of this method is of the order of $O(\epsilon^2)$.

Computationally, the infinitesimal time evolution is achieved by the
multiplication of a state vector with the $4\times 4$ matrices $O(N)$ times.
All the operators $\hat{U}_{l,m}$'s conserve the total number of up spins.
Thus, one can work within the subspace consistent with the conservation,
which make numerical calculations efficient.
In the numerical calculations, we used the exact time evolution when $N\leq
16$, and the decomposition method for $N>16$ with $\epsilon = 0.01$. We
examined whether $\epsilon = 0.01$ is small enough. As can be seen in
Fig.~\ref{sfig2}, the data with $\epsilon=0.1, 0.05$, and $0.01$ are
indistinguishable in the plot. More quantitatively, the values of
$Q_{\rm L}$ at time $t=10^3$ are $6.004072339$, $6.003279566$, and
$6.003025718$ for $\epsilon=$ $0.1$, $0.05$, and $0.01$, respectively.
They are well fitted to the function $Q_{\rm L} = Q_0 + a \epsilon^{1.99926}$ with
$Q_0 = 6.003015128$. The relative error of the value at $\epsilon=0.01$ is
$\simeq 1.7\times 10^{-6}$. It proves that $\epsilon=0.01$ is small enough.

\begin{figure}
\includegraphics*[width=0.5\columnwidth]{sfig2.eps}
\caption{Numerical data for the energy density and the particle density for
the subsystem L using the same condition as Fig.~3 of the main text. The
data obtained by using $\epsilon=0.1, 0.05,$ and $0.01$ cannot be
distinguished.}
\label{sfig2}
\end{figure}

\subsection{Singular value decomposition for reduced density matrix}
Suppose that one separates the $N$ spins into two sets $A$ of $N_A$ spins
and $B$ of $N_B = N-N_A$ spins. We assume that $N_A \leq N_B$ without loss
of generality.
Then, the state vector can be written as
\begin{equation}
|\Psi\rangle = \sum_{a,b} \psi_{a,b} |a\rangle\otimes |b\rangle ,
\end{equation}
where $\{|a\rangle\}$ and $\{|b\rangle\}$ span the Hilbert spaces of
dimension $d_A=2^{N_A}$ and $d_B = 2^{N_B}$ for spins in $A$ and $B$,
respectively. The reduced density operator for the spins in $A$ is
given by
\begin{equation}
\hat{\rho}_A = {\rm Tr}_B |\Psi\rangle\langle \Psi| = \sum_{b}
\psi_{a,b}\psi_{a',b}^* |a\rangle \langle a'|
\end{equation}
and represented by the $({d_A}\times {d_A})$ reduced density matrix
is given by
\begin{equation}
(\rho_A )_{a,a'} = \sum_b \psi_{a,b} \psi_{a',b}^* .
\end{equation}
The reduced density matrix is necessary when one calculates the local
observables and the probability distribution for the spins in $A$.
The eigenvalues are necessary
for the entanglement entropy. The singular value decomposition~(SVD) is a
convenient method to handle the reduced density matrix.

The probability amplitudes $\psi_{a,b}$ can be regarded as the
elements of a $({d_A}\times {d_B})$ positive semidefinite complex matrix
$M$, $M_{ab} =
\psi_{a,b}$. According to the SVD, such a matrix $M$ can be written as
\begin{equation}
M = V_1 \Sigma V_2^\dagger ,
\end{equation}
where $V_1$ is a $({d_A}\times {d_A})$ unitary matrix, $\Sigma$ is a $({d_A}
\times {d_B})$ rectangular diagonal matrix with non-negative real diagonal
elements, and $V_2$ is a $({d_B}\times {d_B})$ unitary matrix.
The SVD can be done with a numerical analysis software.

Once the SVD is done, the reduced density matrix is obtained from the
relation
\begin{equation}
\rho_A = M M^\dagger = V_1 \Sigma V_2^\dagger V_2 \Sigma^\dagger V_1^\dagger
       = V_1 \Sigma \Sigma^\dagger V_1^\dagger .
\end{equation}
That is, the reduced density matrix is given by the similarity transformation
of the diagonal $(d_A\times d_A)$ matrix
\begin{equation}
\Sigma \Sigma^\dagger = {\rm
diag}(|\lambda_1|^2,|\lambda_2|^2,\cdots,|\lambda_{d_A}|^2) .
\end{equation}
The SVD allows one to calculate the entanglement entropy easily.
It is given by
\begin{equation}
S_{ent} = -\sum_{n=1}^{d_A} |\lambda_n|^2 \ln |\lambda_n|^2 .
\end{equation}

\subsection{Temperature-energy relation in the $Q=1$ sector}
Consider the XXZ Hamiltonian of $N$ spins under the periodic boundary
condition. In the $Q=0$ sector where all spins are down~(or all sites are
empty), there is a single eigenstate with the energy
\begin{equation}
E_{0,1} = -\frac{\Delta}{2} N .
\end{equation}
The Hilbert space of the $Q=1$ sector is spanned by the states $\{|1\rangle,
\cdots, |x\rangle,\cdots,|N\rangle\}$ where $|x\rangle$ represents the state
where the particle is at site $x$. The plane
wave $|k_n\rangle = \frac{1}{\sqrt{N}} \sum_x e^{ik_n x} |x\rangle$
is the eigenstate and the energy eigenvalue is given by
\begin{equation}
E_{1,k_n} = -\frac{\Delta}{2}(N-4) - \frac{2}{1+\lambda} \left( \cos k_n +
\lambda \cos 2k_n \right) = -\frac{\Delta}{2}(N-4) +
\varepsilon(k_n)
\end{equation}
with $k_n = 2\pi n / N$~($n=0,\cdots,N-1$). The partition function in the
$Q=1$ sector is given by
\begin{equation}
Z_1(\beta) = 2^{\beta \Delta (N-4)/2}\sum_{n=0}^{N-1}  e^{-\beta
\varepsilon(k_n)}
\end{equation}
with the mean energy
\begin{equation}
\langle E\rangle = - \frac{\partial}{\partial \beta} \ln Z_1(\beta)
= -\frac{\Delta}{2}(N-4) + \frac{\sum_n \varepsilon(k_n) e^{-\beta
\varepsilon(k_n)}}{\sum_n e^{-\beta \varepsilon(k_n)}} .
\end{equation}

We estimate the location of the inversion temperature $T_{\rm I}$
by requiring that
the mean energy in the $Q=1$ sector should be equal to the energy in the $Q=0$
section. It yields that
\begin{equation}
-2 \Delta =  \frac{\sum_n \varepsilon(k_n) e^{-\beta_{\rm I}
\varepsilon(k_n)}}{\sum_n e^{-\beta_{\rm I} \varepsilon(k_n)}}
\end{equation}
with $\beta_{\rm I} = 1/T_{\rm I}$. In the $N\to\infty$ limit, the summation is replaced
by the integration to yield the expression
\begin{equation}
-2\Delta = \frac{ \int_{-\pi}^\pi dk \ \varepsilon(k) e^{-\beta_{\rm I}
\varepsilon(k)}}{ \int_{-\pi}^\pi dk \ e^{-\beta_{\rm I} \varepsilon(k)}}  \equiv
\langle \varepsilon(k)\rangle_{\beta_{\rm I}}.
\end{equation}

Note that $\langle \varepsilon(k)\rangle_{\beta_{\rm I}}$ increases as
$\beta_{\rm I}
\to 0^+$ or $T_{\rm I} \to +\infty$. For small $\beta_{\rm I}$, we can approximate
\begin{equation}
\langle \varepsilon(k)\rangle_{\beta_{\rm I}} \simeq -\beta_{\rm I} \int \frac{dk}{2\pi}
\varepsilon(k)^2  = -2 \beta_{\rm I} \frac{1+\lambda^2}{(1+\lambda)^2} .
\end{equation}
Thus, the inversion temperature has the limit form
\begin{equation}
T_{\rm I} \simeq \frac{1+\lambda^2}{(1+\lambda)^2} \Delta^{-1}
\end{equation}
as $\Delta \to 0^{+}$.
In the opposite limit where $\beta_{\rm I} \to \infty$~($T_{\rm I}\to 0^+$), $\langle
\varepsilon(k)\rangle_{\beta_{\rm I}}$ is dominated by the minimum value of
$\varepsilon(k)$. For positive $\lambda$, $\varepsilon(k)$ takes the minimum
value $-2$ at $k=0$. It yields that $T_{\rm I} \simeq (1-\Delta)/4$.

\subsection{Joule expansion}
In the main text, we present the Joule expansion of a gas into the
vacuum. We have performed the numerical simulation of the Joule expansion
in a mixture of a dense gas and a dilute gas.
The subsystems ${\rm L}$ and ${\rm R}$ is prepared in the energy eigenstates
$|Q_{\rm L}^0 =
q_{\rm L} N, n_{\rm L}\rangle$ and $|Q_{\rm R}^0 = q_{\rm R} N, n_{\rm R}\rangle $, respectively,
with $q_{\rm L} = q = 1/4$ and $q_{\rm R} = 1-q = 3/4$ with the same quantum
number $n_{\rm L}
= n_{\rm R} = n_0$.
The two states are tied to each other by the particle-hole symmetry.
We select those two states in order to keep them in the same initial
temperature.

When the total system reaches a stationary state, we measure the probability
distribution $P_{\rm L}(Q,n)$ that the subsystem ${\rm L}$ is in the state
$|Q,n\rangle$. Then, we fit the probability distribution to the function
$P_{\rm L}(Q,n) \propto \exp[ -(E_{Q,n}-\overline{E_{\rm L}})/T_f - a_{11}
(E_{Q,n}-\overline{E_{\rm L}})^2]$ at the most probable sector $Q=N/2$ to find the
final temperature $T_f$. The numerical result is presented in
Fig.~\ref{sfig3}. We also observe that the gas heats up (cools down) when
the initial temperature is lower (higher) than a threshold temperature.
The finite size effect is also seen in the plot.

\begin{figure}
\includegraphics*[width=0.5\columnwidth]{sfig3.eps}
\caption{Heating/cooling curve when the dense gas of density $\rho_{\rm R}=3/4$ is
mixed with the dilute gas of density $\rho_{\rm L} = 1/4$. The filled symbols
represent the numerical results obtained by fitting the probability
distribution. The open symbols represent the expected result from
the isotherm curve.}
\label{sfig3}
\end{figure}

\end{document}


\section{Introduction}

In finite Group Theory, the nilpotency of a group can be analyzed looking at the behavior of its elements of coprime orders. Indeed, in \cite{baubau} B.\,Baumslag and J.\,Wiegold proved the following result.

\begin{theorem}[\cite{baubau}]\label{teobau} Let $G$ be a finite group in which $|ab|=|a||b|$ whenever the elements $a,b$ have coprime orders. Then $G$ is nilpotent.
\end{theorem}

 Here the symbol $|x|$ stands for the order of the element $x$ in a group $G$. In this direction, in 2016, R.\,Bastos and P.\,Shumyatsky established a characterization for the nilpotency of the commutator subgroup  of a finite group \cite{BS}.

\begin{theorem}[\cite{BS}]\label{teoBS} Let $G$ be a finite group in which $|ab|=|a||b|$ whenever the elements $a,b$ are commutators of coprime orders. Then $G'$ is nilpotent.
\end{theorem}

 This theorem represents a first successful attempt to broaden Theorem~\ref{teobau} to the world of group-words. Recall that a {\it word} (or {\it group-word}) $w=w(x_1, \ldots, x_k)$ is a nontrivial element of the free group $F=F(x_1, \ldots, x_k)$ on free generators $x_1, \ldots, x_k$. A word is a {\it commutator word} if it belongs to the commutator subgroup $F'$ of $F$.

 Given a word $w$ and a group $G$, we denote by $G_w$ the set of all $w$-values, and the subgroup of $G$ generated by $G_w$ is called the {\it verbal subgroup of $G$ corresponding to $w$}, and it is usually denoted by $w(G)$.  An example of commutator word is the $k$th lower central word $\gamma_k=~\gamma_k(x_1, \ldots, x_k)$ defined inductively by the formulae
\[
\gamma_1=x_1
\qquad \text{and} \qquad
\gamma_{k+1}=[\gamma_{k},x_{k+1}]
\quad
\text{for $k\ge 1$.}
\]

 As one can aspect, the verbal subgroup corresponding to $\gamma_k$ is the familiar $k$th term of the lower central series $\gamma_k(G)$, and the $\gamma_k$-values in $G$ are called $\gamma_k$-commutators.

Using the language of group-words, Theorem~\ref{teobau} and Theorem~\ref{teoBS}  show that the  nilpotency of the verbal subgroups $G$ and $G'$ can be obtained imposing the condition about orders on the sets $G_x$ and $G_{[x,y]}$, respectively. This suggests the question whether a similar criterion holds for other group-words.

\begin{qst}\label{qGen}
Let $w$ be a group-word and let $G$ be a finite group in which $|ab|=|a||b|$ whenever the $w$-values $a,b$ have coprime orders. Is then $w(G)$ necessarily nilpotent?
\end{qst}

Easy counterexamples show that Question \ref{qGen} has negative answer in general. Indeed, one can consider the symmetric group $G$ of degree $3$ and the word $w=x^3$. Then the nontrivial $w$-values have order $2$, but $w(G)=G$.

Also in the case of commutator words Question \ref{qGen} has negative answer since one can consider the alternating group $\mathbb{A}_5$ of degree $5$ and the word $[x,y^{10},y^{10},y^{10}]$ all of whose nontrivial values in $\mathbb{A}_5$ have order $2$ (see \cite{MT}).

Despite all the underlined negative answers, C.\,Monetta and A.\,Tortora proved that there is a family of commutator words which answers in the affermative to Question~\ref{qGen}.

\begin{theorem}[\cite{MT}]\label{teoMT}
Let $w=w(x_1, \ldots, x_k)$ be the simple commutator word defined by
\[w=[x_{i_1}, \ldots, x_{i_n}],\]
where $i_{1} \neq i_j$ for every $j \in \{2, \ldots, n\}$ and $i_j \in \{1, \ldots, k\}$ for every $j$. If $G$ is a finite group, then $w(G)$ is nilpotent if and only if $|ab|=|a||b|$ for any $w$-values $a,b\in G$ of coprime orders.
\end{theorem}

Notice that $\gamma_n$-commutators are particular instances of simple commutator words.
Indeed, Theorem~\ref{teoMT} extends the following result proved by R.\,Bastos, C.\,Monetta and P.\,Shumyatsky in \cite{bms}.

\begin{theorem}[\cite{bms}]\label{metanilp} Let $k$ be a positive integer. The $k$th term of the lower central series of a finite group $G$ is nilpotent if and only if $|ab|=|a||b|$ for any $\gamma_k$-commutators $a,b\in G$ of coprime orders.
\end{theorem}

This result contributes towards a positive answer to the following conjecture made in \cite{BS}, which remains still an open problem.

\begin{conjecture}\label{shumy} Let $w$ be a multilinear commutator word and let $G$ be a finite group in which $|ab|=|a||b|$ whenever the $w$-values $a,b$ have coprime orders. Then $w(G)$ is nilpotent.
\end{conjecture}

Recall that a multilinear commutator word is a word obtained nesting commutators but using always different variables. However, not all simple commutator words are multilinear. Indeed, another example of simple commutator words is the $k$th Engel word $\varepsilon_k=~\varepsilon_k(x,y)$ defined inductively by the formulae
\[
\varepsilon_1=[x,y]
\qquad \text{and} \qquad
\varepsilon_{k+1}=[\varepsilon_k,y]=[x, \overset{k}{y, \ldots, y}]
\quad
\text{for $k\ge 1$.}
\]

Therefore, from Theorem~\ref{teoMT} it follows that Question~\ref{qGen} has a positive answer for Engel words, too.

\bigskip

A similar question could be asked for commutators of elements of coprime orders, namely for the coprime commutators $\gamma^*_k$ and $\delta^*_k$.

The coprime commutators $\gamma^*_k$ and $\delta^*_k$  were introduced by P.\,Shu\-myatsky in \cite{pavel} as a tool to study the nilpotency and the solubility of finite groups in terms of commutators of elements of coprime orders. They are defined as follows.

Let $G$ be a finite group. Every element of $G$ is both a $\gamma^*_1$-commutator and a $\delta^*_0$-commutator. Now let $k \geq 2$ and let $X$ be the set of all elements of $G$ that are powers of $\gamma^*_{k-1}$-commutators. An element $g \in G$ is a $\gamma^*_k$-commutator if there exist $a \in X$ and $b \in G$ such that $g=[a,b]$ and $(|a|,|b|)=1$. For $k \geq 1$ let $Y$ be the set of all elements of $G$ that are powers of $\delta^*_{k-1}$-commutators. An element $g \in G$ is a $\delta^*_k$-commutator if there exist $a, b \in Y$ such that $g = [a, b]$ and $(|a|, |b|) = 1$. The subgroups of G generated by all $\gamma^*_k$-commutators and all $\delta^*_k$-commutators are denoted by $\gamma^*_k(G)$ and $\delta^*_k(G)$, respectively.

One can easily see that if $N$ is a normal subgroup of $G$ and $x$ is an element of $G$ whose image in $G/N$ is a $\gamma^*_k$-commutator (respectively a $\delta^*_k$-commutator), then there exists a $\gamma^*_k$-commutator (respectively a $\delta^*_k$-commutator) $y$ in $G$ such that $x \in yN$.

The commutators $\gamma^*_k$ and $\delta^*_k$ are closely related to the commutators $\gamma_k$ and $\delta_k$, that is, they control the nilpotency and the solubility of a group as $\gamma_k$ and $\delta_k$ do. Indeed, in \cite{pavel}, it was shown that $\gamma^*_k(G)=1$ if and only if $G$ is nilpotent, while $\delta^*_k(G)=1$ if and only if $G$ is soluble and the Fitting height of $G$ is at most $k$.

Coprime commutators have been studied by many authors (for example see \cite{as, ast, FC, monak, ps}).  In \cite{monak} it is established that the nilpotent residual $\gamma_{\infty}(G)$ of a finite group $G$ is generated by commutators of primary elements of coprime orders, where a primary element is an element of prime power order. Moreover, in the same paper, V.\,S.\,Monakhov proved that the nilpotent residual $\gamma_{\infty}(G)$ of a finite soluble group $G$ is nilpotent if and only if $|ab| \geq |a| |b|$ for every $a$ and $b$ commutators of primary elements of coprime orders.

Similarly, in \cite{FC} A.\,Freitas de Andrade and A.\,Carrazedo Dantas studied the nilpotency of $\gamma_{\infty}(G)$ looking at the orders of powers of $\delta^*_1$-commutators, where $G$ is a finite group. More precisely they proved that the nilpotent residual $\gamma_{\infty}(G)$ of a finite group $G$ is nilpotent provided that $|ab| = |a||b|$ whenever $a,b$ are powers of $\delta^*_1$-commutators of coprime orders. In \cite{FC}, they also posed the following question.

\begin{qst}\label{alexi}
Let $k$ be a non-negative integer and let $G$ be a finite group in which $|ab| = |a||b|$ whenever the elements $a,b$ are (powers of) $\delta^*_k$-commutators of coprime orders. Is then the subgroup $\delta^*_k(G)$ nilpotent?
\end{qst}

\noindent When $k=0$, the answer is positive by Theorem~\ref{teobau}.
At the same time, one may think to restate Question~\ref{alexi} in terms of $\gamma_k^*$-commutators.

Our first main result is the following, which gives a complete answer for $\gamma^*_k$-commutators.

\begin{thmA}\label{mainGk}
Let $G$ be a finite group and let $k \geq 2$. Then $\gamma_k^*(G)$ is nilpotent if and only if $|ab|=|a| |b|$ whenever $a$ and $b$ are $\gamma_k^*$-commutators such that $(|a|,|b|)=1$.
\end{thmA}

 In particular, Theorem A applies to $\gamma_{\infty}(G)$, showing that the nilpotent residual $\gamma_{\infty}(G)$ of a finite group $G$ is nilpotent if and only if $|ab|=|a| |b|$ whenever $a$ and $b$ are $\gamma_k^*$-commutators of coprime orders, for some $k \geq 2$. Moreover, since $\gamma^*_2$-commutators and $\delta^*_1$-commutators coincide, Theorem A shows that it is enough to impose the condition on the orders of $\delta^*_1$-commutators without considering their powers, and that a group satisfying such a condition is necessarily soluble, improving both results in \cite{monak} and \cite{FC}.

Concerning the case when $k \geq 2$, we answer in the positive to Question~\ref{alexi} imposing the condition on the orders of powers of $\delta^*_k$-commutators. Indeed we succeed in proving the following result.

\begin{thmB}\label{mainDk}
Let $G$ be a finite group and let $k \geq 2$. Then $\delta^*_k(G)$ is nilpotent if and only if $|ab|=|a| |b|$ whenever $a$ and $b$ are powers of $\delta^*_k$-commutators such that $(|a|,|b|)=1$.
\end{thmB}

\section{Preliminaries}
As usual, if $\pi$ is a set of primes, we denote by $\pi'$ the set of all primes that do not belong to $\pi$. For a group $G$ we denote by $\pi(G)$ the set of primes dividing the order of $G$, and the maximal normal $\pi$-subgroup of $G$ is denoted by $O_{\pi}(G)$. For elements $x,y$ of a group $G$ we write $[x, \ _0y]=x$ and $[x,\ _{i+1}y]=[[x,\ _{i}y],y]$ for $i\geq0$.
We denote by $F(G)$ the Fitting subgroup of~$G$.

Let $G$ be a soluble group. The {\it upper Fitting series} of $G$ is the series
\[
1=F_0(G) \leq F_1(G) \leq \cdots \leq F_t(G)=G
\]
where $F_1(G)=F(G)$ and $F_{i+1}(G)/F_i(G)=F(G/F_i)$ for every $i \geq 1$.
Similarly, the  {\it lower Fitting series} of $G$ is the series
\[
G=N_0(G) \geq N_1(G) \geq \cdots \geq N_t(G)=1
\]
where $N_{i}(G)=\gamma_{\infty}(N_{i-1}(G))$ for every $i \geq 1$. The number $t$ is called the {\it Fitting height} of the group and is denoted by $h(G)$.

Throughout the article we use without special references the well-known properties of coprime actions:  if $\alpha$ is an automorphism of a finite group $G$ of coprime order, $(|\alpha|,|G|)=1$, then $C_{G/N}(\alpha)=C_G(\alpha)N/N$ for any $\alpha$-invariant normal subgroup $N$, the equality $[G,\alpha]=[[G,\alpha],\alpha]$ holds, and if $G$ is in addition abelian, then $G=[G,\alpha]\times C_G(\alpha)$. Here $[G,\alpha]$ is the subgroup of $G$ generated by the elements of the form $g^{-1}g^\alpha$, where $g\in G$.

For convenience of the reader we state some results which will be used later.

\begin{lemma}[\cite{pavel}, Lemma 2.4]\label{uuu}
Let $G$ be a finite group and $y_1, \ldots, y_k$ be $\delta^*_k$-commutators in $G$. Suppose the elements  $y_1, \ldots, y_k$  normalize a subgroup $N$ of $G$ such that $(|y_i|, |N|)=1$ for every $i=1,\ldots, k$. Then for every $x \in N$ the element $[x,  y_1, \ldots, y_k ]$ is a $\delta^*_{k+1}$-commutator.
\end{lemma}

Let us call a subgroup $H$ of a group $G$ a tower of height $h$ if $H$ can be written as a product $H=P_1\cdots P_h$, where

\begin{enumerate}
\item $P_i$ is a $p_i$-group ($p_i$ a prime) for $i=1,\dots,h$.
\item $P_i$ normalizes $P_j$ for $i<j$.
\item $[P_i,P_{i-1}]=P_i$ for $i=2,\dots,h$.
\end{enumerate}
It follows from (3) that $p_i\neq p_{i+1}$ for $i=1,\dots,h-1$. In \cite{turull}, it was proved that a finite soluble group $G$ has Fitting height at least $h$ if and only if $G$ possesses a tower of height $h$.

\begin{lemma}[\cite{pavel}, Lemma 2.6]\label{coprimetower}
Let $P_1 \cdots P_h$ be a tower of height $h$. For every $1 \leq i \leq h$ the subgroup $P_i$ is generated by $\delta^*_{i-1}$-commutators contained in $P_i$.
\end{lemma}
The following theorem shows how $\delta^*_k$-commutators control the solubility of a finite group.

\begin{theorem}[\cite{pavel}, Theorem 2.7]\label{fittingdelta}
Let $G$ be a finite group and let $k$ be a positive integer. Then $\delta^*_k(G)=1$ if and only if $G$ is soluble with Fitting height at most $k$.
\end{theorem}

As an immediate consequence of Theorem~\ref{fittingdelta}, the next result shows the relationship between the terms of the lower Fitting series and  $\delta^*_k$-commutators.

\begin{lemma}\label{deltagen}
If $G$ is a finite group, then the $k$th term $N_k(G)$ of the lower Fitting series of $G$ equals $\delta^*_{k}(G)$ for every $k \geq 0$.
\end{lemma}

We conclude this section with a result obtained in \cite{pavel} which generalizes the famous Burnside $p^{\alpha}q^{\beta}$-Theorem (see \cite[Theorem 4.3.3]{go}).

\begin{theorem}[\cite{pavel}, Theorem 2.5]\label{prime} Let $k$ be a positive integer, let $\pi$ be a set consisting of at most two primes and let $G$ be a finite group in which all $\delta^*_k$-commutators are $\pi$-elements. Then $G$ is soluble and $\delta^*_k(G) \leq O_{\pi}(G)$.
\end{theorem}

\section{Proof of Theorem A}

We start with a  proposition which clarifies the connection between $\gamma^*_k(G)$ and the terms of the lower central series of a finite group $G$.

\begin{proposition}\label{gammainfk}
If $G$ is a finite group, then $\gamma^*_k(G)=\gamma_{\infty}(G)$ for every $k \geq 2$.
\end{proposition}

\begin{proof}
We argue by induction on $k$. Let $k=2$.
Since $G/\gamma_{\infty}(G)$ is nilpotent, $[a,b] \in \gamma_{\infty}(G)$ for every $a,b \in G$ of coprime orders. Hence $\gamma^*_2(G) \leq \gamma_{\infty}(G)$. On the other hand, if $p,q \in \pi(G)$ with $p \neq q$, then $p$-elements and $q$-elements of $G/\gamma^*_2(G)$ commute. Therefore $G/\gamma^*_2(G)$ is nilpotent and $\gamma_{\infty}(G) \leq \gamma^*_2(G)$.

Now, notice that $\gamma_{\infty}(G)=[\gamma_{\infty}(G),G]$ and that with a similar argument as before, it can be proved that
\[[\gamma_{\infty}(G),G]= \langle [a,b] \ | \ a \in \gamma_{\infty}(G), b \in G, (|a|,|b|)=1 \rangle.\]
If $k \geq 2$, it follows that
\begin{align*}
\gamma^*_{k+1}(G)&=\langle [a,b] \ | \ a \in \gamma^*_k(G), b \in G, (|a|,|b|)=1 \rangle \\
&= \langle [a,b] \ | \ a \in \gamma_{\infty}(G), b \in G, (|a|,|b|)=1 \rangle \\
&= \gamma_{\infty}(G),
\end{align*}
and we are done.
\end{proof}

The next lemma shows the role played by the condition imposed on the orders of $\gamma^*_k$-commutators.

\begin{lemma}\label{ddd}
Let $k\geq 2$ and let $G$ be a finite group in which $|ab|=|a||b|$ whenever $a$ and $b$ are $\gamma_k^*$-commutators of coprime orders. Assume that $N$ is a subgroup of $G$ normalized by  an element $x$ such that $(|N|,|x|)=1$. Then $[N, x]=1$.
\end{lemma}

\begin{proof}
Let $y \in N$. Then $g=[y, \ _{k-1} x^{-1}]$ is a  $\gamma_k^*$-commutator in $N$. Since $(|g|,|x|)=1$, it follows that
\[|gx|=|g| |x|.\]
However,
\[gx=[y, \ _{k-2} x^{-1}]^{-1}x[y, \ _{k-2} x^{-1}],\]
hence $|g|=1$. Since $[N, \ _{k-1} x^{-1}]=[N, x^{-1}]$, it results that $[N,x]=~1$.
\end{proof}

The soluble case follows.

\begin{proposition}\label{soluGk}
If $G$ is a finite soluble group in which $|ab|=|a||b|$ whenever $a$ and $b$ are $\gamma_k^*$-commutators of coprime orders, for some $k\geq~2$, then $\gamma_k^*(G)$ is nilpotent. Moreover $G$ is metanilpotent.
\end{proposition}

\begin{proof}
Let $h$ be the Fitting height of $G$. Since $\gamma_k^*(G)=\gamma_{\infty}(G)$ by Theorem \ref{gammainfk}, it follows that $\gamma^*_k(G)$ is nilpotent when $h \leq 2$.

Assume $h\geq 3$. Then, there exists a tower $P_1P_2P_3\ldots P_h$ of height $h$ in $G$. Since $P_2=[P_2,P_1]$, it follows that $$P_2=[P_2,\underbrace{P_1, \ldots,P_1}_{(k-1)\ times}].$$
Combining Lemma \ref{ddd} with the fact that $P_2$ is generated by $\gamma_k^*$-com\-mutators of $G$ of $p_2$-orders, we deduce that $P_3$ commutes with $P_2$. On the other hand, $[P_3,P_2]=P_3$, because $P_1P_2P_3\ldots P_h$ is a tower. This is a contradiction, and the proof is complete.
\end{proof}

To succeed in completing the proof we need the following lemma.

\begin{lemma} \label{gk}
Let $G$ be a finite group such that $G=G'$ and let $q\in\pi(G)$. Then $G$ is generated by $\gamma^*_k$-commutators of $p$-power order for primes $p\neq q$.
\end{lemma}

\begin{proof}
For each prime $p\in\pi(G)\setminus\{q\}$, denote by $N_p$ the subgroup generated by all $\gamma^*_k$-values of $p$-power order. Let us show that for each $p$ the Sylow $p$-subgroups of $G$ are contained in $N_p$. Suppose that this is false and choose $p$ such that a Sylow $p$-subgroup of $G$ is not contained in $N_p$. Since $N_p \lhd G$, we can pass to the quotient $G/N_p$ and we can assume that $N_p=1$. Now $G=G'$, so $G$ does not possess a normal $p$-complement. Therefore, by \cite[Theorem 7.4.5]{go}, it follows that $G$ has a $p$-subgroup $H$ and a $p'$-element $a\in N_G(H)$ such that $[H,a]\neq1$.
Then, $$1\neq[H,a]=[H,\underbrace{a,\ldots,a}_{k-1\ times}]\leq N_p,$$ a contradiction. Therefore $N_p$ contains the Sylow $p$-subgroups of $G$. Let $T$ be the product of all $N_p$ for $p\neq q$. Then $G/T$ is a $q$-group, and so $G=T$ because $G=G'$. We conclude that $G$ can be generated by $\gamma^*_k$-values of $p$-power order for $p \neq q$, as desired.

\end{proof}

To give the proof of Theorem A we will need some remarks on {minimal simple groups}. A minimal simple group is a nonabelian simple group whose proper subgroups are soluble. They have been classified by Thompson in his famous paper \cite{thompson}.

\begin{theorem}[\cite{thompson}, Corollary 1] \label{thomClass} Every finite minimal simple group is isomorphic
to one of the following groups:
\begin{enumerate}
\item[(1)] $\Psl(2, 2^p)$, where $p$ is any prime;
\item[(2)] $\Psl(2, 3^p)$, where $p$ is any odd prime;
\item[(3)] $\Psl(2, p)$, where $p > 3$ is any prime such that $p^2 + 1 \equiv 0 \mod 5$;
\item[(4)] $\Psl(3, 3)$;
\item[(5)] $\Sz(2^p)$, where $p$ is any odd prime.
\end{enumerate}
\end{theorem}

In light of Theorem \ref{thomClass}, we recall the following result from \cite{MT}.

\begin{proposition}[\cite{MT}]\label{maj}
Every finite minimal simple group $G$ contains a subgroup $H= A \rtimes T$ where $A$ is an elementary abelian $2$-group and $T$ is a subgroup of odd order such that $C_A(T)=1$.
\end{proposition}

We are now in position to give the proof of our first main result.

\begin{proof}[\bf Proof of Theorem A] Since the necessary condition for $\gamma^*_k(G)$ to be nilpotent is obviously satisfied, we only need to prove the sufficient one.

Suppose that the theorem is false and let $G$ be a counterexample of minimal order. In view of Proposition \ref{soluGk} $G$ is not soluble while all proper subgroups of $G$ so are. Indeed, if $H <G$ then $\gamma^*_k(H)$ is nilpotent by minimality of $G$; moreover, $H/\gamma^*_k(H)=H/\gamma_{\infty}(H)$ is nilpotent, thus $H$ is soluble.
Therefore $G=G'$. Let $R$ be the soluble radical of $G$. It follows that $G/R$ is a nonabelian simple group. Moreover, by Proposition \ref{soluGk}, $R$ is metanilpotent. We want to show that $R$ coincides with $Z(G)$, which is obvious if $R=1$. Suppose $R\neq1$ and choose $q\in\pi(F(G))$. According to Lemma \ref{gk}, $G$ is generated by the $\gamma^*_k$-commutators of $p$-power order for primes $p\neq q$. Let $Q$ be the Sylow $q$-subgroup of $F(G)$. By Lemma \ref{ddd} $[Q,x]=1$, for every $\gamma^*_k$-commutator $x$ of $q'$-order. Therefore $Q\leq Z(G)$. This happens for each choice of $q\in\pi(F(G))$, so that $F(G)= Z(G)$.

Now, since $R$ is metanilpotent, there exists a positive integer $s$ such that $\gamma_s(R)$ is nilpotent. Therefore, if $x \in R$ and $y \in G$, we have
\[[[y,x], \underbrace{x,\ldots,x}_{s-1 \ times}]\in\gamma_s(R)\leq F(G)=Z(G).\]
Consequently, $[y, \ _{s+1} x]=1$, and so every element $x$ is Engel. Hence $R\leq F(G)$ by Theorem \cite[12.3.7]{rob}. Thus $R=Z(G)$ and $G$ is quasisimple.

We claim that $G$ contains a $\gamma^*_k$-commutator  $a$ such that $a$ is a $2$-element and $a$ has order 2 modulo $Z(G)$. Since $G/Z(G)$ is minimal simple, Proposition \ref{maj} implies the existence of a subgroup \[H/Z(G) = A/Z(G) \rtimes T/Z(G)\] where $A/Z(G)$ is an elementary abelian $2$-group and $T/Z(G)$ is a group of odd order such that $C_{A/Z(G)}(T/Z(G))=1$.
Since $A$ is nilpotent, its Sylow $2$-subgroup $P$ is normal in $H$ and $A=PZ(G)$. Therefore, \[A/Z(G) = [A/Z(G), \ _{k-1} T/Z(G)],\] and there exists $[x, \ _{k-1} y]Z(G)$ in $A/Z(G)$ with $x$ in $P$ and $yZ(G)$ in $T/Z(G)$. Moreover, since $\langle y \rangle Z(G)$ is nilpotent, we have $\langle y \rangle Z(G)=WZ(G)$ with $W$ of odd order. Hence we may suppose that $y$ has odd order. If $g=[x, \ _{k-1} y]z$, for some $z \in Z(G)$, then $[g,y]$ is a $\gamma_k^*$-commutator of $2$-power order with order $2$ modulo the center.

Finally, fix an element $a$ with the above properties. Since $G/Z(G)$ is a nonabelian simple group, it follows from \cite[Theorem 3.8.2]{go} that there exists an element $t\in G$ such that the order of $[a,t]$ is odd. On the one hand, it is clear that $a$ inverts $[a,t]$. On the other hand, by Lemma \ref{ddd}, $a$ centralizes $[a,t]$. This is a contradiction.
\end{proof}

 In particular, Theorem A applies to $\gamma_{\infty}(G)$, as showed by the following two corollaries.

\begin {corollary}
The nilpotent residual $\gamma_{\infty}(G)$ of a finite group $G$ is nilpotent if and only if $|ab|=|a| |b|$ whenever $a$ and $b$ are $\gamma_k^*$-commutators of coprime orders, for some $k \geq 2$.
\end {corollary}

\begin {corollary}\label{corollarioE}
Let $G$ be a finite group. Then the nilpotent residual $\gamma_{\infty}(G)$ is nilpotent if and only if $|ab|=|a| |b|$ whenever $a$ and $b$ are $\delta^*_1$-commutators of coprime orders.
\end {corollary}

 As an immediate consequence of Theorem~\ref{teoMT} and Corollary~\ref{corollarioE} we obtain the following corollary.

\begin {corollary}
Let $G$ be a finite group. The following properties are equivalent:
\begin{enumerate}
\item[(1)] $|ab|=|a| |b|$ whenever $a$ and $b$ are $\delta^*_1$-commutators of coprime orders;
\item[(2)] the nilpotent residual of $G$ is nilpotent;
\item[(3)] the Fitting height of $G$ is at most $2$;
\item[(4)] there exists a simple commutator word $w=[x_{i_1}, \ldots, x_{i_k}]$ with $i_1 \neq i_j$ for every $j \in \{2, \ldots, k\}$ such that $|xy|=|x||y|$ for any $w$-values $x,y\in G$ of coprime orders.
\end{enumerate}
\end {corollary}

\section{Proof of Theorem B}

We start with the following lemma.
\begin{lemma}\label{commfond}
Let $k$ be a positive integer and let $G$ be a finite group in which $|ab|=|a| |b|$ whenever $a$ and $b$ are $\delta^*_k$-commutators such that $(|a|,|b|)=1$. Let $x$ be a $\delta^*_k$-commutator of $G$ and let $N$ be a subgroup of $G$ normalized by $x$. If $(|N|,|x|)=1$, then $[N, x]=1$.
\end{lemma}

\begin{proof}
By Lemma \ref{uuu}, if $y \in N$ the element $[y,{}_k x]$ is a $\delta^*_k$-commutator in $N$, whose order is prime to the order of $x$. Hence, since $x^{-1}$ is still a $\delta^*_k$-commutator,
\[\displaystyle{|[y,{}_k x] x^{-1}|=|[y,{}_k x]| |x^{-1}|}.\]
However,  $[y,{}_k x] x^{-1}$ is a conjugate of $x^{-1}$, and so $[y,{}_k x]=1$. Therefore, $[N,x]=[N,{}_k x]=1$, and the lemma is proved.
\end{proof}

The following result shows that in the soluble case we achieve the nilpotency of $\delta^*_k(G)$ in weaker hypothesis; in fact, we impose the condition on the orders of $\delta^*_k$-commutators, without considering the powers.

\begin{proposition}\label{main3}
Let $G$ be a finite soluble group and let $k \geq 1$. If $|ab|=|a| |b|$ whenever $a$ and $b$ are $\delta^*_k$-commutators such that $(|a|,|b|)=1$, then $\delta^*_k(G)$ is nilpotent.
\end{proposition}

\begin{proof}
Let $h=h(G)$. By Theorem \ref{fittingdelta}, when $h \leq k+1$ then $\delta^*_{k+1}(G)=1$. Therefore, Theorem \ref{deltagen} implies that $\delta^*_{k+1}(G)=\gamma_{\infty}(\delta^*_k(G))$ and so $\delta^*_k(G)$ is nilpotent.

Now  assume $h \geq k+2$. Thus there exists a tower \[ P_1 \cdots P_{k+2} \cdots P_h\] of $G$ of height $h$. By Lemma \ref{coprimetower},  $P_{k+1}$ is generated by $\delta^*_k$-commutators in $P_{k+1}$. Combining Lemma \ref{commfond} with the fact that $P_{k+1}$ is generated by $\delta^*_k$-commutators of $p_{k+1}$-orders, we deduce that $P_{k+2}$ commutes with $P_{k+1}$. On the other hand, $[P_{k+2},P_{k+1}]=P_{k+2}$, because $P_1\cdots P_h$ is a tower. This is a contradiction, and the proof is complete.
\end{proof}

Now, we need a result concerning the generators of a minimal nonsoluble group. Recall that a group $G$ is {\it minimal nonsoluble} if every proper subgroup of $G$ is soluble, but $G$ is not soluble.

\begin{lemma}\label{genDk}
Let $G$ be a finite minimal nonsoluble group. If $q \in \pi(G)$, then $G$ can be generated by powers of $\delta^*_k$-commutators of $q'$-order.
\end{lemma}

\begin{proof}
Let $T$ be the subgroup of $G$ generated by all powers of $\delta^*_k$-commutators of $q'$-order. We want to prove that $T=G$. Suppose that $T < G$. Then $T$ is soluble and, being normal, we can consider the quotient $G/T$. If $gT$ is a $\delta^*_k$-commutator of $G/T$, we have $|gT|=q^{\alpha}m$, where $q$ does not divide $m$. Hence $(g)^{q^{\alpha}} \in T$, and so every $\delta^*_k$-commutator of $G/T$ has $q$-order. Therefore, by Theorem \ref{prime}, $G/T$ is soluble. In particular $G$ is soluble, a contradiction.
\end{proof}

In addition, we will need the following fact, already stressed in \cite{pavel}, about elements of order $2$ in minimal simple groups. Recall that an element $g$ of a group $G$ is strongly real if there exists an element $x\in G$ of order $2$ such that $g^x=g^{-1}$.

\begin{proposition}\label{invstar}
In a minimal simple group all elements of order $2$ are $\delta^*_n$-commutators,  for every $n \geq 0$.

\end{proposition}

\begin{proof}
Let $G$ be a minimal simple group. It is well-known that all all elements of order $2$ are conjugate in $G$. By Theorem~\ref{thomClass} we have to analyze five cases.\\

Assume that $G=\Psl(2,3^p)$. By Proposition~\ref{maj}, $G$ contains the alternating group of degree $4$ with Frobenius complement $T$ and Frobenius kernel $A$. Then $T$ is normalized by a strongly real element $u$ of prime order dividing either $(3^p-1)/2$ or $(3^p+1)/2$. If $K$ is the subgroup generated by $u$, then $T \rtimes K$ and $C_T(K)=1$. Let $x$ be a generator of $T$. Then
\[x=[x_1, \ _n u],\]
for some $x_1 \in T$ and for every $n \geq 1$. Assume that $j$ is an element of $G$ of order $2$ such that $u^j=u^{-1}$. Hence, if $y \in K$ we have
\[y=[y_1, \ _n j],\]
 for some $y_1 \in K$ and for every $n \geq 1$. The same happens for $A$, that is, when $a \in A$ there exists $a_1 \in A$ such that
\[a=[a_1, \ _n x],\]
for any $n \geq 1$. We prove that $a,x,y$  are  $\delta^*_n$-commutators for every $n \geq 0$. Indeed, since the case $n=0$ is obvious, suppose that $n >0$ and assume that $a, x, y$ are $\delta^*_{n-1}$-commutators. By Lemma \ref{uuu}, $x$ is a $\delta^*_n$-commutator because $x=[x_1, \ _n u]$. Hence, applying again Lemma \ref{uuu}, $a$ is a $\delta^*_{n+1}$-commutator and so are all elements of order $2$ of $G$. Finally, Lemma~\ref{uuu} implies that $y$ is a  $\delta^*_n$-commutator.  \\

Now assume that $G$ is not isomorphic to $\Psl(2,3^p)$. By Proposition~\ref{maj} $G$ contains a subgroup $H=A \rtimes T$ with $A$ elementary abelian $2$-group and $T$ a group of odd order such that $C_A(T)=1$. We can assume $T$ cyclic. If $T$ contains a strongly real element, then the elements of $T$ and all the elements of order $2$ are $\delta^*_n$-commutators for every $n \geq 0$. Assume that this is the case and let $T= \langle x \rangle$ where $x$ is a strongly real element such that $C_A(x)=1$. Since $x$ is strongly real, there exists $j \in G$ of order $2$ such that $x^j=x^{-1}$. Hence $T=[T, j]$ and for every $y\in T$ there exists $y_1 \in T$ such that $y=[y_1, \ _n j]$. Then, for every $n \geq 1$ and every involution $a \in A$ we have
\[a=[b, \ _{n-1}x]\]
for some $b \in A$. We show that $a$ and $x$ are $\delta^*_n$-commutators for every $n \geq 0$. Indeed, suppose that $n \geq 1$ and $x$ is a $\delta^*_{n-1}$-commutator, being clear the case $n=0$. Since $a=[b, \ _{n-1} x]$, Lemma \ref{uuu} implies that $a$ is a $\delta^*_n$-commutator. Since all elements of order $2$ in $G$ are conjugate,  they are $\delta^*_n$-commutators. Therefore, applying again Lemma \ref{uuu} it follows that $x$ is a $\delta^*_{n+1}$-commutator.\\

Let $G$ be isomorphic to the Suzuki group $\Sz(q)$ with $q=2^p$ and $p$ an odd prime. From \cite[Theorem 9]{suzuki2}, $G$ contains a Frobenius group $H= A \rtimes T$ of order $q^2(q-1)$, where A is a Sylow $2$-subgroup of $G$ and $T$ is a cyclic subgroup of order $q-1$ normalizing $A$. Let $x$ be a generator of $T$. By \cite[Proposition 3]{suzuki2}, the normalizer of $T$ in $G$ is a dihedral group $T \rtimes K$ of order $2(q-1)$, and a generator $j \in K$ is such that $x^j=x^{-1}$. \\

Let $G$ be isomorphic to $\Psl(2,2^p)$, with $p$ any prime. By \cite[Theorem 6.25]{suzuki} $G$ contains a subgroup $H=Q \rtimes L$ which is the semidirect product of a Sylow $2$-subgroup $Q$ order $2^p$  and a cyclic group $T$ of order $2^p-1$. Moreover, the normalizer of $T$ in $G$ is dihedral of order $2(2^p-1)$, and $G$ contains strongly real elements of order dividing $2^p-1$ and $2^p+1$.\\

If $G$ is isomorphic to $\Psl(3,3)$, then $G$ contains a copy of the alternating group of degree $4$, with cyclic Frobenius complement $T= \langle x \rangle$, and $x$ strongly real. \\

Finally, if $G$ is isomorphic to $\Psl(2,p)$, with $p$ odd prime, then $G$ contains a copy of the alternating group  of degree $4$, with cyclic Frobenius complement $T$ of order $3$. A generator of $T$ is strongly real. The proof is concluded.
\end{proof}

Now we are ready to prove our second main result.

\begin{proof}[\bf Proof of Theorem B]
Clearly if $\delta^*_k(G)$ is nilpotent, then $|ab| = |a| |b|$ for any power of $\delta^*_k$-commutators $a,b \in G$ of coprime orders. So we only need to prove the converse.

Suppose that the theorem is false and let $G$ be a counterexample of minimal order. In view of Proposition \ref{main3}, $G$ is not soluble while all proper subgroups of $G$ so are. In fact, if $H < G$ then $\delta^*_k(H)$ is nilpotent, and so $\delta_{k+1}^*(H)=1$. Thus, $H$ is soluble by Theorem \ref{fittingdelta}. It follows that $G = G'$.

If $R$ is the soluble radical of $G$, then $G/R$ is a nonabelian simple group. Now we want to prove that $G/ Z(G)$ is simple. We will do it proving that $R= Z(G)$, which is obvious if $R=1$. Suppose $R \neq 1$ and choose $q \in \pi (F(G))$. According to Lemma \ref{genDk} $G$ is generated by power of $\delta^*_k$-commutators of $q'$-order. Let $Q$ be the Sylow $q$-subgroup of $F(G)$. By Lemma \ref{commfond} $[Q, x] = 1$, for every power of $\delta^*_k$-commutator $x$ of $q'$-order. Therefore $Q \leq Z(G)$, and this happens for each choice of $q \in \pi(F(G))$. Hence, we have $F(G) = Z(G)$.

Now suppose that $R$ is not nilpotent,  that is $\gamma_{\infty}(R) \neq 1$. Since $R$ is soluble, the lower Fitting series of $R$,
\[
R=N_0(R) \geq N_1(R) \geq \cdots \geq N_t(R)=1,
\]
has length $t  \geq 2$. It follows that $N_{t-1}(R)$ is a characteristic nilpotent subgroup of $R$, and so $N_{t-1}(R) \leq F(G)=Z(G)$. Therefore $N_{t-2}(R)$ is nilpotent and $N_{t-1}(R)=\gamma_{\infty}(N_{t-2}(R))=1$, a contradiction. Thus $R$ is nilpotent and $R\leq F(G) = Z(G)$. Therefore $R=Z(G)$ and $G$ is quasisimple.

We claim that $G$ contains a power of $\delta^*_k$-commutator $a$ of $2$-order such that $a$ has order $2$ modulo $Z(G)$. Since $G/Z(G)$ is a minimal simple group, by Lemma \ref{invstar}, there exists a $\delta^*_k$-commutator $gZ(G)$ of order $2$. Then there exists a $\delta^*_k$-commutator $y$ in $G$ such that $gZ(G)=yZ(G)$. Hence, we can choose $a$ to be a power of $y$ such that $a$ is a $2$-element and $a$ has order 2 modulo $Z(G)$.

Let $a$ be an element of $G$ with the above properties. Since $G/Z(G)$ is a nonabelian simple group, it follows from \cite[Theorem 3.8.2]{go} that there exists an element $t\in G$ such that the order of $[a,t]$ is odd. On the one hand, since $1=[a^2,t]=[a,t]^a[a,t]$, $a$ inverts $[a,t]$. On the other hand, by Lemma \ref{commfond}, $a$ centralizes $[a,t]$. This is a contradiction.
\end{proof}

Finally, the following result is an immediate consequence of Theorem~B.

\begin {corollary}
Let $G$ be a finite group and let $k \geq 1$. Assume that $|ab|=|a| |b|$ whenever $a$ and $b$ are powers of $\delta^*_k$-commutators such that $(|a|,|b|)=1$. Then $\delta^*_k(G)$ is nilpotent and the Fitting height of $G$ is at most $k+1$.
\end {corollary}

\section{Acknowledgment}
 The authors wish to thank professor Pavel Shumyatsky for interesting discussions. The work of the first author was partially supported by FAPDF/Brazil, while the second author was supported by the ``National Group for Algebraic and Geometric Structures, and their Applications" (GNSAGA - INdAM). Moreover, this study was carried out during the second author's visit to the University of Brasilia. He wishes to thank the Department of Mathematics for the excellent hospitality.