Investigation of two photon emission in strong field QED using channeling in a crystal

TL;DR

This paper explores two-photon emission processes in strong field QED within crystal channeling, distinguishing between one-step and two-step mechanisms, and provides calculations relevant for high-energy electron experiments in crystalline environments.

Contribution

It offers the first detailed analysis of two-photon emission, including interference effects, in strong field QED channeling, extending understanding beyond single-photon emission studies.

Findings

Two-step emission probability equals the product of single-photon probabilities.

One-step contributions involve significant interference effects.

Thick crystals allow simplified multiple photon emission calculations.

Abstract

We investigate the 2nd order process of two photons being emitted by a high-energy electron dressed in the strong background electric field found between the planes in a crystal. The strong crystalline field combined with ultra relativistic electrons is one of very few cases where the Schwinger field can be experimentally achieved in the electron's rest frame. The radiation being emitted, the so-called channeling radiation, is a well studied phenomenon. However only the first order diagram corresponding to emission of a single photon has been studied so far. We elaborate on how the 2 photon emission process should be understood in terms of a two-step versus a one-step process, i.e., if one can consider one photon being emitted after the other, or if there is also a contribution where the two photons are emitted 'simultaneously'. From the calculated full probability we see that the two-step contribution is simply the product of probabilities for single photon emission while the additional one-step terms are, mainly, interferences due to several possible intermediate virtual states. These terms can contribute significantly when the crystal is thin. Therefore, in addition, we see how one can, for a thick crystal, calculate multiple photon emissions quickly by neglecting the one-step terms, which represents a solution of the problem of quantum radiation reaction in a crystal beyond the usually applied constant field approximation. We explicitly calculate an example of 180 GeV electrons in a thin Silicon crystal and argue why it is, for experimental reasons, more feasible to see the one-step contribution in a crystal experiment than in a laser experiment.

Equations

$$\chi =e\sqrt{(F_{\mu \nu }{p}^{\nu }{)}^2}/m^3,$$

$$S_{fi}=⟨f|U(\infty ,-\infty )|i⟩$$

$$\Psi =\sum _{s=1}^2\int \frac{d^{3}p}{(2\pi {)}^3}\left[b_p^s{\psi }_{p,s}^{-}(x)+c_p^{s†}{\psi }_{p,s}^{+}(x)\right],$$

$$/A=\int \frac{d^{3}k}{(2\pi {)}^3}\sqrt{\frac{4\pi }{2\omega }}\sum _{r=1}^2\left[/{ϵ}_r{a}_k^r{e}^{-ikx}+/{ϵ}_r^{\ast }{a}_k^{r†}{e}^{ikx}\right],$$

$${\psi }^{-}(x)=\frac{1}{\sqrt{2\epsilon }}{e}^{i(p_{x}x+p_{z}z-\epsilon t)}U(y),$$

$${\psi }^{+}(x)=\frac{1}{\sqrt{2\epsilon }}{e}^{-i(p_{x}x+p_{z}z-\epsilon t)}V(y)$$

$$U(y)=\sqrt{\varepsilon+m}\left(\begin{array}[]{c}\boldsymbol{s}^{-}\\ \frac{\boldsymbol{\sigma}\cdot\tilde{\boldsymbol{p}}}{\varepsilon+m}\boldsymbol{s}^{-}\end{array}\right)I^{-}(y)$$

$$V(y)=\sqrt{\varepsilon+m}\left(\begin{array}[]{c}\frac{\boldsymbol{\sigma}\cdot\tilde{\boldsymbol{p}}}{\varepsilon+m}\boldsymbol{s}^{+}\\ \boldsymbol{s}^{+}\end{array}\right)I^{+}(y)$$

$$\left[-\frac{1}{2\epsilon }\frac{d^2}{d{y}^2}+q\phi (y)\right]I(y)=\frac{{\epsilon }^2-p_x^2-p_z^2-m^2}{2\epsilon }I(y),$$

$$I^{-}(y)=e^{i{k}_B^{-}y}{u}_{k_B}^{-}(y),$$

$$\left(-\frac{1}{2\epsilon }\left[\frac{d^2}{d{y}^2}+2i{k}_B\frac{d}{dy}-k_B^2\right]+q\phi (y)\right)u_{k_B}^{-}(y)$$

$$=\frac{{\epsilon }^2-p_x^2-p_z^2-m^2}{2\epsilon }{u}_{k_B}^{-}(y).$$

$$u_{k_B}^{-}(y)=\sum _j{c}_j{e}^{ij{k}_{0}y},$$

$$E_n=\frac{{\epsilon }^2-p_x^2-p_z^2-m^2}{2\epsilon },$$

$$\sum _j\frac{1}{2\epsilon }{\left[j{k}_0+k_B\right]}^2{\delta }_{j,l}{c}_j$$

$$+\sum _j{c}_j\frac{1}{d_p}\int q\phi (y)e^{i(j-l)k_{0}y}dy$$

$$=\sum _j\frac{{\epsilon }^2-p_x^2-p_z^2-m^2}{2\epsilon }{\delta }_{j,l}{c}_j.$$

$$U(y)=\sum _j{c}_j^{-}{\mathbf{S}}_j^{-}{e}^{i(j{k}_0+k_B)y},$$

$$V(y)=\sum _j{c}_j^{+}{\mathbf{S}}_j^{+}{e}^{-i(j{k}_0+k_B)y},$$

$$\boldsymbol{S}_{j}^{-}=\sqrt{\varepsilon+m}\left(\begin{array}[]{c}\boldsymbol{s}^{-}\\ \frac{\boldsymbol{\sigma}\cdot\boldsymbol{p}_{j}}{\varepsilon+m}\boldsymbol{s}^{-}\end{array}\right),$$

