This wikilog article is a draft, it was not published yet.

Statistics of multiphoton emission

⇠ Back to Blog:Science

Our last paper with Guillermo Díaz Camacho et al. on the dynamics of emission of multiphotons, in particular, of photon bundles (?!) involves a lot of complicated combinatoric formulas, such as the average time of detection for the $k$th photon from a $N$-photon Fock state that is spontaneously emitted with a radiative rate $\gamma_a$ and measured by a detector with bandwidth $\Gamma$:

$$\begin{multline} \langle t_{k}^{(N)} \rangle = 2 N! \left(\frac{\Gamma_+ \Gamma \gamma_a}{\Gamma_-^2}\right)^N \sum_{k_1 + k_2 + k_3\atop = N-k}\sum_{k_4 + \cdots+ k_9\atop = k-1} \sum_{k_{10} + k_{11}\atop = 2}\prod_{j=1}^{11}{1\over k_j!}\times\\ \frac{(-1)^{k_3 + k_6 + k_7 + k_8 + k_{11}}}{\gamma_a^{k_1 +k_4+k_7} \Gamma^{k_2+k_5+k_8} \left( \frac{\Gamma_+}{4} \right)^{k_3 + k_6 + k_9} \big( \gamma_a (k_1 + k_7+{k_{11}\over2}) +\Gamma( k_2 +k_8 +{k_{10}\over2}) + \frac{\Gamma_+}{2} (k_3+k_9)\big)^2 }\,. \end{multline}$$

Here I provide a Mathematica code to turn this general result into particular and limiting cases, and discuss briefly the whole thing. The formula is obtained by integration over the joint probability density function $\phi_\Gamma$ for emitting the photons at times $t_i$, which we show in this paper is given by:

