Symmetries in the collective excitations of an electron gas in core-shell nanowires

We study the collective excitations and inelastic light scattering cross-section of an electron gas confined in a GaAs/AlGaAs coaxial quantum well. These system can be engineered in a core-multi-shell nanowire and inherit the hexagonal symmetry of the underlying nanowire substrate. As a result, the electron gas forms both quasi 1D channels and quasi 2D channels at the quantum well bents and facets, respectively. Calculations are performed within the RPA and TDDFT approaches. We derive symmetry arguments which allow to enumerate and classify charge and spin excitations and determine whether excitations may survive to Landau damping. We also derive inelastic light scattering selection rules for different scattering geometries. Computational issues stemming from the need to use a symmetry compliant grid are also investigated systematically.


I. INTRODUCTION
Semiconductor nanowire (NW) lateral heterostructuressuch as coaxial heterojunctions or quantum wells-represent an important new class of nanomaterials with promising versatile properties for future applications in nanotechnology [1]. Many of their advantages derive from the high precision and reproducibility of "bottom-up" NW growing techniques [2], which allows for near-ideal atomically sharp interfaces to be engineered both in the axial and in the radial direction. Therefore, heterostructured NWs provide the possibility to tune quantum confinement properties by band structure engineering in the radial direction, while using the extended axis for facile transport and device integration, including with Si substrates thanks to the strain release in mismatched NW interfaces [3][4][5][6][7].
Coaxial heterostructures may also host a high-mobility electron gas (EG) [8]. The remote doping technique has been recently demonstrated in GaAs/Al x Ga 1−x As core-shell NWs, with the EG confined at the NW heterojunctions [9][10][11]. Hole gas in unintentionally doped structures can also be realized [12]. This opens up the realization of a variety of NW-based electronic devices [13][14][15]. On the other hand, this allows the investigation of fundamental properties of the EG in novel nanoscopic morphologies [10,11].
Traditional probes of the EG based on transport measurements, such as the Hall mobility, are difficult in NWs due to the difficulty of creating the required Ohmic contacts [16,17]. On the other hand, optical spectroscopies are nondestructive and * miguel.royovalls@unimore.it † andrea.bertoni@nano.cnr.it ‡ guido.goldoni@unimore.it contactless probes. The dynamics of photoexcited electronhole plasmas was studied by photoluminescence in single bare NWs [18] and core-shell NWs [19]. More recently, pump-THz probe spectroscopy allowed the determination of mobilities, lifetimes, and surface recombination rates of photoexcited carriers in different III-V NWs [20] and coreshell NWs [21]. Inelastic light scattering (ILS) has been used to extract carrier density and mobility data from the plasmon-phonon coupling modes of photoexcited NWs [22] and multilayered NWs [23]. Recently, we have used mean-field simulations combined with ILS experiments to demonstrate that remote doping induces high-mobility EG, and to assign ILS resonances to separate quasi-one-dimensional (q1D) and quasi-two-dimensional (q2D) channels in the sample [8].
In this paper, we study the collective excitations of an EG confined in a hexagonal coaxial QW (coQW), as engineered in a core-multishell GaAs/Al x Ga 1−x As NW. Channels showing q1D and q2D character are subsequently populated by varying the Fermi energy. We adapt the multisubband random-phase approximation (RPA) and time-dependent local-density approximation (TDLDA) formalisms, well established from studies in lower-dimensional systems, to perform a full 3D modeling of the NW electronic excitations, tracing the calculated CDEs and SDEs to the hexagonal symmetry of the system. Group theory is used to classify the complex set of collective excitations and it is shown that Landau damping into single-particle excitations takes place only for excitations of matching symmetry. Finally, we obtain ILS cross sections and predict ILS spectra under different scattering geometries, showing that the anisotropy of the system may be clearly exposed in ILS spectroscopy. The need for proper calculations of single-particle states and Coulomb matrix elements in a symmetry-compliant grid is emphasized.
The paper is organized as follows. In Sec. II, we sketch the theoretical model. We use the time-dependent density functional theory (TDDFT) formalism with one invariant direction to describe the coQW system (Sec. II A) and we formulate the nonresonant formalism employed to compute the ILS spectra (Sec. II B). Finally, in Sec. II C, the details of the computational procedure are summarized. In Sec. III, we illustrate CDEs and SDEs at various carrier densities. In Sec. III A, we describe the single-particle states used as the basis set to compute the excitations. In Sec. III B, we report our RPA and TDDFT results and, finally, in Sec. III C, we report ILS spectra computed for two relevant scattering geometries. To conclude, in Sec. IV, we summarize and discuss our results. In the Appendix, A we discus the calculation of Coulomb integrals on a triangular grid.

