\section{Introduction}\label{Sect:1}

In the last $25$ years, the \acrfull*{RV} technique has been highly successful in the search of planets orbiting other stars. Technical advancements have enabled the discovery of a wide range of exoplanetary types, in stellar systems with architectures far different from our own Solar System. However, the strive to find smaller, Earth-like planets is currently impeded by stellar variability drowning out potential planetary RV signals. The precision required to detect an Earth twin around a solar-like star is of the order of ${\sim}\SI{0.1}{\meter/\second}$, about one order of magnitude smaller than present-day attainability using the best activity mitigation techniques \citep[e.g.][]{Zhao+22}.

By cross-correlating an observed stellar spectrum with a template of the laboratory wavelengths and expected strengths of absorption lines (depending on the spectral type of the star), RVs have traditionally been derived using the entire spectrum at once. Previous studies \citep{Dumusque18,Cretignier+20a} show that a careful selection of \acrfull*{LBL} \acrshortpl*{RV} can help mitigate the effects of stellar activity. \cite{Cretignier+20a} showed that the RV of some lines are uncorrelated with respect to the mean RV of all available lines, and thus seem to be less affected by stellar activity.

The effects of stellar variability on RVs manifest themselves in many different ways. In the structure of low-mass stars, convection in the envelope leads to stellar oscillations and granulation patterns on the stellar surface. Stellar oscillations are due to the propagation of acoustic waves in the convective envelope of late-type stars, which locally modify the stellar surface and can be measured using RVs \citep[e.g.][]{Kjeldsen&Bedding95,Arentoft+08}. The granulation patterns, consisting of ascending hot gas cells and descending cool ones, form blueshifted granules and redshifted intergranulation lanes. Considering the spectral line contributions across a stellar disk, a component formed at the location of a granule would be slightly blueshifted and deeper (due to the higher intensity of granules) than if observed elsewhere \citep{Nordlund+09}. Due to a larger coverage in area and a higher brightness, the granules lead to an observed net \acrfull*{CB} of disk-integrated lines, as well as asymmetries induced by the weaker redshifted contributions. These shifts and asymmetries, known as the third signature of granulation \citep{Gray09}, differ for different spectral lines, with weak lines being affected the most \citep{Reiners+16,Meunier+17a}.

However, stellar oscillations and the variation of the granulation pattern induce RV perturbations on timescales shorter than a few days. These phenomena are therefore not responsible for the dominant RV perturbations which occur on timescales of several days and longer. In fact, both surface oscillations and granulation have been shown to be relatively easy to mitigate with the choice of an appropriate observation strategy \citep{Dumusque+11,Meunier+17b}. This is done by choosing a longer exposure time and observing the target star multiple times per night with sufficient temporal gaps. On timescales longer than a few days, the dominant RV signal is induced by the impact of active regions, which are stellar surface features caused by concentrated magnetic fields \citep{Schrijver&Zwaan:2000}. These active regions can either be spots, which have a lower contrast compared to the inactive stellar surface, or faculae, which have a slightly higher contrast.

Active regions influence the RV precision in mainly two ways \citep{Bauer+18}. The first effect arises because spots and faculae break the flux balance between the two hemispheres which have Doppler velocities of opposite signs due to rotation. This introduces quasi-sinusoidal variations with a periodicity of half the stellar rotation, imitating the Rossiter-McLaughlin effect caused by a transiting planet, but somewhat dissimilar because of spherical projection and varying filling factor at any given time as spots and faculae grow and decay while moving in and out of sight. The second effect arises when the same magnetically active regions inhibit underlying convective motion, which introduces a modulation with a period of the stellar rotation. \cite{Meunier+10} showed that the suppression of convective blueshift, which unlike the flux effect cannot balance itself out when the approaching and receding hemispheres are equally covered, is the dominating contribution in solar simulations, reaching peak-to-peak variations of ${\sim}\SI{10}{\meter/\second}$.

This paper further investigates the selection of spectral lines used to derive RVs, applying a diagnostic tool obtained through spectral synthesis. We implement a novel approach in which line segments formed at different temperature ranges are believed to be more or less susceptible to stellar activity compared to their LBL averages. To test this, we explore whether RV perturbations are stronger in spectral intervals formed deeper in the photosphere, where radial flow velocities and temperature contrasts are larger.

In Sect.~\ref{Sect:2}, we describe the data that we use in the following analysis. In Sect.~\ref{Sect:3}, we describe our method of obtaining the average formation temperature of each sampled point of a stellar spectrum, and then explain how to measure the RV for different temperature ranges. We apply our method to the Sun and $\alpha\,\mbox{Cen}\,\mbox{B}${} in Sect.~\ref{Sect:4} and demonstrate that the stellar activity signal changes as a function of temperature in the stellar atmosphere. We finally discuss the obtained results in Sect.~\ref{Sect:5}.

\section{Data}\label{Sect:2}

\subsection{Observed spectra}\label{Sect:2.1}

In this paper, we analyse spectra of the Sun and $\alpha\,\mbox{Cen}\,\mbox{B}${} which have been obtained with HARPS-N and HARPS, respectively. Both instruments have very similar resolution, $R\,{=}\,115,000$. More precisely, we analyse the echelle-order merged spectra produced by the HARPS and HARPS-N data reduction software. Those spectra are interpolated on a common wavelength grid to facilitate the comparison of all spectra.

Since we are interested in the stellar activity signal that induces RV variations on the rotational period timescale, i.e. dozens of days for Sun-like stars, we daily bin the available spectra to increase the \acrfull*{S/N}, which will be beneficial to derive more precise information at the spectral line level. The median S/N values for the Sun and $\alpha\,\mbox{Cen}\,\mbox{B}${} daily-binned spectra are $890$ and $590$, respectively. We then correct for the color of those daily-binned spectra using \texttt{RASSINE} \citep{Cretignier+20b}, to get the same continuum level for all the spectra, and post-process them using \texttt{YARARA} \citep{Cretignier+21} to remove telluric contamination and known instrumental systematics. To reach the highest S/N possible, a master spectrum is built by combining all available observations for each target. We note that the master spectrum is shifted to the stellar rest frames, to be able to easily compare the spectral lines with line lists from the VALD3\footnote{VALD stands for \acrlong*{VALD}.} database \citep{Piskunov+95,Kupka+00,Ryabchikova+15}.

\subsection{Synthesised spectra}\label{Sect:2.2}

