text:
\section*{Introduction}

	The electroluminesence (EL) of individual quantum emitters coupled to metallic electrodes has been investigated extensively
	since the first scanning tunneling microscope induced luminescence (STML) experiments were performed \cite{gimzewski1988photon,shamai2011spectroscopy}.  Through precise control of an STM probe's position, both the emitter--probe coupling and the resonant frequency of collective motion of the metallic electrons (plasmons)
	confined in the nanocavity formed between the probe and substrate can be tuned. By adjusting the plasmon frequency, the coupling between particle-hole excitations (excitons) on the quantum emitter and the plasmons can be controlled, leading, for instance, to the observation of plasmon enhanced photon emission
	\cite{gimzewski1989enhanced,berndt1991inelastic}.
	%

	This unprecedented control
	makes STML systems ideal both
	for exploring fundamental aspects of the nonequilibrium electro-optical response of quantum emitters
	and as a testbed to develop quantum-enhanced device technology, e.g. those related to biological sensing \cite{lee2007exciton}, photovoltaic energy conversion \cite{ferry2010design}, or non-classical light generation \cite{chang2006quantum, tame2013quantum,monroe2002quantum, khitrova2006vacuum,raimond2001manipulating}.

	We focus on molecular emitters in particular

	since they can be engineered with the specific emission profiles, dipole moments, wavelengths, and symmetries necessary to harness uniquely quantum resources which may be useful in the development of novel opto-electronic devices  \cite{nothaft2012electrically, lounis2005single,aradhya2013single}.

	When a voltage bias causes the source and drain electrodes' chemical potentials to align with unoccupied and occupied molecular states, respectively,
	a tunnel current and subsequent molecular exciton are produced.  If the exciton decays radiatively, the resulting EL encodes the specific electronic and vibrational state of the molecule.  Molecular vibrational states have been observed in STML experiments of single porphryin molecules \cite{qiu2003vibrationally,wu2008intramolecular, dong2004vibrationally,dong2010generation,chen2010viewing} and fullerene C$_{60}$ and C$_{70}$ clusters \cite{cavar2005fluorescence,rossel2009plasmon}.

	%

	Dong {\em et al.} observed molecular hot-luminescence (HL) from excited vibrational modes in tetraphenylporphyrin (TPP) molecules weakly coupled to metallic electrodes \cite{dong2010generation}.  Their data are a direct observation of the strong dependence of the EL on the resonant frequency of the localized nanocavity plasmons.    In addition, their report of a violation of Kasha's rule, which states that the lowest vibrational transitions should dominate the molecular flourescence, indicates a strong enhancement of the

	spontaneous emission rate (i.e. a Purcell enhancement \cite{purcell1946spontaneous}) caused by the formation of the nanocavity \cite{le2007surface, chen2015molecular}.
	Interestingly, above threshold HL (i.e. $eV < \hbar \omega$) was also observed in TPP junctions.\cite{dong2010generation}  In the weak coupling limit, this effect doesn't appear to be described using a classical plasmonic field,\cite{tian2011electroluminescence} although it may be explained when higher-order  electron-plasmon scattering processes are included\cite{kaasbjerg2015theory}.

	In the study of quantum electrodynamics, Purcell enhancement is a signature of the {\em weak coupling} regime between coupled quantum emitters and optical modes.  As the coupling strength is increased
	there is a transition into the {\em strong coupling} regime, where energy transfers coherently between the emitter and field modes, giving rise to an observable Rabi splitting between the joint emitter-field states.  Systems operating in the strong coupling regime allow for the observation of quantum effects, including single-atom lasing, single photon generation, and all-optical single photon switching \cite{khitrova2006vacuum, mckeever2003experimental, solomon2001single,kimble1998strong,raimond2001manipulating}.
	%

	Coupled molecular excitons and plasmons form joint states known as plexcitons \cite{savasta2010nanopolaritons,chang2006quantum,tame2013quantum,manjavacas2011quantum}.  Molecular plexcitonic states with Rabi splittings up to several hundred meV have been observed  \cite{bellessa2004strong,dintinger2005strong,hakala2009vacuum,schwartz2011reversible, guebrou2012coherent,fofang2008plexcitonic,vasa2013real,melnikau2016rabi,wang2016role}, motivating a detailed
	investigation into the influence of quantum dynamics, chemical structure, many-body interactions, plexcitonic dynamics, and loss mechanisms on the optoelectronic response of these system.

	In this article, we develop a state-based quantum master equation (QME) approach to investigate the EL of molecules in both the weak and strong plexcitonic coupling regimes.
	We first derive an effective multi-state Jaynes-Cummings model for the molecule and quantum plasmon modes, and use the QME framework to describe finite tunneling currents, radiative and non-radiative exciton decay paths, vibrational damping, and finite plasmon lifetimes.   Although similar methods have been used to investigate plasmon-enhanced EL and transport-induced EL in STML systems before \cite{tian2011density,tian2011electroluminescence,zhang2013plasmon,chen2015molecular}, we extend these works to describe the quantum optical regime including a full quantum many-body description of the molecule, plasmon modes, electrodes and their couplings.  As a first application, we simulate the EL of a TPP molecule coupled to a single quantum plasmon mode for several voltages, plexcitonic coupling strengths, and detunings.

	\begin{figure}[tb]
		\centering
		\includegraphics[width=2.8in]{schematic_alone_fig1_v3.eps}

		\caption{A schematic diagram of an STML experiment including the STM probe, molecule, and substrate.  Molecular excitons generated by a finite tunnel current couple to nanocavity plasmons and may decay radiatively as EL.}

		\label{fig:schematic}
	\end{figure}

	\section*{Theoretical Model}

	%

	%

	%

	We consider open quantum systems
	composed of a molecule coupled to electromagnetic field modes, metallic leads (e.g. the substrate and STM probe), and vibrational modes subject to applied voltages and temperature gradients.
	A schematic of the STM-based experiments we consider is shown in Fig~\ref{fig:schematic}.
	The Hamiltonian corresponding to this system may be partitioned as

	\begin{equation}
	H = {H_{\rm mol}} +  {H_{\rm leads}} + H_{\rm vib} + H_{\rm EM} +
	  H_{\rm tun}  + H_{\rm mol - vib} + {H_{\rm mol - EM}},
	\label{eq:Hamiltonian}
	\end{equation}

	where $H_{\rm mol}$ is the molecular Hamiltonian, and the independent lead, vibrational, and electromagnetic baths are described by
	\begin{align}
	{H_{\rm leads}} &= \sum_\alpha \sum_{k \in \alpha } {{\varepsilon _{k\sigma }}c_{k\sigma }^\dag {c_{k\sigma }}},  \\
	{H_{\rm vib}} &= \sum_\alpha \sum_{l \in \alpha } {{\hbar\Omega_{l }}b_{l\sigma }^\dag {b_{l\sigma }}}, \\
	H_{\rm EM} & = \sum_j \hbar \omega_j a^\dagger_j a_j
	\end{align}
	respectively, where $c_{k}$ annihilates an electron in lead mode $k$ with dispersion $\varepsilon_{k\sigma}$ and spin $\sigma$, $b_l$ annihilates a vibrational excitation (phonon) in mode $l$ with energy $\hbar \Omega_l$, and $a_j$ annihilates a photon in mode $j$ with energy $\hbar \omega_j$.

	The tunnel coupling, vibrational mode couplings, and electromagnetic field couplings are given by
	\begin{align}
	H_{\rm tun} & = \sum_\alpha \sum\limits_{k \in \alpha,n\sigma } {{V_{nk}}d_{n\sigma }^\dag {c_{\sigma k}} + H.c.}, \nonumber \\
	H_{\rm mol-vib} &= \sum_\alpha \sum\limits_{l \in \alpha,n\sigma } {{W_{nk}}d_{n\sigma }^\dag {b_l} + H.c.},   \nonumber \\
	H_{\rm mol-EM} &=  \int d{\bf r} J({\bf r})\cdot {\bf A}({\bf r})

	\end{align}
	respectively, where $V_{nk}$ is the tunneling matrix element, $d_{n\sigma}$ annihilates an electron in molecular state $n$ with spin $\sigma$, and $W_{nk}$ is the vibrational coupling between molecular orbital $n$ and lead mode $k$.

	The current density $J({\bf r})$ couples to the vector potential ${\bf A}({\bf r})$, which for quantized electromagnetic modes is given by \cite{scully1997quantum}
	\begin{equation}
	{\bf A}_j({\bf r}) = \sqrt{\frac{\hbar}{2\varepsilon_0 \omega_j V}} {\eta}_j({\bf r}) \left(a_j + a_j^\dagger\right)
	\end{equation}
	where ${\bf A} = \sum_j {\bf A}_j$, $\eta({\bf r})$ is the product of the polarization vector and a function describing the spatial profile of the field, $V$
	is the effective mode volume, and $\hbar\omega_j$ is the energy of mode $j$.

	%

	%

	%

	\subsection*{Quantum Master Equation}
	In general, the system described by Eq.~\ref{eq:Hamiltonian} cannot be solved exactly.  To proceed, we utilize a state-based quantum master equation (QME) approach, where the quantum dynamics of the joint molecule and plasmon system are treated exactly, while the other macroscopic degrees of freedom are traced over using a coarse-graining procedure.  Within the QME framework, the Liouville equation for the reduced density matrix of the system is given by
	\begin{equation}
	\dot{\rho} =  - \frac{i}{\hbar}\left[ {{H_0},\rho } \right] +  \left({\cal L}_{\rm tun} + {\cal L}_{\rm damp} + {\cal L}_{\rm deph} \right) \rho,
	\label{eq:Liouvillian}
	\end{equation}
	where $\rho$ is the
	density operator,

	and the Liouvillian superoperators $\SO_{\rm tun}$, $\SO_{\rm damp}$, and $\SO_{\rm deph}$ describe the nonhermitian evolution of the system due to quantum tunneling, damping of the populations and coherences of states, and pure dephasing, respectively. Once the density matrix is determined, expectation values of observables may be calculated using $\langle {\cal O} \rangle = \Tr{\rho(t) {\cal O}(t)}$.

	%

	%

	After coarse graining, our free-system Hamiltonian $H_0$
	is composed of three terms: the molecular Hamiltonian, quantum plasmon modes, and their couplings.  The Hamiltonian of each term is given by

	\begin{align}
	\label{eq:Hfree_start}

	H_{\rm mol} &= \sum\limits_{n,m,\sigma}  {H^{(1)}_{nm}}d_n^\dag {d_m}  + \!\!\! \sum_{ijnm, \atop{\sigma \sigma'}} \!\!\! {{\frac{U_{ijnm}}{2}}d_{i\sigma }^\dag d_{j\sigma '}^\dag {d_{m\sigma '}}{d_{n\sigma }}}  \nonumber \\ & \phantom{abc} + \sum_{l} \hbar \tilde{\Omega}_l \tilde{b}_l^\dagger \tilde{b}_l + \sum_{l,n} \lambda \hbar\tilde{\Omega}_l (\tilde{b}_l + \tilde{b}_l^\dagger) d_n^\dagger d_n,   \\
	H_{\rm plas} &= \sum\limits_j {\hbar \tilde{\omega}_j \tilde{a}_j^\dag {\tilde{a}_j}},    \\
	H_{\rm int} &= \sum_{n,m,j} \hbar g_{nm}^j d^\dagger_n d_m \left(\tilde{a}_j + \tilde{a}_j^\dagger \right),
	\label{eq:Hfree_end}
	\end{align}
	respectively,

	%
	where $H^{(1)}_{nm}$ is the one-body portion of the molecular Hamiltonian, which may be renormalized by classical electrostatic\cite{barr2012effective} (e.g. image charge) or vibrational effects induced by the electrodes; $U$ is the Coulomb integral, $\tilde{b}_l$ annihilates a (renormalized) phonon in mode $l$ with energy $\hbar \tilde{\Omega}_l$; $\lambda$ is the electron-phonon coupling; and $\tilde{a}_j$ annihilates a plasmon in mode $j$ with energy $\hbar \tilde{\omega}_j$.

	For systems we consider, the dipole approximation of the electromagnetic coupling is sufficient. In this approximation

	the plexcitonic coupling parameter is given by \cite{loudon2000quantum}

	\begin{equation}

	\hbar g^j_{nm}  = \sqrt{\frac{\hbar \omega_j}{2\varepsilon_0 V}} \mu_{nm} u_j(x_0)
	\label{eq:g}
	\end{equation}
	where $j$ is the plasmon mode index, $n$ and $m$ are level indices, $V$ is the mode volume, $\mu_{nm} =-e\langle n \left| \vec{r}  \right| m\rangle$
	is the transition dipole matrix element, $\omega_j$ is the mode's angular frequency, $\varepsilon_0$ is the permittivity of free space, and $u_j(x_0)$ is the mode function evaluated at the emitter's position $x_0$.

	Although we have expressed the molecular Hamiltonian in terms of electron and phonon operators, Eqs.~\ref{eq:Hfree_start}-\ref{eq:Hfree_end} are essentially a multi-state Jaynes-Cummings model where the state energies and matrix elements can be found using a variety of methods (e.g. via exact diagonalization, density functional theory, etc.).

	\subsubsection*{Quantum Transport}

	Electron transport involves the addition and removal of charges which maintain a degree of phase coherence as they traverse a junction. We consider systems
	in which the molecules and electrodes are deliberately decoupled (e.g. via the growth of insulating layers on the metallic  substrate\cite{qiu2003vibrationally}, or by depositing several molecular monolayers  \cite{dong2004vibrationally}) such that the individual molecule's emission is not quenched by interactions with the metallic electrodes\cite{baffou2008molecular,marinica2013plexciton}.

	In this

	regime, the coherence time of electrons on the molecule are short compared to tunneling time, allowing us to neglect the excitation of coherent superposition states and instead describe the transport as a simple kinetic process.
	\cite{beenakker1991theory,harbola2006quantum}.

	Following an expansion of the Liouville equation for the time evolution of the density matrix to second-order in $H_{\rm tun}$,
	the master equation for tunneling is given by \cite{seldenthuis2010electroluminescence, tian2011density,chen2015molecular}
	\begin{equation}
	{\cal L}_{\rm tun}\rho = \sum_{ij} \left[R_{i \rightarrow j}\sigma_{ji}\rho \sigma_{ij} - R_{j \rightarrow i}\sigma_{jj}\rho \sigma_{jj} \right],
	\end{equation}
	where $R_{i\to j}$ is the charging rate between the $N$-particle state $i$ and the N+1-particle state $j$,  $R_{j\to i}$ is the discharge rate between states $j$ and $i$, and
	$\sigma_{ij}$ is a matrix in the free system's state space with element $(i,j)=1$ and all other elements equal to $0$.
	%

	%

	The electronic tunneling rates are given by
	\begin{align}
	R_{i \to j} &= F_{ij} \sum\limits_{\alpha} {\tilde \Gamma _{ij}^\alpha {f_\alpha }\left( {\left[ {E_j - E_i} \right] } \right)},  \\ \nonumber
	{R_{j \to i}} &= F_{ji}\sum\limits_{\alpha} \tilde \Gamma _{ji}^\alpha \left\{ {1 - {f_\alpha }\left( {\left[ {E_j - E_i} \right]} \right)} \right\},
	\end{align}
	where $F_{ij}$ are the Franck-Condon factors (i.e. the overlap between nuclear wave functions), and $f_\alpha(E)=(1+\exp[(E-\mu_\alpha)/k_{\rm B}T_\alpha])^{-1}$ is the Fermi-Dirac distribution for lead $\alpha$ with temperature $T_\alpha$ and chemical potential $\mu_{\alpha}$.  The chemical potential of the source and drain leads are given by $\mu_S=E_f-e\alpha V_b$ and $\mu_D=E_f-e(1-\alpha)V_b$, respectively, with Fermi energy $E_f$, electron charge magnitude $e$, voltage symmetry $\alpha$, and bias voltage $V_b$. We assume a symmetric potential drop, where $\alpha=0.5$.

	The effective tunnel coupling between electrode $\alpha$ and the molecule is given by \cite{bergfield2009many}
	\begin{equation}
	{\tilde \Gamma}_{ij}^\alpha  = \Tr{ {{\Gamma ^\alpha }\left( {{E_j} - {E_i}} \right)C\left( {i,j} \right)}},
	\label{eq:eff_gamma}
	\end{equation}
	where the bare electron tunneling rate matrix
	$\Gamma ^\alpha_{nm}\left( E \right) = 2\pi/\hbar \sum_{k\sigma  \in \alpha } {{V_{nk}}V_{mk}^*\delta \left( {E - {\varepsilon _{k\sigma }}} \right)}$
	is dressed by the many-body renormalization factors \cite{bergfield2009many}
	\begin{equation}
	{\left[ {C\left( {i,j} \right)} \right]_{n\sigma ,m\sigma '}} = \left\langle j \right|d_{n\sigma }^\dag \left| i \right\rangle \left\langle i \right|{d_{m\sigma '}}\left| j \right\rangle.
	\label{eq:manybody}
	\end{equation}
	As indicated by Eq.~\ref{eq:eff_gamma}, both the relative phase and magnitude of the many-body factors influence the effective tunneling rates and therefore the transport and optical response of systems with multiple states.  In addition to the many-body wave function normalization, where the total resonance width of a molecular state is reduced by a factor of $1/N$ ($N$ is the number of atomic orbitals), strong correlations can also lead to an exponential suppression of these terms.

	\subsubsection*{Damping and Dephasing}

	When a free system interacts with the environment,
	an initially excited state can decay via a number of irreversible damping processes.  We account for these loss mechanisms with the
	composite Liouvillian operator $\SO_{\rm damp} = \SO_{\rm rad} + \SO_{\rm cav}+\SO_{\rm vib} $, which describes radiative decay processes,
	the finite lifetime of the nanocavity plasmons, and vibrational relaxation processes, respectively.
	%

	Assuming Markovian baths, $\SO_{\rm rad}$
	can be expressed as a Lindblad master equation \cite{scully1997quantum, seldenthuis2010electroluminescence, tian2011density,tian2011electroluminescence}
	\begin{equation}

	{\cal L}_{\rm rad}\rho  =   \! - \! \sum_{i,j} \!\! \frac{\gamma^{j \to i}_{\rm rad}}{2} \left( \sigma _{ji}\sigma _{ij}\rho  - 2\sigma_{ij}\rho \sigma_{ji} + \rho \sigma_{ji}\sigma_{ij}  \right),
	\end{equation}
	where $\gamma^{j \to i}_{\rm rad}$ is the
	radiative coupling rate between electronic levels $i$ and $j$.  The finite plasmon lifetime is included via the phenomenological decay rate $\kappa_j$ and

	master equation

	\begin{equation}
	\SO_{\rm cav} \rho =  - \sum_j {\kappa_j  \over 2 }\left( {{\tilde{a}_j^\dag }\tilde{a}_j\rho  - 2\tilde{a}_j\rho {\tilde{a}_j^\dag } + \rho {\tilde{a}_j^\dag }\tilde{a}_j} \right),
	\end{equation}
	where $\tilde{a}_j$ annihilates a plasmon in mode $j$.

	The intraband vibrational damping may be described by \cite{seldenthuis2010electroluminescence}
	\begin{align}
	{\SO_{\rm vib}}\rho
	& =  - \gamma_{\rm vib} \sum\limits_i  \left[{ {{\sigma _{ii}}\rho {\sigma _{ii}} - P\left( {\Omega_{\rm vib}^i} \right)\sum\limits_j {{\sigma _{ij}}\rho } {\sigma _{ji}}} } \right],
	\end{align}
	where $\gamma_{\rm vib}$ is the vibrational coupling rate, and $P(\Omega) = e^{-\hbar \Omega / k_B T}/{\cal Z}$, with the partition function ${\cal Z}=\sum_k {e^{-\hbar \Omega^k_{\rm vib} / k_BT}}$.  The states labeled $i$ and $j$ belong to the same electronic manifold.  Although $\gamma_{\rm vib}$ is typically several orders of magnitude larger than the radiative relaxation rate, $\gamma_{\rm rad}$  may be enhanced (e.g. by placing an emitter in a cavity) to exceed the vibrational decay rate, resulting, for instance, in the observation of HL from excited vibrational states
	\cite{dong2010generation, tian2011electroluminescence,chen2015molecular}.

	In addition to relaxation, where excitations are transferred from the system to modes of the environment, pure dephasing is also possible, in which the populations are unaffected but their coherence is reduced.

	The Liouvillian term

	describing
	dephasing

	is given by \cite{johansson2005surface}
	\begin{equation}
	\SO_{\rm deph} = - \gamma_{\rm deph} \sum_{i,j} \left(\sigma_{ij} \rho_{ij} + \sigma_{ji} \rho_{ji} \right)
	\end{equation}
	where $\gamma_{\rm deph}$ is the dephasing rate.

	%
	%

	%

	The relationship between $\kappa$, the total state coupling $\Gamma_{\rm total}$, and the plexcitonic coupling $g$ distinguishes the strong ($g \gg \Gamma_{\rm total}, \kappa$) and weak ($g \ll \Gamma_{\rm total}, \kappa$) coupling regimes.  According to Eq.~\ref{eq:g} the strong coupling regime may be achieved by increasing the transition energy, increasing the molecular dipole moment, or decreasing the effective mode volume accessible to the plasmons  (e.g. by patterning the substrate to exhibit a reduced density of modes).

	%

	%
	%

	%
	%

	\begin{figure}[tb]
		\centering

		\includegraphics[width=3.5in]{tpp_energy_level_v2.eps}
		\caption{The energy level diagram for the charge neutral manifold of tetraphenylporphyrin (TPP) molecule and coupled quantum plasmon modes.  The relevant (Q-band) electronic states $S_0$ and $S_1$ are shown with their associated vibrational levels.  The observed optical gap of TPP is 1.89eV, with a vibrational state spacing of 0.16eV.\cite{dong2010generation}  Based on comparison with experiment, the source and drain tunnel couplings are set to $\hbar \Gamma^{\rm S}=\hbar \Gamma^{\rm D}$=16.4$\mu$eV, while the vibrational damping $\gamma_{\rm vib}$, and radiative decay  $\gamma_r$, are consistent with a vibrational lifetime of 2ps and a radiative lifetime of 2ns.\cite{dong2010generation}

			We set the plasmonic decay to $\hbar \kappa=10meV$ and plexcitonic coupling $g$ is taken as an adjustable parameter.}
		\label{fig:energy_level}

	\end{figure}

	\begin{figure*}[htb]
		\centering

		\includegraphics[width=.44\linewidth]{EL_g_0_kappa_10_delta_0_T_80.eps}
\includegraphics[width=.44\linewidth]{EL_g_10_kappa_10_delta_0_T_80.eps}
\includegraphics[width=.44\linewidth]{EL_g_20_kappa_10_delta_0_T_80.eps}
\includegraphics[width=.44\linewidth]{EL_g_30_kappa_10_delta_0_T_80.eps}
\vspace{-0.5cm}
		\caption{The calculated STML of a TPP molecule coupled to a single plasmon mode is shown as a function of photon energy for four plexcitonic coupling strengths $\hbar g_{00}$=0meV, 10meV, 20meV, and 30meV at three bias voltages $V_b$=1.8V, 1.9V, and 2.5V. In the decoupled cavity limit ($\hbar g_{00}=0$meV), the molecular vibrational spectrum is recovered.  As $g_{11}$ is increased plexcitonic states form giving rise to the characteristic split peaks in the EL.  Although all radiative transitions have non-zero coupling to the plasmonic field, the resonant $Q(0,0)$ transition dominates and the plexcitonic splitting is $\sim2\hbar g_{00}$.  The coherent mixing of multiple states gives rise to the asymmetric peak structure and the above threshold hot HL observed when $V_b$=1.8V; even for junctions operating in the weak coupling regime.  For this system, the strong coupling boundary is $\hbar g_{00} > \sim$15meV.  Calculations are for junctions operating at $T$=80K, to be consistent with measured STML  spectra.\cite{dong2010generation}}
		\label{fig:resonant_TPP_vs_g}
	\end{figure*}

	\section*{Results and Discussion}

	As a first test of our theory, we investigate the STML of a tetraphenylporphyrin (TPP) molecule coupled to a single quantum plasmon mode.  In general, the spontaneous emission spectrum may be found from the density matrix though the use of the quantum regression theorem \cite{gardiner2004quantum,laussy2009luminescence}.  Since our model for the TPP molecule is constructed {\em ad hoc} from experimental data, we consider it to be a sum of Lorentzians given by
	\cite{tian2011density,tian2011electroluminescence,chen2015molecular}

	\begin{equation}
	I(\omega) =  \frac{1}{2\pi} \sum_{i,j} \frac{\gamma_{\rm rad}^{j\rightarrow i} \rho^{jj}_{ss}}{(\omega - \omega_{ji})^2 +\tau_{ji}^2}
	\end{equation}
	where $\gamma_{\rm rad}^{j\rightarrow i}$ is the radiative decay rate between states $j$ and $i$,  $\rho_{ss}$ is the steady-state solution of Eq.~\ref{eq:Liouvillian}, and $\tau_{ji}$ is the full width at half maximum (FWHM) of the EL.  We assume that the plexcitonic EL linewidths are well approximated by the zero-detuned single atom vacuum Rabi splitting linewidths $\tau_{ji}=(2\tau_{ji}^0 + \kappa) / 2$, where $\tau_{ji}^0$=0.05eV is the FWHM extracted from experiments operating in the weak coupling limit \cite{dong2010generation,tian2011density}.

	A schematic of the TPP molecule and the energy level energy diagram for the TPP system are shown in Fig.~\ref{fig:energy_level}.  The $Q$-band
	energy gap of TPP (i.e. the $Q(0,0)$ transition, $S_1(0)\rightarrow S_0(0)$) is 1.89eV, with a vibrational level spacing of $0.16$eV. \cite{dong2010generation, gouterman1978optical}  We consider the TPP molecule coupled to a single quantum plasmon mode with energy $\hbar \tilde{\omega}$ = 1.89eV -  $\Delta$, where $\Delta$ is the detuning parameter. The source and drain tunnel coupling rates ($\Gamma^{\rm S}=\Gamma^{\rm D}=0.4\times10^{10}s^{-1}$), vibrational decay rates ($\gamma_{\rm vib}=0.5\times10^{12}s^{-1}$), and radiative decay rates ($\gamma^{j \rightarrow i}_{\rm rad} = F_{ji} / 2ns$)
	are established by comparison with experiment. \cite{tian2011density,tian2011electroluminescence,dong2010generation}
	The transition rates between vibrational levels and the transition dipole moments are scaled by the appropriate Franck-Condon factors $F_{ij}$, which are found using the harmonic approximation \cite{keil1965shapes} with a Huang-Rhys parameter $S=0.61$ (See Supporting Information).\cite{tian2011electroluminescence}  We assume that the radiative lifetime of all states are equal (i.e. $\gamma_{\rm rad}^{\rm HL}=\gamma_{\rm rad}$), and calculations were performed using a modified version of QuTIP. \cite{johansson2012qutip}

	\begin{figure*}[htb]
		\centering

		\includegraphics[width=.48\linewidth]{EL_g_10_kappa_10_delta_20_T_80.eps}
		\includegraphics[width=.48\linewidth]{EL_g_20_kappa_10_delta_20_T_80.eps}

		\caption{The calculated STML of a TPP molecule coupled to a single plasmon in the weak (left panel) and strong (right panel) coupling regimes for several detuning values.  The EL is enhanced by detuning the plasmon and molecular excitation resonances since the plasmonic nanocavity lifetime is limited by $\kappa$.  In the strong coupling regime ($|g_{00}|/\kappa > 0.25$), detuning can be used to favor the upper or lower plexcitonic state.  Simulations are for junctions operating at $T$=80K.
		}
		\label{fig:EL_detune}
	\end{figure*}

	For a resonant plasmon, the plasmonic decay rate $\kappa$ may be expressed in terms of the mode energy $\omega$ and quality factor $Q$ as $\kappa = \omega / Q$.  Using reports for other metal-insulator-metal nanostructures, \cite{yang2012ultrasmall} we find that values of $Q \approx 100$ are reasonable.  Given that a plasmonic mode's lifetime (and the cavity quality factor) can vary significantly for different junction designs, probe positions, and substrate materials, we consider $Q$=189 in our calculations, such that $\hbar\kappa$=10meV.

	The calculated STML of a TPP junction for four plexcitonic couplings and three bias voltages are shown as a function of photon energy in  Fig.~\ref{fig:resonant_TPP_vs_g}.

	In the decoupled cavity limit (top left panel, $\hbar g_{00}$=0meV), the experimentally observed spontaneous emission peak structure \cite{dong2010generation} is recovered, where the three peaks corresponding to the $Q(0,0)$, $Q(0,1)$, and $Q(0,2)$ transitions of the TPP molecule.  As $g_{00}$ is increased, the molecular exciton and cavity plasmon states mix, forming plexcitonic states separated in energy by $\sim 2 \hbar g_{ij}$.  In the weak coupling regime the radiative transition rate is enhanced via the Purcell effect, making HL possible, while in the strong coupling regime energy exchanges coherently between plasmonic and molecular states resulting in a characteristic split peak structure of the EL.  The boundary between weak and strong coupling regimes is defined by $|g_{00}|/|\kappa_t| = 0.25$, where for the TPP junction $\kappa_t \approx \kappa + \tau_{ij}^0$ = 60meV and therefore $\hbar g_{00} = \sim$15meV.

	The STML of a TPP junction operating in the weak coupling regime is shown in the top right panel of Fig.~\ref{fig:resonant_TPP_vs_g}, where a HL peak at 2.05eV is visible for all bias voltages.  Although the peak was observed experimentally at 1.8V, where $eV_b$ is less than the excitation energy of the molecule,\cite{dong2010generation} it was not seen at this voltage in previous calculations using classical plasmonic fields.\cite{tian2011electroluminescence}  This implies that Purcell enhancement alone can't explain the measured HL.

	The strong coupling STML of the TPP junction
	is shown in the lower two panels of Fig.~\ref{fig:resonant_TPP_vs_g}, where the $Q(0,0)$ and $Q(0,1)$ peaks have split into two peaks.  The asymmetry of these plexcitonic peaks stems from the influence of multiple detuned resonances, where the plexcitonic coupling between levels is reduced by the appropriate Franck-Condon factor. The $Q(0,2)$ and HL peaks are also split but can't be identified with $\kappa=10meV$ due the reduced effective couplings.

	Our simulations show that HL is suppressed when the off-diagonal coupling terms $g_{i \neq j}$ are reduced or, as expected,

	 when the vibrational relaxation rate is increased.  This suggests that the HL peaks are a consequence of the (weak) coherent dynamics between off-resonant states.  Although the above-threshold emission has been explained in terms of higher-order many-body processes \cite{tobiska2006quantum,schneider2010optical, kaasbjerg2015theory}, our calculations support an additional physical explanation in which the tunnel current pumps energy into cavity modes via the nascent plexcitonic states.  As shown in the lower panels of Fig.~\ref{fig:resonant_TPP_vs_g}, when the plexcitonic coupling is increased the detuning between states is reduced and, consistent with our argument, the above-threshold emission is enhanced.

	Finally, we consider the influence of the detuning between molecular transitions and  cavity plasmon resonance energies on the EL.  Physically, detuning can be controlled by adjusting the STM probe's height above the substrate.  As evidenced by the peak at 2.05eV shown in the left panel of Fig.~\ref{fig:EL_detune}, Purcell enhancement in the weak coupling regime results in HL for all detunings.  When the cavity is blue or red detuned relative to the $Q(0,0)$ transition of TPP, the STML spectral weight is shifted towards higher or lower energy peaks, respectively.  Since the TPP junction supports a finite current,
	energy is constantly (albeit weakly) pumped into off resonant cavity plasmon modes, giving the observed shift to the spectrum.  The EL is {\em increased} by detuning since detuning reduces the molecule's effective coupling to non-radiative plasmonic decay paths.

	In the strong coupling regime, shown in the right panel of the same figure, the blue and red detuned plasmon modes again shift the STML spectral weight up or down in energy, respectively.  However, in this regime the strongest peaks are split into distinguishable plexcitonic resonances.  In addition to the characteristic split peak EL, plexcitonic states and the onset of the strong coupling regime can also be identified via this distinct spectral weight shift with detuning.

	\section*{Conclusions}

	We develop a QME approach to investigate the STML of molecules coupled to quantized electromagnetic modes. Within our method we include the effects of electronic tunneling, vibrational damping, and environmental dephasing, and can describe both weak and strong plexcitonic coupling regimes. Our approach extends existing methods and includes a full quantum description of the coherent state dynamics.  Our method is valid in both single-particle and many-body representations, allowing future studies to balance computational effort with chemical accuracy.

	Motivated by the observation of HL in the STML of TPP,\cite{dong2010generation} and the argument that it was a consequence of STM-induced Purcell enhancement \cite{chen2015molecular,tian2011electroluminescence,le2007surface}, we calculated the EL of a TPP molecule coupled to a single quantum plasmon mode.  In the weak coupling regime, we recover the experimentally observed spectra, including the above-threshold HL.  Using a fully quantum plasmon theory, we conclude that the low-bias HL peak may be a consequence of the weakly coherent energy exchange dynamics.  Finally, we identify several signatures of the formation of plexcitonic states: a split peak structure of the EL, and the shifted spectral weight as the plasmon resonance is tuned.

	Although the strong coupling regime has not yet been observed in STML systems, molecular systems with coupling strengths of hundreds of meV have been fabricated\cite{bellessa2004strong,dintinger2005strong,hakala2009vacuum,schwartz2011reversible, guebrou2012coherent,fofang2008plexcitonic,vasa2013real,melnikau2016rabi,wang2016role}.  For the TPP system investigated here, we find  that it is physically plausible to achieve the strong coupling regime if the nanocavity losses are reduced slightly, e.g. via careful material selection or patterning of the substrate to reduce the effective plasmon mode volume (See Supporting Information).  If realized, the ability to measure the spatial distribution of the electro-optical response of molecules operating in the strong coupling regime would be invaluable in the development of myriad quantum information applications, and would herald a new phase in the study of QED and molecular dynamics.

\section*{Acknowledgments}
We acknowledge support from the Air Force Office of Scientific Research (Program Manager Dr. Gernot Pomrenke) under contract number 15RYCOR159.

\section*{Author contributions statement}
J.P.B and J.R.H conceived the idea and performed the theoretical analysis together.  J.P.B. developed the quantum transport theory and code.  Both J.P.B. and J.R.H. wrote the manuscript.

\section*{Additional information}
\textbf{Competing financial interests}: The authors declare that they have no competing interests. \\
\textbf{Supplementary Information}: Details of the Franck-Condon calculations, estimates of the Rabi splitting in TPP, and room temperature STML of TPP can be found in the supporting information


text:
\section{Introduction\label{intro}}

The use and manipulation of entanglement is central to the exploitation of quantum computation (see, e.g., [1-8]). The quantum system obviously ``knows'' what its own entanglement is, though extraction of that information is not obvious; thus, we use dynamic learning methods\cite{lecun,werbos} to map this information onto a single experimental measurement which is our entanglement indicator\cite{behrmanqic}. Our method does not require prior state reconstruction or lengthy optimization\cite{vedral, tamaryan}, nor must the system be ``close'' to a given entangled state \cite{filip}.  An entanglement witness emerges from the learning process. We use knowledge of the smaller two-qubit system as a means of ``bootstrapping'' to larger systems \cite{efron}. As the size of the system grows the amount of additional training necessary diminishes\cite{nabic,behrmanmulti}, unlike other methods, e.g., which require knowledge or reconstruction of the density matrix \cite{bennett2, wootters, ghz}; thus, our method potentially may be of general applicability even to large-scale quantum computers, once they are built.

In any experimental implementation, though, we need also to consider that no setup is perfect: there will always be some uncertainty due to extraneous effects.  In quantum systems there is also the problem of decoherence. Classically learning systems such as neural networks have proven fault tolerant and robust to noise; they are also famously used for noise reduction in signals\cite{wasserman}. A machine learning  approach would seem to be an excellent one for issues like noise, decoherence, or missing or damaged data. Here, we show that this is in fact the case, using as a test bed our entanglement indicator on the simple two-qubit system.

\section{Dynamic learning of an entanglement indicator\label{qnn}}

In previous work, we showed we could successfully train a quantum system to estimate its own degree of entanglement, by mapping a measureable output at the final time, to give an indicator of the entanglement of the prepared, initial state. Briefly the method was as follows; for full details the reader is referred to\cite{behrmanqic,nabic,behrmanmulti,singapore,behrmanieee}

We begin with the Schr\"{o}dinger equation for the time evolution of the density matrix $\rho$\cite{peres}:
\begin{equation}
\frac{d \rho}{dt} = \frac{1}{i \hbar}[H, \rho]
\label{schr}
\end{equation}
where  $H$ is the Hamiltonian. We consider an N-qubit system whose Hamiltonian is:
\begin{equation}
H   = \sum_{\alpha=1}^{N}  K_{\alpha} \sigma_{x\alpha} + \varepsilon_{\alpha} \sigma_{z\alpha} + \sum_{\alpha\neq\beta=1}^{N}\zeta_{\alpha\beta} \sigma_{z\alpha} \sigma_{z\beta}
\label{ham}
\end{equation}
where $\{ \sigma \}$ are the Pauli operators corresponding to each of the qubits, $\{K \}$ are the tunneling amplitudes, $\{ \varepsilon \}$  are the biases, and $\{ \zeta \}$, the qubit-qubit couplings. We choose the usual ``charge basis '', in which each qubit's state is given as 0 or 1; for a system of N qubits there are $2^{N}$ states, each labelled by a bit string each of whose numbers corresponds to the state of each qubit, in order. The amplitude for each qubit to tunnel to its opposing state (i.e., switch between the 0 and 1 states) is its $K$ value; each qubit has an external bias represented by its $\varepsilon$ value; and each qubit is coupled to each of the other qubits, with a strength represented by the appropriate $\zeta$ value. Note that, for example, the operator $\sigma_{x1} = \sigma_{x}\otimes I \: ... \: \otimes I$, where there are (N-1) outer products, acts nontrivially only on qubit 1.

The parameter functions $\{ K(t),\varepsilon(t),\zeta(t) \}$ direct the time evolution of the system in the sense that, if one or more of them is changed, the way a given initial state will evolve in time will also change, because of Eqs.~\ref{schr}-\ref{ham}. This is the basis for using our quantum system as a neural network. The role of the input vector is played by the initial density matrix $\rho(0)$, the role of the output by (some function of) the density matrix at the final time, $\rho(t_{f})$, and the role of the ``weights'' of the network  by the parameter functions of the Hamiltonian, $\{ K,\varepsilon,\zeta  \}$, all of which can be adjusted experimentally\cite{yamamoto}. By adjusting these parameters using a machine learning algorithm we can train the system to evolve in time from an input state to a set of particular final states at the final time $t_{f}$. Because the time evolution is quantum mechanical (and, we assume, coherent), a quantum mechanical function, like an entanglement witness of the initial state, can be mapped to an observable of the system's final state, a measurement made at the final time $t_{f}$. Complete details, including a derivation of the quantum dynamic learning paradigm using quantum  backpropagation\cite{lecun} in time\cite{werbos}, are given in \cite{behrmanqic}. We call this quantum system a ``quantum neural network'' (QNN).

We  found \cite{behrmanqic, behrmanieee}  a set of parameter functions that successfully map the input (initial) state of a two-qubit system to a good approximation of the entanglement of formation of that initial state, using as the output the qubit-qubit correlation function at the final time, $\langle \sigma_{z1}(t_{f})\sigma_{z2}(t_{f})\rangle^{2}$. This set of parameter functions was relatively easily generalized to three- four- and five-qubit systems\cite{behrmanmulti}. Here, we will consider the effect of noise on only the simplest case, of two qubits ($N=2$), and for ease of notation we will call the two qubits $A$ and $B$. With continuum\cite{singapore,behrmanieee} rather than piecewise constant parameter functions\cite{behrmanqic,behrmanmulti} training is more rapid and complete; see Figure~\ref{trzero}. The parameter functions are found by training with a set of just four initial quantum states (``inputs''), as shown in Table~\ref{enttraining}: a fully entangled state (``Bell''), a ``Flat'' state (equal amounts of all basis states), a product state``C''  whose initial $(t=0)$ correlation function $\langle \sigma_{zA}\sigma_{zB}\rangle^{2}$ is nonzero, and a partially entangled state ``P''. The parameter functions are shown in Figures~\ref{K0}~\ref{ep0}~\ref{z0}. Because of the symmetry in the Hamiltonian and in the training set, $K_{A}=K_{B}=K$, and $\epsilon_{A}=\epsilon_{B}=\epsilon$. Each is a relatively simple function, well parametrized by only one frequency ($\epsilon, \zeta$) or two ($K$), as shown in Table~\ref{entparamfit} and plotted in Figures ~ref{K0},\ref{ep0},\ref{z0}. The $(input, output)$ pairs are:

\begin{eqnarray}
input = \rho(0)= |\Psi(0)\rangle \langle \Psi(0)| \\ \nonumber
output = \langle \sigma_{zA}(t_{f}) \sigma_{zB}(t_{f}) \rangle ^{2} \rightarrow target
\end{eqnarray}

\noindent with prepared input states at zero time, and corresponding targets, given in Table \ref{enttraining}. This table also shows the QNN indicated entanglement values after training has finished and, for comparison, the entanglement of formation, calculated using the analytic formula \cite{wootters} for comparison. Entanglement of formation is not the only measure of entanglement, of course, but any one that we chose would have qualitatively similar behavior, which we would like our entanglement indicator to imitate. That is, we seek here not exact \underline{agreement} with $E_{F}$ (in which case we would train the state $\frac{1}{\sqrt{3}}(|00\rangle + |01\rangle + |10\rangle)$ to a target value of 0.55), but a robust and internally  self consistent measure, which we would hope would track well with an analytic measure like $E_{F}$. The QNN indicator systematically underestimates $E_{F}$ for partially entangled pure states; this is because we found through simulation that the net naturally trained to the target value of 0.44. See \cite{behrmanqic} for details.

\begin{figure}
\includegraphics[height=2in]{Fig1-eps-converted-to.pdf}
\caption{Total root mean squared error for the training set as a function of epoch (pass through the training set), for the 2-qubit system, with zero noise. Asymptotic error is $1.6 \times 10^{-3}$. For comparison: with piecewise constant functions a similar level of error required 2000 epochs.}
\label{trzero}
\end{figure}

\begin{table}[!t]
\caption{Training data for QNN entanglement witness.}
\label{enttraining}
\begin{tabular}{l l l l}\\
\hline
Input state $|\Psi(0)\rangle$  & Target  & Trained & $E_{F}$ \\
\hline
Bell = $\frac{1}{\sqrt{2}}(|00\rangle + |11\rangle) = |\Phi_{+}\rangle$ & 1.0 & 0.998 & 1.0 \\
Flat = $\frac{1}{2}(|00\rangle + |01\rangle + |10\rangle + |11\rangle)$    & 0.0  &  $1.2 \times 10^{-5}$ & 0.0 \\
C = $\frac{1}{\sqrt{1.25}}(0.5|10\rangle + |11\rangle) $ & 0.0  &  $1.8 \times 10^{-4}$  & 0.0 \\
P = $\frac{1}{\sqrt{3}}(|00\rangle + |01\rangle + |10\rangle)$   & 0.44  & 0.44 & 0.55 \\
\hline
RMS Error   &   { }  & $1.6 \times 10^{-3}$ & { }  \\
Epochs  &  { }   &  100  &  {} \\
\hline \\
\end{tabular}
\end{table}

\begin{figure}
\includegraphics[height=2in]{Fig2-eps-converted-to.pdf}
\caption{Parameter function $K_{A}=K_{B}$ as a function of time (data points), as trained at zero noise for the entanglement indicator, and plotted with the Fourier fit (solid line).}
\label{K0}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig3-eps-converted-to.pdf}
\caption{Parameter function $\epsilon_{A}=\epsilon_{B}$ as a function of time (data points), as trained at zero noise for the entanglement indicator, and plotted with the Fourier fit (solid line).}
\label{ep0}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig4-eps-converted-to.pdf}
\caption{Parameter function $\zeta$ as a function of time (data points), as trained at zero noise for the entanglement indicator, and plotted with the Fourier fit (solid line).}
\label{z0}
\end{figure}

\begin{table}[!t]
\caption{Curvefit coefficients for parameter functions $K$, $\epsilon$, $\zeta$, for QNN entanglement witness.}
\label{entparamfit}
$f(t) =  a_{0} + a_{1}\cos(\omega t) + b_{1}\sin(\omega t) + a_{2}\cos(2\omega t) +  b_{2}\sin(2\omega t) $ \\
\begin{tabular}{l l l l }\\
Coefficient  & K & $\epsilon$ & $\zeta$ \\
\hline
 $a_{0}$  &   $  0.00248 $ & $9.89\times 10^{-5}$ & $9.89 \times 10^{-5}$ \\
 $a_{1}$  &   $3.68\times 10^{-6}$ &$ -4.96\times 10^{-6}$ & $ -1.46\times 10^{-5}$ \\
 $b_{1}$  &   $ 1.95\times 10^{-6}$ & $4.55\times 10^{-6}$  & $ -2.34\times 10^{-5} $ \\
 $a_{2}$  &  $  3.70\times 10^{-6}$ &  ---  &  ---  \\
 $b_{2}$  &  $8.68\times 10^{-7}$ &   ---  &  ---  \\
 $\omega$   &   $  0.0250  $   &$ 0.0250 $ &$ 0.0575 $ \\
\hline
RMS Error   &  $9.77\times 10^{-7} $ & $1.84\times 10^{-6}$ &  $7.29\times 10^{-6}$  \\
\hline \\
\end{tabular}
\end{table}

Once trained, the parameter functions that are found can be tested, by using the Hamiltonian so defined to calculate the QNN indicator for  other initial states. Testing was therefore done on a large number of states not represented in the training set, including fully entangled states, partially entangled states, product (unentangled) states, and even \underline{mixed} states. (Note that only \underline{pure} states were present in the training sets.) The interested reader is referred to our previous work for the (extensive) testing results \cite{behrmanqic, behrmanmulti, behrmanieee}. Note also that the testing was done using the fitted functions only.

\section{ Learning with noise}

Physical systems contain noise, meaning that there is some uncertainty in each of the elements of the density matrix (though, of course, it must remain hermitian and positive semidefinite, with unit trace to conserve probability.) What if the system on which we train is somewhat noisy?

Let us first define our terms. Because we are working with simulations, we can isolate the different effects of ``noise'' and ``decoherence'': here, we will use ``noise'' to refer to random (uncorrelated) magnitudes, of a given size, added to the density matrix elements; and ``decoherence'' to refer to random phases so added. In general, of course, both effects will be present; we will call that ``complex'' noise. Recall that our entanglement indicator is a mapping to a time correlation function, after evolution in time according to the Hamiltonian of Equation \ref{ham}; that is, the entanglement of the initial state is approximated by a measurement performed at the final time. To simulate white noise (zero mean and specified amplitude), random numbers were added at each timestep $\Delta t=0.8$ ns of that time evolution. These numbers have zero time correlation themselves. The level of noise we report is the amplitude, that is, the root-mean-square-average size of these random numbers, that are imposed at each time step. All the simulations were done for the same total time of 317 timesteps, or about 253 ns. Because the system evolves in time, the noise itself propagates, so that by the time the correlation function is measured, at the final time, numbers that seem quite small can build up to destroy a significant amount of entanglement. For comparison with our testing results, we will always therefore include the entanglement of formation for the noisy density matrix, calculated using the Bennett-Wootters formula\cite{bennett2,wootters} and marked ``BW'' on our testing graphs (Figures 12-17, 25-30, and 38-43) below.

We consider first the case of  magnitude noise only. Figure~\ref{magtrn} shows a typical rms error training curve for a fairly large level of noise. As expected, asymptotic error for the training set does increase with increasing noise. The parameter functions also become ``noisy'': See Figures~\ref{kmag},\ref{epmag},\ref{zmag}.

\begin{figure}
\includegraphics[height=2in]{Fig5-eps-converted-to.pdf}
\caption{Total root mean squared error for the training set as a function of epoch (pass through the training set), for the 2-qubit system, with a (magnitude) noise level of 0.014 at each (of 317 total) timestep. Asymptotic error is $3.1 \times 10^{-3}$, about double what it was with no noise.}
\label{magtrn}
\end{figure}

\begin{figure}
\includegraphics[height=2in]{Fig6-eps-converted-to.pdf}
\caption{Parameter function $K_{A}=K_{B}$ as a function of time, as trained at 0.0089 amplitude noise at each of the 317 timesteps, for the entanglement indicator (data points), and plotted with the Fourier fit (solid line). Note the change in scale from Figure~\ref{K0}, because of the (much larger) spread of the noisy data: the Fourier fit is actually almost the same on this graph.}
\label{kmag}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig7-eps-converted-to.pdf}
\caption{Parameter function $\epsilon_{A}=\epsilon_{B}$ as a function of time, as trained at 0.0089 amplitude noise at each of the 317 timesteps,  for the entanglement indicator (data points), and plotted with the Fourier fit (solid line).}
\label{epmag}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig8-eps-converted-to.pdf}
\caption{Parameter function $\zeta$ as a function of time, as trained at 0.0089  noise at each of the 317 timesteps,  for the entanglement indicator (data points), and plotted with the Fourier fit (solid line).}
\label{zmag}
\end{figure}

