Multichannel decay law

It is well known, both theoretically and experimentally, that the survival probability for an unstable quantum state, formed at $t=0,$ is not a simple exponential function, even if the latter is a good approximation for intermediate times. Typically, unstable quantum states/particles can decay in more than a single decay channel. In this work, the general expression for the probability that an unstable state decays into a certain $i$-th channel between the initial time $t=0$ and an arbitrary $t>0$ is provided, both for nonrelativistic quantum states and for relativistic particles. These partial decay probabilities are also not exponential and their ratio turns out to be not a simple constant, as it would be in the exponential limit. Quite remarkably, these deviations may last relatively long, thus making them potentially interesting in applications. Thus, multichannel decays represent a promising and yet unexplored framework to search for deviations from the exponential decay law in quantum mechanical systems, such as quantum tunnelling, and in the context of particle decays.

Decay processes are ubiquitous in the realm of Quantum Mechanics (QM) and Quantum Field Theory (QFT). Quite interestingly, they comprise utterly different time scales, which range from the very short lifetime of the ρ meson (τ ∼ 10 −23 s) and the Higgs boson (τ ∼ 10 −22 s ), to the extremely long lifetimes of some nuclei, as for instance the double-β decays of 130 Te with τ ∼ 10 20 y and 124 Xe with τ ∼ 10 22 y [1]. It is then remarkable that 'in nuce' a very similar theory of decay can be applied to such different physical systems.
The survival probability p(t) is the probability that an unstable state in QM or an unstable particle in QFT created at t = 0 is still undecayed at the time t > 0. It can be expressed by the following famous formula (e.g. Ref [2] and the following for an explicit QM derivation): where d S (E) is the the so-called energy (or mass) probability density and E th,1 its minimal low-energy threshold. Namely, an unstable state is not an energy eigenstate of the underlying Hamiltonian of the system and dEd S (E) can be interpreted as the probability that it has an energy between E and E + dE. Note, for a stable state d S (E) = δ(E − M ): one recovers p(t) = 1 for each t.
The decay law is typically quadratic for short times, p(t) ≃ 1 − t 2 /τ 2 Z , where τ Z is the so-called Zeno time. As a consequence, for frequently repeated measurements at very small time intervals, the Zeno effect, that is the freezing of state in its unstable initial configuration, takes place [21][22][23][24][25][26][27], as confirmed in experiments on atomic transitions [28][29][30][31] as well as on quantum tunnelling [16,17]. For measurements repeated for slightly larger time intervals, also the inverse Zeno effect, that is an increased decay rate by measurements, is possible [11,12,32,33], as verified in Ref. [17]. Finally, at large times a power law sets in, p(t) ≃ t −α where α > 0 depends on the details of the interaction that governs that particular decay [19]. All these well-known features deal with the survival probability of the unstable state expressed in Eq. (1) and apply to basically each unstable state.
Next, a simple question might be asked: Which is the probability that the decay of the unstable state occurs between t = 0 and the time t? The answer is trivial, since the probability that the decay has actually occurred, denoted as w(t), must be Similarly, the quantity h(t) = w ′ (t) = −p ′ (t) is the probability decay density, with h(t)dt being the probability that the decay occurs between t and t + dt. Yet, in many practical cases that apply to both unstable quantum states and particles, there is not a single decay channel, but N distinct ones. Typically, excited atoms can decay into different low-lying energy levels and elementary particles, such as the τ lepton or the Higgs boson as well as composite hadrons, display various decay channels. In these cases, in the exponential (or Breit-Wigner (BW)) limit the decay width Γ is the sum of N distinct terms, Γ = N i=1 Γ i , where Γ i is the partial decay width in the i-th channel and Γ i /Γ the corresponding branching ratio.
Then, a natural, but less easy question is the following: How to calculate, in a general fashion, the probability, denoted as w i (t), that the decay occurs in the i-th channel between 0 and t?
In the BW limit, the expected result is simply w i (t) = Γi Γ w(t), thus each decay probability w i (t) is a fraction of the total decay probability w(t), but, as we shall see, this is not valid in general. In the precursory work of Ref. [13] a partial approximate solution for the function w i (t) was put forward in specific models and in the recent Ref. [34] the two-channel decay was studied via tunnelling to the 'left' and to the 'right' in an asymmetric double-delta potential. The aim of this work is to present the general form of the i-th channel decay probability w i (t).
Let us start with the QM case. We consider a system that contains an unstable state |S that can decay into final states of the type |E, i ,where i = 1, ..., N enumerates the decay channels (E being the energy of the final state, see below). By denoting with H the Hamiltonian of the system, the time-evolution operator e −iHt/ is clearly of crucial importance for our purposes, since it controls how the unstable state |S evolves in time. First, we rewrite it as (t > 0): out of which the survival probability amplitude is obtained by the expectation value The survival probability of Eq. (1) emerges as The self-energy function Π(E) is the one-state irreducible contribution. Intuitively, it represents the transition |S → |E, i → |S where an integral/sum over E and i is taken. (With no loss of generality, we adopt the choice Re Π(M ) = 0 by a suitable subtraction, thus M is the nominal energy of the state or, more formally, the value at which Re[G −1 S (E = M )] = 0). As a consequence of the optical theorem (e.g. Ref. [35]), the decay width function Γ(E) = 2 Im Π(E) depends on the energy (only in the exponential limit it is a simple constant). The 'on shell' BW decay width is recovered as an approximation by setting E = M , thus Γ = Γ(M ).
The spectral function (or mass distribution) d S (E), already introduced in the very first Eq. (1), is obtained from the last equality of Eq. (5) as: Intuitively, Eq. (5) shows that the propagator G S (E) is rewritten as the 'sum' of free propagators with 'weight' d S (E ′ )dE ′ , which as mentioned above is naturally interpreted as the probability that the state |S has an energy contained in the range (E, E + dE). As a consequence, the probability decay amplitude in Eq. (4) is just the integral over d S (E) weighted by the corresponding time-evolution factor e −iEt/ℏ . The normalization ∞ E th,1 dEd S (E) = 1 is guaranteed for any physical problem (for an explicit proof, see e.g. Refs. [36,37]).
When N decay channels are available, Π(E) is given by the sum over them: where Π i (E) encodes the loop |S → |E, i → |S for a fixed decay channel and Γ i (E) is the i-th decay width function, that vanishes below the corresponding threshold E th,i . Here, we assume for definiteness that E th,1 ≤ E th,2 ≤ ... ≤ E th,N . The functions Γ i (E) are -in principle-obtainable for a given problem, once the particular Hamiltonian is known. Of course, in practice they are often not known exactly and their explicit determination may be confined to a certain energy range and within the validity of the used approach to evaluate them (as, for example, perturbation theory). The real part Re Π i (E) can be obtained by dispersion relations 1 . The partial BW widths are We also rewrite the spectral function as the sum over partial spectral functions where dEd S (E) is the probability that the state |S has an energy between (E, E + dE) and decays in the i-th channel. Then, the integral is easily interpreted as the asymptotic branching ratio for the decay in the i-th channel. In models, the BW branching ratio Γ i /Γ represents a good (although not exact) estimate of the value r i , see later on for an example. Once the functions Γ i (E) are (assumed to be) known, we can map the quantum decay of the unstable state |S onto an effective Lee Hamiltonian, see e.g. [13,36,[38][39][40][41][42][43][44][45][46], which reproduces the propagator of Eq. (5) exactly: The ket |E, i is the state describing the i-th channel decay product with E being the eigenvalue of the non-interacting part of H L . Clearly, the difficulty lies in the previous determination of the decay where C is a subtraction constant. Re Π i (E) is obtained by taking the principal part. In some models, it may be convenient to calculate directly the loop function Π i (E).
functions Γ i (E), which are used here as an input. The searched probability w i (t) that the unstable state decays in the i-th channel between 0 and t is formally defined as: Namely, it is the probability that the original state |S has evolved into the decay product |E, i at the time t > 0. An explicit calculation of the matrix element in the integrand delivers [13,36] thus w i (t) is obtained as: Upon introducing the spectral function as in Eq. (5), w i (t) can be expressed as: or equivalently The expressions w i (t) in Eqs. (13), (14), (15) represent the main result of this work and can be easily evaluated numerically, since they involve only calculable quantities, such as the Γ i (E) and M (which fix also d S (E) and a(t)). Moreover, the limiting case t → ∞ can be calculated by taking into account the latter being the asymptotic branching ratio in the i-th channel of Eq. (9). Moreover, from Eq. (14) one has w i (0) = 0 and the small t expansion reads w i (t) ≃ c i t 2 with 2π 2 c i = ∞ E th,i dEΓ i (E) (one recovers also the Zeno time as 2π 2 τ −2 Z = ∞ E th,1 dEΓ(E)). The ratio w 1 /w 2 ≃ c 1 /c 2 applies for short times, but w 1 /w 2 is approximated by Γ 1 /Γ 2 at intermediate times, showing that in general it cannot be a constant, as confirmed below in a numerical example. Of course, the equality w = N i=1 w i = 1 − p(t) holds for each t. Another quantity of interest is the partial i-th probability decay density h i (t), evaluated as h i (t) = w ′ i (t): h i (t)dt is the probability that the unstable quantum state decays in the i-th channel between t and t + dt. Clearly, h(t) = −p ′ (t) = N i=1 h i (t) also applies. In the renowned BW limit [47][48][49] with Γ = N i=1 Γ i , thus r i = Γ i /Γ is exact in this limit. The temporal evolution is exactly exponential with a(t) = e −i M ℏ t− Γ 2ℏ t and p(t) = e − Γ ℏ t as well as: As already anticipated, the limit w i (∞) = r i is clearly fulfilled. In the BW limit, one has: This is however not true in general, see the following numerical examples and Refs. [13,34]. A useful (although not exact) approximation presented in Ref. [13] is given by It fulfills all the required limits (w i (t) = 0, w i (∞) = r i , p + N i=1 w i = 1, as well as the BW one) and, as numerical tests show, is close to the full result of Eq. (14). The functions w appr i (t) are easier to evaluate because they involve only one-dimensional integrals and can be useful in first approximation. The straight line corresponds to the BW limit Γ 1 /Γ 2 , see Eq. (19).
For illustrative purposes, a numerical example that contains two decay channels is presented in Figs.

1-4 for Γ
This is a rather simple model that however contains the main needed features that each decay should have, see e.g. the study of the decay of an excited hydrogen atom in Ref. [50]. In particular, the square root E − E th,i appears quite naturally in various applications as a phase-space term.
the decay in the first channel is enhanced (the opposite is true for h 1 /h 2 < Γ 1 /Γ 2 ).
It should be stressed that the precise form of the plotted functions depends on the chosen physical system; moreover, in the presented example the functions Γ i (E) are assumed to be perfectly known, what in actual physical systems is generally not the case. Nevertheless, the most relevant features, such as the departures from the simple BW behavior, with special attention on the ratios w 1 /w 2 and h 1 /h 2 depicted in Fig. 2 and Fig. 4, are expected to be valid in general.
Next, we briefly present the main features that concern the extension to QFT. In practical terms, when going from QM to QFT, one needs to replace the energy E with the Mandelstam variable s = E 2 . (Intuitively, in QFT there are positive and negative energy solutions and the QM free propagator (E − M ) −1 becomes (s − M 2 ) −1 ). The general QFT propagator (in natural units) of an unstable particle S reads (see e.g. Ref. [35]): where Π(s) = N i=1 Π i (s) is the sum of the self energies for N distinct decay channels. The partial decay function reads now 2 Γ i (s) = Im Π i (s)/ √ s. The partial (on shell) BW decay widths are Γ i = Im Π i (M 2 )/M and the total one Γ = Im Π(M 2 )/M . Just as in the QM case, the spectral function is the imaginary part of the propagator, d S (s) = − 1 π Im G S (s), thus as function of the energy d S (E) = − 2E π Im G S (s = E 2 ) [51,52], where the factor 2E arises by the variable change. The normalization ∞ E i,th dEd S (E) = 1 holds also in QFT [37] and the survival probability p(t) takes the same form of Eq. (1) [53,54]. A complete study of the QFT case with all technical details is left as an important outlook of the present work.
In QFT, the determination of Γ i (s) (and then Π(s)) is, in general, a highly nontrivial problem. Perturbative approaches are valid for small s, but fail for large energies (this is the main reason why perturbative approach do not capture deviations from the exponential [55]), where nonperturbative methods (such as unitarization schemes [56][57][58]) are needed. Just as in QM, we assume that a suitable form of Γ i (s) is known for a certain specific problem. Then, analogously to Eq. (14), the final result for w i (t) is: which can be easily evaluated numerically. The approximate expression w appr is also applicable in QFT. A useful testing case is realized for Γ i (s) = g 2 i (s − s th,i )/s with s th,i = E 2 th,i (the real part is a constant that can be set to zero [59]), which gives a similar qualitative behavior of the curves depicted in Figs. 1-4. Interestingly, this choice for Γ i (s) has been recently successfully applied to the description of some strongly decaying resonances (ρ-meson, a 1 -meson, ∆-baryon) in Ref. [59], delivering better fits than (non)relativistic BW functions.
In conclusion, for the case of an unstable state with N different available decay channels, we have provided the expressions for the decay probability in each channel in QM (Eqs. (14) and (15)) and QFT (Eq. (22)) and we have presented an illustrative numerical example in Figs. 1-4. These equations can be tested in various experimental setups, such as an asymmetric tunneling potential, along the line of the already performed experiment in Ref. [16] and in agreement with the recent simulation of Ref. [34]. Presently, the possibility to manipulate potentials, e.g. Refs. [60][61][62][63][64][65][66][67], opens up interesting developments. Alternatively, the experiment described in Ref. [20] is also very promising: the non-exponential decay is mapped into arrays of single-mode optical devices, in which the time evolution corresponds to spatial evolution; it is therefore conceivable to realize a similar experiment by engineering two decay channels, in such a way to 'measure' the decay probabilities w i (t).
The study of multiple decay channels for elementary particles is also an important and ambitious goal, even though it is typically difficult to measure deviations from the exponential decay due to the short lifetimes involved. The ratios w i /w j and h i /h j are expected to deviate longer from the BW limit and could therefore be the key to see such deviations in such fundamental systems.