3

𝑃1 scattering

D. R. Entema , J. A. Ollerb , and M. S. Sáncheza

a Grupo de Física Nuclear and IUFFyM, Universidad de Salamanca, E-37008 Salamanca, Spain
b Departamento de Física, Universidad de Murcia, E-30071 Murcia, Spain

June 2, 2025

This is our abstract.

1 Introduction
ind This is our introduction.

2 Schrödinger equation
ind The (uncoupled) 𝐿-wave radial Schrödinger equation given a central two-nucleon (𝑁 𝑁) potential 𝑉 (𝑟), with 𝑟 the
inter-nucleon radial distance, reads
 2
− 𝜕𝑟 + 𝑚 𝑁 𝑉 (𝑟) + 𝐿 (𝐿 + 1)/𝑟 2 − 𝑘 2 𝑢 𝐿 (𝑟; 𝑘) = 0 ,

(1)

where 𝑚 𝑁 = 939 MeV is the nucleon mass and 𝑢 𝐿 (𝑟; 𝑘) is the 𝐿-wave reduced wave function, 𝑘 standing for the
magnitude of the 𝑁 𝑁 center-of-mass momentum. We consider the one-pion exchange (OPE) interaction

𝑚 3𝜋
V𝜋 (®
𝑟) = 𝝉1 · 𝝉2 [ 𝜎
®1 · 𝜎
® 2 Y (𝑚 𝜋 𝑟) + 𝑆12 (𝑟)
ˆ T (𝑚 𝜋 𝑟)] , (2)
3 𝑚 𝑁 𝑀𝑁 𝑁
where 𝑚 𝜋 = 138 MeV is the pion mass, which sets the range of the interaction. In turn, its strength is set by the
low-momentum scale [1, 2] 𝑀 𝑁 𝑁 = (4𝜋/𝑚 𝑁 ) (2 𝑓 𝜋 /𝑔 𝐴) 2 ≈ 290 MeV, with 𝑓 𝜋 = 92.4 MeV and 𝑔 𝐴 = 1.26 the pion
decay and nucleon axial-coupling constants, respectively. Besides,
𝑒 −𝑥
Y (𝑥) = ; (3)
𝑥 
𝑒−𝑥

3 3
T (𝑥) = 1+ + 2 (4)
𝑥 𝑥 𝑥
are dimensionless functions, and
  
3 −3 𝑆 = 0,
®1 · 𝜎
𝜎 ®2 = 2 𝑆(𝑆 + 1) − = (5)
2 +1 𝑆 = 1;
  
3 −3 𝑇 = 0 ,
𝝉1 · 𝝉2 = 2 𝑇 (𝑇 + 1) − = (6)
2 +1 𝑇 = 1

are spin and isospin factors. Since nucleons are fermions, 𝐿 + 𝑆 + 𝑇 has to be an odd number. Finally,

𝑆12 (𝑟)
ˆ = 3 (𝜎 ˆ (𝜎
® 1 · 𝑟) ® 2 · 𝑟)
ˆ −𝜎
®1 · 𝜎
®2 (7)

is the tensor operator. This vanishes in the scalar (singlet) waves, so that 𝑉 (𝑟) collapses to a regular (Yukawa-type)
potential, but induces the interaction in the vector (triplet) waves to diverge as ±1/𝑟 3 at short distances, where the sign
corresponds to that of 𝑆12 (𝑟)
ˆ in a particular wave.

1
 Let us assume that 𝑘 2 ⩾ 0 (scattering case). The potential becomes negligible at large distances, 𝑉 (𝑟 → ∞) → 0,
and Eq. (1) has the general free solution
h i  𝜋  𝜋
𝑢 𝐿 (𝑟 → ∞; 𝑘) = 𝛼 𝐿 (𝑘) 𝚥ˆ𝐿 (𝑘𝑟) + 𝛽 𝐿 (𝑘) 𝑦ˆ 𝐿 (𝑘𝑟) = 𝛼 𝐿 (𝑘) sin 𝑘𝑟 − 𝐿 − 𝛽 𝐿 (𝑘) cos 𝑘𝑟 − 𝐿 , (8)
𝑟→∞ 2 2
where 𝛼 𝐿 (𝑘) and 𝛽 𝐿 (𝑘) are arbitrary functions of the momentum, and 𝚥ˆ𝐿 (𝑥) ≡ 𝑥 𝑗 𝐿 (𝑥) and 𝑦ˆ 𝐿 (𝑥) ≡ 𝑥 𝑦 𝐿 (𝑥) are
(𝐿th-order) reduced first- and second-kind spherical Bessel functions with well-known asymptotic series and behavior
near the origin,
 𝜋  𝜋
𝚥ˆ𝐿 (𝑥 → ∞) = sin 𝑥 − 𝐿 + O (𝑥 −1 ) , 𝑦ˆ 𝐿 (𝑥 → ∞) = − cos 𝑥 − 𝐿 + O (𝑥 −1 ) ; (9)
2 2
𝑥 𝐿+1 (2𝐿 − 1)!! 
1 + O (𝑥 2 ) , 1 + O (𝑥 2 ) .
  
𝚥ˆ𝐿 (𝑥 → 0) = 𝑦ˆ 𝐿 (𝑥 → 0) = − (10)
(2𝐿 + 1)!! 𝑥𝐿
Besides, they verify the general property

𝚥ˆ𝐿 (𝑥) 𝑦ˆ 𝐿′ (𝑥) − 𝑦ˆ 𝐿 (𝑥) 𝚥ˆ𝐿′ (𝑥) = 1. (11)

One thus finds

 
𝑢 𝐿 (𝑟; 𝑘) 𝜕𝑟 𝚥ˆ𝐿 (𝑘𝑟) − 𝚥ˆ𝐿 (𝑘𝑟) 𝜕𝑟 𝑢 𝐿 (𝑟; 𝑘) 𝛽 𝐿 (𝑘)
= − . (12)
𝑢 𝐿 (𝑟; 𝑘) 𝜕𝑟 𝑦ˆ 𝐿 (𝑘𝑟) − 𝑦ˆ 𝐿 (𝑘𝑟) 𝜕𝑟 𝑢 𝐿 (𝑟; 𝑘) 𝑟→∞ 𝛼 𝐿 (𝑘)

It is customary to take 𝛼 𝐿 (𝑘) = cot 𝛿 𝐿 (𝑘) and 𝛽 𝐿 (𝑘) = −1, 𝛿 𝐿 (𝑘) being the 𝐿-wave phase shift. 1 The latter is thus
obtained from the large-distance reduced wave function as
 
−1 𝜕𝑟 𝑦ˆ 𝐿 (𝑘𝑟) − 𝑦ˆ 𝐿 (𝑘𝑟) L 𝐿 (𝑟; 𝑘)
𝛿 𝐿 (𝑘) = cot lim , (13)
𝑟→∞ 𝜕𝑟 𝚥ˆ𝐿 (𝑘𝑟) − 𝚥ˆ𝐿 (𝑘𝑟) L 𝐿 (𝑟; 𝑘)

where the logarithmic derivative

L 𝐿 (𝑟; 𝑘) ≡ 𝑢 𝐿 (𝑟; 𝑘) −1 𝜕𝑟 𝑢 𝐿 (𝑟; 𝑘) (14)

fulfills the first-order inhomogeneous differential equation

𝜕𝑟 L 𝐿 (𝑟; 𝑘) + L 𝐿 (𝑟; 𝑘) 2 = 𝑚 𝑁 𝑉 (𝑟) + 𝐿 (𝐿 + 1) /𝑟 2 − 𝑘 2 , (15)

as found from Eq. (1). To solve Eq. (15), one must impose one initial condition. It is customary to introduce an
ultra-violet momentum cutoff Λ ≫ 𝑘, thus a cutoff radius 𝑟 𝑐 ≡ 1/Λ such that the potential vanishes at short distances
in virtue of the replacement

𝑉 (𝑟) → 𝑉 (𝑟) 𝜃 (𝑟 − Λ1 ) . (16)

For the reduced second-kind spherical Bessel function does not tend to zero around the origin [cf. Eq. (10)], Eq. (16)
leads to the initial condition
𝜕𝑟 𝚥ˆ𝐿 (𝑘𝑟)
L 𝐿 ( Λ1 ; 𝑘) = (𝐿 + 1) Λ 1 + O (𝑘 2 /Λ2 ) ,
 
= (17)
𝚥ˆ𝐿 (𝑘𝑟) 𝑟= Λ1

which is actually energy-independent in the cutoff-removal limit Λ → ∞. Besides, one has the low-energy expansion

L 𝐿 (𝑟; 𝑘) = L 𝐿(0) (𝑟) + L 𝐿(2) (𝑟) 𝑘 2 + L 𝐿(4) (𝑟) 𝑘 4 + O (𝑘 6 ) , (18)

where the first three coefficients can be found successively by solving the differential equations

𝜕𝑟 L 𝐿(0) (𝑟) + L 𝐿(0) (𝑟) 2 = 𝑚 𝑁 𝑉 (𝑟) + 𝐿 (𝐿 + 1) /𝑟 2 ; (19)

𝜕𝑟 L 𝐿(2) (𝑟) + 2 L 𝐿(0) (𝑟) L 𝐿(2) (𝑟) +1 = 0; (20)

