\section{Introduction}
\label{sec1}

AI systems are expected to work in an incremental manner when the amount of knowledge increases over time. They should be capable of learning new concepts while maintaining the ability to recognize previous ones.
However, deep-neural-network-based systems often suffer from serious forgetting problems (called ``catastrophic forgetting'') when they are continuously updated using new coming data. This is due to two facts: (i) the updates can override the knowledge acquired from the previous data~\cite{mccloskey1989catastrophic, McRae1993Catastrophic, Ratcliff1990catastrophic, shin2017continual, KemkerMAHK18}, and (ii) the model can not replay the entire previous data to regain the old knowledge.

To encourage solving these problems, \cite{rebuffi2017icarl} defined a class-incremental learning (CIL) protocol for image classification where the training data of different classes gradually come phase-by-phase.
In each phase, the classifier is re-trained on new class data, and then evaluated on the test data of both old and new classes.
To prevent trivial algorithms such as storing all old data for replaying,
there is a strict memory budget due to which a tiny set of exemplars of old classes can be saved in the memory.
This memory constraint causes a serious data imbalance problem between old and new classes, and indirectly causes the main problem of CIL -- the stability-plasticity dilemma~\cite{mermillod2013stability}.
In particular, higher plasticity results in the forgetting of old classes~\cite{mccloskey1989catastrophic},
while higher stability weakens the model from learning the data of new classes (that contain a large number of samples).
Existing CIL works try to balance stability and plasticity using data strategies. For example, as illustrated in Figure~\ref{fig_new_teaser} (a) and (b), some early methods train their models on the imbalanced dataset where there is only a small set of exemplars for old classes~\cite{rebuffi2017icarl,Li18LWF}, and recent methods include a fine-tuning step using a balanced subset of exemplars sampled from all classes~\cite{Castro18EndToEnd,hou2019lucir,douillard2020podnet}.
However, these data strategies are still limited in terms of effectiveness.
For example, when using the models trained after $25$ phases, LUCIR~\cite{hou2019lucir} and Mnemonics~\cite{liu2020mnemonics} ``forget'' the initial $50$ classes by $30\%$ and $20\%$, respectively, on the ImageNet dataset~\cite{russakovsky2015imagenet}.

In this paper, we address the stability-plasticity dilemma by introducing a novel network architecture called Adaptive Aggregation Networks (AANets). Taking the ResNet~\cite{He_CVPR2016_ResNet} as an example of baseline architectures, we explicitly build two residual blocks (at each residual level) in AANets: one for maintaining the knowledge of old classes (i.e., the stability) and the other for learning new classes (i.e., the plasticity), as shown in Figure~\ref{fig_new_teaser}~(c).
We achieve these by allowing these two blocks to have different levels of learnability, i.e., less learnable parameters in the stable block but more in the plastic one.
We apply aggregation weights to the output feature maps of these blocks, sum them up, and pass the result maps to the next residual level.
In this way, we are able to dynamically balance the usage of these blocks by updating their aggregation weights. To achieve auto-updating, we take the weights as hyperparameters and optimize them in an end-to-end manner~\cite{finn2017model, Wu2019LargeScale, liu2020mnemonics}.

\input{figures/1_framework}

Technically,
the overall optimization of AANets is bilevel. Level-1 is to learn the network parameters for two types of residual blocks, and level-2 is to adapt their aggregation weights.
More specifically, level-1 is the standard optimization of network parameters, for which we use all the data available in the phase.
Level-2 aims to balance the usage of the two types of blocks, for which we optimize the aggregation weights using a balanced subset (by downsampling the data of new classes), as illustrated in Figure~\ref{fig_new_teaser}~(c).
We formulate these two levels in a bilevel optimization program (BOP)~\cite{SinhaMD18bilevelreview} that solves two optimization problems alternatively, i.e., update network parameters with aggregation weights fixed, and then switch.
For evaluation, we conduct CIL experiments on three widely-used benchmarks, CIFAR-100, ImageNet-Subset, and ImageNet. We find that many existing CIL methods, e.g., iCaRL~\cite{rebuffi2017icarl}, LUCIR~\cite{hou2019lucir}, Mnemonics Training~\cite{liu2020mnemonics}, and PODNet~\cite{douillard2020podnet}, can be directly incorporated in the architecture of AANets, yielding consistent performance improvements. We observe that a straightforward plug-in causes memory overheads, e.g., $26\%$ and $15\%$ respectively for CIFAR-100 and ImageNet-Subset. For a fair comparison, we conduct additional experiments under the settings of zero overhead (e.g., by reducing the number of old exemplars for training AANets), and validate that our approach still achieves top performance across all datasets.

\textbf{Our contribution} is three-fold:
1) a novel and generic network architecture called AANets specially designed for tackling the stability-plasticity dilemma in CIL tasks;
2) a BOP-based formulation and an end-to-end training solution for optimizing AANets;
and 3) extensive experiments on three CIL benchmarks by incorporating four baseline methods in the architecture of AANets.
\section{Related Work}
\label{sec2}

\myparagraph{Incremental learning} aims to learn efficient machine models from the data that gradually come in a sequence of training phases.
Closely related topics are referred to as continual learning~\cite{de2019continualsurvey,lopez2017gradient} and lifelong learning~\cite{chen2018lifelong,aljundi2017expert,li2019online}.
Recent incremental learning approaches are either task-based, i.e., all-class data come but are from a different dataset for each new phase~\cite{Li18LWF,shin2017continual,hu2018overcoming,chaudhry2018efficient,riemer2018learning,zhao2020maintaining,davidson2020sequential,chaudhry2018riemannian}, or class-based i.e., each phase has the data of a new set of classes coming from the identical dataset~\cite{rebuffi2017icarl,hou2019lucir,Wu2019LargeScale,Castro18EndToEnd,liu2020mnemonics,hu2021cil,Zhang_2021_CVPR}.
The latter one is typically called class-incremental learning (CIL), and our work is based on this setting.
Related methods mainly focus on how to solve the problems of forgetting old data. Based on their specific methods, they can be categorized into three classes: regularization-based, replay-based, and parameter-isolation-based~\cite{de2019literaturereview,prabhu12356gdumb}.
\\
\myparagraph{\emph{Regularization-based}} methods introduce regularization terms in the loss function to consolidate previous knowledge when learning new data.
Li et al.~\cite{Li18LWF} proposed the regularization term of knowledge distillation~\cite{HintonVD15}.
Hou et al.~\cite{hou2019lucir} introduced a series of new regularization terms such as for less-forgetting constraint and inter-class separation to mitigate the negative effects caused by the data imbalance between old and new classes.
Douillard et al.~\cite{douillard2020podnet} proposed an effective spatial-
based distillation loss applied throughout the model and also a representation comprising multiple proxy vectors for each object class.
Tao et al.~\cite{Tao2020topology} built the framework with a topology-preserving loss to maintain the topology in the feature space.
Yu et al.~\cite{yu2020CVPRsemantic} estimated the drift of previous classes during the training of new classes.
\\
\myparagraph{\emph{Replay-based}} methods store a tiny subset of old data, and replay the model on them (together with new class data) to reduce the forgetting.
Rebuffi et al.~\cite{rebuffi2017icarl} picked the nearest neighbors to the average sample per class to build this subset.
Liu et al.~\cite{liu2020mnemonics} parameterized the samples in the subset, and then meta-optimized them automatically in an end-to-end manner taking the representation ability of the whole set as the meta-learning objective.
Belouadah et al.~\cite{belouadah2019il2m} proposed to leverage a second memory to store  statistics of old classes in rather compact formats.
\\
\myparagraph{\emph{Parameter-isolation-based}} methods are used in task-based incremental learning (but not CIL). Related methods dedicate different model parameters for different incremental phases, to prevent model forgetting (caused by parameter overwritten).
If no constraints on the size of the neural network is given, one can grow new branches for new tasks while freezing old branches.
Rusu et al.~\cite{rusu2016progressive} proposed ``progressive networks'' to integrate the desiderata of different tasks directly into the networks.
Abati et al.~\cite{abati2020conditional} equipped each convolution layer with task-specific gating modules that select specific filters to learn each new task.
Rajasegaran et al.\cite{NIPS2019_9429} progressively chose the optimal paths for the new task while encouraging to share parameters across tasks.
Xu et al.~\cite{xu2018reinforced} searched for the best neural network architecture for each coming task by leveraging reinforcement learning strategies.
\textbf{Our differences} with these methods include the following aspects. We focus on class-incremental learning, and more importantly, our approach does not continuously increase the network size. We validate in the experiments that under a strict memory budget, our approach can surpass many related methods and its plug-in versions on these related methods can bring consistent performance improvements.

\myparagraph{Bilevel Optimization Program}
can be used to optimize hyperparameters of deep models.
Technically, the network parameters are updated at one level and the key hyperparameters are updated at another level~\cite{stackelberg1952theory,Wang2018Distillation,goodfellow2014generative,zhang2020deepemd,Liu2020E3BM,li2019learning}.
Recently, a few bilevel-optimization-based approaches have emerged for tackling incremental learning tasks.
Wu et al.~\cite{Wu2019LargeScale} learned a bias correction layer for incremental learning models using a bilevel optimization framework.
Rajasegaran et al.~\cite{rajasegaran2020itaml} incrementally learned new tasks while learning a generic model to retain the knowledge from all tasks.
Riemer et al.~\cite{riemer2018learning2learn} learned network updates that are well-aligned with previous phases, such as to avoid learning towards any distracting directions.
In our work, we apply the bilevel optimization program to update the aggregation weights in our AANets.
\section{Adaptive Aggregation Networks (AANets)}
\label{sec_method}
\input{figures/2_f_connection}

Class-Incremental Learning (CIL) usually assumes $(N+1)$ learning phases in total, i.e., one initial phase and $N$ incremental phases during which the number of classes gradually increases~\cite{hou2019lucir,liu2020mnemonics,douillard2020podnet,hu2021cil}.
In the initial phase, data $\mathcal{D}_0$ is available to train the first model $\Theta_0$.
There is a strict memory budget in CIL systems, so after the phase, only a small subset of $\mathcal{D}_0$ (exemplars denoted as $\mathcal E_0$) can be stored in the memory and used as replay samples in later phases.
Specifically in the $i$-th ($i\geq1$) phase, we load the exemplars of old classes $\mathcal E_{0:i-1}=\{\mathcal E_0, \dots, \mathcal E_{i-1}\}$ to train model $\Theta_i$ together with new class data $\mathcal{D}_i$.
Then, we evaluate the trained model on the test data containing both old and new classes.
We repeat such training and evaluation through all phases.

The key issue of CIL
is that the models trained at new phases easily ``forget'' old classes.
To tackle this, we introduce a novel architecture called AANets.
AANets is based on a ResNet-type architecture, and each of its residual levels is composed of two types of residual blocks: a plastic one to adapt to new class data and a stable one to maintain the knowledge learned from old classes.
The details of this architecture are elaborated in Section~\ref{subsec_MANets_arch}.
The steps for optimizing AANets are given in Section~\ref{subsec_optimization_MANets}.

\subsection{Architecture Details}
\label{subsec_MANets_arch}