$$\boldsymbol{S}_{j}^{+}=\sqrt{\varepsilon+m}\left(\begin{array}[]{c}\frac{\boldsymbol{\sigma}\cdot\boldsymbol{p}_{j}}{\varepsilon+m}\boldsymbol{s}^{+}\\ \boldsymbol{s}^{+}\end{array}\right),$$

$${\bar{\mathbf{S}}}_f^{-}/{ϵ}^{\ast }{\mathbf{S}}_i^{-}=-{\mathbf{s}}_f^T\left[\mathbf{ϵ}\cdot \mathbf{A}+i\mathbf{B}\cdot \mathbf{\sigma }\right]{\mathbf{s}}_i$$

$$\mathbf{A}=\sqrt{\frac{{\epsilon }_f+m}{{\epsilon }_i+m}}{\mathbf{p}}_i+\sqrt{\frac{{\epsilon }_i+m}{{\epsilon }_f+m}}{\mathbf{p}}_f,$$

$$\mathbf{B}={\mathbf{ϵ}}^{\ast }\times \left(\sqrt{\frac{{\epsilon }_f+m}{{\epsilon }_i+m}}{\mathbf{p}}_i-\sqrt{\frac{{\epsilon }_i+m}{{\epsilon }_f+m}}{\mathbf{p}}_f\right).$$

$$G(x_2,x_1)$$

$$\theta (t_2-t_1)e^{-i\epsilon (t_2-t_1)}{\psi }_{p,n,s}^{-}({\mathbf{x}}_2){\bar{\psi }}_{p,n,s}^{-}({\mathbf{x}}_1)$$

$$-\theta (t_1-t_2)e^{i\epsilon (t_2-t_1)}{\psi }_{p,n,s}^{+}({\mathbf{x}}_2){\bar{\psi }}_{p,n,s}^{+}({\mathbf{x}}_1).$$

$$d{W}_{i\to f}^{(1)}$$

$${\mathcal{M}}_{i\to f}=e\sqrt{\frac{4\pi }{2\omega }}\frac{1}{2\sqrt{{\epsilon }_f{\epsilon }_i}}\sum _j{c}_{n_B+j,f}^{\ast }{c}_{j,i}{\bar{\mathbf{S}}}_{n_B+j,f}^{-}/{ϵ}^{\ast }{\mathbf{S}}_{j,i}^{-},$$

$$f(\theta )$$

$$\simeq E_{n_f}-E_{n_i}+\frac{m^2}{2{\epsilon }_f}-\frac{m^2}{2{\epsilon }_i}+\frac{\omega {\theta }^2}{2}\left(1+\frac{\omega {\text{sin}}^2\phi }{{\epsilon }_f}\right).$$

$$\frac{d{P}_{i\to f}^{(\text{cascade})}}{d{\omega }_{1}d{\omega }_2}$$

$$+\frac{d{W}_{i\to v}^{(1)}}{d\omega }({\omega }_2)\frac{d{W}_{v\to f}^{(1)}}{d\omega }({\omega }_1)].$$

$$S_{i\to f}^{(2)}$$

Full text

Strong field QED is the study of physical processes that take place in a strong background field and nonlinear effects of quantum nature arise when the size of the Lorentz invariant parameter

[TABLE]

is on the order of unity, which is the ratio of the electromagnetic field experienced in the electron’s rest frame compared to the Schwinger field strength $E_{\text{Sch}}=1.32\times 10^{18}\,\text{V/m}$. Here $e$ is the elementary charge, $m$ the electron mass, $F_{\mu\nu}$ the electromagnetic field tensor of the background field and $p^{\nu}$ the electron 4-momentum. We use natural units such that $\hbar=c=1$, $\alpha=e^{2}$. Lindhard was one of the first to realize that when high energy charged particles are aimed close to the direction along an axis or plane in a crystal, the charged particle can become transversely trapped Lindhard (1965). Later it was studied how this motion leads to radiation emission called channeling radiation, especially relevant for electrons and positrons. This is well-studied both experimentally Bak et al. (1985, 1988); Swent et al. (1979); Andersen et al. (1981, 1982); Klein et al. (1985); Alguard et al. (1979); Andersen et al. (2012); Uggerhøj (2005) and theoretically Kumakhov (1976, 1977); Andersen et al. (1981); Sáenz et al. (1981); Kimball and Cue (1985). Crystal channeling represents one of the only phenomena where the Schwinger field can be experimentally achieved in the electron’s rest frame Belkacem et al. (1986); Andersen et al. (1982); Esberg et al. (2010); Wistisen et al. (2018), with the only other example being the famous E-144 SLAC experiment on non lienar Compton scattering Bula et al. (1996) using relativistic electrons colliding with a laser beam. Crystals with ultra relativistic electrons or positrons therefore present a unique possibility to study physics in such strong fields. However a calculation from first principles of emission of more than 1 photon has not been carried out for crystal channeling. The recent studies of 2 photon emission in the collision of relativistic electrons with a laser pulse Seipt and Kämpfer (2012); Mackenroth and Di Piazza (2013); Dinu and Torgrimsson (2018a); King (2015) show that the emission of 2 photons is not exactly the product of probabilities for each emission, however under certain conditions it is an acceptable approximation. The experimental verification of such results are however complicated in the case of the laser pulse colliding with an electron bunch because any two (or more) emitted photons cannot be known to be emitted by the same electron. In crystal experiments as in e.g. Wistisen et al. (2018), it is standard that each incoming particle is recorded as a separate event, and therefore the measured outgoing photons are sure to stem from the single incoming particle. Therefore, in this paper, we will calculate the emission of 2 photons during electron channeling in a crystal, which could potentially be studied experimentally in an experiment similar to the one seen in Wistisen et al. (2018), however with a modified setup to allow for the detection of an additional photon. For the theory of channeling radiation, in particular the development of the semi-classical operator method by Baier et. al. Baier and Katkov (1968) stands out, and has been extensively applied to the phenomenon of channeling Baier et al. (1998).

