Non-Hermitian approach for quantum plasmonics

Plasmonics is the study of light interactions with generally metallic structures that can be resonantly excited to yield intense and localized electromagnetic responses of interest both for fundamental and practical reasons (see, e.g., Refs. 1 and 2). Of interest is the coupling of plasmonic systems with molecules, nanostructures, or materials that exhibit quantum mechanical responses that could conceivably be enhanced or coupled into the plasmonic structure in some fashion to achieve interesting outcomes, i.e., quantum plasmonics.^3,4

We have previously modeled quantum dot/plasmon interactions^5–7 using a cavity quantum electrodynamics approach based on the Lindblad master equation⁸ for the quantum mechanical density matrix. This is a reasonably rigorous approach that can incorporate important environmental effects such as dissipation and dephasing but can be computationally intensive. Furthermore, it is desirable to develop simpler models for analysis purposes that can still convey some of the correct dynamics.

Non-Hermitian quantum mechanics generally corresponds to the study of time-dependent Schrödinger equations (TDSEs) or time-independent Schrödinger equations with Hamiltonian operators that are not Hermitian.^9,10 The non-Hermitian terms in the Hamiltonian are designed to describe processes such as interaction with an environment that are not explicitly included as degrees of freedom in the Hamiltonian. While much recent work has focused on parity-time symmetric non-Hermitian Hamiltonians that have real eigenvalues and can exhibit interesting properties such as exceptional points,^10,11 our focus here is on non-Hermitian Hamiltonians with complex eigenvalues that can mimic to some degree dissipation and dephasing relevant to quantum dot or exciton interactions with plasmonic systems such as metal nanoparticles or arrays of metal nanoparticles. Recently, several papers have, in fact, studied the aspects of quantum plasmonics with such non-Hermitian Hamiltonians.^12–14 The non-Hermitian Hamiltonian results in systems of equations that qualitatively (and sometimes quantitatively) capture the true dynamics while being significantly more efficient to solve, compared to the full density matrix approach. Using such non-Hermitian Hamiltonians can allow the study of systems far larger than otherwise possible, such as the investigation of entanglement dynamics of large numbers of quantum dots coupled to a plasmonic system. The principle contributions of our work are (i) to assess the validity of this approach via a comparison with the full Lindblad dynamics and (ii) to demonstrate how interesting results can be generated for very challenging problems involving large numbers of quantum dots.

Section II outlines both the Lindblad (Sec. II A) and non-Hermitian (Sec. II B) approaches, Sec. III outlines the results we have obtained for one and two quantum dot systems coupled to a plasmon, as well as an example exploring the role of coupling disorder involving 50 quantum dots. Section IV presents concluding remarks.

The Lindblad master equation⁸ for the time-dependent density matrix, ρ(t), is

where H(t) is the system Hamiltonian (which may or may not include external driving in time) and L{} is the Lindblad superoperator that acts on the density matrix as follows:

where the sum is over the number of relevant dissipation or dephasing terms and C_k are collapse operators to be specified below.

The Hamiltonians we study are those N quantum dots coupled to a plasmonic mode,^5–7 

where σ_j is a 2-state quantum dot lowering operator for quantum dot j, b is a bosonic annihilation operator for the plasmonic mode, and the g_j are quantum dot/plasmon coupling rates. (Specific values for these and other parameters entering into the model will be given in Sec. III A.) If an external field, E(t), is being considered, then the dipole operator is taken to be

For simplicity, we take the quantum dots to have the same excitation frequency and dipole moment parameters; in general, this is not necessary.

The collapse operators {C_k} correspond to spontaneous emission and pure dephasing for the quantum dots and plasmon damping (e.g., Ref. 15) and are, respectively,

With specification of a finite basis corresponding to two states per quantum dot and N_pl plasmon states, the density matrix ρ is either 2N_pl × 2N_pl for the case of one quantum dot or 4N_pl × 4N_pl for the case of two quantum dots, and its elements can be written in operator form via (for the N = 2 case)

where each q_j = 0 or 1 for ground or excited quantum dot states and s = 0, 1, 2, …, (N_pl − 1).

We implement the density matrix approach described here using the convenient and freely available Quantum Toolbox in Python (QuTiP).^16,17 (For the low-intensity results here, N_pl = 5 suffices; for the high-intensity results, N_pl = 15.)