In Figure~\ref{fig_framework_connection}, we provide an illustrative example of our AANets with three residual levels.
The inputs $x^{[0]}$ are the images and the outputs $x^{[3]}$ are the features used to train classifiers.
Each of our residual ``levels'' consists of two parallel residual ``blocks'' (of the original ResNet~\cite{He_CVPR2016_ResNet}): the \textcolor[rgb]{0.93,0.47,0.32}{orange} one (called plastic block) will have its parameters fully adapted to new class data, while the \textcolor[rgb]{0.2,0.475,0.71}{blue} one (called stable block) has its parameters partially fixed in order to maintain the knowledge learned from old classes.
After feeding the inputs to Level 1, we obtain two sets of feature maps respectively from two blocks, and aggregate them after applying the aggregation weights $\alpha^{[1]}$.
Then, we feed the resulted maps to Level 2 and repeat the aggregation.
We apply the same steps for Level 3.
Finally, we pool the resulted maps obtained from Level 3 to train classifiers.
Below we elaborate the details of this dual-branch design as well as the steps for feature extraction and aggregation.

\myparagraph{Stable and Plastic Blocks.}
We deploy a pair of stable and plastic blocks at each residual level, aiming to balance between the plasticity, i.e., for learning new classes, and stability, i.e., for not forgetting the knowledge of old classes.
We achieve these two types of blocks by allowing different levels of learnability, i.e., less learnable parameters in the stable block but more in the plastic.
We detail the operations in the following.
In any CIL phase, Let $\eta$ and $\phi$ represent the learnable parameters of plastic and stable blocks, respectively.
$\eta$ contains all the convolutional weights,
while $\phi$ contains only the neuron-level scaling weights~\cite{sun2019meta}. Specifically, these scaling weights are applied on the model $\theta_{\mathrm{base}}$ obtained in the $0$-th phase\footnote{Related work~\cite{hou2019lucir,douillard2020podnet,liu2020mnemonics} learned $\Theta_0$ in the $0$-th phase using half of the total classes. We follow the same way to train $\Theta_0$ and freeze it as $\theta_{\mathrm{base}}$.}.
As a result, the number of learnable parameters $\phi$ is much less than that of $\eta$.
For example, when using $3\times3$ neurons in $\theta_{\mathrm{base}}$, the number of learnable parameters $\phi$ is only $\frac{1}{3\times 3}$ of the number of full network parameters (while $\eta$ has the full network parameters). We further elaborate on these in the following paragraph.

\myparagraph{Neuron-level Scaling Weights.}
For stable blocks, we learn its neuron parameters in the $0$-th phase and freeze them in the other $N$ phases.
In these $N$ phases,  we apply a small set of scaling weights $\phi$ at the neuron-level, i.e., each weight for scaling one neuron in $\theta_\mathrm{base}$.
We aim to preserve the structural pattern within the neuron and slowly adapt the knowledge of the whole blocks to new class data.
Specifically, we assume the $q$-th layer of $\theta_\mathrm{base}$ contains $R$ neurons, so we have $R$ neuron weights as $\{W_{q,r}\}_{r=1}^R$. For conciseness, we denote them as $W_{q}$.
For $W_{q}$, we learn $R$ scaling weights denoted as $\phi_{q}$
Let $X_{q-1}$ and $X_q$ be the input and output feature maps of the $q$-th layer, respectively.
We apply $\phi_q$ to $W_q$ as follows,
\begin{equation}\label{eq_T_operation_layer}
     X_q =(W_q\odot\phi_{q}) X_{q-1},
\end{equation}
where $\odot$ donates the element-wise multiplication. Assuming there are $Q$ layers in total, the overall scaling weights can be denoted as $\phi=\{\phi_{q}\}_{q=1}^Q$.

\myparagraph{Feature Extraction and Aggregation.}
We elaborate on the process of feature extraction and aggregation across all residual levels in the AANets, as illustrated in Figure~\ref{fig_framework_connection}.
Let $\mathcal{F}^{[k]}_{\mu}{(\cdot)}$ denote
the transformation function of the residual block parameterized as $\mu$ at the Level $k$.
Given a batch of training images $x^{[0]}$, we
feed them to AANets to compute the
feature maps at the $k$-th level (through the stable and plastic blocks respectively) as follows,
\begin{equation}\label{eq_kth_two_branch}
     x_{\phi}^{[k]} = \mathcal{F}^{[k]}_{\phi\odot \theta_{\mathrm{base}}}(x^{[k-1]});
     \quad x_{\eta}^{[k]} = \mathcal{F}^{[k]}_{\eta}(x^{[k-1]}).
\end{equation}
The transferability (of the knowledge learned from old classes) is different
at different levels of neural networks~\cite{yosinski2014transferable}.
Therefore, it makes more sense to apply
different aggregation weights for different levels of residual blocks.
Let $\alpha^{[k]}_{\phi}$ and $\alpha^{[k]}_{\eta}$
denote the
aggregation weights of the stable and plastic blocks, respectively, at the $k$-th level.
Then, the
weighted sum of $x_{\phi}^{[k]}$ and $x_{\eta}^{[k]}$ can be derived as follows,
\begin{equation}\label{eq_kth_alpha_fusion}
     x^{[k]} = \alpha^{[k]}_{\phi} \cdot x_{\phi}^{[k]} + \alpha^{[k]}_{\eta} \cdot x_{\eta}^{[k]}.
\end{equation}
In our illustrative example in Figure~\ref{fig_framework_connection}, there are three pairs of weights to learn at each phase. Hence, it becomes increasingly challenging to choose these weights manually if multiple phases are involved.
In this paper, we propose an learning strategy to automatically adapt these weights, i.e., optimizing the weights for different blocks in different phases, see details in Section~\ref{subsec_optimization_MANets}.

\subsection{Optimization Steps}
\label{subsec_optimization_MANets}

In each incremental phase, we optimize two groups of learnable parameters in AANets: (a) the neuron-level scaling weights $\phi$ for the stable blocks and the convolutional weights
$\eta$ on the plastic blocks;
(b) the feature aggregation weights $\alpha$.
The former is for network parameters and the latter is for hyperparameters.
In this paper, we formulate the overall optimization process as a bilevel optimization program (BOP)~\cite{goodfellow2014generative,liu2020mnemonics}.

\myparagraph{The Formulation of BOP.}
In AANets, the network parameters $[\phi, \eta]$ are trained using the aggregation weights $\alpha$ as hyperparameters.
In turn, $\alpha$ can be updated when temporarily fixing network parameters $[\phi, \eta]$.
In this way, the optimality of $[\phi, \eta]$
imposes a constraint on $\alpha$ and vise versa.
Ideally, in the $i$-th phase, the CIL system aims to learn the optimal $\alpha_i$ and $[\phi_i, \eta_i]$ that minimize the classification loss on all training samples seen so far, i.e., $\mathcal{D}_i\cup\mathcal{D}_{0:i-1}$, so the ideal BOP can be formulated as,
\begin{subequations}
\label{eq_global_bilevel_program_original}
\begin{align}\label{eq_global_bilevel_upper_A}
    &\min_{\alpha_i} \mathcal{L}(\alpha_i, \phi^*_i, \eta^*_i; \mathcal{D}_{0:i-1}\cup \mathcal{D}_{i}) \\
\label{eq_global_bilevel_lower_A}
    &\ \text{s.t.} \ [\phi_i^*, \eta^*_i] = \argmin_{[\phi_i, \eta_i]} \mathcal{L}(\alpha_i, \phi_i, \eta_i; \mathcal{D}_{0:i-1}\cup \mathcal{D}_i),
\end{align}
\end{subequations}
where $\mathcal{L}(\cdot)$ denotes the loss function, e.g., cross-entropy loss. Please note that for the conciseness of the formulation, we use $\phi_i$ to represent $\phi_i\odot \theta_{\mathrm{base}}$ (same in the following equations).
We call Problem~\ref{eq_global_bilevel_upper_A} and Problem~\ref{eq_global_bilevel_lower_A} the \emph{upper-level} and
\emph{lower-level} problems, respectively.

\myparagraph{Data Strategy.}
To solve Problem~\ref{eq_global_bilevel_program_original}, we need to use $\mathcal{D}_{0:i-1}$. However, in the setting of CIL~\cite{rebuffi2017icarl,hou2019lucir,douillard2020podnet}, we cannot access $\mathcal{D}_{0:i-1}$ but only a small set of exemplars $\mathcal{E}_{0:i-1}$, e.g., $20$ samples of each old class.
Directly replacing $\mathcal{D}_{0:i-1}\cup \mathcal{D}_{i}$ with $\mathcal{E}_{0:i-1}\cup \mathcal{D}_{i}$ in Problem~\ref{eq_global_bilevel_program_original} will lead to the forgetting problem for the old classes.
To alleviate this issue, we propose a new data strategy in which we use different training data splits to learn different groups of parameters:
1) in the \emph{upper-level} problem, $\alpha_i$ is used to balance the stable and the plastic blocks, so we use the balanced subset to update it, i.e., learning $\alpha_i$ on $\mathcal{E}_{0:i-1}\cup \mathcal{E}_{i}$ adaptively;
2) in the \emph{lower-level} problem, $[\phi_i, \eta_i]$ are the network parameters used for feature extraction, so we
leverage all the available data
to train them, i.e., base-training $[\phi_i, \eta_i]$ on $\mathcal{E}_{0:i-1}\cup \mathcal{D}_{i}$.
Based on these, we can reformulate the ideal BOP in Problem~\ref{eq_global_bilevel_program_original} as a solvable BOP as follows,
\begin{subequations}
\label{eq_global_bilevel_program}
\begin{align}\label{eq_global_bilevel_upper}
    &\min_{\alpha_i} \mathcal{L}(\alpha_i, \phi^*_i, \eta^*_i; \mathcal{E}_{0:i-1}\cup \mathcal{E}_{i}) \\
\label{eq_global_bilevel_lower}
    &\ \text{s.t.} \ [\phi_i^*, \eta^*_i] = \argmin_{[\phi_i, \eta_i]} \mathcal{L}(\alpha_i, \phi_i, \eta_i; \mathcal{E}_{0:i-1}\cup \mathcal{D}_i),
\end{align}
\end{subequations}
where Problem~\ref{eq_global_bilevel_upper} is the \emph{upper-level} problem and Problem~\ref{eq_global_bilevel_lower} is the \emph{lower-level} problem we are going to solve.

\myparagraph{Updating Parameters.}
We solve Problem~\ref{eq_global_bilevel_program} by alternatively updating two groups of parameters ($\alpha_i$ and $[\phi, \eta]$) across epochs, e.g., if $\alpha_i$ is updated in the $j$-th epoch, then $[\phi, \eta]$ will be updated in the $(j+1)$-th epoch, until both of them converge.
Taking the $i$-th phase as an example, we initialize ${\alpha}_i, \phi_i, \eta_i$ with ${\alpha}_{i-1}, \phi_{i-1}, \eta_{i-1}$, respectively.
Please note that $\phi_0$ is initialized with ones, following~\cite{sun2019meta,sun2019mtlj}; $\eta_0$ is initialized with $\theta_{\mathrm{base}}$; and ${\alpha}_0$ is initialized with $0.5$.
Based on our \textbf{Data Strategy}, we use all available data in the current phase to solve the \emph{lower-level} problem, i.e., training [$\phi_i$, $\eta_i$] as follows,
\begin{equation}\label{eq_model_update_1}
    [\phi_i, \eta_i] \gets  [\phi_i, \eta_i] - \gamma_1 \nabla_{[\phi_i, \eta_i]}\mathcal{L}(\alpha_i, \phi_i, \eta_i; \mathcal{E}_{0:i-1}\cup \mathcal{D}_i).
