Automatic Width Computation with Small Mass Gap

Hi,

If the phase-space is so small for a given decay, then it is hard to believe that the narrow-width approximation still holds.
And therefore it is difficult to see a reason to not make the auto-width to crash.
(in other words it is crashing due to numerical precision, but it would make sense to make the code to crash for physical inconsistency between the benchmark and the hypothesis used to compute the width).

> My guess is that for some reason using mmaj = 1.925000e+01 MG fails to check that the PS is to tiny for the process N5 > N4 nu nu and runs into issues. If that is the case that would be a bug.

It is technically not a bug since
1) the paper clearly state that the code should not be used in such limit.
2) the code clearly warns the user that the code should not be used in such limit.

Actually, numerical issue is not something that we can really fully prevent.
It is true that one can try to replace the expression (a**2-b**2) by (a-b)(a+b)
Obviously this might fix the crash for your value and postponed the issue to smaller mass gap but it will still happens at some stage or just move the crash to another expression since this code was never designed in such limit.

But let me stress that this not only a numerical issue, this is a real physics issue. If NWA is not valid (which should be the case for such small mass gap) then they are no way to ensure that our algorithms is valid (i.e. it might discard real contribution, under/over estimate some decay channel/...).

In your mail, you also comment that you expect destuctive interference between two decay channel, this clearly shows that NWA is not valid and here I would need to go back to books in order to know how one can define the width.

Cheers,

Olivier

> On 7 Jun 2022, at 16:45, Johannes Rosskopp <email address hidden> wrote:
>
> New question #702099 on MadGraph5_aMC@NLO:
> https://answers.launchpad.net/mg5amcnlo/+question/702099
>
> Hi Olivier,
>
> There is a problem for the automatic computation of the decay width in MG version 2_9_10 as well as 3_4_0.
>
> With the model file and param card I sent to you already by Email.
> As written there, the mass splitting of the heavy neutrinos is very small.
> The process I'm using (although it should not be relevant for automatic width calculation) is:
> generate ud ud > ww, (ww > mu nn, (nn > mu ud ud))
>
> when setting the mass and yukawa to
> set mmaj = 1.925000e+01
> set yvn2 = 1.035916e-04
> set decay 8000018 = Auto
> set decay 8000020 = Auto
>
> the automatic width computation returns an math domain out of bounds error.
> I tracked the problem down to the function calculate_apx_psarea in the file decay_objects.py
> The problem is the part
>
> return math.sqrt((M ** 2+mass_n ** 2-M_eff_mean**2) ** 2-\
> (2*M *mass_n) ** 2)
>
> that takes the sqrt of smth negative (sqrt(-2.3283064365386963e-10)).
>
> It appears in the computation of the decay of a heavy_nu (N5) into the slighly lighter heavy_nu (N4) + 2 light nu's.
> The phase space for this is obviously very tiny and the function might have trouble with (numerics?).
> I would expect MG to just neglect this diagram and indeed using different parameter points (e.g. mmaj = 14 or mmaj = 20, rest as above) everything goes through fine and there is no decay of N5 > N4 nu nu.
>
> My guess is that for some reason using mmaj = 1.925000e+01 MG fails to check that the PS is to tiny for the process N5 > N4 nu nu and runs into issues. If that is the case that would be a bug.
>
>
> more info (paramters in calculate_apx_psarea that cause the error):
> M = 19.25000844903333
> M_eff_mean = 1.0302869668521453e-13
> mass_list = [0.0, 0.0]
> mass_n = 19.250008449033125
> self.calculate_apx_psarea(M_eff_mean, mass_list) = 0.019894367886486918
> (M ** 2+mass_n ** 2-M_eff_mean**2) ** 2-(2*M *mass_n) ** 2 = -2.3283064365386963e-10
>
>
> Please let me know if you need anything else.
>
> cheers
> Johannes
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi Olivier,

thanks for the fast reply. I'm a bit confused. The comment you mentioned about the interference is about the full process
 p p > w, w > N > mu mu j j.

But the automatic width computation should be independent of the process right? After all in the process I only consider the decay to a muon whereas in the width all the possible decays are considered.

I don't see where the NWA is used, as the decay is a 3-Body decay with an off-shell Z i.e. N5 > N4 v v
I thought the NWA is only used for N > 3 body decays or am I wrong?

cheers
Johannes

Hi,

One think is 100% clear to me is that the algorithm used in the automatic width is not designed for compressed spectra since it is likely missing some (potentially important) contribution. It is possible (likely?) that some contribution to the total width will not be include while they should be.