1 The solution is thus normalized so that its asymptotic behavior (8) yields 𝑢𝐿 (𝑟 → ∞; 𝑘 ) = sin[ 𝑘𝑟 − 𝐿 𝜋/2 + 𝛿 𝐿 (𝑘 ) ]/sin 𝛿 𝐿 (𝑘 ).

2
 𝜕𝑟 L 𝐿(4) (𝑟) + 2 L 𝐿(0) (𝑟) L 𝐿(4) (𝑟) + L 𝐿(2) (𝑟) 2 = 0, (21)

supplemented by the initial conditions

L 𝐿(0) ( Λ1 ) = (𝐿 + 1) Λ ; L 𝐿(2) ( Λ1 ) = 0 ; L 𝐿(4) ( Λ1 ) = 0 . (22)

On the other hand, the 𝐿-wave effective-range expansion is

1 𝑟𝐿 2
𝑘 2𝐿+1 cot 𝛿 𝐿 (𝑘) = − + 𝑘 + 𝑣 𝐿 𝑘 4 + O (𝑘 6 ) . (23)
𝑎𝐿 2
Given 𝐿 = 0, 1, . . . , 𝑎 𝐿 is the so-called scattering length, scattering volume, . . . ; 𝑟 𝐿 is the effective range, effective
momentum, . . . ; the first shape parameter 𝑣 𝐿 ≡ 𝑣 2,𝐿 has mass dimension 2𝐿 − 3. Then, expanding Eq. (13) around
𝑘 = 0 and matching yields

𝑟 2𝐿+1  1 1 
𝑎𝐿 = lim 2
− (0)
; (24)
𝑟→∞ (2𝐿 − 1)!! 1 + 2𝐿 𝐿 + 𝜆 (𝑟)
𝐿
ℭ 𝐿(1) ℭ 𝐿(2)
𝑟𝐿 = ℭ 𝐿(0) + + 2 ; (25)
𝑎𝐿 𝑎𝐿
𝔇 𝐿(1) 𝔇 𝐿(2) 𝔇 𝐿(3)
𝑣 2,𝐿 = 𝔇 𝐿(0) + + 2 + 3 , (26)
𝑎𝐿 𝑎𝐿 𝑎𝐿

where we introduced the dimensionless shorthand (𝑛 = 0, 1, 2)

𝜆 𝐿(2𝑛) (𝑟) = 𝑟 1−2𝑛 L 𝐿(2𝑛) (𝑟) , (27)

together with

2 (2𝐿 − 1)!!2  (2) 1 

ℭ 𝐿(0) = lim 𝜆 (𝑟) + , (28)
𝑟→∞ 𝑟 2𝐿−1 𝐿 1 − 2𝐿
4𝑟 2  (2) 1
ℭ 𝐿(1) = − lim 𝜆 𝐿 (𝑟) + , (29)
𝑟→∞ 1 + 2𝐿 2
2𝑟 3+2𝐿  (2) 1 
ℭ 𝐿(2) = lim 𝜆 (𝑟) + ; (30)
𝑟→∞ (2𝐿 + 1)!!2 𝐿 3 + 2𝐿
 𝜆 𝐿(2) (𝑟) 2 + 𝜆 𝐿(2) (𝑟) − 2𝐿−3
1 
𝔇 𝐿(0) = lim (2𝐿 − 1)!! 𝑟 2 3−2𝐿 (4)
𝜆 𝐿 (𝑟) + , (31)
𝑟→∞ 1 + 2𝐿
2𝑟 4  (4) 3𝜆 (2) (𝑟) 2 (2 + 𝐿) [1 + 4𝜆 𝐿(2) (𝑟)] 
𝔇 𝐿(1) = − lim 𝜆 𝐿 (𝑟) + 𝐿 + , (32)
𝑟→∞ 1 + 2𝐿 2 + 4𝐿 (2 + 4𝐿) (3 + 2𝐿)
𝑟 5+2𝐿  (4) 3𝜆 𝐿(2) (𝑟) 2 + 7+2𝐿 (2) 1
3+2𝐿 [𝜆 𝐿 (𝑟) + 5+2𝐿 ]

𝔇 𝐿(2) = lim 𝜆 (𝑟) + , (33)
𝑟→∞ (2𝐿 + 1)!!2 𝐿 1 + 2𝐿
𝑟 6+4𝐿  (2) 1 2
𝔇 𝐿(3) = − lim 𝜆 (𝑟) + . (34)
𝑟→∞ (2𝐿 + 1)!!4 𝐿 3 + 2𝐿

For the particular case of the 3 𝑃1 channel, 𝐿 = 𝑆 = 𝑇 = 1. Besides, 𝐽 = 𝐿 so this is indeed an uncoupled wave
with 𝑆12 = 2. Hence, OPE (2) becomes

𝑚 3𝜋
𝑉 (𝑟) = [Y (𝑚 𝜋 𝑟) + 2T (𝑚 𝜋 𝑟)] , (35)
3𝑚 𝑁 𝑀 𝑁 𝑁
which is singular and repulsive. This latter feature is what prevents the collapse rendering an ill-defined solution one
encounters when cutoffs Λ ≳ 1 GeV are applied to waves where OPE is singular and attractive and no counterterm is
included [3,4], such as the 3 𝑃0 wave. Indeed, rather quick cutoff convergence is observed for the 3 𝑃1 wave. Conversely,
when one or plus counterterms are added the phase shifts are not well-defined anymore.

3
3 Lippmann–Schwinger equation
ind Now we consider the partial-wave projected Lippmann–Schwinger equation without mixing (𝐿 = 𝐿 ′ ), provided
with a (local) sharp-cutoff regulator
∫ Λ
′ ′ 𝑚𝑁 𝑣 𝐿 ( 𝑝 ′ , 𝑞) 𝑞 2 𝑡 𝐿 (𝑞, 𝑝; 𝑘)
𝑡 𝐿 ( 𝑝 , 𝑝; 𝑘) = 𝑣 𝐿 ( 𝑝 , 𝑝) + 𝑑𝑞 , (36)
2𝜋 2 0 𝑘 2 − 𝑞 2 + 𝑖 0+
with the partial-wave projected interaction given by the Fourier transform
∫ ∞
′ 4𝜋
𝑣 𝐿 ( 𝑝 , 𝑝) = 𝑑𝑟 𝚥ˆ𝐿 ( 𝑝 ′ 𝑟) 𝑉 (𝑟) 𝚥ˆ𝐿 ( 𝑝𝑟) . (37)
𝑝′ 𝑝 0
Once the off-shell amplitude 𝑡 𝐿 ( 𝑝 ′ , 𝑝; 𝑘) is found, the phase shift is obtained as
 
𝑖 𝑖 2
𝛿 𝐿 (𝑘) = − log 1 − 𝑚 𝑁 𝑘 𝑇𝐿 (𝑘 ) , (38)
2 2𝜋
with the on-shell amplitude

𝑇𝐿 (𝑘 2 ) ≡ 𝑡 𝐿 (𝑘, 𝑘; 𝑘) . (39)

Specializing for the 3 𝑃1 interaction (35), Eq. (37) yields

 ′2
𝑝 + 𝑝 2 − 𝑚 2𝜋 ( 𝑝 ′ 2 − 𝑝 2 ) 2 − 𝑚 4𝜋
 ′
( 𝑝 + 𝑝) 2 + 𝑚 2𝜋

4𝜋
𝑣 1 ( 𝑝 ′ , 𝑝) = − log . (40)
𝑚 𝑁 𝑀𝑁 𝑁 4𝑝 ′ 𝑝 16𝑝 ′ 2 𝑝 2 ( 𝑝 ′ − 𝑝) 2 + 𝑚 2𝜋
Again, the repulsive nature of this interaction allows one to safely take Λ → ∞ in Eq. (36), so that 𝑡 1 ( 𝑝 ′ , 𝑝; 𝑘), thus
𝛿1 (𝑘), will be well convergent.

4 Unsubtracted N/D method

ind Consider an 𝑁/𝐷 decomposition of the 𝐿-wave projected on-shell amplitude
𝑁 𝐿 ( 𝐴)
𝑇𝐿 ( 𝐴) = , (41)
𝐷 𝐿 ( 𝐴)

where 𝑁 𝐿 ( 𝐴) has only a left-hand cut (LHC) for 𝐴 < −𝑚 2𝜋 /4 ≡ 𝐿, and 𝐷 𝐿 ( 𝐴) has only a unitarity or right-hand cut
(RHC) for 𝐴 > 0. Namely,

ℑ [𝑁 𝐿 ( 𝐴)] = 𝐷 𝐿 ( 𝐴) Δ 𝐿 ( 𝐴) 𝜃 (𝐿 − 𝐴) (42)

in terms of the amplitude discontinuity defined for 𝐴 ∈ LHC

Δ 𝐿 ( 𝐴) ≡ (2𝑖) −1 [𝑇𝐿 ( 𝐴 + 𝑖0+ ) − 𝑇𝐿 ( 𝐴 − 𝑖0+ )] , (43)

and
𝑚𝑁 √
ℑ [𝐷 𝐿 ( 𝐴)] = − 𝐴 𝑁 𝐿 ( 𝐴) 𝜃 ( 𝐴) , (44)
4𝜋
respectively. In this approach, the functions 𝑁 𝐿 ( 𝐴) and 𝐷 𝐿 ( 𝐴) are obtained through solving a system of two integral
equations (IEs) provided with 𝔫 and 𝔪 subtractions, respectively — this is what we call the 𝑁 𝐷 𝔫𝔪 equations.
Since there is always the ambiguity of an overall factor in the definitions of 𝑁 𝐿 ( 𝐴) and 𝐷 𝐿 ( 𝐴) [cf. Eq. (41)], the
normalization must be fixed by hand. This translates into the fact that our IEs must subsume at least one subtraction.
Most commonly, one adopts the zero-energy prescription