Instead of solving the Lindblad master equation [Eq. (1)], the non-Hermitian model involves solving a time-dependent Schrödinger equation (TDSE),

where H_c(t) is a complex Hamiltonian matrix or operator given by

where H(t) is the system Hamiltonian [Eq. (2)] and the C_k are the collapse operators appearing in the Lindblad superoperator [Eq. (3)], which take on the more explicit forms for our problems of interest in Eq. (5).

We have written the complex Hamiltonian in the form of Eq. (8) in order to draw its connection with the Lindblad master equation for the density matrix, and in this form, it can also be recognized as the effective Hamiltonian that enters into stochastic Schrödinger equation approaches as the first stage before any probabilistic collapses.^18,19 However, one can easily re-express Eq. (8) as

where $Γ=2γ2*+γ1$⁠. This connection between a non-Hermitian Hamiltonian and the stochastic Schrödinger equation has been pointed out before (see, e.g., Refs. 12–14).

If the wave packet is written as

then Eq. (6) yields first-order differential equations for the time derivatives $dαs,q1,q2(t)/dt$ and the contributions from the imaginary terms in Eqs. (8) or (9), i.e., the loss terms in these equations, are easily seen to be $−Γ2q1+q2+sγplαs,q1,q2(t)$⁠. There is, therefore, no loss for number states q₁ = 0, q₂ = 0, or s = 0; in addition, the effective loss rate for each s > 0 plasmon state is $sγpl2$⁠, i.e., it increases with s. The non-Hermitian model results in a wave packet of size 2N_pl for the case of one quantum dot. This is compared to the density matrix size of 2N_pl × 2N_pl for the full Lindblad master equation. Thus, the non-Hermitian model results in significant computational savings while still including the effects of dissipation and dephasing.

We first consider the case of one quantum dot or exciton interacting with a plasmonic system. For an earlier, pioneering quantum mechanical treatment of such coupled exciton-plasmon systems, see Ref. 20. Figure 1(a) is a schematic of the type of system being considered. We employ the parameters of Shah et al.,⁵ which we list in Table I, along with other parameters used in our calculations. The parameters in this system are chosen to be consistent with the optical properties of two ellipsoidal gold nanoparticles of dimensions in the 20–30 nm range with a small (4 nm diameter) CdSe quantum dot placed within the 6 nm gap between the particles, all embedded within a medium of refractive index n_med = 1.5. Illumination is along the common, long axis and under such conditions, it is possible to have a dipole or exciton-induced transparency effect^5,21 wherein a sharp dip is superimposed on a broader, plasmonic line shape. This phenomenon arises from Fano interferences.^5,20,21 In this sense, it is similar to a different and better known phenomenon, electromagnetically induced transparency (EIT).²² However, EIT is a nonlinear optical effect requiring both an intense pump and a probe laser.²² 

We compute optical absorption spectra by considering the system in its ground state initially and exposing it to a short Gaussian pulse with sufficient energy content to describe the spectral region of interest,

The optical spectrum is then computed from⁵ (SI units)

where the polarizability is given by

with the time integrals extending over times consistent with the system’s response, the expectation value of the dipole moment [see Eq. (4)], ⟨μ(t)⟩, rising and falling to zero after the pulse.

Figure 1(b) displays the result of the non-Hermitian model (solid curve) and the full density matrix result based on solving the Lindblad density matrix (symbols). The agreement is quite good, with there being some small discrepancies in the narrowest region of the dipole-induced transparency. The dip is quite sensitive to quantum dot dephasing, as is evidenced in Fig. 2 where we show how varying the dephasing rate from $ℏγ2*=0$ to 0.005 08 eV, keeping all other parameters as in Table I, can dramatically reduce the transparency. Figure 2(a) is the non-Hermitian model result, and Fig. 2(b) is the density matrix result, again showing that the non-Hermitian model agrees well with the density matrix one.