\end{equation}
Then, we use a balanced exemplar set
to solve the \emph{upper-level} problem, i.e., training ${\alpha}_i$ as follows,
\begin{equation}\label{eq_model_update_2}
    {\alpha}_i \gets {\alpha}_i - \gamma_2 \nabla_{{\alpha}_i}\mathcal{L}({\alpha}_i, \phi_i, \eta_i; \mathcal{E}_{0:i-1}\cup \mathcal{E}_i),
\end{equation}
where $\gamma_1$ and $\gamma_2$ are the \emph{lower-level} and \emph{upper-level} learning rates, respectively.

\input{pseudocode/algo}

\subsection{Algorithm}
\label{subsec_algo}

In Algorithm~\ref{alg:manet}, we summarize the overall training steps of the proposed AANets in the $i$-th incremental learning phase (where $i \in [1,...,N]$). Lines 1-4 show the preprocessing including loading new data and old exemplars (Line 1), initializing the two groups of learnable parameters (Lines 2-3), and selecting the exemplars for new classes (Line 4). Lines 5-12 optimize alternatively between the network parameters and the Adaptive Aggregation weights. In specific, Lines 6-8 and Lines 9-11 execute the training for solving the \emph{upper-level} and \emph{lower-level} problems, respectively. Lines 13-14 update the exemplars and save them to the memory.
\section{Experiments}
\label{sec_exp}

We evaluate the proposed {AANets}
on three CIL benchmarks, i.e., CIFAR-100~\cite{krizhevsky2009learning}, ImageNet-Subset~\cite{rebuffi2017icarl} and ImageNet~\cite{russakovsky2015imagenet}.
We incorporate AANets into four baseline methods and boost their model performances consistently for all settings.
Below we describe the datasets and implementation details (Section~\ref{subsec_datasets}), followed by the results and analyses (Section~\ref{subsec_results}) which include a detailed ablation study, extensive comparisons to related methods, and some visualization of the results.

\subsection{Datasets and Implementation Details}
\label{subsec_datasets}
\input{tables/ablation_2}

\myparagraph{Datasets.}
We conduct CIL experiments on two datasets, CIFAR-100~\cite{krizhevsky2009learning} and ImageNet~\cite{russakovsky2015imagenet}, following closely related work~\cite{hou2019lucir,liu2020mnemonics,douillard2020podnet}. CIFAR-100 contains $60,000$ samples of $32\times32$ color images for $100$ classes. There are $500$ training and $100$ test samples for each class. ImageNet contains around $1.3$ million samples of $224\times224$ color images for $1000$ classes. There are approximately $1,300$ training and $50$ test samples for each class.
ImageNet is used in two CIL settings: one based on a subset of $100$ classes (ImageNet-Subset) and the other based on the full set of $1,000$ classes.
The $100$-class data for {ImageNet-Subset}
are sampled from ImageNet in the same way as~\cite{hou2019lucir,douillard2020podnet}.

\myparagraph{Architectures.}
Following the exact settings in~\cite{hou2019lucir,liu2020mnemonics}, we deploy a $32$-layer ResNet as the baseline architecture (based on which we build the AANets) for CIFAR-100.
This ResNet consists of $1$ initial convolution layer and $3$ residual blocks (in a single branch). Each block has $10$ convolution layers with $3\times3$ kernels. The number of filters starts from $16$ and is doubled every next block. After these $3$ blocks, there is an average-pooling layer to compress the output feature maps to a feature embedding. To build AANets, we convert these $3$ blocks into three levels of blocks and each level consists of a stable block and a plastic block, referring to Section~\ref{subsec_MANets_arch}.
Similarly, we build AANets for ImageNet benchmarks but taking an $18$-layer ResNet~\cite{He_CVPR2016_ResNet} as the baseline architecture~\cite{hou2019lucir,liu2020mnemonics}.
Please note that there is no architecture change applied to the classifiers, i.e., using the same FC layers as in~\cite{hou2019lucir,liu2020mnemonics}.

\myparagraph{Hyperparameters and Configuration.}
The learning rates $\gamma_1$ and $\gamma_2$ are initialized as $0.1$ and $1\times10^{-8}$, respectively.
We impose a constraint on each pair of $\alpha_{\eta}$ and $\alpha_{\phi}$ to make sure $\alpha_{\eta}+\alpha_{\phi}=1$.
For fair comparison, our training hyperparamters are almost the same as in~\cite{douillard2020podnet,liu2020mnemonics}.
Specifically, on the CIFAR-100 (ImageNet), we train the model for $160$ ($90$) epochs in each phase, and the learning rates are divided by $10$ after $80$ ($30$) and then after $120$ ($60$) epochs.
We use an SGD optimizer with the momentum $0.9$ and the batch size $128$ to train the models in all settings.

\myparagraph{Memory Budget.}
By default, we follow the same data replay settings used in~\cite{rebuffi2017icarl,hou2019lucir,liu2020mnemonics,douillard2020podnet}, where each time reserves $20$ exemplars per old class.
In our ``strict memory budget'' settings, we strictly control the memory budget shared by the exemplars and the model parameters.
For example, if we incorporate AANets to LUCIR~\cite{hou2019lucir}, we need to reduce the number of exemplars to balance the additional memory used by model parameters (as AANets take around $20\%$ more parameters than plain ResNets).
As a result, we reduce the numbers of exemplars for AANets from $20$ to $13$, $16$ and $19$, respectively, for CIFAR-100, ImageNet-Subset, and ImageNet, in the ``strict memory budget'' setting.
For example, on CIFAR-100, we use $530k$ additional parameters, so we need to reduce $530k \mathrm{floats}\times4\mathrm{bytes/float}\div(32\times32\times3\mathrm{bytes}/\mathrm{image})\div100\mathrm{classes}\approx7\mathrm{images}/\mathrm{class}$.

\myparagraph{Benchmark Protocol.}
We follow the common protocol used in \cite{hou2019lucir,liu2020mnemonics,douillard2020podnet}.
Given a dataset, the model is trained on half of the classes in the $0$-th phase.
Then,
it learns
the remaining classes evenly in the subsequent $N$ phases.
For $N$, there are three options as $5$, $10$, and $25$, and the corresponding settings are called ``$N$-phase''.
In each phase, the model is evaluated on the test data for all seen classes. The average accuracy (over all phases) is reported.
For each setting, we run the experiment three times and report averages and $95\%$ confidence intervals.

\input{tables/sota}
\subsection{Results and Analyses}
\label{subsec_results}

Table~\ref{table:ablation} summarizes the statistics and results in $8$ ablative settings.
Table~\ref{table_sota} presents the results of $4$ state-of-the-art methods \emph{w/} and
\emph{w/o} AANets as a plug-in architecture, and the reported results from some other comparable work.
Figure~\ref{figure_gradcam} compares the activation maps (by Grad-CAM~\cite{selvaraju2017gradcam}) produced by different types of residual blocks and for the classes seen in different phases.
Figure~\ref{figure_values_plots} shows the changes of values of $\alpha_{\eta}$ and $\alpha_{\phi}$ across $10$ incremental phases.

\myparagraph{Ablation Settings.}
Table~\ref{table:ablation} shows the ablation study.
By differentiating the numbers of learnable parameters, we can have $3$ block types: 1) ``all'' for learning all the convolutional weights and biases; 2) ``scaling'' for learning neuron-level scaling weights~\cite{sun2019meta} on the top of a frozen base model $\theta_{\mathrm{base}}$; and 3) ``frozen'' for using only $\theta_{\mathrm{base}}$ (always frozen). In Table~\ref{table:ablation}, the pattern of combining blocks is A+B where A and B stands for the plastic and the stable blocks, respectively.
Rows 1 is the baseline method LUCIR~\cite{hou2019lucir}.
Row 2 is a double-branch version for LUCIR without learning any aggregation weights.
Rows 3-5 are our AANets using different combinations of blocks.
Row 6-8 use ``all''+``scaling'' under an additional setting as follows. 1) Row 6 uses imbalanced data $\mathcal{E}_{0:i-1}\cup \mathcal{D}_{i}$ to train $\alpha$ adaptively. 2) Row 7 uses fixed weights $\alpha_{\eta}=\alpha_{\phi}=0.5$ at each residual level. 3) Row 8 is under the ``strict memory budget'' setting, where we reduce the numbers of exemplars to $14$ and $17$ for CIFAR-100 and ImageNet-Subset, respectively.

\myparagraph{Ablation Results.}
In Table~\ref{table:ablation}, if comparing the second block (ours) to the first block (single-branch and double-branch baselines), it is obvious that using AANets can clearly improve the model performance, e.g., ``scaling''+``frozen'' gains an average of $4.8\%$ over LUCIR for the ImageNet-Subset, by optimizing $1.4M$ parameters during CIL --- only $12.6\%$ of that in LUCIR.
Among Rows 3-5, we can see that for the ImageNet-Subset, models with the fewest learnable parameters (``scaling''+``frozen'') work the best.
We think this is because we use shallower networks for learning larger datasets (ResNet-32 for CIFAR-100; ResNet-18 for ImageNet-Subset), following the \textbf{Benchmark Protocol}.
In other words, $\theta_{\mathrm{base}}$ is quite well-trained with the rich data of half ImageNet-Subset ($50$ classes in the $0$-th phase), and can offer high-quality features for later phases.
Comparing Row 6 to Row 3 shows the efficiency of using a balanced subset to optimize $\alpha$.
Comparing Row 7 to Row 3 shows the superiority of learning $\alpha$ (which is dynamic and optimal) over manually-choosing $\alpha$.
\\
\textbf{About the Memory Usage.}
By comparing Row 3 to Row 1, we can see that AANets can clearly improve the model performance while introducing small overheads for the memory, e.g., $26\%$ and $14.5\%$ on the CIFAR-100 and ImageNet-Subset, respectively.
If comparing Row~8 to Row~3, we find that though the numbers of exemplars are reduced (for Row~8), the model performance of AANets has a very small drop, e.g., only $0.3\%$ for the $5$-Phase CIL models of CIFAR-100 and ImageNet-Subset. %
Therefore, we can conclude that AANets can achieve rather satisfactory performance under strict memory control --- a desirable feature needed in class-incremental learning systems.

\input{figures/5_gradcam}
\input{figures/3_plot}