Much of this variation is not meaningful, though. In fact, the numbers in Table~\ref{entparamfit} are slightly different from the ones we found in \cite{behrmanieee}. Is that difference significant? To investigate this, we tried testing with all but a selected one of the parameter functions' Fourier coefficients set to the trained values, but assigning random numbers (of the right order of magnitude) to the Fourier coefficients of the one chosen. The system was remarkably insensitive to this procedure when the random function was the tunneling function $K$: as long as $K$ was the right order of magnitude, the indicator still tested extremely well - in fact for most of the testing the indicator results were identical. This was not true if $\epsilon$ or $\zeta$ were randomized, however: errors were substantial, particularly when the Fourier frequency $\omega$ is randomized. Still, exact agreement with the trained functions is not, apparently, necessary.

We can, as before, fit each function to a Fourier series, and test using the fitted functions. Figures~\ref{FKmag},\ref{Fepmag},\ref{Fzmag} show the Fourier coefficients for $K_{A}=K_{B}$, $\epsilon_{A}=\epsilon_{B}$ , and $\zeta$, respectively, as functions of increasing noise level. The Fourier components of the tunneling parameter function $K$ are clearly the least sensitive to noise level, while for $\epsilon$ they are a bit more so and for $\zeta$ the most. This is in accordance with the observed insensitivity of the entanglement indicator to $K$: that is, it is true both that the indicator is relatively insensitive to the $K$ function, and that the $K$ function's Fourier fit is insensitive to environmental noise.

