Supplementary Information - Theory of photoluminescence spectral line shapes of semiconductor nanocrystals

The synthesis of the CdSe cores and the first fast-injection shell growth was carried out by following a previously published procedure by Carbone et al.^Carbone2007 For the CdS shell growths, to a 250mL round bottom flask was added 100nmol of the CdSe cores dissolved in hexane, 3mL of ODE, 3mL of oleylamine, and 3mL of oleic acid. The solution was degassed at r.t. for 2 hr. and then for 5 min at 100${}^{\circ}C$ to remove the hexane and water. The solution was then stirred under N2 and the temperature was set to rise to 310${}^{\circ}C$. Once the temperature reached 200${}^{\circ}C$, a solution of Cd(OA)2 in ODE (0.08M) and a separate solution of octanethiol dissolved in ODE (0.096M) were injected at a rate of 2.5mL/hr. Finally, the CdSe/CdS NCs were overgrown with a ZnS shell. For this, 0.1M Zn(OA)2 and 0.12 M octanethiol/ODE were injected at 3 ml/hr to the CdSe/CdS NCs for the appropriate amount of time to achieve desired shell thickness. After injection, the solution is further annealed at 310${}^{\circ}C$ for 15 min.

Samples for single molecule spectroscopy were prepared by diluting the stock-solution of nanocrystals in a PMMA solution (3% by weight in toluene), and spincoating (5000 rpm, 1 min) the solution on quartz slides (MTI, optical grade). The samples were cooled to liquid helium temperatures with a closed cycle liquid helium cryostat (Cryostation, Montana Instruments). To adjust the temperature, Montana Instrument’s Agile Temperature Sample Mount (ATSM), which allows fast and reliable adjustment of the temperature without too much spatial drift. Single particle spectroscopy was performed using a home-built confocal fluorescence microscope. Single nanocrystals were excited with a 488nm CW laser (Coherent Sapphire) for the recording of the temperature-dependent emission spectra and single nanocrystal temperature-dependent lifetimes. The single particle emission was filtered spatially (telescope (8 cm focal length with a 30 $\mu m$ pinhole) and spectrally (dichroic notch (488nm, Edmund Optics) and longpass filters (500nm, Thorlabs)). Single particle spectra were recorded by diffracting the mission in a monochromator (Acton 2050, Princeton Instruments) and detection with a cooled EMCCD camera (ProEM512B, Princeton Instruments). For the PL lifetime measurements, single nanocrystals were excited non-resonantly with a picosecond diode laser (PicoQuant, LDH-C-400 - 400 nm). The emission was collected in time-tagged mode (T3 mode, Picoquant Hydraharp). Photon-detection was performed with an APD (Excelitas) with a nominal timing resolution of FWHM 200ps.

The configurations of the core-shell CdSe/CdS NCs in theoretical calculations were constructed by first cutting the desired core size from bulk wurtzite CdSe geometry, and then adding several monolayers of shell (CdS). The core-shell geometries were relaxed using the conjugate gradient minimization with the classical Stillinger-Weber force field, which was parametrized for the bulk CdSe and CdS systems.^Zhou2013 Upon geometry relaxation, the outermost layer of atoms were removed and replaced with ligand pseudopotentials representing the passivation layer.^{Rabani1999, Jasrasaria2022a} The molecular formulae, core diameters, shell thicknesses and overall diameters are summarized below in Table 1 for all the CdSe/CdS core–shell NCs calculated in this work.

As stated in the main text, we use a model Hamiltonian to describe the excitonic states, phonon modes, and exciton-phonon interactions within a semiconductor nanocrystal, which is also weakly perturbed by an external electric field (See Eq.(1) from the main text). To parameterize this Hamiltonian within reasonable computational cost, we resort to the atomistic semiempirical pseudopotential method ^{Wang1995, Wang1996, Rabani1999} combined with the Bethe-Salpeter Equation (BSE).^{Rohlfing2000, Eshet2013} We also use a Stillinger-Weber force field ^Zhou2013 to describe the equilibrium geometry of the NC and to obtain the phonon modes. This pseudopotential-BSE approach and the set of parameters for CdSe and CdS have been shown to yield accurate descriptions of the exciton fine structure in CdSe/CdS core-shell nanocrystals systems.^{Rabani1999, Philbin2018, Jasrasaria2022a, Jasrasaria2021}

The atomistic semi-empirical pseudopotential approach, adapted from Ref. Jasrasaria2022a, was used to calculate the exciton fine structure of the CdSe/CdS nanocrystals. A non-interacting electron picture was used. The pseudopotential was assumed to have a local form in the reciprocal space, with a pre-factor expanded up to cubic order in the trace of the local strain tensor to account for deformation potential coupling:^{Jasrasaria2022a, Rabani1999, Wang1999}

$$\left\{-\frac{1}{2}\nabla^{2}+\sum_{\mu}\hat{v}_{\mu}\left(\mathbf{r}-\mathbf{R}_{\mu}\right)\right\}\psi_{i}=E_{i}\psi_{i}\,,$$

(1)

$$\hat{\tilde{v}}_{\mu}\left(\mathbf{k}\right)=\left[1+a_{4}Tr\epsilon_{\mu}+a_{5}\left(Tr\epsilon_{\mu}\right)^{3}\right]\cdot\frac{a_{0}\left(k^{2}-a_{1}\right)}{a_{2}e^{a_{3}k^{2}}-1}\,,$$

(2)

where $\hat{v}_{\mu}\left(\mathbf{r}\right)$ and $\hat{v}_{\mu}\left(\mathbf{k}\right)$ are the pseudopotential around atom $\mu$ in the real and reciprocal space, respectively. $\epsilon_{\mu}$ is the local strain tensor around atom $\mu$, and its trace is calculated as the tetrahedron volume formed by the nearest neighbors. Parameters $a_{0}$ through $a_{3}$ were fitted to the band structures of bulk CdSe and CdS obtained empirically or through high-accuracy electronic structure calculations, such as DFT+GW.^{Cohen1966, Hybertsen1986} The local strain-dependent prefactors ($a_{4}$ and $a_{5}$) were fitted to the bulk deformation potentials of the CBM and VBM for each semiconductor material.^{Jasrasaria2022a, Li2006} After the single-electron Hamiltonian is constructed, the filter diagonalization technique was used to obtain the quasi-particle states (electron or hole states) near the band edge of the nanocrystal.^{Neuhauser2000, Toledo2002} The pseudopotential parameters used in the work for Cd, Se and S atoms are given in Table 2.^{Jasrasaria2022a}

Excitonic states were expressed as a linear combination of electron-hole product wavefunctions, where the single-particle states are correlated. The coefficients ($c_{a,i}^{n}$ in the equation below) were obtained by solving the Bethe-Salpether equation within the static screening approximation.^Rohlfing2000

$$\psi^{n}\left(\mathbf{r}_{e},\mathbf{r}_{h}\right)=\sum_{a,i}c_{a,i}^{n}\phi_{a}\left(\mathbf{r}_{e}\right)\phi_{i}\left(\mathbf{r}_{h}\right)\,,$$

(3)

where $\psi^{n}\left(\mathbf{r}_{e},\mathbf{r}_{h}\right)$ is the exciton wavefunction. $\phi_{a}\left(\mathbf{r}_{e}\right)$ and $\phi_{i}\left(\mathbf{r}_{h}\right)$ are the single-particle electron and hole states obtained from the filter diagonalization technique, as discussed above. To converge the ground exciton state energy for the NC structures in this work, 60 electron and 60 hole states near the band edge were used.

The exciton-phonon coupling matrix element to first order was calculated non-adiabatically as the derivative of semi-empirical pseudopotential, following the methods and developments shown in Ref. Jasrasaria2021. The exciton-atomic-coordinate coupling $V_{n,m}^{\mu k}$ has contributions from the electron and the hole channels, where the coefficients from the Bethe-Salpether equation are factored in:

$$V_{n,m}^{\mu k}=\sum_{abi}c_{ai}^{n}c_{bi}^{m}v_{ab,\mu}^{\prime}-\sum_{aij}c_{ai}^{n}c_{aj}^{m}v_{ij,\mu}^{\prime}\,,$$

(4)

where $v_{rs,\mu}^{\prime}=\left\langle\phi_{r}\left|\left(\frac{\partial v_{\mu}\left(\left|\mathbf{r}-\mathbf{R}_{\mu}\right|\right)}{\partial R_{\mu k}}\right)_{\mathbf{R}_{0}}\right|\phi_{s}\right\rangle$, $c_{ai}^{n}$ and $c_{bj}^{m}$ are coefficients from the Bethe-Salpether equation, $\phi_{r}$ and $\phi_{s}$ are the single-particle states. The coupling matrix $V$ was then transformed from the atomic-coordinate $\mu k$ to the phonon mode coordinates $\alpha$ through the eigenvectors of the mass-weighted Hessian matrix, yielding the exciton-phonon coupling matrix.

$$V_{n,m}^{\alpha}=\sum_{\mu k}\frac{1}{\sqrt{m_{\mu}}}E_{\mu k,\alpha}V_{n,m}^{\mu k}\,,$$

(5)

where $E_{\mu k,\alpha}$ is the $\mu k$ element of the $\alpha$ eigenvector of the mass-weighted Hessian matrix D, which was constructed by taking the second-order derivative of the Stillinger-Weber potential.

Our photoluminescence linewidth model combines the information from exciton fine structure, phonon modes, and exciton-phonon couplings to calculate the emission spectrum within linear response and the dipole approximation.^{mukamel1995principles} The emission spectrum is calculated as the Fourier transform of the dipole auto-correlation function in the time domain. We also adopt the Condon approximation, which assumes that the dipole operator couples the ground and excited states. We first consider only the transition between the lowest excitonic state $\left|\psi_{1}\right\rangle$ and the ground state $\left|\psi_{g}\right\rangle$, and then extend our model to all emission transitions.

$$\hat{\bm{\mu}}=\bm{\mu}\left(\left|\psi_{g}\right\rangle\left\langle\psi_{1}\right|+\left|\psi_{1}\right\rangle\left\langle\psi_{g}\right|\right)\,,$$

(6)

The dipole auto-correlation function was then calculated by taking a partial trace over the electronic (system) degrees of freedom, resulting in an average over the nuclear/phonon (bath) degrees of freedom. The last term in Eq. (7) can be rewritten in the form of a time-ordered exponential, also known as the dephasing function: ^{Skinner1986, Egorov1997}

	$\displaystyle C_{\mu\mu}\left(t\right)$	$\displaystyle=\left\langle\hat{\bm{\mu}}_{H}\left(t\right)\hat{\bm{\mu}}_{H}\left(0\right)\right\rangle$
		$\displaystyle=\text{$\left|\bm{\mu}\right|$}^{2}\left\langle e^{i\left(E_{1}+H_{b}+\Delta\right)t/\hbar}e^{i\left(E_{g}+H_{b}\right)t/\hbar}\right\rangle_{B}$
		$\displaystyle=\text{$\left|\bm{\mu}\right|$}^{2}e^{i\left(E_{1}-E_{g}\right)t/\hbar}\cdot\left\langle F_{1}\left(t\right)\right\rangle\,,$		(7)

where the dephasing function between excitonic state $\left|\psi_{1}\right\rangle$ and the ground state is $\left\langle F_{1}\left(t\right)\right\rangle\equiv\left\langle\exp_{T}\left\{\frac{i}{\hbar}\int_{0}^{t}d\tau\Delta\left(\tau\right)\right\}\right\rangle$. And both $H_{b}=\sum_{\alpha}\hbar\omega_{\alpha}b_{\alpha}^{\dagger}b_{\alpha}$ and $\Delta=\sum_{\alpha}V_{1,1}^{\alpha}q_{\alpha}$ are operators in the nuclear (phonon) Hilbert space. $\left\langle\cdot\right\rangle=Tr_{b}\left[\rho_{b}\cdot\right]$. $\Delta\left(\tau\right)=e^{i\left(H_{b}\right)t/\hbar}\Delta e^{-i\left(H_{b}\right)t/\hbar}$ is the Heisenberg representation of the coupling operator. The subscript $T$ keeps the operators in decreasing time order.

A cumulant expansion on the dephasing function is then performed.^{Skinner1986, Egorov1997a, Egorov1996, Egorov1995} Since the operator $\Delta$ is linear in the phonon coordinates, only the second cumulant is non-zero, and is expressed as:

$$K_{2}\left(t\right)=\left(-\frac{i}{\hbar}\right)^{2}\int_{0}^{t}d\tau_{2}\int_{0}^{t-\tau_{2}}d\tau_{1}\sum_{\alpha,\beta}V_{1,1}^{\alpha}V_{1,1}^{\beta}\left\langle q_{\alpha}\left(\tau_{1}\right)q_{\beta}\left(0\right)\right\rangle\,,$$

(8)

Assuming all the phonons are harmonic oscillators and treating them quantum mechanically, the bath correlation function can be evaluated analytically.

$$\left\langle q_{\alpha}\left(t\right)q_{\beta}\left(0\right)\right\rangle=\frac{\hbar}{2\omega_{\alpha}}\left[\left(n\left(\omega_{\alpha}\right)+1\right)e^{-i\omega_{\alpha}t}+n\left(\omega_{\alpha}\right)e^{i\omega_{\alpha}t}\right]\delta_{\alpha\beta}\,,$$

(9)

where $n\left(\omega_{\alpha}\right)=\frac{1}{\exp\left(\beta\hbar\omega_{\alpha}\right)-1}$ is the Bose-Einstein distribution.

Thus, our model gives an analytical expression for the emission spectrum between the lowest exciton state to the ground state.^{Skinner1986, Egorov1997a, Egorov1996, Egorov1995}

$$I_{1\rightarrow g}\left(\omega\right)=\frac{\text{$\left|\bm{\mu}_{g1}\right|$}^{2}}{2\pi}\int_{-\infty}^{\infty}dte^{-i\left(\omega-\frac{E_{1}-E_{g}}{\hbar}\right)t}\exp\Bigg{\{}\frac{1}{2\hbar}\sum_{\alpha}\frac{\left(V_{1,1}^{\alpha}\right)^{2}}{\left(\omega_{\alpha}\right)^{3}}\left[C_{\alpha}^{\Re}(t)+iC_{\alpha}^{\Im}(t)\right]\Bigg{\}}\,,$$

(10)

where $C_{\alpha}^{\Re}(t)$ and $C_{\alpha}^{\Im}(t)$ are given in the main text. Our model is then easily extended to multiple electronic transitions, the effect of which only becomes noticeable when the temperature is high. In PL emission spectroscopy, the initial excitonic states are populated according to a Boltzmann thermal distribution, and the strength of the spectrum is proportional to the transition dipole moment. Therefore, the spectra of each exciton-ground-state transition are summed up with the corresponding weights, as expressed below:

$$I\left(\omega\right)=\frac{1}{Q_{n}}\sum_{n}e^{-\beta E_{n}}\left|\bm{\mu}_{gn}\right|^{2}\int_{-\infty}^{\infty}dte^{i\left(\omega-\omega_{ng}\right)t}\exp\Bigg{\{}\frac{1}{2\hbar}\sum_{\alpha}\frac{\left(V_{n,n}^{\alpha}\right)^{2}}{\left(\omega_{\alpha}\right)^{3}}\left[C_{\alpha}^{\Re}(t)+iC_{\alpha}^{\Im}(t)\right]\Bigg{\}}\,,$$

(11)

where a sum over all excitonic states $n$ is performed, $E_{n}$ is the energy of state $n$, and $\omega_{gn}=\frac{E_{n}-E_{g}}{\hbar}$. $\left|\bm{\mu}_{gn}\right|$ is the transition dipole moment between exciton state $n$ and the ground state, calculated by the following integral on the single-particle wavefunctions. $e$ is the electron charge.

$$\left|\bm{\mu}_{gn}\right|^{2}=\left\langle\psi_{g}\left|e\mathbf{r}\right|\psi_{n}\right\rangle=e\sum_{ai}c_{a,i}^{n}\int d\mathbf{r}\phi_{a}\left(\mathbf{r}_{e}\right)\mathbf{r}\phi_{i}\left(\mathbf{r}_{h}\right)\,,$$

(12)

Obtaining the FWHM from the analytical calculations of the spectra is tricky, since our model Hamiltonian lacks a dephasing time ^Skinner1986 and exhibits recurrences due to the finite number of phonon modes, resulting in infinitely narrow vibronic peaks (see main text Fig.1(b)). In order to reveal the intrinsic dynamics of the electronic transition and calculate experimentally-relevant FWHM, we convolved the calculated vibronic spectra with a time-domain Gaussian filter function with a width of $\sim 4$meV before comparing with experimental single-molecule PL measurements. This width is also consistent with the instrument resolution limit of the experimental setup, and the timescale of the decay of this filter function is faster than the recurrence time, $t_{r}\approx 2\pi/\omega_{\rm min}$. Furthermore, to match ensemble measurements, we convolved the calculated spectrum with another (broader) Gaussian function to account for the nearly-normal distribution of inhomogeneity of local environments in ensemble PL experiments.

We investigated whether anharmonicity of atomic motion at high temperatures, manifested through the finite lifetimes of phonons, affects the width of the PL spectrum. We concluded that the finite lifetimes of phonons only broadens the FWHM of PL spectra by around $10$meV at $300$K, and is thus insufficient to explain the mismatch in FWHM - temperature relationship between our linear exciton-phonon coupling model and the experimental single-NC measurements. We then empirically correct the dephasing rate by accounting for higher-order expansion terms in the exciton-phonon couplings, and analyze these additional decay channels (as discussed in the main text). Here, the details of anharmonicity are given.

The finite lifetime of the phonons were calculated following the methods outlined in Ref. Guzelturk2021. The relaxation lifetimes $1/T_{1}$ of each phonon mode in a $3.0$ nm core, $3$ML shell CdSe/CdS nanocrystal were calculated from the velocity auto-correlation function $C(t)=\langle v(0)v(t)\rangle$ using the following relation:^{Guzelturk2021}

$$\frac{1}{T_{1}}=\tilde{\xi}^{\prime}(\omega)=\Re\left[\frac{C(0)-i\omega\tilde{C}(\omega)}{\tilde{C}(\omega)}+i\omega\right]\,,$$

(13)

where $\tilde{\xi}^{\prime}(\omega)$ is the real part of the Laplace transform of the friction kernel from the generalized Langevin equation at the frequency of the mode. The velocity auto-correlation functions were calculated directly by simulating molecular dynamics (MD) trajectories which were sampled from a microcanonical ensemble at $300$K. The same Stillinger-Weber force field parameterized for bulk CdSe and CdS (as described in the main text and earlier sections of the SI) was used. The transformation from atomic coordinates to phonon mode coordinates was done through the eigenvectors of the mass-weighted Hessian matrix. The lifetimes of phonons that couple most strongly to the ground bright excitonic state have the most significant effects on the dephasing dynamics. As shown in Fig.2e of Ref. Guzelturk2021, the MD simulations reveal that the most strongly coupled acoustic phonons at $\sim 0.5$ THz have lifetimes of about $1$ ps, while the CdSe optical phonon modes at $\sim 8$ THz have lifetimes of around $100$ fs.

The finite phonon lifetime is incorporated into the bath correlation function by considering its damping effect on each phonon motion. The equations of motion of the creation and annihilation operators in these quantum mechanical, damped harmonic oscillators can be rewritten as:

$$\dot{\hat{b}}_{\alpha}=-i\omega_{\alpha}\hat{b}_{\alpha}-\frac{1}{T_{1}}\hat{b}_{\alpha}\,,\quad h.c.\,,$$

(14)

where $\hat{b}_{\alpha}$ and $\hat{b}^{\dagger}_{\alpha}$ are the annihilation and creation operators for phonon mode $\alpha$ with frequency $\omega_{\alpha}$ and lifetime $T_{1}$. For simplicity, we assume that all phonons have the same lifetime and thus, the bath correlation function for the damped harmonic oscillator is exponentially damped.

$$\langle q_{\alpha}\left(t\right)q_{\beta}\left(0\right)\rangle=e^{-\frac{t}{T_{1}}}\langle q_{\alpha}\left(t\right)q_{\beta}\left(0\right)\rangle_{HO}\,,$$

(15)

From here, we can plug the damped bath correlation function into Eq. (8) to obtain the analytical PL spectra:

$$I\left(\omega\right)=\frac{1}{Q_{n}}\sum_{n}e^{-\beta E_{n}}\left|\mu_{gn}\right|^{2}\int_{-\infty}^{\infty}dte^{i\left(\omega-\omega_{ng}\right)t}\exp\Bigg{\{}\frac{1}{2\hbar}\sum_{\alpha}\frac{\left(V_{n,n}^{\alpha}\right)^{2}}{\left(\omega_{\alpha}\right)^{3}}C_{\alpha}^{DHO}(t)\Bigg{\}}\,,$$

(16)

where for the each damped harmonic oscillators,

	$\displaystyle C_{\alpha}^{DHO}(t)=\omega_{\alpha}^{2}\cdot\Bigg{[}$	$\displaystyle\left(n(\omega_{\alpha})+1\right)\frac{e^{-t/T_{1}-i\omega_{\alpha}t}+\frac{t}{T_{1}}+i\omega_{\alpha}t-1}{(\frac{1}{T_{1}}+i\omega_{\alpha})^{2}}$
		$\displaystyle+n(\omega_{\alpha})\frac{e^{-t/T_{1}+i\omega_{\alpha}t}+\frac{t}{T_{1}}-i\omega_{\alpha}t-1}{(\frac{1}{T_{1}}-i\omega_{\alpha})^{2}}\Bigg{]}\,,$		(17)

Using the phonon lifetime $T_{1}\approx 100$fs at $300$K for $3$nm CdSe / $3$ML CdS core-shell nanocrystal, we obtained the anharmonic correction to the PL spectrum. At $300$K, incorporating anharmonicity increases the model PL FWHM from $28.2$meV (for harmonic phonons) to $37.8$meV. This correction of $\approx 10$meV in the single-NC FWHM only accounts for a very small fraction of the discrepancy between the experimental linewidth ($68.10$meV) and our model with linear exciton-phonon coupling($28.2$meV). Thus, we ruled out anharmonicity as the primary mechanism for the temperature-activated emission channel, and investigated high-order exciton-phonon coupling expansions, as discussed in the main text.