Sign of the terms in the summation of the second order rdts

[oscarxblanco]

Dear all,
during the review of the the function get_rdt that calculates the resonant driving terms in pyat (see #821) I saw a discrepancy in the sign of the terms inside the summation, Eq. (46) in Chun-xi Wang's article Explicit formulas for 2nd-order driving terms due to sextupoles and chromatic effects of quadrupoles, March 10, 2012, ANL/APS/LS-330 and the implementation coded in pyat, AT matlab and Elegant.

In the article
https://www.aps.anl.gov/files/APS-sync/lsnotes/files/APS_1429490.pdf

the pair $i,j$ are indexes to represent a pair of sextupoles along the ring where $j$ is strictly downstream wrt to $i$ if we use the Poisson bracket definition. This expression is rewritten by Wang as a sum over all $i,j$ combinations of a function $f(i,j)$ that he defines along his article for every resonant term. I remark that the terms in this summation are added when $j>i$ and substracted when $i>j$.

However, this summation is reversed in all code implementations (added when $j<i$ and substracted when $i<j$. See the python lines here below.

at/pyat/at/physics/rdt.py

Lines 228 to 229 in 3777663

	bsign = np.array([np.sign(sm[i] - sm) * 1j * b3lm[i] * b3lm
	for i in range(nelem)])

In the expression np.sign(sm[i] - sm) , sm contains the positions along the $s$-coordinate for each sextupole ordered from minimum to maximum, and the index i is the same index $i$ in Wang's equation. Then, consider the case when $i=1$, the sign expression should return zero for the first element in sm and positive for all others, as I wrote on the table above. However, np.sign(sm[i] - sm) will be negative for all elements except the first (I'm a little bit concerned about this floating point substraction to get the zero when $i=j$, but, it works).

i\j 1 2 3 ...
1   0 - - - -

You could say that maybe $i$ and $j$ are swapped, but, that is not the case. Compare, for example, $f(i,j)$ for $h_{31000}$ in the article and the implementation to be sure i is $i$ and all $j$'s are the elements in the numpy array.

at/pyat/at/physics/rdt.py

Lines 240 to 243 in 3777663

	rdts2["h31000"] += (1.0 / 32) * np.sum([
	rbbx[i] * pxm3[i] * cpxm
	for i in range(nelem)]
	)

the term pxm3[i] is the third power of $e$ to the horizontal phase advance multiplied by the imaginary number $i$, and cpxm is the complex conjugate which brings the minus in front of the phase for all $j$'s.

This is repeated in the C implementation RDTelegantAT.cpp used by the AT matlab function computeRDT

at/atmat/atphysics/NonLinearDynamics/RDTelegantAT.cpp

Line 328 in 3777663

termSign = SIGN(ed[iE].s - ed[jE].s);

at/atmat/atphysics/NonLinearDynamics/RDTelegantAT.cpp

Lines 334 to 336 in 3777663

	h31000 += (1./32)*termSign*ii*ed[iE].b3L*ed[jE].b3L*
	ed[iE].rbetax*ed[jE].rbetax*ed[iE].betax*ed[jE].betax*
	ed[iE].px[3]*conj(ed[jE].px[1]);

And this is also the case of Elegant, as it is reported in this link:

https://www3.aps.anl.gov/forums/elegant/viewtopic.php?f=17&t=1402

This is no coincidence. The code written in python was benchmarked using numerical results of the matlab+C implementation, which was at some point adapted to AT from Elegant. The Elegant developers reported in their forum (see the link above) that they benchmarked their code wrt OPA https://ados.web.psi.ch/opa/ , but, there is no further detail into which version of OPA, source code or the test.

We could assume that it is a mistake in Eq. (46) by Wang. I tried to reconstruct step by step his calculation and got the same expression, for me it is correct. I attach to this message some notes in a pdf focused on checking the sign of the terms in the summation. If you find yourself willing to spend some time reading it, please, do not hesitate to ask me for more details.
RDT2ndordersummsignElegant_addkets.pdf

There are aspects that I checked on Wang's citations and didn't find any error. But, maybe a look by a more expert reviewer on symplectic maps operations could shed some light. For example, the order of map operations used by Wang, or the direction on which the propagation of phase advance should be taken into account.

@swhite2401 mentioned that RDTs could be compared to particle tracking, and from literature their effect could be compared with real machine measurements.

I start this discussion in order to document the status of my findings on this topic. I'm on the way to go for the next steps: understand more on symplectic maps, and also try to obtain some evidence in tracking.

Best regards,
o