\myparagraph{Comparing to the State-of-the-Art.}
Table~\ref{table_sota} shows that taking our AANets as a plug-in architecture for $4$ state-of-the-art methods~\cite{rebuffi2017icarl,hou2019lucir,liu2020mnemonics,douillard2020podnet} consistently
improves their model performances. E.g., for CIFAR-100,  LUCIR \emph{w/} AANets and Mnemonics \emph{w/} AANets respectively gains $4.9\%$ and $3.3\%$ improvements on average.
From Table~\ref{table_sota}, we can see that our approach of using AANets achieves top performances in all settings.
Interestingly, we find that AANets can boost more performance for simpler baseline methods, e.g., iCaRL. iCaRL \emph{w/} AANets achieves mostly better results than those of LUCIR on three datasets, even though the latter method deploys various regularization techniques.

\myparagraph{Visualizing Activation Maps.} Figure~\ref{figure_gradcam} demonstrates the activation maps visualized by Grad-CAM for the final model (obtained after $5$ phases) on ImageNet-Subset ($N$=5).
The visualized samples from left to right are picked from the classes coming in the $0$-th, $3$-rd and $5$-th phases, respectively.
For the $0$-th phase samples, the model makes the prediction according to foreground regions (right) detected by the stable block and background regions (wrong) by the plastic block.
This is because, through multiple phases of full updates, the plastic block forgets the knowledge of these old samples while the stable block successfully retains it.
This situation is reversed when using that model to recognize the $5$-th phase samples. The reason is that the stable block is far less learnable than the plastic block, and may fail to adapt to new data.
For all shown samples, the model extracts features as informative as possible in two blocks. Then, it aggregates these features using the weights adapted from the balanced dataset, and thus can make a good balance of the features to achieve the best prediction.

\myparagraph{Aggregation Weights
($\alpha_{\eta}$ and $\alpha_{\phi}$).
} Figure~\ref{figure_values_plots} shows the values of $\alpha_{\eta}$ and $\alpha_{\phi}$ learned during training $10$-phase models.
Each row displays three plots for three residual levels of AANets respectively.
Comparing among columns, we can see that Level 1 tends to get larger values of $\alpha_{\phi}$, while Level 3 tends to get larger values of $\alpha_{\eta}$. This can be interpreted as lower-level residual blocks learn to stay stabler which is intuitively correct in deep models.
With respect to the learning activity of CIL models, it is to continuously transfer the learned knowledge to subsequent phases.
The features at different resolutions (levels in our case) have different transferabilities~\cite{yosinski2014transferable}.
Level 1 encodes low-level features that are more stable and shareable among all classes.
Level 3 nears the classifiers, and tends to be more plastic such as to fast to adapt to new classes.

\section{Conclusions}
\label{sec_conclusion}

We introduce a novel network architecture AANets specially for CIL.
Our main contribution lies in addressing the issue of stability-plasticity dilemma in CIL by a simple modification on plain ResNets --- applying two types of residual blocks to respectively and specifically learn stability and plasticity at each residual level, and then aggregating them as a final representation.
To achieve efficient aggregation, we adapt the level-specific and phase-specific weights in an end-to-end manner.
Our overall approach is \emph{generic} and can be easily incorporated into existing CIL methods to boost their performance.

\vspace{0.1cm}
\myparagraph{Acknowledgments.} {This research is supported by A*STAR under its AME YIRG Grant (Project No. A20E6c0101), the Singapore Ministry of Education (MOE) Academic Research Fund (AcRF) Tier 1, Alibaba Innovative Research (AIR) programme, Major Scientific Research Project of Zhejiang Lab (No. 2019DB0ZX01), and Max Planck Institute for Informatics.}

\section{Results for Different CIL Settings.}
\label{subsec_ablation_same_cls}
We provide more results on the setting with the same number of classes at all phases~\cite{rajasegaran2020itaml} in the second block (\textbf{``same \# of cls''}) of Table~\ref{table_supp_same_cls}. For example, $N$=$25$ indicates $100$ classes evenly come in $25$ phases, so $4$ new classes arrive in each phase (including the $0$-th phase).
Further, in this table, each entry represents an accuracy of the \textbf{last phase} (since all-phase accuracies are not comparable to our original setting) averaged over $3$ runs;
and ``update $\theta_{\mathrm{base}}$'' means that $\theta_{\mathrm{base}}$ is updated as $\theta_{\mathrm{base}}\gets\phi_i \odot \theta_{\mathrm{base}}$ after each phase.
All results are under ``strict memory budget'' and ``all''+``scaling'' settings, so $\phi_i$ indicate the meta-learned weights of SS operators.
The results show that
1)~``\emph{w/} AANets'' performs best in all settings and brings consistent improvements; and 2)~``update $\theta_{\mathrm{base}}$'' is helpful for CIFAR-100 but harmful for ImageNet-Subset.
\input{supplementary/tables/same_num_of_cls}

\section{Strict Memory Budget Experiments}
\label{subsec_stric_memory}
\input{supplementary/tables/sota_strict_memory}

In Table~\ref{table_sota_stric_memory},
we present the results of $4$ state-of-the-art methods \emph{w/} and
\emph{w/o} AANetss as a plug-in architecture, under the ``strict memory budget'' setting which strictly controls the total memory shared by the exemplars and the model parameters.
For example, if we incorporate AANetss to LUCIR~\cite{hou2019lucir}, we need to reduce the number of exemplars to balance the additional memory introduced by AANetss (as AANetss take around $20\%$ more parameters than the plain ResNets used in LUCIR~\cite{hou2019lucir}).
As a result, we reduce the numbers of exemplars for AANetss from $20$ to $13$, $16$ and $19$, respectively, for CIFAR-100, ImageNet-Subset, and ImageNet, in the ``strict memory budget'' setting.
For CIFAR-100, we use $530k$ additional parameters, so we need to reduce $530k \mathrm{floats}\times4\mathrm{bytes/float}\div(32\times32\times3\mathrm{bytes}/\mathrm{image})\div100\mathrm{classes}\approx6.9\mathrm{images}/\mathrm{class}$, and $\lceil 6.9 \rceil=7\mathrm{images}/\mathrm{class}$.
For ImageNet-Subset, we use $12.6M$ additional parameters, so we need to reduce $12.6M \mathrm{floats}\times4\mathrm{bytes/float}\div(224\times224\times3\mathrm{bytes}/\mathrm{image})\div100\mathrm{classes}\approx3.3\mathrm{images}/\mathrm{class}$, and $\lceil 3.3 \rceil=4\mathrm{images}/\mathrm{class}$.
For ImageNet, we use $12.6M$ additional parameters, so we need to reduce $12.6M \mathrm{floats}\times4\mathrm{bytes/float}\div(224\times224\times3\mathrm{bytes}/\mathrm{image})\div100\mathrm{classes}\approx0.3\mathrm{images}/\mathrm{class}$, and $\lceil 0.3 \rceil=1\mathrm{image}/\mathrm{class}$.
From Table~\ref{table_sota_stric_memory}, we can see that our approach of using AANetss still achieves the top performances in all CIL settings even if the ``strict memory budget'' is applied.

\section{More Ablation Results}
\label{subsec_ablation}
\input{supplementary/tables/more_ablation}
In Table~\ref{table:ablation_add}, we supplement the ablation results obtained in more settings. ``$4\times$'' denotes that we use $4$ same-type blocks at each residual level.
Comparing Row 7 to Row 2 (Row 5) shows the efficiency of using different types of blocks for representing stability and plasticity.

\section{Additional Plots}
\label{supp_sec_plots}

In Figures~\ref{figure_acc_plots}, we present the phase-wise accuracies obtained on CIFAR-100, ImageNet-Subset and ImageNet, respectively. ``Upper
Bound'' shows the results of joint training with all previous data accessible in every phase. We can observe that our method achieves the highest accuracies in almost every phase of different settings.
In Figures~\ref{figure_supp_values_cifar} and \ref{figure_supp_values_imgnetsub}, we supplement the plots for the values of $\alpha_{\eta}$ and $\alpha_{\phi}$ learned on the CIFAR-100 and ImageNet-Subset ($N$=$5$, $25$).
All curves are
smoothed with a rate of $0.8$ for a better visualization.

\section{More Visualization Results}
\label{subsec_more_visualization}
Figure~\ref{figure_gradcam_supp} below
shows the activation maps of a ``goldfinch'' sample (seen in Phase 0) in different-phase models (ImageNet-Subset, $N$=$5$).
Notice
that the plastic block gradually loses its attention on this sample (i.e., forgets it), while the stable block retains it. AANets benefit from its stable blocks.
\input{supplementary/figures/more_visualization}

\section{Source Code in PyTorch}
\label{suppsec_code}