Figure 3 displays optical spectra results for two quantum dots interacting with the plasmonic mode. Once again, very good agreement is seen between the non-Hermitian and density matrix approaches. In this case, the dip is a little more pronounced than in the single quantum dot case of Fig. 1. In fact, these results are consistent with a Tavis–Cummings picture²³ wherein the effect of N two-level systems interacting with a cavity (the plasmon) can be described by a single system interacting with the cavity with an effective system–cavity coupling factor of $Ng$⁠. Indeed, if we carry out a calculation with one quantum dot interacting with the plasmon but employ dot–plasmon coupling $2g$⁠, the corresponding spectrum is virtually superimposable on the full two-dot results.

The optical spectra discussed above were inferred from Fourier transformation of the results from excitation with short pulses. The non-Hermitian model, however, can fail when the system is pumped with continuous wave (cw) light, e.g., when the exponential term in Eq. (11) is set to unity. With both plasmon dissipation and dephasing operative, the system cannot achieve steady state populations owing to the loss terms leading to eventually complete loss of wave packet amplitude. Thus, the model cannot describe the saturation effect of Ref. 5 that results when cw light is applied with ever increasing magnitudes of E_L. In this case, the Lindblad density matrix formalism leads to diminishing of the transparency effect, i.e., the transparency dips become ever less deep to eventually being absent for sufficiently large magnitude E_L.⁵ Figure 4 illustrates the failure of the non-Hermitian model for cw excitation and E_L = 1.4 × 10⁻⁶ a.u. (corresponding to an intensity on the order of 0.1 MW/cm², i.e., 100 times larger than that used in the calculations to generate Figs. 1–3). In this case, we are continually applying the driving on the frequency or energy of the transparency dip, 2.042 eV. Whereas the low-intensity optical cross section from Fig. 1 is ≈3 × 10⁻¹¹ cm², the resulting cross section in this case is ≈6 × 10⁻¹¹ cm², i.e., the dip has nearly been eliminated. In Fig. 4, we see that the density matrix calculations (symbols) yield a steady-state quantum dot number state expectation value, ⟨σ^†σ⟩, which is also the probability of excitation of the quantum dot, of slightly over 0.4, and the plasmon number state average, ⟨b^†b⟩, attains a value just under 0.1. In contrast, the non-Hermitian model (solid curves) is completely decayed away. The wave packet norm is shown as a green curve in Fig. 4, and the failure of the model can be correlated with it becoming significantly less than unity. The inset in Fig. 4 focuses on the shorter time behavior and shows that at least in this limit there is some agreement between the two approaches. We have experimented with a variety of “fixes,” including renormalization after each time step and the introduction of gain terms. However, we have not yet found a suitable fix that would maintain the simplicity of a single wave packet propagation.

Rather than drive the system with a pulsed or cw laser, as in Sec. III A, one can instead imagine the system to be already excited in some particular way and follow the subsequent time evolution in the absence of driving. For the case of two quantum dots coupled to a plasmon, we have found that if one dot is excited, that energy transfer, mediated by the plasmon, can occur to the other quantum dot and that during the course of this process, the two quantum dots can exhibit a reasonable degree of entanglement.^6,7,24 The non-Hermitian model can reproduce this type of behavior. This behavior does occur for the model parameters we have been using so far corresponding to those of Ref. 5. However, the dynamics is somewhat uninteresting, exhibiting no transitory coherences. Somewhat more interesting behavior occurs when one uses the parameters of Ref. 24, which are designed to be consistent with a more complex plasmonic structure that exhibits gap plasmons. Employing these parameters (also listed in Table I), and now with no driving (E_L = 0) but setting the initial condition to be one of the quantum dots is in its excited state, we obtain the results in Figs. 5(a) and 5(b).

The upper panel of Fig. 5(a) shows that excited quantum dot probability comes down from unity as the second quantum dot excitation probability rises and there is also a subsequent, secondary coherence. However, while qualitatively correct, the non-Hermitian model (solid curves) does show discrepancies with the density matrix results (symbols). The lower panel of Fig. 5(b) shows the associated concurrence²⁵ of the two quantum dots, calculated by tracing out the plasmon quantum numbers to obtain a reduced two quantum dot density matrix and then performing the required computations.^6,7 The concurrence is a measure of the degree of entanglement between the quantum dots, taking on a value of unity for maximal entanglement and 0 if the system is completely unentangled. Moderate values of entanglement, particularly at short times, are achieved and the level of agreement between non-Hermitian and Lindblad density matrix concurrences is good, indicating (especially at the short times) that the non-Hermitian model can describe entanglement.