\begin{figure}
\includegraphics[height=2in]{Fig9-eps-converted-to.pdf}
\caption{Fourier coefficients for the tunneling parameter functions $K$, as functions of noise level.}
\label{FKmag}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig10-eps-converted-to.pdf}
\caption{Fourier coefficients for the bias parameter functions $\epsilon$, as functions of noise level. }
\label{Fepmag}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig11-eps-converted-to.pdf}
\caption{Fourier coefficients for the coupling parameter function $\zeta$, as functions of noise level.}
\label{Fzmag}
\end{figure}

All of the parameters look fairly stable at these levels of noise, from this point of view. But a more important question is: how much does noise interfere with the net's ability to detect entanglement?

Consider a pure state, specifically, the state $P(\gamma) = \frac{|00 \rangle + |11 \rangle + \gamma |01 \rangle}{\sqrt{2+|\gamma|^{2}}}$. For $\gamma=0$ this is of course the (fully entangled) Bell state; as $\gamma$ increases the entanglement decreases. $P(\gamma)$ is one of the states we tested on in the 2008 paper\cite{behrmanqic}, and, as we found then, the entanglement as computed by the quantum neural network, without noise, tracks the entanglement of formation very well.  How much noise or uncertainty can the network tolerate? We answer this question in two ways. First, we suppose the QNN was trained on the perfect (zero noise) system of the four-pair training set, and we add increasing amounts of noise to test the system for various nonzero values of $\gamma$, and compare the results to the entanglement of formation both at zero noise and at 0.0069 noise. See Figure~\ref{P0mag}. Then, we suppose the QNN was trained on a noisy system, and, again, test with increasing levels of noise. Figures~\ref{P89mag} and \ref{P13mag} show, respectively, an intermediate level and a high level of noise present during training. Recall that the noise should be understood as occurring over the total time evolution interval: that is, an independent (uncorrelated) noise at the given rms level was added at \underline{each} timestep. The two curves showing the entanglement of formation can therefore be thought of as a kind of ``error bars'': the correct entanglement of the system should lie somewhere between the zero noise result and the result at maximum noise, insofar as the QNN tracks well with the entanglement of formation. Presumably the presence of noise does destroy entanglement, but, since the measurement itself is noisy, it is not certain how much is destroyed and how much is simply a bad measurement. Still, it is obvious from the results that, indeed, the QNN technique does an excellent job of remaining robust to pure noise.