This method allowed to include quantum effects such as the electron spin and the photon recoil, which are important when $\chi$ is no longer small, while needing only the classical trajectory of the electron/positron in the external field. The authors of this method, seeking analytical results, in most applications to channeling, applied the approximation of the local constant field which greatly simplifies calculations. The constant field approximation means that while a particle moves in an external field, which is not constant, one applies the result of the constant field formula locally, i.e. in a small time step. Effectively this means neglecting that the radiation emitted before or after can interfere with this radiation. This is valid only for certain parameters of fields and particle energies. However the semi-classical operator method can be used to calculate the radiation emission under general circumstances without much effort, also when the constant field approximation is no longer valid Wistisen (2014, 2015), which with modern computing power makes it one of the most powerful methods to calculate the radiation emitted by ultra-relativistic electrons in a general field configuration. There are caveats however, which are two-fold. Firstly, the notion of a classical trajectory should make sense. Or, in other words, the quantum numbers associated with the motion should be large, a subject recently studied in Wistisen and Di Piazza (2019, 2018). Secondly, the derivation starts out from the first-order diagram of a dressed electron emitting a single photon. Therefore the emission rate of two, or more, photons can not be predicted by this method without approximations. The emission of a single photon yields a rate, an emission probability per unit time, and as such one can construct the probability for emitting several photons by applying this rate for each consecutive emission. In this way, the probability to emit, e.g., two photons would be proportional to time, or thickness of the crystal, squared, and so on. We will call this process the ’cascade’ process. Herein lies an approximation, where interference between different emissions is neglected. We show that the two-photon emission probability contains the cascade along with one-step terms which scale linearly with the crystal thickness. Therefore, for sufficiently thin crystals, these one-step terms will become important. This phenomenon is also discussed in pair production of electron/positron pairs from high energy photons in a strong field where one also distinguishes between the two-step and the one-step, or ’trident’ process. This has been investigated in crystals in Esberg et al. (2010) and has received renewed interest with the prospect of studying such phenomena in high-intensity laser fields Hu et al. (2010); Ilderton (2011); Acosta and Kämpfer (2019); King and Fedotov (2018); Dinu and Torgrimsson (2018b); Mackenroth and Di Piazza (2018). In this paper we make quantitative calculations of the angularly integrated probability, differential in photon energies, of emission of two photons by an electron in the planar Doyle-Turner potential Baier et al. (1998); Doyle and Turner (1968); Avakian et al. (1982); Møller (1995). We do this by finding numerical solutions of the Dirac equation by solving the problem in a basis of plane waves, which is possible due to the periodicity of the transverse potential in a crystal, as shown in Wistisen and Di Piazza (2019). If the cascade terms are enough to properly describe the radiation emission is a highly relevant question as it closely relates to the phenomenon of quantum radiation reaction, the emission of multiple photons when $\chi$ is large, Di Piazza et al. (2010), recently studied using channeling radiation and in laser experiments Wistisen et al. (2018); Poder et al. (2018); Cole et al. (2018). In the crystal experiment it was seen that even for energies as high as $180$ GeV positrons, where it could be expected that the constant field approximation would be acceptable, it was shown that discrepancies arise due to this, and therefore a more general theory was called for. The current theory of quantum radiation reaction in lasers relies on the local constant field approximation Di Piazza et al. (2010); Neitz and Di Piazza (2013); Blackburn et al. (2014); Baier et al. (1998); Vranic et al. (2016); Li et al. (2014), and it is unknown if one can calculate the emission of many photons in a way that avoids calculating all the corresponding higher order diagrams, when going beyond the constant field approximation. This question will be addressed in the case of a crystal, in the current paper.

We use the Feynman slash notation such that $\not{a}=a_{\mu}\gamma^{\mu}$, where $\gamma^{\mu}$ are the Dirac gamma matrices and $a^{\mu}$ an arbitrary four-vector. We adopt the metric tensor $\eta^{\mu\nu}=\text{diag}(+1,-1,-1,-1)$.

I Formalism

In QED the transition amplitude from a given initial state $\left|i\right\rangle$ to a final state $\left|f\right\rangle$ is given by

[TABLE]

where $U$ is the time evolution operator, often written as $U(\infty,-\infty)=\mathcal{T}\text{exp}\left(-i\int_{-\infty}^{\infty}V(t)dt\right)$ where $\mathcal{T}$ is the time-ordering operator and $V(t)=\int e\bar{\Psi}\not{A}\Psi d^{3}x$ is the quantized interaction. We then write our quantized fields as

[TABLE]

[TABLE]

where $\psi_{p,s}^{-}(x)$ and $\psi_{p,s}^{+}(x)$ are an orthonormal and complete set of electron and positron solutions, respectively, in the background field. $\int\frac{d^{3}p}{(2\pi)^{3}}$ denotes a summation over all states, and $p$ the relevant quantum numbers which we will find later. The $b$, $c$ and $a$ operators are the annihilation operators of the electron, positron and photon field respectively, obeying the relations, that the only non-zero (anti-)commutators are $\left\{b_{p}^{r},b_{q}^{s\dagger}\right\}=\left\{c_{p}^{r},c_{q}^{s\dagger}\right\}=[a_{p}^{r},a_{q}^{s\dagger}]=(2\pi)^{3}\delta^{rs}\delta^{(3)}(\boldsymbol{p}-\boldsymbol{q})$, where the $\{\}$ brackets denote the anti-commutator and $[]$ the commutator.