What is not clear to me is if the theoretical formula that we use for the computation of the width
(by using 1 to 2 and 1 to 3 body decay) is still valid or not. I know taht to use such formula we rely on the optic theorem.
But I need to check if the NWA is not also used there to pass from the pole formalism to 1>2 amplitude.
( I guess it does not apply if you have two particles so close in mass that they can technically mix when slightly off-shell)

So in short, I would not trust myself (today/now) to know how to compute such case, so I would certainly not trust the automatic width computation (even if it was returning a number)

Cheers,

Olivier

> On 7 Jun 2022, at 17:55, Johannes Rosskopp <email address hidden> wrote:
>
> Question #702099 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/702099
>
> Status: Answered => Open
>
> Johannes Rosskopp is still having a problem:
> Hi Olivier,
>
> thanks for the fast reply. I'm a bit confused. The comment you mentioned about the interference is about the full process
> p p > w, w > N > mu mu j j.
>
> But the automatic width computation should be independent of the process
> right? After all in the process I only consider the decay to a muon
> whereas in the width all the possible decays are considered.
>
> I don't see where the NWA is used, as the decay is a 3-Body decay with an off-shell Z i.e. N5 > N4 v v
> I thought the NWA is only used for N > 3 body decays or am I wrong?
>
> cheers
> Johannes
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi,

I discussed again with my colleges today and we are quite sure that it's not a physical problem but rather a numerical issue of MG.
After reading about the optical theorem etc. the only place where a NWA enters is when having a intermediate particle in a decay. However, that is not the typical NWA as I know it, since this intermediate particle has to be off-shell.
E.g. in the decay N5 > N4 Z*, Z* > v v, the Z is off shell.

Rather the propagator is replaced from
1/(p^2 - m^2 - i Im(Sig(p^2))) -> 1/(p^2 - m^2 + i m Gamma)

but this only effect the particles that are intermediate in the N body decays, which in our case are only the W and Z Boson and NOT the heavy neutrinos, so everything should be fine. Please correct me if I'm mistaken. I also don't see any other possible issues regarding the compressed spectrum.

In any case, to be sure I computed the decay width once for a big mass splitting (1Gev) and once for a very small one (10^-13 GeV), keeping the mass if the lighter neutrino and its couplings fixed. The result is that the widths are equal in both cases. This proves that the effect of N5 on the width of N4 is negligible (in case there is any effect at all which I doubt). Additionally even with the big mass splitting the decay N5 > N4 v v, does not show up as a branching ratio in the param card. This is because in the model it is strongly suppressed. So also the effect of N4 on the width of N5 is negligible (as expected).

Following this argument, I decided to just fix the numerical issue by hand, replacing Sqrt(x) by Sqrt(abs(x)).
It seems that everything works fine that way.

Thanks for taking the time to discuss the issue. I'll mark the Problem as solved but feel free to reply :)

Hi,

> Following this argument, I decided to just fix the numerical issue by hand, replacing Sqrt(x) by Sqrt(abs(x)).

You can obviously do all the change that you want to the code if that fit your situation.

This being said, this is a quite weird fix, that I would certainly not recomend.
If "x" can be negative this means that the computation of "x" is not reliable.
moving to abs(x) is obviously not reliable. It does avoid the crash but the associated value is still non-sense.

Now obviously if that channel does not significantly contributes to the width and that the code discard it anyway later then it should indeed not matter that "x" is a random number.

I discussed again with my colleges today and we are quite sure that it's not a physical problem but rather a numerical issue of MG.
After reading about the optical theorem etc. the only place where a NWA enters is when having a intermediate particle in a decay. However, that is not the typical NWA as I know it, since this intermediate particle has to be off-shell.
E.g. in the decay N5 > N4 Z*, Z* > v v, the Z is off shell.

As said before this is a topic that I do not master enough to comment on. My issue is not on Z (this is standard) but on the N5/N4 where some assomptotic state hyppothesis should also be taken which might not fully makes sense since the mass gap is so small that the two particles are actually kind of mixing. As I said, I would need to check the litterature to see how this is handle since I do not know.

Now even if they are no theoretical issue of using partial width, the selection algorithm is certainly not general to handle all the case. In particular the concept of open and close channel will not be valid anymore (or maybe it does but this is also not clear to me and if so why it does).

Now if it works for you --potentially with the above patch--, please go ahead.
I'm just sorry that I can not offer you more in this case.

Cheers,

Olivier

PS: Personally, I do not consider this as a bug since you are outside the regime where the code is designed.
Would be nice to see what can be done in such regime but this will require some large investment (likely too many for the associated impact)