\begin{figure}
\includegraphics[height=2in]{Fig12-eps-converted-to.pdf}
\caption{Entanglement of the state $P$ as a function of $\gamma$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero noise (blue) and at 0.0069 noise (orange). In each case the QNN was trained at zero noise, but tested at the given level.}
\label{P0mag}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig13-eps-converted-to.pdf}
\caption{Entanglement of the state $P$ as a function of $\gamma$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero noise (blue) and at 0.0069 noise (orange). In each case the QNN was trained at a noise level of 0.0089, and then tested at the given level.}
\label{P89mag}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig14-eps-converted-to.pdf}
\caption{Entanglement of the state $P$ as a function of $\gamma$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero noise (blue) and at 0.0069 noise (orange). In each case the QNN was trained at 0.013 noise, and then tested at the given level.}
\label{P13mag}
\end{figure}

Second, we consider testing on mixed states. We might expect the QNN to perform significantly less well with these kinds of states, because the training set (Table \ref{enttraining}) contained no mixed states. In fact with zero noise the QNN tested well on several classes of mixed states \cite{behrmanqic}; we need to see if that success is maintained with noisy conditions. Figures~\ref{M0mag},\ref{M89mag},\ref{M13mag} show results for the mixed states $M(\delta) = \frac{\delta|11 \rangle \langle 11| + |\Phi^{+} \rangle \langle \Phi^{+} |}{\delta + 1}$, where $|\Phi^{+} \rangle $ is the Bell state, given in Table \ref{enttraining}, as functions of $\delta$, for QNN trained at zero, 0.0089, and 0.013 noise amplitude, respectively.  Again, the entanglement of formation results can serve as an approximate error bound, and we see that the QNN's entanglement indicator is robust to noise. Indeed, comparison of these figures with the ones for pure states shows that performance on mixed states is even better. Because of this we might expect the QNN indicator to show greater robustness to decoherence than to magnitude noise; we will see in the next section that this is, indeed, the case.