In Wistisen and Di Piazza (2019, 2018) we discussed the Dirac equation with the potential found in the crystal, but we will here repeat the results we need in order to calculate the emission of 2 photons. It was found in Wistisen and Di Piazza (2019) that the electron solution can be written as follows

[TABLE]

and the positron solutions can then be written as (see appendix A)

[TABLE]

and $U$ and $V$ are given by

[TABLE]

[TABLE]

where $\tilde{\boldsymbol{p}}=\left(p_{x}+q\varphi(y),\text{sign}(q)i\frac{d}{dy},p_{z}\right)$, $\varphi(y)$ is the electrostatic potential, $q=\pm e$ is the charge, the superscript on $I(y)$ refers to the charge sign, $\boldsymbol{s}$ is a two component vector describing the spin, which we can choose as either \left(\begin{array}[]{cc}1&0\end{array}\right)^{T} or \left(\begin{array}[]{cc}0&1\end{array}\right)^{T} , corresponding to spin-up and spin-down respectively for the electron, and opposite for the positron. From the choice of the form of the spinors $U$ and $V$, it is also clear that $\varepsilon$ positive should be used (see appendix A). $I(y)$ is the solution to the equation

[TABLE]

For $\varphi(y)$ we will use the Doyle-Turner model Baier et al. (1998); Doyle and Turner (1968); Avakian et al. (1982); Møller (1995), chosen as symmetric around 0. In a crystal this potential $\varphi(y)$ is periodic with the period of the inter planar distance which we will denote as $d_{p}$. Because of this, the solution (for the electron) can be written as a Bloch wave such that

[TABLE]

and where $u_{k_{B}^{-}}^{-}(y)$ is also periodic with period $d_{p}$ and $k_{B}^{-}$ is the Bloch momentum, which can be taken to be in the interval $0\leq k_{B}^{-}<k_{0}$, $k_{0}=\frac{2\pi}{d_{p}}$. It then follows from Blochs theorem that these solutions form an orthogonal and complete set of solutions of Eq. (6). Inserting $I^{-}(y)$ of Eq. (7) into Eq. (6) gives us the equation governing $u_{k_{B}}^{-}(y)$

[TABLE]

The periodicity of $u_{k_{B}}(y)$ means it can be written as a Fourier series,

[TABLE]

To ensure normalization we should have $\sum_{j}\left|c_{j}\right|^{2}=1$ (see appendix B). It is now clear that this is an eigenvalue problem for each $k_{B}$ where the quantized eigenvalue is

[TABLE]

where $n$ is the quantum number corresponding to the value of this energy in ascending order and where [math] is the ground state. This equation leads to a quantization of e.g. $p_{x}$. The coefficients $c_{j}$ are found by solving the matrix eigenvalue problem obtained by inserting Eq. (9) in Eq. (8) and multiply with $\frac{1}{d_{p}}e^{-ilk_{0}y}$ and integrate over $y$ from [math] to $d_{p}$ to exploit orthogonality

[TABLE]

This was done with the electron function $I^{-}(y)$ in mind, but the positron coefficients can be obtained just by changing $q$. With these things taken into consideration, we now see that we can write the $U(y)$ and $V(y)$ functions in terms of the coefficients $c_{j}$ such that

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

where $\boldsymbol{p}_{j}=(p_{x}+E_{n}-\frac{(jk_{0}+k_{B})^{2}}{2\varepsilon},jk_{0}+k_{B},p_{z})$. For the calculation of radiation emission from electrons we will need the quantity $\bar{\boldsymbol{S}}_{f}^{-}\not{\epsilon}^{*}\boldsymbol{S}_{i}^{-}$, where we have put labels for the initial state $i$ and final state $f$, however these still each depend on the index $j$. This quantity can then be written as

[TABLE]

where

[TABLE]

[TABLE]

Now since we have an orthonormal complete set of solutions, we can write the propagator in terms of these states as Beresteckij et al. (2008)

[TABLE]

This expression can be simplified due to the simple expression for the wave functions in all coordinates but the $y$ coordinate. However, we will not carry this out, as it is easier to see how the cascade part of the radiation emission arises by starting from the above expression.

II Single photon emission and cascade

We will now briefly mention some results obtained in Wistisen and Di Piazza (2019) on the single photon emission probability which is relevant to build the expected cascade contribution. We found that the rate of emission is given by

[TABLE]

where we defined

[TABLE]

where $n_{B}$ is the integer such that $0\leq k_{B,f}<k_{0}$, where $k_{B,f}=k_{B,i}-k_{y}-n_{B}k_{0}$, $\boldsymbol{S}_{j,i}$ corresponds to the initial state and $c_{j,i}$ is coefficient with index $j$ corresponding to the initial state $i$. See the appendix of Wistisen and Di Piazza (2019) for the details on why $\mathcal{M}$ reduces to a single sum over $j$. As shown in Wistisen and Di Piazza (2019) there are large terms in $\varepsilon_{f}+\omega-\varepsilon_{i}$ which cancel, leaving behind the relevant small terms, because the relevant transverse energies $E_{n}$, comparable to the potential depth, are much smaller than the whole particle energy i.e. eV versus GeV. We could rewrite the content of the delta function as

[TABLE]

