\section{Supplemental Material for \\ ``Athermal Fracture of Elastic Networks:\\ How Rigidity Challenges the Unavoidable Size-Induced Brittleness''}

\renewcommand{\thefigure}{S\arabic{figure}}
\noindent

\section{Simulation methods}
We generate networks by diluting regular lattices consisting of $L \times L$ nodes with spacing $l_0$ and different topologies. A fraction $p$ of all possible nearest-neighbor bonds is present. The associated average network connectivity is therefore $\langle z \rangle = p \, z_{\text{max}}$, with $z_{\text{max}}$ depending on network topology (see also Fig.~\ref{SfigB}). All bonds are harmonic springs with stiffness $\mu$ and rest length $l_0$ and break irreversibly when exceeding the threshold $\lambda$ (same for all springs). By setting $l_0=1$ and $\mu=1$, we can express all quantities in reduced units, and $\lambda$ can be equivalently considered as deformation or force threshold. The energy of the system reads $\mathcal{H}=\frac{\mu}{2l_0}\sum_{\langle ij \rangle} {(l_{ij}-l_0)}^2$, where $\langle ij \rangle$ indicates bonded pairs of nodes and $l_{ij}$ the distance between them. We deform the networks uniaxially under athermal and quasistatic loading conditions, using the FIRE algorithm for the energy minimization~\cite{bitzek2006}. The tolerance $F_{\text{RMS}}$ is typically set at $10^{-5}$ for $L \leq 256$, but it is varied for larger sizes as discussed in the text. In our case, $F_{\text{RMS}}$ is the maximum numerical tolerance for both the system root-mean-square force~\cite{bitzek2006} and the maximum force on each single node. The networks are deformed uniaxially in the $y$-direction by employing a very small strain increment $\Delta \epsilon=0.001$, after which the energy minimization procedure is performed to obtain the equilibrium (off-lattice) node positions. Based on the mechanically equilibrated configuration, we identify the bonds that exceed $\lambda$ (if any) and remove the weakest bond defined as the one that exceeds the threshold the most, i.e. the one with the largest $(l_{ij}-l_0)/\lambda$. In most of the simulations, we additionally allow the strain increment to be reduced by up to two order of magnitude ($\Delta \epsilon=0.00001$) if the most stressed bond is very close to $\lambda$, to ensure that we are in the quasistatic condition and only one bond at a time can be broken. Typically, we use periodic boundary conditions in all directions. However, following Ref.~\cite{zhang2017}, we also perform simulations in which the nodes on the top and bottom boundaries are free to move only in the $x$-direction and move affinely in the $y$-direction, without noticing any appreciable qualitative difference. Furthermore, as detailed below, we also perform simulations by employing different boundary conditions in the $x$-direction. We confirmed that the conclusions of our study are robust with respect to all the variations on the simulation protocol.

\section{Observables and data analysis}
We define the mechanical response of the network by considering the stress along the direction of the deformation. In particular, we use the $yy$-component of the virial stress tensor $\sigma \equiv \sigma_{yy}^{\text{vir}} = \frac{1}{d V} \sum_{\langle ij \rangle} f_{ij,y} \, r_{ij,y}$ , where the sum runs over all the bonded pairs of nodes $\langle ij \rangle$, $\mathbf{f}_{ij}$ is the force acting on node $i$ due to $j$, $\mathbf{r}_{ij}$ is the vector connecting the two nodes, $d$ is the dimensionality of the network and $V$ is its (instantaneous) volume. Using a different definition for the macroscopic stress (e.g. $\sigma_{xx}^{\text{vir}}+\sigma_{yy}^{\text{vir}}$) does not change the conclusions of the paper. A network is defined rigid if $\sigma$ is greater than a small value ($10^{-4}$) and $\epsilon_{\text{R}}$ is the first strain for which this condition is satisfied. $\sigma_{\text{c}}$ is the maximum in the stress response and $\epsilon_{\text{p}}$ the associated strain. The work of extension $W_{\text{f}}$ is calculated as the integral of the stress-strain curve until fracture (network broken in two parts). $\Delta \sigma_{\text{max}}$ is the maximum stress drop and it is calculated following a procedure highlighted in Ref.~\cite{bonfanti2018}. It consists of (i) calculating the derivative of the stress-strain curve; (ii) identifying the strain interval for which the derivative is greater than a certain threshold (here we use zero, but we checked that our results are robust using different thresholds); (iii) for each strain interval, calculate the stress drop as the difference in the stress associated to the extremal values of the interval; (iv) the maximum stress drop is $\Delta \sigma_{\text{max}}$. In most of the cases, this procedure gives exactly the same result as just considering the difference in stress between two consecutive strains; whereas in few cases the values are slightly different but the overall trends are consistent with both definitions.

All quantities are averaged over many configurations to ensure proper statistical sampling of the disorder. Errors are calculated as standard deviation divided by the number of sampled configurations (standard error on the mean) and are shown when they are larger than the symbols displayed in the graphs. For example, for the size scaling analysis for networks close to the isostatic point ($p=0.65$ and $p=0.70$) and for small $\lambda$ we employ more than 2000 configurations for small system sizes ($L \leq 32$); 1000 for $L=50$; 500 for $L=64$; 200 for $L=128,256$; 50 for $L=512$; and $\sim 10$ configurations for larger system sizes ($L=640,768,1024$) that are statistically sufficient due to the self-averaging nature of the fracture process in such large systems. In general, for all networks we average quantities over at least 200 configurations when $L=128$ and a larger number (typically 500) for $L=64$ or smaller sizes. Results for different network geometries are also averaged over approximately 200 configurations for 2D networks and over 20 to 100 configurations for 3D networks, depending on system size.

\section{Additional results}

\subsection{Fixed system size}
In Fig.~\ref{SfigA}, some additional results for triangular networks with fixed system size $L=128$ and different $p$ and $\lambda$ are shown. In particular, we note in Fig.~\ref{SfigA}(a) that for fixed $\lambda$ the work of extension $W_f$, i.e. the energy needed to break the material (calculated as the integral of the stress-strain curve), shows a maximum around the isostatic point ($p_{iso}\simeq0.66$) for small $\lambda$, analogously to the trend for $n_f$ (shown in Fig.~1(c) of the main text). This has been observed also in Ref.~\cite{zhang2017} and implies that networks close to the isostatic point can be counterintuitively considered tougher than networks with more bonds (when $\lambda \leq 0.10$). However, as shown in Fig.~\ref{SfigA}(b), when considering the maximum stress $\sigma_{\text{c}}$ as a function of connectivity we observe that there is no maximum and the more intuitive monotonic dependence on $p$ --- the more the bonds, the stronger the material --- is observed. This discrepancy between the dependence on connectivity of two quantities describing resistance to fracture is an effect of the size-dependent brittleness demonstrated in this study: at moderate system size, the networks are ductile and significant energy is dissipated in bond-breaking events occurring after the peak stress. To further quantify the tunability of the fracture process by varying $p$ and $\lambda$ at fixed system size, we show the peak strain $\epsilon_p$ in panel (c) and the maximum stress drop in panel (d).

\subsection{Different topologies}
In Fig.~\ref{SfigB}, we show $n_{\text{f}}$ as a function of the network average connectivity $\langle z \rangle=pz_{\text{max}}$ for different network topology and dimensionality, and small $\lambda$. We observe that the position of the maximum in $n_{\text{f}}$ depends not only on dimensionality but also on topology and orientation. For example, despite the fact that the diluted square (SQ) and the 2D diamond (D2) networks have the same connectivity properties, SQ is more brittle (fewer bonds are broken) because it features more bonds already aligned with the direction of deformation. Analogously, because of geometric effects (a large number of system-spanning bonds at small connectivity) in the honeycomb (HC) and in the 3D diamond (D3) networks, a larger number of bonds is broken when approaching the limit of undiluted lattices. Only for the diluted triangular (TR) and face-centred-cubic (FCC) networks, that are sufficiently isotropic, the maximum of $n_{\text{f}}$ roughly corresponds to the isostatic point $z_{\text{iso}}$. Our results enrich the data reported in Refs.~\cite{zhang2017,driscoll2016,berthier2019}. We can summarize that damage is maximized close to rigidity points, either of isostatic or geometric nature, when $\lambda$ is small.

\subsection{Size scaling}
In Fig.~\ref{SfigC}, additional results from the size scaling are shown. To extract the damage fractal dimension $d_{\text{f}}$ and generate the plot of Fig.~2(b) of the main text, we simulate at least six different system sizes (typically seven, up to $L=256$) and we fit $N_{\text{f}} \sim N_{\text{in}}^{d_{\text{f}}/2}$ where $N_{\text{in}}$ is the number of intact initial springs. This is analogous to $N_{\text{f}} \sim L^{d_{\text{f}}}$ but sometimes better fits are obtained with the previous expression since the values on the x-axis span more orders of magnitude. As shown in the main text, we notice that the fractal dimension $d_{\text{f}}$ is (much) greater than 1 for a large range of $p$ and $\lambda$. This is a distinct feature of our elastic network model compared to the scalar random fuse models~\cite{alava2006} or models assuming linear elasticity~\cite{nukala2005}, where the number of broken bonds at peak load scales \emph{at most} linearly with the system size. In our case $d_f$ refers to the exponent of the total broken bonds at failure and not to the exponent of broken bonds at peak load. However, since the fracture becomes more brittle (i.e. bonds break only before and at peak load) the exponent for the broken bonds at peak load is actually larger than $d_f$. We also show how the maximum stress $\sigma_{\text{c}}$ and peak strain $\epsilon_p$ decrease with increasing system size. Inspired by previous works~\cite{duxbury1987,manzato2012,roy2017b}, we fitted the decay of $\sigma_{\text{c}}$ with three functional forms: (i) logarithmic decay $\propto 1/\ln(L)$; (ii) Duxbury-Leath-Bale (DLB)~\cite{duxbury1987,manzato2012} $\propto 1/\sqrt{\ln(L^2)}$; (iii) power-law decay $\propto 1/L^{\beta}$. We found that the power-law decay fits our data the best, with an exponent $\beta$ that clearly depends on $p$ and $\lambda$. For example, we obtain $\beta=0.73$ for $p=0.65$, $\lambda=0.03$ (Fig.~2(d) of main text) and $\beta=0.89$ for $p=0.70$, $\lambda=0.03$ (Fig.~\ref{SfigC}(d)). Differently from previous models where it is assumed that the material strength (i.e. $\sigma_c$) is zero in the thermodynamic limit~\cite{alava2006}, in our case we have an offset value corresponding to a finite network strength $\sigma_c^{(\infty)} \neq 0$ that depends on $p$ and $\lambda$. Our data suggest that the spring network model studied here is qualitatively different from scalar models that have been so far considered as the golden standard for statistical physics studies of fracture.