A. The linear-response TDDFT approach for q1D systems
In linear-response theory, the excitation energies of an interacting electron system can be obtained from the poles of the density-density response function [40]. In the Lehman representation, this so-called irreducible response function is written as Here, 0 and n are the many-electron ground-and excitedstate wave functions, respectively, n = E n − E 0 are the excitations energies,n(R) is the density operator expressed in the spatial coordinate R, and η is a positive small damping parameter. In the TDDFT formalism [41],˜ (R,R ,ω) is obtained from the response function of the noninteracting Kohn-Sham (KS) system 0 (R,R ,ω). The latter is formally obtained from Eq. (1) evaluating the matrix elements in the numerator assuming single Slater determinants built from the KS orbitals. For a NW that is translationally invariant along the z direction, the latter can be factorized as where φ n (r) is an envelope function over the NW in-plane directions r ≡ (x,y), and k z is the momentum along the inwire direction z. Correspondingly, the energy of the KS states ε n (k z ) is parabolic in k z . Likewise, the density operator is conveniently Fourier transformed along the invariant direction, which yields,n(r,q z ) = N l=1 δ(r − r l )e −iq z z l , with N being the total number of electrons. Altogether, the KS response function reads where g = 2 accounts for electron-spin degeneracy, f n (k z ) is the Fermi occupation function, and q z is the change in the in-wire momentum.
To obtain the response function of the interacting system, we expand it in terms of the KS orbitals as follows: The matrix elements˜ ij lm (q z ,ω) are then related to the elements of the KS response function (4) through the following Dyson equation: Here, v ij kn (q z ) and u XC ij kn are the direct Coulomb and exchangecorrelation matrix elements, respectively, which describe the dynamic interaction of two electrons, one of which gets scattered from state i to j and the other from k to n, with an exchange of momentum q z .
The direct Coulomb matrix elements read whereV (r − r ,q z ) is the Fourier transform of the Coulomb operator in the z direction. The exchange-correlation matrix elements read withf XC (r,r ) being the dynamic exchange-correlation kernel defined in the adiabatic local-density approximation as the derivative of the static exchange-correlation potential with respect to the ground-state density,f XC (r,r ) = dV XC [n0](r) dn 0 (r ) δ(r − r ). In particular, using the LDA functional conceived by Gunnarsson and Lundqvist [42], the exchange-correlation kernel reads [43] f XC (r,r ) = −1.704 a * 0 (r) 3 r s (r) 2 × 1 + 0.6213 r s (r) 11.4 + r s (r) Ry * (r) δ(r − r ).
Here, a * 0 (r) and Ry * (r) are the material-dependent effective Bohr radius and Rydberg energy, respectively, and r s (r) is the Wigner-Seitz radius.
The imaginary part of the irreducible response function (5) is proportional to the dynamic structure factor according to the fluctuation-dissipation theorem. As such, it gives a direct measure of the spectral strength of various elementary excitations of the electron system. In the long-wavelength limit (q z → 0), collective charge-density excitations (CDEs) or plasmons carry all the spectral weight and show up as narrow spectral peaks. However, as q z increases, electron-hole single-particle excitations (SPEs) start gaining spectral weight and show up as weak broad bands.
The TDDFT formalism also describes spin-density excitations (SDEs). The latter are only coupled through indirectexchange interactions. Consequently, they always appear redshifted with respect to their plasmonic counterparts. The SDEs are obtained from the poles of the so-called reducible response function, (R,R ,ω). The latter is calculated within a procedure similar to the above one, setting to zero the direct Coulomb integrals in Eq. (6). Furthermore, by deactivating the exchange-correlation integrals in Eq. (6), one recovers the CDEs within the random-phase approximation (RPA).

B. ILS cross section
In a nonresonant formalism, the ILS cross section is estimated from the imaginary part of the appropriate complete momentum-dependent response functions [30]. For CDEs, we Fourier transform the irreducible response function to obtaiñ Here, q is the in-plane component of the total momentum Q ≡ (q,q Z ) exchanged in the scattering process, i.e., Q = Q i − Q s , with Q i and Q s being the momenta of the incident and scattered photons, respectively. From the response functions, we may also calculate the density fluctuation induced by the electromagnetic field at a given energy and momentum, the so-called induced density distribution (IDD) from Kubo's correlation formula [44], δn(r,q,q z ,ω) = dr ˜ (r,r ,q z ,ω) e i qr .
The scattering cross section and IDD for SDEs are obtained analogously from the reducible response function, (R,R ,ω).