Now we may use that $\delta(\varepsilon_{f}+\omega-\varepsilon_{i})=\frac{1}{|f^{\prime}(\theta_{0})|}\delta(\theta-\theta_{0})$ where $\theta_{0}$ is the positive solution to $f(\theta)=0$. From the formula for single photon emission, Eq. (19), we can construct the cascade contribution to two photon emission. We wish to know the probability of finding a photon in the momentum interval $d^{3}k_{1}$ around $k_{1}$ while also finding a photon within another interval $d^{3}k_{2}$ around $k_{2}$. This can happen in two ways, either the particle emits $k_{1}$ while transitioning from the initial state, and then subsequently $k_{2}$ or vice versa. We are however interested in the angular integrated spectrum, that is $dP_{i\rightarrow f}^{(\text{cascade})}/d\omega_{1}d\omega_{2}$ and therefore an additional factor of $\frac{1}{2}$ must be added due to counting the same point in phase space twice Feynman (1965), and so we obtain

[TABLE]

III Two photon emission

Expanding the time evolution operator to second order, allowing for two photon emission we have that the S-matrix element is

[TABLE]

When specifying the final state as $\left\langle p_{f},k_{1},k_{2}\right|$, an electron and two photons and the initial state as just an electron, $\left|p_{i},0,0\right\rangle$, $S_{fi}^{(2)}$ can be rewritten in terms of the wave functions and the propagator. In Beresteckij et al. (2008) this is done for the Compton scattering matrix element, which is the same diagram as here, except that an incoming photon is instead outgoing. The matrix element is therefore

[TABLE]

Now we define

[TABLE]

where $\mathcal{M}$ is defined as in Eq. (20) where $v$ is used to denote the virtual state from the propagator, and is shorthand for the dependence on $p_{x,v}$, $k_{B,v}$, $p_{z,v}$, $n_{v}$ and $s_{v}$. The superscript $-$ on $M^{-}$ and $\mathcal{M}^{-}$ denotes that the virtual state is the electron state $\psi_{v}^{-}$, and $M^{+}$, $\mathcal{M}^{+}$ is the same but with the positron virtual state. The matrix element may then be written as

[TABLE]

Therefore the term in the second line is seen as the electron first emits a photon with momentum $\boldsymbol{k}_{1}$ at $t_{1}$ and then propagates to a later time $t_{2}$ and emits a second photon with momentum $\boldsymbol{k}_{2}$. The term in the third line is then the electron emitting the photon with momentum $\boldsymbol{k}_{1}$ at a time $t_{1}$ turning the electron into a positron going into the past and emitting the photon with momentum $\boldsymbol{k}_{2}$ at the earlier time $t_{2}$. This last term is heavily suppressed in our case which we can see as follows. Denote $a=\varepsilon_{v}+\omega_{1}-\varepsilon_{i}$ and $b=\varepsilon_{f}+\omega_{2}-\varepsilon_{v}$, then we may use that

[TABLE]

where $\epsilon$ is a small real number for which one in the end should take the limit $\epsilon\rightarrow 0$ and therefore we have

[TABLE]

We have also that $-\theta(t_{2}-t_{1})=\frac{i}{2\pi}\int\frac{1}{\varepsilon_{V}-i\epsilon}e^{i(t_{2}-t_{1})\varepsilon_{V}}d\varepsilon_{V}$ and therefore we have the term from the third line of Eq. (26) carries the factor of

[TABLE]

and therefore this term will always be very far off-shell, as the virtual particle on-shell condition can never be met as it corresponds to the spontaneous production of an electron, positron and photon from the crystal field, where the produced positron is subsequently annihilated with the incoming electron to emit another photon. Having carried out the integrations over time we obtain that

[TABLE]

Now we may integrate over $p_{v}$ to obtain

[TABLE]

and then $v_{1}$ denotes the virtual state with momentum given by $p_{x,v_{1}}=p_{x,i}-k_{x,1}$, $p_{z,v_{1}}=p_{z,i}-k_{z,1}$ and $k_{B,v_{1}}^{-}=k_{B,i}^{-}-k_{y,1}-n_{B,1}^{-}k_{0}$ and $-k_{B,v_{1}}^{+}=k_{B,i}^{-}-k_{y,1}-n_{B,1}^{+}k_{0}$, i.e. that photon with label $1$ is emitted at the vertex connected with the initial particle. From the amplitude we get the transition probability according to

[TABLE]

where we have added a factor of $1/2$ in front due to identical particles in the final state, and that we in the end want to integrate over all angles, and would therefore, again, be counting double Feynman (1965). From this full result, it is seen that the result can diverge when $\epsilon\rightarrow 0$ because $\varepsilon_{i}-\varepsilon_{v}-\omega_{1}=0$ is possible. The nature of the divergence is however different for some of the terms, namely the ones which are the norm square of each term underneath the sum, $\left|\mathcal{M}_{i\rightarrow v_{1}}^{-}\mathcal{M}_{v_{1}\rightarrow f}^{-}/\left(\varepsilon_{i}-\varepsilon_{v_{1}}-\omega_{1}-i\epsilon\right)\right|^{2}$, where the limit of $\epsilon\rightarrow 0$ will yield an infinite result, even after integration over one of the angles $\theta_{1}$ or $\theta_{2}$. On the other hand, while the remaining terms, of the interference type, still diverge, they can be integrated over $\theta_{1}$ or $\theta_{2}$ to yield a convergent result. To learn the meaning of this divergence due to the denominator, see also Hu (2011), we may write

[TABLE]

and note that

[TABLE]

if we evaluate the integrals of $\mathcal{M}_{i\rightarrow v_{1}}^{-}\mathcal{M}_{v_{1}\rightarrow f}^{-}$ with the factor $\delta(a+b)\delta(b)$ we get well defined results, as this just amounts to the product of two 1.st order emissions. It is therefore useful to write

[TABLE]

