Number of active nf in Singlet FONLL about yadism

[giacomomagni]

On a similar point for the heavy SF can you confirm that also here:

yadism/src/yadism/coefficient_functions/kernels.py

Lines 315 to 325 in 46ca667

	ns_partons[ihq] = w
	ns_partons[-ihq] = w if kind != "F3" else -w
	ch_av = w / (nf) if kind != "F3" else 0.0
	for pid in range(1, nf):
	s_partons[pid] = ch_av
	s_partons[-pid] = ch_av
	return (
	Kernel(ns_partons, light_cfs.NonSinglet(esf, nf)),
	Kernel({21: ch_av}, light_cfs.Gluon(esf, nf)),
	Kernel(s_partons, light_cfs.Singlet(esf, nf)),
	)

we have the desired behaviour ?

I.e the nf-th parton is not coupling to the boson in singlet diagrams but singlet and gluon coefficients are called with nf.

And not:

 ns_partons[ihq] = w 
 ns_partons[-ihq] = w if kind != "F3" else -w 
 ch_av = w / (nf) if kind != "F3" else 0.0 
 for pid in range(1, nf + 1): 
     s_partons[pid] = ch_av 
     s_partons[-pid] = ch_av 
 return ( 
     Kernel(ns_partons, light_cfs.NonSinglet(esf, nf)), 
     Kernel({21: ch_av}, light_cfs.Gluon(esf, nf)), 
     Kernel(s_partons, light_cfs.Singlet(esf, nf)), 
 ) 

[alecandido]

When skip_heavylight is used to generate the nc_weights the heavy parton is not coupling to the boson
as correctly said, but still the singlet (and the gluon) charges are divided by the number of light flavors plus the heavy quark nl+1.

I'm not sure if I'm fully understanding, all this business would have required proper nomenclature and notation.
But if I remember correctly, the meaning should be the following:

• when you use skip_heavylight you should be dropping the contributions of the former heavy flavor (so the nl + 1 in the lower patch, that in the higher patch it would be just light as the others, and it is called just $h$ and $\bar{h}$ in the FONLL paper)
• this should be done to compute e.g. bottom contributions when FONLL is used for charm
• the idea is that bottom pretends that charm does not exist, in order to decouple the two mass effects (including only one kind of contributions would be inconsistent, you need to compute bottom mass effects in a consistent scheme, and thus you just pretend charm not to exist)

However, this is part of the very complex business of that infamous diagram I posted many times. skip_heavylight in some of the fancy terms in 3 or 4 flavors. But from the diagram I do not actually see which (above I was mainly describing Fheavy(3f:b), but that is in 3 flavors, where charm is actually heavy...).

                             ─────────────────────────────────────────────────────────
                              mc: FONLL    mb: FONLL    mt: ZM-VFNS    intrinsic: [4]
                                                Observable: Ftotal
                             ─────────────────────────────────────────────────────────


                                                ┌───────────────────┐
                                                │ active flavors: 3 │
                                                └───────────────────┘
                                                Flight = Flight(zm:3) + Fmiss(3f:c, /c, ~~3) + Fmiss(3f:b, /b, ~~3)
                                                Fcharm = Fheavy(3f:c) + Fheavy(3f:c, I4) + Fmiss(3f:b, /b, ~~c)
                                                Fbottom = Fheavy(3f:b)
                                                Ftotal = Flight + Fcharm + Fbottom
                           ┌───────────────────┐
                           │ active flavors: 4 │
                           └───────────────────┘
    Flight = Flight(FO:c) + Fmiss(4f:b, /b, ~~4)
Fcharm = Fheavy(FO:c, I4) + Fmiss(4f:b, /b, ~~c)
                          Fbottom = Fheavy(4f:b)
              Ftotal = Flight + Fcharm + Fbottom
                                                ┌───────────────────┐
                                                │ active flavors: 5 │
                                                └───────────────────┘
                                                Flight = Flight(FO:b)
                                                Fbottom = Fheavy(FO:b)
                                                Ftotal = Flight + Fbottom
                           ┌───────────────────┐
                           │ active flavors: 6 │
                           └───────────────────┘
                           Flight = Flight(zm:6)
                                 Ftotal = Flight