\begin{figure}
\includegraphics[height=2in]{Fig15-eps-converted-to.pdf}
\caption{Entanglement of the state $M$ as a function of $\delta$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero noise (blue) and at 0.0069 noise (orange). In each case the QNN was trained at zero noise, but tested at the given level.}
\label{M0mag}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig16-eps-converted-to.pdf}
\caption{Entanglement of the state $M$ as a function of $\delta$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero noise (blue) and at 0.0069 noise (orange). In each case the QNN was trained at a noise level of 0.0089, and then tested at the given level.}
\label{M89mag}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig17-eps-converted-to.pdf}
\caption{Entanglement of the state $M$ as a function of $\delta$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero noise (blue) and at 0.0069 noise (orange). In each case the QNN was trained at 0.013 noise, and then tested at the given level.}
\label{M13mag}
\end{figure}

\section{ Decoherence}

We now turn to the case of ``pure'' decoherence, that is, we introduce random phases to the elements of the density matrix, without changing their magnitudes. Figure~\ref{phtrn} shows a typical rms error training curve for a fairly large level of phase noise. Again, asymptotic error for the training set does increase with increasing phase noise (decoherence). The parameter functions also become ``noisy'': See Figures~\ref{Kph},\ref{epph},\ref{zph}.