We provide our PyTorch code on  \href{https://class-il.mpi-inf.mpg.de/}{https://class-il.mpi-inf.mpg.de/}.
To run this repository, we kindly advise you to install Python 3.6 and PyTorch 1.2.0 with Anaconda.

\input{supplementary/figures/acc_plot}
\input{supplementary/figures/3_value_plot_cifar}
\input{supplementary/figures/4_value_plot_imgnet}


\section{Introduction}
The physics of atomic clusters has gained broad importance as a domain of study of new physical objects, often termed as super-atoms, over the past few decades~\cite{jena2018super}. Such clusters are aggregates of atoms, ranging from two to a few thousand of atoms, and  provide the platform to study a new phase of matter intermediate to atoms or molecules and solids ~\cite{RevModPhys.58.533}. These systems also offer valuable opportunities to understand how the bulk properties emerge from their individual constituent atoms~\cite{stace1988small}. Possibilities of designing new class of materials with tailor made characteristics using clusters as building blocks, instead of atoms, have rendered this field diverse and vibrant~\cite{castleman2009clusters}. The extensive research on clusters over the years has revealed many unusual properties which are of interests to the field of physical and chemical sciences, material sciences, and biological sciences, making these explorations interdisciplinary in nature~\cite{dinh2013introduction}.

Many intriguing phenomena are associated with the photoresponse spectrum of atomic clusters~\cite{kreibig2013optical}. These include phenomena, such as, plasmon resonances~\cite{xia2009,madjet2009,ekardt1985}, Fano type resonances~\cite{miroshnichenko2010fano}, and diffractive modulations~\cite{jaenkaelae2011} in the photoelectron signal due to the largely well-defined cluster edges~\cite{madjet2001}. One primary focus of the present work is the photoionizing response of the plasmon resonances - resonances that form due to the collective oscillations of the valence electron cloud. Excitation of a lower energy giant surface plasmon resonance is known to be a prominent feature in the photoresponse of alkali metal clusters below the ionization threshold. This feature in metal clusters has applications in nano-optical devices, chemical and biological sensing, and bio-medicine etc~\cite{mirin2009,loo2005,liao2006}, besides its eminent role as ``spectral laboratories" to assess many-electron effects. Therefore, developing detailed and accurate understanding of the origin, the underlying dynamics, and related observable effects of this resonance are a matter of significant priority.

In addition, the presence of a higher-energy volume like plasmon makes the photo-spectrum of atomic clusters~\cite{xia2009} more robust compared to the spectrum of the corresponding bulk metal. It is well known that due to the translational invariance, the coupling of the volume oscillation in bulk to light is not feasible~\cite{ferrell1958predicted}. However, the situation is different for the case of finite systems like metal clusters due to the broken translational invariance that may cause boundary reflections of the plasma wave. As a result, optically active volume plasmon resonance can also become possible. The presence of such a volume plasmon excitation is observed for the first time in the experimental work of Xia et al.~\cite{xia2009}. Their photon depletion measurements of Na$_{20}$ and Na$_{92}$ clusters showed a broad volume plasmon resonance mounted on the decay ridge of giant surface plasmon and peaking slightly above 4 eV. The tail of this combined structure was found extended to the ionization region. The present calculation is motivated to address this spillover part of the plasmon structure closely above the ionization threshold. In studying this, we also predict a universal minimum at the decay end of the structure from an interchannel coupled phase-coherence effect arising from many-electron correlations. Furthermore, an analysis of the many-body induced self-consistent field potential at the plasmon spillover energies point to a correlation driven attractive force that suggests possible time delays of the emerging photoelectron.

The dynamical response of the Na clusters to the external electromagnetic radiation is calculated using a linear response density functional theory (DFT) approach. A fairly competent method belonging to this class is the jellium based time-dependent local density approximation (TDLDA)~\cite{onida2002electronic}. In the past, application of TDLDA to atomic clusters is found to be successful in explaining the collective phenomena occurring at low energies~\cite{madjet2009} and also the diffractive oscillations present at high energies~\cite{madjet2001}. The ease and transparency of a jellium based model enable physicists to interpret the key physics that determines the dynamics. A recent study on fullerenes indicated that a better agreement with experiment can be obtained by using one such DFT approach with Leeuwen-Baerends (LB94) exchange-correlation functional that produces the correct ground state asymptotic properties~\cite{choi2017effects}. Another study has tested the efficacy of the method for Na clusters~\cite{shaik2020}. Hence, this method is chosen in the present work.

The paper is organized in the following way. A basic description of the theoretical methodology is provided in Section II which has two parts. The details of ground state structure in the spherical jellium formalism along with a brief account of the LB94 parametrization scheme is provided in Sec. II A. A brief description of the method that incorporates electron correlations in response to the radiation is given in Sec. II B . Section III constitutes the results and discussions of: photoionization cross sections and its comparison with experiment (Sec. III A), the photoelectron oscillator strength (Sec. IIIB), the many-body correlation induced potential in the spillover region of plasmon (Sec. III C), and the prediction of a ``correlation minimum" at the high-energy end of the plasmon structure (Sec. III D). Finally, Section IV concludes the study.

\section{Theoretical Methodology}

\subsection{Ground states of Na$_{20}$ and Na$_{92}$}

To investigate the ground state electronic structure of Na$_{n}$($n$=20 and 92), a DFT approach is adopted. For such systems of closed-shell configurations the jellium model serves as a very good approximation. In this model, the jellium potential ($V_{jel}(\bm r)$) replaces the ionic core of 20 and 92 Na$^{10+}$ ions, respectively for Na$_{20}$ and Na$_{92}$, by potentials constructed after homogeneously smearing their positive charges into jellium spheres. The radius of each cluster is determined by the number $n$ of ions present in the system. The radius of Na$_{20}$ is calculated to be 10.67 a.u.\ and that of Na$_{92}$ to be 17.74 a.u.\ using the formula  $R_{c}$=$r_{s}n^{1/3}$, where $r_{s}$ is the Wigner-Seitz radius (3.93 a.u.) of a Na atom. The Kohn-Sham equations for $n$ delocalized valence electrons, a 3s$^1$ electron from each Na atom, are solved to obtain the ground state structures of Na$_{20}$ and Na$_{92}$. It is to be noted that to match the valence ionization thresholds with the experimental values~\cite{Chandezon}, suitable constant pseudo potentials are added.
\begin{figure}[ht]
\includegraphics[width=\linewidth]{ground.eps}
\caption{a) Ground-state radial wavefunctions for HOMO and HOMO-1 levels calculated for Na$_{20}$ and Na$_{92}$. b) Ground-state radial potentials and energy levels calculated for Na$_{20}$ and Na$_{92}$.}
\label{gpot}
\end{figure}

In terms of the single-particle density $\rho(\bm r)$, the ground state self-consistent field LDA potential reads as,
\begin{equation}\label{lda-pot}
V_{\scriptsize \mbox{LDA}}(\mathbf{r}) = V_{\mbox{\scriptsize jel}}(\mathbf{r}) + \int d\mathbf{r}'\frac{\rho(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|} + V_{\scriptsize \mbox{xc}}[\rho(\mathbf{r})],
\end{equation}
where the second and third terms on the RHS are respectively the direct and exchange-correlation (xc) potentials. Approximate form for $V_{\scriptsize \mbox{xc}}$ has to be used since its exact form is unknown. In the present work, we approximate the $V_{\scriptsize \mbox{xc}}$ with the LB94 functional which provides an accurate asymptotic description of the ground-state properties. This functional belongs to a gradient corrected class of approximation by van Leeuwen and Baerends, parameterized in terms of the reduced density and its gradient $\nabla\rho(\bm r)$~\cite{van1994}, to be given by,\\
\begin{equation}\label{lb94}
V_{\mbox{\scriptsize LB}} = -\beta{[\rho({\bf r})]^{1/3}}\frac{(\xi X)^2}{1+3\beta\xi X \sinh^{-1}(\xi X)},
\end{equation}
where $\beta$=0.01 is a fitting parameter and $X = [\nabla\rho]/\rho^{4/3}$. The parameter $\xi$ arises due to the change of system from spin-polarized to spin-unpolarized form ~\cite{oliver79}. This scheme of gradient correction of density is derived from a gradient expansion series that naturally eliminates self-interactions and produces the correct asymptotic behavior.

\begin{table}[ht]
\caption{Binding energies (BE) of HOMO and HOMO-1 levels of Na$_{20}$ and Na$_{92}$ in the harmonic oscillator notation.}
\label{tab1}
\begin{tabular}{c|c cccc}
\hline\noalign{\smallskip}
    &    & BE$_{\mbox{\scriptsize HOMO}}$ (eV) &    & {\smallskip} BE$_{\mbox{\scriptsize HOMO-1}}$ (eV)\\
\noalign{\smallskip}\hline\noalign{\smallskip}
Na$_{20}$  &  & $2s$ {\smallskip} $-3.75$ & &$1d$ {\smallskip} $-4.20$      \\
Na$_{92}$ &     & $1h$ {\smallskip}  $-3.47$ & & $3s$  {\smallskip} $-3.57$  \\
\noalign{\smallskip}\hline
\end{tabular}
\end{table}

The calculated binding energies of the highest occupied molecular orbitals (HOMO) and (HOMO-1) of these systems,  of 2s and 1d character respectively for Na$_{20}$ and 1h and 3s respectively for Na$_{92}$, are given in Table~\ref{tab1} and their radial wavefunctions are shown in Fig.~\ref{gpot}(a). It can be seen that HOMO level of  Na$_{20}$ has a single node whereas Na$_{92}$ is nodeless. The situation reverses for HOMO-1. As a result, the radial structures of HOMO and HOMO-1 wavefunctions significantly differ in these two systems. Fig.~\ref{gpot}(b) shows the ground-state radial potentials of Na$_{20}$ and Na$_{92}$ clusters and the level energies.