𝐷 𝐿 (0) = 1, (45)

so that around threshold

𝑁 𝐿 ( 𝐴 → 0) ∝ 𝐴𝐿 . (46)

4
Hence, the simplest version of the 𝑁 𝐷 approach we consider here consists of the 𝑁 𝐷 01 equations
∫ 𝐿
1 Δ 𝐿 (𝜔 𝐿 ) 𝐷 𝐿 (𝜔 𝐿 )
𝑁 𝐿 ( 𝐴) = 𝑑𝜔 𝐿 ; (47)
𝜋 −∞ 𝜔𝐿 − 𝐴
𝑚𝑁 𝐴 ∞
∫
𝑁 𝐿 (𝜔 𝑅 )
𝐷 𝐿 ( 𝐴) = 1− 2
𝑑𝜔 𝑅 √ . (48)
4𝜋 0 𝜔 𝑅 (𝜔 𝑅 − 𝐴)

Substituting Eq. (47) into Eq. (48), permuting the RHC and LHC integrals, and employing the algebraic result 2
∫ ∞
𝑑𝜔 𝑅 𝑖𝜋
√ = − √ √ , (49)
0 𝜔 𝑅 (𝜔 𝑅 − 𝐴) (𝜔 𝑅 − 𝜔 𝐿 ) 𝐴 𝜔𝐿 + 𝜔𝐿 𝐴

we get a self-consistent IE for the 𝐷 function that reads

√ ∫ 𝐿
𝑚𝑁 𝐴 𝐷 𝐿 (𝜔 𝐿 ) Δ 𝐿 (𝜔 𝐿 )
𝐷 𝐿 ( 𝐴) = 1 − 𝑖 2
𝑑𝜔 𝐿 √ √ √ . (50)
4𝜋 −∞ ( 𝐴 + 𝜔𝐿 ) 𝜔𝐿

We recall that — provided a local interaction such as OPE — the LHC discontinuity of the amplitude Δ 𝐿 ( 𝐴) (43) can
be found without having to explicitly compute 𝑇𝐿 ( 𝐴). This is done by means of the strategy pioneered in Refs. [5, 6],
where the input is the LHC discontinuity of the on-shell potential

𝑉𝐿 (𝑘 2 ) ≡ 𝑣 𝐿 (𝑘, 𝑘) . (51)

For the particular case of the 3 𝑃1 wave,

1 𝐿 𝐿2
  
4𝜋 𝐿
𝑉1 ( 𝐴) = + + 2 log 1 − (52)
𝑚 𝑁 𝑀𝑁 𝑁 2 𝐴 𝐴 𝐴

[c.f. Eq. (40)], thus (𝐴 ∈ R)

1  2𝜋𝐿  2
(2𝑖) −1 𝑉1 ( 𝐴 + 𝑖0+ ) − 𝑉1 ( 𝐴 − 𝑖0+ )
 
= 𝜃 (𝐿 − 𝐴) . (53)
𝑚 𝑁 𝑀𝑁 𝑁 𝐴

5 Comparison among the three methods and with the data

ind Let 𝛿1Sch (𝑘), 𝛿1LS (𝑘), and 𝛿1𝑁 𝐷 (𝑘) the 3 𝑃1 phase shifts obtained by means of the three non-perturbative methods
that were described above — the Schrödinger equation, the Lippmann–Schwinger equation, and the unsubtracted 𝑁/𝐷
method, respectively. Then, introduce the relative differences

𝛿1Sch (𝑘)
ΔSch
LS (𝑘) = 1− ; (54)
𝛿1LS (𝑘)
𝛿1Sch (𝑘)
ΔSch
𝑁 𝐷 (𝑘) = 1− . (55)
𝛿1𝑁 𝐷 (𝑘)

In Fig. 1 we give a log plot of both functions in the whole elastic regime. It is verified that ΔSch Sch
LS (𝑘) and Δ 𝑁 𝐷 (𝑘) take
−7 −6
typical values of O (10 ) and O (10 ), respectively. This confirms the consistency and equivalence between the three
methods. (Why does the agreement get worse for 𝑘 → 0?)
Besides, we note that Eqs. (24)–(26) give for the 3 𝑃1 -wave first three ERE parameters

𝑎𝑉 = 1.7088 fm3 ; 𝑟 𝑘 = −8.1965 fm−1 ; 𝑣 2 ∼ 1 fm , (56)

where the lack of numerical stability in the extraction of the shape parameter prevents from giving an accurate value.
Alternatively, one can compute these parameters following the strategy proposed in Ref. [7], thus getting

𝑎𝑉 = 1.70889820646 fm3 ; 𝑟 𝑘 = −8.19198434500 fm−1 ; 𝑣 2 = 0.956394207030 fm . (57)

5
 Sch Sch
10-4 ΔLS (k) ΔND (k)

Relative discrepancy
10-5

10-6

10-7

10-8
0 100 200 300 400
k (MeV)
Figure 1: Relative discrepancies (54) and (55) between 3 𝑃1 phase shifts as obtained from the first and second method (solid red
line), and from the first and third method (dotted blue line), as a function of the center-of-mass momentum 𝑘 (in MeV).

It is observed that the numerical accuracy of this method is clearly superior. However, the agreement between both
sets of values can be deemed satisfactory.
We can also compare the 3 𝑃1 phase shifts coming out from any of the three methods with the phenomenology given
by the Granada database [8, 9] for 3 𝑃1 neutron-proton scattering, see Fig. 2. We provide as well the curve resulting
from the truncated ERE (shape-parameter approximation) using the values of Eq. (57) on the one hand, and the values
extracted by the Granada group on the other hand. The latter are [10]

𝑎 𝑉(Gr) = 1.536 fm3 ; 𝑟 𝑘(Gr) = −8.50 fm−1 ; 𝑣 2(Gr) = 0.02 fm . (58)

The somewhat small relative deviations (≈ 10% and ≈ 4%, respectively) between the scattering volume and the
effective momentum as predicted from the 𝑁 𝐷 01 solution (57) and the corresponding phenomenological numbers (58)
attest for the fact that OPE alone actually renders a pretty decent agreement with the partial-wave analysis (PWA) 3 . Yet
it is noteworthy that the PWA-extracted shape parameter is unnaturally small [𝑣 2(𝐺𝑟 ) ≈ (70 𝑚 𝜋 ) −1 ], actually two orders
of magnitude below the 𝑁 𝐷 01 prediction. We remark though that the 𝑁/𝐷 method, unlike the two first methods that
were detailed above, offers us the chance to do better by increasing the number of subtractions in our IEs, thus fixing
in principle an arbitrarily large number of physical inputs in a non-perturbative fashion. [Discuss the fact that fixing
only 𝑎 𝑉 , or even 𝑎 𝑉 and 𝑟 𝑘 , seems not to work here.] In what follows, we explore the scenario where the number of
subtractions 𝔫 and 𝔪 is taken so as to reproduce by construction the physical values of Eq. (58).

6 The N/D23 solution

ind The form of the 𝑁 𝐷 23 equations is found from specializing Eqs. (99) and (100) for 𝔪 = 3. (From now on, we
omit the 𝐿 dependence in our equations since we fix 𝐿 = 1 as we do in the “appendix”.) For convenience, we give
them here explicitly:

𝐴2
∫ 𝐿
4𝜋 𝑎 𝑉 Δ(𝜔 𝐿 ) 𝐷 (𝜔 𝐿 ) 1
𝑁 ( 𝐴) = − 𝐴 + 𝑑𝜔 𝐿 ; (59)
𝑚𝑁 𝜋 −∞ 𝜔2𝐿 𝜔𝐿 − 𝐴
𝑎𝑉 𝑟 𝑘 √
𝐷 ( 𝐴) = 1− 𝐴 − 𝑎 𝑉 𝑣 2 𝐴2 + 𝑖 𝑎 𝑉 𝐴 𝐴
2

2 It is understood that both 𝜔 𝐿 and 𝐴 have positive (arbitrarily small) imaginary parts, preventing any ambiguity when taking the square roots.
√
Nevertheless, Eq. (49) is valid in the whole complex 𝐴 plane provided that 𝐴 is defined with arg( 𝐴) ∈ [0, 2 𝜋 ).
3 This is in contrast with what happens for other low-𝐿 𝑁 𝑁 channels such as the 1 𝑆 , 3 𝑆 , or 3 𝑃 waves, where short-distance physics happens to
0 1 0
play a significantly more relevant role. Of course, these are all waves where OPE is attractive, unlike what happens with the 3 𝑃1 wave.

6
 0

-5

-10
δ1 (deg) -15
ND01
-20
ERE
-25 ERE Granada
Granada
-30
0 100 200 300 400
k (MeV)
Figure 2: 3𝑃
1 phase shift (in degrees) versus center-of-mass momentum 𝑘 (in MeV). The solid (red) line is the curve obtained
from the unsubtracted 𝑁/𝐷 method. The black dots are taken from the Granada partial-wave analysis [8, 9]. The dashed (blue) and
dot-dashed (green) curves stand for the shape-parameter approximation using as inputs the values of Eqs. (57) and (58), respectively.