\begin{figure}
\includegraphics[height=2in]{Fig18-eps-converted-to.pdf}
\caption{Total root mean squared error for the training set as a function of epoch (pass through the training set), for the 2-qubit system, with a phase noise level of 0.014 at each (of 317 total) timestep. Asymptotic error is $1.3 \times 10^{-3}$, approximately the same as with no noise.}
\label{phtrn}
\end{figure}

\begin{figure}
\includegraphics[height=2in]{Fig19-eps-converted-to.pdf}
\caption{Parameter function $K_{A}=K_{B}$ as a function of time, as trained at 0.0089 phase noise at each of the 317 timesteps, for the entanglement indicator (data points), and plotted with the Fourier fit (solid line)}
\label{Kph}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig20-eps-converted-to.pdf}
\caption{Parameter function $\epsilon_{A}=\epsilon_{B}$ as a function of time, as trained at 0.0089 phase noise at each of the 317 timesteps,  for the entanglement indicator (data points), and plotted with the Fourier fit (solid line).}
\label{epph}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig21-eps-converted-to.pdf}
\caption{Parameter function $\zeta$ as a function of time, as trained at 0.0089 phase noise at each of the 317 timesteps,  for the entanglement indicator (data points), and plotted with the Fourier fit (solid line).}
\label{zph}
\end{figure}