\subsection{Universality of fracture abruptness trend}
In Figs.~\ref{SfigD} and~\ref{SfigE}, we show that the non-monotonic size-dependence of fracture abruptness $\Delta \sigma_{\text{max}}$ is universal with respect to network topology and robust with respect to the simulation procedure. In Fig.~\ref{SfigD}(a), we compare the trend for the same set of initial configurations simulated with a protocol where the small strain increment is kept fixed at $\Delta \epsilon=0.001$ and one where it can be adjusted to much smaller values (up to $\Delta \epsilon=0.00001$) if the most stressed bond is very close to $\lambda$. A minimum is observed in both cases for similar values of $L$. A slightly smaller $\Delta \sigma_{\text{max}}$ is observed for large $L$ when $\Delta \epsilon$ can be reduced by two orders of magnitude. Next, we check if the choice of boundary conditions in the $x$-direction, transverse to the deformation direction, has some qualitative effect on the fracture abruptness. We tested three different boundary conditions: (i) periodic (mimicking an infinitely large system); (ii) free, i.e. without bonds crossing the $x$-direction (the system is finite in that direction); (iii) with periodic bonds and with adjustable size in the $x$-direction such that the area is kept constant during the uniaxial deformation. As shown in Fig.~\ref{SfigD}(b), the qualitative trends are robust and independent of the fundamentally different boundary conditions employed (the first one was chosen for the rest of our study since it seems the most computationally efficient). Furthermore, we confirm that the non-monotonic trend persists when changing network topology. As already highlighted in the main text, the minimum and the maximum of $\Delta \sigma_{\text{max}}$ as a function of $L$ depend on network properties. This can be appreciated in Fig.~\ref{SfigD}(c) where we report selected cases for various 2D and 3D networks. Finally, preliminary results confirm the non-monotonic trend also for shear deformation, supporting the conclusion that it is a universal size-dependent behavior of diluted networks. As shown in Fig.~\ref{SfigE}, we also succeeded in manually collapsing the curves of $\Delta \sigma_{\text{max}}$ obtained for networks close to (just below and slightly above) the isostatic point and various $\lambda$ into a single master curve. The rescaling $L \rightarrow \alpha L$ and $\Delta \sigma_{\text{max}} \rightarrow \beta \Delta \sigma_{\text{max}}$ was performed by imposing that the different curves approximately share the same x-value for the minimum. We obtained values for $\alpha$ that are inversely proportional to $\lambda$, suggesting that the ``critical'' system size $L^* \sim 1/\alpha$ at which brittleness kicks in decreases by making the springs stronger. On the other hand, $\beta$ increases more than linearly with $\lambda$, consistently with the idea that due to inhomogeneous stress distribution, the increase of the springs strength has non-linear consequences on the overall network.

\subsection{Effect of tolerance in energy minimization}
Finally, in Fig.~\ref{SfigF} we report the size scaling of the ductility interval $\Delta \epsilon_{\text{duct}}$, defined as the strain difference between $\epsilon_p$ and the strain corresponding to the maximum stress drop, for two different values of the tolerance $F_{\text{RMS}}$ and various networks. Brittle fracture (vanishing $\Delta \epsilon_{\text{duct}}$) is observed in the athermal limit, whereas a non-monotonic size-dependence of the ductility interval is observed for larger $F_{\text{RMS}}$. The system size at which ductility is increasing again (the minimum in $\Delta \epsilon_{\text{duct}}$ as a function of $L$) does strongly depend on $\lambda$ and therefore seems to be coupled with the network rigidity. If a variation of the tolerance in the energy minimization protocol $F_{RMS}$ is interpreted as a variation of the thermal effects in the system (that would keep the system not exactly at the energy minimum), these results might suggest that by increasing temperature it is possible to observe a more ductile fracture occurring via many large cracks. Future research should further investigate the intriguing possibility that the minimum temperature at which such ductile behavior is first observed, and the material rigidity are strongly coupled. However, to precisely assess the role of temperature, more traditional simulation methods for thermal systems should be employed.

\clearpage

\begin{figure*}[h!t]
	\begin{center}
		\includegraphics[width=0.9\textwidth]{./SfigA.pdf}
	\end{center}
	\caption{Additional results for triangular networks with $L=128$ as a function of $p$ and various $\lambda$ (see legend in panel (b)). (a) Work of extension $W_{\text{f}}$, defined as the integral of the stress-strain curve. (b) Maximum stress $\sigma_{\text{c}}$. (c) Strain $\epsilon_{\text{p}}$ corresponding to the stress peak as a function of $p$ for various $\lambda$. Inset: strain $\epsilon_{\text{R}}$ at which the network is rigid upon uniaxial extension as a function of $p$. (d) Maximum stress drop $\Delta \sigma_{\text{max}}$ versus $p$ for different $\lambda$.  }
	\label{SfigA}
\end{figure*}

\begin{figure*}[h!t]
	\begin{center}
		\includegraphics[width=0.75\textwidth]{./SfigB.pdf}
	\end{center}
	\caption{Fraction of broken bonds as a function of averaged network connectivity $\langle z \rangle$ for (a) 2D and (b) 3D networks with the different architectures shown on the right, same $\lambda$ and similar system size. Maxima are observed close to rigidity points that depend on topology, geometry and dimensionality.}
	\label{SfigB}
\end{figure*}

\begin{figure*}[h!t]
	\begin{center}
		\includegraphics[width=0.9\textwidth]{./SfigC.pdf}
	\end{center}
	\caption{(a) Number of broken bonds at failure $N_{\text{f}}$ as a function of number of initial intact springs $N_{\text{in}}$ for various combination of $p$ and $\lambda$. Dashed lines are best power-law fits from which the damage fractal dimension $d_{\text{f}}$ is extracted. (b) Size scaling of broken bonds for triangular network close from above to the isostatic point with small $\lambda$. Black circles indicate single-run values, red squares averages and dashed line best power-law fit. (c) Size scaling of strain at peak stress $\epsilon_{\text{p}}$ for networks in the proximity of the isostatic point with small $\lambda=0.03$. Symbols represent averages and lines are best fits. Data are well fitted by a power-law decay with an offset $\epsilon_{\text{p}} = \epsilon_{\text{p}}^{(\infty)} + A/L^{\beta}$, with $\epsilon_{\text{p}}^{(\infty)}$ the strain at peak in the thermodynamic limit. For $p=0.65$, we obtain $\epsilon_{\text{p}}^{(\infty)}=0.050$ and $\beta=0.47$. For $p=0.70$, we obtain $\epsilon_{\text{p}}^{(\infty)}=0.016$ and $\beta=0.24$. (d) Size scaling of the network strength $\sigma_{\text{c}}$ (log-log plot). Symbols correspond to averages, error bars show the standard deviations and lines are the best fits. Three functional forms were used: logarithmic decay (dotted black line), DLB decay (dashed red line), and power-law decay $\sigma_{\text{c}} = \sigma_{\text{c}}^{(\infty)} + A/L^{\beta}$ (green solid line). The latter fits the data the best, we obtain $\sigma_{\text{c}}^{(\infty)}=0.0053$ and $\beta=0.89$. }
	\label{SfigC}
\end{figure*}

\begin{figure*}[h!t]
	\begin{center}
		\includegraphics[width=0.8\textwidth]{./SfigD.pdf}
	\end{center}
	\caption{(a) Size scaling of maximum stress drop for two different simulation protocols. Black circles represent data from simulations where the strain step can be adjusted (reduced) to smaller steps when there are bonds very close to their thresholds (see section on simulation methods); whereas blue squares are data from simulations where the strain step is kept fixed at a small value. (b) Size scaling of maximum stress drop for four networks with $p$ and $\lambda$ as indicated in the legend, and three different simulation protocols for which periodic boundary conditions in the $x$ direction (normal to the deformation direction) are present (solid lines) or not (dashed lines), or for which the $x$ direction is adjusted to keep the system area constant (dotted lines). (c) Size scaling of fracture abruptness for diluted 2D honeycomb (HC), 2D square (SQ) and 3D FCC (inset) networks with averaged connectivity $\langle z \rangle$ and $\lambda$ as indicated.}
	\label{SfigD}
\end{figure*}

\begin{figure*}[h!t]
	\begin{center}
		\includegraphics[width=0.9\textwidth]{./SfigE.pdf}
	\end{center}
	\caption{Rescaled maximum stress drop $\Delta \sigma_{\text{max}}$ as a function of rescaled system size $L$ for (a) $p=0.65$ (see also Fig.~3(d) of main text) and (b) $p=0.70$, color-coded according to $\lambda$. Symbols represent averages and error bars are standard deviations. (c) Both data sets combined (circles for $p=0.65$, squares for $p=0.70$). The values used for the rescaling parameters $\alpha$ and $\beta$ are shown in the insets as a function of $\lambda$. $\beta$ grows (more than linearly) with increasing $\lambda$. $\alpha$ is related to the transition system size $L^*$ discussed in the main text and decreases by two order of magnitudes in the simulated range of $\lambda$. }
	\label{SfigE}
\end{figure*}

\begin{figure*}[h!t]
	\begin{center}
		\includegraphics[width=0.9\textwidth]{./SfigF.pdf}
	\end{center}
	\caption{Size scaling of the ductility interval $\Delta \epsilon_{\text{duct}} = \epsilon_{\Delta\sigma_{\text{max}}} - \epsilon_{\text{p}}$, where $\epsilon_{\Delta\sigma_{\text{max}}}$ is the strain at which the maximum stress drop is observed (see also Fig.~4 of the main text), for simulations with (a) $F_{\text{RMS}}=10^{-5}$ and (b) $F_{\text{RMS}}=10^{-6}$ for different $p$, $\lambda$ as indicated in the legend. For small $F_{\text{RMS}}$ (panel (b)), it is evident that fracture becomes more brittle upon increasing system size. For larger $F_{\text{RMS}}$ (panel (a)), that are nevertheless three orders of magnitude smaller than $\lambda$, a non-monotonic trend is observed and networks with large $L$ reacquire ductility. At a first glance, the minima in $\Delta \epsilon_{\text{duct}}$ seem to roughly correspond with the minima in $\Delta \sigma_{\text{max}}$ (see Fig.~3(a) of main text), but further studies are needed to confirm this correspondence. }
	\label{SfigF}
\end{figure*}

\end{document}