∫ 𝐿
"  √ #
𝑚𝑁 𝐴 𝐷 (𝜔 𝐿 ) Δ(𝜔 𝐿 ) 1 𝑟𝑘 1 𝐴 1 𝐴 𝐴
+ 𝑑𝜔 𝐿 − + 𝐴 − −𝑖 √ √ √ . (60)
4𝜋 2 −∞ 𝜔2𝐿 𝑎𝑉 2 𝜔 𝐿 𝑎 𝑉 𝜔2𝐿 𝜔 𝐿 ( 𝜔 𝐿 + 𝐴)

Fortunately, Δ( 𝐴) keeps being the same as before, since it only accounts for the large-distance (pion) physics that is
the only source of the LHC discontinuity. By construction, the resulting amplitude (41) fulfills at low energies
√
 
𝐴 𝑚𝑁 1 𝑟𝑘 2 3
lim = − + 𝐴 + 𝑣2 𝐴 − 𝑖 𝐴 𝐴 + O ( 𝐴 ) , (61)
𝐴→0 𝑇 ( 𝐴) 4𝜋 𝑎𝑉 2

where 𝑎 𝑉 , 𝑟 𝑘 , and 𝑣 2 are free parameters.

In Fig. 3 we plot the phase shifts coming out from the new 𝑁 𝐷 23 solution with the Granada parameters (58)
together with the Granada curve. To provide a more complete visual impression, we include as well the phase shifts
due to the original 𝑁 𝐷 01 solution. It is apparent that 𝑁 𝐷 23 improves the agreement with the phenomenology up to
𝑘 ≲ 𝑚 𝜋 . However, at larger momenta there seems to be a lack of repulsion that actually makes the phenomenological
description a bit worse than it was for the 𝑁 𝐷 01 phase shifts. It remains to be seen whether this inconvenience can be
amended by further increasing the number of subtractions. Also, what about TPE?
Note that the 𝑁 𝐷 23 solution that takes the values of the first three ERE parameters from Ref. [10] predicts the
fourth and fifth ERE parameters to be

𝑣3 = 3.3082 fm3 ; 𝑣 4 = −63.6565 fm5 (!) , (62)

according to the algorithm to extract them given in Ref. [7]. In contrast, Ref. [10] provides more natural values

𝑣 3(𝐺𝑟 ) = −1.05 fm3 ; 𝑣 4(𝐺𝑟 ) = 0.56 fm5 . (63)

6.1 𝑁 𝐷 01 as a particular case of 𝑁 𝐷 23

ind A compelling cross-check of the consistency of our calculations is that one must recover the original 𝑁 𝐷 01 solution
from what we call the 𝑁 𝐷 ★23 solution. The latter is found from the 𝑁 𝐷 23 equations simply by fixing the three ERE
parameters that are taken as inputs to the corresponding values predicted by the 𝑁 𝐷 01 solution (57). We quantify in
Fig. 4 the relative deviations between both phase-shift curves for different values of 𝑛max (number of Gauss points
employed along the LHC numerical integration) used to obtain the 𝑁 𝐷 ★23 solution. Consistently, it is found that the
agreement is improved provided that 𝑛max is increased, eventually reaching a satisfactory typical accuracy.

7
 0

-5

-10
δ1 (deg)

-15
ND01
-20
ND23
-25 ERE Granada
Granada
-30
0 100 200 300 400
k (MeV)
Figure 3: 3𝑃
1 phase shift (in degrees) versus center-of-mass momentum 𝑘 (in MeV). The dashed (red) and solid (blue) lines are
the curves corresponding to the 𝑁 𝐷 01 and 𝑁 𝐷 23 solutions, respectively. The black dots are taken from the Granada partial-wave
analysis [8,9]. The dot-dashed (green) curve stands for the shape-parameter approximation using as inputs the Granada values [10].

nmax=1000
1
nmax=2000


0.01 nmax=3000
δND01 (k; 4000)
(k; nmax)

nmax=4000
10-4
23
δND

10-6
*

 1-

10-8

0 100 200 300 400 500

k (MeV)

Figure 4: Relative deviation (in absolute value) of the 𝑁 𝐷★23 phase shift (see the text) with respect to the 𝑁 𝐷 01 phase shift as a
function of the center-of-mass momentum (in MeV) for different numbers of Gauss points 𝑛max = 1000, 2000, 3000, 4000. (The
𝑁 𝐷 01 phase shift has been obtained by fixing 𝑛max = 4000.) [The agreement appears to be worse for 𝑘 → 0, but quickly gets better
as 𝑘 is increased.]

8
7 The N/D34 solution
7.1 From the 𝑁 𝐷 23 to the 𝑁 𝐷 34 IE
√
ind We start from Eqs. (78) and (79) particularized for 𝔫 = 2 and 𝔪 = 3, which can be rewritten as [𝜌( 𝐴) ≡ − 𝑚4 𝜋𝑁 𝐴]

𝐴2 𝐴3
∫ 𝐿 ∫ 𝐿
Δ(𝜔 𝐿 ) 𝐷 (𝜔 𝐿 ) Δ(𝜔 𝐿 ) 𝐷 (𝜔 𝐿 )
𝑁 ( 𝐴) = 𝜈2 𝐴 + 𝑑𝜔 𝐿3
+ 𝑑𝜔 𝐿 3 ; (64)
𝜋 −∞ 𝜔𝐿 𝜋 −∞ 𝜔 𝐿 (𝜔 𝐿 − 𝐴)
∫ ∞
2 𝐴2 ( 𝐴 − 𝐶) 𝜌(𝜔 𝑅 ) 𝑁 (𝜔 𝑅 )
𝐷 ( 𝐴) = 1 + 𝛿2 𝐴 + 𝛿3 𝐴 + 𝑑𝜔 𝑅
𝜋 0 𝜔2𝑅 (𝜔 𝑅 − 𝐶) 2
𝐴2 ( 𝐴 − 𝐶) 2 ∞
∫
𝜌(𝜔 𝑅 ) 𝑁 (𝜔 𝑅 )
+ 𝑑𝜔 𝑅
𝜋 0 𝜔 𝑅 (𝜔 𝑅 − 𝐶) 2 (𝜔 𝑅 − 𝐴)
2
" ∫ 𝐿 #
2 2 1 Δ(𝜔 𝐿 ) 𝐷 (𝜔 𝐿 )
= 1 + 𝛿2 𝐴 + 𝛿3 𝐴 + 𝐴 ( 𝐴 − 𝐶) 𝜈2 𝐺 (0) − 𝑑𝜔 𝐿 𝐺 (𝜔 𝐿 )
𝜋 −∞ 𝜔2𝐿
" ∫ 𝐿 #
2 2 1 Δ(𝜔 𝐿 ) 𝐷 (𝜔 𝐿 )
+ 𝐴 ( 𝐴 − 𝐶) 𝜈2 G(0; 𝐴) − 𝑑𝜔 𝐿 G(𝜔 𝐿 ; 𝐴) , (65)
𝜋 −∞ 𝜔2𝐿

with the auxiliary integrals

∫ ∞
1 𝜌(𝜔 𝑅 )
𝐺 (𝜔 𝐿 ) = 𝑑𝜔 𝑅
𝜋 0(𝜔 𝑅 − 𝐶) 2 (𝜔 𝑅 − 𝜔 𝐿 )
𝑚𝑁 1
= − √ √ √ ; (66)
8𝜋 −𝐶 ( −𝐶 + −𝜔 𝐿 ) 2
∫ ∞
1 𝜌(𝜔 𝑅 )
G(𝜔 𝐿 ; 𝐴) = 𝑑𝜔 𝑅
𝜋 0 (𝜔 𝑅 − 𝐶) 2 (𝜔 𝑅 − 𝐴) (𝜔 𝑅 − 𝜔 𝐿 )
" #
𝑚𝑁 1 1
= − √ √ √ − √ √ . (67)
8𝜋 −𝐶 (𝜔 𝐿 − 𝐴) ( −𝐶 + −𝜔 𝐿 ) 2 ( −𝐶 + −𝐴) 2

It can be checked that Eq. (65) matches Eq. (83) when the latter is specialized for 𝔫 = 2 and 𝔪 = 3. To be continued?
The 𝑁 𝐷 34 equations, which are recovered from taking 𝔪 = 4 in Eqs. (102) and (103), read
4𝜋 𝑎 𝑉  𝑣4 
𝑁 ( 𝐴) = − 𝐴 1− 𝐴
𝑚𝑁 𝑣3
" #
2 ∫ 𝐿  
𝐴 Δ(𝜔 𝐿 ) 𝐷 (𝜔 𝐿 ) 1 1 𝑟𝑘 2 3 √ 1
+ 𝑑𝜔 𝐿 − − + 𝜔 𝐿 + 𝑣 2 𝜔 𝐿 + 𝑣 3 𝜔 𝐿 − 𝑖𝜔 𝐿 𝜔 𝐿 + (;68)
𝜋 −∞ 𝜔2𝐿 𝑣 3 𝜔4𝐿 𝑎𝑉 2 𝜔𝐿 − 𝐴
𝑎𝑉 𝑟 𝑘  2𝑣 4   𝑟 𝑘 𝑣4   𝑣2𝑣4  √  𝑣4 
𝐷 ( 𝐴) = 1− 𝐴 1+ − 𝑎 𝑉 𝑣 2 𝐴2 1 − − 𝑎 𝑉 𝑣 3 𝐴3 1 − 2 + 𝑖 𝑎 𝑉 𝐴 𝐴 1 − 𝐴
2 𝑎𝑉 𝑟 𝑘 𝑣3 2𝑣 2 𝑣 3 𝑣3 𝑣3
∫ 𝐿 ( "  √ #
𝑚𝑁 𝐴 𝐷 (𝜔 𝐿 ) Δ(𝜔 𝐿 ) 𝐴 1 𝑟𝑘 1 𝐴 1 𝐴 𝐴
+ 𝑑𝜔 𝐿 − + 𝐴 − −𝑖 √ √ √
4𝜋 2 −∞ 𝜔2𝐿 𝜔𝐿 𝑎𝑉 2 𝜔 𝐿 𝑎 𝑉 𝜔2𝐿 𝜔 𝐿 ( 𝜔 𝐿 + 𝐴)
)
√
 