\subsection{Dynamical Response}
The dynamical response of the cluster subjected to the external dipole field `$z$' can be calculated by employing a linear-response DFT scheme, known as time-dependent LDA (TDLDA) modified by LB94~\cite{madjet2008photoionization}. The dipole field can induce a frequency-dependent complex change in the electron density $\delta \rho$~\cite{petersilka1996}, which, within a linear-response perturbative framework, can be written as:
\begin{equation}\label{cor-dens}
\delta \rho (\mathbf{r}^{\prime}; \omega) = \int \chi (\mathbf{r}, \mathbf{r}^{\prime}; \omega)
z  d\mathbf{r},
\end{equation}
where $\chi$ is the full susceptibility of the system which incorporates all the dynamical electron correlations. The above equation can be re-written in terms of independent particle (IP) susceptibility, $\chi_0$, and the complex field $\delta V$ as,
\begin{equation}\label{cor-dens-1}
\delta \rho (\mathbf{r}^{\prime}; \omega) = \int \chi_0 (\mathbf{r}, \mathbf{r}^{\prime}; \omega)
\delta V(\mathbf{r}; \omega)  d\mathbf{r},
\end{equation}
where $\chi_0$ is constructed by the ground-state single-electron orbitals $\phi_{nl}$ and energies $\epsilon_{nl}$~\cite{feibelman1975microscopic}. The IP susceptibility $\chi_0$ is related to $\chi$ by the matrix equation:
\begin{eqnarray}\label{chi-chi0}
\chi = \chi_0[1-(\partial V/\partial \rho)\chi_{0}]^{-1}.
\end{eqnarray}
The total field, $\delta V$, is given by the following relation:
\begin{equation}\label{comp-pot}
\delta V(\mathbf{r}; \omega)=z+V_{ind}(\mathbf{r}; \omega),
\end{equation}
where
\begin{equation}\label{V_ind}
V_{ind}(\mathbf{r}; \omega) = \int \frac{\delta \rho (\mathbf{r}^{\prime};
\omega)}{|\mathbf{r}-\mathbf{r}^{\prime}|}d\mathbf{r}^{\prime}+\Big[\frac{\partial V_{XC}}{\partial \rho}\Big]_{\rho = \rho_{0}} \delta \rho(\mathbf{r}; \omega).
\end{equation}

Using the matrix inversion method ~\cite{BERTSCH1990}, Eq.~\ref{chi-chi0} is finally solved for $\chi$ which in turn is used for obtaining $\delta \rho$ and thereby $\delta V$ by making use of Eq.~\ref{cor-dens} and Eq.~\ref{comp-pot} in a self-consistent way.

In this formalism, the photoionization cross section corresponding to a bound-to-continuum dipole transition $n\ell\rightarrow k\ell^\prime$ is then calculated using
\begin{equation}\label{cross-pi}
\sigma_{n\ell\rightarrow k\ell'} \sim |\langle k\ell'|\mbox{$\delta V$}|n\ell\rangle|^2.
\end{equation}
It is clear from the above equation that in addition to the external perturbation $z$, the calculation involves the complex induced field $V_{ind}$ driven by the many-electron correlations. Obviously, setting $\delta V=z$ yields the independent particle LDA cross section  that ignores correlations. A comparison of LDA and TDLDA, therefore, easily possible in this method to study the role of many-electron effects in the photoionization process.
\begin{figure}[ht]
\includegraphics[width=\linewidth]{fig1.eps}
\caption{TDLDA and LDA cross sections of Na$_{20}$ (a) and Na$_{92}$ (b).}
\label{lda-tdlda}
\end{figure}

\section{Results and Discussion}

\subsection{Photoionization cross-sections and comparison with measurements}

Total photoionization cross-section of Na$_{20}$ and Na$_{92}$ clusters, calculated in TDLDA+LB94 and the corresponding LDA results, are shown in Fig.~\ref{lda-tdlda}(a) and Fig.~\ref{lda-tdlda}(b) respectively. As seen, the TDLDA cross section profiles, as a function of photon energies, strongly differ from the LDA profiles by (i) significant enhancements near the ionization threshold region, (ii) the presence of autoionization resonances, and (iii) the appearance of a minimum at about 7.20 eV for Na$_{20}$ and 6.22 eV for Na$_{92}$. The near-threshold enhancement in TDLDA cross sections combines the remnant of the giant surface plasmon resonance, that flourishes below the ionization threshold, and a major portion of the broader volume plasmon resonance, both emerging from the collective electronic motions when subjected to an external electromagnetic field. %
This resonance spillover from the discrete to continuum spectrum is embedded by a host of narrow spikes, which are the autoionization resonances resulting from the degeneracy of ionization channel with the inner-level single-electron discrete excitations~\cite{fano1961effects}. As seen, these features are completely missing in the LDA predictions which neglect the electron correlations. Structures in the LDA profiles in this near-threshold region are due to the gradual openings of inner-level ionization channels.
\begin{figure}[ht]
\includegraphics[width=\linewidth]{fig2.eps}
\caption{TDLDA cross sections along with experimental data from ~\cite{xia2009} for Na$_{20}$ (a) and Na$_{92}$ (b). Lorentzian fits to the experimental cross sections are also shown.}
\label{tdlda-expt}
\end{figure}

The existence of the volume plasmon resonance and the extended tail of the surface resonance to the ionization part is experimentally seen in the work of~\cite{xia2009}. Comparisons of our calculations with these measured data is shown in Fig.~\ref{tdlda-expt}(a) for Na$_{20}$ and Fig.~\ref{tdlda-expt}(b) for Na$_{92}$. The data covers the range of 2 eV to 5.64 eV. Following Ref.\ \cite{xia2009}, fits to the experimental data as a sum of four Lorentzian profiles for Na$_{20}$ and a sum of three for Na$_{92}$ are also shown in Fig.~\ref{tdlda-expt}. The lower energy Lorentzians in both cases represent the giant portion of the plasmon below the first ionization threshold that our photoionization results can not access. The present TDLDA calculations span the ionization part of the spectrum only and, as noted, overlap very well with the higher energy Lorentizian fits used for both the clusters. The peak energy position (E$_0$) and width ($\Gamma$) of the fourth Lorentizian curve for the higher energy plasmon resonance (HEP) for Na$_{20}$ and the corresponding position and width for the third Lorentzian curve for Na$_{92}$ are shown in Table~\ref{tab2}. It is clear from the figure that the steady background parts of TDLDA cross sections agree well with the experiment describing the spillover of the plasmon resonances. As the size of the cluster increases the peak positions of the volume plasmon shifts to higher energies (blue shifted). Curiously, this size effect is found to be opposite in the case of fullerenes (C$_{60}$ and C$_{240}$) where peak positions shifted to lower energies as the size increased~\cite{choi2017effects}. It is to be noted that the narrow autoionization resonances are missing in the experimental data. This is due to the finite temperature effects of the metal clusters in experimental conditions which lead to the coupling of electronic motion with the temperature induced vibrational and rotational modes~\cite{bertsch1989thermal}, as was also found and discussed earlier for C$_{60}$~\cite{madjet2008photoionization}.

\begin{table}
\caption{Resonance positions ($E_o$), FWHM ($\Gamma$) of the higher energy plasmon (HEP)~\cite{xia2009}.}
\label{tab2}
\begin{tabular}{c  c  c}
\hline\noalign{\smallskip}
              & $E_o$ (eV)   & $\Gamma$ (eV) \\
\noalign{\smallskip}\hline\noalign{\smallskip}

Na$_{20}$ HEP  & 4.04  & 1.19     \\
\noalign{\smallskip}\hline\noalign{\smallskip}
Na$_{92}$ HEP & 4.20  & 1.16   \\
\noalign{\smallskip}\hline
\end{tabular}
\end{table}

\begin{figure}[ht]
\includegraphics[width=\linewidth]{fig3.eps}
\caption{Accumulative oscillator strengths for Na$_{20}$ (a) and Na$_{92}$ (b) where the theoretical results are calculated using the present TDLDA+LB94 cross sections and the experimental OS are calculated using fitted values to the data in~\cite{xia2009}.}
\label{OS}
\end{figure}

\subsection{Oscillator strength}

From the experimental data shown in Fig.~\ref{tdlda-expt}, it is quite clear that the bulk of the absorption oscillator strength (OS) is consumed in the excitation part of the spectrum below the ionization threshold. In order to get a quantitative measure of the fraction of the OS exhausted in the ionization process, we introduce an ``accumulative" OS as a function of photon energy given by the following relation:
\begin{equation}\label{OSt}
OS_{t}(E_{f})=OS_{b}+\int_{E_{th}}^{E_{f}} \sigma(E)dE,
\end{equation}
where OS$_{b}$ is defined as the baseline OS that is calculated by integrating the experimental cross section~\cite{xia2009} from the starting photon energy of the measurement to the theoretical ionization thresholds (E$_{th}$). OS$_{b}$ therefore corresponds to the oscillator strength exhausted by the plasmon resonances below the threshold. The second term in Eq.~\ref{OSt} corresponds to the incremental addition to the OS$_{b}$, cumulatively with energy, due to the plasmon spill to the ionization region from E$_{th}$ to E$_{f}$. Here E$_{f}$ will vary from E$_{th}$ onwards allowing us to calculate OS$_{t}$ as a function of E$_{f}$. Note that until 5.64 eV, the photon energy corresponding to the last experimental datum, E$_{f}$'s are the various energies for which experimental measurements are available. In Eq.~\ref{OSt}, $\sigma(E)$ corresponds to TDLDA total cross section for the calculations of theoretical OS$_{t}$. For calculating the experimental OS$_{t}$, on the other hand, fitted cross section data between two successive measurements are used.

Figure~\ref{OS} shows OS on a zero-to-one scale, because it represents OS per atom, that is the total OS divided by $n$. For Na$_{20}$, the OS$_{b}$=0.5230 and for Na$_{92}$, OS$_{b}$=0.4352. These suggest that about 52\% of OS is exhausted below the ionization threshold of Na$_{20}$, while it is about 44\% for Na$_{92}$. The implication is that more fraction of electrons are available to participate in the ionization process for the larger cluster. Hence, the size of plasmon spill increases for  Na$_{92}$. The rise of the theoretical OS curve above the experimental curve is due primarily to the narrow autoionizing resonances in the theoretical spectrum which add strengths. This accounts roughly about 20\% of OS for Na$_{20}$, while it is 30\% for Na$_{92}$. The experimental measurements cannot access such narrow single particle resonances due to the finite temperature effects~\cite{bertsch1989thermal}; similar distinction was also noted for the plasmon resonances in C$_{60}$~\cite{madjet2008photoionization}. However, the missing part of the oscillator strengths in the measurements may be distributed over higher photon energies, in which case it is expected that the two curves may eventually intersect at higher energies. One comment should be made about the measurements though. The energy deposited into the cluster will not only dissipate through ionization that we calculate, but also be shared with thermalization and evaporation channels as well. Even though the corresponding branching ratio is unknown, the ionization is known to be a much faster process (in attoseconds) compared to thermalization/evaporation (in femtoseconds). Therefore, one may reasonably assume that the main contribution to the measured data above the threshold is likely from the ionization process.
\begin{figure*}[!ht]
\includegraphics[width=16cm]{scf.eps}
\caption{The real and imaginary part of radial self-consistent field potential, $V_{ind}(r)$, for Na$_{20}$ (left panels) and Na$_{92}$ (right panels).}
\label{scf}
\end{figure*}

\subsection{Collectivity induced field}

TDLDA dipole matrix element $ \langle k\ell'|z+V_{ind}|n\ell \rangle$ in \eq{cross-pi} requires knowledge of $V_{ind}$, which is the complex induced field driven by electron correlations and is singularly responsible for the plasmonic enhancement in the cross section. The behavior of the real and imaginary parts of $V_{ind}$ across the collective resonance region is well known~\cite{wendin73,zangwill80}. Im($V_{ind}$) characteristically shows a well-type shape across the energy range of the resonance where the minimum of the well occurs near the energy of the resonance peak. On the other hand, Re($V_{ind}$) executes an oscillation by switching the sign over this range where it sluices through the zero at the resonance peak. These two distinct behaviors can be combined in a unified picture. Im($V_{ind}$) has a predominant collective character, while Re($V_{ind}$) represents effects of the external field. Therefore, as the resonance builds with increasing energy and approaches its peak, the effect of external field reduces while the collective motion grows. At the peak, the effect of external field is negligible, where the collective response dominates. At the waning part of the resonance this trend reverses.

In the current study, however, we access the remnant of the plasmon structure spilling over to the photoionization channel which is fragmented by many single electron narrower resonances mixing coherently with the plasmon effect. This coherence indicates that there are interferences between the single electron Auger and collective plasmon. So the simple mechanism expressed above becomes complex, the induced field becomes structured, but an interesting general trend should remains that we now explain.

Figure~\ref{scf} shows the 3D plots of real and imaginary part of radial induced field, $V_{ind}(r; \omega)$, for Na$_{20}$ and Na$_{92}$. All the plots show many  smaller structures, which are because of single-electron resonances that become relatively stronger as the plasmon spill gradually weakens with increasing energy.  However, there appears wider well-shapes for Im ($V_{ind}$) of both systems. This broader range of negative fields suggests the emergence of an attractive field as a result of the electrons' collective dynamics in these plasmonic energies. Consequently, the liberating photoelectrons will experience a resistance against their exit and will likely slow down. This can lead to observable time-delayed emissions of electrons~\cite{schultze} in the plasmon spillover range. Particularly, this may find relevance in the context of the electron's intrinsic Wigner-type delay properties~\cite{magrakvelidze2015attosecond}. Time delay studies in photoemission belong to a contemporary field of interest based on methods of attosecond photoemission measurements in RABITT or streaking spectrometry~\cite{Pazourek}. Similar attractive force, driven by electronic collective interactions, has earlier been predicted at plasmon photoionization energies of C$_{60}$~\cite{magrakvelidze2020}.
\begin{figure*}[ht]
\includegraphics[width=\linewidth]{fig4.eps}
\caption{Subshell cross-sections for Na$_{20}$, (a) and (c), and Na$_{92}$, (b) and (d), compared with corresponding total cross-section obtained using LDA (left-panels) and TDLDA (right panels).}
\label{subcross}
\end{figure*}

\subsection{Correlation minimum}

As seen in Fig.~\ref{lda-tdlda}, the total photoionization cross sections undergo oscillations both in LDA and TDLDA starting from energies where the plasmon structure begins to fizzle. The oscillations far above the plasmon region are well understood and are attributed to the diffraction of the photoelectron waves from the edges of the cluster~\cite{madjet2001}. Since the feature is associated with the cluster geometry, it is seen even in the independent particle LDA results. The merging of the oscillations at higher energies in TDLDA and LDA reflects on the fact that the oscillations have nothing to do with the electron correlations. However, the TDLDA cross sections reveal a rather strong additional minimum occurring at 7.20 eV for Na$_{20}$ and 6.22 eV for Na$_{92}$ which are clearly missing in their LDA counterparts. Incidentally, this feature was also revealed in earlier calculations~\cite{ekardt1985}, but was never properly interpreted. If it were a diffractive feature, it would be present in the LDA results as well. Obviously, that is not being the case suggests that this lower energy minimum must be induced by an electron correlation effect.

Further insights into this minimum can be harnessed by looking at the individual subshell cross sections. In Fig.~\ref{subcross} are shown the cross sections from the three outermost subshells of Na$_{20}$ and Na$_{92}$ along with the total cross sections obtained from LDA (Fig.~\ref{subcross}(a) and Fig.~\ref{subcross}(b)) and TDLDA(Fig.~\ref{subcross}(c) and Fig.~\ref{subcross}(d)). It is clear from the TDLDA results that in each of the subshell cross sections for a cluster this particular minimum in consideration occurs dramatically at the same low photon energy to create a rather emphatic minimum in the total cross section. But this feature is totally absent in the LDA curves.

For a subshell, while the creation of a minimum derives directly from the ionization amplitude, its position in energy traces in to the ionization phase. The positions of the minima at high energies are seen to be very sensitive to the ionizing subshells. It was shown earlier that the energy-positions of high energy minima are a function of the phase of the individual subshell LDA transition amplitudes~\cite{frank97}, since the contribution of $\delta V$ [\eq{comp-pot}] is virtually zero at these energies. In addition, at very high energies the contributions to the scattering-phase from even the short-range and Coulomb potentials are negligible. As a result, the final state wavefunction is approximately of the form $\psi_f \approx \cos(kr-\frac{l^{\prime} \pi}{2})$, where $l^\prime = l \pm 1$ being the dipole selected final state angular momentum, `$l$' is the initial state and `$k$' is the photoelectron momentum. Consequently, the squared transition amplitudes differ by a phase of $l^\prime \pi$. This explains why the oscillations in the cross sections of two subshells of angular momenta differing by an odd integer, such as, s and p, p and d, s and h etc.\ are roughly out-of-phase, while those differing by an even integer, such as, s and d, p and f etc.\ are roughly in-phase~\cite{mccune2008unique}. This can be seen in both LDA and TDLDA subshell results at very high energies. At not-so-high energies, these patterns are not seen to be exactly followed due to non-negligible short-range and Coulomb phases.