\section{Introduction}
    Lung cancer has the highest incidence and mortality rates worldwide \cite{bray2018global}. Early diagnosis and treatment of pulmonary nodules can increase the survival rate of patients. Computed tomography (CT) has been widely used and proved effective for detecting pulmonary nodules. However, manually identifying nodules in CT scans is often time-consuming and tedious, because a radiologist needs to read the CT scans slice by slice, and a chest CT may contain over 200 slices. Accurate and precise nodule segmentation can provide more in-depth assessment of the shape, size and change rate of the nodule. When nodule is identified, a follow up scan in 3 - 12 months is usually required to assess its growth rate \cite{kalpathy2016comparison}. The growth of the lung tumor may be an indicator for malignancy, and an accurate nodule segmentation can be used for measuring the growth rate of the nodule.

In recent years, deep convolutional neural network has emerged as a leading method for automatically detecting and segmenting pulmonary nodules and have achieved great success. State-of-the-art frameworks for nodule detection often ustilize the 3D region proposal network (RPN) \cite{ren2015faster} for nodule screening \cite{tang2018automated,zhu2017deeplung,setio2017validation,liao2019evaluate}, followed by a 3D classifier for false positive reduction \cite{ding2017accurate,tang2019end}. Although single stage detector has also been proposed in \cite{khosravan2018s4nd}, their hit criteria was different from what was more commonly adopted \cite{setio2017validation}. Moreover,
the refinement provided by the extra classifiers may correct some errors made by the detectors. In terms of nodule segmentation, U-Net \cite{ronneberger2015u} and V-Net \cite{milletari2016v} like structure is predominantly used \cite{wang2017central,wu2018joint,aresta2018iw}.
In practice, a computer aided diagnosis (CAD) system for pulmonary nodule detection and segmentation often consists of several independent subsystems, optimized separately.

There are some limitations on handling each task completely independent. First, it is time-consuming and resource intensive to train several deep convolutional neural networks. Although each component is designed for different purposes, they share the common procedure of extracting feature representations that characterize pulmonary nodules. Second, the performance of the whole system may not be optimal, because separately training several systems prevents communication between each other and learning intrinsic feature representations. Intuitively, the segmentation mask of the nodule should provide a strong guide for the neural network to learn discriminative features, which may in turn improve the performance of nodule detection.

Although multi-task learning (MTL) and feature sharing offer an attractive solution to combine different tasks, a naive implementation may cause other problems \cite{cheng2018revisiting}. First, because of the mismatched goals of localization and classification, it may be sub-optimal if these two tasks are performed using the same feautre map. Second, a large receptive field may integrate irrelevant information from other parts of the image, which may negatively affect and confuse the classification of nodules, especially small ones. \cite{cheng2018revisiting} decoupled localization and classification to address the problem in natural imaging. However, completely separating the two tasks without sharing any feature extraction backbone, still prevents cross-talk between two networks and may not be the most efficient. Therefore, a decoupled false positive reduction, that pools features from early scales of the feature extraction backbone, is proposed to address this problem, which allows learning both task-independent and task-dependant features.

Here, we propose a new end-to-end framework, called NoduleNet, for solving pulmonary nodule candidate screening, false positive reduction and segmentation jointly. NoduleNet consists of three parts: nodule candidate screening, false positive reduction and segmentation refinement (\Cref{fig:model}). These three components share the same underlying feature extraction backbone and the whole network is trained in an end-to-end manner.

\begin{figure}[!h]
\centering
 \includegraphics[width=0.9\textwidth]{model.pdf}
\caption{\bf Overview of NoduleNet. NoduleNet is an end-to-end framework for pulmonary nodule detection and segmentation, consisting of three sequential stages: nodule candidate screening, false positive reduction and segmentation refinement. $k$ is the number of anchors. FC is short for fully connected layer.}
\label{fig:model}
\end{figure}

Our main contributions are summarized as follows:
\begin{itemize}
    \item We propose a unified model to integrate nodule detection, false positive reduction and nodule segmentation within a single framework, trained end-to-end in a multi-task fashion.
    \item We demonstrate the effectiveness of the model, improving nodule detection accuracy by 10.27\% compared to the baseline model trained only for nodule detection, and achieving a state-of-the-art nodule segmentation accuracy of 83.10\% on Dice-Sørensen coefficient (DSC).
    \item We carry out systematic ablation studies to verify the contributions of several design tricks underlying NoduleNet, including decoupled features maps, segmentation refinement subnet, and multi-task training.
\end{itemize}

\section{NoduleNet}

\subsubsection{Nodule candidate screening (NCS)}
To generate nodule candidates, a $3 \times 3 \times 3$ 3D convolutional layer is applied to the feature map (feature\_map\_4 in Figure~\ref{fig:model}), followed by two parallel $1 \times 1 \times 1$ convolutional layers to generate classification probability and six regression terms associated with each anchor at each voxel on the feature map.
An anchor is a 3D box, which requires six parameters to specify: central z-, y-, x- coordinates, depth, height and width. We chose cube of size 5, 10, 20, 30 and 50 as the 5 anchors in this work. Then, we minimize the same multi-task loss function as \cite{ren2015faster}.

\subsubsection{Decoupled false positive reduction (DFPR)}
Unlike \cite{ren2015faster} that performs classification using features pooled from the same feature map as RPN (feature\_map\_4). Learning using coupled feature map may lead to sub-optimal solutions of the two tasks. Instead, we use 3D region of interest (ROI) pooling layer to pool features from early feature map that has a small receptive field (down\_4). This not only ensures the false positive reduction network has a small receptive field and can learn feature representations that are substantially different from nodule candidate screening network, but also allows sharing of a few feature extraction blocks. The false positive reduction network minimizes the same multi-task in loss function as the NCS
.

\subsubsection{Segmentation refinement (SR)}
As shown in \Cref{fig:model}, segmentation is performed at the same scale of the original input CT image, by progressively upsampling the cropped high-level feature map (feature\_map\_4) and concatenating them with low-level semantically strong features.

This approach is fundamentally different from the mask branch proposed in \cite{he2017mask}. In \cite{he2017mask}, the authors perform segmentation by only using downsampled feature map and then resize the predicted mask back to the original image scale, which may lose precision due to bounding box regression errors and loss of more fine-grained local features.

Another advantage is that, only the regions have nodules are upsampled to the original image scale, which only accounts for a small area of the whole input image. This saves a large amount of GPU memory, making whole volume input feasible during training and testing, as compared to upsampling the whole feature map to original input scale in \cite{milletari2016v}.

The segmentation refinement network minimizes the soft dice loss of the predicted mask sets $\{m\}$ and the ground truth mask sets $\{g\}$ of the input image.

\section{Results}
\subsubsection{Data and experiment configurations}
We used LIDC-LDRI \cite{armato2011lung} for evaluting the performance of NoduleNet. LIDC-LDRI is a large-scale public dataset for studying lung cancers, which contains 1018 sets of CT scans collected from multiple sites with various slice thickness. Nodules with diameter equal or greater than 3 mm in this dataset have contour outlined by up to four radiologists. We included only those CT scans met the selection criteria of LUNA16 \cite{setio2017validation} in this work. If the two segmentation masks provided by two radiologists have an intersection over union (IoU) larger than 0.4, we consider the two masks are referring to the same nodule. We consider nodules annotated by at least 3 out of 4 radiologists the ground truth, resulting in a total number of 586 CT scans with 1131 nodules.
Note that the number of CT scans and nodules included in this work may be different from previous work \cite{wang2017central,wu2018joint,aresta2018iw}, due to different inclusion criteria.

A six-fold cross validation was performed to demonstrate the performance of NoduleNet. All models in the experiment were trained using stochastic gradient descent (SGD) with initial learning rate 0.01, momentum 0.9 and $l2$ penalty 0.0001, for 200 epochs. The learning rate was scheduled to decrease to 0.001 after 100 epochs and to 0.0001 after another 60 epochs.

Free-Response Receiver Operating Characteristic (FROC) \cite{kundel2008receiver} analysis was adopted for evaluating the performance of nodule detection. We used the same hit criterium and competition performance metric (CPM) as in the LUNA16 \cite{setio2017validation}. Intersection over union (IoU) and S${\o}$rensen-Dice coefficient (DSC) were used for evaluating the performance of nodule segmentation.

\subsubsection{Nodule detection performance}
In order to fully verify and understand our aforementioned assumptions, we conducted extensive experiments using different network architectures and design choices. We use \textbf{$\textrm{N}_1$} to represent network that only has NCS branch, \textbf{$\textrm{N}_2$} for network has both NCS and FPR branches, and \textbf{$\textrm{N}_3$} for network has all NCS, FPR and SR branches. \textbf{$\textrm{F}_c$} represents the FPR branch is built on the same feature map as NCS, and \textbf{$\textrm{F}_d$} means the FPR branch is built on the decoupled feature map mentioned in previous section. \textbf{R} means the training data is extraly augmented with $xy$ - plane rotation. \textbf{NCS} means the predicted probability comes from NCS branch, \textbf{FPR} means the predicted probability comes from FPR branch, and \textbf{FU} means the predicted probability is fused from NCS and FPR. Note that \textbf{$\textrm{N}_1$} is the widely used 3D RPN for nodule detection \cite{tang2018automated,tang2019end,zhu2017deeplung,liao2019evaluate}, which was served as a strong baseline for evaluating the performance of each added component. The results are summarized in \Cref{table:performance}.

As seen from \Cref{table:performance}, the sensitivity at 8 false positives per patient rate has a consistent improvement of 1.0\% to 1.5\% by adding the segmentation refinement network (\textbf{$\textrm{N}_3$}), which demonstrates the effectiveness of using the extra nodule segmentation information.

The average sensitivity of the NoduleNet using decoupled false positive reduction (\textbf{$\textrm{F}_d$}) has around 3\% to 4\% improvement over the NoduleNet using coupled false positive (\textbf{$\textrm{F}_c$}). Moreover, by adding rotation in data augmentation (\textbf{R}), the performance of \textbf{FPR} branch is further improved by around 2.5\% while the performance of \textbf{NCS} branch remains almost the same. This verifies our assumption that classification should learn invariant features, while localization may learn co-variant features. Those findings demonstrate the importance of decoupling modules that are essentially learning different tasks.

By fusing the predicted probability from NCS and FPR, the performance was consistently improved by 0.7\% - 1.0\%, demonstrating that combining predictions from branches that perceive different level of context information is important.

By adding false positive reduction and segmentation refinement network, the performance of the baseline detector (\textbf{NCS}) is correspondingly improved, showing the effectiveness of multi-task learning and feature sharing.