1 1 𝑟𝑘 2 1 𝑟𝑘 2 √
− − + 𝐴 + 𝑣2 𝐴 − 𝑖 𝐴 𝐴 − + 𝜔 𝐿 + 𝑣 2 𝜔 𝐿 − 𝑖𝜔 𝐿 𝜔 𝐿 . (69)
𝑣 3 𝜔4𝐿 𝑎𝑉 2 𝑎𝑉 2

Hence, the amplitude around threshold fulfills

√
 
𝐴 𝑚𝑁 1 𝑟𝑘 2 3 4 5
lim = − + 𝐴 + 𝑣2 𝐴 + 𝑣3 𝐴 + 𝑣4 𝐴 − 𝑖 𝐴 𝐴 + O ( 𝐴 ) , (70)
𝐴→0 𝑇 ( 𝐴) 4𝜋 𝑎𝑉 2

where 𝑎 𝑉 , 𝑟 𝑘 , 𝑣 2 , 𝑣 3 , and 𝑣 4 can now be simultaneously fixed to the data, cf. Eqs. (58) and (63).

9
 1 nmax=1000
nmax=2000
0.01 nmax=3000


δND23 (k; 4000)
(k; nmax)
nmax=4000
10-4

10-6
34
δND
*

10-8
 1-

10-10

0 100 200 300 400 500

k (MeV)

Figure 5: Same as in Fig. 4, but replacing 𝛿 𝑁 𝐷01 → 𝛿 𝑁 𝐷23 , 𝛿 𝑁 𝐷★ → 𝛿 𝑁 𝐷★ . Funnily enough, in this case the agreement is better
23 34
at low scattering momenta, unlike Fig. 4.

7.2 𝑁 𝐷 23 as a particular case of 𝑁 𝐷 34

ind Recall that at the end of the previous section (see Fig. 4) we verified that our results are consistent given the very
small relative discrepancy between the 𝑁 𝐷 01 and 𝑁 𝐷 ★23 phase shift along the whole elastic regime. In the same spirit,
now (see Fig. 5) we plot the relative discrepancy between, on the one hand, the 𝑁 𝐷 23 phase shift obtained by fixing
the three first ERE parameters to the Granada values (58) and, on the other hand, the 𝑁 𝐷 ★34 phase shift found by fixing
in Eqs. (68) and (69) the three first ERE parameters to the Granada prescription (58), and the fourth and fifth ERE
parameters to the predictions of the 𝑁 𝐷 23 solution (62). Again, we perform this cross-check for different numbers of
Gauss points 𝑛max , and observe that the agreement tends to improve as 𝑛max is increased. Overall, the correspondence
between both solutions is satisfactory.

7.3 Sensitivity of 𝐷 ( 𝐴) along the LHC on 𝑎 𝑉 and 𝑟 𝑘

ind Here we explore how little changes in the input values of the two first ERE parameters affect the solution of Eq.
(69) for 𝐴 < 𝐿. We take the reference values of the five first ERE parameters

𝑎 𝑉 = 0.525 . . . 𝑚 −3 −1 −3 −5
𝜋 ; 𝑟 𝑘 = −12.154 . . . 𝑚 𝜋 ; 𝑣 2 = 0.014 . . . 𝑚 𝜋 ; 𝑣 3 = 1.132 . . . 𝑚 𝜋 ; 𝑣 4 = −10.648 . . . 𝑚 𝜋 , (71)

obtained from the 𝑁 𝐷 23 solution by means of the procedure explained in Ref. [7]. Then we introduce the solution
along the LHC

𝐷 𝔦,𝔧 ( 𝐴) ≡ 𝐷 ( 𝐴; 𝑎 𝑉 (1 + 𝔦𝜆), 𝑟 𝑘 (1 + 𝔧𝜆)) , (72)

with the small number 𝜆 = 10−4 , thus the relative deviations

𝐷 𝔦,𝔧 ( 𝐴)
Δ𝔦,𝔧 ( 𝐴) ≡ 1− . (73)
𝐷 0,0 ( 𝐴)

In Fig. 6 we represent Δ𝔦,𝔧 as a function of −𝐿 < | 𝐴| < 16𝑚 2𝜋 (with 𝐿 = −𝑚 2𝜋 /4) for 𝔦, 𝔧 = −1, 0, +1. The discrepancies
are surprisingly large, but this does not affect the physical solution, namely the results along the RHC. We confirm this
element in what follows.

7.4 Sensitivity of 𝐷 (𝑘) along the RHC on 𝑎 𝑉 and 𝑟 𝑘

ind Now we show how small shifts in 𝑎 𝑉 and 𝑟 𝑘 render changes in the values taken by the 𝐷 function along the unitarity
cut. We take the 𝐷 function fulfilling (69) corresponding to the 𝑁 𝐷 ★34 case for 𝐴 ⩾ 0, then introduce

𝐷 𝔦,𝔧 (𝑘) ≡ 𝐷 (𝑘; 𝑎 𝑉 (1 + 𝔦𝜆), 𝑟 𝑘 (1 + 𝔧𝜆)) , 𝑘 ⩾ 0, (74)

10
Figure 6: [Generally speaking, it is observed that the deviation is smaller (larger) in the two cases when the variations on 𝑎 𝑉 and
𝑟 𝑘 have opposite (the same) sign, while the remaining four cases where either 𝑎 𝑉 or 𝑟 𝑘 remain unchanged lie somewhere in the
middle.]

with the small number 𝜆 = 10−2 , thus the relative deviations

𝐷 𝔦,𝔧 (𝑘)
Δ𝔦,𝔧 (𝑘) ≡ 1− . (75)
𝐷 0,0 (𝑘)
In Fig. 7 we plot such deviations for 𝔦, 𝔧 = −1, 0, +1 in the interval 0 ⩽ 𝑘 ⩽ 150 MeV, where the effective-range
approximation is expected to hold reasonably. (Discuss the OPE cut at 𝑚 𝜋 /2?) The good news is that, as we
immediately see, such discrepancies are much less drastic than they were along the LHC, c.f. Fig. 6.
Given the unitarity condition ℑ[𝐷 (𝑘)] = 𝜌(𝑘) 𝑁 (𝑘) 𝜃 (𝑘), it is manifest that the 𝑁 function is real along the
unitarity cut, hence the phase shift 𝛿(𝑘) = tan−1 (ℑ[𝑇 (𝑘)]/ℜ[𝑇 (𝑘)]) actually depends only on the 𝐷 function, that
is to say, 𝛿(𝑘) = − tan−1 (ℑ[𝐷 (𝑘)]/ℜ[𝐷 (𝑘)]) — this is of course a general result. In Fig. 8 we plot the relative
differences
𝛿𝔦,𝔧 (𝑘)
Δ𝔦,𝔧 (𝑘)
e ≡ 1− (76)
𝛿0,0 (𝑘)
for 𝔦, 𝔧 = −1, 0, +1 and 0 ⩽ 𝑘 ⩽ 150 MeV. The following elements can be remarked here:

1
Δ-1,0
Δ+1,0
0.100
Relative difference

Δ0,-1
Δ0,+1
0.010 Δ-1,-1
Δ-1,+1
Δ+1,-1
0.001
Δ+1,+1

10-4
0 20 40 60 80 100 120 140
k (MeV)

Figure 7: Δ𝔦,𝔧 (𝑘) function, c.f. (73).

11
 1
Δ-1,0
Δ+1,0

Relative difference
0.100
Δ0,-1
Δ0,+1
0.010 Δ-1,-1
Δ-1,+1
0.001 Δ+1,-1
Δ+1,+1
10-4
0 20 40 60 80 100 120 140
k (MeV)

Figure 8: e
Δ𝔦,𝔧 (𝑘) function, c.f. (76).

ND01
ND*23
8
10
ND23
|N | (m-π2)

104 ND*34

10-4

1 10 100 1000 104 105 106

k (MeV)

Figure 9: Absolute value of the 𝑁 function (in 𝑚 𝜋 units) versus center-of-mass momentum 𝑘 ⩽ 106 MeV. The different curves
correspond to the several 𝑁/𝐷 solutions considered up to now (see the text).

• When 𝑘 → 0, all the six curves coming out from varying 𝑎 𝑉 with respect to the Granada value (𝔦 = ±1) match
at 𝜆. This is consistent with the expectation because the near-threshold phase shift becomes proportional to the
scattering volume, c.f. Eq. (23).
• Actually, for the two cases where the variations upon the scattering volume and the effective momentum have
opposite signs (𝔦 + 𝔧 = 0), one sees that e
Δ𝔦,𝔧 ≈ 𝜆 in the whole momentum regime in consideration — this seems
reasonable as well due to middle-range cancellations. In contrast, the agreement is somewhat worse for the two
curves with |𝔦 + 𝔧| = 2.
• For the two cases where 𝑎 𝑉 remains unchanged but 𝑟 𝑘 is shifted (𝔦 = 0, 𝔧 = ±1) the coincidence is almost perfect
when 𝑘 → 0, but gets worse as 𝑘 is increased. Again, this matches the expectation.