where then the factor $\epsilon/(b^{2}+\epsilon^{2})$ acts like a delta-function for small enough $\epsilon$, yielding a finite value when we perform the integrals in Eq. (31), and then it is clear that this is divergent as $\epsilon\rightarrow 0$ due to the factor of $1/\epsilon$. However this should be understood in terms of an additional factor of $T$ for this term. To see this, consider the origin of this expression from Eq. (28), but consider instead that we had a finite time, and integrate over $a$ and $b$

[TABLE]

and we also have that

[TABLE]

and so we see that we must replace $\left|\frac{1}{b-i\epsilon}\right|^{2}\rightarrow\pi T\delta(b)$, and therefore these terms turn out to give us the cascade contribution. To see how the probability from Eq. (31) splits up into this cascade along with additional terms, we will denote the quantity underneath the norm-square as $R^{-}=\sum_{n_{v}}J_{1}^{-}(n_{v})+J_{2}^{-}(n_{v})$ corresponding to the terms with the virtual electron and similarly $R^{+}=\sum_{n_{v}}J_{1}^{+}(n_{v})+J_{2}^{+}(n_{v})$, where

[TABLE]

[TABLE]

and $J_{2}$ is $J_{1}$ with $(\epsilon_{1},k_{1})\leftrightarrow(\epsilon_{2},k_{2})$ and we then define $R=R^{+}+R^{-}$. The quantity we want is then $|R|^{2}=|R^{+}|^{2}+|R^{-}|^{2}+2\text{Re}\left[R^{+}\left(R^{-}\right)^{*}\right]$. In the $R^{+}$ term, it is never possible for the denominator to become [math] and therefore it can be directly calculated (see appendix C). For

[TABLE]

the product of the terms with the same subscript and where $n_{v}=n^{\prime}_{v}$ are the cascade which are the only problematic terms and so need special attention as described above. Therefore it is useful to employ that

[TABLE]

and so the terms in the first two lines are convergent contributions to the one-step process and the terms on the last line are the cascade terms, except that the spin sum is still underneath the norm-square. In appendix D we show that the interference due to spin will be [math] when the photon polarization can be taken as real and that either the sum over initial or final spins (we will do both) is carried out. And so we can write the differential probability of emission, with a given initial state, as

[TABLE]

IV Choice of regularization

Consider the terms proportional to $T^{2}$ from the above equation

[TABLE]

which by comparison with Eq. (22) and Eq. (19) is seen to be in agreement with the expected cascade result. Above we chose a certain way to regularize the divergence, by recognizing that the divergent terms correspond to the cascade terms, and that in taking the time limit from $\pm\infty$, some information about the duration of interaction was lost, which we put back in, in a way that is correct when $T$ is large enough i.e. larger than the photon formation length roughly estimated by $l_{f}=2\gamma^{2}(1-\omega/\varepsilon)/\omega$ Baier et al. (1998), which in our case is roughly $\gamma/m\sim 0.8$ µm, because $\omega$ is on the order of $\varepsilon$. Another way often found in literature Oleinik (1967, 1968); Roshchupkin (1996); Lötstedt et al. (2007); Gonthier et al. (2014) is to say that the virtual state is unstable and therefore replace the energy of the virtual particle according to $\varepsilon_{v}\rightarrow\varepsilon_{v}-i\Gamma_{v}/2$ where $\Gamma_{v}=\sum_{f}W_{v\rightarrow f}$ is the total decay width of the virtual state from all processes. This is equivalent to adding the effect of the line width in atomic Raman scattering Bransden et al. (2003). Effectively this corresponds to replacing the $\epsilon$ in the denominator with $-\Gamma_{v}/2$ which lifts the divergence. However one can see that with this substitution, see Eq. (34), one would obtain that

[TABLE]

where $f(b)$ is a function peaked around $b=0$ which obeys $\int f(b)db=1$ and therefore resembles the delta-function $\delta(b)$, but with a non-zero width $\Gamma_{v}$. If we then again calculate the cascade part according to this we would obtain

[TABLE]

and if we assume that the dominant contribution to the decay width is due to radiation emission we have that the total width is

[TABLE]

therefore if we approximate $f(b)\simeq\delta(b)$ and integrate over $d^{3}k_{2}$ and sum over $n_{f}$ we will obtain a factor of the total rate $\Gamma_{v_{1}}$, which cancels out, and so we have that

[TABLE]

which is just the single photon emission probability. Therefore this approach leads to the prediction that it is just as likely to emit 2 photons as it is 1. This is not a meaningful result and the reason is that the integration over time has been carried out over all times, i.e. it is assumed that $T\gg 1/\Gamma_{v}$ which means it is guaranteed that the virtual state decays. However in that case not only 2 photon emission is likely, also larger number of photons, which we do not take into account. For Raman scattering the approach is reasonable when $T\gg 1/\Gamma_{v}$ such that it is guaranteed that an excited state will decay before the observation is made. However if the interaction time is very short $T\ll 1/\Gamma_{v}$, it is also expected that Raman scattering should have a dependence as $T^{2}$, as each sub process, excitation and decay, is characterized by a rate, and the probability is therefore the product of $\left(W^{\text{excite}}T\right)\left(W^{\text{decay}}T\right)$. The substitution $\varepsilon_{v}\rightarrow\varepsilon_{v}-i\Gamma_{v}/2$ therefore corresponds to the replacement $W^{\text{decay}}T\rightarrow 1$ and then combines the processes corresponding to the first order diagrams of excitation first, and subsequently decay, with the second order diagram which allows for off-resonant excitation and decay. We are interested in the case when $T<1/\Gamma_{v}$ such that 2 photon emission is unlikely compared to 1 photon emission, and therefore higher number of photon emissions can be neglected. In this case one can also think of the previously obtained result for the cascade contribution, as the contribution of the finite crystal length to the line width, which corresponds to setting $\Gamma_{v}/2=1/T$, which will be the dominant contribution to the line width when $T\ll 1/\Gamma_{v}$.