All together, NoduleNet outperforms a strong baseline single stage detector by 10.27\%. Note that performance reported in LUNA16 may not be directly comparable to this work, because of different nodule selection criteria, and training and testing data splits. Also, this work focuses on the joint learning of nodule detection and segmentation, whereas the LUNA16 focuses only on nodule detection.

\begin{table}[!ht]
\centering
\begin{tabular}{l c c c c c c c c}
\hline
\textbf{Method} & 0.125 & 0.25 & 0.5 & 1.0 & 2.0 & 4.0 & 8.0 & Avg.\\ \hline
\textbf{\textbf{$\textrm{N}_1$} (\textbf{NCS}) \cite{tang2018automated,tang2019end,zhu2017deeplung,liao2019evaluate}} & 52.17 & 62.51 & 71.09 & 80.46 & 87.27 & 91.07 & 94.43 & 77.00 \\
\textbf{\textbf{$\textrm{N}_2$} + \textbf{$\textrm{F}_c$} (\textbf{NCS})} & 53.85 & 62.07 & 71.09 & 79.22 & 86.74 & 90.98 & 93.28 & 76.75 \\
\textbf{\textbf{$\textrm{N}_2$} + \textbf{$\textrm{F}_c$} (\textbf{FPR}) \cite{tang2019end}} & 55.79 & 66.93 & 75.77 & 82.40 & 88.68 & 91.78 & 93.10 & 79.21 \\
\textbf{\textbf{$\textrm{N}_3$} + \textbf{$\textrm{F}_c$} (\textbf{NCS})} & 53.67 & 63.84 & 74.62 & 83.20 & 88.51 & 92.04 & 94.96 & 78.69 \\
\textbf{\textbf{$\textrm{N}_3$} + \textbf{$\textrm{F}_c$} (\textbf{FPR})} & 57.38 & 65.96 & 77.19 & 84.97 & 89.92 & 93.28 & 95.40 & 80.59 \\ \hline
\textbf{\textbf{$\textrm{N}_2$} + \textbf{$\textrm{F}_d$} (\textbf{NCS})} & 56.15 & 66.93 & 74.54 & 82.23 & 88.59 & 92.22 & 95.05 & 79.39 \\
\textbf{\textbf{$\textrm{N}_2$} + \textbf{$\textrm{F}_d$} (\textbf{FPR})} & 61.98 & 71.26 & 78.78 & 85.41 & 89.30 & 92.22 & 95.31 & 82.04 \\
\textbf{\textbf{$\textrm{N}_3$} + \textbf{$\textrm{F}_d$} (\textbf{NCS})} & 61.45 & 70.20 & 78.16 & 84.62 & 90.27 & 93.63 & 96.20 & 82.08 \\
\textbf{\textbf{$\textrm{N}_3$} + \textbf{$\textrm{F}_d$} (\textbf{FPR})} & 68.08 & 73.56 & 81.70 & 85.94 & 90.80 & 93.90 & 96.55 & 84.36 \\
\textbf{\textbf{$\textrm{N}_3$} + \textbf{$\textrm{F}_d$} (\textbf{FU})} & 68.70 & 75.60 & 82.23 & 87.36 & 92.04 & 94.96 & 96.46 & 85.34 \\ \hline
\textbf{\textbf{$\textrm{N}_3$} + \textbf{$\textrm{F}_d$} + \textbf{R} (\textbf{NCS})} & 62.78 & 70.65 & 78.43 & 84.44 & 89.74 & 93.10 & 95.49 & 82.09 \\
\textbf{\textbf{$\textrm{N}_3$} + \textbf{$\textrm{F}_d$} + \textbf{R} (\textbf{FPR})} & 69.23 & 77.01 & 84.70 & 89.48 & 93.37 & 95.23 & $\textbf{96.55}$ & 86.51 \\
\textbf{\textbf{$\textrm{N}_3$} + \textbf{$\textrm{F}_d$} + \textbf{R} (\textbf{FU})} & $\textbf{70.82}$ & $\textbf{78.34}$ & $\textbf{85.68}$ & $\textbf{90.01}$ & $\textbf{94.25}$ & $\textbf{95.49}$ & 96.29 & $\textbf{87.27}$ \\ \hline
\end{tabular}
\caption{\bf CPM of different methods on the LIDC dataset based on six-fold cross validation. Shown are nodule detection sensitivities (unit: \%) with each column denoting the threshold false positive rate per CT scan (FPs/scan). The last column denotes the average sensitivities across the seven pre-defined FPs/scan thresholds. }
\label{table:performance}
\end{table}

\subsubsection{Nodule segmentation performance}
In \Cref{table:performance_IOU}, we compared the segmentation performance of NoduleNet to other deep learning based methods trained and tested on LIDC dataset \cite{wang2017central,wu2018joint,aresta2018iw}. NoduleNet outperformed previous state-of-the-art deep learning based method by 0.95\% on DSC, without the need to train a separate and dedicated 3D DCNN for nodule segmentation.
We randomly selected several nodules for visualizing the segmentation quality (\Cref{fig:visualization}).