Interestingly, if we maintain the plasmon loss and spontaneous emission terms but neglect the quantum dot dephasing, i.e., set $γ2*$ = 0, we obtain the result in Fig. 5(b), with excellent agreement between the non-Hermitian and density matrix models, both in terms of probabilities and concurrences. Note that a steady state is established in the non-Hermitian model and it agrees with the density matrix one. In the Lindblad master equation approach, dephasing only affects the coherences while leaving the populations unchanged. Since the non-Hermitian Schrödinger equation approach deals directly with the probability amplitudes, it will affect both the coherences and populations simultaneously. It is only in the limit that dephasing goes to zero and that both approaches will agree quantitatively.

Dephasing is included in the non-Hermitian model [Eq. (9)] as an imaginary part in some of the diagonal elements. This effectively acts as an additional source of dissipation, removing norm from system. In the full Lindblad master equation, however, dephasing does not just dissipate the number states; it also mixes the states and the effective dissipation due to dephasing is less than what the non-Hermitian model predicts. This is evidenced in Fig. 5(a) where the non-Hermitian dynamics are slightly below the density matrix dynamics.

We should point out that all the non-Hermitian model results of this subsection do not actually require numerical integrations since there is no driving and the dynamics is restricted to the one-excitation manifold. If the states in this one-excitation manifold are taken to be |s = 1, q₁ = 0, q₂ = 0⟩, |s = 0, q₁ = 1, q₂ = 0⟩, |s = 0, q₁ = 0, q₂ = 1⟩, the relevant 3 × 3 Hamiltonian matrix is

which is a generalization of the approach in the Appendix of Ref. 7 (see also Ref. 26) to include quantum dot dephasing and spontaneous emission. The eigenvalues and eigenvectors can readily be obtained, and the dynamics of a particular initial condition within the one-excitation manifold can easily be obtained. It is also clear that this can easily be extended to studying the dynamics of many more quantum dots interacting with a plasmon mode. Indeed, it is possible to write down the analytical solution for the N quantum dot extension of Eq. (14), also allowing for inhomogeneous couplings (see the  Appendix for details). As an example, we consider 50 quantum dots in resonance with a single plasmonic mode, as schematically illustrated in Fig. 6(a), with (i) homogeneous couplings [all dot/plasmon couplings equal to the value given in the lower parameter set of Table I (0.0167 eV)] and (ii) inhomogeneous couplings (dot/plasmon couplings randomly drawn from a normal distribution with mean 0.0167 eV and standard deviation 0.004 175 eV). We choose this example to study the effects of disorder on bi-partite entanglement properties of quantum plasmonic systems. Figure 6(b) displays the average bipartite concurrence²⁵ for both cases with the initial condition being one quantum dot excited. The magnitudes of these concurrences are much smaller than the two quantum dot case. This is due to the phenomenon of monogomy of entanglement,²⁷ where the amount of bipartite concurrence between any two systems of a set of systems decreases as the size of the set increases. Such states can still be maximally entangled (in a multi-partite sense), as has been shown for the W state.²⁸ In addition, Fig. 6 shows differences between the homogeneous and inhomogeneous coupling cases, with the transitory, maximum concurrences being somewhat larger for the homogeneous case and the inhomogeneous case exhibiting somewhat higher frequency oscillations. Nonetheless, the non-negligible concurrences seen in the presence of disorder are encouraging for practical studies, which inevitably involve some degree of disorder. Figures 6(c) and 6(d) display the excited state populations for quantum dots 2–50 (with quantum dot 1 being the initially excited one) for the homogeneous and inhomogeneous cases, respectively. Of course, all quantum dots act the same for the homogeneous case, whereas Rabi oscillations are only significant for quantum dots with sufficiently large coupling strengths in the inhomogeneous case. Similar trends are evident for the pair-wise concurrences in Figs. 6(d) and 6(f).