V Discussion of results

In the figures in this paper we show the calculations made for a 180 GeV electron in the Doyle-Turner potential Baier et al. (1998); Doyle and Turner (1968); Avakian et al. (1982); Møller (1995) for the (110) planes in Silicon and for the state $n=25$. This is a quite low lying state which for electrons will have a high radiation power Wistisen and Di Piazza (2019). Electrons were chosen for this reason as it is not as numerically heavy when the quantum numbers are relatively small, as opposed to the positron case, which would require large quantum numbers to obtain an appreciable value of the quantum non-linearity parameter $\chi$, which means that quantum effects such as spin and recoil are important in the emission process. To compare with an experiment one should average over the distribution of the initial states which depends on the particle beam angular mean and divergence. In (41) the integrals over $\varphi$ and $\theta$ are carried out numerically over the intervals $0<\varphi<2\pi$ and $0<\theta<\frac{1.5}{\gamma}\times(1+\xi)$, and therefore includes nearly all emitted radiation. From the result of Eq. (41) we see that the part scaling with $T^{2}$ is the cascade, obtained by simple multiplication of probabilities, and will dominate unless the crystal is very thin, due to the remaining terms being proportional to $T$. Therefore, if one made a Monte Carlo approach using the single photon emission rate using the quantum numbers of the current state, instead of using the constant field approximation with the current value of the field, one would obtain the dominant (cascade) contribution, which will be accurate also when the constant field approximation is no longer valid. In figure (2) we show the result from the cascade process. In figure (3) we show the one-step terms and finally in figure (4) we show the ratio of these one-step terms to the cascade terms for $T=20\text{\textmu m}$. From this figure we see that the one-step terms can become significant compared to the cascade terms for short crystals. This ratio scales as $1/T$. Therefore one needs a thin crystal for the one-step contribution to be significant, so thin that the probability to emit more than 1 photon becomes small. One may rightfully ask based on these figures, if one picks a very small value of $T$, the total probability could seemingly become negative, however the results shown are only valid when $T\gg l_{f}\sim 0.8$ µm as estimated earlier. For the $180$ GeV case calculated here, the probability to emit a photon with energy above $1$ GeV from a 20 µm crystal is roughly $7\%$ and therefore the probability corresponding to the cascade for two-photon emission above this photon energy is $0.25\%$, and as can be seen in figure (4) the spectrum in the region where the radiation is most abundant, the ratio is around $\pm 20\%$. This number serves as an upper limit to the size of the effect, because under experimental conditions one would obtain the average from a population of many different levels with different quantum number $n$, and this averaging would likely reduce the size of the effect. If we assume the size of the effect to be this upper limit, one would need enough events such that one would have enough statistics to see an effect of such a size from only $0.25\%$ of the events. If this setup was realized by adding a calorimeter to a setup as the one used in Wistisen et al. (2018) we can estimate the number of particles required to see this.

Making a histogram of 20 bins in each direction of $\omega_{1}$ and $\omega_{2}$ and assuming 100 counts on average in each bin, one would need roughly $3.2\times 10^{8}$ electrons and assuming an electron rate of $10^{4}/\text{min}$ this translates into roughly 22 days of measuring time. This would therefore be a challenging experiment and having in mind that there would likely also be systematic uncertainties, the realistic outcome of such an experiment would be to put a constraint on the size of such one-step terms, rather than their direct observation.

VI Conclusion

In conclusion, we have shown how to accurately calculate the two photon emission rate for a high energy electron (or positron) channeled in a crystal. This calculation shows that the full probability contains what is known as the cascade, which could have been obtained multiplying probabilities of single photon emissions, as well as additional interference terms, called the one-step contribution. The one-step contribution scales only linearly with the crystal length, and therefore one needs a thin crystal to see the effect of these terms. We have calculated the size of all contributions to the emission probability for 180 GeV electrons in Silicon and found that with a long measuring time, the one-step contribution could possibly be seen. Since these effects are however small, we also see how to solve the problem of quantum radiation reaction, under general circumstances, in a crystal, by using the single photon emission rate in consecutive emissions, corresponding to the particle’s current state.

VII Acknowledgments

The author gratefully acknowledges useful discussions with Antonino Di Piazza and Karen Z. Hatsagortsyan. This work was partially supported by a research grant (VKR023371) from VILLUM FONDEN and later by the Alexander von Humboldt-Stiftung. In addition the author acknowledges the support of NVIDIA Corporation with the donation of the Titan V GPU used for this research.

Appendix A

The general (unnormalized) solution to the Dirac equation with potential energy $V(\boldsymbol{r})=-e\varphi(\boldsymbol{r})$ can be written as

[TABLE]

The Dirac equation then becomes

[TABLE]

[TABLE]

The electron solution is then

[TABLE]

We then obtained an equation for $\phi(\boldsymbol{r})$ by isolating $\chi(\boldsymbol{r})$ in Eq. (48) and inserting in Eq. (47). This solution has the property that it is well defined when $\varepsilon=m$. Another solution can be found by isolating $\chi(\boldsymbol{r})$ in Eq. (47) and inserting in (48). However this solution is not well defined when $\varepsilon=m$ and therefore one must use the negative energy solution $\varepsilon=-\sqrt{p^{2}+m^{2}}=-\varepsilon_{p}$ therefore we have

[TABLE]

where $\varepsilon_{p}$ is the positive energy of the positron. The equation for $\chi$ we can now be obtained by using

[TABLE]

which is equivalent with

[TABLE]

