# Show $\lim\limits_{t\to\infty}\Bigg|\sum\limits_{n=0}^\infty\Theta(t-nR)\frac{(\Gamma(t-nR))^n}{n!}e^{-\Gamma(t-nR)}\Bigg|^2=\frac1{(1+\Gamma R)^2}$

I have encountered the following problem while studying non-Markovian effects in real-time dynamics of open quantum systems.

In particular, I was studying a system comprised of two qubits (qubit is a standard shorthand for two level quantum system) separated by distance (in configuration space e.g. a laboratory) $$R$$ from one another and coupled to a one-dimensional Bosonic reservoir hosting a pair of boson species corresponding to right and left moving photon fields with linear dispersion relation and fixed propagation speed (equal to $$1$$ in my units (Planck's constant is set to be $$2\pi$$)). It is really one of the simplest system one can imagine which possesses the so-called delayed coherent quantum feedback, a property that dynamics of the quantum system at any time $$T$$ depends on the entire history of its dynamics for all times $$t\in[0, T]$$ in deterministic (i.e. classical, not quantum, and thus very human-being controllable) fashion.

You can easily guess that without any knowledge of physics whatsoever, you have quantum information that propagates with finite velocity $$v=1$$ over given distance $$R$$, so there is an intrinsic time delay $$\tau=vR$$ in the system and due to Lorentz covariance this cannot be altered by any quantum effects. The Hamiltonian operator of such a system can be written as $$H=H_{0}+V$$, where \begin{align}H_{0}&=\sum_{n=1, 2}\frac{\Delta_{n}}{2}\sigma_{3}^{(n)}+\sum_{\mu=1, 2}\int_{-\infty}^{\infty}dk\omega_{\mu}(k)a^{\dagger}_{\mu}(k)a_{\mu}(k),\\V&=\sum_{n=1, 2}\sum_{\mu=1, 2}\int_{-\infty}^{\infty}dkg_{\mu, n}(k)a^{\dagger}_{\mu}(k)\sigma_{+}^{(n)}+g_{\mu, n}^{*}(k)\sigma_{-}^{(n)}a_{\mu}(k).\end{align} Here $$\Delta_{n}$$ is the detuning of the qubit number $$n$$ from some reference energy $$k_{0}$$, i.e. $$\Delta_{n}=\Omega_{n}-k_{0}$$, where $$\Omega_{n}$$ is the transition frequency of the qubit number $$n$$ and $$k_{0}$$ is the momentum around which the spectrum of the bath of boson fields was linearised (we of course expect photons of energy ($$hvk_{0}/(2\pi)=k_{0}$$) $$k_{0}$$ to couple most strongly to some system with transition frequency $$\min[\Omega_{1}, \ \Omega_{2}]$$). Further, $$\omega_{\mu}(k)=hvk/(2\pi)=k$$ is the energy of the photon of flavour $$\mu$$ and momentum $$k$$. Such a photon is created by operator $$a^{\dagger}_{\mu}(k)$$ and destroyed by $$a_{\mu}(k)$$, these obey non-zero commutation $$[a_{\mu}(k), a^{\dagger}_{\mu'}(k')]=\delta(k-k')\delta_{\mu, \mu'}$$. $$\sigma_{j}^{(1)}=\sigma_{j}\otimes\sigma_{0}, \ \sigma_{j}^{(2)}=\sigma_{0}\otimes\sigma_{j}$$ , where $$\sigma_{\pm}=(\sigma_{1}+i\sigma_{2})/2$$, and $$\sigma_{1,2, 3}$$ are usual Pauli-matrices ($$\sigma_{0}$$ is an identity on $$\mathbb{C}^{2}$$). The coupling constants are defined as $$g_{\mu, n}(k)=\sqrt{\frac{\Gamma_{n}}{2\pi}}e^{-ic_{\mu}c_{n}(k_{0}+k)R/2},$$ where $$\Gamma_{n}$$ is the bare decay rate of a single qubit into the continuum, and $$c_{s}=(-)^{s+1}, \ s=\mu, \ n$$ distinguishes the coupling to right/left (with index $$\mu$$) photons of an atom number $$n=1, 2$$, located at $$\pm R/2$$.

As a mock up problem I consider the following. Consider the system described above Hamiltonian prepared at time $$t_{0}=0$$ in the following 2-parameter family of states
$$|\psi(0)\rangle=(\cos\vartheta|1\rangle\otimes|0\rangle+e^{i\varphi/2}\sin\vartheta|0\rangle\otimes|1\rangle)\otimes|\Omega\rangle.$$ Here the quantum states are ordered as qubit 1, qubit 2, bath, i.e. $$|a\rangle\otimes|b\rangle\otimes|c\rangle$$ = qubit one is in state $$a$$, qubit $$2$$ is in $$b$$ and bosons are in $$c$$. Here $$|\Omega\rangle$$ is the "vacuum" state of bosons defined by $$a_{\mu}(k)|\Omega\rangle=0,\forall \mu=1,2 , \ k\in\mathbb{R}$$. In physics you'll call such a setup a spontaneous emission problem. I was able to deduce with the help of diagrammatic techniques (see our recent paper https://arxiv.org/pdf/2101.07603.pdf related to the $$T\rightarrow\infty$$ limit of non-Markovian systems) that the exact survival probability amplitude for an initial state defined above has the form (this result is exact in all parameter regime iff $$\Theta(t)$$ the Heaviside step function is defined equal to $$1$$ at $$t=0$$, it is indeed ambiguous since we use use both Hille–Yosida theorem and Sokhotski-Plemelj theorem which are good until $$T=t_{0}=0$$ where divergencies happen, then you have to maintain the order of limits with some care) $$P(t)=|a(t)|^{2},\quad a(t)=\oint_{C_{+\eta}}\langle{\psi(0)|G(z)|\psi(0)\rangle}e^{-izt}\frac{dz}{2\pi i}, \ t>0.$$ Here the integration contour $$C_{+\eta}$$ is a Bromiwich style contour suspending itself above the real axis with positive imaginary part $$\eta\searrow0$$. $$G(z)$$ is the operator valued function known in physics as the retarded Green's function. One thus sees that all important information about dynamics is contained in the poles and branch cuts of $$G(z)$$. The projection of the retarded Green's function onto the single atomic excitation subspace can be determined analytically in the exact form: $$G^{(1)}(z)=\frac{1}{1+g_{1}(z)g_{2}(z)\Gamma_{1}\Gamma_{2}e^{2i(z+k_{0})R}}(g_{1}(zc|10\rangle\langle{10}|+g_{2}(z)|01\rangle\langle{01}|$$ $$-ig_{1}(z)g_{2}(z)\sqrt{\Gamma_{1}\Gamma_{2}}e^{i(z+k_{0})R}(\sigma_{+}^{(1)}\sigma_{-}^{(2)}+\sigma_{-}^{(2)}\sigma_{+}^{(1)})$$ $$=G_{e1}(z)|10\rangle\langle{10}|+G_{e2}(z)|01\rangle\langle{01}|+G_{o}(\sigma_{+}^{(1)}\sigma_{-}^{(2)}+\sigma_{-}^{(2)}\sigma_{+}^{(1)}),$$ where \begin{align} G_{ej}(z)=\frac{g_{j}(z)}{1+g_{1}(z)g_{2}(z)\Gamma_{1}\Gamma_{2}e^{2i(z+k_{0})R}}, \quad G_{o}(z)=\frac{-ig_{1}(z)g_{2}(z)\sqrt{\Gamma_{1}\Gamma_{2}}e^{i(z+k_{0})R}}{1+g_{1}(z)g_{2}(z)\Gamma_{1}\Gamma_{2}e^{2i(z+k_{0})R}}, \end{align} and $$g_{j}(z)=(z-\Delta_{j}+i\Gamma_{j})^{-1}$$ are the single-quit green's functions (note that their Laplace transform is trivial $$e^{-i\Delta_{j}t}e^{-\Gamma_{j}t}$$ due to locality (like you've learned in high school-an exponential decay)). The complete answer is a bit involved but can be clearly expressed in terms of Fourier images of $$G_{ej}(z), \ G_{o}(z)$$, these are

\begin{align} G_{e1}(t)&=\oint_{C_{+\eta}}\frac{dz}{2\pi i}\frac{g_{1}(z)}{1+g_{1}(z)g_{2}(z)\Gamma_{1}\Gamma_{2}e^{2i(z+k_{0})R}}e^{-izt}\\ &=e^{-i\Delta_{1}t}e^{-\Gamma_{1}t}+\sum_{m=1}^{\infty}(-\Gamma_{1}\Gamma_{2}e^{2ik_{0}R})^{m}I(2m, m+1, m, t),\\ G_{e2}(t)&=\oint_{C_{+\eta}}\frac{dz}{2\pi i}\frac{g_{2}(z)}{1+g_{1}(z)g_{2}(z)\Gamma_{1}\Gamma_{2}e^{2i(z+k_{0})R}}e^{-izt}\\ &=e^{-i\Delta_{2}t}e^{-\Gamma_{2}t}+\sum_{m=1}^{\infty}(-\Gamma_{1}\Gamma_{2}e^{2ik_{0}R})^{m}I(2m, m, m+1, t),\\ G_{o}(t)&=\oint_{C_{+\eta}}\frac{dz}{2\pi i}\frac{-ig_{1}(z)g_{2}(z)\sqrt{\Gamma_{1}\Gamma_{2}}e^{i(z+k_{0})R}}{1+g_{1}(z)g_{2}(z)\Gamma_{1}\Gamma_{2}e^{2i(z+k_{0})R}}e^{-izt}\\ =&-i\sqrt{\Gamma_{1}\Gamma_{2}}e^{ik_{0}R}\sum_{m=0}^{\infty}(-\Gamma_{1}\Gamma_{2}e^{2ik_{0}R})^{m}I(2m+1, m+1, m+1, t) \end{align} where \begin{align} &I(a, b, c, t)=\oint_{C_{+\eta}}\frac{dz}{2\pi i}e^{-iz(t-aR)}g_{1}^{b}(z)g_{2}^{c}(z) =\int_{C_{+}}\frac{dz}{2\pi{i}}\frac{e^{-iz(t-aR)}}{(z-\Delta_{1}+i\Gamma_{1})^{b}(z-\Delta_{2}+i\Gamma_{2})^{c}}\\ &=\Theta(t-aR)\Bigg(\sum_{k=0}^{b-1}\frac{(-1)^{k}(c+k-1)!(-i(t-aR))^{b-k-1}}{k!(b-k-1)!(c-1)!}\frac{e^{-i\Delta_{1}(t-aR)}e^{-\Gamma_{1}(t-aR)}}{(\Delta_{2}-\Delta_{1}+i(\Gamma_{2}-\Gamma_{1}))^{c+k}}\\ &+\sum_{k=0}^{c-1}\frac{(-1)^{k}(b+k-1)!(-i(t-aR))^{c-k-1}}{k!(c-k-1)!(b-1)!}\frac{e^{-i\Delta_{2}(t-aR)}e^{-\Gamma_{2}(t-aR)}}{(\Delta_{1}-\Delta_{2}+i(\Gamma_{1}-\Gamma_{2}))^{b+k}}\Bigg). \end{align} So far so good. These functions are looking physically correct, etc. You can see an infinite number of revival peaks in survival probability which happen on all integer multiples of delay time $$vR$$, i.e. when the photon emitted by one atom reaches the other, isn't it cool you can control what atom does by just moving it around?=)

Many physicists though are recently discussing the possibility of so called DARK states, the states for which $$p(t)=1, \ \forall t>0$$. Their argument is based on Markov approximation though, it is the limit where $$\max{\Gamma_{n}}_{n=1, \ 2}\times vR\ll1$$. The first assumption to be done is that the quits are equivalent and "bright" (no detuning) $$\Delta_{1}=\Delta_{2}=0, \ \Gamma_{n}=\Gamma$$. The second is $$\theta=\pi/4, \mod 2\pi, \varphi=4\pi, \mod 2\pi$$ for bright and $$\theta=\pi/4, \mod 2\pi, \varphi=2\pi, \mod 2\pi$$ for dark states. My analysis shows \begin{align} \label{eq: DE58} p_{B}(t)=&\Bigg{|}G_{e}(t)+G_{o}(t)\Bigg{|}^{2}=\Bigg{|}\sum_{n=0}^{\infty}\Theta(t-nR)\frac{(-\Gamma(t-nR))^{n}}{n!}e^{-\Gamma(t-nR)}\Bigg{|}^{2},\\ \label{eq: DE59} p_{D}(t)=&\Bigg{|}G_{e}(t)-G_{o}(t)\Bigg{|}^{2}=\Bigg{|}\sum_{n=0}^{\infty}\Theta(t-nR)\frac{(\Gamma(t-nR))^{n}}{n!}e^{-\Gamma(t-nR)}\Bigg{|}^{2}. \end{align} I.e. for the dark state we obtain $$p_{D}(t)=1, \forall t>0$$ at $$R=0$$. For $$R\neq0$$ when $$t\rightarrow \infty$$ the sereies converges to a finite value $$p(t)\rightarrow\text{constant}.$$ By numerical exercise I'd like to claim that the $$t$$ infinity limit of $$p_{D}(t)$$ is precisely $$1/(1+\Gamma R)^{2}$$. Can you help me to prove or disprove that? And if you happen to know the closed form the above series do not hesitate to share.

By the long-time limit I understand you take above series with all $$\Theta$$ equal to $$1$$ (do not bother about how rigorous this is please). To be specific: if you happen to know something about series like $$x^{n}(c-nx)^n/n!$$ give us a clue.

Another clue from numerics: at super large times $$\Gamma t\gg1$$ we deduce that the quantum decay of bright state follows a sub-exponential power law $$1/t$$ in the long time regime.

Numerics show the $$R/(2\pi t)$$ law, see picture.

• "Since I've noticed that questions coming from theoretical physics on MSE attract more attention when some context is exposed" — I think you may have gone a little overboard. It takes a little effort to find your actual question. Commented Jul 28, 2021 at 23:48
• @Jacob Manaker true, sorry for that! Commented Jul 28, 2021 at 23:53
• Is $\Gamma$ a constant in your series?
– Gary
Commented Jul 29, 2021 at 11:56
• This gamma is actually the dipole moment transition element if we consider a real atom from periodic table, but may be anything really, since experimentalists nowadays can mock up that situation in 1000 different setup Commented Jul 29, 2021 at 12:18
• @TheSimpliFire can you recommend any desmos handbook? It seems to about to become a big thing today! Commented Jul 29, 2021 at 14:57

We begin by proving a simple lemma that will aid the limit calculations.

Lemma 1. For all $$S>0$$, the real parts of the roots of $$z+1-e^{-zS}$$ over $$\Bbb C$$ are non-positive.

Proof: Evidently one root is $$z=0$$. Now let $$z=u+iv$$ so that $$z+1-e^{-zS}=0$$ yields the equations $$(u+1)e^{uS}=\cos vS$$ and $$ve^{uS}=-\sin vS$$. If $$u>0$$, then $$(u+1)e^{uS}>1$$ which exceeds the bounds of $$\cos vS$$, so $$u$$ must be negative. $$\quad\square$$

As mentioned in the post, we can drop the Heaviside step function as it has no effect when $$t$$ is positive. Thus, the conjecture that $$p_D(+\infty)=1/(1+\Gamma R)^2$$ is equivalent to the following.

Proposition 2. Let $$\Gamma,R$$ be positive constants. Then $$\displaystyle\lim\limits_{t\to+\infty}\sum_{n=0}^\infty\frac{(\Gamma(t-nR))^n}{n!e^{\Gamma(t-nR)}}=\frac1{1+\Gamma R}$$.

Proof: Let $$S=\Gamma R$$ and $$x=\Gamma t$$ so that the statement becomes $$\displaystyle\lim\limits_{x\to+\infty}\sum_{n=0}^\infty\frac{(x-Sn)^n}{n!e^{x-Sn}}=\frac1{1+S}$$. We now apply a Feynman-like trick on the sum to convert the summand into an integral. Since $$(x-Sn)^ne^{Sn-x}$$ is the $$n$$th derivative of $$e^{z(x-Sn)}$$ at $$z=-1$$, we obtain $$\sum_{n=0}^\infty\frac{(x-Sn)^n}{n!e^{x-Sn}}=\sum_{n=0}^\infty\frac{D_z^ne^{z(x-Sn)}\mid_{z=-1}}{n!}=\sum\operatorname{Res}\left(\sum_{n=0}^\infty\frac{e^{z(x-Sn)}}{(z+1)^{n+1}}\right)$$ by Cauchy's integral formula. The inner sum is a standard geometric progression given by $$\frac{e^{zx}}{z+1}\cdot\frac1{1-e^{-zS}/(z+1)}=\frac{e^{zx}}{z+1-e^{-zS}}$$ so we seek singularities of the denominator, which occur when $$z=0$$ or $$z=W_k(Se^S)/S-1$$ where $$W_k$$ denotes the $$k$$th branch of the Lambert-$$W$$ function and $$k\ne0$$ is an integer. By taking any circular contour that includes $$z=-1,0$$, we have \begin{align}\sum_{n=0}^\infty\frac{(x-Sn)^n}{n!e^{x-Sn}}&=\operatorname{Res}\left(\frac{e^{zx}}{z+1-e^{-zS}},0\right)+\sum_{0<|k|<\ell}\operatorname{Res}\left(\frac{e^{zx}}{z+1-e^{-zS}},\frac{W_k(Se^S)}S-1\right)\\&=\frac1{1+S}+\sum_{0<|k|<\ell}\left[e^{(W_k(Se^S)/S-1)x}\lim\limits_{z\to W_k(Se^S)/S-1}\frac{z+1-W_k(Se^S)/S}{z+1-e^{-zS}}\right]\end{align} where $$\ell$$ is a finite integer and the limit is a finite constant dependent on $$k,S$$ — note WolframAlpha evaluates it as $$1/(1+W_k(Se^S))$$. From Lemma 1, the real part of $$W_k(Se^S)/S-1$$ must be negative, so everything in the sum converges to zero as $$x\to+\infty$$. $$\quad\square$$

Remark 3. An argument akin to that of Proposition 2 can be used to show that $$p_B(+\infty)=0$$ where we replace the derivative at $$z=-1$$ by the derivative at $$z=1$$.

• Is there a preferred way to quote you? Commented Jul 29, 2021 at 16:24
• We actually new about the solution of pole equation in terms of Lambert function, but did not have any idea about it and thus expanded integral in the series as you can see. You have basically done opposite) Commented Jul 29, 2021 at 17:03
• @KirylPesotski if you are going to use the answer for a paper that you are preparing to submit to a journal, you can click the "cite" to cite this answer Commented Jul 29, 2021 at 19:27
• @KirylPesotski I now have a complete proof. Commented Jul 29, 2021 at 19:38
• @TheSimpliFire I'm very excited by this proof, similar series now appear 2 times a week in Physical Review A. About $10$ years ago nobody used to treat the finite $R$ regime seriously since it was not possible to experimentally assess any finite size effect. All chips in waveguide QED are very tiny, and the speed of propagation of information is so large that $R=0$ approximation is very good. Now experiments allow to actually study dynamics in real time, so this thing became ultra hot! Commented Jul 30, 2021 at 12:43