7.5 𝑁 (𝑘) along the RHC

ind In Fig. 9 we plot the dependence of the 𝑁 function respectively associated to the solutions 𝑁 𝐷 01 , 𝑁 𝐷 ★23 , 𝑁 𝐷 23 ,
and 𝑁 𝐷 ★34 that were explored and discussed above.
From inspection of Fig. 9, it seems consistent that the 𝑁 𝐷 23 and 𝑁 𝐷 ★34 curves display great mutual agreement for
𝑘 ≲ 105 MeV, since it matches the fact that both solutions are physically equivalent as we checked above, c.f. Fig. 5.
However, it may look suspicious that the 𝑁 𝐷 01 and 𝑁 𝐷 ★23 curves start to differ significantly for 𝑘 ≳ 100 MeV — much
below the inelastic threshold — as this can be presumed to contradict the equivalence between the two former physical
solutions that was verified in Fig. 4. As another cross-check, I have computed again the 𝑁 function associated to the
𝑁 𝐷 ★23 solution in an indirect way. Generically, the 𝑁 function can be expressed as

4𝜋 𝐷 (𝑘) 4𝜋
𝑁 (𝑘) = = sin 𝛿(𝑘) 𝑒 𝑖 𝛿 (𝑘 ) 𝐷 (𝑘) . (77)
𝑚 𝑁 𝑘 cot 𝛿(𝑘) − 𝑖𝑘 𝑚𝑁 𝑘
It has been verified that both ways to obtain the 𝑁 function differ by an error of numerical origin only, typically of
O (10−13 ) or less.
Besides, in Figs. 10 and 11 we plot
Blah. . .

12
 ND01
1010
ND*23

|Re D|
106
100

1 10 100 1000 104 105 106

k (MeV)

Figure 10: |ℜ 𝐷 (𝑘)| functions for the 𝑁 𝐷 01 and 𝑁 𝐷★23 solutions.

105
|Re T -1| (m-π2)

1000
10 ND01
0.1 ND*23
10 104
k (MeV)

Figure 11: |ℜ 𝑇 −1 (𝑘)| functions for the 𝑁 𝐷 01 and 𝑁 𝐷★23 solutions.

ND01 ND*23
0.100

0.001
θ

10-5

10-7

10 20 30 40 50
k (MeV)

13
 A zero of 𝑁 renders a vanishing phase shift, thus a pole in 𝑇 −1 , often referred as a CDD pole in the literature [11]
(cite to be updated). I obtained that the 𝑁 function crosses the RHC at the following locations 4 :
• 𝑁 𝐷 01 𝑘 0 = 4.2 × 104 MeV
• 𝑁 𝐷 ★23 𝑘 0 = 3.5 × 104 MeV

• 𝑁 𝐷 23 𝑘 0 = 3.5 × 104 , 8.4 × 103 , 535 MeV

• 𝑁 𝐷 ★34 𝑘 0 = 3.5 × 104 , 8.4 × 103 , 535 MeV
Surprisingly to me, from the above it seems that the number of CDD poles does not depend on the number subtractions
of a given 𝑁/𝐷 solution, but rather on the physics coming out from it. Remarkably, when there is a CDD pole, the
pair of consecutive values of the momentum such that 𝑁 changes its sign is always the same. Comment on discrete
versus continuous invariance.

A N/D equations for an arbitrary number of subtractions

ind [J.A.: Presentar notación 𝑁 𝐷 𝔫𝔪 .] Consider an 𝑁/𝐷 decomposition of the 𝑃-wave amplitude matrix encoding
𝔫 ⩾ 2 and 𝔪 ⩾ 2 subtractions, with 𝔪 ⩾ 𝔫, that are applied on the integral equations fulfilled by 𝑁 ( 𝐴) and 𝐷 ( 𝐴),
respectively. All the 𝑁 subtractions are taken at threshold 5 ; only two of the 𝐷 subtractions are applied at 𝐴 = 0, while
the remaining 𝔪−2 subtractions are taken at 𝐴 = 𝐶 < 0. That is to say,
𝔫−1
𝐴𝔫
∫ 𝐿
∑︁
𝑗 Δ(𝜔 𝐿 ) 𝐷 (𝜔 𝐿 )
𝑁 ( 𝐴) = 𝜈 𝑗+1 𝐴 + 𝑑𝜔 𝐿 ; (78)
𝑗=1
𝜋 −∞ 𝜔𝔫𝐿 (𝜔 𝐿 − 𝐴)
𝔪−1 √
𝑚 𝐴2 ( 𝐴 − 𝐶) 𝔪−2 ∞ 𝜔 𝑅 𝑁 (𝜔 𝑅 )
∑︁ ∫
𝑗
𝐷 ( 𝐴) = 1+ 𝛿 𝑗+1 𝐴 − 𝑑𝜔 𝑅 , (79)
𝑗=1
4𝜋 2 0 𝜔2𝑅 (𝜔 𝑅 − 𝐶) 𝔪−2 (𝜔 𝑅 − 𝐴)

where 𝜈2 , . . . , 𝜈𝔫 , 𝛿2 , . . . , 𝛿𝔪 constitute a set of a priori unknown coefficients one can determine from the knowledge
of the first 𝔫+𝔪−2 𝑃-wave effective-range expansion parameters, namely 6 𝑣 0 ≡ −1/𝑎 𝑉 , 𝑣 1 ≡ 𝑟 𝑘 /2, . . . , 𝑣 𝔫+𝔪−3 [J.A.:
cuesta distinguir los dos tipos de nu’s para el desarrollo de N(A) y para el ERE], by means of the matching equation [7]
𝔫+𝔪−3
√

4𝜋 𝐴  ∑︁ 𝑗 𝔫+𝔪−2 

𝐷 ( 𝐴) + i 𝐴 𝐴 𝑁 ( 𝐴) =  𝑣 𝑗 𝐴 + O ( 𝐴 )  𝑁 ( 𝐴) . (80)
𝑚 
 𝑗=0 
 
Substituting Eq. (78) into (79) yields
𝔪−1 𝔫−1
∑︁ 𝑚 𝐴2 ( 𝐴 − 𝐶) 𝔪−2 ∑︁
𝐷 ( 𝐴) = 1+ 𝛿 𝑗+1 𝐴 𝑗 − 𝜈 𝑗+1 𝐺 𝑗 −1,𝔪−2 (0; 𝐴)
𝑗=1
4𝜋 𝑗=1

𝑚 𝐴2 ( 𝐴 − 𝐶) 𝔪−2
∫ 𝐿
Δ(𝜔 𝐿 ) 𝐷 (𝜔 𝐿 )
+ 𝑑𝜔 𝐿 𝐺 𝔫−2,𝔪−2 (𝜔 𝐿 ; 𝐴) , (81)
4𝜋 2
−∞ 𝜔𝔫𝐿
with the auxiliary RHC integral
∞
√
𝜔 𝑎𝑅 𝜔 𝑅
∫
1
𝐺 𝑎,𝑏 (𝜔 𝐿 ; 𝐴) = 𝑑𝜔 𝑅 . (82)
𝜋 0 (𝜔 𝑅 − 𝐶) 𝑏 (𝜔 𝑅 − 𝐴) (𝜔 𝑅 − 𝜔 𝐿 )
Eventually one finds
𝔪−1
∑︁
𝐷 ( 𝐴) = 1 + 𝛿 𝑗+1 𝐴 𝑗
𝑗=1

4 I have found the precise numbers by linear extrapolation between each pair of points where 𝑁 flips sign, but I guess one can refine such estimate.
5 There is no risk of infrared divergences since 𝐿 < 0.
6 Sometimes 𝑎 and 𝑟 are called the scattering volume and the effective momentum, respectively.
𝑉 𝑘

14
 𝔫−1
" 3 3
#
𝑚 𝐴2 ( 𝐴 − 𝐶) 𝔪−2 ∑︁ 𝐴 𝑗− 2 𝑗 − 23 5 𝐶 𝐶 𝑗 −𝔪+ 2
 
1
− 𝜈 𝑗+1 − 𝐹 (1, 𝑗 − ; 𝑗 − 𝔪 + ; )
4𝜋 ( 𝐴 − 𝐶) 𝔪−2 𝔪−3 2 1 2 2 𝐴 𝐴
𝑗=1
" 3 𝔫− 3
i 𝑚 𝐴2 ( 𝐴 − 𝐶) 𝔪−2 𝐴𝔫− 2 𝜔𝐿 2
∫ 𝐿
Δ(𝜔 𝐿 ) 𝐷 (𝜔 𝐿 )
− 𝑑𝜔 𝐿 𝔫 −
4𝜋 2 −∞ 𝜔 𝐿 (𝜔 𝐿 − 𝐴) ( 𝐴 − 𝐶) 𝔪−2 (𝜔 𝐿 − 𝐶) 𝔪−2
3 3
!#
𝔫 − 32 5 𝐶 𝐶 𝔫−𝔪+ 2 5 𝐶 𝐶 𝔫−𝔪+ 2
 