C. Computational methods
To evaluate the system response functions, we need the static properties of the system, that is, the single-particle energies and corresponding envelope functions of the q1D subbands, as well as the total electron density. At the meanfield level, these should be obtained from a self-consistent density functional theory (DFT) calculation [8]. However, since in this paper we aim at the dynamic properties of the EG, we clear up the calculation of the static properties by neglecting mean-field effects. In other words, the energy subbands are only determined by the unscreened band-offset modulation of the coQW. Since single-particle states would be different within different screening schemes, using a common unscreened confinement allows us to expose the difference in the dynamic properties within different formalisms, namely, TDLDA and RPA. On the other hand, mean-field effects would be dominated by band-offset confinement in narrow QWs as the present one. Mean-field calculations of static properties would be obviously required to study doped heterojunctions [10].
The envelope functions and subband energies are obtained within a single-band effective-mass approximation. Assuming in-wire spatial invariance along z and factorizing the envelope functions as in Eq. (2), the 2D envelope functions φ n (r) are given by where V (r) is the spatial confinement potential determined by the core-shell band offset. Equation (12) is numerically integrated in a hexagonal domain delimited by the NW surface using a symmetry-compliant triangular grid with 27.55 points/nm 2 , and assuming Dirichlet boundary conditions.
In the following section, we shall investigate the dynamical properties of the EG at different densities and subband occupations. This is performed by fixing the position of the Fermi energy to our interest, and then calculating the Fermi occupation of the subbands (we assume zero temperature throughout this paper) and the electron density, which is used in the dynamic exchange-correlation kernel, given by Eq. (9).
The computation of the Coulomb matrix elements entering the Dyson equation (6) is the numerically most intensive part of the procedure, as it requires the calculation of the 4D integrals (7) in a dense grid. To speed up the calculation, we have adapted the Fourier convolution theorem, which has been widely used to calculate Coulomb integrals in rectangular grids [45], to the case of a triangular grid. The method is described in the Appendix. We also take advantage of the system symmetries to reduce the number of Coulomb integrals to be computed as follows: (i) Since the electron wave functions φ n (r) are real, the following symmetries hold: (ii) As a consequence of the existence of a center of inversion in the system, we can classify the wave functions as gerade or ungerade, which allows us to discard a priori the computation of integrals with odd integrand.
(iii) Coulomb elements, in which the first two indexes (ij ) correspond to empty subbands, vanish. This is because such elements are multiplied in the Dyson equation (6) by the KS response function 0 ij (q z ,ω), which is zero for empty ij subbands [see Eq. (4)].
These very same considerations can be adopted to reduce the calculation of the exchange-correlation terms. However, in that case, the integrals are easier to compute due to the locality of the LDA potential, which reduces the dimensionality to 2D, and standard numerical integration methods can be used.
We do not consider effects from polarization charges accumulated at the NW surface and interfaces due to dielectric constant mismatch. In the present modeling, the only contribution from the dielectric confinement in the groundstate calculation would be the self-energy. This term has been shown to produce a small and almost constant shift of the energies, leaving unaffected the wave functions [46]. Dielectric confinement effects also enter the dynamic Coulomb and exchange-correlation matrix elements. This might affect the depolarization and excitonic shift of the collective modes, as it has been shown for NWs embedded in low dielectric constant matrices [37]. However, this quantitative effect does not break the hexagonal symmetry of the present system and, therefore, it will not alter the qualitative conclusions drawn here.
To compute the interacting response functions, we rewrite the Dyson equation (6) in tensor notation as or, equivalently,˜ In the above equation, I is the identity matrix, and ε = I − 0 (v + u XC ) is the dielectric tensor of the electron gas, whose inverse yields the required solution,˜ = 0 I ε −1 .
Recovering the subband notation leads tõ Therefore, in order to calculate a single element of the response function, we first build up the complete dielectric tensor of the EG and then invert it by means of efficient routines [47].