We now explore the role of disorder on the average bipartite concurrence for 50 quantum dots coupled to a single plasmonic mode. Figures 7(a) and 7(b) show the average bipartite concurrence for mean coupling strengths of $ℏg¯$ = 6.7 meV and 83.4 meV. The disorder parameter is characterized by the standard deviation σ_g = (0, 0.13, 0.25) $×g¯$⁠, varied from top to bottom in Figs. 7(a) and 7(b), with all results averaged over 100 instantiations. Figure 7(a) shows that the average concurrence does not change much as a function of disorder, demonstrating a surprising robustness to the coupling strength inhomogeneity. On the other hand, for a large enough mean coupling strength and relative disorder, as shown in Fig. 7(b), we observe a noticeable difference in the ensemble-averaged time-dependent concurrence. The average bipartite concurrence has different behavior in short and long times due to the dark and bright modes of the system.⁶ In the limit of zero dephasing, the dark mode persists indefinitely and may be considered as a quasi-steady-state of the non-Hermitian Hamiltonian, where N − 1 quantum dots form an approximate W many-body entangled state, where a W state is an equal superposition of all single-excited states.^7,28 The short time dynamics is dominated by bright modes, which exhibit fast oscillations that dampen out as a function of disorder. This result is highlighted in Fig. 7(c) where we plot the equilibration time as a function of the normalized disorder parameter. Here, we define the equilibration time as the time when the difference between the minimum in maximum inside a window including at least two maxima is below a cutoff of 0.000 23. This value was chosen by examining the difference between the minimum and maximum throughout the whole evolution and choosing a value where the decay of the difference dramatically decreases.

The origin of the short-time oscillations in Figs. 7(b) and 7(c) is due to two bright modes that are the eigenvectors of the coupled N-quantum dot and single plasmon system. In the short-time limit where the fast oscillations are observed, the wavefunction may be expanded as a superposition of two “bright” states (see the  Appendix), resulting in the approximate form for the concurrence,

where G is the square root of the sum of the squares of all the coupling strengths. By taking an ensemble average of this quantity, it can be shown that oscillations dampen out as a function of disorder, qualitatively consistent with our results in Fig. 7(b). It is very interesting to note that despite there being only one decay constant, (γ_pl + Γ)/2, when G is appropriately averaged over, we find what appears to be an effective increase in the decay rate as disorder is increased. Our results may be useful for sensing applications or other quantum applications that could make use of ensemble averages, e.g., arrays of sensors.

In general, there is a computational advantage to the non-Hermitian approach in the sense that for a problem with Hilbert space dimension M, the effort and storage scale as M for the non-Hermitian approach as opposed to M² for the Lindblad approach. In the above 50 quantum dot example, however, we made use of a symmetry and were able to focus on the single-excitation manifold leading to an effective Hilbert space of dimension N + 1 = 51; a corresponding Lindblad treatment while scaling as (51)² (i.e., a factor of 51 more effort/storage) is still relatively easy to carry out for N = 50. If we went much higher in N or could not invoke any symmetries, computational advantages would become more useful.

We have investigated a simple non-Hermitian model to describe quantum dot/plasmon interactions. We found that it yielded generally very good results in the linear optical excitation regime for models of one and two quantum dots coupled to a plasmonic structure. It led to poor results in the limit of cw driving as intensity was increased, however, and could not describe saturation effects in this limit. Nonetheless, we were also able to show that the model could describe non-trivial coherences and entanglement in the un-driven case, i.e., scenarios wherein one imagines an optical process has already excited a quantum dot and energy transfer can occur via the plasmon to excite the other quantum dot. Such non-Hermitian models will be useful for studying cases involving many quantum dots coupled to plasmonic structures that cannot be easily simulated with more complete density matrix approaches.

In the future, we plan to investigate approaches to remedying the failure of the non-Hermitian model in the high intensity cw case. One avenue is to exploit the fact that the model represents the first stage of the stochastic Schrödinger equation.^18,19 Another avenue is to attempt to incorporate non-linearities explicitly into the equations that include aspects of gain, as in the optical Bloch equation work in Refs. 29 and 30.

Based on our disorder study involving 50 quantum dots coupled to a plasmon, which showed promising properties can be achieved with increasing disorder, we also plan to explore such systems in the realms of quantum sensing and in the generation of novel quantum states. Further work will explore the role of exceptional points in these systems, as well as conditioned or post-selected measurements.^11,31

This work was performed at the Center for Nanoscale Materials, a U.S. Department of Energy Office of Science User Facility, and supported by the U.S. Department of Energy, Office of Science, under Contract No. DE-AC02-06CH11357.