$$\phi_\Gamma(t_1,\dots,t_N)=(-1)^{N}N!\gamma_a^N\left(\frac{\Gamma}{\Gamma_-}\right)^{2N }\prod_{i=1}^{N}\left[{d\over dt_i}({e^{-\gamma_a t_i}\over\gamma_a}+{e^{-\Gamma t_i}\over\Gamma}-4{e^{-\Gamma_+ t_i/2}\over\Gamma_+})\right]{\mathbf{1}}_{[t_{i-1}, t_{i+1}[}(t_{i})$$ where $\Gamma_\pm\equiv\Gamma\pm\gamma_a$. When tracing out all the photons to keep only the $k$th of interest, through some clever math hackery, we arrive to the computation of $\int_0^\infty\left({e^{-\gamma_a t}\over\gamma_a}+{e^{-\Gamma t}\over\Gamma}-4{e^{-\Gamma_+ t/2}\over\Gamma_+}\right)^N\,dt$ which is easy to do with the multinomial theorem, which is where all the mad coefficients come from. The nice thing, of course, is that it covers all cases with a single formula, i.e., it can tell you about the 1st photon from a biphoton or the 7th photon from a 15-photon bundle. Now the direct calculation of these particular cases is a bit tedious. We can ask a formal computer language, like Mathematica, to do it. This is the generic expression to be summed over:

tkN[k1_, k2_, k3_, k4_, k5_, k6_, k7_, k8_, k9_, k10_, k11_] := 
 2 \[CapitalNu]! (((\[CapitalGamma] + \[Gamma]a) \[CapitalGamma] \
\[Gamma]a)/(\[CapitalGamma] - \[Gamma]a)^2)^\[CapitalNu] 1/(
  k1! k2! k3! k4! k5! k6! k7! k8! k9! k10! k11!) (-1)^(
  k3 + k6 + k7 + k8 + k11)/(\[Gamma]a^(k1 + k4 + k7) \[CapitalGamma]^(
   k2 + k5 + k8) ((\[CapitalGamma] + \[Gamma]a)/4)^(
   k3 + k6 + k9) (\[Gamma]a (k1 + k7 + k11/2) + \[CapitalGamma] (k2 + k8 + k10/
        2) + (\[CapitalGamma] + \[Gamma]a)/2 (k3 + k9))^2)

and defining the following module to provide all the ways to write an integer $s$ as the sum of exactly $k$ terms (possibly zero):

PartitionForMultinomialSum[n_, k_] := Module[{},
  Flatten[
   Permutations /@ (PadRight[#, k] & /@ 
      Select[IntegerPartitions[n], Length[#] <= k &]), 1]]

for instance,

In:= PartitionForMultinomialSum[3, 2]
 
Out= {{3, 0}, {0, 3}, {2, 1}, {1, 2}}

we can now compute, e.g., the average time of detection for the k=2nd photon from a N=3-photon bundle as:

In:= k = 2; \[CapitalNu] = 3;
In:= listkN = Flatten[Table[Flatten[{lk1, lk2, lk3}],
   {lk1, PartitionForMultinomialSum[\[CapitalNu] - k, 3]},
   {lk2, PartitionForMultinomialSum[k - 1, 6]},
   {lk3, PartitionForMultinomialSum[2, 2]}
   ], 3 - 1]
Out= {{1, 0, 0, 1, 0, 0, 0, 0, 0, 2, 0}, {1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 
  2}, {1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1}, {1, 0, 0, 0, 1, 0, 0, 0, 0, 
  2, 0}, {1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 2}, {1, 0, 0, 0, 1, 0, 0, 0, 
  0, 1, 1}, {1, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0}, {1, 0, 0, 0, 0, 1, 0, 
  0, 0, 0, 2}, {1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1}, {1, 0, 0, 0, 0, 0, 
  1, 0, 0, 2, 0}, {1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2}, {1, 0, 0, 0, 0, 
  0, 1, 0, 0, 1, 1}, {1, 0, 0, 0, 0, 0, 0, 1, 0, 2, 0}, {1, 0, 0, 0, 
  0, 0, 0, 1, 0, 0, 2}, {1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1}, {1, 0, 0, 
  0, 0, 0, 0, 0, 1, 2, 0}, {1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 2}, {1, 0, 
  0, 0, 0, 0, 0, 0, 1, 1, 1}, {0, 1, 0, 1, 0, 0, 0, 0, 0, 2, 0}, {0, 
  1, 0, 1, 0, 0, 0, 0, 0, 0, 2}, {0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 
  1}, {0, 1, 0, 0, 1, 0, 0, 0, 0, 2, 0}, {0, 1, 0, 0, 1, 0, 0, 0, 0, 
  0, 2}, {0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1}, {0, 1, 0, 0, 0, 1, 0, 0, 
  0, 2, 0}, {0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 2}, {0, 1, 0, 0, 0, 1, 0, 
  0, 0, 1, 1}, {0, 1, 0, 0, 0, 0, 1, 0, 0, 2, 0}, {0, 1, 0, 0, 0, 0, 
  1, 0, 0, 0, 2}, {0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1}, {0, 1, 0, 0, 0, 
  0, 0, 1, 0, 2, 0}, {0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 2}, {0, 1, 0, 0, 
  0, 0, 0, 1, 0, 1, 1}, {0, 1, 0, 0, 0, 0, 0, 0, 1, 2, 0}, {0, 1, 0, 
  0, 0, 0, 0, 0, 1, 0, 2}, {0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1}, {0, 0, 
  1, 1, 0, 0, 0, 0, 0, 2, 0}, {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2}, {0, 
  0, 1, 1, 0, 0, 0, 0, 0, 1, 1}, {0, 0, 1, 0, 1, 0, 0, 0, 0, 2, 
  0}, {0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 2}, {0, 0, 1, 0, 1, 0, 0, 0, 0, 
  1, 1}, {0, 0, 1, 0, 0, 1, 0, 0, 0, 2, 0}, {0, 0, 1, 0, 0, 1, 0, 0, 
  0, 0, 2}, {0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1}, {0, 0, 1, 0, 0, 0, 1, 
  0, 0, 2, 0}, {0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2}, {0, 0, 1, 0, 0, 0, 
  1, 0, 0, 1, 1}, {0, 0, 1, 0, 0, 0, 0, 1, 0, 2, 0}, {0, 0, 1, 0, 0, 
  0, 0, 1, 0, 0, 2}, {0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1}, {0, 0, 1, 0, 
  0, 0, 0, 0, 1, 2, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 2}, {0, 0, 1, 
  0, 0, 0, 0, 0, 1, 1, 1}}
 
In=Simplify[Total[
  tkN[#[[1]], #[[2]], #[[3]], #[[4]], #[[5]], #[[6]], #[[7]], #[[
      8]], #[[9]], #[[10]], #[[11]]] & /@ listkN]]
Out[200]= (150 \[CapitalGamma]^8 + 2345 \[CapitalGamma]^7 \[Gamma]a + 
 14493 \[CapitalGamma]^6 \[Gamma]a^2 + 
 41371 \[CapitalGamma]^5 \[Gamma]a^3 + 
 58786 \[CapitalGamma]^4 \[Gamma]a^4 + 
 41371 \[CapitalGamma]^3 \[Gamma]a^5 + 
 14493 \[CapitalGamma]^2 \[Gamma]a^6 + 
 2345 \[CapitalGamma] \[Gamma]a^7 + 
 150 \[Gamma]a^8)/(6 \[CapitalGamma] \[Gamma]a (\[CapitalGamma] + \
\[Gamma]a) (2 \[CapitalGamma] + \[Gamma]a) (3 \[CapitalGamma] + \
\[Gamma]a) (5 \[CapitalGamma] + \[Gamma]a) (\[CapitalGamma] + 
   2 \[Gamma]a) (\[CapitalGamma] + 3 \[Gamma]a) (\[CapitalGamma] + 
   5 \[Gamma]a))

$$\langle t_2^{(3)} \rangle = \frac{\Big(150 \gamma_a^8 + 2345 \gamma_a^7 \Gamma +14493 \gamma_a^6 \Gamma^2 + 41371 \gamma_a^5 \Gamma^3 + 58786 \gamma_a^4 \Gamma^4 + 41371 \gamma_a^3 \Gamma^5 + 14493 \gamma_a^2 \Gamma^6 +2345 \gamma_a \Gamma^7 + 150 \Gamma^8\Big)}{6 \Gamma \gamma_a (\Gamma +\gamma_a) (2 \Gamma + \gamma_a) (\Gamma + 2 \gamma_a)(3 \Gamma + \gamma_a)(\Gamma + 3 \gamma_a) (5 \Gamma +\gamma_a) (\Gamma + 5\gamma_a)}$$

Maybe not particularly enlightening, but that is the result. This can also be written as a sum of inverse combinations of the two emission rates:

$$\langle t_1^{(1)} \rangle =\frac{2}{\Gamma +\gamma_a}+\frac{1}{\Gamma }+\frac{1}{\gamma_a}$$

$$\langle t_1^{(2)} \rangle =\frac{1}{8} \left(\frac{8}{\Gamma +\gamma_a}+\frac{9}{3 \Gamma +\gamma_a}+\frac{9}{\Gamma +3 \gamma_a}+\frac{4}{\Gamma }+\frac{4}{\gamma_a}\right)$$

$$\langle t_2^{(2)} \rangle =\frac{3}{8} \left(\frac{8}{\Gamma +\gamma_a}-\frac{3}{3 \Gamma +\gamma_a}-\frac{3}{\Gamma +3 \gamma_a}+\frac{4}{\Gamma }+\frac{4}{\gamma_a}\right)$$

$$\langle t_1^{(3)} \rangle =\frac{1}{243} \left(2 \left(\frac{62}{2 \Gamma +\gamma_a}+\frac{125}{5 \Gamma +\gamma_a}+\frac{62}{\Gamma +2\gamma_a}+\frac{125}{\Gamma +5 \gamma_a}+\frac{81}{\Gamma +\gamma_a}\right)+\frac{81}{\Gamma}+\frac{81}{\gamma_a}\right)$$

$$\langle t_2^{(3)} \rangle =\frac{1}{1944}\Bigg(\frac{3240}{\Gamma +\gamma_a}-\frac{1984}{2 \Gamma +\gamma_a}+\frac{6561}{3 \Gamma +\gamma_a}-\frac{4000}{5 \Gamma +\gamma_a}-\frac{1984}{\Gamma +2 \gamma_a}+\frac{6561}{\Gamma +3 \gamma_a}+{}\\{}-\frac{4000}{\Gamma +5 \gamma_a}+\frac{1620}{\Gamma }+\frac{1620}{\gamma_a}\Bigg)$$

$$\begin{multline} \langle t_3^{(3)} \rangle =\frac{1}{1944}\Bigg(\frac{7128}{\Gamma +\gamma_a}+\frac{992}{2 \Gamma +\gamma_a}-\frac{6561}{3 \Gamma +\gamma_a}+\frac{2000}{5 \Gamma +\gamma_a}+\frac{992}{\Gamma +2 \gamma_a}+{}\\{}-\frac{6561}{\Gamma +3 \gamma_a}+\frac{2000}{\Gamma +5 \gamma_a}+\frac{3564}{\Gamma }+\frac{3564}{\gamma_a}\Bigg) \end{multline}$$

This is the formula for the average of the squared time of the $k$th photon from a $N$-photon bundle:

$$\begin{multline} \langle \big(t_k^{(N)}\big)^2\rangle = \frac{2\gamma_a^N}{P(N,N)}\left( \frac{\Gamma}{\Gamma_-} \right)^{2N}k {N\choose k} \sum_{k_1+k_2+k_3=N-k}{N-k\choose k_1, k_2, k_3}\\{}\times \sum_{k_4+\cdots+k_9=k-1}{k-1\choose k_4, k_5, k_6, k_7,k_8,k_9} \sum_{k_{10}+k_{11}=2}{2\choose k_{10}, k_{11}}\\{}\times \frac{(-1)^{k_3 + k_6 + k_7 + k_{8}+k_{11}}}{\gamma_a^{k_1+k_4+k_7} \Gamma^{k_2+k_5+k_8} \left( \frac{\Gamma_+}{4} \right)^{k_3 + k_6 +k_9} ( \gamma_a (k_1+k_7+k_{11}/2) +\Gamma(k_2+k_8 +k_{10}/2) + \frac{\Gamma_+}{2} (k_3+k_9))^3} \end{multline}$$

$$\begin{tabular}{r|cccccccccc} &1&&&&&\hskip-.5cm $k$&&&&10\\ \hline 1&1 & & & & & & & & & \\ &${1\over2}$ & ${3\over2}$ & & & & & & & & \\[.1cm] &${1\over3}$ & ${5\over6}$ & ${11\over6}$ & & & & & & & \\[.1cm] &${1\over4}$ & ${7\over12}$ & ${13\over12}$ & $25\over12$ & & & & & & \\[.1cm] &${1\over5}$ & ${9\over20}$ & ${47\over60}$ & $77\over60$ & $137\over60$ & & & & & \\%[.1cm] \raisebox{.2cm}{$N$}&${1\over6}$ & ${11\over30}$ & ${37\over60}$ & $19\over20$ & $29\over20$ & $49\over20$ & & & & \\[.1cm] &${1\over7}$ & ${13\over42}$ & ${107\over210}$ & $319\over420$ & $153\over140$ & $223\over140$ & $363\over140$ & & & \\[.1cm] &${1\over8}$ & ${15\over56}$ & ${73\over168}$ & $533\over840$ & $743\over840$ & $341\over280$ & $481\over280$ & $761\over280$ & & \\[.1cm] &${1\over9}$ & ${17\over72}$ & ${191\over504}$ & $275\over504$ & $1879\over2520$ & $2509\over2520$ & $3349\over2520$ & $4609\over2520$ & $7129\over2520$ & \\[.1cm] 10&${1\over10}$ & ${19\over90}$ & ${121\over360}$ & $1207\over2520$ & $1627\over2520$ & $2131\over2520$ & $2761\over2520$ & $3601\over2520$ & $4861\over2520$ & $7381\over2520$ \\[.1cm] \end{tabular}$$

— • 16:48, 18 November 2021