A. Single-particle states
The reference system for the following sections is a coremultishell NW (CSNW), as the one outlined in Fig. 1(a), with a hexagonal Al 0.3 Ga 0.7 As central core with 100 nm of diameter, a 20-nm-wide shell of GaAs, and a 10 nm Al 0.3 Ga 0.7 As capping layer. The 20 nm GaAs shell is a coQW for the conduction electrons. Note that this is somehow different from the usual samples, particularly in that the core is typically grown from GaAs. However, in doped samples, the core is usually depopulated due to band bending [8,10]. As we do not perform self-consistent calculations here, we use a Al 0.3 Ga 0.7 As core to exclude the formation of low-energy core states, which we are not interested in.
In the calculation, the GaAs/Al 0.3 Ga 0.7 As band offset is taken as 0.284 eV [48] and the origin of energies is placed at the GaAs conduction-band edge. The position-dependent effective mass m * (r) is 0.067 in GaAs and 0.092 in Al 0.3 Ga 0.7 As regions [49].
In Figs. 1(b) and 1(c), we show the normalized DOS and the in-plane envelope functions φ n (r), respectively, for the lowest  single-particle states. Two types of states can be identified, i.e., those preferentially localized in the corners and those which are delocalized over the coQW. The former tend to have lower energy [8]. The envelope functions show an increasing number of nodal planes normal to the coQW plane, which can be interpreted as the discretized momentum of a QW which is wrapped around the axis. States with nodal surfaces parallel to the coQW plane, analogous to excited subbands in planar QWs, lie at higher energies and are not shown here.
For the discussion of collective excitations and their ILS cross section, it will be useful to classify the single-particle states on the basis of their symmetry. Due to the overall hexagonal symmetry of the system, the eigenstates form a basis of the irreducible representations of the D 6h group. In Fig. 1(c), we label the symmetry irreducible representation (irrep) of each state, with the E type being doubly degenerate representations. The symmetry-induced degeneracies show up in the DOS as the higher peaks [see Fig. 1(b)].

B. Elementary excitation dispersion
We now look into the collective excitations of the EG in the coQW system at different density regimes, corresponding to the occupation of an increasing number of subbands. Here we are interested in classifying the lowest-energy excitations of the system, which correspond to CDE and SDE along the NW and around the NW section. Accordingly, to calculate the response functions (5), we use a basis set of singleparticle states restricted to the 12 states shown in Fig. 1(c) to simplify our analysis. This clearly neglects higher-energy excitations in the radial direction, since the QW higher states are not included. While we are not interested in quantitative predictions here, the excitations discussed in the following results are well converged with respect to the number of single-particle states except for the highest transitions.
The imaginary part of the response functions, given by Eq. (5), have poles at energies corresponding to electronic excitations in which the ground state is coupled with excited states through the density operator. We recall that within the mean-field formalism employed here, the wave functions of the effective noninteracting system are single Slater determinants. Hence, since the density operator is a one-body operator, the only accessible excited states are described by Slater determinants differing in a single excitation from the determinant describing the ground state. All excitations are coupled by the dynamic Coulomb and exchange-correlation matrix in the response function. In the following, we will adopt a widely used notation which labels the collective excitations with the single-particle transition of the final state with the largest contribution in the response function. We will also indicate the irrep of such final state as evaluated by multiplying the irreps of the two subbands participating in the single-particle transition.

One occupied subband
We first consider the case of the Fermi energy midway between subbands n = 1 and n = 2,3, which corresponds to a linear electron density of ∼0.007 × 10 7 cm −1 . In such a low-density regime, the linear-response formalism has been shown to fail to calculate the interaction of a q1D electron gas with an impurity [50]; however, its use to study elementary excitations of q1D electron systems is well grounded [51]. The dispersion of CDEs calculated within the RPA is shown in Fig. 2 as a function of q z /k F , with k F being the Fermi wave vector. One may recognize (i) one intrasubband CDE with vanishing energy at q z = 0; (ii) seven intersubband CDEs associated with transitions from n = 1 to each set of higher subbands, as indicated.
Since n = 1 has irrep A 1g , the symmetry of the excitation, which is the product of the irreps of the involved states, coincides with the irrep of the final state. Therefore, excitations to final states which are twofold degenerate are degenerate as well, and have, in general, larger spectral weight.
The dispersion of the intersubband excitations in q z is characteristic of 1D electron systems [32]. At small q z , the CDEs are blueshifted by the depolarization field from their analogous single-particle excitation, and approach the upper bound of single-particle continuum as q z increases, without entering it, and are not Landau damped. The intrasubband CDE shows a typical dispersion for q1D systems [36][37][38]52,53] proportional to q z √ |ln(q z )| in the long-wavelength limit.