\begin{table}[!ht]
\centering
\begin{tabular}{l c c c c c}
\hline
\textbf{Approach} &  \multicolumn{2}{c}{\textbf{\# Nodules}} & \textbf{\# Consensus} & \textbf{IoU} (\%) & \textbf{DSC} (\%)\\
 & \textbf{train} & \textbf{test} & \\ \hline
Wu \textit{et al}.  \cite{wu2018joint} & 1404 & 1404 & 3 & N\textbackslash A & 73.89 $\pm$ 3.87 \\
Aresta \textit{et al}.  \cite{aresta2018iw} & 1593 & 1593 & 3 & 55.00 $\pm$ 14.00 & N\textbackslash A \\
Wang \textit{et al}.  \cite{wang2017central} & 350 & 493 & 4 &71.16 $\pm$ 12.22 & 82.15 $\pm$ 10.76 \\ \hline
NoduleNet & 1131 & 1131 & 3 & 69.98 $\pm$ 10.80 & 81.80 $\pm$ 8.65 \\
NoduleNet & 1131 & 712 & 4 & \textbf{71.85 $\pm$ 10.48} & \textbf{83.10 $\pm$ 8.85} \\ \hline
\end{tabular}
\caption{\bf IoU (\%) and DSC (\%) performance of nodule segmentation between different methods. ``\# Consensus'' means each method includes nodules that are annotated by at least ``\# Consensus'' experts.}
\label{table:performance_IOU}
\end{table}

\begin{figure}[!h]
\centering
 \includegraphics[width=0.9\textwidth]{visualization2.pdf}
\caption{\bf Examples of nodule segmentation generated by NoduleNet.}
\label{fig:visualization}
\end{figure}

\vspace{-0.5cm}
\section{Conclusion}
In this work, we propose a new end-to-end 3D DCNN, named NoduleNet, for solving pulmonary nodule detection, false positive reduction and segmentation jointly. We performed systematic analysis to verify the assumptions and intuitions behind the design of each component in the architecture. Cross validation results on LIDC dataset demonstrate that our model achieves a final CPM score of 87.27\% on nodule detection and DSC score of 83.10\% on nodule segmentation, representing current state-of-the-arts on this dataset. The techniques introduced here are general, and can be readily transferred to other models.

\bibliographystyle{splncs04}


\section{Introduction}

Metallic nanoparticles are an important optoelectronic platform because of their
strong plasmon resonance at visible energies, which can be tuned by size, shape,
and environment\cite{Willets2007,Odom2011}; however, the accessible carrier
densities are limited to those of the parent metals and are typically
$10^{22}$~cm$^{-3}$ or higher.  Recently, the doping of semiconductor
nanoparticles has enabled access to
much lower electron densities~\cite{Mocatta2011,Luther2011,Scholl2012,Schimpf2014,Schimpf2015,Shen2016}.
These doped nanoparticles, with their tunable charge carrier density, exhibit
strong intraband absorption in a wide range of the THz
regime, enabling promising new infrared and plasmonic applications~\cite{Scotognella2013,Faucheaux2014,Kriegel2017,Lhuillier2017}.

Chemical techniques have enabled semiconductor nanoparticles to be doped with
as few as 1-100 electrons. In this regime, the classical electrostatic picture of
plasmons breaks down, demanding a theory of so-called ``quantum
plasmons''~\cite{Pitarke2007,Morton2011,Bernadotte2013,Tame2013,Krauter2014,Casanova2016,Varas2016,Bursi2016,Zhang2017}.
Furthermore, the collective plasmon picture
becomes dubious for systems containing only 1-10 excess electrons, suggesting
a transition to ``single-particle'' excitations~\cite{Shen2016} or -- as we will argue --
excitonic transitions.
Here, we aim to present a detailed quantum mechanical understanding of the
microscopic nature of these intraband excitations over a range of
experimentally-relevant sizes and densities.

The classical Mie theory of plasmons in nanoparticles predicts a localized surface plasmon resonance when $\varepsilon_{1}(\Omega_\mathrm{LSPR})=-2\varepsilon_{\mathrm{m}}$, where
$\varepsilon_1(\omega)$ is the real part of the complex bulk dielectric function for the nanoparticle, $\varepsilon(\omega) = \varepsilon_1(\omega) + i\varepsilon_2(\omega)$
and $\varepsilon_{\mathrm{m}}$ is the dielectric constant of the medium\cite{Willets2007}.
The Drude plasmon pole approximation to the dielectric function,
$\varepsilon(\omega)
    = \varepsilon_\infty - \Omega_\mathrm{p}^2/(\omega^2+i\gamma\omega)$,
leads to the expression for the plasmon frequency
$
\Omega_{\mathrm{LSPR}} = \sqrt{
    \Omega_\mathrm{p}^2/(\varepsilon_{\infty}+2\varepsilon_{\mathrm{m}})-\gamma^2}$,
where $\Omega_\mathrm{p}=\sqrt{4\pi \rho}$ is the bulk plasma frequency, $\rho$ is the
free charge carrier density, $\varepsilon_\infty$ is the high-frequency
dielectric constant, and $\gamma$ is the scattering rate of the electrons.
For a given material and medium, this Mie+Drude plasmon frequency only
depends on the density and not on the size, therefore failing to
account for quantum confinement effects in small nanoparticles. The Drude dielectric function can be replaced with
a microscopic dielectric function that accounts for quantum confinement effects, but this approach typically neglects
interparticle interactions\cite{Kraus2017,Genzel1975,Scholl2012,Jain2014,Schimpf2014}.

The prevailing quantum mechanical theory of plasmons in metals is the
random-phase approximation (RPA)~\cite{Bohm1953,GellMann1957}. As a theory of the ground-state energy density
of bulk metals, the RPA famously removes the divergences encountered in
finite-order perturbation theory. As a theory of the dynamical response, the RPA
predicts the collective plasmon excitation, including its dispersion and strong
oscillator strength~\cite{PinesNozieres,GiulianiVignale}. Despite its success for simple bulk metals, the RPA (by which
we mean the time-dependent Hartree approximation) is not
an accurate theory of excitation energies in molecules, casting doubt
on its applicability to quantum plasmonics. In particular, the RPA fails to
describe bound states such as excitons. To go beyond the RPA requires the tools
of higher-level many-body theory or quantum chemistry.

In this paper, we investigate a model of interacting electrons confined to the interior of a sphere of radius $R$.  This model generalizes the uniform electron gas
(UEG), sometimes referred to as ``jellium'', which is the canonical model of
bulk metals and their plasmonic excitations. In the $R\rightarrow \infty$ limit,
our model approaches the UEG (up to a background charge density).  At finite
$R$, quantum confinement produces a one-particle spectrum that is gapped, which
alters the nature of the dominant excitations. We develop the requisite
machinery, especially the two-electron integrals, and evaluate the performance
of the RPA, time-dependent Hartree-Fock (TDHF), and configuration interaction with
single excitations (CIS), which is the Tamm-Dancoff approximation (TDA) to TDHF.
For accessible system sizes, we compare our results to higher-level
equation-of-motion coupled-cluster theory with single and double excitations
(EOM-CCSD)~\cite{Emrich1981,Koch1990,Bishop1991,Stanton1993,Bartlett2007,Krylov2008},
which has recently been used to characterize plasmons in the bulk
UEG~\cite{Lewis2019}.

The intraband excitation in our model evolves from confinement-dominated
and single-particle in character, to excitonic, to plasmonic, when the number of electrons is increased.
By analyzing the underlying physics and comparing to higher-level
EOM-CCSD, we argue that the evolution of
the excitation is most accurately described by TDHF; the RPA and CIS are
distinct approximations to TDHF and are only able to correctly describe plasmons
and excitons, respectively. The transition from a confinement-dominated to
plasmon-like excitation can be driven by both the number and density of
electrons, and we find that for a fixed number of electrons, increasing the
density actually \textit{decreases} the plasmonic character of the excitation,
opposite to the prediction of noninteracting models~\cite{Jain2014}, due to
increasing quantum confinement. The combination of local density functional
theory (DFT) and the RPA, a popular technique in
literature\cite{prodan2003structural,Zhang2014,Zhang2017,Ipatov2018}, predicts results that are
similar to TDHF, but due to a cancellation of errors.
The character
of the excitation across the entire range of number of electrons and density, in
particular the intermediate excitonic state, can only be described by properly
accounting for the attractive electron-hole (exchange) interaction, similar
to the situation in molecules or semiconductors.

\section{Theory}
\subsection{Model}
\label{ssec:theory_model}

As a model of a doped nanoparticle, we treat the conduction band electrons as
particles in an infinitely deep spherical well, where the atomic details of the
nanoparticle are represented by the effective mass, radius, and
dielectric constant. When used, the dielectric constant approximates the
effect of the ignored valence electrons and higher excitations, which screen
the Coulomb interaction.  A more sophisticated model of the surface would
include a finite or stepped barrier, however we do not
expect our qualitative conclusions to be sensitive to the details of the
surface.  Furthermore, we neglect the dielectric
contrast with the environment, which alters both the single particle band gap
and the optical gap, although these two effects partially cancel in low-dimensional
semiconductors~\cite{Cho2018}. We will treat the two-body
Coulomb interactions between the conduction
band electrons, which represents the focus of this work and goes beyond simple
models of noninteracting electrons under confinement.
We note that this ``jellium sphere'' model and various levels of theory have also been used
to describe the structure and excitations of nuclei (including the giant dipole resonance~\cite{ring2004nuclear})
and especially the optical properties of alkali clusters~\cite{Puska1985,Guet1992,deHeer1993,Brack1993,Koskinen1994PRB,Madjet1995,Patterson2019}.

In first quantization and atomic units, the total $N$-electron Hamiltonian is
\begin{equation}
\begin{split}
H &= \sum_{n=1}^{N} \left[-\frac{1}{2m^*}\nabla_n^2 + v({\bm{r}}_n)\right]
    + \sum_{n=1}^{N}\sum_{m<n} \frac{1}{|{\bm{r}}_m-{\bm{r}}_n|},
\end{split}
\end{equation}
where $m^*$ is the effective mass of the conduction band electrons and
$v(r<R) = 0$ and $v(r\ge R) = \infty$. Although charge neutrality would imply an
additional $r$-dependent harmonic potential, here we neglect this potential
because many experimental procedures for nanoparticle doping (e.g.~photodoping
with hole scavengers~\cite{Schimpf2014}) do not preserve charge neutrality.

We use an orthogonal one-particle basis of eigenfunctions of the
one-electron part of the above Hamilonian, corresponding to the well-known
particle-in-a-sphere (PIS),
\begin{equation}
    \phi_{nlm}(r,\theta,\phi) = N_{nl}^{-1} j_{l}(\alpha_{nl}r) Y_{lm}(\theta,\phi) = R_{nl}(r) Y_{lm}(\theta,\phi),
\end{equation}
where $r$ is the radial coordinate, $\alpha_{nl} = k_{nl} /R$,
$k_{nl}$ is the $n$th zero of the spherical Bessel function $j_l(r)$, and the
normalization constant is $N_{nl} = \sqrt{R^3/2} \left|j_{l+1}(k_{nl})\right|$.
Each orbital is characterized by three quantum numbers $n,l,m$, with the limits
$n\geq 1$, $l\geq 0$, and  $m=-l,...,l$. The noninteracting orbital energies are given by
\begin{equation}
    \varepsilon_{nlm} = \frac{k_{nl}^2}{2m^*R^2},
\label{eq:oneeen}
\end{equation}
and are $(2l+1)$-fold degenerate. Unlike the hydrogen atom, the orbital energies
of the PIS are not degenerate with respect to the principle
quantum number $n$, but are $m$-fold degenerate for a given $n$ and $l$.

In this orthogonal single-particle basis, the interacting, second-quantized
Hamiltonian is
\begin{equation}
H = \sum_{p} \varepsilon_p a_p^\dagger a_p
    + \frac{1}{2}\sum_{pqrs} \langle pq|rs\rangle a_p^\dagger a_q^\dagger
        a_s a_r,
\end{equation}
where the indices $pqrs$ run over spin-orbitals, i.e.~$p=(n,l,m,\sigma)$
and the two-electron integrals are given by
\begin{equation}
\langle pq|rs \rangle
    = \int d{\bm{x}}_1 \int d{\bm{x}}_2 \phi_p^*({\bm{x}}_1) \phi_q^*({\bm{x}}_2) r_{12}^{-1} \phi_r({\bm{x}}_1) \phi_s({\bm{x}}_2),
\end{equation}
where ${\bm{x}}=({\bm{r}},\sigma)$ is a combined space and spin variable.

We study nanoparticles containing $N=2$, 8, 18, 32, 50, 72,
and 98 electrons, which correspond to closed-shell solutions of restrictedf
Hartree-Fock (RHF). For these system sizes, we find that the RHF solution
occupies the orbitals 1s,1p,1d,...~up to 1$l_\mathrm{max}$, and therefore
these closed-shell fillings correspond to
$N=2\sum_{l=0}^{l_\mathrm{max}}(2l+1)=2(l_\mathrm{max}+1)^2$, where
the factor of 2 accounts for spin.

The naive way to grow the basis set is to add PIS orbitals based on increasing
energy; however, the RHF orbitals are pure eigenfunctions of $l$, so increasing
the number of basis functions $n$ for each $l=0...l_\mathrm{max}$ is sufficient
to converge the ground state calculation. In order to capture singly excited
states we add an
additional shell $l_\mathrm{max}+1$ based on the dipole selection rule $\Delta
l=\pm 1$. The rapidly increasing degeneracy of the basis
functions limits the number of electrons to 98, which we converge with 483 basis
functions ($n_\mathrm{max}$ of $[10,9,9,8,8,7,7,7]$ for
$l=0...l_\mathrm{max}+1=0...7$) requiring about 50 GB to store
the two-electron integrals.

Expressions for the two-electron integrals $\langle pq|rs\rangle$, which are not
analytic but can be reduced to two-dimensional quadrature along the radial axis,
are given in the Appendix. All
electronic structure calculations are performed by defining a custom Hamiltonian
for use in the PySCF software package~\cite{Sun2018}.

\subsection{Excited states}

We focus on quantum chemical single-excitation theories due to their
favorable $O(N^4)$ scaling with system size, which makes them practical for future
atomistic studies.
Specifically, we consider excited states of the form
\begin{equation}
|\Psi_n\rangle = \sum_{ai}
    \left[X_{ai} a_a^\dagger a_i + Y_{ai} a_i^\dagger a_a \right]
    |\Psi_0\rangle
\label{eq:singleexitedstate}
\end{equation}
where here and throughout $i,j,k,l$ and $a,b,c,d$ index occupied and unoccupied
HF orbitals, and $X_{ai}$ and $Y_{ai}$ correspond to coefficients for the
excitation and deexcitation of an electron from orbital $i$ to $a$,
respectively. The deexcitation operator implies that the ground state
$|\Psi_0\rangle$ is potentially correlated, though unspecified.
The amplitudes
$X_{ai}$ and $Y_{ai}$  are obtained from the eigenvalue
problem~\cite{ring2004nuclear,Scuseria2008}
\begin{equation}
\label{eq:rpa}
\left(
\begin{array}{cc}
\mathbf{A}   & \mathbf{B} \\
-\mathbf{B}^* & -\mathbf{A}^*
\end{array}
\right)
\left(
\begin{array}{c}
\mathbf{X} \\
\mathbf{Y}
\end{array}
\right)
=
\left(
\begin{array}{c}
\mathbf{X} \\
\mathbf{Y}
\end{array}
\right) \mathbf{\Omega},
\end{equation}
where $\mathbf{\Omega}$ is a diagonal matrix of excitation energies.
The single excitation theories considered here correspond to specific
choices of the $\mathbf{A}$ and $\mathbf{B}$ matrices.
The most ``complete'' theory is TDHF, for which
\begin{subequations}
\begin{align}
A_{ai,bj} &= (\varepsilon_a - \varepsilon_i) \delta_{ab}\delta_{ij} + \langle ib || aj \rangle, \\
B_{ai,bj} &= \langle ij || ab \rangle,
\end{align}
\end{subequations}
and the antisymmetrized integrals are defined as $\langle pq||rs\rangle =
\langle pq|rs\rangle - \langle pq|sr\rangle$. The RPA is obtained by neglecting
this antisymmetrization (consistent with time-dependent Hartree theory), and retains only
the ``direct'' Coulomb matrix elements.
The neglected ``exchange'' Coulomb matrix elements (also sometimes called ``direct electron-hole
interactions'') are responsible for the
electron-hole attraction and the formation of bound excitons.
The Tamm-Dancoff approximation (TDA) corresponds to $\mathbf{B}=0$,
which neglects potential correlation in the ground state.
When applied to TDHF and the RPA, the TDA
leads to theories we will refer to
as CIS and the RPA(TDA).

The RPA is the minimal theory necessary for the description
of plasmons. For the UEG, the RPA predicts a plasmon dispersion $\Omega(q)$ that
has the correct
long-range limit, $\Omega(q\rightarrow 0) = \Omega_\mathrm{p}$, where
$\Omega_\mathrm{p} = \sqrt{4\pi \rho}$ is the classical plasma frequency~\cite{PinesNozieres,GiulianiVignale}. The TDA
(including CIS) predicts a collective excitation, but one whose energy
unphysically diverges in the long-range limit,
$\Omega^\mathrm{TDA}(q\rightarrow 0) \rightarrow \infty$.
Here, we will see the same behavior in the $R\rightarrow \infty$ limit,
highlighting the failure of the TDA for large plasmonic nanoparticles.

For comparison, we also present results obtained at lower and higher levels
of theory, where available. At a low level, we consider simply the orbital
energy differences from the noninteracting (PIS) and mean-field (HF) theories,
i.e.~the difference in energies of the occupied and unoccupied orbitals with the largest
transition dipole matrix element.
For HF, these orbitals are always the highest occupied and lowest unoccupied molecular orbitals (HOMO and LUMO).
At a high level, we present results from
equation-of-motion coupled-cluster theory with single and
double excitations (EOM-CCSD) up to 50 electrons. Importantly, we note that EOM-CCSD
includes ground-state correlation and excited-state double excitations,
which both contribute to screening in an effective single-excitation
theory like the Bethe-Salpeter equation. We use a basis set that adds PIS
orbitals based on increasing energy to converge the EOM-CCSD calculations.

Since we are interested in the absorption properties of doped nanoparticles, we calculate the
dynamical polarizability
\begin{equation}
\alpha(\omega) = \sum_{m} |\langle \Psi_0|\hat{\mu}|\Psi_m\rangle|^2
\delta(\omega - \Omega_m)
\label{eq:suscept}
\end{equation}
where $\hat{\mu} = \sum_{n=1}^{N} \hat{r}_n = \sum_{pq} r_{pq} a_p^\dagger a_q$
is the dipole operator.

\section{Results and Discussion}

We study a model nanoparticle with $m^* = 0.28$; which is characteristic
of the conduction band of ZnO, whose plasmonic properties under doping have
been experimentally studied~\cite{Schimpf2014}.  As discussed in
Sec.~\ref{ssec:theory_model}, we consider systems containing 2, 8, 18, 32, 50,
72, and 98 electrons and study experimentally
relevant densities from $\rho=1.4\times 10^{20}$ to $1\times 10^{22}$~cm$^{-3}$,
which defines the radius $R$ for a given number of electrons $N$.  In this range,
the radii of the nanoparticles are on the order of 1-10~nm.

\subsection{Spectral properties and peak position}

\begin{figure}
\includegraphics[scale=1]{Manuscript_Fig1.png}
\caption{Optical and electronic structure properties of doped nanocrystals, for
$m^* = 0.28$ at density $\rho=1.4\times 10^{20}$~cm$^{-3}$ (left) and $1\times
10^{22}$~cm$^{-3}$ (right). \textbf{Top:} The noninteracting (left column) and
mean-field (RHF, right column) density of states for 2, 8, 32, and 98 electrons.
The manifold of unoccupied states is truncated due to
the finite basis set. \textbf{Middle:} The TDHF absorption spectrum (red)
density of excited states (grey, filled); the energy of the HF gap is indicated by a black
solid peak.
\textbf{Bottom:} Energy of
the peak absorption versus radius, for 2, 8, 18, 32, 50, 72, and 98 electrons.
The shaded region marks the bulk $\Omega_\mathrm{p}=\sqrt{4\pi \rho/m^*}$ and
Mie+Drude theory $\Omega_\mathrm{LSPR}=\Omega_\mathrm{p}/\sqrt{3}$ plasmon frequency.}
\label{fig:1}
\end{figure}

First we analyze the single-particle orbital energies.  The top two panels of Figure~\ref{fig:1} show the the density of
states (DOS) predicted by the noninteracting (PIS, left) and mean-field (HF, right) level of theory,
where blue lines indicate occupied orbitals and red lines indicate unoccupied orbitals.  Results are
shown for
$N=2,8,32,$ and 98 electrons (bottom to top) at the smallest and largest densities considered
here, $1.4\times 10^{20}$~cm$^{-3}$ and $1\times 10^{22}$~cm$^{-3}$ (left and right), the former
of which corresponds to the ZnO system.
In both cases, the energy spacings are reduced with increasing $R$ due to the
scaling of the kinetic energy term, $1/(2m^*R^2)$.  In the nonineracting
results, the energy of the HOMO (roughly the Fermi energy) decreases with
increasing $R$, which is typical particle-in-a-box type physics.  However in the
HF results, the spectrum is further compressed and the energy of the HOMO is
shifted to \textit{higher} energies due to the mean-field effects of the Coulomb
interaction.

In the middle two panels of Figure~\ref{fig:1}, we show the RHF (black) and TDHF absorption
spectrum (red), and the density of excited states $D(\omega) = \sum_m
\delta(\omega-\Omega_m)$  (grey, filled) for the same systems as above. From
the spectrum of excited states, which is becoming gapless in the limit of large
$N$, we observe a dominant peak in the absorption spectrum that arises from a
state that is typically \textit{not} the lowest in energy.  This redistribution
of oscillator strength from a low-energy continuum into a single high-energy
state is reminiscent of the plasmon peak in the dynamical structure factor of
the UEG~\cite{ring2004nuclear,Lewis2019}; the bright state is only the lowest in energy for $N=2$.
Henceforth, we focus on the excitation energy of the dominant bright peak, which
is the energy of the excited state with the largest transition dipole matrix
element.  The spectra in Figure~\ref{fig:1} also highlight a computational
challenge, as the strongly absorbing state of interest is buried in the interior
of the eigenvalue spectrum of the Hamiltonian.

In the bottom two panels of Figure~\ref{fig:1}, we compare the energy of this dominant absorption  peak versus radius, for these same two densities, as predicted by various theories.
We present results for TDHF, the RPA, CIS, RPA(TDA), and EOM-CCSD, along with
the noninteracting and HF transitions for comparison.
The EOM-CCSD result, the most accurate solution here, demonstrates the qualitative
evolution of the excitation, which we separate into three regions. At small $R$, the excitation energy
is dominated by the kinetic (confinement) energy, and scales with the band gap. At intermediate
$R$, the EOM-CCSD result is \textit{below} the HF gap and exhibits a minimum,
which we attribute to the formation of intraband excitons.  At large $R$, the
excitation energy increases with increasing $R$ and goes \textit{above} the HF
gap before reaching a plateau near the classical plasma energy.  This latter
behavior is consistent with the $R\rightarrow \infty$ limit of our model,
i.e.~the UEG, which has a well-known $q\rightarrow 0$ plasmon
at the classical plasma frequency.

We observe from Figure~\ref{fig:1} that TDHF is the most accurate
single excitation theory. For a given density, as the
number of electrons (or radius) increases, the energy of the TDHF excitation
crosses from below to above the HF gap, which follows the EOM-CCSD result. This behavior is consistent with the physics
embodied in TDHF, which contains the ingredients necessary to capture
the three regimes described above, i.e.~confinement-dominated excitations, excitons
(due to exchange integrals, i.e.~direct electron-hole interactions), and plasmons (due to the nonzero $\mathbf{B}$
matrix). These claims are validated by comparison with the ``approximations''
to TDHF, which contain only a subset of these ingredients: the RPA correctly predicts the evolution of the excitation to the
classical plasmon (always above the HF gap), but cannot lower the energy of
the excitaton at any radius due
to the lack of electron-hole interactions, and CIS contains the
exchange Coulomb interaction but neglects the $\mathbf{B}$ matrix,
so it is accurate in
describing bound excitons at small radii, but fails to correctly describe the
plasmonic state at large radii. In particular,
CIS predicts an excitation energy which goes above the HF gap and diverges at
large radius, which is precisely analogous to its $q\rightarrow 0$ behavior
in the UEG.
Finally, the RPA(TDA) result is not accurate at any radius.

Despite the good
overall agreement of TDHF with EOM-CCSD at high electron densities, the TDHF solution suffers at low
density and large $R$ (always underpredicting the EOM-CCSD result), which is
precisely where electron correlation is expected to be strongest.
In this regime, one may anticipate the onset of Wigner crystallization and multireference
character~\cite{Thompson2002,Jung2003,Thompson2004,Thompson2005}. While this physics could be approximately addressed via
spin-symmetry breaking, we do not pursue this approach here.

\subsection{Characterization of the excited state}

\begin{figure}
\includegraphics[scale=1]{Manuscript_Fig2.png}
\caption{
Summary of TDHF and RPA results for 2, 8, 18, 32, 50, 72, and 98
electrons at densities logarithmically spaced from
\SI{1.4E20}{cm^{-3}} to \SI{1E22}{cm^{-3}}. \textbf{Top:} Energy of the peak
absorption versus radius. The shaded regions mark the bulk
$\Omega_\mathrm{p}=\sqrt{4\pi \rho/m^*}$ and Mie theory
$\Omega_\mathrm{LSPR}=\Omega_\mathrm{p}/\sqrt{3}$ plasmon frequency for each
density. Dashed lines connect results obtained with 2, 18, 50, and 98 electrons.
\textbf{Bottom:} The sum of the deexcitation coefficients,
$\sum_{ai}|Y_{ai}|^2$, for the peak absorption state versus radius, calculated with
TDHF (left), and the RPA (right). The data are connected by lines of constant
density but changing particle number (solid) and by lines of constant particle
number but changing density (dashed).}
\label{fig:2}
\end{figure}

We now discuss a microscopic characterization of the wavefunction of the
strongly absorbing excited state. The top two plots of Figure~\ref{fig:2}
summarize the energy of the peak absorption for TDHF and the RPA, and the bottom
two plots of Figure~\ref{fig:2} plot the sum of the $Y$ coefficients for the
peak absorption, which is a measure of the plasmonic character of the
excitation. The plasmonic character is sensitive to both the absolute number of
electrons and the density, and the emergence of a plasmon is not solely a
density-driven transition. For a fixed number of electrons, the plasmonic
character decreases with increasing density, opposite to the noninteracting
result \cite{Jain2014}; at higher densities, a greater number of electrons is
needed to reach the same plasmonic character as a fewer number of electrons in a
lower density system. This decrease in plasmonic character at smaller $R$ is due
to the quantum confinement of the one-electron kinetic energy that scales as
$1/(2m^*R^2)$, and further exacerbated by the strong electron-hole interaction
(compare TDHF to the RPA). When the energy of the peak absorption is compared to
the plasmonic character, we find the surprising result that even though the
energy lies within the range of classical plasma frequencies, the excitation can
be far from plasmonic, and is instead excitonic or single-particle-like.

Figure~\ref{fig:3} is a series of plots of
the induced charge density,
\begin{equation}
\begin{split}
\delta \rho({\bm{r}},\omega)&=\int d{\bm{r}}' \chi({\bm{r}},{\bm{r}}',\omega)v_\mathrm{ext}({\bm{r}}',\omega) \\
&= \sum_n \frac{2\Omega_n}{\omega^2-\Omega_n^2}\sum_{ai}|X_{ai}^n+Y_{ai}^n|^2\langle i|\hat{z}|a\rangle\phi_a({\bm{r}})\phi_i({\bm{r}}),
\end{split}
\end{equation}
where $\chi$ is the density-density linear response function, for an external electric field
oriented along the $z$-axis.
We evaluate the induced charge density at the peak absorption energy for the
noninteracting, HF, TDHF, and CIS levels of theory. The doping density of
\SI{1.4e20}{cm^{-3}} and
2 (top), 8 (middle), and 98 (bottom) electrons correspond to nanoparticle radii
of 1.5~nm, 2.4~nm, and 5.5~nm, and excitation character of confined, excitonic,
and plasmonic. As the number of electrons increases, the induced density
concentrates at the surface, in agreement with classical theory.  The addition
of mean-field interactions (HF) repels the non-interacting induced density to
the surface. The induced density at the mean-field level is not changed under configuration
mixing, because the HOMO-LUMO transition has the largest dipole matrix
element $\langle i | \hat{z} | a\rangle$ and this configuration contributes with
a large weight in the bright state.  In other words, the excited-state wavefunctions are all
qualitatively similar, despite predicting very different energies.

\begin{figure}
\includegraphics[scale=1]{Manuscript_Fig3}
\caption{Characterization of the strongly absorbing excited state of nanoparticles for a doping
density of density \SI{1.4e20}{cm^{-3}}. \textbf{Top:} The induced charge
density for 2 ($R=1.5$~nm), 8 ($R=2.4$~nm), and 98 electrons ($R=5.5$~nm), which
roughly correspond to a single-particle, excitonic, and plasmonic excitations.
Results are presented at four levels of theory (noninteracting, RHF, TDHF, and CIS).
\textbf{Bottom:} The many-body contributions to the energy of the peak
absorption versus radius as given in Eq.~\ref{eq:ep}.
}
\label{fig:3}
\end{figure}

In light of the above observation, we next quantify the contributions of the various singles theories to
the total excitation energy,
\begin{align}
\begin{split}
\Omega_n &= 2\sum_{ai} (\varepsilon_a-\varepsilon_i) |X_{ai}+Y_{ai}|^2 \\
&\hspace{1em} +4\sum_{abij}\langle ib|aj\rangle\left(X_{bj}+Y_{bj}\right)\left(X_{ai}+Y_{ai}\right) \\
&\hspace{1em} -2\sum_{abij}\langle ib|ja\rangle\left(X_{ai}X_{bj}+Y_{ai}Y_{bj}\right) \\
&\hspace{1em} -2\sum_{abij}\langle ij|ba\rangle\left(X_{ai}Y_{bj}+Y_{ai}X_{bj}\right).
\end{split}
\label{eq:ep}
\end{align}
The second line in Eq.~(\ref{eq:ep}), containing the ``direct'' $\langle ib|aj\rangle$ integrals that describe the interaction energy between single-particle excitations, has been termed the ``plasmonicity''~\cite{Zhang2017} and is nonzero for TDHF, RPA, CIS, and RPA(TDA). The third line, containing the ``exchange'' $\langle ib|ja\rangle$ integrals from the $\mathbf{A}$ matrix, provides
a measure of the exciton binding energy and is nonzero for TDHF and CIS.  Finally, the fourth line,
containing the ``exchange'' $\langle ij|ba\rangle$ integrals from the $\mathbf{B}$ matrix has no
classical interpretation and is only nonzero for TDHF elements.

The bottom of Figure~\ref{fig:3} plots the above contributions versus radius for a doping
density of \SI{1.4e20}{cm^{-3}}, for all four singles theories. CIS and the
RPA(TDA) do not account for ground state correlations ($\mathbf{B}=0$) and
predict a diverging plasmonicity. The RPA adds a nearly constant plasmonicity at
\textit{all} $R$, such that the excitation energy within the RPA is only modulated
by the mean-field gap, which is consistent with the prediction of a higher energy
collective excitation that is a combination of degenerate single-particle
excitations\cite{ring2004nuclear}. The RPA and RPA(TDA) overpredict the excitation
energy at small radii due to the lack of the excitonic interaction (see Figure~\ref{fig:2}, bottom).
CIS fails at large $R$, but when quantum confinement
dominates at small $R$, both CIS and TDHF correctly describe the bound excitonic
state: notably, for the most confined case of two electrons, the exciton binding
energy lowers the excitation energy by roughly 1.5~eV.
The exciton binding energy decreases with increasing $R$
as the electron-hole spatial overlap becomes smaller.
The nonclassical exchange contribution from the $\mathbf{B}$ matrix makes a small
positive contribution at all values of $R$, never exceeding 0.2~eV.

TDHF mixes the correct large $R$ limit of the RPA and the small $R$ limit of
CIS. Importantly, TDHF inherits the RPA prediction of a relatively constant
plasmonicity at all radii, Figure~\ref{fig:3} bottom left. Thus, at least at this density,
there is no clearly distinguishable single-particle excited state, and the
exciton binding energy and plasmonicity make non-negligible contributions at all
values of $R$.

\subsection{Schematic model}

We can understand the dependence of the energy of the excitation on radius by
considering a schematic model of the RPA and RPA(TDA) theories, i.e.~those without
antisymmetrized integrals. This treatment is motivated by the presentation in Ref.~\onlinecite{ring2004nuclear}.
The results of the schematic model at all levels of theory will be compared to the
numerical results at four densities in Figure~\ref{fig:4}.

We consider the
factorization of the direct two-electron integrals
$\langle ib|aj\rangle \approx \lambda \rho_{ai} \rho_{bj}$, with $\lambda > 0$.
In this approximation, the RPA(TDA) equation for an excited state $n$ can be
written as
\begin{equation}
\left[ \Omega_n-(\varepsilon_a-\varepsilon_i) \right] X_{ai}^n
= \lambda \rho_{ai} \sum_{bj} \rho_{bj} X_{bj}^n,
\end{equation}
which leads to~\cite{ring2004nuclear}
\begin{align}
\frac{1}{\lambda} = \sum_{ai}
    \frac{\rho_{ai}^2}{\Omega_n - (\varepsilon_a-\varepsilon_i)}.
\end{align}
The latter equation can be
solved graphically for $\Omega_n$, leading to a number of single-particle
excitations with energies approximately given by $\varepsilon_a - \varepsilon_i$
and one higher-energy collective excitation (the plasmon).
For illustrative purposes, we consider the subspace containing only the
HOMO and LUMO (each potentially degenerate), such that $\varepsilon_a-\varepsilon_i = \varepsilon$, which yields
for the plasmon state
\begin{align}
\Omega = \varepsilon + \lambda \sum_{ai} \rho_{ai}^2 \approx \varepsilon + \lambda N_\mathrm{trans} \bar{\rho}^2,
\end{align}
where $N_{\mathrm{trans}}$ is the number of transitions
and $\bar{\rho}$ is an average quantity.

We now seek to understand the $R$ dependence of this excitation energy at
\textit{fixed} density.  For all of the system sizes studied here, the
closed-shell RHF solution is only stable for those configurations for which all
occupied orbitals have $n=1$, i.e.~1s$^2$, 1p$^6$, 1d$^{10}$, 1f$^{14}$, ...,
and so on.  The angular momentum of the HOMO is thus defined by the number of
electrons, $N = 2(l_\mathrm{max}+1)^2$ or
$l_\mathrm{max} = \sqrt{N/2}-1$ at fixed density.  To a good approximation, we
find that the zeros of the spherical Bessel functions can be written as $k_{nl}
= n\pi + 1.32l$ (in particular, the value 1.32 is empirically better than the
asymptotic value $\pi/2$).  This yields a noninteracting band gap from $1l$ to
$1(l+1)$ of
\begin{equation}
\varepsilon_\mathrm{g}^{\mathrm{NI}}(R)
    = \frac{6.55}{2m^*R^2} + \frac{5.04 \sqrt{\rho}}{2m^* \sqrt{R}}.
\end{equation}
At the HF level, we find that the form of the one-electron
contribution to the band gap is very similar and most significantly modified by
the exchange contribution, which we model with the form $1/R$ for all densities,
$\varepsilon_\mathrm{g}^\mathrm{HF}(R)
= \varepsilon_\mathrm{g}^{\mathrm{NI}}(R) + 1/R$.

The HOMO and LUMO each have a degeneracy proportional
to $l(l+1)$ which, combined with the dipole selection rule
$\Delta m = 0,\pm1$, leads to a number of transitions
$N_\mathrm{trans} \propto l_\mathrm{max} \propto \sqrt{\rho R^3}$. The
Coulomb interaction has a scaling $\lambda(R) \propto R^{-1}$.
Therefore, within the RPA(TDA), the schematic model predicts an excitation
energy
\begin{equation}
\Omega^{\mathrm{RPA(TDA)}}(R)
    = \varepsilon_\mathrm{g}^{\mathrm{HF}}(R) + c \sqrt{\rho R},
\end{equation}
where we take $c=1$. At small $R$, the excitation energy of a doped
nanoparticle is given by the kinetic-energy-determined band gap; at large $R$,
the excitation energy diverges due to the Coulomb interaction. This divergence
is unphysical and analogous to the behavior of the TDA in the $q\rightarrow 0$
limit of the three-dimensional uniform electron gas.

The divergence is fixed in the full RPA, which is given in the schematic model by
\begin{align}
\frac{1}{\lambda} = \sum_{ai} \rho_{ai}^2
    \frac{2(\varepsilon_a-\varepsilon_i)}
        {\Omega_n^2 - (\varepsilon_a-\varepsilon_i)^2}.
\end{align}
The same approximations as above for a spherical nanoparticle leads to
\begin{equation}
\Omega^{\mathrm{RPA}}(R)
    = \sqrt{\left(\varepsilon_\mathrm{g}^\mathrm{HF}(R)\right)^2
                + 2 \sqrt{\rho R} \varepsilon_\mathrm{g}^{\mathrm{HF}}(R)}.
\end{equation}
This excitation energy has the same kinetic-energy-determined band gap at small
$R$, but now has a finite $R\rightarrow \infty$ limit,
$\Omega^{\mathrm{RPA}}(R\rightarrow \infty)
= \sqrt{5.04\rho/m^*}$.  Importantly, the schematic model recovers the exact
limiting form of the classical plasmon energy, up to constants of order 1,
i.e.~$\sqrt{5.04} \approx 2.24$ compared to $\sqrt{4\pi}\approx 3.54$.

\begin{figure}[t]
\includegraphics[scale=1]{Manuscript_Fig4.png}
\caption{The results of the approximate solutions of the schematic model (solid lines) compared
to the numerical results (symbols) at four different
densities.}
\label{fig:4}
\end{figure}

Before continuing on to theories with exchange (CIS and TDHF), we first
turn to an analysis of the excitation coefficients. Again within the degenerate
schematic model, the $X$ coefficients in RPA(TDA) are
\begin{align}
    X_{ai}^{\mathrm{RPA(TDA)}} = \frac{1}{\sqrt{C}} \times \rho_{ai}
\end{align}
where $C = \sum_{ai} |\rho_{ai}|^2$ is a normalization constant,
and in the RPA are
\begin{align}
X_{ai}^{\mathrm{RPA}} &= \frac{1}{\sqrt{C}}
    \times \frac{\rho_{ai}}{\Omega - \varepsilon_\mathrm{g}} \\
Y_{ai}^{\mathrm{RPA}} &= -\frac{1}{\sqrt{C}}
    \times \frac{\rho_{ai}}{\Omega + \varepsilon_\mathrm{g}},
\end{align}
where $\Omega$ is the energy of the collective state and
\begin{equation}
C = 4\sum_{ai}|\rho_{ai}|^2 \varepsilon_\mathrm{g} \Omega_{\mathrm{P}}
\end{equation}
We note that as $\varepsilon_\mathrm{g}\rightarrow 0$ and the
model becomes more metallic, the excited
state is a plasmon and the $X$ and $Y$ coefficients become equal in
magnitude.
In contrast, as the gap $\varepsilon_\mathrm{g}$ increases, the $X$
coefficients dominate and $Y_{ai}\rightarrow 0$.
This behavior is observed numerically, as shown in Figure~\ref{fig:2}.
The sum of deexcitation coefficients for the collective state,
$\sum_{ai} |Y_{ai}|^2$, is therefore another measure of plasmonic character.
From Figure~\ref{fig:2}, it is clear to see that this character increases
with increasing $R$ at fixed $N$ or increasing $N$ at fixed $R$, because
both situations correspond to reducing the gap $\varepsilon_\mathrm{g}$. By this
measure, a low plasmonic character does not imply that the excitation is
single-particle-like, because the excitation can still be delocalized over the
$X_{ai}$ coefficients.

The inclusion of antisymmetrized integrals in CIS and TDHF prevents an
analytic treatment of the schematic model because the integral factorization
$\langle ia|jb\rangle \approx \lambda \rho_{ij} \rho_{ab}$ does not facilitate
the solution of the eigenvalue problem.
However, the largest-in-magnitude element $\langle ia|ia\rangle$ can be included
exactly as it just shifts the HF band gap
$\varepsilon_\mathrm{g}\rightarrow \varepsilon_\mathrm{g} - \langle ia|ia\rangle$.
Within the schematic model, we take all such excitonic Coulomb integrals to be equal
and obeying the scaling $\langle ia|ia\rangle = c^\prime/R$, with $c^\prime = 1.4$.
This approximation gives the CIS and TDHF excitation energies as
\begin{equation}
\Omega^{\mathrm{CIS}}(R)
    = \varepsilon_\mathrm{g}^{\mathrm{HF}}(R) + \sqrt{\rho R} - \frac{1.4}{R}
\end{equation}
\begin{equation}
\Omega^{\mathrm{TDHF}}(R)
    = \sqrt{\left(\varepsilon_\mathrm{g}^\mathrm{HF}(R)-\frac{1.4}{R}\right)^2
                + 2 \sqrt{\rho R} \left(\varepsilon_\mathrm{g}^{\mathrm{HF}}(R)-\frac{1.4}{R}\right)}.
\end{equation}
We note that the remaining excitonic Coulomb integrals $\langle ia|jb\rangle$ could be
included via perturbation theory, though we do not pursue this here.

In Figure~\ref{fig:4}, we show the performance of the analytic solutions of this schematic model
compared to our numerical calculations. The noninteracting, HF, RPA, and
RPA(TDA) schematic models fit the data remarkably well, with the only
parameterizations being the fit of the spherical Bessel function zeros and the
proposed $1/R$ form of the HF exchange contribution to the gap. The CIS and TDHF
schematic models fit the data reasonably well, with the main difficulty being
the estimation of the contribution of the antisymmetrized two-electron
integrals. Importantly, the schematic model captures the overall behaviour of
our calculations based upon simple scaling arguments on the density and radius,
and supports the physical interpretations given to the various single-excitation theories.

\subsection{Effect of screening and comparison to DFT}

\begin{figure}
\includegraphics[scale=1]{Manuscript_Fig5.png}
\caption{
The energy of the peak absorption versus radius, parameterized to the
ZnO system ($\rho=1.4\times10^{20}$, $m^*=0.28$,
$\varepsilon=3.72$). \textbf{Left:} theories with
exchange: RHF gap and TDHF excitation energy
compared to the experimental results of Ref.~\onlinecite{Schimpf2014}.
\textbf{Right:} theories without exchange: LDA gap and
the DFT+RPA excitation energy, again compared to to experiment. The inset at right
is the plasmonicity energy from TDHF (red) and DFT+RPA (yellow).}
\label{fig:5}
\end{figure}

In our model, the effect of screening from the ignored valence electrons can be approximately
accounted for by a scaling of the Coulomb interaction,
$r_{12}^{-1}\rightarrow(\varepsilon r_{12})^{-1}$. We can approximate ZnO nanoparticles
by setting $\varepsilon=3.72$. Figure~\ref{fig:5} plots the energy of the
intraband absorption as a function of radius for this ``ZnO'' nanoparticle, at
the experimental doping density of \SI{1.4E20}{cm^{-3}}, against the
experimental results from Ref.~\onlinecite{Schimpf2014}. We compare TDHF (left) with DFT plus the RPA (right), where
the RPA excited state is calculated from the solution of the restricted
Kohn-Sham equations with the local density approximation (LDA). The LDA+RPA
approach has been applied previously to our model~\cite{Varas2016,Zhang2017} and has found success in describing the
experimental ZnO results~\cite{Goings2014,Ipatov2018,Gerchikov2018}.

From Figure~\ref{fig:5}, both TDHF and DFT+RPA compare favorably with the
experimental result, although they embody different physical effects as explored
in this paper. As usual, the LDA gap is smaller than the RHF gap, due to the
latter's treatment of exchange. The inset on the right of Figure~\ref{fig:5}
plots the sum of all contributions given in Eq.~(\ref{eq:ep}) \textit{except}
for the orbital-energy contribution, from LDA+RPA and RHF+TDHF. The RPA only acts to increase
the small LDA gap at all $R$; by contrast, the TDHF lowers the large RHF gap at
small $R$ (due to excitonic effects) and increases the RHF gap at large $R$ (due
to plasmonic effects). Therefore, ignoring exchange in the ground state and
excited state, as in the LDA+RPA approach, can produce an accurate result due to
cancellation of errors\cite{Jung2004}.

\section{Conclusions}

In summary, we provide a fully quantum mechanical study of a confined, interacting electron gas
as a model for doped semiconductor nanoparticles.  We observe strongly absorbing excited states
whose wavefunction character can be classified as single-particle-like (confinement dominated),
excitonic, or plasmonic.  Within the framework of the most computationally affordable single-excitation
theories, only TDHF is capable of capturing the qualitative behavior at all studied densities and particle
sizes.  We also present a schematic model of the strongly absorbing excited state that reproduces
the $R$-dependence observed in our simulations.

Our model is simple in order to focus on the essential features of electronic
interactions in the excited states of confined systems.  The model neglects atomistic details as well as surface,
ligand, or solvent effects.  The model is also ignorant of the doping mechanism and neglects the
atomic defect potential that is introduced by impurity doping (but not by electron transfer or photodoping).
Nonetheless, the results of our calculations argue strongly against the interpretation or simulation of
doped nanoparticle spectra based on single-particle transitions between orbitals, and we propose an interpretation
of intraband excitons as the primary excitations at low doping or small nanoparticles.

Looking forward towards atomistic or tight-binding~\cite{Pi2013} simulations, our work has two important ramifications. First, the
TDA fails spectacularly and should be avoided in all simulations seeking to address
the possibility of plasmonic excitations.  Second, the retention of attractive electron-hole ``exchange'' integrals
is essential for an accurate wavefunction description of excitonic states.
With these criteria in mind,
we suggest that the most promising and affordable ab initio methods are TDHF (as explored here), TDDFT with hybrid
functionals, or the $GW$+Bethe-Salpeter equation approach without the TDA~\cite{Sander2015}.

\section*{Acknowledgments}

T.C.B.~thanks Philippe Guyot-Sionnest for early conversations related to this work.
All calculations were performed with the PySCF software package~\cite{Sun2018},
using resources provided by the University of Chicago Research Computing Center.
This work was supported by the Air Force Office of Scientific Research under
AFOSR Award No.~FA9550-18-1-0058 and by the National Science Foundation CAREER
program under Award No.~CHE-1848369.
The Flatiron Institute is a division of
the Simons Foundation.