In conclusion, I actually do not see where something like skip_heavilight would make sense.
Brief legend:

• 3f:b is "3 light flavors + a massive b"
• zm:6 is "6 light flavors"
• FO:c is "FONLL for charm quark"
• ~~x means that x is coupling
• I4 includes intrinsic contributions for charm
• /b is special for "Fmissing", and states that the missing diagram is associated to a massive bottom (even though it seems redundant with the FNS indication...)

However, some time is passed, and it was also confusing at the time we worked out, so I should spend some time to understand from scratch.
My current proposal is to supersede all this business with the so-called "numerical FONLL". I'd rather call it "straightforward/basic FONLL", "FONLL at the level of observables" instead of coeff funcs, or just the "original FONLL" (since it's what they define before working out the very complicated coeff funcs). They are both computed numerically, and that's just the computing tool, not the abstract algorithm. But in the end is just a name (even if I'd choose a telling one for a paper).

In that case, there won't be any missing, skip_heavylight or anything like that. Just simple coefficient functions computed in a very well-defined FFNS (the most exotic thing will be asy, that's just a limit).

[giacomomagni]

Thanks @alecandido for the summary. My questions are a bit more specific. Take EM for simplicity, no I4.
So starting from your syntax the questions I'm wondering are:

1. skip_heavylight is only used in generate_light_diff. Take Flight(FO:c), nl=3 as an example: here you use skip_heavylight when you need to include the contributions of the charm initiated diagrams in which the boson is coupled to a light quark (ie Singlet diagrams). Now you remove the charm from the couplings using skip_heavylight: the averaged coupling is done summing on u,d,s, but dividing by nf=nl+1. Then you call the coefficient function Singlet with nl+1.
   Is this what you actually want?
   I thought you should remove charm from couplings; build the averaged charge summing on u,s,d and divide by nl and then calling Singlet with nl.
   As said the 2 methods are identical up to NNLO where the nf dependency factors out S(x,nf) = nf * S(x) and no matter what you do you will end up with a coupling summed on u,d,s times S(x). However at N3LO this is not the case anymore since the nf dependency do not factors out completely from the singlet cf so S(x,nf) != nf * S(x) and the 2 methods can lead to different results

2. The second point is related to generate_single_flavor_light and in particular when is called in generate_heavy_diff, so
   Fheavy(FO:c), nl=3.
   Is it correct that the heavy quark is included in the NS diagrams but is not included in the Singlet couplings although both light_cfs.Gluon and light_cfs.Singlet are called with nf = nl +1 ?
   

   yadism/src/yadism/coefficient_functions/kernels.py

   Lines 315 to 325 in 46ca667

   	ns_partons[ihq] = w
   	ns_partons[-ihq] = w if kind != "F3" else -w
   	ch_av = w / (nf) if kind != "F3" else 0.0
   	for pid in range(1, nf):
   	s_partons[pid] = ch_av
   	s_partons[-pid] = ch_av
   	return (
   	Kernel(ns_partons, light_cfs.NonSinglet(esf, nf)),
   	Kernel({21: ch_av}, light_cfs.Gluon(esf, nf)),
   	Kernel(s_partons, light_cfs.Singlet(esf, nf)),
   	)

   
   I thought here you should include the heavy quark also in the Singlet couplings so having range(1, nf+1) since you have ´nl+1´ light flavors.

Hope this explain my doubts a bit better... Thanks for any clarification.

[alecandido]

Ok, concerning your point 1., the additional context helped: I still don't remember exactly why we do it, but I can tell you that if you put something in the coeff func call, it has to be nl + 1.