Three occupied subbands
Next we consider a case with the Fermi energy midway between n = 2,3 and n = 4,5, with a linear electron density ∼0.039 × 10 7 cm −1 . The RPA result is plotted in Fig. 3(a) and in an enlarged scale in Fig. 3(c). The low-energy CDEs consist of (i) two intrasubband CDEs corresponding to the occupied subbands; (ii) the same intersubband CDEs from subband n = 1 as in the single occupied subband case with an additional CDE for the (1→2,3) transition, with a negative dispersion in the long-wavelength limit, lying between the two branches of the (1→2,3) single-particle continuum; (iii) multiple intersubband transitions from n = 2,3 to higher empty states.
Several additional issues may be discussed in this case. CDE degeneracies. Since n = 2,3 belong to the degenerate E 1u set, intersubband excitations to higher subbands lead to different number of CDEs, depending on the symmetry of the final state. Based on the product of the irreps of the involved subbands, we can identify the following cases: (i) If the final subband is nondegenerate, then the CDE is twofold degenerate. For instance, (2,3→7) yields the following degenerate irrep: E 1u ⊗ B 1u = E 2g .   The bounds of the single-particle excitation continua for excitations from n = 1 (n = 2,3) to higher subbands are delimited by blue solid (dashed) lines, which can be calculated analytically from energy and momentum conservation. Dark blue rectangles highlight selected CDEs which are Landau damped. Panels (c) and (d) show low-energy CDEs between the occupied subbands calculated within the RPA and TDLDA, respectively. Black dots: CDEs; gray (red) area: single-particle continuum for transitions from n = 1 (n = 2,3). The labels show the subbands involved in the excitation and the symmetry of the excited state.
(ii) If the final subband is twofold degenerate and of the same symmetry as n = 2,3, then two CDEs appear, one being doubly degenerate; for instance, (2,3→10,11) where A 2g is the basis of the antisymmetric representation of the permutation group. CDEs to such type of excited states are forbidden in the present spin-independent calculation since we always deal with singlet, and hence antisymmetric, spin wave functions which require symmetric orbital parts.
(iii) If the final subband is twofold degenerate and it is not of the same symmetry as n = 2,3, then three CDEs arise, one being doubly degenerate; for instance, for the excitation (2,3→8, 9) Depolarization shift. Intersubband CDEs couple through off-diagonal elements of the dielectric tensor [see, e.g., Eqs. (14) and (15)] only if they belong to the same symmetry. This leads to assorted depolarization shifts for different intersubband CDEs in Fig. 3(a). In general, however, these are larger than in Fig. 2 due to the increase of electron density.
Low-energy excitations. In Fig. 3(c), we observe two intersubband CDEs associated with the transition (1→2,3). The one with low energy shows an anomalous negative dispersion in the long-wavelength limit. Such type of plasmon is exclusive from q1D systems with more than one occupied subband and is free of Landau damping for large ranges of q z as it lies between the bounds of the (1→2,3) single-particle continuum [37,54,55]. The two intrasubband plasmons (1→1) and (2,3→2,3) are highlighted in Fig. 3(c). The (1→1)A 1g excitation is more dispersive than the corresponding excitation in Fig. 2. This is due to the higher-density regime and the coupling with the (2,3→2,3) intrasubband plasmon. The two CDEs couple because ( ] ⊕ E 2g has a symmetry component in common with the (1→1)A 1g plasmon. On the other hand, the (2,3→2,3)E 2g component appears as a slender acoustic plasmon [56] lying between the two intrasubband single-particle continua [see Fig. 3(c)]. This type of plasmon, free of Landau damping, is exclusive of q1D systems with more than one occupied subband [37,38,[54][55][56][57]. The antisymmetric component (2,3→2,3)A 2g is dark in the charge-density channel for the same aforementioned reasons.
Exchange and correlation effects. Figure 3(b) shows the CDEs calculated within the TDLDA. These are in oneto-one correspondence with the RPA CDEs, but redshifted by the dynamic exchange-correlation matrix elements. The exchange-correlation vertex correction may also overcome the direct Hartree term, bringing the plasmon below the corresponding single-particle excitations, as predicted by Das Sarma et al. [58] and later observed by Ernst et al. [59] in very dilute QWs. The symmetry-selective Landau damping phenomena observed in the RPA spectrum are also present within the TDLDA, marked by dark blue squares in Fig. 3(c), although they are not as well resolved as in Fig. 3(a).

Five occupied subbands
In Fig. 4, we show the CDEs for a case with the Fermi energy midway between n = 4,5 and n = 6 with a linear electron density ∼0.089 × 10 7 cm −1 . Calculations are shown for TDDFT only. In this higher-density regime, all of the intersubband plasmons appear blueshifted from their corresponding single-particle excitations in spite of the exchangecorrelation corrections. The same symmetry arguments used above to assign the excitation spectra apply also in this case and in higher subband occupations. We do not include labeling of all CDEs appearing in Fig. 4 for the sake of conciseness. Landau damping of selected CDEs is marked with dark blue rectangles for plasmons (2,3→7)E 2g , (2,3→6)E 2g , and (4,5→6)E 1u , which are clearly Landau damped in the single-particle continua with the same symmetry, (1→8,9)E 2g , (2,3→7)E 2g , and (4,5→7)E 1u , respectively.

Spin-density excitations
SDEs are computed from the imaginary part of the TDLDA reducible response function. In Fig. 5, we show SDEs for three occupied subbands. SDEs appear, in general, redshifted with respect to their corresponding single-particle excitations. This is due to the so-called excitonic shift caused by the exchangecorrelation matrix elements. SDEs are less dispersive than CDEs (compare Figs. 5 and 3). This is originated from the fact that the intersubband collective excitations become dispersive as they approach their corresponding single-particle continua. However, the bottom bound of the continua is less dispersive than the top one. Hence, the SDEs merge in the single-particle continua at larger q z . Once there, the slope of the dispersion is also lower. Additional differences with respect to CDEs are observed in the low-energy region: (i) Less peaks are visible which, nevertheless, show higher configuration mixing. For example, at ∼0.38 and ∼0.2 meV, there are two SDEs which have contributions from transitions (2,3→6) plus (1→4,5)E 2g and (2,3→4,5) plus (1→2,3)E 1u , respectively.
(ii) A single intrasubband SDE associated with transition (2,3→2,3)E 2g is observed. It disperses almost linearly in q z between the two intrasubband single-particle continua as the slender acoustic plasmon observed in Fig. 3.
SDEs associated with transitions (2,3→10,11) and (2,3→2,3) lead to excited states of symmetry A 1g ,E 2g , and E 2g , respectively. Such excited states resulting from transitions between degenerate states of the same E 1u symmetry are a basis of the symmetric representation of the permutation group. Notice that here these are not proper excited states for SDEs, which implicate a triplet, and hence symmetric, spin function requiring an antisymmetric orbital function. Therefore, for the two mentioned SDEs, one would expect to obtain the antisymmetric [A 2g ] state. This physically incorrect 155416-7 result is a shortcoming consequence of the widely used spin-independent formalism employed here [27,43,60].

C. Inelastic light-scattering spectra
We now focus on the ILS spectra in the nonresonant formalism discussed in Sec. II B, which has been widely used in spectroscopy of q2D and q1D electron systems (see, e.g., Refs. [24] and [25], and references therein). The authors have recently used this formalism to successfully assign ILS resonances of a high-mobility EG in modulation-doped core-multishell NWs [8].
The ILS cross section of CDEs is obtained from the imaginary part of the complete Fourier transform of the density-density response function, given by Eq. (10), which has a clear physical interpretation: the coefficients˜ ij lm (q z ,ω) give the spectral intensity of the CDEs or SDEs, while the matrix elements φ i (r)|e −i qr |φ j (r) φ m (r )|e i qr |φ l (r ) represent the coupling of light with electron charge density, which depend on the setup geometry. Therefore, the geometry of the experiment, defined by Q = Q i − Q s , sets specific selection rules. This is well known in q2D systems, where intrasubband or intersubband excitations can be selectively excited by choosing the exchanged momentum along or perpendicular to the QW plane, respectively. In coQW, the situation is clearly more complex. It is easy to realize that the photon momentum always has both in-plane and vertical components with respect to some of the facets, and one never recovers the ideal QW geometry.
Below, we will consider two backscattering configurations (see left insets in Fig. 6), both with the incoming and scattered photons perpendicular to the NW axis: (i) photons are perpendicular to the top/bottom NW facets; (ii) photons are parallel to the top/bottom NW facets. Accordingly, we set q z = 0 and we assume a typical excitation energy of 1.92 eV, which yields |Q i | = 6.9 × 10 5 cm −1 , and approximate |Q s | = |Q i |.
In the presence of the photons, the symmetry of the system is reduced from D 6h to D 2h . In configuration (i), the real and imaginary parts of the kernel, e ±i qr , have the irreps A 1g and B 2u of D 2h group, respectively: only those CDEs which have a A 1g or B 2u component in D 2h will be observed. Similarly, in configuration (ii), the kernel has the irreps A 1g and B 3u . This implies that some CDEs that are observed in one geometry are dark in the other, and vice versa.
As an example, in Figs. 6(b) and 6(c), we show the calculated ILS spectra for a NW with three occupied subbands computed from the TDLDA density-density response function, with a damping parameter η = 0.01 meV. In Fig. 6(a), we also show for comparison the imaginary part of the response function calculated at q z = 0, which has one peak for each CDE of the system. Peaks appear in both configurations if the corresponding excitations have irreps in D 6h , which reduces to A 1g or to B 2u ⊕ B 3u in D 2h . This is the case, for instance, of the peak labeled (2,3→10,11)A 1g , which has a symmetric representation also under the effect of the photons, as it can also be observed in the IDDs shown in the insets. Some of the resonances, however, can only be observed in one of the scattering geometries. Peaks (1→6)B 2u and (1→7)B 1u , for example, are only observed in configurations (ii) and (i) [Figs. 6(c) and 6(b)], respectively. Indeed, the irreps of these transitions in D 2h are B 3u and B 2u , respectively. The corresponding IDD in the insets clearly shows the same symmetries. Due to the same argument, CDEs associated with the same transition involving degenerate subbands can selectively be observed in the different geometries. For instance, the peak (2,3→8,9)B 1u , which has a very low spectral weight [see Fig. 6

IV. SUMMARY AND CONCLUSIONS
We have used TDLDA and RPA methodologies to study electron collective excitations and ILS cross section of GaAs/Al 0.3 Ga 0.7 As core-multishell NWs. This is a complex system from the electronic point of view, where q1D and q2D channels coexist. The large dimensions of the target system have required a 3D computational scheme. We have shown how to make calculations manageable for such a complex geometry by exploiting the symmetries of the system and, particularly, we have developed a fast and accurate approach to calculate Coulomb matrix elements in hexagonal grids. The latter are shown to be necessary to obtain convergence and correct degeneracies with a reasonable grid density (see details in the Appendix).
We have studied CDEs and SDEs of the EG at different density regimes, i.e., different number of occupied subbands of different localization (bents or facets) and degeneracy. We have identified a number of features ensuing from the discrete, here D 6h , symmetry of the system: (i) We have observed intersubband collective excitations, both CDEs and SDEs associated with transitions between twofold degenerate subbands split in different peaks in the spectra. The number of peaks is in agreement with the number of accessible excited states predicted by the symmetry group theory.
(ii) We have observed symmetry-selective Landau damping, namely, collective excitations are only Landau damped in single-particle continua associated with transitions of the same symmetry as the collective mode.
(iii) We have observed intrasubband slender acoustic plasmons and intersubband plasmons with negative dispersion exclusive of q1D systems with multiple subband occupation.
We have calculated the ILS spectra for two relevant scattering geometries in a backscattering configuration. As a result of the discrete symmetry of the system, the spectra are substantially anisotropic as the photon momentum is rotated around the NW axis. Some of the collective modes can be only observed in one of the geometries, being dark in the other. Selection rules are shown which explain the observed features. Although not included in the present study, ILS experiments may access excitations at larger energies where intersubband excitations between excited QW states, i.e., with a radial nodal plane, are involved. Simulations in this higherenergy range must also include coupling with GaAs phonon resonances [8].

ACKNOWLEDGMENTS
We acknowledge partial financial support from APOSTD/2013/052 Generalitat Valenciana Vali+d Grant, EU-FP7 Initial Training Network INDEX, and University of Modena and Reggio Emilia through grant "Nanoand emerging materials and systems for sustainable technologies."

APPENDIX: COULOMB MATRIX ELEMENTS CALCULATION IN HEXAGONAL GRIDS
The Coulomb matrix elements that we have to calculate are These, by taking g kn (r ) = φ * k (r ) φ n (r ) and rearranging the integrals, can be written as In the above equation, h(r) = d r V C (r − r ,q z ) g kn (r ) is the convolution ofV C (r − r ,q z ) and g kn (r ). Therefore, according to the Fourier convolution theorem, the Fourier transform of h(r) can be obtained ash(q) =Ṽ C (q,q z )g kn (q), where the tilde means a Fourier transformed function in momentum space.
Since h(r) can be equivalently obtained by performing the inverse Fourier transform ofh(q), i.e., h(r) = F −1 [h(q)], it is possible to calculate the Coulomb matrix elements as With this approach, the dimensionality of the Coulomb integrals is reduced from 4D to 2D at the expense of performing three discrete Fourier transforms (DFTs). The outcome in terms of computation time is highly favorable thanks to the existing fast Fourier transform (FFT) algorithms. Besides, V C (q,q z ) needs to be only calculated once and it has a well-known analytical expression, where ε ∞ is the high-frequency dielectric constant and q D is the momentum derived from the Debye length λ D as q D = 1 λ D . For the calculations shown in the present study, we have employed the GaAs dielectric constant, ε ∞ = 10.86, and a Debye length λ D = 1 μm. Available libraries with implemented FFT algorithms exclusively work with numerical sampling on rectangular grids. We have used this type of algorithm to calculate our Coulomb matrix elements by extrapolating our functions onto a rectangular grid beforehand. However, such a procedure leads to qualitative errors in the Coulomb integrals. Such errors, which do not vanish as the grid is made denser, are due to the inefficacy of rectangular sampling to capture the implicit hexagonal symmetry of the present system. For instance, in the first column of Table I, we show two diagonal Coulomb matrix elements between the ground and the first twofold degenerated excited states calculated with a FFT algorithm in a rectangular grid. They both should have the same value since the degeneracy is induced by the D 6h symmetry; however, we observe a discrepancy of ∼8 meV nm. There exist alternative algorithms [61] in which the input data is sampled in a hexagonal grid but the output of the Fourier transform is sampled on rectangular grids. Their use does not overcome the discrepancy in the Coulomb integrals calculation, though, as can be seen in the second column of Table I. A proper calculation of the Coulomb integrals requires the use of an algorithm that performs the DFT completely in a hexagonal grid. Here we have adapted to our interest the formalism reported by Mersereau [62] initially devoted to signal processing.
The formulas appearing in Ref. [62] to perform the hexagonal discrete Fourier transform (HDFT) of a 2D function g(n 1 ,n 2 ) defined in a hexagonal domain are reported here for completeness: g(k 1 ,k 2 ) = Here, n 1,2 (k 1,2 ) are integer coordinates denoting points of the grid in the position (momentum) space, and the sum is restricted to those points inside the domain limiting the grid R H (N,N). The choice of the latter is not trivial since a DFT assumes that the input function is periodically replicated in the space. Moreover, the replication pattern will determine the sampling of the output function of the DFT. Thus, in order to obtaing(k 1 ,k 2 ) sampled on a hexagonal grid, one has to assume that g(n 1 ,n 2 ) is periodically replicated with a hexagonal pattern. This implies that the following relations have to be accomplished: g(n 1 ,n 2 ) = g(n 1 − 3N,n 2 ) = g(n 1 − 2N,n 2 − N ) = g(n 1 − N,n 2 − 2N ), with N being the replication period. A regular hexagonal domain as the one shown in Fig. 7(a), delimited by the solid line, which is equivalent to the one in which our functions are sampled, is not a proper domain to perform the HDFT. The reason is that its periodic replication according to Eqs. (A7) leads to empty gaps and overlapping between the replicas, which is known to produce aliasing. Instead, we use as R H (N,N) the deformed hexagonal region depicted by the gray area in Fig. 7(a), which can be exactly hexagonally replicated, as illustrated in Fig. 7(b). Proceeding in this way, we obtain an output functiong(k 1 ,k 2 ) hexagonally sampled on a equivalent domain as the one illustrated by the gray area in Fig. 7(c). The dimensions of the domain in the momentum space, shown as W 1 and W 2 in Fig. 7(c), are fixed by the sampling theorem as W 1 = 4π 3T 1 and W 2 = π T 2 , where T 1 and T 2 are the sampling intervals in the position space [see Fig. 7(a)].
As shown in the third column of Table I, the use of the HDFT solves the symmetry-induced numerical discrepancies in the calculation of the Coulomb matrix elements. Its formal computation, though, is computationally more demanding than in a rectangular sampling formalism, mainly due to the impossibility to separate the Fourier kernel in 1D-like DFTs. However, there exists a fast Fourier implementation of the algorithm, fairly described in Ref. [62], that will allow one to increase the efficiency of the computation.