Here, we generalize the analytical treatment of Ref. 7 to allow for effective quantum dot decays. The underlying single excitation manifold is given by the set of states {|s = 1, q₁ = 0, q₂ = 0, …, q_N = 0⟩, |s = 0, q₁ = 1, q₂ = 0, …, q_N = 0⟩, …, |s = 0, q₁ = 0, q₂ = 0, …, q_N = 1⟩}. The N + 1 × N + 1 non-Hermitian Hamiltonian matrix within this manifold is then

corresponding to a plasmon mode (decay constant γ_pl) interacting with N two-level quantum dot systems with couplings g₁, g₂, …, g_N and effective decay constant $Γ=γ1+2γ2*$⁠. Considering successive solutions of the corresponding determinants for the eigenvalues for N = 2, 3, etc., cases, it is easy to see that the eigenvalues, {λ_k}, of Eq. (1) are given by

where $G=g12+g22+⋯+gN2$⁠. The corresponding eigenvectors of Eq. (1), v^k$=v1k,…,vN+1kT$⁠, represent the eigenstates

The first N − 1 eigenvectors can be obtained with the procedure introduced in the Appendix of Ref. 7. That is, the matrix eigenvalue equation for these eigenvectors, even with quantum dot decays, leads to all these eigenvectors being orthogonal to the vector g = (0, g₁, g₂, …, g_N)^T. Therefore, a Gram–Schmidt procedure applied to g and a set of N − 1 vectors with random coefficients can lead to a suitable realization of the degenerate vectors, v^k, for k = 1, 2, …, N − 1. While Eq. (A2) shows that the corresponding eigenvalues have negative imaginary parts corresponding to loss, the magnitude of this loss in the plasmonics context is less than the loss associated with λ_N and λ_N+1, which involves γ_pl. For this reason, we term these first N − 1 states “dark” states. Regarding the two remaining “bright” states, the eigenvectors are

All the eigenvectors so-constructed are orthogonal. A given initial state, |ψ(t = 0)⟩, then evolves in time according to

where $nk=∑jvjk2$ (assuming some of the v^k have not been otherwise normalized) and b_k = ⟨ϕ_k, ψ(0)⟩$≡∑j=1N+1vjkdj$⁠, where d is the vector representation |ψ(t = 0)⟩in the single excitation manifold basis. Note that there are no complex conjugates involved in the normalization condition or in the definition of ⟨w, x⟩, which is a simple way of correctly incorporating the bi-orthogonal nature of the problem since the left eigenvectors are the transposes of the right eigenvectors for such a complex symmetric Hamiltonian matrix.

To calculate the average concurrence, we use the definition for concurrence for pure states: C = 2|ab − cd|, for a general two-qubit (or quantum dot) wavefunction, |ψ(t)⟩= a|00⟩+ b|01⟩+ c|10⟩+ d|11⟩. In the single excitation manifold, the concurrence is given by C = 2|cd|, where “c” and “d” are the probability amplitudes of qubit 1 in the excited-state and qubit 2 in the excited-state. For an N-qubit system, the concurrence between any two qubits i and j may be written as C_ij = 2|c_i c_j|, where c_i (c_j) corresponds to the probability amplitude of qubit i (j) in the excited-state. It can also be shown that the average concurrence may be written as C_avg(t) = $2N(N−1)∑ijcitcjt$⁠. Noting that disorder primarily affects the short-time behavior that is dominated by the bright states, we make the approximation that probability amplitude for the jth qubit is given by (assuming qubit 1 is initially excited)

where e_j is a unit vector corresponding to the jth qubit (one at the qubit location and zero everywhere else). This quantity can be shown to be proportional to

Note that the above approximate relation for the amplitude of qubit j, due to the two bright states, is the same for all j, and thus, the average concurrence inferred from it is simply proportionate to |c_jc_j|, which leads to Eq. (15). It is worth noting that the t = 0 limit is incorrect, since it suggests a high concurrence and we know that the concurrence is zero at t = 0. To obtain the correct behavior, the dark state contributions need to be included. Nonetheless, we find that this relation describes the correct oscillation dynamics that leads to the relaxation of the bright states to the dark states.