However, the energy-concurrent minima in the TDLDA profiles occur at a low enough energy. Even though the diffractive mechanism forming oscillations is emerging, this range is still in the waning part of the plasmon structure. Hence, the contribution from the correlation phase, that is the phase of the amplitude with significant strength of $\delta V$, still dominates over the $l^\prime \pi$ in this range. The correlation phase, being collective in nature, is subshell independent, leading to all subshells experiencing effectively the same phase~\cite{magrakvelidze2015attosecond}. This ensures that the minimum for all the subshells to appear coherently at the same energy. Since this coherence originates from a direct influence of the collective dynamics, we call this feature a correlation minimum, which can, in principle, be accessible by the photoelectron spectroscopy (PES).

\section{Conclusions}
In summary, the present jellium based linear response DFT calculation describes the remnant of the giant surface and the bulk of the volume plasmon in the photoionization cross section of Na$_{n}$ ($n$=20 and 92) clusters.  Results show that an appreciable amount of plasmon spillover into the ionization continuum occurs in these systems. The steady background part of the cross section exhibits reasonably good agreement with previous experimental results of the absorption of these clusters. However, a detailed scrutiny of the  oscillator strength calculated from the cross section and its comparison with the strength extracted from the measurements enables us to quantify the contribution of the single electron Auger-type resonances that the theory predicts. Such, rather narrow, resonances shown in the theory, which does not incorporate the temperature and vibro-rotational effects, are missing in the measurements that include these effects. Furthermore, a deeper scrutiny of the many-body induced self-consistent DFT field reveals the presence of an attractive force in the plasmon spillover energy region. We speculate that this force can cause an observable delay in the emission of photoelectrons at these energies, enabling these clusters as interesting candidates for time-delay measurements. The current study further uncovers the presence of a correlation minimum in the cross sections which appears at the waning range of the plasmon resonance structure and which may also attract experimental effort. To this end, the collective-motion spectral properties divulged in the current DFT study of Na clusters are likely quite general and may apply at various degrees of prominence in other metallic clusters.
\begin{acknowledgments}
The research is supported by the SERB, India, Grant No.\ EMR/2016/002695 (HRV) and by the US National Science Foundation Grant No.\ PHY-1806206 (HSC).
\end{acknowledgments}


\section{Introduction}

In hypernuclear physics the lambda-nucleon-nucleon ($\Lambda N N$) three-body problem plays a very important role, akin to that of the NN and NNN systems in the non-strange sector of nuclear physics. In particular, the spin-isospin channel $\Lambda n p$ ($J^{\pi}=1/2^+$) is traditionally used for constraining lambda-nucleon forces. After a bound state (or resonance?) was inferred for $\Lambda n n$ ($J^{\pi}=1/2^+$) \cite{rap2013}, these two systems together now serve to provide more constraints on lambda-nucleon forces. The lambda hypertriton ($\Lambda n p$) is important because it can be used to constrain both the $\Lambda p$ and $\Lambda n$ forces, while $\Lambda n n$ is crucial for constraining the $\Lambda n$ force. It is therefore of uppermost importance that the same set of $\Lambda p$ and $\Lambda n$ forces, with charge symmetry breaking, is used to study both three-body problems. Through this approach, the following important \textit{double} question can be answered: what $\Lambda n p$ binding energy is consistent with a $\Lambda n n$ bound state or resonance?

More than half a century after the first observation of a hypernucleus in nuclear emulsions, there is still no published data for $\Lambda n$ scattering. This is because there are no lambda hyperon targets and, in the case of inverse kinematics, lambda hyperon beams are very difficult to prepare. The lack of $\Lambda n$ data is a serious setback for all methods used in the development of $\Lambda N$ potentials: Meson-Exchange theory, Chiral Effective Field theory, quantum chromodynamics and, recently, inverse scattering theory. As a result, basic structural properties like binding energy, size and lifetime of the simplest hypernucleus cannot be \textit{accurately} computed.

Since 1958, a number of calculations on $\Lambda n n$ predicted that a bound state was not possible \cite{dal1958a, dal1958b, dow1959, gar1987,bel2008}. After evidence for a possible bound state was inferred from heavy-ion experiments by the HypHI Collaboration at GSI, a number of recent theoretical studies based on improved $\Lambda N$ forces from Meson-Exchange theory, Chiral Effective Field theory and quantum chromodynamics have also argued against the existence of such a bound state \cite{gal2014, gar2014, hiy2014, afn2015, kam2016, fil2016, gib2019}. In order to provide a new perspective on the $\Lambda N$ force, inverse scattering theory has recently been applied on theoretical $\Lambda p$ and $\Lambda n$ scattering data \citep{meo2019}. Spin averages of the recovered $\Lambda p$ and $\Lambda n$ potentials (GLM-YN0) were tested in \cite{meo2020} by computing the binding energy and root-mean-square radius of $\Lambda n p$ ($J^{\pi}=1/2^+$), using Malfliet-Tjon V potential for the NN interaction. The aim of this paper is to use GLM-YN0 potentials with the semi-realistic spin-dependent Malfliet-Tjon I/III NN potential that correctly binds the deuteron, in exploring the $\Lambda n p$ ($J^{\pi}=1/2^+$) and $\Lambda n n$ ($J^{\pi}=1/2^+$) three-body problems.

\section{Lambda-neutron and lambda-proton interactions: GLM-YN0 potentials}

The charge-asymmetric pair of $\Lambda p$ and $\Lambda n$ potentials used in this paper were developed through Gel'fand-Levitan-Marchenko theory \citep{meo2019}. These potentials were restored through the application of this method on sub-threshold theoretical $s$-wave scattering phases. Based on the spin multiplicities of the $^1\text{S}_0$ and $^3\text{S}_1$ $\Lambda N$ states, effective potentials are constructed from one-quarter spin singlet and three-quarters spin triplet potentials from inverse scattering theory \cite{meo2020}:

\begin{align}
\label{eq:spinlp}
V_{\Lambda p} = \frac{1}{4}V_{\Lambda p}(^1\text{S}_0) + \frac{3}{4}V_{\Lambda p}(^3\text{S}_1)
\end{align}

\begin{align}
\label{eq:spinln}
V_{\Lambda n} = \frac{1}{4}V_{\Lambda n}(^1\text{S}_0) + \frac{3}{4}V_{\Lambda n}(^3\text{S}_1)
\end{align}
The data obtained from these spin-averages, which will henceforth be referred to as GLM-YN0 potentials, was fitted with the following three-range Gaussians:
\begin{align}
\label{eq:three_gaussians}
V_{\Lambda N}(r) = \sum_{i=1}^{3} V_{i} \exp \left \{ \frac{-(r - \mu_i)^2}{\sigma^2_i} \right \}.
\end{align}

Figure \ref{fig:both_effective1} illustrates the features of GLM-YN0 potentials, especially a noticeable charge symmetry breaking. The theoretical scattering phases used for inverse scattering theory in \cite{meo2019} carry in them not only this feature of charge asymmetry, but also $\Lambda - \Sigma$ coupling. Therefore, it can be stated that even though GLM-YN0 potentials only explicitly have a central term, the effects of charge symmetry breaking and $\Lambda - \Sigma$ coupling are \textit{built} into this central term, since phase equivalence with NSC1997f was verified in \cite{meo2019}. The parameters $V_{i}, \mu_i$ and $\sigma_i$ in Equation \eqref{eq:three_gaussians} were determined through a nonlinear least squares fit. Estimated values are reported in Table \ref{fit}.

\begin{figure}[!h]
\centering
\includegraphics[width=1.0\linewidth]{both_effective1.eps}
\caption[]{GLM-YN0 effective $\Lambda p$ and $\Lambda n$ potentials: these potentials are constructed from one-quarter of the spin singlet state and three-quarters of the triplet state potential from inverse scattering theory.}
\label{fig:both_effective1}
\end{figure}

\begin{table}[h!]
\centering
\caption[]{Parameters for GLM-YN0 effective $\Lambda p$  and $\Lambda n$ potentials.}
\begin{tabular}{p{1.0cm}p{3.0cm}p{3.5cm}p{3.6cm}}
           &                  &   $\Lambda p$             &                               \\  \hline
           &  $V_i$ /MeV      &   $\mu_i$/ fm             &  $\sigma_i$/fm                \\ \hline
$i=1$      & $45.88      $    &   $0.1148 \pm 0.0006601$  &  $-0.3932 \pm 0.0008502$      \\
$i=2$      & $8.106e+07  $    &   $-1.193 \pm 0.001948$   &  $0.3575 \pm 0.0005306$       \\
$i=3$      & $-47.04     $    &   $0.3748 \pm 0.0001386$  &  $0.1667 \pm 0.0002179$       \\   \hline
           &                  &                           &                               \\
           &                  &                           &                               \\
           &                  &    $\Lambda n$            &                               \\ \hline
           &  $V_i$ /MeV      &   $\mu_i$/ fm             &  $\sigma_i$/fm                \\ \hline