Figures~\ref{FKph},\ref{Fepph},\ref{Fzph} show the Fourier coefficients for $K_{A}=K_{B}$, $\epsilon_{A}=\epsilon_{B}$ , and $\zeta$, respectively, as functions of increasing decoherence level. It is clear that in the case of decoherence the parameter functions are even more stable than in the previous case of noise.

\begin{figure}
\includegraphics[height=2in]{Fig22-eps-converted-to.pdf}
\caption{Fourier coefficients for the tunneling parameter functions $K$, as functions of decoherence level.}
\label{FKph}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig23-eps-converted-to.pdf}
\caption{Fourier coefficients for the bias parameter functions $\epsilon$, as functions of decoherence level.}
\label{Fepph}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig24-eps-converted-to.pdf}
\caption{Fourier coefficients for the coupling parameter function $\zeta$, as functions of decoherence level.}
\label{Fzph}
\end{figure}

Testing results for the pure state $P$ subject to decoherence are shown in Figures~\ref{P0ph},\ref{P89ph},\ref{P13ph}, as trained, respectively, at zero, 0.0089, and 0.013 decoherence (phase noise), and tested with various levels of phase noise, and compared with the entanglement of formation at zero and at 0.0069 phase noise. The results are extremely good, better even than the ones in Figures~\ref{P0mag},\ref{P89mag},\ref{P13mag}; clearly, the QNN is even better at dealing with decoherence than with ``pure'' noise. Possibly the various (random) phases tend to cancel each other; but since for \underline{definite} phase shifts the QNN underestimates the entanglement, sometimes drastically \cite{behrmanqic,behrmanieee}, this was an unexpectedly good result.

\begin{figure}
\includegraphics[height=2in]{Fig25-eps-converted-to.pdf}
\caption{Entanglement of the state $P$ as a function of $\gamma$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero phase noise (blue) and at 0.0069 phase noise (orange). In each case the QNN was trained at zero phase noise, but tested at the given level.}
\label{P0ph}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig26-eps-converted-to.pdf}
\caption{Entanglement of the state $P$ as a function of $\gamma$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero phase noise (blue) and at 0.0069 phase noise (orange). In each case the QNN was trained at a phase noise level of 0.0089, and then tested at the given level.}
\label{P89ph}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig27-eps-converted-to.pdf}
\caption{Entanglement of the state $P$ as a function of $\gamma$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero phase noise (blue) and at 0.0069 phase noise (orange). In each case the QNN was trained at 0.013 phase noise, and then tested at the given level.}
\label{P13ph}
\end{figure}