This is the same equation as the one we obtained for $\phi(\boldsymbol{r})$, except with the sign of $e$ changed such that, after making the same approximations as we did in Wistisen and Di Piazza (2018):

[TABLE]

We therefore make the ansatz in line with the usual approach (the sign on the momenta is changed):

[TABLE]

[TABLE]

Then

[TABLE]

with $\boldsymbol{p}=(p_{x}+q\varphi(\boldsymbol{r}),i\frac{d}{dy},p_{z})$, inserting $I^{+}(y)$, this becomes

[TABLE]

with $\boldsymbol{p}=(p_{x}+E_{n}-\frac{(jk_{0}+k_{B})^{2}}{2\varepsilon},jk_{0}+k_{B},p_{z})$.

Appendix B

The electron state can be written as (putting back in the volume factor)

[TABLE]

where

[TABLE]

where $\boldsymbol{p}_{j}=(p_{x}+E_{n}-\frac{(jk_{0}+k_{B})^{2}}{2\varepsilon},jk_{0}+k_{B},p_{z})$ and then

[TABLE]

Explicitly we have that $c_{j}=c_{j}(p_{x},k_{B},p_{z},n)$. Now since both $k_{B}$ and $k_{B}^{\prime}$ obey that $0\leq k_{B}<k_{0}$ we have that $-k_{0}<k_{B}-k_{B}^{\prime}<k_{0}$ and therefore $k_{B}-k_{B}^{\prime}$ can never be an integer value of $k_{0}$ unless $k_{B}-k_{B}^{\prime}=0$, and therefore we can write

[TABLE]

[TABLE]

However the vector $\boldsymbol{c}$ is a normalized ($|\boldsymbol{c}|=1$), eigenvector of a hermitian matrix and the vectors corresponding to $n$ and $n^{\prime}$ have different eigenvalues of this matrix, and are therefore orthogonal, so

[TABLE]

Now consider

[TABLE]

and therefore

[TABLE]

There $\simeq$ refers only to the normalization. The states are exactly orthogonal, but in the normalization we neglect corrections which are suppressed by at least $\xi/\gamma$ compared to leading order. So finally

[TABLE]

Appendix C

Even though we consider the radiation from electrons, the propagator contains terms from the positron $\psi_{p,n,s}^{+}(\boldsymbol{x}_{2})\bar{\psi}_{p,n,s}^{+}(\boldsymbol{x}_{1})$. Therefore we will need to calculate

[TABLE]

and $\mathcal{M}_{v\rightarrow f}^{+}$ so we need

[TABLE]

Then

[TABLE]

[TABLE]

Here we may use that $\boldsymbol{p}_{v}\times\left(\boldsymbol{\epsilon}^{*}\times\boldsymbol{p}_{i}\right)=\boldsymbol{\epsilon}^{*}(\boldsymbol{p}_{i}\cdot\boldsymbol{p}_{v})-\boldsymbol{p}_{i}(\boldsymbol{p}_{v}\cdot\boldsymbol{\epsilon}^{*})$ and so

[TABLE]

Now consider the other part for $\mathcal{M}_{v\rightarrow f}^{+}(\boldsymbol{k},\epsilon)$

[TABLE]

This is the same as before except with $i\rightarrow v$ and $v\rightarrow f$ . And now we want the quantity

[TABLE]

where now $n_{B,1}^{+}$ is chosen such that $k_{B,v}=k_{y}-k_{B,i}+n_{B,1}^{+}k_{0}$ is in the FBZ. Note that $-k_{B,v}=k_{B,i}-k_{y}-n_{B,1}^{+}k_{0}$ for which we already have the solution, called $k_{B,v}^{-}=k_{B,i}-k_{y}-n_{B,1}^{-}k_{0}$, and therefore

[TABLE]

therefore $n_{B,1}^{+}=n_{B,1}^{-}+1$. For the $\mathcal{M}_{v\rightarrow f}^{+}$ term one obtains that $k_{B,f}=-k_{B,v}-k_{y}-n_{B,2}^{+}k_{0}$ and for this term one has that $n_{B,2}^{+}=n_{B,2}^{-}-1$, in terms of the $n_{B,2}^{-}$ value for the corresponding electron term in the propagator. And that the $l$ index is given by $l=n_{B,2}^{+}-j$.

Appendix D

We need to consider $|\sum_{s_{v}}\mathcal{M}_{2}\mathcal{M}_{1}|^{2}$, in particular we would like to show that $\text{Re}\left(\left[\mathcal{M}_{2,\uparrow}\mathcal{M}_{1,\uparrow}\right]\left[\mathcal{M}_{2,\downarrow}\mathcal{M}_{1,\downarrow}\right]^{\dagger}\right)$ is 0, where the arrows denote the spin state of the virtual particle. This we may rearrange and consider therefore the product $\mathcal{M}_{2,\downarrow}^{\dagger}\mathcal{M}_{2,\uparrow}$. Now we may use that $\mathcal{M}$ can be written as

[TABLE]

Now for simplicity we define

[TABLE]

[TABLE]

and then we have that

[TABLE]

Therefore

[TABLE]

We assume that $\boldsymbol{\epsilon}^{*}=\boldsymbol{\epsilon}$, which is possible if we choose linear polarization as our basis, and we will perform the summation of final spins and therefore $\boldsymbol{s}_{f}\boldsymbol{s}_{f}^{\dagger}$ is the identity

[TABLE]

where we used that $\tilde{\boldsymbol{B}}$ is a real vector. For the other term, $\mathcal{M}_{1,\uparrow}\mathcal{M}_{1,\downarrow}^{\dagger}$ the same can be done, and here the argument hinges upon summation over initial spins, therefore, if either a summation is carried out over initial or final spins, the spin interference terms will be 0.