While I have written many posts on this blog about various topics in physics or math unrelated to my graduate work as well as posts promoting papers from my graduate work, it is rare that I've written direct technical posts about my graduate work. It is even more unusual that I should be doing so 2 years after leaving physics as a career. However, I felt compelled to do so after meeting again with my PhD advisor (a day before the Princeton University 2020 Commencement, which was held in person after a delay of 2 years due to this pandemic), as we had a conversation about the problem of infinite local density of states (LDOS) in a lossy medium.
Essentially, the idea is the following. Working in the frequency domain, the electric field produced by a polarization density in any EM environment is \( E_{i}(\omega, \vec{x}) = \int G_{ij}(\omega, \vec{x}, \vec{x}')P_{j}(\omega, \vec{x}')~\mathrm{d}^{3} x' \) which can be written in bra-ket notation (dispensing with the explicit dependence on frequency) as \( |\vec{E}\rangle = \hat{G}|\vec{P}\rangle \). The LDOS is proportional to the power radiated by a point dipole and can be written as \( \mathrm{LDOS}(\omega, \vec{x}) \propto \sum_{i} \mathrm{Im}(G_{ii}(\omega, \vec{x}, \vec{x})) \). This power should be finite as long as the power put into the dipole to keep it oscillating forever at a given frequency \( \omega \) is finite. However, there appears to be a paradox in that if the position \( \vec{x} \) corresponds to a point embedded in a local lossy medium, the LDOS diverges there.
I wondered if an intuitive explanation could be that loss should properly imply the existence of energy leaving the system by traveling out of its boundaries, so the idea of a medium that is local everywhere (in the sense that the susceptibility operator takes the form \( \chi_{ij}(\omega, \vec{x}, \vec{x}') = \chi_{ij}(\omega, \vec{x})\delta^{3} (\vec{x} - \vec{x}') \) at all positions) and is lossy at every point in its domain may not be well-posed as energy is somehow disappearing "into" the system instead of leaving it. Then, I wondered if the problem may actually be with locality and whether a nonlocal description of the susceptibility could help. This is where my graduate work could come in. Follow the jump to see a very technical sketch of how this might work (as I won't work out all of the details myself).
Scalar response
Scalar response implies the simplification that the vacuum Green's function and susceptibility are scalars. For the purpose of this derivation, the vacuum Green's function can be approximated as \( G^{\mathrm{vac}}(\omega, \vec{x}, \vec{x}') = \left(\frac{\omega}{c}\right)^{2} \frac{\exp(\mathrm{i}\omega |\vec{x} - \vec{x}'|/c)}{4\pi|\vec{x} - \vec{x}'|} \); in operator form, it will be written as \( \hat{G}^{\mathrm{vac}} \). Meanwhile, the susceptibility will be taken to schematically represent a single electron with a given density, so that can be written as \( \chi(\omega, \vec{x}, \vec{x}') = \alpha(\omega) f(\omega, \vec{x})f^{\star}(\omega, \vec{x}') \) for some polarizability \( \alpha \) and some density function \( f \). In bra-ket notation, this is written as \( \hat{\chi} = \alpha|f\rangle\langle f| \), making clear that \( \hat{\chi} \) has a rank of 1. Additionally, the identity operator \( \hat{1} \) has the representation \( \delta^{3}(\vec{x} - \vec{x}') \).
In the scalar case, the LDOS can be defined using bra-ket notation as \[ \mathrm{LDOS}(\omega, \vec{x}) \propto \mathrm{Im}(\langle \delta_{\vec{x}}, \hat{G} \delta_{\vec{x}} \rangle) \] where \( |\delta_{\vec{x}}\rangle \) has the representation \( \delta^{3}(\vec{x} - \vec{x}') \) using \( \vec{x}' \) as the dummy variable. Furthermore, the definition of the total Green's function is \[ \hat{G} = \hat{G}^{\mathrm{vac}} + \hat{G}^{\mathrm{vac}} \hat{T} \hat{G}^{\mathrm{vac}} \] after defining the T-operator \[ \hat{T} = \hat{\chi} (\hat{1} - \hat{G}^{\mathrm{vac}} \hat{\chi})^{-1} \] in terms of the desired operators. Plugging in the definition of \( \hat{\chi} \) and using the fact that it has a rank of 1 gives \( \hat{T} = \frac{\alpha}{1 - \alpha \langle f, \hat{G}^{\mathrm{vac}} f\rangle} |f\rangle\langle f| \). This therefore gives \( \mathrm{LDOS}(\omega, \vec{x}) \propto \mathrm{Im}(\langle \delta_{\vec{x}}, \hat{G}^{\mathrm{vac}} \delta_{\vec{x}} \rangle) + \mathrm{Im}\left(\frac{\alpha \langle \delta_{\vec{x}}, \hat{G}^{\mathrm{vac}} f\rangle\langle f, \hat{G}^{\mathrm{vac}} \delta_{\vec{x}} \rangle}{1 - \alpha \langle f, \hat{G}^{\mathrm{vac}} f\rangle} \right) \). The first term is evaluated as \( \mathrm{Im}(\langle \delta_{\vec{x}}, \hat{G}^{\mathrm{vac}} \delta_{\vec{x}} \rangle) = \frac{1}{4\pi} \left(\frac{\omega}{c}\right)^{3} \) regardless of \( \vec{x} \). The second term will not only depend on \( \vec{x} \) but will depend on the specific choices of the functions \( \alpha(\omega) \) and \( f(\omega, \vec{x}) \).
Consider a generic complex isotropic \( \alpha(\omega) \), representing loss, but a specific isotropic Gaussian density \( f(\omega, \vec{x}) = \frac{1}{(2\pi\sigma^{2})^{3/2}} \exp\left(-\frac{|\vec{x}|^{2}}{2\sigma^{2}}\right) \) where the width parameter is defined as \( \sigma(\omega) = (|\alpha(\omega)|/3)^{1/3} / (2\sqrt{\pi}) \) in terms of the polarizability. In such a case, the relevant matrix elements can be calculated as \( \langle f, \hat{G}^{\mathrm{vac}} f\rangle = \frac{1}{8\pi\sigma} \left(\mathrm{i}q - q~\mathrm{erfi}(q/2) + \frac{4}{\sqrt{\pi}}\exp(q^{2} / 4)\right) \) after defining \( q \equiv \frac{2\omega\sigma}{c} \) and \( \langle f, \hat{G}^{\mathrm{vac}} \delta_{\vec{x}}\rangle = \frac{\exp(-q^{2} / 8)}{8\pi\rho \times 2\sigma} (\exp(\mathrm{i}\rho q) \mathrm{erfc}(-\mathrm{i}q/(2\sqrt{2}) - \sqrt{2}\rho) - \) \( \exp(-\mathrm{i}\rho q) \mathrm{erfc}(-\mathrm{i}q/(2\sqrt{2}) + \sqrt{2}\rho)) \) using the same definition of \( q \) as before and defining \( \rho = \frac{|\vec{x}|}{2\sigma} \). The analytical expression for the LDOS becomes quite complicated & cumbersome, but what is clear is that neither the real nor imaginary parts of the terms inside of the definition of the LDOS diverge. Therefore, I think it is fair to claim that under this model of lossy nonlocal response with a single dipolar nanoparticle, the LDOS doesn't diverge for any \( \vec{x} \).
Vector response
The case of vector response is harder to write analytically just because the susceptibility has a rank of 3 instead of 1, so instead of taking the reciprocal of a scalar, one must take the inverse of a 3-by-3 matrix. That said, the spatial dependence of the susceptibility, if isotropic, is still assumed to come only through the single function \( f(\omega, \vec{x}) \) (and the rank of 3 just comes from the 3 independent Cartesian directions). Furthermore, the problem with the LDOS diverging in a uniform lossy medium arises because the imaginary part of the susceptibility affects the term in the Green's function that goes as \( |\vec{x} - \vec{x}'|^{-3} \) (which is a non-integrable singularity) as \( \vec{x}' \to \vec{x} \), but multiplying that term (and any singular terms that are more integrable) by \( f^{\star}(\omega, \vec{x})f(\omega, \vec{x}')~\mathrm{d}^{3} x~\mathrm{d}^{3} x' \) and integrating yields a finite result. Therefore, I see no reason to believe that the LDOS would diverge in that case either. Thus, I am more convinced that nonlocality is needed for previously problematic results to become mathematically sensible in lossy systems.
A note about rank for infinite-dimensional spaces
The notion of rank for operators in infinite-dimensional vector spaces is complicated. For a finite-dimensional vector space of size \( N \), an operator can have a whole number rank of at most \( N \). However, for an operator on an infinite-dimensional vector space in which the labels for the singular values are countable, meaning \( \hat{A} = \sum_{n = 0}^{\infty} \sigma_{n} |u^{(n)}\rangle\langle v^{(n)}| \), it is possible for the operator to have infinite but not full rank, for example if \( \sigma_{n} \) is nonzero only for odd labels \( n \). Weirder things can happen when labels for the singular values are uncountable; for simplicity, the label \( \theta \) will be assumed to lie in the interval \( [0, 1] \), so \( \hat{A} = \int_{0}^{1} \sigma(\theta) |u(\theta)\rangle\langle v(\theta)|~\mathrm{d}\theta \). If \( \sigma(\theta) \) vanishes only for \( \theta > 1/2 \), then \( \hat{A} \) is not of full rank, but its rank is "proportional" to that of the identity operator represented as \( \delta(\theta - \theta') \) on the full interval \( [0, 1] \). If \( \sigma(\theta) \) vanishes only for \( \theta \) that lies in the Cantor set fractal defined on the interval \( [0, 1] \), then \( \hat{A} \) has "almost full" rank. If \( \sigma(\theta) \) vanishes only for \( \theta \) that lies in the complement of the Cantor set fractal defined on the interval \( [0, 1] \), then \( \hat{A} \) has "almost zero" rank but is not the same as the zero operator that maps every vector to the null vector. I don't have formal training in how operator ranks should be rigorously defined for infinite-dimensional vector spaces, but these are some thoughts that came to my mind about the subject.