1 1
− 𝐹 (1, 𝔫 − ; 𝔫 − 𝔪 + ; ) − 2 𝐹1 (1, 𝔫 − ; 𝔫 − 𝔪 + ; ) . (83)
𝔪−3 2 1 2 2 𝐴 𝐴 2 2 𝜔𝐿 𝜔𝐿

Through Eqs. (78) and (79), Taylor expanding the left-hand side of Eq. (80) around 𝐴 = 0 one finds

√ 𝔫+𝔪−2
4𝜋 𝐴 ∑︁
𝐷 ( 𝐴) + i 𝐴 𝐴 𝑁 ( 𝐴) = 𝔏 𝑗 𝐴 𝑗 + O ( 𝐴𝔫+𝔪−1 ) , (84)
𝑚 𝑗=1

with
 𝔫−1  3
𝔪 − 3 ∑︁ 𝑘 − 23 𝔪 − 𝑘 − 2 𝜈 𝑘+1
 
4𝜋

𝔏𝑗 = 𝜃 1− 𝑗 + 𝜃 𝑗 −2 𝜃 𝔪− 𝑗 𝛿 𝑗 + i 𝜃 𝑗 −3 (−1) 𝑗+𝔪
𝑚 𝑗 −3 𝔪 − 3 𝑗 − 𝑘 − 3 𝐶 𝑗 −𝑘− 32
𝑘=1 2
𝑗 −3  𝑘 "
𝐶 3− 𝑗
∫ 𝐿 
i ∑︁ 𝔪−2 Δ(𝜔 𝐿 ) 𝐷 (𝜔 𝐿 ) 𝐶
+ 𝜃 𝑗 −3 (−1) 𝑗 𝑑𝜔 𝐿 − −
𝜋 𝑗 − 𝑘 − 3 −∞ 5
𝜔𝐿 (1 − 𝜔𝐶𝐿 ) 𝔪−2
𝑘=0 𝜔 𝐿2
 𝔫− 𝑗+ 3 !#
3 (𝔪 − 𝔫 − 32 ) (𝔪 + 𝑘 − 𝑗 + 1)

𝔪 𝔫− 2 𝐶
2 1 5 𝐶 𝐶
+ (−1) 2 𝐹1 (1, 𝔫 − ; 𝔫 − 𝔪 + ; ) + ,(85)
𝔪−3 𝔫− 3 2 2 𝜔𝐿 𝜔𝐿 ( 𝑗 − 𝑘 − 𝔫 − 23 ) (𝔪 − 2)
𝜔𝐿 2
where 𝜃 𝑖 is the Heaviside function of discrete argument — 0 (1) for 𝑖 = . . . , −2, −1 (𝑖 = 0, 1, . . .). Similarly, through
Eq. (78) it turns out for the right-hand side of Eq. (80)
𝔫+𝔪−3  𝔫+𝔪−2
 ∑︁ 𝔫+𝔪−2 
∑︁
ℜ 𝑗 𝐴 𝑗 + O ( 𝐴𝔫+𝔪−1 ) ,
𝑗


 𝑣 𝑗 𝐴 + O(𝐴 )  𝑁 ( 𝐴) = (86)
 𝑗=0  𝑗=1
 
with
𝑗
" ∫ #
𝐿
∑︁ 𝜃 𝑗 −𝑘−𝔫 Δ(𝜔 𝐿 ) 𝐷 (𝜔 𝐿 )
ℜ𝑗 = 𝜃 𝑗 −𝑘−1 𝜃 𝔫+𝑘− 𝑗 −1 𝜈 𝑗 −𝑘+1 + 𝑑𝜔 𝐿 𝑗 −𝑘+1
𝑣𝑘 . (87)
𝑘=0
𝜋 −∞ 𝜔𝐿
Hence, the set of equalities

𝔏1 = ℜ1 ; (88)
𝔏1 = ℜ2 ; (89)
..
. (90)
𝔏𝔫+𝔪−2 = ℜ𝔫+𝔪−2 (91)

is a system of linear equations from where 𝜈2 , . . . , 𝜈𝔫 , 𝛿2 , . . . , 𝛿𝔪 can be solved in terms of 𝑣 0 , 𝑣 1 , . . . , 𝑣 𝔫+𝔪−3 , as

anticipated above. In turn, such system is encoded by the matrix identity

𝑨𝑿 = 𝑩, (92)

where
𝜈2
© . ª
𝑎 1,1 𝑎 1,2 ... 𝑎 1,𝔫+𝔪−2 ­ . ® 𝑏1
© ª ­ . ® © ª
­ 𝑎 2,1 𝑎 2,2 ... 𝑎 2,𝔫+𝔪−2 ­ ®
­𝜈 ® ­ 𝑏2 ®
𝑿 = ­ 𝔫®;
®
𝑨 = ­­ .. .. .. .. ®; 𝑩 = ­ . ®® ,
­ (93)
­ 𝛿2 ® ­ .. ®
. ®
­ . . . ® ­ . ®
𝑎 𝑎 𝔫+𝔪−2,2 ... 𝑎 𝔫+𝔪−2,𝔫+𝔪−2 ¬ ­ . ®
« 𝔫+𝔪−2,1 ­ . ® «𝑏 𝔫+𝔪−2 ¬
«𝛿 𝔪 ¬

15
 with the matrix elements
" #
3
(−1) 𝑖+𝔪 𝔪 − 3 𝑗 − 23 𝔪 − 𝑗 − 2
  
4𝜋
𝑎 𝑖, 𝑗 = 𝑣 𝑖− 𝑗 𝜃 𝑖− 𝑗 − i 𝜃
3 𝑖−3
𝜃 𝔫−1− 𝑗 − 𝛿𝔫+𝑖−2, 𝑗 𝜃 𝑖−2 ; (94)
𝐶 𝑖− 𝑗 − 3
2 𝑖 − 3 𝔪 − 3 𝑖 − 𝑗 − 𝑚
2
∫ 𝐿 " 𝑖 𝑖−3   𝑘
4𝜋 1 Δ(𝜔 𝐿 )𝐷 (𝜔 𝐿 ) ∑︁ 𝜃 𝑖−𝑘−𝔫 𝑣 𝑘 𝜃 𝑖−3 ∑︁ 𝔪 − 2 𝐶
𝑏𝑖 = 𝜃 1−𝑖 + 𝑑𝜔 𝐿 − + i −
𝑚 𝜋 −∞ 𝜔𝐿 𝑘=0
𝜔𝑖−𝑘
𝐿
(−𝐶) 𝑖 𝑘=0 𝑖 − 𝑘 − 3 𝜔𝐿
3 
!#
3 𝐶 𝔫+ 2 3 3
(𝔪 − 𝔫 − ) (𝔪 + − + 1)
 
𝔫 − 1 5 𝐶 𝐶 2 𝑘 𝑖  𝐶
× (−1) 𝔪 2 𝐹 (1, 𝔫 − ; 𝔫 − 𝔪 + ; ) + − .
𝔪 − 3 𝜔𝔫𝐿 2 1 2 2 𝜔𝐿 𝜔𝐿 (𝑖 − 𝑘 − 𝔫 − 32 ) (𝔪 − 2)
3
(1 − 𝜔𝐿 ) 𝔪−2 𝜔 2
𝐶 𝐿
(95)

In this way, one can express the free parameters contained in the column matrix 𝑿 = 𝑨 −1 𝑩 in terms of physical
magnitudes. Such relationships are then employed in the full DRs (78) and (83), and those enter the expression of the
on-shell amplitude as

𝑁𝔫,𝔪 ( 𝐴)
𝑇𝔫,𝔪 ( 𝐴) = , (96)
𝐷 𝔫,𝔪 ( 𝐴)

where the number of subtractions has been made explicit for better clarity. Now, in virtue of Eq. (80), at low energies
one has
4𝜋 𝐴 (𝔫)
= 𝜙𝔫+𝔪−3 ( 𝐴) + e
𝑣 𝔫+𝔪−2 𝐴𝔫+𝔪−2 + O ( 𝐴𝔫+𝔪−1 ) , (97)
𝑚 𝑇𝔫,𝔪 ( 𝐴)

where we introduced the abbreviation

𝑖
∑︁ √
𝜙𝑖 ( 𝐴) = 𝑣𝑗 𝐴𝑗 − i 𝐴 𝐴 , (98)
𝑗=0

𝑣 𝑖 (𝔫) stands for the first ERE parameter that is not treated as an input, but comes out as a prediction. The 𝑁𝔫,𝔪 ( 𝐴)
and e
and 𝐷 𝔫,𝔪 ( 𝐴) IEs, together with the corresponding lowest-order predictions for the ERE parameters, for the first few
values of 𝔫 and arbitrary 𝔪 = 𝔫, 𝔫+1, . . . are the following:
• For 𝔫 = 2,