$i=1$      & $186.9      $    &   $- 0.3476 \pm 0.001364$ &  $-0.5469 \pm 0.001125$       \\
$i=2$      & $6.74e+04   $    &   $- 0.383 \pm 0.001433$  &  $0.191 \pm 0.0005638$        \\
$i=3$      & $- 52.14    $    &   $0.3243 \pm 0.0001977$  &  $0.2013 \pm 0.0002466$       \\   \hline
\end{tabular}
\label{fit}
\end{table}

\section{Neutron-proton and neutron-neutron interactions: Malfliet-Tjon I/III potential}

Malfliet-Tjon I/III potential is used to represent the NN interaction. In the spin-triplet states, the $(S,I)=(1,0)$ potential is

\begin{align}
\label{eq:mtiii_np}
V_{np}(r) =  [V^t_A \text{exp}(-\mu_a r ) + V^t_R \text{exp}(-\mu_r r )]/r
\end{align}

In the spin-singlet state, $(S,I)=(0,1)$, the potential is give by

\begin{align}
\label{eq:mti_nn}
V_{nn}(r) =  [V^s_A \text{exp}(-\mu_a r ) + V^s_R \text{exp}(-\mu_r r )]/r
\end{align}

The parameters used for these NN potentials, shown in Table \ref{mt_iii_i}, are taken from \cite{fri1990}. These parameters were established so that the Malfliet-Tjon I/III potential can be used as a tool for benchmarking few-body computations in nuclear physics. With these parameters, the binding energy of the deuteron was found to be -2.2307 MeV.

\begin{table}[h!]
\centering
\caption[]{Parameters for Malfliet-Tjon I/III NN potentials, from \cite{fri1990}.}
\begin{tabular}{p{3.0cm}p{3.0cm}p{3.0cm}p{3.0cm}}
$V^t_A$    & $V^t_R$          &   $\mu_a$                 &  $\mu_r $                \\ \hline
-626.885   & 1438.72          &   1.55                    &  3.11                    \\
$V^t_A$    & $V^t_R$          &                           &                          \\ \hline
-513.968   & 1438.72          &                           &                           \\   \hline
\end{tabular}
\label{mt_iii_i}
\end{table}

\section{Hyperspherical method}

Let $\{\vec{x}_i, \vec{y}_i \}$ be Jacobi coordinates, after elimination of centre-of-mass motion \cite{tho2004}. These coordinates are transformed into six-dimensional hyperspherical coordinates $\{\rho, \Omega_5\}$, where $\rho = \arctan (y_i/x_i)$ is the hyperradius and $\Omega_5$ represents one hyperangle ($\theta_i$) and two sets of polar angles ($\nu_{x_i}, \nu_{y_i}, \omega_{x_i}, \omega_{y_i}$). Faddeev amplitudes are given by
\begin{align}
\label{eq:faddeev}
 \psi_{\alpha_i}^{i,J}(x_i, y_i) =  \sum_{K_i=K_{min}}^{K_{max}} \rho^{-5/2} \chi_{\alpha_i, K_i}^{i,J}(\rho) \phi_{K_i}^{\ell_{xi}, \ell_{yi}}(\theta_i),
\end{align}
where $K_i$ are hyperangular momenta and $\alpha_i = \{ (\ell_{xi}, \ell_{yi})L_i, (s_j,s_k)S_{xi} \}J_i;s_i$ is the coupling order for Faddeev amplitudes. The $\Lambda NN $ three-body problem is simulated through the following system of coupled hyperradial equations \cite{smi1977}:

\begin{align}
\label{eq:hyperradial}
 \left \{ -\frac{\hbar^2}{2m}\frac{d^2}{d \rho^2} + \frac{\hbar^2}{2 m \rho^2}\mathcal{L}_{K_i}(\mathcal{L}_{K_i} +1)  - E  \right \} \chi_{\alpha_i, K_i}^{i,J}(\rho) =- \sum_{ j \alpha_j K_j } V_{\alpha_i K_i, \alpha_j K_j}^{ij} (\rho)       \chi_{\alpha_j, K_j}^{j,J}(\rho),
\end{align}
where $m$ is a unit mass (taken here to be the mass of a nucleon, 939.0 MeV), $\mathcal{L}_{K_i}=K_i+3/2$ and $V_{\alpha_i K_i, \alpha_j K_j}^{ij}$ are coupling potentials. This coupled system of hyperradial equations is converted to the following eigenvalue problem \cite{tho2004}:
\begin{align}
\label{eq:eigenvalue}
H \mathbf{a}=E\mathbf{a}
\end{align}
where $\mathbf{a}$ is the eigenvector and $E$ the eigenvalue.

\section{Results and discussion}

\subsection{$\Lambda n p$ ($J^{\pi}=1/2^+$)}

The eigenvalue problem in Equation \ref{eq:eigenvalue} was solved using the computer code \texttt{FaCE} \cite{tho2004}. For $\Lambda n p$ ($J^{\pi}=1/2^+$), the $np$ potential in Equation \ref{eq:mtiii_np} was used, together with GLM-YN0 $\Lambda p$ and $\Lambda n$ potentials presented in this paper. The binding energy and rms matter radius computed are found to be (-3.0759 MeV, 7.7 fm). This gives a lambda separation energy of $B_{\Lambda}=0.85$ MeV, relative to the deuteron binding energy for the $np$ potential used in the calculation (-2.2307 MeV). In a recent experiment it was found that $B_{\Lambda}=0.41$ MeV \cite{ada2020a}, which is a consequential increase from the widely used value $B_{\Lambda}=0.13$ MeV \cite{jur1973}. A more accurate experiment, that is \textit{expected} to measure an even higher $\Lambda n p$ binding energy, is in the planning phase at JLab \cite{gog2020}. A higher lambda separation energy in $\Lambda n p$ has many consequences in hypernuclear physics, some of which are discussed within chiral effective field theory in \cite{le2020} for $B_{\Lambda}=0.23$ to 0.51 MeV.

The rms radius of $\Lambda n p$ is another property that is dependent on the $\Lambda n p$ wavefunction, even though it has never been measured in any experiment. As expected, the value computed here is smaller than that from models in which $\Lambda n p$ has a lower binding energy. For example, a $\Lambda$-deuteron rms radius of 10.3 fm was recently computed in pionless effective field theory \cite{hil2019} while a value of 10.6 fm is reported in \cite{bra2019}.

A recent high-precision measurement by the STAR Collaboration at RHIC-BNL \cite{ada2019} found the lifetime of $\Lambda n p$ ($J^{\pi}=1/2^+$) to be about 50\% shorter than that of a free $\Lambda$ hyperon: this clearly implies that $\Lambda p$ and $\Lambda n$ forces are much stronger than currently used forces that compute lambda hypertriton lifetimes nearly equal to that of a free $\Lambda$ hyperon. Our potentials are therefore in accord with the STAR Collaboration lifetime measurements that challenge the notion of $\Lambda n p$ being a loose $\Lambda$ hyperon around a deuteron core. Therefore, an important test for GLM-YN0 $\Lambda p$ and $\Lambda n$ potentials lies in the computation of the lifetime of $\Lambda n p$ ($J^{\pi}=1/2^+$).

\subsection{$\Lambda n n$ ($J^{\pi}=1/2^+$)}

For $\Lambda n n$ ($J^{\pi}=1/2^+$), the $nn$ potential in Equation \ref{eq:mti_nn} was used, together with GLM-YN0 $\Lambda n$ potential. By examining all the computed eigenvalues of Equation \ref{eq:eigenvalue}, two possible sets of values for the binding energy and rms matter radius are found: (-0.0695 MeV, 12.3 fm) and (-0.6112 MeV, 10.1 fm). The HypHI Collaboration, in their announcement of a possible $\Lambda n n$ bound state, did not state any value for its binding energy. Neither was a spin and parity assigned. There is therefore no experimental results to compare with. Notwithstanding, one can immediately suspect that the second value is too high because there are no $\Lambda n$ or $n n$ bound states, compared to $\Lambda n p$ for which there is a $n p$ bound state.

The most likely binding energy from our calculation (-0.0695 MeV) is close to that reported in \cite{afn2015, gib2019} as -0.069 MeV. The result from \cite{afn2015, gib2019} is obtained by scaling up meson theory $\Lambda N$ potentials to about 35\% in a few-body calculation carried out through the Complex Energy Method. In \cite{bel2008}, a bound state was found by scaling up their $\Lambda N$ potential by 50\%, though the value of the binding energy computed was not stated. Motivated by the charge asymmetry relation between $\Lambda p$ and $\Lambda n$ forces, an important question arises from the studies in \cite{bel2008, afn2015, gib2019}, and all other studies that computed a $\Lambda n n$ bound state by scaling up the attraction in the $\Lambda n$ force: what is the binding energy of $\Lambda n p$ after both $\Lambda p$ and $\Lambda n$ are scaled up by the same amount? In \cite{hiy2014}, this question was explored using a lambda-nucleon force which is phase-equivalent with the NSC1997f force: by scaling up the tensor component in this force (which acts in the $^3\text{S}_1$ state) by 1.20, a bound state is found for $\Lambda n n$ at -0.054 MeV, and the corresponding lambda separation energy for $\Lambda n p$ computed is $B_{\Lambda}=0.83$ MeV. This value of $B_{\Lambda}$ is nearly the same as that presented in this paper for GLM-YN0 potentials ($B_{\Lambda}=0.85$ MeV).

Recently, a new experiment to search for a bound state or resonance in $\Lambda n n$ ($J^{\pi}=1/2^+$) was carried out at JLab \cite{tan2018}. At the time of writing this paper, data from this experiment was still being analysed.

\section{Conclusions}

In this paper, charge asymmetric lambda-nucleon potentials have been proposed. These potentials were obtained by averaging spin-dependent lambda-nucleon potentials recovered through the Gel'fand-Levitan-Marchenko method on theoretical scattering data. The potentials (GLM-YN0) are constructed from one-quarter spin singlet plus three-quarters spin triplet potentials. GLM-YN0 potentials were used in solving the $\Lambda NN$ three-body problem, together with Malfliet-Tjon I/III nucleon-nucleon potential. $\Lambda n p$ ($J^{\pi}=1/2^+$) binding energy and root-mean-square matter radius computed were found to be (-3.0759 MeV, 7.7 fm). These two results reveal that $\Lambda p$ and $\Lambda n$ forces are much stronger, a fact that is strongly evident from the 2018 $\Lambda n p$ ($J^{\pi}=1/2^+$) lifetime measurements by the STAR Collaboration. GLM-YN0 $\Lambda n$ forces gave rise to a bound $\Lambda n n$ ($J^{\pi}=1/2^+$) state, whose binding energy and root-mean-square radius are reported as (-0.0695 MeV, 12.3 fm) or (-0.6112 MeV, 10.1 fm). This $\Lambda n n$ ($J^{\pi}=1/2^+$) bound state arises \textit{naturally} within GLM-YN0 potentials. Computing $\Lambda n p$ ($J^{\pi}=1/2^+$) lifetime and the size of charge symmetry breaking in isospin multiplets are decisive tests on the accuracy of our lambda-nucleon forces from inverse scattering theory.

\section{References}

    \bibliographystyle{elsarticle-num-names}