2

I have a question about an Integrals which I can do by hand but I want to implement this Integral in Mathematica. Actually it is the definition of a 2 Particle Phase Space. The Expression is the following:

$$\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\theta (\text{E1}) \theta (\text{E2}) \delta \left(\text{E1}^2-\text{k1}^2-\text{m1}^2\right) \delta \left(\text{E2}^2-\text{k2}^2-\text{m2}^2\right) \delta (-\text{E1}-\text{E2}+\text{Ea}+\text{Eb}) \delta (-\text{k1}-\text{k2}+\text{ka}+\text{kb}) T(\text{E2},\text{k2},\text{E1},\text{k1})d\text{k2}d\text{E2}d\text{k1}d\text{E1}$$

When I try to integrate this with Mathematica it does nothing. (BTW T here is just a function e.g. the matrix element squared).I mean, when I do it by hand I treat e.g. the last delta dist. as a "function" of k1 and replace all the k1s in the other distributions and functions. So I kill step by step most of the deltas, at least two of them. But Mathematica does nothing. Can't Mathematica handle expressions like that? Or do I have to use some tricks? I would be very happy for some help:)

I mean, I know that it could be tricky for Mathematica to treat a specific delta as only a dist dependent on one specific argument for the specific integration.

Cheers,

Marcel

Marcel
  • 101
  • 4
  • $\theta$ is HeavisideTheta[]? – J. M.'s missing motivation Jun 18 '15 at 08:58
  • yep, sorry. Just copied from mathematica as LaTeX output:) – Marcel Jun 18 '15 at 09:04
  • I have run into this before and decided that the best way was not to use Integrate but my own Head, that behaves like Integrate when the integrand contains no DiracDeltas but does something smarter when it does exist. I think it's just a few simple replacement rules. – evanb Sep 16 '15 at 21:37

2 Answers2

1

Delta functions are equivalent to a system of equations. They have several solutions. I am not quite sure, but it seems that it is this that is the reason for Mma doubts. 

sl = Solve[{-k1 - k2 + ka + kb == 0, e1^2 - k1^2 - m1^2 == 0, 
    e2^2 - k2^2 - m2^2 == 0, -e1 - e2 + ea + eb == 0} , {k1, k2, e1, 
    e2}];

You could help Mma as follows substituting the first and then the second solution:

HeavisideTheta[e1]*HeavisideTheta[e2]*T[e2, k2, e1, k1] /. sl[[1]]

I cannot imagine, however, what you can do with the answer. Have fun!

Alexei Boulbitch
  • 39,397
  • 2
  • 47
  • 96
  • This is partly true, but it omits the transformation that needs to be done in order to get rid of the squared integration variables inside the delta functions, which introduces additional denominators. – Jens Jul 18 '15 at 17:03
0

I am sure there are more elegant ways to solve this, but one possibility might be the following module. As noted by Alexei Boulbitch it is necessary to transform the arguments of the delta functions as equations to be solved by Solve. The following implementation includes the necessary jacobian factors and allows in principle to start from the complete integrand. All integrations (which are evaluated through delta functions) are assumed to run from minus to plus infinity. Integration bounds can in principle included through corresponding HeavisideTheta-functions. The current implementation only applies to a scalar integrand (an extension to vector integrals shouldn't be too hard). The result must be then used as further input for Integrate or NIntegrate to evaluate the remaining integrations. Note that the module makes use of the function removeTimes defined at Convert head Times to List

deltaevalv[expr_, intlist_] := 
Module[{tmp, deltaprod, arglist, arglisteq0, jac, resS}, 
tmp = expr /. delta[a_] :> 1; deltaprod = expr/tmp; 
arglist = removeTimes[deltaprod] /. delta[a_] :> a; 
jac = Abs[
 Det[Table[
   D[arglist[[i]], intlist[[j]]], {i, 1, 
    Dimensions[intlist][[1]]}, {j, 1, 
    Dimensions[intlist][[1]]}]]]^-1; 
arglisteq0 = 
Table[arglist[[i]] == 0, {i, 1, Dimensions[intlist][[1]]}]; 
resS = Solve[arglisteq0, intlist]; 
Sum[jac *tmp /. resS[[j]], {j, 1, Dimensions[resS][[1]]}]]

The arguments are defined as follows: expr_ denotes the complete integrand which includes delta functions \delta(a) as delta[a] and intlist_ are the variables which should be evaluated through the delta functions. It is implicitly assumed that the number of intlist agrees with the number of delta functions. An example would be

deltaevalv[f[x1] g[x2] h[x3] delta[x1^2 - m12] delta[x2^2 - m22] delta[
x2 + x1 - x3 - 4], {x1, x2, x3}]

which yields the result

(f[-Sqrt[m12]]*g[-Sqrt[m22]]*h[-4 - Sqrt[m12] - Sqrt[m22]])/
(4*Abs[Sqrt[m12]*Sqrt[m22]]) + 
(f[Sqrt[m12]]*g[-Sqrt[m22]]*h[-4 + Sqrt[m12] - Sqrt[m22]])/
(4*Abs[Sqrt[m12]*Sqrt[m22]]) + 
(f[-Sqrt[m12]]*g[Sqrt[m22]]*h[-4 - Sqrt[m12] + Sqrt[m22]])/
(4*Abs[Sqrt[m12]*Sqrt[m22]]) + 
(f[Sqrt[m12]]*g[Sqrt[m22]]*h[-4 + Sqrt[m12] + Sqrt[m22]])/
(4*Abs[Sqrt[m12]*Sqrt[m22]])
Community
  • 1
MM8778
  • 56
  • 3