𝐴2
∫
4𝜋 𝐴 𝐿
Δ(𝜔 𝐿 ) 𝐷 2,𝔪 (𝜔 𝐿 ) 1
𝑁2,𝔪 ( 𝐴) = + 𝑑𝜔 𝐿 2
; (99)
𝑚 𝑣0 𝜋 −∞ 𝜔𝐿 𝜔𝐿 − 𝐴
𝔪−2 √ √ !
𝜔𝐿 𝜔𝐿 − 𝐴 𝐴
∫ 𝐿
𝜙𝔪−1 ( 𝐴) 𝑚𝐴 Δ(𝜔 𝐿 ) 𝐷 2,𝔪 (𝜔 𝐿 ) ∑︁ 𝐴𝑖
𝐷 2,𝔪 ( 𝐴) = + 𝑑𝜔 𝐿 𝜙𝑖 (𝜔 𝐿 ) 𝑖+1 + i (100)
;
𝑣0 4𝜋 2 −∞ 𝜔2𝐿 𝑖=0
𝜔𝐿 𝜔𝐿 − 𝐴
∫ 𝐿
𝑚 𝑣0 Δ(𝜔 𝐿 ) 𝐷 2,𝑖 (𝜔 𝐿 )
𝑣 𝑖 (2)
e = − 2 𝑑𝜔 𝐿 𝜙𝑖−1 (𝜔 𝐿 ) . (101)
4𝜋 −∞ 𝜔𝑖+2
𝐿

• For 𝔫 = 3,

𝐴2
∫ 𝐿
4𝜋 𝐴  𝑣𝔪  Δ(𝜔 𝐿 ) 𝐷 3,𝔪 (𝜔 𝐿 )  1 𝜙𝔪−1 (𝜔 𝐿 ) 
𝑁3,𝔪 ( 𝐴) = 1− 𝐴 + 𝑑𝜔 𝐿 − ; (102)
𝑚 𝑣0 𝑣 𝔪−1 𝜋 −∞ 𝜔2𝐿 𝜔𝐿 − 𝐴 𝜔𝔪𝐿 𝑣 𝔪−1
1  𝑣𝔪  𝑚𝐴∫ 𝐿 Δ(𝜔 𝐿 ) 𝐷 3,𝔪 (𝜔 𝐿 )
𝐷 3,𝔪 ( 𝐴) = 𝜙𝔪−1 ( 𝐴) − 𝐴 𝜙𝔪−2 ( 𝐴) + 𝑑𝜔 𝐿
𝑣0 𝑣 𝔪−1 4𝜋 2 −∞ 𝜔2𝐿
𝔪−2 √ √ !
∑︁ 𝐴𝑖 𝜙𝔪−2 ( 𝐴) 𝜙𝔪−1 (𝜔 𝐿 ) 𝜔𝐿 𝜔𝐿 − 𝐴 𝐴
× 𝜙𝑖 (𝜔 𝐿 ) 𝑖+1 − +i ; (103)
𝑖=0
𝜔𝐿 𝜔𝔪𝐿 𝑣 𝔪−1 𝜔𝐿 − 𝐴
𝑣 2𝑖−1
∫ 𝐿
𝑚 𝑣0 Δ(𝜔 𝐿 ) 𝐷 3,𝑖−1 (𝜔 𝐿 ) h 𝑣 𝑖−1  i
𝑣 𝑖 (3) = − 𝑑𝜔 𝐿 𝜙 𝑖−2 (𝜔 𝐿 ) − 𝜔 𝐿 𝜙 𝑖−3 (𝜔 𝐿 ) . (104)
4𝜋 2 −∞
e
𝑣 𝑖−2 𝜔𝑖+2
𝐿
𝑣 𝑖−2

16
 • For 𝔫 = 4,
!
(𝑣 𝔪−2 𝑣 𝔪+1 − 𝑣 𝔪−1 𝑣 𝔪 ) 𝐴 + 𝑣 2𝔪 − 𝑣 𝔪−1 𝑣 𝔪+1 𝐴2


𝐴2
∫ 𝐿
4𝜋 𝐴 Δ(𝜔 𝐿 )𝐷 4,𝔪 (𝜔 𝐿 )
𝑁4,𝔪 ( 𝐴) = 1+ 2
+ 𝑑𝜔 𝐿
𝑚 𝑣0 𝑣 𝔪−1 − 𝑣 𝔪−2 𝑣 𝔪 𝜋 −∞ 𝜔2𝐿
" #
1 (𝑣 𝔪−2 − 𝑣 𝔪−1 𝐴) 𝜙𝔪 (𝜔 𝐿 ) − 𝜔 𝐿 (𝑣 𝔪−1 − 𝑣 𝔪 𝐴) 𝜙𝔪−1 (𝜔 𝐿 )
× +   ; (105)
𝜔𝐿 − 𝐴 𝜔𝔪+1 𝑣 2 −𝑣 𝑣
𝐿 𝔪−1 𝔪−2 𝔪
!
(𝑣 𝔪−2 𝑣 𝔪+1 − 𝑣 𝔪−1 𝑣 𝔪 ) 𝐴 𝜙𝔪−2 ( 𝐴) + 𝑣 2𝔪 − 𝑣 𝔪−1 𝑣 𝔪+1 𝐴2 𝜙𝔪−3 ( 𝐴)


1
𝐷 4,𝔪 ( 𝐴) = 𝜙𝔪−1 ( 𝐴) +
𝑣0 𝑣 2𝔪−1 − 𝑣 𝔪−2 𝑣 𝔪
" 𝔪−2
!
𝑚𝐴
∫ 𝐿
Δ(𝜔 𝐿 ) 𝐷 4,𝔪 (𝜔 𝐿 ) ∑︁ 𝐴𝑖 1 𝜙𝔪−2 ( 𝐴) 𝜙𝔪−3 ( 𝐴) 𝜔𝐴𝐿
+ 2 𝑑𝜔 𝐿 𝜙𝑖 (𝜔 𝐿 ) 𝑖+1 − +
4𝜋 −∞ 𝜔2𝐿 𝑖=0
𝜔𝐿 𝜔𝐿 𝑣2
1 − 𝔪−1 1 − 𝑣𝔪−2
𝑣2
𝑣𝔪
𝑣𝔪−2 𝑣𝔪 𝔪−1

(𝑣 𝔪−2 − 𝑣 𝔪−1 𝜔 𝐿 ) 𝜙𝔪−1 (𝜔 𝐿 ) 𝜙𝔪−2 ( 𝐴) − 𝐴 (𝑣 𝔪−1 − 𝑣 𝔪 𝜔 𝐿 ) 𝜙𝔪−2 (𝜔 𝐿 ) 𝜙𝔪−3 ( 𝐴)

+  
𝜔𝔪+1
𝐿 𝑣 2
𝔪−1 − 𝑣 𝔪−2 𝑣 𝔪
√ √ #
𝜔𝐿 𝜔𝐿 − 𝐴 𝐴
+i ; (106)
𝜔𝐿 − 𝐴
𝑣 2𝑖−2 − 𝑣 𝑖−3 𝑣 𝑖−1
∫ 𝐿
𝑣 𝑖−3 𝑣 𝑖−2 − 𝑣 𝑖−4 𝑣 𝑖−1 𝑚 𝑣0 Δ(𝜔 𝐿 ) 𝐷 4,𝑖−2 (𝜔 𝐿 )
𝑣 𝑖 (4)
e = 2
𝑣 𝑖−1 − 2
𝑣 𝑖−2 − 2
𝑑𝜔 𝐿
𝑣 𝑖−3 − 𝑣 𝑖−4 𝑣 𝑖−2 𝑣 𝑖−3 − 𝑣 𝑖−4 𝑣 𝑖−2 4𝜋 −∞ 𝜔𝑖+2
𝐿
" #
2 −
𝑣 𝑖−3 𝑣 𝑖−2 − 𝑣 𝑖−4 𝑣 𝑖−1   𝑣 𝑖−2 𝑣 𝑣
𝑖−3 𝑖−1

2

× 𝜙𝑖−3 (𝜔 𝐿 ) − 𝜔 𝐿 𝜙𝑖−4 (𝜔 𝐿 ) + 2 𝜔 𝐿 𝜙𝑖−5 (𝜔 𝐿 ) (107)
.
𝑣 2𝑖−3 − 𝑣 𝑖−4 𝑣 𝑖−2 𝑣 𝑖−3 − 𝑣 𝑖−4 𝑣 𝑖−2

J.A.: Nosotros obtuvimos la 23 como nueva solución. Si pensamos en la ambigüedad de CDD’s, se debería ir
aumentado también el índice n, puesto que se están introduciendo ceros, y no solo aumentar el índice m. Así de 23
pasaría a 34.

References
[1] D.B. Kaplan, M.J. Savage, and M.B. Wise. Phys. Lett. B 424 (1998) 390.
[2] D.B. Kaplan, M.J. Savage, and M.B. Wise. Nucl. Phys. B 534 (1998) 329.
[3] A. Nogga, R.G.E. Timmermans, and U. van Kolck. Phys. Rev. C 72 (2005) 054006.

[4] M. Pavón Valderrama and E. Ruiz Arriola, Phys. Rev. C 74 (2006) 064004; 75 (2007) 059905(E).
[5] D.R. Entem and J.A. Oller, Phys. Lett. B 773 (2017) 498.
[6] J.A. Oller and D.R. Entem, Ann. Physics 411 (2019) 167965.
[7] J.A. Oller, Phys. Rev. C 93 (2016) 024002.

[8] R. Navarro-Pérez, J.E. Amaro, and E. Ruiz Arriola, Phys. Rev. C 88 (2013) 024002.
[9] R. Navarro-Pérez, J.E. Amaro, and E. Ruiz Arriola, Phys. Rev. C 88 (2013) 064002.
[10] R. Navarro-Pérez, J.E. Amaro, and E. Ruiz Arriola, J. Phys. G 43 (2016) 114001.

[11] F.J. Dyson, Phys. Rev. 106 (1957) 157.

17