Figures~\ref{M0ph},\ref{M89ph},\ref{M13ph} show the performance of the QNN on the mixed state $M$, as trained, respectively, at zero, 0.0089, and 0.013 decoherence (phase noise), and tested with various levels of phase noise, and compared with the entanglement of formation at zero and at 0.0069 phase noise. Again, we see that the QNN entanglement indicator is robust to decoherence.

\begin{figure}
\includegraphics[height=2in]{Fig28-eps-converted-to.pdf}
\caption{Entanglement of the state $M$ as a function of $\delta$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero phase noise (blue) and at 0.0069 phase noise (orange). In each case the QNN was trained at zero noise, but tested at the given level.}
\label{M0ph}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig29-eps-converted-to.pdf}
\caption{Entanglement of the state $M$ as a function of $\delta$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero phase noise (blue) and at 0.0069 phase noise (orange). In each case the QNN was trained at a phase noise level of 0.0089, and then tested at the given level.}
\label{M89ph}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig30-eps-converted-to.pdf}
\caption{Entanglement of the state $M$ as a function of $\delta$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero phase noise (blue) and at 0.0069 phase noise (orange). In each case the QNN was trained at 0.013 phase noise, and then tested at the given level.}
\label{M13ph}
\end{figure}

\section{ Noise plus decoherence}

Finally, we consider the case of noise plus decoherence, that is, what we are calling random complex noise. For this case, we add both magnitude and phase noise. Figure~\ref{cxtrn} shows a typical rms error training curve for a the same level of complex noise as in Figures ~\ref{magtrn} and \ref{phtrn}. Again, asymptotic error for the training set does increase with increasing complex noise. The parameter functions again become ``noisy'': See Figures~\ref{Kcx},\ref{epcx},\ref{zcx}.

\begin{figure}
\includegraphics[height=2in]{Fig31-eps-converted-to.pdf}
\caption{Total root mean squared error for the training set as a function of epoch (pass through the training set), for the 2-qubit system, with a complex noise level of 0.014 at each (of 317 total) timestep. Asymptotic error is $3.0 \times 10^{-3}$, approximately the same as with only magnitude noise.}
\label{cxtrn}
\end{figure}

\begin{figure}
\includegraphics[height=2in]{Fig32-eps-converted-to.pdf}
\caption{Parameter function $K_{A}=K_{B}$ as a function of time, as trained at 0.0089 complex noise at each of the 317 timesteps, for the entanglement indicator (data points), and plotted with the Fourier fit (solid line). Note the change in scale from Figure~\ref{K0}, because of the (much larger) spread of the noisy data: the Fourier fit is actually almost the same on this graph.}
\label{Kcx}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig33-eps-converted-to.pdf}
\caption{Parameter function $\epsilon_{A}=\epsilon_{B}$ as a function of time, as trained at 0.0089  complex noise at each of the 317 timesteps,  for the entanglement indicator (data points), and plotted with the Fourier fit (solid line).}
\label{epcx}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig34-eps-converted-to.pdf}
\caption{Parameter function $\zeta$ as a function of time, as trained at 0.0089 complex noise at each of the 317 timesteps,  for the entanglement indicator (data points), and plotted with the Fourier fit (solid line).}
\label{zcx}
\end{figure}

Again, we test to see how much the Fourier fit changes, this time with complex noise. Figures~\ref{FKcx},\ref{Fepcx},\ref{Fzcx} show the Fourier coefficients as a function of complex noise level, for $K$, $\epsilon$, and $\zeta$, respectively. We can see that the indicator is, again, relatively stable.

\begin{figure}
\includegraphics[height=2in]{Fig35-eps-converted-to.pdf}
\caption{Fourier coefficients for the tunneling parameter functions $K$, as functions of complex noise level.}
\label{FKcx}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig36-eps-converted-to.pdf}
\caption{Fourier coefficients for the bias parameter functions $\epsilon$, as functions of complex noise level.}
\label{Fepcx}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig37-eps-converted-to.pdf}
\caption{Fourier coefficients for the coupling parameter function $\zeta$, as functions of complex noise level.}
\label{Fzcx}
\end{figure}

Again we use the Fourier fitted functions to test on both pure and mixed states. Figures~\ref{P0cx},\ref{P89cx},\ref{P13cx} show performance of the QNN on the entanglement of the state $P$, as trained, respectively, at zero, 0.0089, and 0.013 amplitude complex noise, and tested with various levels of complex noise, and compared with the entanglement of formation at zero and at 0.0069 complex noise. Figures~\ref{M0cx}~\ref{M89cx}~\ref{M13cx} show the performance of the QNN on the mixed state $M$, as trained, respectively, at zero, 0.0089, and 0.013 complex noise, and tested with various levels of complex noise, and compared with the entanglement of formation at zero and at 0.0069 complex noise. Results are excellent; in fact, they are somewhat better than for the case of only ``magnitude'' noise. In some sense allowing for decoherence makes the indicator even more robust.

\begin{figure}
\includegraphics[height=2in]{Fig38-eps-converted-to.pdf}
\caption{Entanglement of the state $P$ as a function of $\gamma$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero noise (blue) and at 0.69\% noise plus decoherence (orange). In each case the QNN was trained at zero noise, but tested at the given level.}
\label{P0cx}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig39-eps-converted-to.pdf}
\caption{Entanglement of the state $P$ as a function of $\gamma$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero noise (blue) and at 0.0069 noise plus decoherence (orange). In each case the QNN was trained at a noise level of 0.0089 complex noise, and then tested at the given level.}
\label{P89cx}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig40-eps-converted-to.pdf}
\caption{Entanglement of the state $P$ as a function of $\gamma$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero complex noise (blue) and at 0.0069 noise plus decoherence(orange). In each case the QNN was trained at 0.013 level complex noise, and then tested at the given level.}
\label{P13cx}
\end{figure}

\begin{figure}
\includegraphics[height=2in]{Fig41-eps-converted-to.pdf}
\caption{Entanglement of the state $M$ as a function of $\delta$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero noise (blue) and at 0.0069 noise plus decoherence(orange). In each case the QNN was trained at zero noise, but tested at the given level.}
\label{M0cx}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig42-eps-converted-to.pdf}
\caption{Entanglement of the state $M$ as a function of $\delta$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero noise (blue) and at 0.0069 noise plus decoherence (orange). In each case the QNN was trained at a complex noise level of 0.0089, and then tested at the given level.}
\label{M89cx}
\end{figure}
\begin{figure}
\includegraphics[height=2in]{Fig43-eps-converted-to.pdf}
\caption{Entanglement of the state $M$ as a function of $\delta$, as calculated by the QNN, and compared with the entanglement of formation (marked ``BW'') at zero noise (blue) and at 0.0069 noise plus decoherence (orange). In each case the QNN was trained at 0.013 complex noise, and then tested at the given level.}
\label{M13cx}
\end{figure}

\section{ Conclusions}

In previous work, we have proposed an entanglement indicator for a general qubit system. This indicator is a quantum system that processes the state whose entanglement is to be estimated. The parameters of the quantum system are adjusted via a supervised learning process using a sparse training set of states whose entanglement is well-defined. The learning is continued until sufficient training is achieved.  The trained parameter functions are well represented by single frequency functions (first order Fourier curve fit.)

We have shown here that those functions are robust to fairly high levels of noise and decoherence. The quantum neural network tests well on ``unknown'' states, both pure and mixed. Performance on decoherence in particular was excellent. We are reasonably confident that our results show that quantum neural networks are well suited for dealing with these types of problems in quantum computing.

We are currently working to extend our results on noise and decoherence to multiple-qubit systems, using the well-known neural network technique of bootstrapping \cite{efron}. We also wish to understand exactly why the QNN is so robust to noise and decoherence. Classical neural networks are robust to noise and single neuron/synapse failure because of the multiple-redundance of parallel computing. Here, we have only a very small number of qubits/neurons, but we have designed our quantum network as operating over propagation in time, which can be written as a superposition of a very large number of definite time paths, using the Feynman path integral representation of quantum mechanics \cite{feynman}. In this picture, the instantaneous states of the quantum system at intermediate times, which are integrated over, play the role of ``virtual neurons''\cite{behrman}. In other words, it is possible that quantum superposition ensures redundance, even when the physical number of qubits is small, and thereby supplies fault-tolerance.

\section{Acknowledgements}
\noindent
This work was supported in part by the National Science Foundation under Grant No. NSF PHY05-51164, through the KITP Scholars program (ECB), at the Kavli Institute for Theoretical Physics, University of California at Santa Barbara, Santa Barbara, CA.

redpajama_set_name: RedPajamaArXiv