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}$$