To obtain the average formation height in the stellar atmosphere for each sampled point of a spectrum, we take advantage of spectral synthesis to evaluate the formation contribution function. We perform \acrfull*{1D} \acrfull*{LTE} spectral synthesis using \texttt{PySME}\footnote{Available at \url{https://github.com/AWehrhahn/SME}.}, a Python implementation of \acrlong*{SME} \citep[SME;][]{Valenti&Piskunov96,Piskunov&Valenti17}. For modeling the spectra of the Sun and $\alpha\,\mbox{Cen}\,\mbox{B}${}, we use \texttt{MARCS} model atmospheres \citep{Gustafsson+08} with different stellar parameters and line lists from VALD3 where the line detection threshold is set to $\SI{1}{\percent}$ in accordance with the chosen parameters.

The synthesized spectra are computed on the same wavelength grid as the master spectra of the targets (see Sect.~\ref{Sect:2.1}). For each target, two synthesized spectra will be produced. The first will be generated at the instrumental resolution, and will include micro- and macroturbulences, and rotational broadening using the projected rotational velocity. The goal of this synthesized spectrum is to model as closely as possible the observed spectra. Another synthesis will be generated without turbulence and rotational broadening at the highest resolution possible for the observed wavelength grid. This synthesis will be used to detect line blends which are unresolved at the HARPS and HARPS-N resolution.

\begin{figure*}[t!]
	\includegraphics[width=\textwidth]{Figures/Fig_01.pdf}
	\caption{\textit{Left}: Normalized cumulative contribution function as a function of temperature at $5$ points across the spectral line shown in the right panel. The dashed line indicates the location of $T_{1/2}$, i.e. where the contribution function is $\SI{50}{\percent}$. \textit{Right}: Line profile of one of the solar lines. The $5$ points for which their contribution function is computed are indicated by big dots color-coded with the $T_{1/2}$ values.}
	\label{Fig:1}
\end{figure*}

\begin{table}[h!]
    \caption{Stellar parameters adopted for the spectral synthesis. The parameters are the effective temperature, $T_{\mathrm{eff}}$, surface gravity, $\log{g}$, metalicity, $[\mathrm{M}/\mathrm{H}]$, micro- and macroturbulences, $v_{\mathrm{mic}}$ and $v_{\mathrm{mac}}$, and projected rotational velocity, $v\sin{i}$.}
    \begin{threeparttable}
	\begin{tabular*}{\linewidth}{c @{\extracolsep{\fill}} c @{\extracolsep{\fill}} c}
		\toprule
		\midrule
		\textbf{Parameter}        & \textbf{Sun}  & \textbf{$\alpha\,\mbox{Cen}\,\mbox{B}${}} \\
		\midrule
		$T_{\mathrm{eff}}$ [K]    & 5770\tnote{a} & 5189\tnote{c}     \\
		$\log{g}$                 & 4.00\tnote{a} & 4.30\tnote{c}     \\
		$[\mathrm{M}/\mathrm{H}]$ & 0.00\tnote{a} & 0.22\tnote{c}     \\
		$v_{\mathrm{mic}}$ [km/s] & 0.85\tnote{b} & 0.95\tnote{c}     \\
		$v_{\mathrm{mac}}$ [km/s] & 3.98\tnote{b} & 4.87\tnote{b}     \\
		$v\sin{i}$ [km/s]         & 1.63\tnote{b} & 1.00\tnote{b}     \\
		\bottomrule
	\end{tabular*}
	\begin{tablenotes}
	\item[a] \footnotesize{\texttt{PySME} default.}
	\item[b] \footnotesize{\cite{Valenti&Fischer05}.}
	\item[c] \footnotesize{\cite{Morel18}.}
	\end{tablenotes}
	\end{threeparttable}
	\label{Tab:01}
\end{table}

\section{Methods}\label{Sect:3}

In this section, we describe how to compute, for each spectral point in the master spectrum, the stellar temperature below which $\SI{50}{\percent}$ of the emergent flux originates (see Sect.~\ref{Sect:3.1}), what we will hereafter call the average formation temperature. We thereafter investigate if the stellar activity signal depends on this average formation temperature by dividing each spectral line into a few average formation temperature bins, and measuring the RV for each of those bins. Finally, to reach the meter-per-second precision level in RV to be able to observe activity signals in the Sun and $\alpha\,\mbox{Cen}\,\mbox{B}${}, we combine the information of each average formation temperature bin obtained for all selected lines.

\subsection{Computing the average formation temperature}\label{Sect:3.1}

\begin{figure*}[t!]
	\includegraphics[width=\textwidth]{Figures/Fig_02.pdf}
	\caption{\textit{Top}: HARPS-N solar spectrum in units of normalized flux. \textit{Bottom}: Same spectrum as above converted to units of average formation temperature.}
	\label{Fig:2}
\end{figure*}

\texttt{PySME} automatically calculates the temperature profile, $T(\tau_0)$, on a grid of optical depths, $\tau$, where the zero subscript signifies a reference wavelength of $\SI{5000}{\angstrom}$. In order to evaluate other wavelengths, $\lambda$, the corresponding optical depth, $\tau_\lambda$, is given by \citep{Gray:2005},
\begin{equation}\label{Eq:1}
	\mathrm{d}\tau_\lambda = \frac{\kappa_\lambda}{\kappa_0}\mathrm{d}\tau_0 \, ,
\end{equation}
where $\kappa$ is the sum of the line and continuum opacities. The cumulative contribution function, $C(\tau_\lambda)$, which provides the emergent flux at a certain depth, is proportional to the source function, $S(\tau_\lambda)$, according to the following integral,
\begin{equation}\label{Eq:2}
	C(\tau_\lambda) \propto \int_{0}^{\tau_\lambda} S(\tau_\lambda')\mathrm{e}^{-\tau_\lambda'} \ \mathrm{d}\tau_\lambda' \, .
\end{equation}
Evaluating the source function, which is also provided by \texttt{PySME}, for a discrete grid in optical depths, the integral in Eq.~\ref{Eq:2} becomes a summation, and the contribution function can be normalized by dividing with its last point. Finally, by having both the temperature profile and contribution function computed on the same grid of optical depths for any given wavelength, a change of variable can be made to obtain the contribution function as a function of temperature. In essence, since we are implementing a 1D plane-parallel model atmosphere, both optical depth and temperature are equivalent variables in the sense that they are both monotonically increasing with geometrical depth (although at different rates). However, we chose to convert to temperature as a more intuitive metric which can also be easily related to the commonly known stellar effective temperatures. We refrain from attempting to convert to geometrical depth, as this involves the inclusion of wavelength-dependent stellar radii and their uncertainties.

As a diagnostic average formation, we extract the temperature at which the cumulative contribution function is equal to $\SI{50}{\percent}$, denoted $T_{1/2}$, which is obtained through a linear interpolation. Fig.~\ref{Fig:1} shows $T_{1/2}$ computed at $5$ different points across one of the solar lines, demonstrating that points close to the continuum are formed in hotter (deeper) parts of the photosphere, whereas the line core is formed in cooler (shallower) regions. Although Fig.~\ref{Fig:1} only demonstrates the derived $T_{1/2}$ value for $5$ line points across a single spectral line, the calculations are in fact repeated for every single observed wavelength point of the HARPS and HARPS-N spectra. The bottom panel of Fig.~\ref{Fig:2} shows the entire HARPS-N solar spectrum converted into units of $T_{1/2}$ rather than flux. It can clearly be seen that the formation temperature roughly traces the shapes of the spectral lines locally. The global decrease at longer wavelengths is due to the increase of the continuum opacity, which in the optical wavelength range is dominated by the bound-free absorption of the negative hydrogen ion \citep{Gray:2005}.

\subsection{Line selection}\label{Sect:3.2}

\begin{figure*}[t!]
	\includegraphics[width=\textwidth]{Figures/Fig_03.pdf}
	\caption{Observed (black) and synthesized (green and blue) spectra shown at the position of three exemplary spectral lines. The synthesis including micro- and macroturbulences, rotational broadening and generated at the instrumental resolution of HARPS and HARPS-N is shown in green. The synthesis without turbulence and rotational broadening, and at the highest resolution for the observational wavelength grid is shown in blue. \textit{Left}: Line considered unblended and included in our study. \textit{Middle}: Line rejected due to a blend seen in the unbroadened synthesis. \textit{Right}: Line rejected due to large discrepancy between the observation and the broadened synthesis.}
	\label{Fig:3}
\end{figure*}

\begin{figure*}[t!]
	\includegraphics[width=\textwidth]{Figures/Fig_04.pdf}
	\caption{Solar convective blueshift as a function of $T_{1/2}$ at the line cores. Only lines which could be cross-matched with the VALD3 line list are displayed. The markers of lines from the final selection are color-coded with the line depth shown in the colorbar, the other lines are shown in gray. The gravitational redshift of $\mathrm{RV}_{\mathrm{grav}}{=}GM/Rc$, where $M$ and $R$ are the solar mass and radius, has been subtracted.}
	\label{Fig:4}
\end{figure*}

\begin{figure*}[t!]
	\includegraphics[width=\textwidth]{Figures/Fig_05.pdf}
	\caption{$\SI{15}{\angstrom}$ spectral window color-coded according to the average formation temperature. The shaded areas indicate whether a spectral line is selected (beige) or rejected (grey). The labels above some spectral lines specify the element and ionization of identified lines.}
	\label{Fig:5}
\end{figure*}

In order to properly asses the average formation temperatures for each spectral point, we need to consider lines that are well modeled by the spectral synthesis. To identify spectral lines that are unblended and symmetric, we initially used the morphological parametrisation described in Appendix C of \cite{Cretignier+20a}. Starting from this line selection, we then consider our two syntheses. For each spectral line, we check if the high-resolution profile contains multiple local minima, in which case we reject the line due to blending. Additionally, we also measure the line-depth weighed $\chi^2$ between the instrumental-resolution profile and the observed profile, and reject all the lines for which the $\chi^2$ is more than $4\sigma$ larger than the median value, which can occur due to faulty atomic parameters in the VALD3 database. Fig.~\ref{Fig:3} demonstrates the comparison between observed and synthesized solar spectra for an accepted line, a line rejected because of hidden blends and a line rejected because of large discrepancy.

For each spectral line in the master spectrum, the CB is computed by fitting a second-order polynomial to the center-most $7$ points. The CB is derived from the difference between the fitted center and the laboratory wavelengths found from a cross-match with the VALD3 line list.

Fig.~\ref{Fig:4} shows the solar CB versus $T_{1/2}$ interpolated at the line cores, with errors estimated from the uncertainties of the polynomial coefficients. The figure showcases how our selection of unblended lines clearly accentuates the dependency of CB on certain line properties, as well as how the relation is transformed from non-linear to linear when plotted against formation temperature as opposed to the more common line depth \citep[see e.g. Fig. 3 in][]{Reiners+16}. The CB of $\alpha\,\mbox{Cen}\,\mbox{B}${}, shown in Fig.~\ref{Fig:B1}, follows the same relation, although with a smaller slope than for the Sun which reflects the expected behaviour of more extensive convective envelopes and thus smaller convective velocity gradients in stars with a later spectral type \citep[e.g.][]{Meunier+17b,Liebing+21}. The CB of $\alpha\,\mbox{Cen}\,\mbox{B}${} does not, however, reach near-zero values at low temperatures, which is likely due to the difficulty of constraining the gravitational redshift and orbital motion of the two main components in the Alpha Centauri system \citep[e.g.][]{Pourbaix+02}.

Although the linearity with temperature appears clear when only considering lines for which the estimated core formation temperature is assumed to be reliable, we are alike previous studies \citep[e.g.][]{Gray09,Reiners+16} unable to provide a physical explanation for its origin. To verify the interpretation of the trend, namely that the convective velocity increases linearly with the stellar temperature profile, one could suggestively compare with more complex, magneto-hydrodynamical (MHD) stellar models which take convection into consideration. Such a model comparison is outside the scope of this study, and we note that the following analysis is ultimately unaffected by the CB since our RV time series are derived with respect to a master spectrum, for which the absolute shift is of lesser importance. The CB, which provides the shift relative to laboratory measurements, is simply presented here to validate our line selection and derivation of average formation temperature. By showing its inverse proportionality with respect to the core formation temperature, which on average increases with decreasing line depth, we are able to demonstrate consistency with the literature.

\subsection{RV per temperature bin for individual spectral lines}\label{Sect:3.3}

Once $T_{1/2}$ is calculated for all points of the master spectrum (see Sect.~\ref{Sect:3.1}), our goal is to measure the RV for different bins in average formation temperature as a function of time, by measuring it for each individual spectrum from which the master is built. The RV is derived by template matching as described in \citet{Dumusque18} and \citet{Cretignier+20a}, using the master spectrum as the template. We construct the temperature bins by heuristically dividing the entire range of $T_{1/2}$ values (for the entire spectrum) into $N$ intervals of equal length, where $N$ is a positive integer. Then, for each line profile we measure the RVs of separate segments, by identifying which points are formed within the unique temperature bins. We remark that each line will at most be divided into $N$ segments, however, the majority of lines will only have points in a subset of bins. Generally, the core will constitute the line segment formed at the coolest temperature bin, whereas hotter temperature bins will have nearly symmetrical contributions from the left and right wings. As an example, if we assume the line in Fig.~\ref{Fig:1} is divided into $3$ segments, then points near dot $3$ will be associated to a low temperature bin, points near dots $2$ and $4$ will be associated to a medium temperature bin, and points nears dots $1$ and $5$ will be associated to a high temperature bin. Line segments with less than three points are discarded, to have a good estimate of the flux derivative needed to compute the RVs \citep[see][]{Bouchy+01}. The end product is an RV matrix with dimensions $N_{\mathrm{spec}}{\times}N_{\mathrm{line}}{\times}N_{\mathrm{temp}}$, where $N_{\mathrm{spec}}$ is the number of spectra in the time series, $N_{\mathrm{line}}$ is the number of spectral lines, and $N_{\mathrm{temp}}$ is the selected number of temperature bins, for which some of the entries are empty due to a lack of line points. A matrix of equal size is created for the corresponding RV uncertainties. Averaging the RV measurement of all lines for a given temperature bin (see Sect.~\ref{Sect:3.4}) allows us to reach high RV precision, and therefore to investigate if the stellar activity signal has a dependence on average formation temperature.

\begin{figure*}[t!]
	\includegraphics[width=\textwidth]{Figures/Fig_06.pdf}
	\caption{\textit{Left}: Sun RV time series for $1$ to $4$ temperature bins. Each point is the weighted average of all spectral lines with RV values for a given bin. The shaded intervals represent the $16\textsuperscript{th}$ to $84\textsuperscript{th}$ percentile of all RVs at any given time point. Times are given in \acrfull*{BJD}. \textit{Right}: GLS periodograms of the time series shown in the left panels. The legend displays the temperature bin symbol (see main text and Table~\ref{Tab:02} for definition). The dashed line indicates the $\SI{1}{\percent}$ FAP level.}
	\label{Fig:6}
\end{figure*}

\begin{table*}[t!]
    \caption{Specifications of the temperature-binned RV time series in Fig.~\ref{Fig:6}. For each configuration, the columns specify the number of temperature bins, the bin symbols as described in the text, the intervals in $T_{1/2}$, the total number of lines considered for all bins, the fractions of lines considered for each bin individually, and the RV RMS for the chosen time interval.}
	\begin{tabular*}{\textwidth}{c @{\extracolsep{\fill}} c @{\extracolsep{\fill}} c @{\extracolsep{\fill}} c @{\extracolsep{\fill}} c @{\extracolsep{\fill}} c}
		\toprule
		\midrule
		\textbf{Nr. of bins} & \textbf{Symbol} & \textbf{$T_{1/2}$ [K]} & \textbf{Nr. of lines} & \textbf{Fraction of lines [\%]} & \textbf{RMS [m/s]}\\
		\midrule
		\multirow{1}{*}{1} & $\Delta T_{1/2}(1,1)$ & 4095-5457 & \multirow{1}{*}{773} & 100 & 1.25\\
		\midrule
		\multirow{2}{*}{2} & $\Delta T_{1/2}(2,1)$ & 4095-4776 & \multirow{2}{*}{771} &  55 & 1.17\\
		                   & $\Delta T_{1/2}(2,2)$ & 4776-5457 &                      & 100 & 1.32\\
		\midrule
		\multirow{3}{*}{3} & $\Delta T_{1/2}(3,1)$ & 4095-4549 & \multirow{3}{*}{781} &  16 & 1.69\\
		                   & $\Delta T_{1/2}(3,2)$ & 4549-5003 &                      &  86 & 1.22\\
		                   & $\Delta T_{1/2}(3,3)$ & 5003-5457 &                      &  67 & 1.39\\
		\midrule
		\multirow{4}{*}{4} & $\Delta T_{1/2}(4,1)$ & 4095-4435 & \multirow{4}{*}{780} &   8 & 2.33\\
		                   & $\Delta T_{1/2}(4,2)$ & 4435-4776 &                      &  54 & 1.13\\
		                   & $\Delta T_{1/2}(4,3)$ & 4776-5116 &                      &  79 & 1.32\\
		                   & $\Delta T_{1/2}(4,4)$ & 5116-5457 &                      &  37 & 1.41\\
		\bottomrule
	\end{tabular*}
	\label{Tab:02}
\end{table*}

Identifying line segments on the master spectrum and then using those for each individual spectrum, could raise concerns about different spectral sampling and Doppler shifts between the spectra. The first problem is not an issue as the master spectrum and the individual spectra are sampled on the same wavelength grid. Regarding the shifts, although the master spectrum and the individual spectra are all in the same restframe, stellar activity will induce at most, in the case of our target stars, peak-to-peak Doppler shifts of a few dozens of meters-per-seconds. Such shifts represent less than a tenth of a pixel (${\sim}\SI{820}{\meter/\second}$ for HARPS and HARPS-N), and will thus not inflict with the identification of line segments.

Although line profiles in units of average formation temperature closely resemble their flux counterparts, the idea behind using temperature bins rather than flux bins is that the two quantities are inherently different both locally for single lines and globally for the entire spectrum. If one were to compare two lines of equal line depths yet unequal atomic properties (such as elemental species, ionization stage, excitation energies, among others), they would not necessarily be formed at the same temperature, which depends on the line opacity and hence on these individual properties. This can be seen in Fig.~\ref{Fig:5}, where the flux profiles of a few lines have been color-coded with their $T_{1/2}$ values. The cores of shallow lines are not colored the same as the wings of deeper lines at equal flux levels. Likewise, as was shown in Fig.~\ref{Fig:2}, when considering the entire spectrum, there are also wavelength-dependent variations of the continuum opacity. Both of these effects motive the usage of formation temperature as a preferred variable for probing photospheric height, and the types of activity occurring at various heights.

\subsection{RV per temperature bin for spectral line average}\label{Sect:3.4}

Before averaging line segments from the same formation temperature bin together to increase the RV precision, we curate the selection of line segments included in the averaging. For each RV time series per line and per temperature bin, we remove outliers in the RVs and their uncertainties through an iterative $4\sigma$ clipping (the clipping is only upper-ended for the uncertainties, i.e. we only remove points for which the uncertainty is $4$ standard deviations larger than the median errorbar). RV time series which contain more than $\SI{5}{\percent}$ outliers are entirely rejected. In addition, time series with a $z$-score, defined as the RV \acrfull*{RMS} divided by the median uncertainty, smaller than $2$ are also excluded before averaging to rejected time series with no significant information. We note that for the $4$ bin configuration, the $\SI{5}{\percent}$ outlier criterion rejects an average $\SI{5.5}{\percent}$ of RV time series, while the $z$-score criterion rejects $\SI{9.9}{\percent}$ of them.

The average RV time series per temperature bin is then obtained by averaging the RVs of all remaining spectral line segments, weighting by their inverse variance at each time point. This results in an RV matrix with dimensions $N_{\mathrm{spec}}{\times}N_{\mathrm{temp}}$.

\section{Results}\label{Sect:4}

\subsection{Short-term RVs of the Sun}\label{Sect:4.1}

Short-term RV measurements of the Sun are performed on $105$ daily-binned spectra taken between Jun. 20 and Oct. 19, 2016, in order to cover a few solar rotational periods during an active solar phase. Fig.~\ref{Fig:6} shows the inverse-variance weighted RV time series of all considered lines for when the formation temperature is increasingly divided into $1$ to $4$ bins of equal length (see Sect.~\ref{Sect:3.3}), denoted $\Delta T_{1/2}(N,n)$, where $N$ is the number of bins and $1 {\leq} n {\leq} N$ is the bin index. The figure also shows the corresponding normalized, \acrlong*{GLS} \citep[GLS;][]{Zechmeister&Kurster09} periodograms for each time series with the $\SI{1}{\percent}$ \acrfull*{FAP} indicated. The uppermost left panel of Fig.~\ref{Fig:6}, for the configuration of just $1$ temperature bin (equivalent to no line segmentation based on formation temperature), is thus showing the same RV time series as other traditional methods of measurements would produce, i.e. by measuring LBL RVs and averaging over all lines, or by cross-correlating with a spectral type mask. The periodogram in the uppermost right panel is also showing two expected peaks at the solar synodic rotation period of ${\sim}\SI{27}{d}$ and half of it (the first harmonic). As the RVs are separated into an increasingly larger number of temperature bins (second, third and fourth rows in Fig.~\ref{Fig:6}), we note the emergence of a signal in the coolest temperature bin which is anti-correlated with respect to the total signal and has a significantly stronger periodogram peak at half the rotation period.

Table~\ref{Tab:02} lists the temperature intervals used for each bin configuration, the number of selected lines (which differs between configurations depending on the sigma-clipping and $z$-score criteria described above), and the RV RMS for each bin. For the configuration with $4$ bins, the RMS verifies that the coolest and hottest bins, ${\Delta}T_{1/2}(4,1)$ and ${\Delta}T_{1/2}(4,4)$, show larger dispersion compared to the intermediate bins. Even though the bin averages are calculated using varying number of lines (due to some bins only existing for a certain subset of lines), the $z$-score criteria ensures that larger uncertainties are only included when the RMS is up-scaled by at least the same amount.

While Fig.~\ref{Fig:6} provides a visualisation of the collective (or rather, average) behaviour of all considered lines over time, it does not inform us about the individual lines which contribute to each bin. Fig.~\ref{Fig:A1} shows a corner plot of various line parameters where each dot is an individual line. The dots are colored according to the temperature bin in which they belong, for the configuration with $4$ bins (same as the bottom panel of Fig.~\ref{Fig:6}). Regarding the distribution of the various line parameters and their co-variances, we state the following remarks:
\begin{itemize}
    \item Central wavelength, $\lambda_{\mathrm{c}}$: All bins except the hottest (red) does not have a dependency on chromaticity. This can be explained by the gradient of the continuum opacity (see Fig.~\ref{Fig:2}) causing the hottest bin to be biased toward shorter wavelengths.
    \item Line depth, $d_{\mathrm{c}}$: All bins include lines of any line depth, except the coolest (purple) bin which primarily consists of the cores of the strongest lines.
    \item Convective blueshift, $\mathrm{RV}_{\mathrm{CB}}$: Increases (i.e. approaches zero due to its negative nature) gradually with decreasing temperature. The line segments from the coolest bin are the least affected by CB.
    \item RV RMS for the time interval in Fig.~\ref{Fig:6}: Similar distributions for all bins. Line segments with the smallest scatter originate from an intermediate (cyan) bin.
    \item Pearson correlation between LBL RV and mean RV (i.e. average total RV without temperature binning, same as the upper panel of Fig.~\ref{Fig:6}), $\mathcal{R}(\mathrm{RV}_i, \langle\mathrm{RV}\rangle)$: Most lines segments are weakly correlated with the average RV, except those from the coolest bin which are evenly spread around zero correlation.
    \item Pearson correlation between LBL RV and \ion{Ca}{II} H\&K $S$-index, $\mathcal{R}(\mathrm{RV}_i, S)$: Most line segments are weakly correlated with the $S$-index, except those from the coolest bin which are mostly weakly anti-correlated.
\end{itemize}

\subsection{Solar simulations with \texttt{SOAP-GPU}}\label{Sect:4.2}

\begin{figure}[t!]
	\includegraphics[width=\linewidth]{Figures/Fig_07.pdf}
	\caption{\textit{Top panels}: Flattened intensitygram (left) and magnetogram (right) of the Sun taken by SDO/HMI on Jul. 17, 2016 at UT12:00:00. \textit{Bottom left panel}: Pixels identified as spots (blue) or faculae (red). \textit{Bottom right panel}: Circles representing the same area for spots (blue) and large faculae (red) as the panel to the left, used as \texttt{SOAP-GPU} input.}
	\label{Fig:7}
\end{figure}

To enable the comparison between the temperature-binned RVs and the impact of different active regions, we implement solar simulations for the same time interval as in Sect.~\ref{Sect:4.1} with the \texttt{SOAP-GPU} code \citep{Zhao&DumusquePrep}, a GPU-based successor to the \texttt{SOAP 2.0} code \citep[\acrlong*{SOAP},][]{Dumusque+14}. \texttt{SOAP-GPU} simulates RVs by modeling the observed solar disk as a grid of pixels. Each pixel is assigned a line-of-sight velocity from the combined contributions of rotation and CB. The code then injects each pixel with one of two spatially resolved solar spectra, either a quiet Sun spectrum or a sunspot spectrum, which are Doppler shifted according to the velocity of each pixel. A disk-integrated spectrum is computed by summing all pixels, after scaling the intensity levels with limb darkening and lower/higher contrast for the pixels selected as spot/facula pixels.

From the \acrfull*{HMI} instrument onboard the \acrfull*{SDO}\footnote{Available at \url{https://sdo.gsfc.nasa.gov}.}, we retrieve continuum intensitygrams and magnetograms with a cadence of $\SI{720}{\second}$. We retrieve one frame at noon for every day available between the dates used in Sect.~\ref{Sect:4.1}. We then follow the methodology in \cite{Haywood+16} to extract a threshold map of spot locations from the flattened intensitygrams and facula locations from the unsigned magnetograms. Spot pixels are identified as those where the intensity is below $\SI{89}{\percent}$ of the quiet Sun intensity, and faculae pixels as those where the unsigned magnetic field is larger than $3$ times its standard deviation. With a clustering algorithm, we search for pixel clusters in the threshold maps to find the most prominent regions in each set of images. Isolated pixels are rejected due to possible fluctuations. Large faculae are defined as regions covering at least $\SI{20}{\micro Hem}$, the rest of the selected pixels being associated to the network covering almost the entire disk. For each spot and large facula, the spherical coordinates and radius of the corresponding circle with equal area (see Fig.~\ref{Fig:7}) are computed and given as the sole input to \texttt{SOAP-GPU}.

\begin{figure}[t!]
	\includegraphics[width=\linewidth]{Figures/Fig_08.pdf}
	\caption{\texttt{SOAP-GPU} generated RVs for the same time points as in Fig.~\ref{Fig:6}, separated into active regions (spots, faculae and both) and effects (flux, convection and both). The green curve in the lower right panel is the observed RVs without temperature binning.}
	\label{Fig:8}
\end{figure}

\texttt{SOAP-GPU} enables the total RV signal to be separated into the contributions of different active regions (i.e. spots, faculae or both combined) and effects (i.e. flux, convection or both combined), shown in Fig.~\ref{Fig:8}. We thereafter linearly interpolate the simulated RVs onto the same time sampling as the observations and cross-correlate the temperature-binned RV time series with the simulated contribution of every combination of region/effect by computing the Pearson correlation coefficient for each configuration while allowing for a phase lag going from $-14$ to $14$ days (covering just over one rotational period) in steps of $1$ day (see Fig.~\ref{Fig:9}). We find that the two hottest temperature bins are strongly correlated with the convective effect of faculae, the second-to-coolest bin to not be strongly correlated with any type of contribution, and the coolest bin to be ambiguously anti-correlated with faculae flux with a lag of $-1$ day or faculae convective inhibition with a lag of $5$ days. It is intuitively difficult to argue why stellar activity signals would propagate from one temperature region to the next on a timescale of several days, therefore, it is more likely that the coolest bin in temperature is related to the flux effect induced by faculae.

\begin{figure}[t!]
	\includegraphics[width=\linewidth]{Figures/Fig_09.pdf}
	\caption{Pearson correlation between observed and simulated RVs as a function of time lag. The $4$ temperature bins are represented by the same colors as in Fig.~\ref{Fig:6}.}
	\label{Fig:9}
\end{figure}

\subsection{Long-term RVs of the Sun}\label{Sect:4.3}

\begin{figure*}[t!]
	\includegraphics[width=0.49\textwidth]{Figures/Fig_10_1.pdf}
	\includegraphics[width=0.49\textwidth]{Figures/Fig_10_2.pdf}
	\includegraphics[width=0.49\textwidth]{Figures/Fig_10_3.pdf}
	\includegraphics[width=0.49\textwidth]{Figures/Fig_10_4.pdf}
	\caption{Sliding periodogram of the temperature-binned RV time series of the Sun. The upper panels show the total RV time series with the short-term interval used in Sect.~\ref{Sect:4.1} indicated by the gray, shaded area. The right panels show the GLS periodogram of the total time series, with the solar rotation period and its first harmonic indicated by the dashed lines. The columns of the central colormaps show the GLS periodogram of every $100$-day window, where darker colors indicate stronger power. The window is shifted by steps of 1 day, and gaps correspond to windows with more than $\SI{25}{\percent}$ of the data points missing.}
	\label{Fig:10}
\end{figure*}

Long-term RV measurements of the Sun are performed on $607$ daily-binned spectra taken between Jul. 29, 2015 and Jul. 15, 2018, spanning almost $3$ years near the end of Solar Cycle 24. In order to simultaneously visualize the variability which was observed on the short-term timescale in Sect.~\ref{Sect:4.1} over several similar intervals, and the overall variability when considering the total available time span, we make use of a sliding periodogram, as shown in Fig.~\ref{Fig:10} for the configuration of $4$ temperature bins. The sliding periodogram consists primarily of a color map, where each column represents the GLS periodogram of a fixed time window (here set to $100$ days to equate a few solar rotation periods alike the short-term interval). The window is slid $1$ day for each column, and windows containing less than $\SI{75}{\percent}$ of the expected number of data points (i.e. $75$ points for our window size of $100$ days) are left as blank columns. The color map will thus always have gaps of $\SI{25}{\percent}$ of the window size on the left- and right-most edges, and its vertical extent (the maximum fitted period) will be equal to the window size itself. The upper panel of the periodogram shows the entire RV time series with its conventional periodogram displayed in the right panel.

The sliding periodograms indicate that the $100$-day windows show consistent power distribution over approximately the first half of the covered time period, corresponding to the near end of the most recent solar cycle during which the Sun showed significant levels of activity and active region coverage. On the overall time span, we find the most extreme temperature bins (the coolest and hottest) to be dominated by the long-term magnetic cycle, and apparently oppositely affected similar to their short-term behaviors. The intermediate bins are less affected by low-frequency signals while still showing relatively strong amplitudes at half the solar rotation period.

\subsection{RVs of $\alpha\,\mbox{Cen}\,\mbox{B}${}}\label{Sect:4.4}

RV measurements of $\alpha\,\mbox{Cen}\,\mbox{B}${} are performed on the 2010 observations which have been extensively studied \citep[see e.g.][]{Dumusque18,Cretignier+20a}. The RV curation of $\alpha\,\mbox{Cen}\,\mbox{B}${} follows the same steps as for the Sun, with the results shown in Appendix~\ref{Sect:B}.

We find a similar behaviour for the temperature-binned RVs, with the coolest bin being even more strongly anti-correlated with the $S$-index. The main difference for $\alpha\,\mbox{Cen}\,\mbox{B}${}, however, is the lack of power in the periodograms at around half the rotation period. This could be explained by the stellar inclination of $\alpha\,\mbox{Cen}\,\mbox{B}${}, believed to be $\SI{45}{\degree}^{+9}_{-19}$ \citep{Dumusque14}, weakening the impact of the flux effect due to a decreased rotational velocity component along the line of sight. Such an inclination also causes active regions near the northern pole to remain visible during the entire rotation, therefore amplifying the signal at the rotation period, and not half of it like for the equator-on Sun.

\section{Discussion}\label{Sect:5}

The analysis performed in this paper implies that the activity-induced, temporal RV signal is dependent on the formation temperature inside the stellar photosphere, by demonstrating that the measured time series differs in amplitude and periodicity when computed on line segments formed at various temperature ranges. From solar simulations on short timescales, we find that convection and its inhibition inside active regions due to strong magnetic fields is responsible for most of the observed effect. This dominant contribution is primarily measured in line segments formed at higher temperature ranges, i.e. deeper parts of the photosphere where convective velocities are expected to be larger and their suppression more noticeable. However, for the lower temperature ranges (high in the photosphere), we detect an effect which does not seem to be explained by convection, and is likely due to a change of the physical conditions. Some indicators allude that the signal in this cooler temperature range could be from the flux effect of active regions (due to its stronger periodogram power at half the solar rotation period). Although, this is not satisfactory supported by the simulated RV flux contributions without considering unexplained lags, we do not rule out that the simulated RV curves might not fully capture line shape variations which affect line segments formed at some temperatures more than other. The explanation by some other, unaccounted physics driving the behavior is outside the scope of this paper, however, we point out to a possible change in the contrast of either the active regions, or the quiet Sun granular pattern \citep[e.g.][]{Janssen&Cauzzi06,Cheung+07}.

With average formation temperature as a diagnostic tool, foreseeable improvements could still be made in its derivation for certain spectral portions. For example, the inclusion of non-LTE departure in the syntheses could enable us to more accurately model the strongest lines, whose line cores are the primary constituents of the coolest temperature ranges. By improving their synthesis and thus also their formation height, one could disentangle their primary source of variability more easily and perhaps better understand the physical processes behind.

We also note that in Sect.~\ref{Sect:4.2}, although we compare the temperature-binned RVs to the \texttt{SOAP-GPU} simulations, which outputs values in agreement with observations when measured across all temperature regimes, it is not straightforward to which extent the comparison is valid when considering only segments of lines. Indeed, even if \texttt{SOAP-GPU} was able to fully model the physics, the simulated velocities would always be averaging over all formation temperatures. Our assumption is that the isolated contributions still affect the entire lines the same way as they would on the line segments where they are the most pronounced, although to a diluted extent when combined with the less sensitive segments. Therefore, despite being unable to compare the magnitudes of the effects, the shapes and phases should be comparable.

It is evident from the existing literature, that it is indeed the interplay between magnetic fields and multi-scale convective motion that is currently impeding RV precision to reach the sub-$\SI{}{\meter/\second}$ level. Although a careful selection of lines result in some temperature bins having higher RV precision, it remains to be determined which partition strategy is the most optimal. The trade-off between the number of bins and the number of lines per bin is so far arbitrarily made and should perhaps be physically motivated by stellar structure models and the instrumental noise level rather than the circumstantial range of temperatures from the sample. In a future paper, we intend to explore this parameter space further by, among other aspects, investigating for which temperature interval a minimum dispersion is attained, and at which break-off temperature the strongest periodogram peak transitions between the fundamental and harmonic periodicities. We will also consider developing a mitigation technique in which a linear combination or gradient of the total RV measured in various temperature ranges can be used to distinguishing the activity signals from a planetary-induced, pure Doppler signal, which should be of equal magnitude and invariant between temperature bins.

\section{Conclusions}\label{Sect:6}

In this paper, we build further upon the technique of measuring radial velocities of individual spectral lines. By employing LTE spectral synthesis, we are able to estimate the average formation temperature of spectral line points as sampled by the HARPS-N and HARPS high-resolution spectrographs, for the Sun and $\alpha\,\mbox{Cen}\,\mbox{B}${} respectively.

With a careful selection of unblended, symmetric lines and with average formation temperate as a diagnostic tool, we investigate how various activity-related properties depend on the formation temperature. Our primary reasoning being that formation temperature maps more consistently with photospheric depth than line depth, and should be a more suitable variable for probing different layers. We begin by showing that the absolute convective blueshift of line cores becomes a linear function of formation temperature, instead of a non-linear function of line depth. This is observed for both the studied target stars.

We are furthermore able to demonstrate that spectral line segments formed at different temperatures in the photosphere can be used to diagnose the impact of stellar activity on measured RVs, at both rotational (see Figs.~\ref{Fig:6} and \ref{Fig:B2}) and magnetic cycle timescales (see Fig.~\ref{Fig:10}). Comparison with RV simulations covering a few rotations during which the Sun experienced heightened levels of active region coverage, indicates that convection and its inhibition due to magnetic fields is the dominant effect and perturbs spectral parts formed deeper in the photosphere (at higher temperatures). For both short and long timescales, the line segments formed at cooler temperature ranges exhibit an inverted RV variation when compared to the hotter temperature ranges. Measuring the RV for different temperatures in the photosphere could therefore be a powerful tool to differentiate between stellar activity signals and a planetary signal.

It remains to be investigated whether different types of stellar activity can be entirely disentangled through partitioning of formation temperature as a proxy for photospheric depth, and if variations in RV with temperature can be used to mitigate the activity patterns and enhance the signal of potential planetary signals.

\text{}\\
\footnotesize
\noindent
\textit{Acknowledgments.} We thank the referee for useful comments which helped improve the clarity of the manuscript. We thank Lucia Kleint for interesting discussions. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement SCORE No. 851555). This work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna. This work has made use of the HARPS-N solar RVs as well as ESO public HARPS data.

\balance

\bibliographystyle{aa}


\section{Introduction}

\IEEEPARstart{F}{or} highly correlated input signals, the normalized
subband adaptive filtering (NSAF)~\cite{lee2004improving} algorithm
provides faster convergence than the normalized least mean square
(NLMS) algorithm and retains comparable complexity. The NSAF
algorithm was proposed based on the multiband structure of subband
filters~\cite{lee2009subband}, which adjusts the fullband filter's
coefficients to remove the aliasing and band edge effects of the
conventional subband structure~\cite{lee2009subband}. However, in
practice the non-Gaussian noise with impulsive samples could
commonly happen such as in echo cancellation, underwater acoustics,
audio processing, and
communications~\cite{nikias1995signal,zimmermann2002analysis}, and
in this scenario, the NSAF performance degrades. To deal with
impulsive noises, several robust subband algorithms based on
different robust criteria were proposed,
see~\cite{kim2017delayless,huang2017combined,hur2016variable,zheng2020robust,yu2016two,yu2019m}
and references therein, and most of them can be unified as the NSAF
update with a specific scaling factor.

Furthermore, it is interesting to improve the adaptive filter
performance by exploiting the system sparsity. For example, the
impulse responses of propagation channels in underwater acoustic and
radio communications are usually
sparse~\cite{radecki2002echo,schreiber1995advanced}, only a few
coefficients of which are non-zero. Aiming at sparse systems,
existing examples are classified into the proportionate type and
sparsity-aware type. The family of proportionate NSAF (PNSAF)
algorithms~\cite{abadi2011family} assigns an individual gain to each
filter coefficient, which has faster convergence than the NSAF
algorithm. Later, robust PNSAF algorithms were also
presented~\cite{zheng2017robust,yu2019m} to deal with impulsive
noises. On the other hand, the family of sparsity-aware algorithms
incorporates the sparsity-aware penalty into the original NSAF's and
PNSAF's cost functions; as a result, sparsity-aware NSAF
(SA-NSAF)~\cite{yu2019sparsity,heydari2021improved} and
sparsity-aware PNSAF~\cite{puhan2019zero} algorithms were developed.
In sparse system identification, these sparsity-aware algorithms can
obtain better convergence and steady-state performance than their
original counterparts.

However, the superiority of sparsity-aware algorithms depends mainly
on the sparsity-penalty parameter, which is often chosen in an
exploratory way thus reducing the practicality of the algorithms.
Besides, they encounter the problem of choosing the step-size, which
controls the tradeoff between convergence rate and steady-state
misadjustment. Hence, adaptation techniques for the sparsity-penalty
and the step-size parameters are necessary. In the literature, they
are rarely discussed simultaneously regardless of the Gaussian noise
or impulsive noise scenarios. In~\cite{ji2020sparsity}, the variable
parameter SA-NSAF (VP-SA-NSAF) algorithm was proposed for the
Gaussian noise, in which these two parameters are jointly adapted
based on a model-driven method, but it requires knowledge of
variances of the subband noises. In~\cite{yu2021sparsity}, by
optimizing the parameters in the sparsity-aware
individual-weighting-factors-based sign subband adaptive filter
(S-IWF-SSAF) algorithm with the robustness in the impulsive noise,
the variable parameters S-IWF-SSAF (VP-S-IWF-SSAF) algorithm was
presented, while it lacks the generality in sparsity-aware subband
algorithms.

In this paper, we propose a unified sparsity-aware robust NSAF
(SA-RNSAF) framework to handle impulsive noises, which can result in
different algorithms by only changing the robustness criterion and
the sparsity-aware penalty. We then devise adaptive schemes for
adjusting the step-size and the sparsity-aware penalty weight, and
develop the alternating optimization of the parameters based
SA-RNSAF (AOP-SA-RNSAF) algorithm, with fast convergence and low
steady-state misadjustment for sparse systems.

\section{Statement of Problem and SA-RNSAF Algorithm}

Consider a system identification problem. The relationship between
the input signal $u(n)$ and desired output signal $d(n)$ at time~$n$
is given by
\begin{align}
\label{001}
d(n) = \bm u^\text{T}(n) \bm w^o+\nu(n),
\end{align}
where the $M\times1$ vector $\bm w^o$ is the impulse response of the
sparse system that we want to identify, $\bm
u(n)=[u(n),u(n-1),...,u(n-M+1)]^\text{T}$ is the $M\times1$ input
vector, and $\nu(n)$ is the additive noise independent of $u(n)$.
For estimating~$\bm w^o$, the SAF with a coefficient vector~$\bm
w(k)$ is used, shown in Fig.~\ref{Fig1} with $N$ subbands, where $k$
denotes the sample index in the decimated domain. The input signal
$u(n)$ and the desired output signal~$d(n)$ are decomposed into
multiple subband signals $u_i(n)$ and $d_i(n)$ via the analysis
filters $\{\bm h_i\}_{i=1}^N$, respectively. For each subband input
signal~$u_i(n)$, the corresponding output of the fullband
filter~$\bm w(k)$ is $y_i(n)$. Then, both $d_i(n)$ and $y_i(n)$ are
critically decimated to yield signals $d_{i,D}(k)$ and $y_{i,D}(k)$,
respectively, with lower sampling rate, namely, $d_{i,D}(k)=d_i(kN)$
and $y_{i,D}(k)=\bm u_i^\text{T}(k)\bm w(k)$, where $\bm
u_i(k)=[u_i(kN),u(kN-1),...,u(kN-M+1)]^\text{T}$. By subtracting
$y_{i,D}(k)$ from $d_{i,D}(k)$, the decimated subband error signals
are obtained:
\begin{align}
e_{i,D}(k) = d_{i,D}(k)-\bm u_i^\text{T}(k)\bm w(k),\;i=1,2,...,N,
\label{002}
\end{align}
which are used to adjust the coefficient vector~$\bm
w(k)$~\footnote{In some applications, we could also eventually need
the output error $e(n)$ in the original time domain. To this end, we
obtain $\bm w(n)$ by copying $\bm w(k)$ for every $N$ input samples,
and then compute the output error by $e(n) = d(n)-\bm u^\text{T}(n)
\bm w(n)$. }.
\begin{figure}[htb]
    \centering
    \includegraphics[scale=0.42] {fig1.eps}

    \vspace{-1em} \caption{Multiband structure of subband adaptive filter. }
    \label{Fig1}
\end{figure}

In practice, the additive noise~$\nu(n)$ can be non-Gaussian
consisting of Gaussian and impulsive components
~\cite{jidf,spa,intadap,mbdf,jio,jiols,jiomimo,sjidf,ccmmwf,tds,mfdf,l1stap,mberdf,jio_lcmv,locsme,smtvb,ccmrls,dce,itic,jiostap,aifir,ccmmimo,vsscmv,bfidd,mbsic,wlmwf,bbprec,okspme,rdrcb,smce,armo,wljio,saap,vfap,saalt,mcg,sintprec,stmfdf,1bitidd,jpais,did,rrmber,memd,jiodf,baplnc,als,vssccm,doaalrd,jidfecho,dcg,rccm,ccmavf,mberrr,damdc,smjio,saabf,arh,lsomp,jrpaalt,smccm,vssccm2,vffccm,sor,aaidd,lrcc,kaesprit,lcdcd,smbeam,ccmjio,wlccm,dlmme,listmtc,smcg},
. Hence, for the identification of a sparse vector~$\bm w^o$ in the
presence of impulsive noise, we define the following minimization
problem:
\begin{equation}
\label{003}
\begin{array}{rcl}
\begin{aligned}
\arg \min_{\bm w(k+1)} &\left[ ||\bm w(k+1) - \bm w(k)||_2^2 +\rho f(\bm w(k+1)) \right],\\
\end{aligned}
\end{array}
\end{equation}
subject to
\begin{subequations} \label{eq:4}
    \begin{align}

e_{p,i}(k) =& \left[ 1 - \mu \phi_i(k) \right] e_{i,D}(k),

    \label{eq:4a}\\
    \phi_i(k)=&\frac{\varphi'(e_{i,D}(k))}{e_{i,D}(k)}, \label{eq:4b}
    \end{align}
\end{subequations}
for subbands $i=1,...,N$, where $e_{p,i}(k)=d_{i,D}(k)-\bm u_i^\text{T}(k)\bm w(k+1)$ denotes the \textit{a posteriori} decimated subband error, $\mu>0$ will be called the step-size in the sequel, and $\phi_i(k)$ is called the scaling factor of the $i$-th subband. In~\eqref{003}, $f(\bm w)$ is a sparsity-aware penalty function and $\rho>0$ is the weight of this penalty term. In~\eqref{eq:4b}, $\varphi'(e) \triangleq \frac{\partial \varphi(e)}{\partial e}$, where $\varphi(e)\geq 0$ is an even function of variable~$e$, defining the robustness to impulsive noise.

By using the Lagrange multiplier method, we obtain the solution of~\eqref{003} subject to~\eqref{eq:4a} as
\begin{equation}
\label{005}
\begin{array}{rcl}
\begin{aligned}
\bm w&(k+1) =\bm w(k) + \mu \sum_{i=1}^{N} \phi_i(k) \frac{e_{i,D}(k)\bm u_i(k)}{||\bm u_i(k)||_2^2} - \\
&\rho \left[ f'(\bm w(k+1)) - \sum_{i=1}^{N} \frac{\bm u_i(k) \bm u_i^\text{T}(k)}{||\bm u_i(k)||_2^2} f'(\bm w(k+1)) \right].
\end{aligned}
\end{array}
\end{equation}
Note that the derivation of~\eqref{005} also uses an approximation
in the SAF, that is~$\bm u_i^\text{T}(k) \bm u_j(k)\approx 0$ for
$i\neq j$~\cite{lee2004improving}. Then, by introducing an
intermediate estimate~$\bm \psi(k)$, we propose to
implement~\eqref{005} in two steps:
\begin{subequations} \label{eq:6}
    \begin{align}

    \bm \psi(k) = &\bm w(k) + \mu \sum_{i=1}^{N} \phi_i(k) \frac{e_{i,D}(k)\bm u_i(k)}{||\bm u_i(k)||_2^2},

    \label{eq:6a}\\
    \bm w(k+1) =& \bm \psi(k) - \rho \bm P(k), \label{eq:6b}
    \end{align}
\end{subequations}
where
\begin{equation}
\label{007}
\begin{array}{rcl}
\begin{aligned}
\bm P(k) = f'(\bm \psi(k)) - \sum_{i=1}^{N} \frac{\bm u_i(k) \bm u_i^\text{T}(k)}{||\bm u_i(k)||_2^2} f'(\bm \psi(k)).
\end{aligned}
\end{array}
\end{equation}
This completes the derivation of the update for the SA-RNSAF
algorithm. In this algorithm, the steps~\eqref{eq:6a}
and~\eqref{eq:6b} have their own roles. The former behaves like the
RNSAF algorithm to obtain a coarse estimate~$\bm \psi(k)$ of~the
sparse vector $\bm w^o$ in impulsive noise. Subsequently, the
step~\eqref{eq:6b} forces the inactive coefficients in $\bm \psi(k)$
to zero, thus obtaining a more accurate sparse estimate~$\bm
w(k+1)$.

It is noteworthy that the parameters~$\mu$ and $\rho$ control the
SA-RNSAF's performance. Specifically, the step-size~$\mu$ controls
the convergence rate and steady-state misadjustment of the
algorithm. Moreover, the SA-RNSAF algorithm can be superior to the
RNSAF algorithm when dealing with sparse systems, but $\rho$ must be
chosen within a theoretically existing range while this range is
unpredictable actually (see Remark~1 below). As such, we will derive
adaptive recursions for adjusting $\mu$ and $\rho$. However, it is
challenging to solve the global optimization problem on $\mu$ and
$\rho$, as~\eqref{eq:6a} and \eqref{eq:6b} depend on each other.
Interestingly, $\mu$ and $\rho$ mainly affect the
steps~\eqref{eq:6a} and \eqref{eq:6b}, respectively, thus we can use
alternating optimization~\cite{hong2016convergence} to solve this
global optimization problem. Accordingly, the adaptations of $\mu$
and $\rho$ will be designed independently according to~\eqref{eq:6a}
and~\eqref{eq:6b}, respectively.

\section{Proposed AOP-SA-RNSAF Algorithm}

By using the band-dependent variable step-size (VSS) $\mu_i(k)$ and
$\rho(k)$ instead of some fixed values, we rearrange~\eqref{eq:6a}
and \eqref{eq:6b} as follows:
\begin{subequations} \label{eq:8}
    \begin{align}

    \bm \psi(k) = &\bm w(k) + \sum_{i=1}^{N} \mu_i(k) \phi_i(k) \frac{e_{i,D}(k)\bm u_i(k)}{||\bm u_i(k)||_2^2},

    \label{eq:8a}\\
    \bm w(k+1) =& \bm \psi(k) - \rho(k) \bm P(k). \label{eq:8b}
    \end{align}
\end{subequations}
\subsection{Adaptation of the step-size}
By subtracting~\eqref{eq:8a} from $\bm w^o$, we obtain
\begin{equation}
\label{009}
\begin{array}{rcl}
\begin{aligned}
\widetilde{\bm \psi}(k) = &\widetilde{\bm w}(k) - \sum_{i=1}^{N} \mu_i(k) \phi_i(k) \frac{e_{i,D}(k)\bm u_i(k)}{||\bm u_i(k)||_2^2},
\end{aligned}
\end{array}
\end{equation}
where $\widetilde{\bm w}(k)=\bm w^o - \bm w(k)$ and $\widetilde{\bm \psi}(k)=\bm w^o - \bm \psi(k)$ define the deviation vectors for the final estimate~$\bm w(k)$ and the intermediate estimate $\bm \psi(k)$ with respect to the true value. By pre-multiplying~$\bm u_i^\text{T}(k)$ on both sides of~\eqref{009} and applying the approximation~$\bm u_i^\text{T}(k) \bm u_j(k)\approx 0$ for $i\neq j$ again, it is established that
\begin{equation}
\label{010}
\begin{array}{rcl}
\begin{aligned}
e_{\varepsilon,i}(k) = \left[ 1 -\mu_i(k) \phi_i(k) \right]  e_{i,D}(k)
\end{aligned}
\end{array}
\end{equation}
for $i=1,2,...,N$, where $e_{\varepsilon,i}(k) = d_{i,D}(k)-\bm u_i^\text{T}(k) \bm \psi(k)$ defines the intermediate \textit{a posteriori} error at the $i$-th subband resulting from the step~\eqref{eq:6a}. By squaring both sides of~\eqref{010} and taking the expectations over all the terms, the following relation is obtained:
\begin{equation}
\label{011}
\begin{array}{rcl}
\begin{aligned}
\text{E}\{e_{\varepsilon,i}^2(k)\} = \left[ 1 -\mu_i(k) \phi_i(k) \right]^2  \text{E}\{e_{i,D}^2(k)\},
\end{aligned}
\end{array}
\end{equation}
where $\text{E}\{\cdot\}$ denotes the mathematical expectation. In~\eqref{011}, a common assumption is used that the step-size $\mu_i(k)$ and the scaling factors $\{\phi_i(k)\}_{i=1}^N$ are deterministic at iteration~$k$ in contrast with the randomness of $e_{i,D}(k)$~\cite{ni2010variable,hur2016variable}.

Motivated by~\cite{ni2010variable}, we wish to compute the subband step-sizes in such a way that $\text{E}\{e_{\varepsilon,i}^2(k)\}=\sigma_{\nu,i}^2$, $i=1,2,...,N$, which means that the powers of the intermediate \textit{a posteriori} subband errors always equal those of the subband noises, where $\sigma_{\nu,i}^2\triangleq\text{E}\{\nu_{i,D}^2(k)\}$ denotes the power of the $i$-th subband noise excluding impulsive interferences. on this requirement, then from~\eqref{011} we can obtain the following equation:
\begin{equation}
\label{012}
\begin{array}{rcl}
\begin{aligned}
\mu_i(k) \phi_i(k) = 1- \sqrt{\frac{\sigma_{\nu,i}^2}{\sigma_{e_{i,D}}^2(k)}},
\end{aligned}
\end{array}
\end{equation}
where $\sigma_{e_{i,D}}^2(k)\triangleq \text{E}\{e_{i,D}^2(k)\}$ indicates the power of $e_{i,D}(k)$ without impulsive noises. For robust adaptive algorithms with the scaling factors, there is a common property~\cite{huang2017combined,hur2016variable} that when impulsive noises happen, the scaling factors $\phi_i(k)$ will become very small (close to zero), thereby preventing the adaptation~\eqref{eq:8a} from the interference caused by impulsive noises. If the impulsive noise is absent, $\phi_i(k)$ will approximately equal one to ensure fast convergence. As such, we can change~\eqref{012} to
\begin{equation}
\label{013}
\begin{array}{rcl}
\begin{aligned}
\mu_i(k) = 1- \sqrt{\frac{\sigma_{\nu,i}^2}{\sigma_{e_{i,D}}^2(k)}}.
\end{aligned}
\end{array}
\end{equation}
To implement \eqref{013}, the statistical quantities $\sigma_{e_{i,D}}^2(k)$ and $\sigma_{\nu,i}^2$ are replaced with their estimates $\hat{\sigma}_{e_{i,D}}^2(k)$ and $\hat{\sigma}_{\nu,i}^2(k)$, respectively. Specifically, $\hat{\sigma}_{e_{i,D}}^2(k)$ is calculated in an exponential window way as
\begin{equation}
\label{014}
\begin{array}{rcl}
\begin{aligned}
\hat{\sigma}_{e_{i,D}}^2(k) = \zeta \hat{\sigma}_{e_{i,D}}^2(k-1) + (1-\zeta)\phi_i^2(k) e_{i,D}^2(k),
\end{aligned}
\end{array}
\end{equation}
where $\zeta$ is a weighting factor often chosen as $\zeta=1-1/(\kappa M)$ with $\kappa \geq1$. Similar to~\cite{ni2010variable}, $\hat{\sigma}_{\nu,i}^2(k)$ is calculated by the following equations:
\begin{subequations} \label{eq:15}
    \begin{align}

    \hat{\sigma}_{u_{i}}^2(k) &= \zeta \hat{\sigma}_{u_{i}}^2(k-1) + (1-\zeta) u_{i}^2(kN),

    \label{eq:15a}\\
    \hat{\bm r}_{ue_{i}}(k) &= \zeta \hat{\bm r}_{ue_{i}}(k-1) + (1-\zeta) \phi_i(k) e_{i,D}(k) \bm u_{i}(k), \label{eq:15b}\\
    \hat{\sigma}_{\nu,i}^2(k) & = \hat{\sigma}_{e_{i,D}}^2(k) - \frac{||\hat{\bm r}_{ue_{i}}(k)||_2^2}{\hat{\sigma}_{u_{i}}^2(k) + \epsilon_1}, \label{eq:15c}
    \end{align}
\end{subequations}
where $\epsilon_1$ is a small positive number (e.g., $10^{-5}$). Note that, we introduce the scaling factor $\phi_i(k)$ in (14) and (15b) for each subband to suppress impulsive noises.

Accordingly, \eqref{013} can be rewritten as
\begin{equation}
\label{016}
\begin{array}{rcl}
\begin{aligned}
\mu_i(k) = 1- \sqrt{\frac{\hat{\sigma}_{\nu,i}^2(k)}{\hat{\sigma}_{e_{i,D}}^2(k) + \epsilon_2}},
\end{aligned}
\end{array}
\end{equation}
where $\epsilon_2$ is also a small positive number. It is stressed that the estimated values of multiple statistical quantities are used in~\eqref{eq:15c}, and thus $\hat{\sigma}_{\nu,i}^2(k)$ could be negative at some iterations. To avoid this, we add the step $\hat{\sigma}_{\nu,i}^2(k)\leftarrow \hat{\sigma}_{\nu,i}^2(k-1)$ after~\eqref{eq:15c}.
\subsection{Adaptation of the sparsity penalty weight}
By subtracting~\eqref{eq:8b} from $\bm w^o$, we obtain
\begin{equation}
\label{017}
\begin{array}{rcl}
\begin{aligned}
\widetilde{\bm w}(k+1) = \widetilde{\bm \psi}(k) + \rho(k) \bm P(k).
\end{aligned}
\end{array}
\end{equation}
By pre-multiplying both sides of~\eqref{017} by their transpose, we obtain
\begin{equation}
\label{018}
\begin{array}{rcl}
\begin{aligned}
||\widetilde{\bm w}&(k+1)||_2^2 = ||\widetilde{\bm \psi}(k)||_2^2 + \triangle(k),
\end{aligned}
\end{array}
\end{equation}
where
 \begin{equation}
 \label{019}
 \begin{array}{rcl}
 \begin{aligned}
 \triangle(k) = 2 \rho(k) \widetilde{\bm \psi}^\text{T}(k)\bm P(k) + \rho^2(k) ||\bm P(k)||_2^2.
 \end{aligned}
 \end{array}
 \end{equation}

Remark~1: \eqref{018} clearly reveals that the proposed SA-RNSAF algorithm will outperform the RNSAF algorithm for identifying sparse systems, if and only if $\triangle(k)<0$ holds\footnote{Following a derivation similar to that in~Appendix~D in \cite{yu2021sparsity}, $\triangle(k)<0$ is likely to be true as long as $\bm w^o$ is sparse.}. It follows that $\rho(k)$ should satisfy the inequality
  \begin{equation}
 \label{020}
 \begin{array}{rcl}
 \begin{aligned}
0<\rho(k) < 2 \frac{ [\bm \psi(k) - \bm w^o]^\text{T} \bm P(k)}{||\bm P(k)||_2^2}.
 \end{aligned}
 \end{array}
 \end{equation}
Moreover, since $\triangle(k)$ is the quadratic function of~$\rho(k)$, there is an optimal $\rho(k)$ such that $\triangle(k)$ arrives at the negative maximum value. Consequently, the optimal $\rho(k)$ is given as
\begin{equation}
\label{021}
\begin{array}{rcl}
\begin{aligned}
\rho_\text{opt}(k) = \frac{ [\bm \psi(k) - \bm w^o]^\text{T} \bm P(k)}{||\bm P(k)||_2^2}.
\end{aligned}
\end{array}
\end{equation}

Although Remark~1 states that the relations~\eqref{020} and~\eqref{021} are existing in sparse systems, they are incalculable due to the fact that the sparse vector~$\bm w^o$ is unknown. To solve this problem, we use the previous estimate $\bm w(k)$ to approximate $\bm w^o$, then~\eqref{021} can be reformulated as
\begin{equation}
\label{022}
\begin{array}{rcl}
\begin{aligned}
\hat{\rho}_\text{opt}(k) = \max \left\lbrace \frac{ [\bm \psi(k) - \bm w(k)]^\text{T} \bm P(k)}{||\bm P(k)||_2^2},0\right\rbrace ,
\end{aligned}
\end{array}
\end{equation}
where $\hat{\rho}_\text{opt}(k)$ is set to zero at $k=0$.

The recursion~\eqref{eq:8} equipped with $\mu_i(k)$ in~\eqref{016} and $\hat{\rho}_\text{opt}(k)$ in~\eqref{022} constitutes the proposed AOP-SA-RNSAF algorithm.

Remark~2: The proposed AOP-SA-RNSAF update generalizes different algorithms, depending on the choice of $\varphi(e)$ in~\eqref{eq:4b} and $f(\bm w)$ in~\eqref{003}. In the literature, several robust criteria against impulsive noises~\cite{zheng2017robust,yu2016two,huang2017combined,hur2016variable,yu2019m} defined by $\varphi(e)$ and sparsity-aware penalties~\cite{yu2021sparsity,heydari2021improved,ji2020sparsity,yu2019sparsity,de2014sparsity} defined by $f(\bm w)$ have been studied, which can be applied in the AOP-SA-RNSAF. Nevertheless, this paper does not consider the effect of different choices of $\varphi(e)$ and/or $f(\bm w)$, which is worth studying in future work. Note that, when setting $\varphi(e)=\frac{1}{2}e^2$, the proposed algorithm is called the alternating optimization of parameters based SA-NSAF (AOP-SA-NSAF) suited for Gaussian noise environments, which is a sparsity-aware variant of the VSS-NSAF algorithm presented in~\cite{ni2010variable}.

Remark~3: By firstly computing the inner product~$\bm u_i^\text{T}(k)f'(\bm \psi(k))$ in~$\bm P(k)$, and then calculating~$\bm P(k)$ only requires $2M$ multiplications, $2M-M/N$ additions, and 1 division. Therefore, the complexity of the proposed AOP-SRNSAF algorithm is still low with $\mathcal{O}(M)$ arithmetic operations per input sample.
\section{Simulation Results}
To evaluate the proposed AOP-SA-RNSAF algorithm, simulations are conducted to identify the acoustic echo paths with $M=512$ taps. The sparsenesses, defined as $\chi(\bm w^o) = \frac{M}{M-\sqrt{M}} \left(1-\frac{||w^o||_1}{\sqrt{M}||w^o||_2} \right)$, of two echo paths are $\chi(\bm w^o_1)=0.9357$  (sparse)~\cite{yu2021sparsity} and $\chi(\bm w^o_2)=0.3663$ (dispersive or non-sparse)~\cite{zheng2017robust}, respectively. The length of the adaptive filters is the same as that of $\bm w^o$. The correlated input signal $u(n)$ is a first-order autoregressive (AR) process with the pole at 0.9, generated by filtering a white Gaussian noise with zero-mean and unit variance. The analysis filters~$\{\bm h_i\}_{i=0}^{N-1}$ for decomposing signals $d(n)$ and $u(n)$ are obtained by cosine-modulated filter banks, where the length of the prototype filter for $N=4$ subbands is 33 to obtain 60 dB stopband attenuation. The high stopband attenuation is to guarantee that adjacent analysis filters have almost no overlap and the cross-correlation of nonadjacent subbands is negligible~\cite{lee2009subband}. The normalized mean square deviation (NMSD), defined as $\text{E}\{||\bm w(n)-\bm w^o||_2^2/||\bm w^o||_2^2\}$, is the performance measure. All the results are the average over 50 independent runs.

For the AOP-SA-RNSAF algorithm, we use the modified Huber~(MH) function for $\varphi(e)$ and the log-penalty for $f(\bm \psi(k))$. The MH function is formulated as $\varphi(e_{i,D}(k))=e_{i,D}^2(k)/2$ if $|e_{i,D}(k)| < \xi_i$ and $\varphi(e_{i,D}(k))=0$  if $|e_{i,D}(k)| \geq \xi_i$~\cite{yu2019m}, where $\xi_i$ is the threshold. Accordingly, when $|e_{i,D}(k)| \geq \xi_i$ (usually impulsive noises occur), then the scaling factor in~\eqref{eq:4b} is $\phi_i(k)=0$, which makes the adaptation step~\eqref{eq:8a} freeze to suppress impulsive interferences; otherwise, $\phi_i(k)=1$, which retains fast convergence. Note that, the threshold $\xi_i$ for each subband $i$ is set to $\xi_i = 2.576 \hat{\sigma}_{\varepsilon,i}(k)$, where $\hat{\sigma}_{\varepsilon,i}^2(k)$ is the variance of $e_{i,D}(k)$ excluding impulsive samples. $\hat{\sigma}_{\varepsilon,i}^2(k)$ is computed by~$\hat{\sigma}_{\varepsilon,i}^2(k) = \lambda \hat{\sigma}_{\varepsilon,i}^2(k-1) + c_\sigma(1-\lambda) \text{med}(\bm a_{\varepsilon,i})$, where $\lambda\in (0.9, 1)$ is the forgetting factor (but $\lambda=0$ at $k=0$), $\text{med}(\cdot)$ denotes the median operator to remove outliers in the data window $\bm a_{\varepsilon,i}=[e_{i,D}^2(k),e_{i,D}^2(k-1),...,e_{i,D}^2(k-N_w+1)]$ with a length of $N_w$, and $c_\sigma=1.483(1+5/(N_w-1))$ is the correction factor. The log-penalty is given as~$f(\bm \psi_k) = \sum_{m=1}^{M} \ln(1+|\psi_{m,k}|/\theta)$~\cite{yu2021sparsity} which characterizes the sparsity of systems, where $\psi_{m,k}$ is the $m$-th element of $\bm \psi_k$, and the shrinkage factor~$\theta>0$ cuts apart inactive and active entries in~$\bm \psi_k$. Thus, $f'(\bm \psi_{m,k})$ in~\eqref{007} is computed element-wise as~$f'(\psi_{m,k}) = \frac{\text{sgn}(\psi_{m,k})}{\theta + |\psi_{m,k}|}, m=1,...,M$. In our simulations, the additive noise $\nu(n)$ is described by the symmetric $\alpha$-stable process, also called the $\alpha$-stable noise, whose characteristic function is formulated as $\phi(t)=\exp(-\gamma \lvert t \lvert^\alpha)$~\cite{nikias1995signal}. The parameter~$\alpha \in (0,2]$ represents the impulsiveness of the noise that for smaller $\alpha$ leads to stronger impulsive noises, and $\gamma>0$ behaves like the variance of the Gaussian density. In particular, it reduces to the Gaussian noise for the case of~$\alpha=2$. In the following simulations, we set $\gamma=0.02$.

Example~1: the impulsive noise is absent, i.e., $\alpha=2$. The proposed AOP-SA-NSAF algorithm in Remark~2 is compared with the NSAF, VP-SA-NSAF~\cite{ji2020sparsity}, VSS-NSAF, and VSS-PNSAF algorithms in Fig.~\ref{Fig3}, where both VSS-NSAF and VSS-PNSAF are obtained from~\cite{yu2019m} but in the Gaussian noise we reset~$\varphi(e)=\frac{1}{2}e^2$ instead of using the MH function. For a fair evaluation, we select the log-penalty parameter $\theta=0.005$ for all the sparsity-aware algorithms. As can be seen, the VSS-NSAF algorithm obtains fast convergence and low steady-state misadjustment, which overcomes the trade-off problem in the NSAF algorithm. By considering the sparsity of the underlying system, both VP-SA-NSAF and VSS-PNSAF algorithms further improve the convergence rate. As compared to the VSS-PNSAF algorithm, the proposed AOP-SNSAF algorithm shows slower initial convergence, but it achieves higher reduction in the steady-state misadjustment.
\vspace{-0.5cm}
 \begin{figure}[htb]
    \centering
    \includegraphics[scale=0.5] {fig2.eps}
    \vspace{-1em} \caption{NMSD performance of NSAF-type algorithms in the Gaussian noise. The parameters of algorithms are listed as follows: $\eta=0.99$, $\lambda=0.95$, and $\mu_{\max}=1$ for VP-SA-NSAF; $\tau=3$ for VSS-NSAF; $\tau=5$ and $\zeta=0$ for VSS-PNSAF; $\kappa=6$ for AOP-SA-NSAF. }
    \label{Fig3}
 \end{figure}

Example~2: $\alpha=1.6$ displays the presence of impulsive noises. Fig.~\ref{Fig4} depicts the NMSD performance of the NSAF, M-NSAF~\cite{yu2019m}, VSS-M-NSAF~\cite{yu2019m}, VSS-M-PNSAF~\cite{yu2019m}, VP-IWF-SSAF with RA~\cite{yu2021sparsity}, and the proposed AOP-SA-RNSAF algorithms\footnote{Since the variance of the $\alpha$-stable noise is nonexistent, here we do not show the performance of the VP-SA-NSAF algorithm.}. For the M-estimate based algorithms, we choose the common M-estimate parameters $\lambda=0.99$ and $N_w =20$. It is seen that the NSAF algorithm shows poor misadjustment in the $\alpha$-stable noise, and other algorithms exhibit robust convergence. Among these robust algorithms, the proposed AOP-SA-RNSAF algorithm is the best choice for identifying sparse systems, due to the fact that it has lower steady-state misadjustment than the VSS-M-PNSAF and VP-S-IWF-SSAF with RA algorithms, even if it has a slightly slower initial convergence than the VSS-M-PNSAF algorithm.
\begin{figure}[htb]
    \centering
    \includegraphics[scale=0.5] {fig3.eps}
    \hspace{2cm}\caption{NMSD performance of NSAF-type algorithms in the $\alpha$-stable noise. The parameters of algorithms are listed as follows: $\mu_{\min}=10^{-5}$, $\tau=2$, and $\chi=1$ for VP-S-IWF-SSAF; $\tau=3$ for VSS-M-NSAF; $\tau=5$ and $\zeta=0$ for VSS-M-PNSAF; $\kappa=6$ for AOP-SA-RNSAF.}
    \label{Fig4}
\end{figure}

It can be seen in Figs.~\ref{Fig3} and~\ref{Fig4} that, after $\bm w^o$ becomes dispersive at the middle of input samples, the proportionate-type (i.e., VSS-PNSAF, VSS-PNSAF) and sparsity-aware type (i.e., VP-SA-NSAF, AOP-SA-NSAF, AOP-SA-RNSAF) algorithms still show almost the same performance as the competing algorithms (i.e., VSS-NSAF and VSS-M-NSAF) in both Gaussian and $\alpha$-stable noise scenarios. In addition, as $\alpha$ decreases from 2 to 1.6, the steady-state misadjustment of the proposed AOP-SA-RNSAF algorithm increases, but this algorithm is still convergent.
\section{Conclusion}
In this paper, a unified SA-RNSAF framework for algorithms was developed for identifying sparse systems in impulsive noise environments. By replacing directly the specified robustness criterion and sparsity-aware penalty, it can yield different SA-RNSAF algorithms. We then developed adaptive techniques for the step-size and the sparsity penalty weight in the SA-RNSAF algorithm, thus arriving at the AOP-SA-RNSAF algorithm with a further performance improvement in terms of the convergence rate and steady-state misadjustment. Simulations for the sparse system identification have demonstrated the effectiveness of the proposed algorithms.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\newpage
\bibliographystyle{IEEEtran}


\section{Introduction}
We will consider several graph coloring problems.
In the \emph{list coloring problem}, we have a graph $G$ and a list $L(v)$ of colors at each vertex $v \in V(G)$. In this setting, we say that an \emph{$L$-coloring} of $G$ is a proper coloring $\varphi:V(G) \rightarrow \bigcup_{v \in V(G)} L(v)$ of $G$ in which $\varphi(v) \in L(v)$ for every vertex $v \in V(G)$. If $G$ always has an $L$-coloring whenever $|L(v)| = f(v)$ for each vertex $v \in V(G)$, then we say that $G$ is \emph{$f$-choosable}. If $f$ is a constant function $f(v) = k$, then we say that the \emph{list-chromatic number} (or \emph{choosability}) of $G$ is at most $k$, and we write $\ch(G) \leq k$.

The \emph{DP-coloring problem} is a generalization of the list coloring problem introduced by Dvo\v{r}\'ak and Postle \cite{DP}, defined as follows.
Given a graph $G$ and a function $f:V(G) \rightarrow \mathbb N$, an \emph{$f$-fold cover} of $G$ is a graph $H$ obtained by the following process:
\begin{itemize}
\item For each vertex $v \in V(G)$, add a clique $K_{f(v)}$ to $H$, and write $L(v)$ for the vertex set of this clique.
\item For each edge $uv \in E(G)$, add a matching between $L(u)$ and $L(v)$.
\end{itemize}
Then, we say that an independent set in $H$ of size $|V(G)|$ is a \emph{DP-coloring} of $G$ with respect to $H$. If $G$ always has a DP-coloring for every $f$-fold cover $H$ of $G$, then we say that $G$ is DP-$f$-colorable.
If $f$ is a constant function $f(v) = k$, then we say that $H$ is a \emph{$k$-fold} cover of $G$, and if $G$ always has a DP-coloring for every $k$-fold cover $H$ of $G$, then we say that the DP-chromatic number of $G$ is at most $k$, and we write $\chi_{DP}(G) \leq k$.
Given a cover $H$ of $G$, we often refer to the vertices of $H$ as \emph{colors}, and if $c \in L(v)$ for a vertex $v \in V(G)$, then we say that the color $c$ is \emph{above} $v$. Note that when $f(v) = k$ is a constant function, if the cliques in $H$ corresponding to vertices in $G$ are replaced with independent sets, and if each matching between sets $L(u)$ and $L(v)$ is a perfect matching, then $H$ is a $k$-sheeted covering space of $G$, and a DP-coloring of $G$ is equivalent to an independent transversal of the fibers in $H$ above the vertices of $G$ (see \cite{Hatcher} for an introduction to graphs as topological spaces).

Every list-coloring problem can be transformed into a DP-coloring problem as follows. Given a graph $G$ with a color list $L'(v)$ at every vertex $v \in V(G)$, we construct a cover $H$ of $G$ by adding a clique with vertex set $L(v)$ for every vertex $v \in V(G)$, with elements of $L(v)$ corresponding to colors in $L'(v)$. Then, we consider each edge $uv \in E(G)$, and we add an edge in $H$ between each pair $(c,c') \in L(u) \times L(v)$ for which $c$ and $c'$ both correspond to a common color from $L'(u) \cap L'(v)$. When $H$ is constructed this way, a DP-coloring of $G$ with respect to $H$ is equivalent to an $L'$-coloring of $G$. Therefore, it holds that $\ch(G) \leq \chi_{DP}(G)$.

We will also consider two online graph coloring problems. The \emph{online DP-coloring} problem takes place in the form of a DP-coloring game between two players, called Lister and Painter. The game is played on a graph $G$ with a function $f:V(G) \rightarrow \mathbb N$. At the beginning of the game, each vertex $v \in V(G)$ has $f(v)$ tokens. On each turn $i$, Lister removes some number $m_i(v)$ (possibly zero) of tokens from each vertex $v \in V(G)$ and then reveals a clique $K_{m_i(v)}$ above each vertex $v$. Furthermore, for each edge $uv \in E(G)$, Lister reveals a matching between the revealed cliques above $u$ and $v$. The cliques and matchings revealed on this turn $i$ form a cover $H_i$ of $G$. After $H_i$ is revealed, Painter chooses an independent set from the vertices of $H_i$. The game ends when $G$ has no more tokens for Lister to remove. Painter wins the game if she manages to choose at least one color above each vertex of $G$ before the game is over; otherwise, Lister wins. If Painter always has a winning strategy in the DP-coloring game on a graph $G$ when each vertex $v \in V(G)$ begins with $k$ tokens, then we say that the \emph{online DP-chromatic number} (or \emph{DP-paintability}) of $G$ is at most $k$, and we write $\chi_{DPP}(G) \leq k$. Observe that if each vertex of $G$ begins with $k$ tokens, then Lister has the option of revealing a $k$-fold cover $H$ of $G$ on the first turn and asking Painter to find a DP-coloring of $G$ with respect to $H$, and therefore $\chi_{DP}(G) \leq \chi_{DPP}(G)$.

If, in the DP-coloring game, Lister is only allowed to remove at most one token from each vertex of $G$ during each turn and must always reveals edges wherever possible, then we call this variant of the game the \emph{list-coloring game}. For the list-coloring game, we may equivalently imagine that on each turn $i$, Lister reveals a single color $c_i$ above each vertex of some induced subgraph $G'_i$ of $G$, and Painter must choose some independent set $I_i$ of $G'_i$ and color each vertex in $I_i$ with $c_i$.
In this equivalent setting, each vertex $v \in V(G)$ still begins with $f(v)$ tokens, and a single token is removed from $v$ whenever Lister reveals a color above $v$.
In this setting, Painter wins the game if and only if she can color every vertex of $G$ before the game ends.
 If Painter always has a strategy to win the list-coloring game on a graph $G$ when each vertex $v \in V(G)$ begins with $f(v)$ tokens, then we say that $G$ is \emph{$f$-paintable}. If $f$ is a constant function $f(v) = k$, then we say that $G$ is \emph{$k$-paintable}, and we write $\chi_P(G) \leq k$. The online list-coloring game was originally invented using this framework of revealing colors above vertices independently by Schauz \cite{Schauz} and Zhu \cite{ZhuP}.

At the end of the list-coloring game on $G$ with a constant function $f(v) = k$, the colors revealed at each vertex $v$ form a set $L(v)$ of $k$ colors, and if Painter wins the game, then Painter completes a proper $L$-coloring of $G$. Since Lister is free to form any list assignment $L$ on the vertices of $G$, it follows that if $G$ is $k$-paintable, then $G$ is also $k$-choosable, and hence $\chi_{\ell}(G) \leq \chi_P(G)$. Also, since the online list-coloring game is at least as difficult for Lister as the DP-coloring game, it also follows that $\chi_P(G) \leq \chi_{DPP}(G)$.

After putting all of the inequalities between these parameters together, we obtain two inequality chains:
\begin{eqnarray*}
 \chi_{\ell}(G)  \leq  & \chi_{DP}(G)  & \leq  \chi_{DPP}(G); \\
 \chi_{\ell}(G)  \leq &  \chi_{P}(G)  & \leq  \chi_{DPP}(G) .
\end{eqnarray*}
Given these inequality chains, it is natural to ask whether the differences between adjacent parameters can be arbitrarily large.
For three out of these four differences, we find an affirmative answer by letting $G = K_{n,n}$.
Indeed, Bernshteyn \cite{Bernshteyn} showed that a graph $G$ of average degree $d$ satisfies $\chi_{DP}(G) = \Omega \left ( \frac{d}{\log d} \right )$, implying that $\chi_{DP}(K_{n,n})  =  \Omega \left (\frac{n}{\log n} \right )$.
Since it is known that $\chi_{\ell}(K_{n,n}) \leq \chi_P(K_{n,n}) = \log_2 n + O(1)$ \cite{Duraj}, this shows that
\[\chi_{DP}(K_{n,n}) - \chi_{\ell}(K_{n,n}) = \Omega \left ( \frac{n}{\log n} \right )\textrm{ \indent and \indent} \chi_{DPP}(K_{n,n}) - \chi_{P}(K_{n,n}) = \Omega \left ( \frac{n}{\log n} \right ).\]
 Duraj, Gutowski, and Kozik \cite{Duraj} also showed that  \[\chi_{P}(K_{n,n}) - \chi_{\ell}(K_{n,n}) = \Omega(\log \log n).\]
Therefore, by letting $G = K_{n,n}$, we achieve an arbitrarily large difference for each adjacent pair of parameters except for $\chi_{DPP}(G) - \chi_{DP}(G)$. For this final difference, Kim, Kostochka, Li, and Zhu \cite{KKLZ} showed there exist graphs $G$ for which $\chi_{P}(G) - \chi_{DP}(G) \geq 1$, implying that this final difference can be positive. However, it has not been shown that this difference can be arbitrarily large.

In this note, we will show that the difference $\chi_{DPP}(G) - \chi_{DP}(G)$ can also be arbitrarily large, answering a question of Kim, Kostochka, Li, and Zhu \cite{KKLZ}. Rather than choosing $G = K_{n,n}$, we will construct a graph $G_t$ for each $t \geq 1$ that satisfies $\chi_{DPP}(G_t) - \chi_{DP}(G_t) \geq t$. We construct our graphs $G_t$ by generalizing an idea from the original paper of Kim, Kostochka, Li, and Zhu \cite{KKLZ}. Our graphs $G_t$ will also satisfy the additional property that $\chi_{P}(G_t) - \chi_{DP}(G_t) \geq t$.
By combining this equality with the already known estimate $\chi_{DP}(K_{n,n}) - \chi_P(K_{n,n}) = \Omega \left ( \frac{n}{\log n} \right )$, we hence see that the difference $\chi_{DP}(G) - \chi_{P}(G)$ can achieve both positive and negative values of arbitrarily large magnitude.

\section{The construction}
For each integer $t \geq 1$, we will construct a graph $G_t$ that satisfies $\chi_P(G_t) - \chi_{DP}(G_t) \geq t$. Since $\chi_{DPP}(G) \geq \chi_P(G)$ for all graphs $G$, our graphs $G_t$ will also satisfy
 $\chi_{DPP}(G_t) - \chi_{DP}(G_t) \geq t$.
Our construction is based on a generalization of an idea of Kim, Kostochka, Li, and Zhu~\cite{KKLZ}.

As we are concerned with showing that the paintability of each graph $G_t$ is large enough, we begin with an observation about the online list-coloring game.
If Lister and Painter play the online list-coloring game on a graph $G$ with some initial token assignment, then Lister wins if and only if he can reach a position in which
each uncolored vertex $v \in V(G)$ has some $g(v)$ remaining tokens, and the uncolored subgraph of $G$ is not $g$-choosable.
In the original paper of Kim, Kostochka, Li, and Zhu~\cite{KKLZ}, the authors take advantage of this idea in order to construct a graph $G$ satisfying $\chi_P(G) \geq \chi_{DP}(G) + 1$.
In order to show that their graph $G$ has a large enough paintability, these authors show that in the online list-coloring game on their graph $G$, Lister always has a strategy to create an uncolored $K_{1,k}$ subgraph of $G$ in which each leaf $\ell$ has $g(\ell) = 1$ token and the center vertex $v$ has $g(v) = k$ tokens. Since $K_{1,k}$ is not $g$-choosable, it follows that Lister has a strategy to win the online list-coloring game on their graph $G$.

In our construction, we will use a similar idea.
We first fix the value $k = 2^{8t^3}$. (With more careful calculation, our proof would work with a smaller value of $k$, but we use this larger value for clearer presentation.)
In each of our graphs $G_t$, we will show that Lister can always create an uncolored $K_{t,k^t}$ subgraph in which each $t$-degree vertex $u$ has $g(u) = t$ tokens and each $k^t$-degree vertex $v$ has $g(v) = k$ tokens.
The following lemma shows that if Lister manages to create such a subgraph of $G_t$, then Lister wins the online list-coloring game.
\begin{lemma}
\label{lem:ktkt}
Given the function $g: V(K_{t,k^t}) \rightarrow \mathbb N$ defined above, $K_{t,k^t}$ is not $g$-choosable.
\end{lemma}
\begin{proof}
We let the $t$ vertices $v_1, \dots, v_t$ of degree $k^t$ have pairwise disjoint color lists $L(v_1), \dots, L(v_t)$. Then, for each of the $k^t$ elements $L \in L(v_1) \times \dots \times L(v_t)$, we let $L$ be the list of a vertex $u$ of degree $t$. Then, for any $L$-coloring of $v_1, \dots, v_t$, some vertex $u$ of degree $t$ has no available colors in its list, and hence $K_{t,k^t}$ is not $g$-choosable.
\end{proof}

The most important piece of our construction of $G_t$ will be the following gadget $H_t$. We construct our gadget $H_t$ along with a function $h:V(H_t) \rightarrow \mathbb N$ as follows. We let $H_t$ contain $(t+1)k^t$ copies $K^1, \dots, K^{(t+1)k^t}$ of the clique $K_{t+1}$, and we write $u^{\ell}_1, \dots, u^{\ell}_{t+1}$ for the vertices of each clique $K^{\ell}$. We write $U$ for the set of all of these vertices of the form $u_j^{\ell}$; in other words, we let $U$ consist of all vertices that we have introduced so far. Then, for each value $1 \leq j \leq t+1$, we add $t$ independent vertices $x_j^1,\dots,x_j^t$, and we make each of these vertices adjacent to $u^1_j, u^2_j,\dots,u^{(t+1)k^t}_j$. We write $X$ for the set consisting of all of these vertices of the form $x_j^i$. For each vertex $u_j^{\ell} \in U$, we let $h(u_j^{\ell}) = t+1$, and for each vertex $x_j^i \in X$, we let $h(x_j^i) = k - t + 1$. We now prove two lemmas that show that under appropriate circumstances, winning the online list-coloring game on $H_t$ as Painter is much harder than finding a DP-coloring $H_t$.

\begin{lemma}
\label{lem:no}
$H_t$ is not $(h+t-1)$-paintable.
\end{lemma}
\begin{proof}
We give each vertex $v \in V(H_t)$ exactly $h(v) + t - 1$ tokens, and we show that Painter cannot win the list-coloring game on $H_t$.

On each of the first $t$ turns, Lister reveals a color at each vertex of each clique $K^{\ell}$ and reveals an edge wherever possible. After these first $t$ turns, each clique $K^{\ell}$ must have an uncolored vertex $u_j^{\ell}$ with exactly $h(u_j^{\ell}) - 1 = t$ tokens. Furthermore, since we have $(t+1)k^t$ cliques $K^{\ell}$, each with at least one uncolored vertex, there must exist some value $1 \leq j^* \leq t+1$ for which at least $k^t$ vertices of the form $u_{j^*}^{\ell}$ are uncolored. However, the $k^t$ vertices of the form  $u_{j^*}^{\ell}$ along with the $t$ vertices $x^1_{j^*}, \dots, x^t_{j^*}$ form an uncolored $K_{t,k^t}$ subgraph in which each $t$-degree vertex $v$ has only $g(v) = t$ remaining tokens, and each $k^t$-degree vertex $v$ has only $g(v) = k$ remaining tokens. Since Lemma \ref{lem:ktkt} shows that this $K_{t,k^t}$ subgraph is not $g$-choosable, Lister has a strategy to win the game.
\end{proof}

\begin{lemma}
\label{lem:yes}
$H_t$ is DP-$h$-colorable.
\end{lemma}
\begin{proof}
Consider an $h$-fold cover $H'$ of $H_t$. Recall that given a vertex $v \in V(H_t)$, we say that $v$ corresponds to a clique $K_{h(v)}$ in $H'$ with a vertex set $L(v)$, and we say that $L(v)$ contains $h(v)$ \emph{colors}. Using this terminology, we say that each color in $H'$ appears above only one vertex of $H_t$, as the sets $L(v)$ forming the cliques of $H'$ are pairwise disjoint.

We will first color the vertices $x^i_j \in X$. For each vertex $x^i_j \in X$ and color $c \in L(x^i_j)$, we assign a set $S_c \subseteq [(t+1)k^t]$ that consists of those indices $\ell$ for which $L(u_j^{\ell})$ contains a color adjacent to $c$. We would like to color the $t(t+1)$ vertices of $X$ using $t(t+1)$ colors $c_1, \dots, c_{t(t+1)}$ that correspond to a family $\mathcal S =\{ S_{c_1}, \dots, S_{c_{t(t+1)}}\}$ such that for each value $1 \leq q \leq t+1$, the following property holds:
\begin{equation}\tag{$\star$}
\label{star}
\textrm{The intersection of any $q$ sets of $\mathcal S$ contains at most $k^{t  - \frac{qt}{t+1}} - 1$ elements.}
\end{equation}
In particular, the intersection of any $t+1$ sets of $\mathcal S$ is empty.

We show that we may greedily color each vertex $x_j^i \in X$ while satisfying (\ref{star}). Suppose we wish to color some vertex $x \in X$ and that we have already colored some subset $Y \subseteq X$ while satisfying (\ref{star}). For each subset $A \subseteq Y$ of size $q \in [0,t]$ whose vertices are colored with colors $a_1, \dots, a_q$, we must choose some color $c \in L(x)$ for which
\begin{equation}\tag{$\bullet$}
\label{eqn:int}
|S_c \cap S_{a_1} \cap \dots \cap S_{a_q}| \leq k^{t - \frac{(q+1)t}{t+1}} - 1.
\end{equation}
Note that since $h(u^{\ell}_j) = t+1$ for each vertex $ u^{\ell}_j \in U$, each value $\ell \in [(t+1)k^t]$ appears in at most $t+1$ sets $S_c$ for colors $c \in L(x)$. Furthermore, since $|S_{a_1} \cap \dots \cap S_{a_q}| < k^{t - \frac{qt}{t+1}}$ whenever $q \geq 1$, the number of colors $c \in L(x)$ that do not satisfy (\ref{eqn:int}) for a given $A \subseteq X$ is at most
\[\frac{(t+1)k^{t - \frac{qt}{t+1}} }{k^{t - \frac{(q+1)t}{t+1}}} = (t+1)k^{\frac{t}{t+1}}.\]
Furthermore, since fewer than $2^{t(t+1)}$ subsets $A \subseteq Y$ can be chosen, the number of colors $c \in L(x)$ that would cause (\ref{star}) to be violated is less than
\begin{eqnarray*}
2^{t(t+1)}(t+1) k^{\frac{t}{t+1}} &= & 2^{t(t+1) + \frac{8t^4}{t+1}}(t+1) \\
 &=&2^{t(t+1) + 8(t^3 - t^2 + t - 1 + \frac{1}{t+1} )}(t+1) \\
 &< & 2^{8t^3 - t} \\
 & < & k - t + 1.
\end{eqnarray*}
Therefore, some color $c \in L(x)$ can be used to color $x$ without violating (\ref{star}).

Now, with every vertex $x_j^i \in X$ colored, and with (\ref{star}) satisfied, no index $\ell$ belongs to the intersection of more than $t$ sets $S_c$, where $c$ is a color used at a vertex $x_j^i$, and hence after coloring the vertices in $X$, at most $t$ colors are unavailable at each clique $K^{\ell}$. Therefore, the vertices of each clique $K^{\ell}$ can be ordered so that the first vertex has at least one available color, the second vertex has at least two available colors, and so forth, until the last vertex has $t+1$ available colors. Therefore, each remaining clique $K^{\ell}$ can be $DP$-colored with its available colors, and the lemma is proven.
\end{proof}

Now, we construct our graph $G_t$. First, we make $k^{k-2t}$ copies of the graph $H_t$, and we index them by the $(k-2t)$-tuples in $[k]^{k-2t}$. We also add $k-2t$ vertices $y_1, \dots, y_{k-2t}$ that are adjacent to all vertices of $U$ in each copy of $H_t$. We write $\tilde{U}$ for this set of neighbors of $y_1, \dots, y_{k-2t}$. The following two theorems show that $\chi_P(G_t) - \chi_{DP}(G_t) \geq t$.

\begin{theorem}
$\chi_{DP}(G_t) \leq k-t+1$.
\end{theorem}
\begin{proof}
We may give $G_t$ a DP-coloring with lists of size $(k-t+1)$ as follows. First, we arbitrarily color the vertices $y_1, \dots, y_{k-2t}$. Next, we observe that
the vertices in $\tilde U$ have lost at most $k-2t$ available colors, so
for each vertex $v$ in a copy of $H_t$, $v$ has at least $h(v)$ available colors remaining. Therefore, by Lemma \ref{lem:yes}, we may finish our DP-coloring of $G_t$ by giving each remaining copy of $H_t$ a DP-$h$-coloring.
\end{proof}

\begin{theorem}
\label{thm:no}
$\chi_P(G_t) > k$.
\end{theorem}
\begin{proof}
Suppose that the online list-coloring game is played on $G_t$ with $k$ tokens at each vertex. We will show that Lister has a winning strategy.
For each pair $(i,j)$ satisfying $1 \leq i \leq k$
and
 $1 \leq j \leq k-2t$, Lister executes the following command. When Lister executes the command for a given pair $(i,j)$, we say that this takes place on \emph{turn $(i,j)$}.
\begin{quote}
 Reveal a color $c_{i,j}$ above $y_j$ and above each vertex of $\tilde{U}$ that belongs to a copy of $H_t$ indexed by a $(k-2t)$-tuple
 with the value $i$ in the $j$th coordinate.
\end{quote}
For each value $j \in [k-2t]$, we write $L(y_j) = \{c_{1,j}, \dots, c_{k,j}\}$ for the set of colors revealed above $y_j$.

For each $j \in [k-2t]$,
we may assume that for some value $i_{j} \in [k]$, Painter colors $y_{j}$
during turn $(i_{j}, j)$
and hence does not color any vertex of $\tilde U$ during turn $(i_j, j)$.
Indeed, if this is not the case, then $y_j$ is never colored, and Painter will not have another chance to color $y_{j}$.
Therefore,
for each value $j \in [k-2t]$, we may assume that
no vertex in a copy of $H_t$ indexed by a $(k-2t)$-tuple with an $i_j$ entry in the $j$th coordinate is colored
using a color in $L(y_j)$.

Now, let $H$ be the copy of $H_t$ indexed by the $(k-2t)$-tuple $(i_1, \dots, i_{k-2t})$. By our observation above, $H$ has not been colored by a color in $L(y_1) \cup \dots \cup L(y_{k-2t})$, and hence no vertex of $H$ has been colored.
Furthermore,
since $k-2t$ tokens have been removed from each vertex in $\tilde{U} \cap V(H)$, it follows that
for each vertex $v \in V(H)$, only $h(v)+t-1$ tokens remain at $v$. Therefore, Lister can follow the strategy in Lemma \ref{lem:no} on $H$ in order to win the game, and thus $\chi_P(G_t) > k$.
\end{proof}

\section{Conclusion}
While we have shown for each $t \geq 1$ the existence of a graph $G_t$ for which $\chi_P(G_t) - \chi_{DP}(G_t) \geq t$, it is still open whether there exists a sequence $\{G_t\}_{t \geq 1}$ of graphs for which
\[\lim_{t \rightarrow \infty} \frac{\chi_{DPP}(G_t)}{\chi_{DP}(G_t)} > 1 \textrm{ \indent or \indent } \lim_{t \rightarrow \infty} \frac{\chi_{P}(G_t)}{\chi_{\ell}(G_t)} > 1.\]
On the other hand, it is unknown whether $\chi_P(G)$ can be bounded above by a linear or even polynomial function of $\chi_{\ell}(G)$, and it is unknown whether $\chi_{DPP}(G)$ can be bounded above by a linear function of $\chi_{DP}(G)$. Duraj, Gutowski, and Kozik \cite{Duraj} have pointed out that currently, the best known bound for $\chi_P(G)$ in terms of $\chi_{\ell}(G)$ comes from the inequality chain
\[\chi_P(G) \leq d(G)+1 \leq 4 {{s^4 } \choose s} \log \left ( 2{{s^4 } \choose s} \right ) + 1,\]
where $d(G)$ is the degeneracy of $G$, $s = \chi_{\ell}(G)$, and the last inequality in the chain is from a result of Alon \cite{Alon}. For $\chi_{DPP}(G)$, the best upper bound in terms of $\chi_{DP}(G)$ comes from the inequality chain
\[\chi_{DPP}(G) \leq d(G) + 1\leq (2 + o(1)) \chi_{DP}(G)  \log \chi_{DP}(G),\]
where $d(G)$ is the degeneracy of $G$, and the last inequality comes from a result of Bernshteyn \cite{Bernshteyn}. It is likely that a deeper understanding of these coloring parameters is necessary to determine tight bounds between these parameters and answer these questions completely.

\section{Acknowledgment}
I am grateful to Alexandr Kostochka for offering helpful feedback on an earlier version of this paper.

\raggedright
\bibliographystyle{plain}