The reason is that you want to remove the contribution of the coupling charm, but only the one coupling and nothing else.
So, all the possible loops have to be computed with nl + 1 flavors, since even the charm is running when decoupled from the EW boson.
Regarding the coupling itself: I believe you're just subtracting a contribution, since it will be accounted for somewhere else. So you do not want to change the average, you just want to remove charm from the sum: everything else has to remain the same of the skip_heavylight=False case.

In practice, I believe that skip_heavylight=True should have the effect of setting the electric charge of the former heavy quark (in the case we're discussing, the charm quark) to 0. But only the electric charge, so charm is still running, i.e. it has the usual strong coupling and mass that you would use in 4FNS.

In short: you should not change your average, because you're not dropping charm as you would do to change FNS, just subtracting a set of diagrams to put them somewhere else.

So, yes, I believe what is done is consistent, simply cumbersome.

And now I also remember why: you're computing the contributions to Flight, but Flight defined in the paper is defined based on electric charges (i.e. a la APFEL).
Our usual definition of Flight in yadism is no heavy quark ever. But this can't be for Flight FONLL, since without heavy quarks there is no FONLL, so we restored APFEL definition to compute Ftotal = Flight + sum(Fheavy for heavy in heavy_quarks), you can see it from the diagram.

[alecandido]

Concerning point 2.: the average in Fheavy is correct for the same reason. The singlet contribution is defined as the one in which the incoming quark does not couple directly to the EW photon (so the quark line incoming and the quark line coupling are not the same, but they are separated by a gluon exchange).

But the singlet and gluon couplings are weird, since they contain an average by definition. In principle, it should be a sum, not an average, because it's the sum over diagrams. And that's why it does not matter at NLO: it is fully artificial, since there are no quark loops, the only quark line is free, because decoupled by the incoming parton, but weighted with the electric charge. So the two dependencies cancel (average in the coupling, number of flavors in coeff funcs).
This is just a result of Vogt definition.

So, wherever you use it, you have to use the average with the number of running flavors, so all of them, because it is actually counteracted by the number of flavors dependency artificially introduced in the coeff func.

I'm still thinking at why we're not including charm in the incoming line.
It is most likely correct: incoming lines are pretty simple, there is mostly no complex behavior there (apart for intrinsic), so intuitively I'd believe the same of you.

However, since they are simple, it would be a rather gross mistake, and I know that we and APFEL are doing the same, because we obtained a perfect benchmark of FONLL implementation (in the comparable region, of course). So, it would be a repeated gross mistake, done by three different people (i.e. me + @felixhekhorn, and Valerio).
For this reason, I have a larger prior to consider it correct, but I'm still missing the why. Any idea @felixhekhorn?

[alecandido]

P.S.: on one side, I'd really love to get rid of all these complications with "numerical FONLL", that's why I'm insisting with @andreab1997 that we should not implement it side-by-side to the current one, but we should use the current version to dump grids, scratch everything from FONLL implementation, and implement the new one in a clean environment, benchmarking with the grids.
This version got sufficiently complicated on its own, so I would not complicate it even more with an additional implementation (we don't have the manpower to maintain both), and just focus on a single, simpler one.

However, we still have to write a yadism paper, and of course all these considerations will fit inside (as we agreed with @felixhekhorn long ago, yadism paper will be an FNSs paper). So it's good to review all these details, and also good to have them written somewhere (even though it should have been in the docs, and the more technical part in the docstrings :P).

[giacomomagni]

Okay, thanks for the detailed explanations:

1. concerning this point, I now understand the logic and I'll fix DIS_TP accordingly.
2. I agree with you that this will make quite a difference, and also in DIS_TP we don't include the initial heavy quark line contribution in the singlet coupling when nf == ihq although singlet cf is called with nf also there. So I'd say this is a feature more than a problem, but I couldn't remember why...
   I'm sure @felixhekhorn or @andreab1997 know it.

So this issue is not a bug anymore, once we recall also that we can safely close.

[giacomomagni]

completed in #183