Transverse-momentum-dependent parton distributions up to N$^3$LL from Drell-Yan data

We present an extraction of unpolarised Transverse-Momentum-Dependent Parton Distribution Functions based on Drell-Yan production data from different experiments, including those at the LHC, and spanning a wide kinematic range. We deal with experimental uncertainties by properly taking into account correlations. We include resummation of logarithms of the transverse momentum of the vector boson up to N$^3$LL order, and we include non-perturbative contributions. These ingredients allow us to obtain a remarkable agreement with the data.


Introduction
The analysis of hard scattering processes involving nucleons in the initial state allows us to obtain information on their internal structure, encoded in parton distribution functions (PDFs).
After decades of studies, we have obtained a detailed knowledge of unpolarised collinear PDFs: they provide information about matter at the subnuclear level and are indispensable in almost any prediction involving high-energy hadrons. Collinear PDFs describe the distribution of partons inside the nucleon as a function of the longitudinal momentum fraction x. Collinear factorisation theorems lead to a precise definition of collinear PDFs based on 2 Theoretical framework In this section we describe the theoretical framework of our analysis. In section 2.1, we review the TMD factorisation formula for the Drell-Yan (DY) process. In section 2.2, we briefly describe the evolution of TMDs and how they can be matched onto the collinear JHEP07(2020)117 P B nucleon < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > P A nucleon < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > k A k B k ⊥A k ⊥B quark < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > quark < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > photon < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > q T q P 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " S F m T F Y 2 v I q 9 F O + D + c z s C m w t A Z b s = " > A A A B 6 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E q M e i F 4 8 V 7 A e 0 o W y 2 m 3 b p 7 i b s T o Q S + h e 8 e F D E q 3 / I m / / G p M 1 B W x 8 M P N 6 b Y W Z e E E t h 0 X W / n d L G 5 t b 2 T n m 3 s r d / c H h U P T 7 p 2 C g x j L d Z J C P T C 6 j l U m j e R o G S 9 2 L D q Q o k 7 w b T u 9 z v P n F j R a Q f c R Z z X 9 G x F q F g F H O p N f Q q w 2 r N r b s L k H X i F a Q G B V r D 6 t d g F L F E c Y 1 M U m v 7 n h u j n 1 K D g k k + r w w S y 2 P K p n T M + x n V V H H r p 4 t b 5 + Q i U 0 Y k j E x W G s l C / T 2 R U m X t T A V Z p 6 I 4 s a t e L v 7 n 9 R M M b / x U 6 D h B r t l y U Z h I g h H J H y c j Y T h D O c s I Z U Z k t x I 2 o Y Y y z O L J Q / B W X 1 4 n n a u 6 5 9 a 9 h + t a 8 7 a I o w x n c A 6 X 4 E E D m n A P L W g D g w k 8 w y u 8 O c p 5 c d 6 d j 2 V r y S l m T u E P n M 8 f A 7 C N j A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " S F m T F Y 2 v I q 9 F O + D + c z s C m w t A Z b s = " > A A A B 6 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E q M e i F 4 8 V 7 A e 0 o W y 2 m 3 b p 7 i b s T o Q S + h e 8 e F D E q 3 / I m / / G p M 1 B W x 8 M P N 6 b Y W Z e E E t h 0 X W / n d L G 5 t b 2 T n m 3 s r d / c H h U P T 7 p 2 C g x j L d Z J C P T C 6 j l U m j e R o G S 9 2 L D q Q o k 7 w b T u 9 z v P n F j R a Q f c R Z z X 9 G x F q F g F H O p N f Q q w 2 r N r b s L k H X i F a Q G B V r D 6 t d g F L F E c Y 1 M U m v 7 n h u j n 1 K D g k k + r w w S y 2 P K p n T M + x n V V H H r p 4 t b 5 + Q i U 0 Y k j E x W G s l C / T 2 R U m X t T A V Z p 6 I 4 s a t e L v 7 n 9 R M M b / x U 6 D h B r t l y U Z h I g h H J H y c j Y T h D O c s I Z U Z k t x I 2 o Y Y y z O L J Q / B W X 1 4 n n a u 6 5 9 a 9 h + t a 8 7 a I o w x n c A 6 X 4 E E D m n A P L W g D g w k 8 w y u 8 O c p 5 c d 6 d j 2 V r y S l m T u E P n M 8 f A 7 C N j A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " S F m T F Y 2 v I q 9 F O + D + c z s C m w t A Z b s = " > A A A B 6 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E q M e i F 4 8 V 7 A e 0 o W y 2 m 3 b p 7 i b s T o Q S + h e 8 e F D E q 3 / I m / / G p M 1 B W x 8 M P N 6 b Y W Z e E E t h 0 X W / n d L G 5 t b 2 T n m 3 s r d / c H h U P T 7 p 2 C g x j L d Z J C P T C 6 j l U m j e R o G S 9 2 L D q Q o k 7 w b T u 9 z v P n F j R a Q f c R Z z X 9 G x F q F g F H O p N f Q q w 2 r N r b s L k H X i F a Q G B V r D 6 t d g F L F E c Y 1 M U m v 7 n h u j n 1 K D g k k + r w w S y 2 P K p n T M + x n V V H H r p 4 t b 5 + Q i U 0 Y k j E x W G s l C / T 2 R U m X t T A V Z p 6 I 4 s a t e L v 7 n 9 R M M b / x U 6 D h B r t l y U Z h I g h H J H y c j Y T h D O c s I Z U Z k t x I 2 o Y Y y z O L J Q / B W X 1 4 n n a u 6 5 9 a 9 h + t a 8 7 a I o w x n c A 6 X 4 E E D m n A P L W g D g w k 8 w y u 8 O c p 5 c d 6 d j 2 V r y S l m T u E P n M 8 f A 7 C N j A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " S F m T F Y 2 v I q 9 F O + D + c z s C m w t A Z b s = " > A A A B 6 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E q M e i F 4 8 V 7 A e 0 o W y 2 m 3 b p 7 i b s T o Q S + h e 8 e F D E q 3 / I m / / G p M 1 B W x 8 M P N 6 b Y W Z e E E t h 0 X W / n d L G 5 t b 2 T n m 3 s r d / c H h U P T 7 p 2 C g x j L d Z J C P T C 6 j l U m j e R o G S 9 2 L D q Q o k 7 w b T u 9 z v P n F j R a Q f c R Z z X 9 G x F q F g F H O p N f Q q w 2 r N r b s L k H X i F a Q G B V r D 6 t d g F L F E c Y 1 M U m v 7 n h u j n 1 K D g k k + r w w S y 2 P K p n T M + x n V V H H r p 4 t b 5 + Q i U 0 Y k j E x W G s l C / T 2 R U m X t T A V Z p 6 I 4 s a t e L v 7 n 9 R M M b / x U 6 D h B r t l y U Z h I g h H J H y c j Y T h D O c s I Z U Z k t x I 2 o Y Y y z O L J Q / B W X 1 4 n n a u 6 5 9 a 9 h + t a 8 7 a I o w x n c A 6 X 4 E E D m n A P L W g D g w k 8 w y u 8 O c p 5 c d 6 d j 2 V r y S l m T u E P n M 8 f A 7 C N j A = = < / l a t e x i t > P 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 V l 7 Y Z 4 5 J 7 9 U a M D E S S X s r g t V 2 z c = " > A A A B 6 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K o M e i F 4 8 V 7 A e 0 o W y 2 m 3 b p 7 i b s T o Q S + h e 8 e F D E q 3 / I m / / G p M 1 B W x 8 M P N 6 b Y W Z e E E t h 0 X W / n d L G 5 t b 2 T n m 3 s r d / c H h U P T 7 p 2 C g x j L d Z J C P T C 6 j l U m j e R o G S 9 2 L D q Q o k 7 w b T u 9 z v P n F j R a Q f c R Z z X 9 G x F q F g F H O p N W x U h t W a W 3 c X I O v E K 0 g N C r S G 1 a / B K G K J 4 h q Z p N b 2 P T d G P 6 U G B Z N 8 X h k k l s e U T e m Y 9 z O q q e L W T x e 3 z s l F p o x I G J m s N J K F + n s i p c r a m Q q y T k V x Y l e 9 X P z P 6 y c Y 3 v i p 0 H G C X L P l o j C R B C O S P 0 5 G w n C G c p Y R y o z I b i V s Q g 1 l m M W T h + C t v r x O O o 2 6 5 9 a 9 h 6 t a 8

a I o w x n c A 6 X 4 M E 1 N O E e W t A G B h N 4 h l d 4 c 5 T z 4 r w H 8 v W k l P M n M I f O J 8 / B
T W N j Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 V l 7 Y Z 4 5 J 7 9 U a M D E S S X s r g t V 2 z c = " > A A A B 6 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K o M e i F 4 8 V 7 A e 0 o W y 2 m 3 b p 7 i b s T o Q S + h e 8 e F D E q 3 / I m / / G p M 1 B W x 8 M P N 6 b Y W Z e E E t h 0 X W / n d L G 5 t b 2 T n m 3 s r d / c H h U P T 7 p 2 C g x j L d Z J C P T C 6 j l U m j e R o G S 9 2 L D q Q o k 7 w b T u 9 z v P n F j R a Q f c R Z z X 9 G x F q F g F H O p N W x U h t W a W 3 c X I O v E K 0 g N C r S G 1 a / B K G K J 4 h q Z p N b 2 P T d G P 6 U G B Z N 8 X h k k l s e U T e m Y 9 z O q q e L W T x e 3 z s l F p o x I G J m s N J K F + n s i p c r a m Q q y T k V x Y l e 9 X P z P 6 y c Y 3 v i p 0 H G C X L P l o j C R B C O S P 0 5 G w n C G c p Y R y o z I b i V s Q g 1 l m M W T h + C t v r x O O o 2 6 5 9 a 9 h 6 t a 8

a I o w x n c A 6 X 4 M E 1 N O E e W t A G B h N 4 h l d 4 c 5 T z 4 r w H 8 v W k l P M n M I f O J 8 / B
T W N j Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 V l 7 Y Z 4 5 J 7 9 U a M D E S S X s r g t V 2 z c = " > A A A B 6 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K o M e i F 4 8 V 7 A e 0 o W y 2 m 3 b p 7 i b s T o Q S + h e 8 e F D E q 3 / I m / / G p M 1 B W x 8 M P N 6 b Y W Z e E E t h 0 X W / n d L G 5 t b 2 T n m 3 s r d / c H h U P T 7 p 2 C g x j L d Z J C P T C 6 j l U m j e R o G S 9 2 L D q Q o k 7 w b T u 9 z v P n F j R a Q f c R Z z X 9 G x F q F g F H O p N W x U h t W a W 3 c X I O v E K 0 g N C r S G 1 a / B K G K J 4 h q Z p N b 2 P T d G P 6 U G B Z N 8 X h k k l s e U T e m Y 9 z O q q e L W T x e 3 z s l F p o x I G J m s N J K F + n s i p c r a m Q q y T k V x Y l e 9 X P z P 6 y c Y 3 v i p 0 H G C X L P l o j C R B C O S P 0 5 G w n C G c p Y R y o z I b i V s Q g 1 l m M W T h + C t v r x O O o 2 6 5 9 a 9 h 6 t a 8 7 a I o w x n c A 6 X 4 M E 1 N O E e W t A G B h N 4 h l d 4 c 5 T z 4 r w 7 H 8 v W k l P M n M I f O J 8 / B T W N j Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 V l 7 Y Z 4 5 J 7 9 U a M D E S S X s r g t V 2 z c = " > A A A B 6 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K o M e i F 4 8 V 7 A e 0 o W y 2 m 3 b p 7 i b s T o Q S + h e 8 e F D E q 3 / I m / / G p M 1 B W x 8 M P N 6 b Y W Z e E E t h 0 X W / n d L G 5 t b 2 T n m 3 s r d / c H h U P T 7 p 2 C g x j L d Z J C P T C 6 j l U m j e R o G S 9 2 L D q Q o k 7 w b T u 9 z v P n F j R a Q f c R Z z X 9 G x F q F g F H O p N W x U h t W a W 3 c X I O v E K 0 g N C r S G 1 a / B K G K J 4 h q Z p N b 2 P T d G P 6 U G B Z N 8 X h k k l s e U T e m Y 9 z O q q e L W T x e 3 z s l F p o x I G J m s N J K F + n s i p c r a m Q q y T k V x Y l e 9 X P z P 6 y c Y 3 v i p 0 H G C X L P l o j C R B C O S P 0 5 G w n C G c p Y R y o z I b i V s Q g 1 l m M W T h + C t v r x O O o 2 6 5 9 a 9 h 6 t a 8 7 a I o w x n c A 6 X 4 M E 1 N O E e W t A G B h N 4 h l d 4 c 5 T z 4 r w 7 H 8 v W k l P M n M I f O J 8 / B T W N j Q = = < / l a t e x i t > k 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " i R + q A F 2 H j 9 H 8 E / B H B Q z w z 7 O R s + k = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K U I 9 F L x 4 r 2 g 9 o Q 9 l s J + 3 S z S b s b o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p Y T K o D c o V t + o u Q N a J l 5 M K 5 G g O y l / 9 Y c z S C K V h g m r d 8 9 z E + B l V h j O B s 1 I / 1 Z h Q N q E j 7 F k q a Y T a z x a n z s i F V Y Y k j J U t a c h C / T 2 R 0 U j r a R T Y z o i a s V 7 1 5 u J / X i 8 1 4 b W f c Z m k B i V b L g p T Q U x M 5 n + T I V f I j J h a Q p n i 9 l b C x l R R Z m w 6 J R u C t / r y O m n X q p 5 b 9 e 6 v K o 2 b P I 4 i n M E 5 X I I H d W j A H T S h B Q x G 8 A y v 8 O Y I 5 8 V 5 d z 6 W r Q U n n z m F P 3 A + f w D 5 V 4 2 U < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " i R + q A F 2 H j 9 H 8 E / B H B Q z w z 7 O R s + k = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K U I 9 F L x 4 r 2 g 9 o Q 9 l s J + 3 S z S b s b o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p Y T K o D c o V t + o u Q N a J l 5 M K 5 G g O y l / 9 Y c z S C K V h g m r d 8 9 z E + B l V h j O B s 1 I / 1 Z h Q N q E j 7 F k q a Y T a z x a n z s i F V Y Y k j J U t a c h C / T 2 R 0 U j r a R T Y z o i a s V 7 1 5 u J / X i 8 1 4 b W f c Z m k B i V b L g p T Q U x M 5 n + T I V f I j J h a Q p n i 9 l b C x l R R Z m w 6 J R u C t / r y O m n X q p 5 b 9 e 6 v K o 2 b P I 4 i n M E 5 X I I H d W j A H T S h B Q x G 8 A y v 8 O Y I 5 8 V 5 d z 6 W r Q U n n z m F P 3 A + f w D 5 V 4 2 U < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " i R + q A F 2 H j 9 H 8 E / B H B Q z w z 7 O R s + k = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K U I 9 F L x 4 r 2 g 9 o Q 9 l s J + 3 S z S b s b o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p Y T K o D c o V t + o u Q N a J l 5 M K 5 G g O y l / 9 Y c z S C K V h g m r d 8 9 z E + B l V h j O B s 1 I / 1 Z h Q N q E j 7 F k q a Y T a z x a n z s i F V Y Y k j J U t a c h C / T 2 R 0 U j r a R T Y z o i a s V 7 1 5 u J / X i 8 1 4 b W f c Z m k B i V b L g p T Q U x M 5 n + T I V f I j J h a Q p n i 9 l b C x l R R Z m w 6 J R u C t / r y O m n X q p 5 b 9 e 6 v K o 2 b P I 4 i n M E 5 X I I H d W j A H T S h B Q x G 8 A y v 8 O Y I 5 8 V 5 d z 6 W r Q U n n z m F P 3 A + f w D 5 V 4 2 U < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " i R + q A F 2 H j 9 H 8 E / B H B Q z w z 7 O R s + k = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K U I 9 F L x 4 r 2 g 9 o Q 9 l s J + 3 S z S b s b o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p Y T K o D c o V t + o u Q N a J l 5 M K 5 G g O y l / 9 Y c z S C K V h g m r d 8 9 z E + B l V h j O B s 1 I / 1 Z h Q N q E j 7 F k q a Y T a z x a n z s i F V Y Y k j J U t a c h C / T 2 R 0 U j r a R T Y z o i a s V 7 1 5 u J / X i 8 1 4 b W f c Z m k B i V b L g p T Q U x M 5 n + T I V f I j J h a Q p n i 9 l b C x l R R Z m w 6 J R u C t / r y O m n X q p 5 b 9 e 6 v K o 2 b P I 4 i n M E 5 X I I H d W j A H T S h B Q x G 8 A y v 8 O Y I 5 8 V 5 d z 6 W r Q U n n z m F P 3 A + f w D 5 V 4 2 U < / l a t e x i t > k 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 m 0 X A h F b 0 T M d l 4 v U E e G z Q t r y S L w = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / o h 6 9 L B b B U 0 l E 0 G P R i 8 e K 9 g P a U D b b T b t 0 s w m 7 E 6 G E / g Q v H h T x 6 i / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v T K U w 6 H n f T m l t f W N z q 7 x d 2 d n d 2 z 9 w D 4 9 a J s k 0 4 0 2 W y E R 3 Q m q 4 F I o 3 U a D k n V R z G o e S t 8 P x 7 c x v P 3 F t R K I e c Z L y I K Z D J S L B K F r p Y d z 3 + 2 7 V q 3 l z k F X i F 6 Q K B R p 9 9 6 s 3 S F g W c 4 V M U m O 6 v p d i k F O N g k k + r f Q y w 1 P K x n T I u 5 Y q G n M T 5 P N T p + T M K g M S J d q W Q j J X f 0 / k N D Z m E o e 2 M 6 Y 4 M s v e T P z P 6 2 Y Y X Q e 5 U G m G X L H F o i i T B B M y + 5 s M h O Y M 5 c Q S y r S w t x I 2 o p o y t O l U b A j + 8 s u r p H V R 8 7 2 a f 3 9 Z r d 8 U c Z T h B E 7 h H H y 4 g j r c Q Q O a w G A I z / A K b 4 5 0 X p x 3 5 2 P R W n K K m W P 4 A + f z B / f T j Z M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 m 0 X A h F b 0 T M d l 4 v U E e G z Q t r y S L w = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / o h 6 9 L B b B U 0 l E 0 G P R i 8 e K 9 g P a U D b b T b t 0 s w m 7 E 6 G E / g Q v H h T x 6 i / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v T K U w 6 H n f T m l t f W N z q 7 x d 2 d n d 2 z 9 w D 4 9 a J s k 0 4 0 2 W y E R 3 Q m q 4 F I o 3 U a D k n V R z G o e S t 8 P x 7 c x v P 3 F t R K I e c Z L y I K Z D J S L B K F r p Y d z 3 + 2 7 V q 3 l z k F X i F 6 Q K B R p 9 9 6 s 3 S F g W c 4 V M U m O 6 v p d i k F O N g k k + r f Q y w 1 P K x n T I u 5 Y q G n M T 5 P N T p + T M K g M S J d q W Q j J X f 0 / k N D Z m E o e 2 M 6 Y 4 M s v e T P z P 6 2 Y Y X Q e 5 U G m G X L H F o i i T B B M y + 5 s M h O Y M 5 c Q S y r S w t x I 2 o p o y t O l U b A j + 8 s u r p H V R 8 7 2 a f 3 9 Z r d 8 U c Z T h B E 7 h H H y 4 g j r c Q Q O a w G A I z / A K b 4 5 0 X p x 3 5 2 P R W n K K m W P 4 A + f z B / f T j Z M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 m 0 X A h F b 0 T M d l 4 v U E e G z Q t r y S L w = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / o h 6 9 L B b B U 0 l E 0 G P R i 8 e K 9 g P a U D b b T b t 0 s w m 7 E 6 G E / g Q v H h T x 6 i / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v T K U w 6 H n f T m l t f W N z q 7 x d 2 d n d 2 z 9 w D 4 9 a J s k 0 4 0 2 W y E R 3 Q m q 4 F I o 3 U a D k n V R z G o e S t 8 P x 7 c x v P 3 F t R K I e c Z L y I K Z D J S L B K F r p Y d z 3 + 2 7 V q 3 l z k F X i F 6 Q K B R p 9 9 6 s 3 S F g W c 4 V M U m O 6 v p d i k F O N g k k + r f Q y w 1 P K x n T I u 5 Y q G n M T 5 P N T p + T M K g M S J d q W Q j J X f 0 / k N D Z m E o e 2 M 6 Y 4 M s v e T P z P 6 2 Y Y X Q e 5 U G m G X L H F o i i T B B M y + 5 s M h O Y M 5 c Q S y r S w t x I 2 o p o y t O l U b A j + 8 s u r p H V R 8 7 2 a f 3 9 Z r d 8 U c Z T h B E 7 h H H y 4 g j r c Q Q O a w G A I z / A K b 4 5 0 X p x 3 5 2 P R W n K K m W P 4 A + f z B / f T j Z M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 m 0 X A h F b 0 T M d l 4 v U E e G z Q t r y S L w = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / o h 6 9 L B b B U 0 l E 0 G P R i 8 e K 9 g P a U D b b T b t 0 s w m 7 E 6 G E / g Q v H h T x 6 i / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v T K U w 6 H n f T m l t f W N z q 7 x d 2 d n d 2 z 9 w D 4 9 a J s k 0 4 0 2 W y E R 3 Q m q 4 F I o 3 U a D k n V R z G o e S t 8 P x 7 c x v P 3 F t R K I e c Z L y I K Z D J S L B K F r p Y d z 3 + 2 7 V q 3 l z k F X i F 6 Q K B R p 9 9 6 s 3 S F g W c 4 V M U m O 6 v p d i k F O N g k k + r f Q y w 1 P K x n T I u 5 Y q G n M T 5 P N T p + T M K g M S J d q W Q j J X f 0 / k N D Z m E o e 2 M 6 Y 4 M s v e T P z P 6 2 Y Y X Q e 5 U G m G X L H F o i i T B B M y + 5 s M h O Y M 5 c Q S y r S w t x I 2 o p o y t O l U b A j + 8 s u r p H V R 8 7 2 a f 3 9 Z r d 8 U c Z T h B E 7 h H H y 4 g j r c Q Q O a w G A I z / A K b 4 5 0 X p x 3 5 2 P R W n K K m W P 4 A + f z B / f T j Z M = < / l a t e x i t > k ?2 < l a t e x i t s h a 1 _ b a s e 6 4 = " w D A g F r s 8 w w b H 7 f P S Y Z F T v y U + S Y g = " > A A A B 8 n i c b V B N S w M x E M 3 W r 1 q / q h 6 9 B I v g q e w W Q Y 9 F L x 4 r 2 A / Y L i W b Z t v Q b B K S W a E s / R l e P C j i 1 V / j z X 9 j 2 u 5 B W x 8 M P N 6 b Y W Z e r A W 3 4 P v f X m l j c 2 t 7 p 7 x b 2 d s / O D y q H p 9 0 r M o M Z W 2 q h D K 9 m F g m u G R t 4 C B Y T x t G 0 l i w b j y 5 m / v d J 2 Y s V / I R p p p F K R l J n n B K w E n h Z J D 3 N T M a N 2 a D a s 2 v + w v g d R I U p I Y K t A b V r / 5 Q 0 S x l E q g g 1 o a B r y H K i Q F O B Z t V + p l l m t A J G b H Q U U l S Z q N 8 c f I M X z h l i B N l X E n A C / X 3 R E 5 S a 6 d p 7 D p T A m O 7 6 s 3 F / 7 w w g + Q m y r n U G T B J l 4 u S T G B Q e P 4 / H n L D K I i p I 4 Q a 7 m 7 F d E w M o e B S q r g Q g t W X 1 0 m n U Q / 8 e v B w V W v e F n G U 0 R k 6 R 5 c o Q N e o i e 5 R C 7 U R R Q o 9 o 1 f 0 5 o H 3 4 r 1 7 H 8 v W k l f M n K I / 8 D 5 / A A M z k Q 8 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " w D A g F r s 8 w w b H 7 f P S Y Z F T v y U + S Y g = " > A A A B 8 n i c b V B N S w M x E M 3 W r 1 q / q h 6 9 B I v g q e w W Q Y 9 F L x 4 r 2 A / Y L i W b Z t v Q b B K S W a E s / R l e P C j i 1 V / j z X 9 j 2 u 5 B W x 8 M P N 6 b Y W Z e r A W 3 4 P v f X m l j c 2 t 7 p 7 x b 2 d s / O D y q H p 9 0 r M o M Z W 2 q h D K 9 m F g m u G R t 4 C B Y T x t G 0 l i w b j y 5 m / v d J 2 Y s V / I R p p p F K R l J n n B K w E n h Z J D 3 N T M a N 2 a D a s 2 v + w v g d R I U p I Y K t A b V r / 5 Q 0 S x l E q g g 1 o a B r y H K i Q F O B Z t V + p l l m t A J G b H Q U U l S Z q N 8 c f I M X z h l i B N l X E n A C / X 3 R E 5 S a 6 d p 7 D p T A m O 7 6 s 3 F / 7 w w g + Q m y r n U G T B J l 4 u S T G B Q e P 4 / H n L D K I i p I 4 Q a 7 m 7 F d E w M o e B S q r g Q g t W X 1 0 m n U Q / 8 e v B w V W v e F n G U 0 R k 6 R 5 c o Q N e o i e 5 R C 7 U R R Q o 9 o 1 f 0 5 o H 3 4 r 1 7 H 8 v W k l f M n K I / 8 D 5 / A A M z k Q 8 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " w D A g F r s 8 w w b H 7 f P S Y Z F T v y U + S Y g = " > A A A B 8 n i c b V B N S w M x E M 3 W r 1 q / q h 6 9 B I v g q e w W Q Y 9 F L x 4 r 2 A / Y L i W b Z t v Q b B K S W a E s / R l e P C j i 1 V / j z X 9 j 2 u 5 B W x 8 M P N 6 b Y W Z e r A W 3 4 P v f X m l j c 2 t 7 p 7 x b 2 d s / O D y q H p 9 0 r M o M Z W 2 q h D K 9 m F g m u G R t 4 C B Y T x t G 0 l i w b j y 5 m / v d J 2 Y s V / I R p p p F K R l J n n B K w E n h Z J D 3 N T M a N 2 a D a s 2 v + w v g d R I U p I Y K t A b V r / 5 Q 0 S x l E q g g 1 o a B r y H K i Q F O B Z t V + p l l m t A J G b H Q U U l S Z q N 8 c f I M X z h l i B N l X E n A C / X 3 R E 5 S a 6 d p 7 D p T A m O 7 6 s 3 F / 7 w w g + Q m y r n U G T B J l 4 u S T G B Q e P 4 / H n L D K I i p I 4 Q a 7 m 7 F d E w M o e B S q r g Q g t W X 1 0 m n U Q / 8 e v B w V W v e F n G U 0 R k 6 R 5 c o Q N e o i e 5 R C 7 U R R Q o 9 o 1 f 0 5 o H 3 4 r 1 7 H 8 v W k l f M n K I / 8 D 5 / A A M z k Q 8 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " w D A g F r s 8 w w b H 7 f P S Y Z F T v y U + S Y g = " > A A A B 8 n i c b V B N S w M x E M 3 W r 1 q / q h 6 9 B I v g q e w W Q Y 9 F L x 4 r 2 A / Y L i W b Z t v Q b B K S W a E s / R l e P C j i 1 V / j z X 9 j 2 u 5 B W x 8 M P N 6 b Y W Z e r A W 3 4 P v f X m l j c 2 t 7 p 7 x b 2 d s / O D y q H p 9 0 r M o M Z W 2 q h D K 9 m F g m u G R t 4 C B Y T x t G 0 l i w b j y 5 m / v d J 2 Y s V / I R p p p F K R l J n n B K w E n h Z J D 3 N T M a N 2 a D a s 2 v + w v g d R I U p I Y K t A b V r / 5 Q 0 S x l E q g g 1 o a B r y H K i Q F O B Z t V + p l l m t A J G b H Q U U l S Z q N 8 c f I M X z h l i B N l X E n A C / X 3 R E 5 S a 6 d p 7 D p T A m O 7 6 s 3 F / 7 w w g + Q m y r n U G T B J l 4 u S T G B Q e P 4 / H n L D K I i p I 4 Q a 7 m 7 F d E w M o e B S q r g Q g t W X 1 0 m n U Q / 8 e v B w V W v e F n G U 0 R k 6 R 5 c o Q N e o i e 5 R C 7 U R R Q o 9 o 1 f 0 5 o H 3 4 r 1 7 H 8 v W k l f M n K I / 8 D 5 / A A M z k Q 8 = < / l a t e x i t > k ?1 < l a t e x i t s h a 1 _ b a s e 6 4 = " c P U x V 5 0 b / 9 g R V T j 2 / R n H S z B e + l 0 = " > A A A B 8 n i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J o Y x n B x E B y h L 3 N X L J k b / f Y 3 R P C k Z 9 h Y 6 G I r b / G z n / j J r l C E x 8 M P N 6 b Y W Z e l A p u r O 9 / e 6 W 1 9 Y 3 N r f J 2 Z W d 3 b / + g e n j U N i r T D F t M C a U 7 E T U o u M S W 5 V Z g J 9 V I k 0 j g Y z S + n f m P j 0 V r y S t m j u E P v M 8 f A a 6 R D g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " c P U x V 5 0 b / 9 g R V T j 2 / R n H S z B e + l 0 = " > A A A B 8 n i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J o Y x n B x E B y h L 3 N X L J k b / f Y 3 R P C k Z 9 h Y 6 G I r b / G z n / j J r l C E x 8 M P N 6 b Y W Z e l A p u r O 9 / e 6 W 1 9 Y 3 N r f J 2 Z W d 3 b / + g e n j U N i r T D F t M C a U 7 E T U o u M S W 5 V Z g J 9 V I k 0 j g Y z S + n f m P j 0 V r y S t m j u E P v M 8 f A a 6 R D g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " c P U x V 5 0 b / 9 g R V T j 2 / R n H S z B e + l 0 = " > A A A B 8 n i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J o Y x n B x E B y h L 3 N X L J k b / f Y 3 R P C k Z 9 h Y 6 G I r b / G z n / j J r l C E x 8 M P N 6 b Y W Z e l A p u r O 9 / e 6 W 1 9 Y 3 N r f J 2 Z W d 3 b / + g e n j U N i r T D F t M C a U 7 E T U o u M S W 5 V Z g J 9 V I k 0 j g Y z S + n f m P j 0 V r y S t m j u E P v M 8 f A a 6 R D g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " c P U x V 5 0 b / 9 g R V T j 2 / R n H S z B e + l 0 = " > A A A B 8 n i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J o Y x n B x E B y h L 3 N X L J k b / f Y 3 R P C k Z 9 h Y 6 G I r b / G z n / j J r l C E x 8 M P N 6 b Y W Z e l A p u r O 9 / e 6 W 1 9 Y 3 N r f J 2 Z W d 3 b / + g e n j U N i r T D F t M C a U 7 E T U o u M S W 5 V Z g J 9 V I k 0 j g Y z S + n f m P In a reference frame in which two colliding nucleons move along the z direction with 4-momenta P 1 and P 2 , a quark with 4-momentum k 1 and transverse momentum k ⊥1 annihilates with a parton with 4-momentum k 2 and transverse momentum k ⊥2 . A (virtual) photon (or Z) is produced with 4-momentum q and transverse momentum q T = k ⊥1 + k ⊥2 .
PDFs. Section 2.3 collects the perturbative ingredients of the factorised formula within the particular choice of the evolution scales adopted in this analysis. In section 2.4, we discuss how these perturbative ingredients are to be combined to achieve a given logarithmic accuracy of the resummation provided by TMD factorisation. In this context, we also review the different logarithmic-counting prescriptions used in the literature, highlighting the possible differences. Finally, in section 2.5 we motivate the introduction of a non-perturbative contribution that needs to be determined from data, and we discuss its particular functional form.

Drell-Yan cross section in TMD factorisation
In the inclusive Drell-Yan process two hadrons h 1 and h 2 with 4-momenta P 1 and P 2 , respectively, collide with center-of-mass energy squared s = (P 1 + P 2 ) 2 and produce a neutral vector boson γ * /Z with 4-momentum q and large invariant mass Q = q 2 . The vector boson eventually decays into a lepton and an antilepton with 4-momenta constrained by momentum conservation, q = l + l . The absolute value of the transverse momentum and the rapidity of the neutral boson (or, equivalently, of the lepton pair) are defined as where the z direction is defined by the hadronic-collision axis (see figure 1). We are specifically interested in the transverse-momentum distribution of the vector boson in the small-q T region (q T Q). In this regime, the (unpolarised) differential cross section factorises and can be expressed in terms of the (unpolarised) TMDs of the two JHEP07(2020)117 where α is the electromagnetic coupling and P is the phase-space reduction factor due to possible kinematic cuts on the final-state leptons (see appendix C). 1 The hard factor H represents the perturbative part of the hard scattering and depends on the hard scale Q and on the renormalisation scale µ. The summation over q in eq. (2.3) runs over the active quarks and antiquarks at the scale Q, and c q are the respective electroweak charges given by where e q , V q , and A q are respectively the electric, vector, and axial charges of the flavour q; V and A are the vector and axial charges of the lepton ; sin θ W is the weak mixing angle; M Z and Γ Z are mass and width of the Z boson. The second line of eq. (2.3) displays the convolution of the TMDs f q 1 and fq 1 of the hadrons h 1 and h 2 , respectively. It describes the annihilation of a quark q, with longitudinal momentum fraction x 1 = Qe y / √ s and transverse momentum k ⊥1 , with the corresponding antiquarkq, with longitudinal momentum fraction x 2 = Qe −y / √ s and transverse momentum k ⊥2 . In the annihilation, the momentum conservation is guaranteed by the presence of δ (2) k ⊥1 + k ⊥2 − q T (see figure 1).
As a consequence of renormalisation and of the removal of the rapidity divergences [27], TMDs acquire a dependence on the renormalisation scale µ and on the so-called rapidity scale ζ. We will discuss our choice for these scales in section 2.3. Here, we just remark that the rapidity scales ζ 1 and ζ 2 in eq. (2.3) must obey the kinematic constraint ζ 1 ζ 2 = Q 4 .
It is convenient to rewrite the convolution in the conjugate position space by using the Fourier transform of each TMD, defined as 2 In the presence of cuts on single lepton variables, an additional parity-violating term contributes to the cross section [26]. However, in appendix C we argue that this contribution is negligible in the experimental conditions considered in this paper. 2 For simplicity, in the rest of the paper we will refer to the bT -dependent functionf1 as to TMD but understanding that this is in fact the Fourier transform of the actual TMD f1. Note that in ref. [21] the variable ξT was used in place of bT . The reason was to avoid confusion with the impact parameter used in the GPD literature for which the symbol bT is typically used. In this paper, we decided to use bT as it is more common in the TMD, qT -resummation, and SCET literature but keeping in mind that this is not the impact parameter but the Fourier conjugate variable of qT . Finally, we notice that in ref. [21] the Fourier transform was defined with an extra 1/(2π) factor.

JHEP07(2020)117
where b T is the absolute value of the vector b T (b T = |b T |). By using eq. (2.7), we can rewrite the convolution of TMDs as where J 0 is the 0-th order Bessel function of the first kind that has the following integral representation By inserting eq. (2.8) into the cross section in eq. (2.3), we finally get which is the formula actually implemented in our analysis of Drell-Yan data.

TMD evolution and matching
In eq. (2.10), the dependence of the TMDsf on the scales µ and ζ arises from the removal of the ultraviolet and rapidity divergences in their operator definition. Each dependence is controlled by an evolution equation: where γ is the anomalous dimension of the Renormalisation Group (RG) evolution in µ, and K is the anomalous dimension of the Collins-Soper evolution in √ ζ [28]. Notice that, for brevity, we have dropped the flavour index q andq. Moreover, since in this section we will only be concerned with the dependence off 1 on the scales µ and ζ, we will also temporarily drop the dependence on x and b T . In addition to the evolution equations in eq. (2.11), the rapidity anomalous dimension K obeys its own RG equation: where γ K is known as cusp anomalous dimension. Since the crossed double derivatives of f 1 must be equal, using eqs. (2.11) and (2.12) we also get

JHEP07(2020)117
Using the point ζ = µ 2 as a boundary condition, the solution of this differential equation is where γ F (α s (µ)) ≡ γ(µ, µ 2 ). If the TMDf 1 is known at some starting scales µ 0 and ζ 0 , the solution of the evolution equations in eq. (2.11) readŝ where the so-called Sudakov form factor R accounts for the perturbative evolution off 1 and it is defined as (2.16) We note that eq. (2.16) can be implemented in various ways [29][30][31][32]. In this work, we follow the standard approach described in [27]. Moreover, we calculate all ingredients involved in eq. (2.16) by adopting a fully numerical approach.
An important property of the TMDf 1 is that at small values of b T it can be matched onto the collinear PDF f 1 . Reinstating for clarity the x and b T dependence and introducing the matching coefficient function C, we can write 3 Then, the actual evolved TMD becomeŝ

Perturbative content
In order to use eq. (2.18) in phenomenological applications, we need to define the values of both the initial and final pairs of scales, (µ 0 , ζ 0 ) and (µ, ζ). It turns out that in the MS renormalisation scheme there exists a particular scale, with γ E the Euler constant, such that the rapidity anomalous dimension K and the matching coefficient C computed at µ 0 = √ ζ 0 = µ b admit a pure perturbative expansion free of explicit logarithms of the scales. Therefore, µ b provides a natural choice for µ 0 and √ ζ 0 . The final renormalisation scale µ must match the one used in the hard factor H in eq. (2.10). Therefore, µ has to be of order Q for avoiding large logarithms in H: we choose µ = Q. Any variation of µ with respect to this choice can be accounted for by expanding the solution of the RG equation for the strong coupling α s . The rapidity scales ζ 1 and ζ 2 in eq. (2.10) are bound to comply with ζ 1 ζ 2 = Q 4 . Therefore, the natural choice is ζ 1 = ζ 2 = Q 2 . However, we stress that any choice that fulfils this constraint leads to the same cross section. In fact, from eq. (2.16) it should be evident that the evolution factors R entering the two TMDs in eq. (2.10) combine in such a way that the result only depends on the product ζ 1 ζ 2 .
After choosing the scales, we discuss the perturbative ingredients that result from this particular choice. We first consider the hard function H. Up to two-loop accuracy, its perturbative expansion is The coefficients H (n) can be read off from, e.g., ref. [33]. When going beyond O(α 2 s ), the hard function acquires a non-trivial flavour structure (see, e.g., ref. [34]). As a consequence, H should in principle be moved inside the flavour sum in eq. (2.10). However, in the present analysis we do not consider corrections beyond O(α 2 s ) and eq. (2.10) is appropriate. Next, we consider the matching function C introduced in eq. (2.17). By making the flavour and x dependences explicit, the C have the following perturbative expansion The coefficient functions C (n) ij up to n = 2 have been computed in refs. [35,36]. They have been reported also in ref. [34], where the authors have verified the consistency of the results. The calculation of the O(α 3 s ) corrections to the quark matching functions appeared very recently in ref. [37].
As for the anomalous dimensions K, γ F , and γ K in the Sudakov form factor in eq. (2.16), their perturbative expansions read, respectively, The coefficients K (n) are listed up to n = 3 in ref. [36] and up to n = 2 in ref. [34]. They differ by a factor −2 due to a different definition of K. Also the coefficients γ (n) F are given in refs. [34,36] up to n = 2, and they differ by a minus sign due to a different definition of the anomalous dimension. Finally, the coefficients γ (n) K were originally computed in ref. [38] and are also given in refs. [34,36] up to n = 2, where they differ by a factor 2. The coefficient γ K has been recently computed in refs. [39][40][41].

Logarithmic ordering
In this section, we discuss how to combine in a consistent way the perturbative ingredients of eqs. (2.20)-(2.22) for the computation of the cross section in eq. (2.10) (see also refs. [42,43]).
As is well known, TMD factorisation provides resummation of large logarithms of Q/q T or, equivalently, of Q/µ b . The resummation is implemented in the Sudakov form factor R in eq. (2.16) whose perturbative expansion reads Because of the inner sum running up to 2n, eq. (2.23) exposes the double-logarithmic nature of the resummation. This structure can be traced back to the evolution equations in eq. (2.11) that resum two different categories of logarithms. However, our particular choice of the scales (µ 0 = √ ζ 0 = µ b and µ = √ ζ = Q) makes the two categories to coincide, producing up to two logarithms for each power of α s . Consequently, eq. (2.23) must include all powers of α s if the scales are such that α s L 2 1.
The expansion (2.23) can be rearranged to define a logarithmic ordering as where [k/2] is the integer part of k/2. According to this definition, the term k = 0 in eq. (2.25) gives the leading-logarithmic (LL) approximation, the term k = 1 gives the next-to-leading-logarithmic (NLL) approximation, and so on. Multiplication of R N k LL by a power p of α s gives 27) where the symbol ∼ means that the left-and right-hand sides have the same logarithmic accuracy. This step is relevant because in the cross section the Sudakov form factor, eq. (2.25), can be multiplied by some power of α s originating from the hard factor H and/or the matching functions C. Equation (2.27) states that, at the cross section level, the inclusion of an additional power of α s in the perturbative expansion of H and/or C implies a contribution two orders higher with respect to the leading term in the logarithmic expansion. For example, at LL and NLL accuracy the functions H and C can be computed JHEP07(2020)117 at O(1), at NNLL and N 3 LL they need to include the O(α s ) corrections, and so on. This logarithmic counting is illustrated in the left panel of figure 2: the diagonal bands represent the terms included in each R N k LL , with H (n) the perturbative coefficients of either H or C or a combination of the two. The counting discussed above generally applies to any process whose amplitude factorises in the appropriate limit, such as DY in the q T Q limit (TMD factorisation). However, in the specific case of DY (i.e., inclusive with respect to soft-collinear QCD radiation) also the phase space for the emission of n real particles in b T space factorises (see, e.g., ref. [44]). This feature, along with the factorisation of the amplitude in the q T Q limit, allows one to exponentiate soft-collinear emissions such that the Sudakov form factor can be written in the following general form (see, e.g., ref. [45] where the functions g (i) are such that g (i) (0) = 0. As compared to the general counting in eq. (2.23), exponentiation relates all the terms in eq. (2.23) of the type α n s L m with n + 1 < m ≤ 2n to the lower-order terms. In eq. (2.28), the logarithmic counting is performed at the level of the argument of the exponential. In this context, the terms Lg (1) , g (2) , α s g (3) , etc., resum, respectively, the LL contributions α n s L n+1 , the NLL contributions α n s L n , the NNLL contributions α n s L n−1 , etc.. Contrary to eq. (2.23), this counting is driven by the condition α s L 1. This extends the validity of the resummed result (truncated at a given level: NLL, NNLL, etc.) to larger values of L (smaller values of q T /Q).
The logarithmic counting applied to the argument of the exponential is equivalent to consider the logarithm of the cross section [33]. In fact, neglecting for simplicity the matching functions, we schematically have

JHEP07(2020)117
The logarithm of H can be expanded as The first term α s H (1) contributes to the tower α n s L n−1 , that is the NNLL contribution. The second term α 2 s H (2) − H (1)2 /2 contributes to the α n s L n−2 tower, thus to the N 3 LL contribution. The same counting applies to the matching functions C. The conclusion is that including O(α s ) contributions in H and C implies introducing NNLL corrections, O(α 2 s ) contributions in H and C contribute to N 3 LL accuracy, and so on. A graphical representation of this counting is sketched in the right panel of figure 2. Again, the bands represent the logarithmic towers, while H (n) are the appropriate coefficients of the expansion of either ln H or ln C or a combination. This logarithmic counting has been used in several papers (see, e.g., refs. [10,33,46,47]). In this work, we will simply denote this counting with the acronyms NLL, NNLL, and so on, and for convenience we will refer to it as to "standard counting". A slightly different counting has also been widely used in the literature (see, e.g., refs. [42,[48][49][50][51]). Expanding the Sudakov form factor (2.28) and multiplying it by the expansion of the hard function in eq. (2.20), we obtain for the cross section where the rightmost term stems from the combination of the first-order terms α s H (1) and Lg (1) in both expansions. As it is clear from the previous discussion, this term has the same form α n s L n as g (2) . Then one can argue that NLL accuracy requires the inclusion not only of g (2) but also of H (1) [48]. This argument works to all orders: at any given logarithmic accuracy, it prescribes to include one more order in the perturbative expansion of H (and/or C) with respect to the standard counting. We will refer to this counting as the to "primed counting", denoting it as NLL , NNLL , and so on. The apparent contradiction between the standard and primed countings is resolved by observing that the first term of the perturbative expansion of α s Lg (1) is proportional to α 2 s L 2 . When considering the general expansion of the cross section given in eqs. (2.25)-(2.27), a term proportional to α 2 s L 2 is of the form α n s L 2n−2 and thus belongs to the NNLL tower. This is formally subleading with respect to the NLL accuracy determined by the g (2) term in the exponent.
Accurate predictions over a wide range in q T require matching resummed calculations (valid at q T Q) to the corresponding fixed-order calculation (valid at q T Q). In this context, the primed ordering turns out to be more advantageous. Indeed, the accuracy of a fixed-order calculation is measured in terms of powers of α s relative to the leading term. In order to produce a Z boson with large q T , it is necessary to produce (at least) a second object with large transverse momentum against which the Z boson recoils, i.e., a jet. As a consequence, the leading-order (LO) contribution to the q T distribution of the Z at fixed order is O(α s ). The NLL prescription correctly reproduces the small-q T limit of the LO fixed-order calculation. It is then possible to realise the matching in an additive way by

JHEP07(2020)117
Accuracy H and C K and γ F γ K PDF and α s evolution combining the NLL resummed calculation with the LO fixed-order one (NLL + LO). The procedure can be extended to higher orders: NNLL + NLO, N 3 LL + NNLO, and so on. Conversely, in the standard counting the matching to the LO fixed-order calculation requires to go further to NNLL accuracy (NNLL + LO), combining in this way a rather accurate calculation at small q T with a poorly accurate calculation at large q T . At higher orders one has N 3 LL + NLO, N 4 LL + NNLO, and so on. We remark that other forms of matching can be used to overcome the limitation of the standard counting [33,52,53]. Finally, table 1 summarises the perturbative ingredients to be used for a consistent computation of the cross section in eq. (2.10) for both the standard and the primed countings. The numbers in table 1 give the maximum power of α s at which the corresponding quantity is to be computed, while the last column reports the corresponding accuracy in computing the evolution of the collinear PDFs and of the coupling α s . 5 In this analysis, we have used the PDF sets of the MMHT2014 family [54] at the appropriate perturbative order accessed through the LHAPDF interface [55].

Non-perturbative content and its parameterisation
In the previous section, we noticed that in the MS scheme the rapidity evolution kernel K and the matching functions C can be made free of logarithms of the scales by introducing the natural scale µ b defined in eq. (2.19). Consistently, in the perturbative expansion of K (see first line of eq. (2.22)) and C (see eq. (2.21)) the strong coupling α s must be computed at µ b . For large values of b T , µ b becomes small such that α s (µ b ) may potentially become very large and eventually diverge when µ b reaches the Landau pole at Λ QCD . As a matter of fact, the integral in eq. (2.10) does require accessing large values of b T . It is then necessary to regularise this divergence by introducing a prescription that avoids integrating over the Landau pole. Different possibilities are available (see, e.g., refs. [53,56]). In this paper, we adopt the prescription originally proposed in ref. [57]: we introduce the arbitrary JHEP07(2020)117 parameter b max that denotes the maximum value of b T at which perturbation theory is considered reliable. Hence, b max must be such that Moreover, we also want to prevent µ b from becoming much larger than the hard scale Q (µ b Q). Despite not strictly mandatory (especially when considering only small values of q T ), this feature makes it possible to expand the cross section integrated in q T , with the lowest-order term reproducing the lowest-order collinear result [58]. To this end, we define and introduce a monotonic function b * (b T ) with the following asymptotic behaviours In this analysis, we adopt for b * (b T ) the same functional form chosen in ref. [21] that guarantees a smooth and rapid convergence towards the asymptotic limits: Now, we simply write the TMDf 1 aŝ (2.36) This separation effectively defines f NP . The advantage is that, due to the behaviour of b * (b T ) for large values of b T ,f 1 (x, b * (b T ), µ, ζ) remains in the perturbative region. The non-perturbative contributions are instead confined into f NP , that has to be determined through a fit to experimental data. However, using eq. (2.36), we can work out some general properties of f NP . First, f NP does not depend on the renormalisation scale µ. To see this, using eqs. (2.15) and (2.16) ). The dependence on µ evidently cancels in the ratio. In addition, for large values of b T µ b * saturates to some minimal value while µ b becomes increasingly JHEP07(2020)117 small. As a consequence of this departure between µ b * and µ b , as well as between √ ζ and µ b , the exponential in eq. (2.37) tends to be suppressed, and so does f NP . Conversely, as b T becomes small b * approaches b min . Using the definition in eq. (2.33), it follows that µ b * saturates to Q while µ b becomes larger and larger. In this limit, we have [58] where p is some positive number. Since TMD factorisation applies to leading-power in q T /Q, we can neglect the power suppressed contribution such that f NP → 1 for b T → 0.
It is important to stress that the separation between perturbative and non-perturbative components of a TMD is arbitrary and depends on the particular choice of b * (or in general on the prescription used to regularise the Landau pole). For any given choice, only the combination in eq. (2.36) is meaningful, and it is misleading to refer to f NP as to the non-perturbative part of TMDs in a universal sense. Following the requirements discussed above, we parameterise f NP as with Q 0 = 1 GeV and with the g 1 (x) and g 1B (x) functions given by (2.40) There are a total of 9 free parameters (λ, g 2 , g 2B , N 1 , σ, α, N 1B , σ B , α B ) to be determined from data. Apart from the logarithmic dependence on ζ, the functional form (2.39) is motivated by empirical considerations. The first line parameterises the "intrinsic" TMD nonperturbative contribution and it only depends on x and b T . The second line accounts for the non-perturbative correction to the perturbative evolution. Therefore, it only depends on b T (on top of the known dependence on ζ).
The intrinsic contribution is a combination of a q-Gaussian (or Tsallis) distribution (first term) and a standard Gaussian distribution (second term). The q-Gaussian has a larger tail than the standard Gaussian, meaning that it gives a bigger contribution to the TMD at small transverse momentum. We found that this combination is able to reproduce the behaviour at very small q T of the experimental distributions from the lowest to the highest energies considered in our analysis.
The functions g 1 and g 1B in eq. (2.40) are related to the width of the TMD distribution. They are expected to depend on x on the basis of model calculations (see ref. [59] and references therein) and more generally from Lorentz invariance constraints on the proton JHEP07(2020)117 light-front wave functions (see, e.g., the discussion in ref. [60]). To best describe experimental data, we found it necessary to have wider TMDs at intermediate x. A log-normal dependence of g 1 and g 1B allowed us to properly describe the datasets differential in the boson rapidity y. In fact, as we will show below, the x dependence of f NP is almost entirely determined by the ATLAS datasets, the only ones differential in y. Our present results are quite different from the ones obtained through fits to semi-inclusive DIS data [21]. We expect that the addition of further datasets from DIS experiments [61,62] will provide more sensitivity to the x dependence and possibly lead to different results.
The non-perturbative components of the TMDs could depend also on flavour [17,25,63]. However, in this work we refrain from including such dependence since DY data are not very sensitive to it. We stress that the fact that we can achieve a good description of data does not exclude the presence of a flavour dependence, which is actually expected on the basis of model calculations [64][65][66][67][68][69], lattice QCD studies [70], and also if QED corrections are taken into account [71,72]. Higher sensitivity to flavour dependence may be provided again by semi-inclusive DIS data with different targets and final-state hadrons and possibly by W -boson production data [73].
Concerning the b T dependence of the non-perturbative evolution in the second line of eq. (2.39), we have used a customary quadratic term [4,8,12,74] with an additional quartic term. The latter contribution appears to be useful to reproduce the energy evolution displayed by the data. Other choices of the functional form have been discussed in, e.g., refs. [20,[75][76][77]. This contribution could be also determined using lattice QCD [78].

Experimental data
In this section we describe the experimental data included in this analysis. We considered q T distributions in DY production from a variety of datasets. Some of these were already included in the analysis of ref. [21], i.e. data from: E605 [79], E288 [80], CDF Run I [81] and Run II [82], and D0 Run I [83] and Run II [84]. We refer the reader to ref. [21] for more details. The new datasets included in the present analysis are: • forward Z-production data from the LHCb experiment at 7 [85], 8 [86], and 13 [87] TeV, • Z-production data from the CMS experiment at 7 [88] and 8 [89] TeV, • Z-production data differential in rapidity from the ATLAS experiment at 7 [88] and 8 [90] TeV, • off-peak (low-and high-mass) DY data from the ATLAS experiment at 8 TeV [90], • preliminary Z-production data from the STAR experiment at 510 GeV. 6 6 We thank the STAR Collaboration for providing us with the data. Finally, we originally considered also measurements from the PHENIX experiment at the center-of-mass energy of 200 GeV [91]. However, due to the cut on q T /Q discussed below, only two data points from this dataset would be included in the fit. Therefore, we decided to exclude it. The breakdown of the entire dataset included in our analysis is reported in table 2. For visualisation purposes, in figure 3 we show the kinematic coverage of each datasets in the x 1 vs. x 2 plane, with x 1,2 = Qe ±y / √ s. The shaded areas are determined considering the corresponding ranges in Q and y, and the center-of-mass energy √ s. 7 As expected, the lowerenergy experiments (E605, E288, and STAR) are placed in the large-x region (x 0.1). Particularly important are the new (preliminary) STAR measurements that cover a kinematic region that is scarcely populated. The Tevatron experiments, CDF and D0, cover a particularly wide kinematic region at intermediate values of x. These experiments (except D0 Run II with muons) provide data extrapolated over the full range in rapidity y, thus extending across the full available phase space. Finally, the LHC experiments (LHCb, CMS, and ATLAS) are placed at lower values of x. The LHCb datasets are in a region in which x 1 is particularly small and x 2 particularly large: this is due to the fact that the data is taken in the forward region, 2 < y < 4.5. The ATLAS datasets are binned in rapidity and thus are expected to be particularly sensitive to the x dependence of the TMDs. Indeed, we will show below that the x dependence of TMDs is mostly constrained by these datasets.
Since our analysis is based on the TMD factorisation formula in eq. (2.10), only data at small q T can possibly be described. Hence, we impose a cut to exclude measurements JHEP07(2020)117 Total 353 ------ Table 2. Breakdown of the datasets included in this analysis. For each dataset, the table includes information on: the number of data points (N dat ) passing the nominal cut on q T /Q, the observable delivered, the center of mass energy √ s, the range(s) in invariant mass Q, the angular variable (either y or x F ), possible cuts on the single final-state leptons, and the public reference (when available). The total number of data points amounts to 353. Note that for E605 and E288 400 GeV we have excluded the bin in Q containing the Υ resonance (Q 9.5 GeV).

JHEP07(2020)117
with large q T by requiring q T /Q < 0.2. Since the measurements are delivered in transversemomentum bins [q T,min : q T,max ] integrated over some range in invariant mass [Q min : Q max ], the cut is conservatively imposed on the ratio q T,max /Q min . The second column in table 2 reports the number of data points (N dat ) for each dataset that pass this cut: the total number of points included in our analysis is 353.
An important feature of all the new datasets listed above is that the cross sections are given within a certain fiducial region. In particular, kinematic cuts on transverse momentum p T and pseudo-rapidity η of the final-state leptons are enforced. The values of the cuts are reported in the next-to-last column of table 2. Our predictions are corrected by means of the phase-space reduction factor P introduced in eq. (2.10), which takes into account these cuts. Details concerning the calculation of P are given in appendix C.
As evident from the "Observable" column of table 2, experimental cross sections are released in different forms. In addition, some of them are normalised to the total (fiducial) cross section while others are not. In our analysis, we expressed all the absolute cross sections in terms of the observable given in eq. (2.10) (details on the transformations between different observables can be found in ref. [21]). When necessary, the total cross section σ required to normalise the differential cross sections is computed using DYNNLO [94,95] with the MMHT2014 collinear PDF sets [54], taking into account the selection cuts and consistently with the perturbative order of the differential cross section. More precisely, the total cross section is computed at LO for NLL accuracy, at NLO for NLL' and NNLL, and at NNLO for NNLL' and N 3 LL. The values of the total cross sections at different orders are reported in table 3. We stress that in this analysis no additional normalisations have been applied, with the consequence that both the shape and the normalisation of the experimental distributions have an impact on the fit.
Most of the considered experimental datasets are released with a set of uncorrelated and correlated uncertainties. As already pointed out in ref. [16], a proper treatment of the experimental uncertainties is crucial to achieve a reliable extraction of TMDs. In other words, the χ 2 , which quantifies the agreement between data and predictions and is minimised during the fit, has to be computed taking into account the nature of the various uncertainties. Particular care has to be taken with the (correlated) normalisation uncertainties. As is well known, an inappropriate description of normalisation uncertainties may lead to underestimate the predictions: that is the so-called D'Agostini bias [96,97]. Different prescriptions have been devised to avoid this problem [98]: in this analysis we adopt the so-called iterative t 0 -prescription [99].
In the presence of correlated uncertainties, the χ 2 can be split as [98] where χ 2 D has an uncorrelated structure (diagonal) while χ 2 λ is a penalty term related to the presence of correlations (see, e.g., appendix B of ref. [16]). For the computation of χ 2 D , theoretical predictions are properly shifted to take into account the effect of the correlated uncertainties. In fact, shifted predictions are a better proxy for visual comparisons to experimental data. Therefore, in the following it is understood that all plots will display shifted predictions.  [54] and required for the computation of the normalised differential cross sections at the different perturbative orders.

JHEP07(2020)117
A further important aspect is the use of collinear PDFs. In order to extract f NP defined in eq. (2.36), it is necessary to assume a given set of collinear PDFs (MMHT2014 in our case). PDF uncertainties reflect the experimental uncertainty of the dataset used for their extraction. It is therefore natural to attribute an experimental nature to this uncertainty and include it in the calculation of the χ 2 . To do so, we computed the PDF errors as relative to the central value 8 and included them in the experimental covariance matrix as uncorrelated uncertainties. The propagation of the resulting experimental uncertainty into the fitted TMDs is achieved through Monte Carlo sampling. Specifically, we generate N rep ( 200) replicas of the original dataset taking into account all the uncertainties and then perform a fit on each single replica. The resulting ensemble of distributions can be used to compute central values and uncertainties as averages and correlations, respectively.
A final remark concerns the integration over the final-state phase space. The basic quantity to be compared to data is where the ranges [y min : y max ], [Q min : Q max ], and [q T,min : q T,max ] define the phase-space integration region and the integrand is given in eq. (2.10). In order to speed up the JHEP07(2020)117 numerical computation of the theoretical predictions, the integration over the bins in q T and Q is often performed approximating the q T -bin integral with its central value and using the narrow-width approximation for the integral over Q around the Z peak. We stress that in this analysis the integrals in eq. (3.2) are computed exactly. While the integrals over y and Q do need to be computed numerically, the integral over q T can be performed (semi)analytically exploiting a property of the Bessel functions J n (see appendix B). This greatly reduces the amount of numerical computations.

Results
In this section, we present the results of our extraction of unpolarised TMDs from a comprehensive set of DY data (see section 3). In section 4.1, we present the quality of the fit at N 3 LL, the best accuracy we can presently reach. In section 4.2 we discuss the TMDs extracted from the nominal fit. In section 4.3, we discuss the convergence of the perturbative corrections. In section 4.4, we focus on the x dependence of the TMDs and we argue that it is mostly constrained by the y-differential ATLAS cross sections. Finally, in section 4.5, we assess the range of validity of TMD factorisation by considering the fit quality as a function of the cut on q T /Q.

Fit quality
In this section, we discuss the quality of the reference fit at N 3 LL with cut q T /Q < 0.  The mean value provides a democratic representative of the ensemble. Other choices are possible, such as the median or the mode of the ensemble. In fact, only the full ensemble of replicas carries the full statistical information. However, the reason for using eq. (4.1) is that quantifying the goodness of our fit becomes easier, as it will be clear in the following. Table 4 reports the breakdown of the χ 2 s normalised to the number of data points, N dat , for each dataset. The uncorrelated (χ 2 D ) and the correlated (χ 2 λ ) contributions to the total χ 2 (see eq. (3.1)) are also reported. The global χ 2 is shown at the bottom of the table.
The value of the global χ 2 is very close to one (1.02), indicating that the fit is able to describe measurements over a wide energy range, from the low-energy fixed-target datasets to the LHC ones. It is important to stress that a substantial contribution to the global χ 2 is given by the correlated penalty term, χ 2 λ /N dat = 0.14. This highlights the importance of a correct treatment of the correlated uncertainties. More specifically, the systematic shifts induced by correlations are often large, indicating that the fit does need to adjust the predictions within the experimentally correlated ranges.
Concerning the single experiments, we observe that the low-energy data (E605, E288, and STAR) have generally lower χ 2 s than the Tevatron (CDF and D0) and LHC (LHCb,

JHEP07(2020)117
CMS, and ATLAS) high-energy data. This is mostly due to the fact that the experimental uncertainties of the former are typically larger than the latter. In particular, the lowenergy data are affected by large normalisation (correlated) uncertainties. Consequently, the relative importance of the correlated contribution χ 2 λ to the total χ 2 is generally larger for the low-energy datasets than for the high-energy ones.
It is interesting to comment on the quality of the fit to the new datasets from RHIC and the LHC that were not included in the analysis of ref. [21] (see section 3). The preliminary measurements from STAR have a χ 2 equal to 0.836. This is particularly encouraging because, as clear from figure 3, this dataset covers a scarcely populated kinematic region and shows no tension with other data. Also the LHC datasets extend the kinematic coverage of the DY data considered in ref. [21]. These measurements are particularly precise and thus very effective in constraining TMDs. We observe that the LHCb datasets are very nicely described with χ 2 s that never exceed 1.3. The CMS data, despite having slightly larger χ 2 , are also well described. The two CMS datasets provide only eight points in total and thus their impact on the fit is modest. The ATLAS datasets, amongst the LHC ones, are by far the most abundant. We observe that the ATLAS 8 TeV datasets are well described, except for the first two low-rapidity bins. The 7 TeV ones present larger values of χ 2 , above 2. Given the extremely high precision of these datasets, even small effects (e.g., power corrections) could give a significant contribution to χ 2 in these conditions. We consider it already a success to obtain a value of χ 2 for these datasets that does not affect too much the global χ 2 . We note that a key feature of these datasets (except the off-peak ones) is that they are differential in the vector-boson rapidity y. As we will see in section 4.4, the x dependence of f NP plays a crucial role in improving the χ 2 .
In order to provide a visual assessment of the fit quality, figure 4 displays the data/theory comparison for a representative selection of datasets. We remind the reader that in each plot theoretical predictions are appropriately shifted to account for correlated uncertainties [16], while the experimental error bars are given by the sum in quadrature of the uncorrelated uncertainties. The upper panel of each plot shows the absolute q T distribution, while the lower panel shows the ratio to data. The plots in the upper row of figure 4 refer to one invariant-mass bin of E605 and CDF Run II already considered in ref. [21]. The remaining plots refer to some of the new datasets, namely STAR, LHCb 8 TeV, ATLAS 8 TeV on-peak at 1.6 < |y| < 2, and ATLAS 8 TeV off-peak at 116 GeV < Q < 150 GeV. As expected, there is a very good agreement between data and theory, for both the old and the new datasets. Finally, it is interesting to observe that the uncertainties of the upper and middle rows of figure 4 are larger than those in the two lower rows. This is due to the fact that the ATLAS distributions are normalised to the total cross section leading to a cancellation of some uncertainties, such as those due to luminosity and collinear PDFs.

TMD distributions
We discuss now the TMD distributions extracted from our reference N 3 LL fit. We stress once again that only the combination in the r.h.s. of eq. (2.36) is meaningful.
In order to assess the sensitivity of the experimental dataset to f NP , it is interesting to look at the values of the free parameters obtained from the fit. In table 5

JHEP07(2020)117
Parameter Value g 2 0.036 ± 0.009 0.012 ± 0.003 Table 5. Average and standard deviation over the Monte Carlo replicas of the free parameters fitted to the data and graphical representation of the correlation matrix.
parameter over the Monte Carlo replicas, along with the respective standard deviation, is reported. All parameters are well constrained. 9 It is interesting to observe that the parameter λ, that measures the relative weight of Gaussian and q-Gaussian in eq. (2.39), is close to 0.5 indicating that these contributions weigh approximately the same. Concerning the values of the parameters g 2 and g 2B associated to the non-perturbative contribution to TMD evolution, we find that the coefficient g 2B of the quartic term is small but significantly different from zero. This seems to suggest that higher-power corrections to the commonly assumed quadratic term g 2 may be required by the data. Further insight concerning the appropriateness of the functional form in eqs. (2.39)-(2.40) can be gathered by looking at the statistical correlations between parameters. In the right panel of table 5, we show a graphical representation of the correlation matrix of the fitted parameters. The first observation is that (off-diagonal) correlations are generally not very large. There is however one exception, i.e. the parameters σ and λ seem to be strongly anti-correlated. This may indicate that the interplay between q-Gaussian and Gaussian may be significantly x dependent. We leave a deeper study of this feature to a future publication.
To conclude this section, in figure 5 we show the down-quark TMD at µ = √ ζ = Q = 2 GeV (left plot) and 10 GeV (right plot) as a function of the partonic transverse momentum k ⊥ for x = 0.001, 0.1, 0.3. The 1-σ uncertainty bands are also shown. As expected, TMDs are suppressed as k ⊥ grows and the suppression becomes relatively stronger as Q increases.

Perturbative convergence
In the previous section we discussed the quality of our fit at N 3 LL, which is the best accuracy presently available. In this section we show how the inclusion of perturbative corrections is crucial to achieve a better description of the experimental data. To this end, we performed fits at NLL , NNLL, and NNLL (see section 2.  Table 6. Values of the global χ 2 of the fits at NLL , NNLL, NNLL , and N 3 LL accuracy. N 3 LL fit. We did not consider LL and NLL accuracies because in both cases the description of the data is very poor (χ 2 20). Table 6 reports the values of the global χ 2 for each of the four accuracies considered. In order to appreciate the significance of the differences, 10 we have reported the absolute values of the χ 2 without dividing by the number of data points N dat . Figure 6 shows a graphical representation of table 6. The global quality of the fit improves significantly as the perturbative accuracy increases. In addition, figure 6 shows that the convergence rate decreases when going to larger perturbative orders. On the one hand, we conclude that it is necessary to include higher perturbative corrections to obtain a good description of the data and that N 3 LL corrections are still significant. On the other hand, it appears that the perturbative series is nicely converging and N 3 LL accuracy seems appropriate within the current experimental uncertainties. In order to quantify the numerical impact of higher-order corrections, in figure 7 we compare the predictions for all the available perturbative orders to the ATLAS 8 TeV data in the bin 66 GeV < Q < 116 GeV and 1.6 < |y| < 2. This plot shows how the inclusion of higher-order corrections improves the shape of the predictions, particularly around the peak region.

Reduced dataset and x dependence
The non-perturbative function f NP , eq. (2.36), accounts for the large-b T behaviour of TMDs. It is in general a function of b T , ζ, and x. While the asymptotic dependence on b T is driven by first-principle considerations (see section 2.5) and the evolution with ζ is determined by the Collins-Soper equation (2.11), the dependence on x is totally unknown. Moreover, a direct access to the x dependence is particularly difficult to achieve because it requires cross-section data finely binned in rapidity y. In the dataset considered here, only the ATLAS experiment delivers data differential in rapidity. Therefore, one would expect that these datasets provide most of the sensitivity to the x dependence of TMDs.
In order to test this conjecture, we employed a particularly simple x-independent parameterisation of the non-perturbative function:  Table 7. The values of the global χ 2 normalised to the number of data points N dat from the fit to the full dataset and to a reduced dataset without the y-differential ATLAS datasets, both using the parameterisation in eq. (4.2). For completeness, we also report the best-fit values of the parameters g 1 and g 2 .
with two free parameters, g 1 and g 2 , and Q 2 0 = 1.6 GeV 2 (inspired by the pioneering work of Davies, Webber, and Stirling. [4]). Using eq. (4.2) we first performed a fit at N 3 LL to the full dataset. Then we excluded the ATLAS datasets differential in rapidity (but we kept the off-peak ATLAS 8 TeV datasets because inclusive in rapidity). The resulting χ 2 s normalised to the number of data points are reported in table 7. For completeness, we also show the best-fit values of the parameters g 1 and g 2 .
Firstly, the χ 2 of the fit to the full dataset using eq. (4.2) (1.339) is significantly larger than that obtained using the parameterisation in eqs. (2.39)-(2.40) (1.020). This suggests that an x-dependent f NP is required to obtain a good description of the data. Secondly, the χ 2 of the fit without the y-differential ATLAS data comes out to be particularly low (0.895). We conclude that at N 3 LL accuracy the x dependence of the TMDs extracted from the currently available DY data is mostly constrained by the ATLAS data differential in the boson rapidity y. We note however that the agreement with the very precise ATLAS data may be influenced also by other small corrections (e.g. power corrections).

Dependence on the cut on q T /Q
As discussed in section 2, our analysis is based on TMD factorisation whose validity is restricted to the region q T Q. As a consequence, we consider only measurements that respect this constraint. More precisely, we require that the maximum value of the ratio q T /Q for a point to be included in the fit be 0.2 (see section 3). Despite this particular value seems to be generally recognised in the literature (see, e.g., ref. [15]), it is interesting to study how the global description of the dataset changes by varying this cut. This will help us assess more quantitatively the validity range of TMD factorisation. Figure 8 displays the behaviour of the global χ 2 /N data for the N 3 LL fit as a function of the q T /Q cut ranging between 0.1 and 0.28 in steps of 0.02. As expected, the quality of the fit tends to degrade as the cut on q T /Q increases. Of course, it is impossible to draw a line between validity and non-validity regions. However, this study gives a quantitative justification for choosing the value 0.2 for the q T /Q cut.

Conclusions
In this paper we presented an extraction of TMDs from Drell-Yan data accurate up to N 3 LL. The dataset used in this analysis includes low-energy data from FNAL (E605 and E288) and RHIC (STAR) and high-energy data from Tevatron (CDF and D0) and the LHC (LHCb, CMS, and ATLAS), for a total of 353 data points.
The fit was performed with a proper treatment of the experimental uncertainties, which were propagated into the fitted TMD distributions by means of the Monte Carlo sampling method. This allowed us to obtain a very good description of the entire dataset (χ 2 /N dat = 1.02) without the need of introducing ad hoc normalisations. A more detailed analysis of the fit quality shows that both low-and high-energy datasets are separately well described. This is a remarkable achievement given the very high precision of the LHC datasets, especially those from ATLAS.
A particularly interesting aspect of our analysis concerns the QCD convergence of the perturbative series. We performed fits at NLL , NNLL, NNLL , and N 3 LL accuracy and showed that the fit quality improves significantly going from NLL to N 3 LL. The difference between the highest orders, i.e. NNLL and N 3 LL, is moderate but still significant. This shows at the same time that the perturbative series is converging, but also that N 3 LL corrections are relevant in relation to the current experimental uncertainties.
We parameterised the non-perturbative contributions by adopting a reasonably flexible functional form: all nine free parameters turned out to be well constrained, with moderate correlations amongst them. An important feature of our parameterisation of the nonperturbative contribution f NP is its explicit x dependence. We proved that the x-dependent part of f NP is mostly constrained by the rapidity-dependent on-peak data at 7 and 8 TeV from ATLAS. While on the one hand, this was to be expected because the x dependence is strictly connected with the rapidity y, on the other hand it also demonstrates that most of the datasets are not sensitive to the x dependence of TMDs.
Finally, we studied the validity range of TMD factorisation in Drell-Yan by varying the cut on q T /Q. In line with the literature, we found that the region q T 0.2 Q is appropriate when working within the TMD factorisation framework.

JHEP07(2020)117
In this paper we set the foundation for a number of future studies. In the first place, we plan to extend the fitted dataset by including the abundant and precise semi-inclusive DIS data from HERMES [61] and COMPASS [62,100], as well as future data from Jefferson Lab at 12 GeV [101]. On top of providing access to TMD fragmentation functions, we expect that the inclusion of semi-inclusive DIS data will have an impact on the determination of the x dependence of TMD PDFs and will make it possible to determine the flavour dependence of the non-perturbative function f NP . We remark that a better knowledge of TMDs will be important not only to obtain a deeper knowledge of hadron stucture and QCD, but also for precision studies in high-energy processes involving hadrons, for instance for the determination of critical Standard Model parameters such as the W mass [25,63].
In the future, the Electron-Ion Collider will provide an unprecedented opportunity to make progress in the determination of TMDs [102,103]. Nevertheless, we are convinced that the era of precision physics with TMDs has already started and it will be beneficial also for studies at higher energies in the perturbative domain of QCD.

JHEP07(2020)117
for this kind of tasks since many years, the second (ceres-solver) is relatively new and typically used for more complex problems such as image recognition, 3D modeling, etc.. Recently, the xFitter Collaboration [111] has used ceres-solver for fitting collinear PDFs [112], showing that this tool is suitable also for this kind of tasks. Having two independent codes within the same framework turned out to be particularly useful to cross check our results.
All the datasets included in this analysis, except the preliminary STAR data, have been taken from the public HEPData repository [113] in YAML format and slightly adapted to fit our needs.
Finally, we mention that the TMDs sets determined in this analysis will be made publicly available also through the TMDplotter interface [114].
B Integrating over q T Experimental measurements of differential distributions are usually delivered as integrated over finite regions of the final-state kinematic phase space (see eq. (3.2)). As a consequence, in order to compare theoretical predictions to data, it is necessary to carry out these integrations. These nested integrals, if evaluated numerically, represent a heavy task that makes an extraction of TMDs from Drell-Yan data computationally very intensive and thus slow. While the integrals over Q and y do need to be computed numerically, the integration in q T can be carried out analytically which substantially reduces the numerical load. To do so, we exploit the following property of the Bessel functions Neglecting for the moment the dependence on q T of the phase-space reduction factor P (which is strictly correct for inclusive observables in the final-state leptons), the differential cross section in eq. (2.10) has the following structure where S is a function that depends on b T (and on the other kinematic variables) but not on q T . Using eq. (B.2), one finds

(B.4)
In conclusion, the quantity

JHEP07(2020)117
is the indefinite integral over q T (the primitive function) of the cross section in eq. (2.10).
Analogously to the unintegrated cross section, K can be computed numerically by performing a Bessel transform of degree one rather than degree zero. Therefore, the integral over a q T bin can be evaluated by taking the difference of K computed at the bin bounds: which is enormously more convenient than computing the integral numerically.

B.1 Kinematic cuts
In the presence of kinematic cuts, such as those on the final-state leptons, the analytic integration over q T discussed above cannot be directly performed. The reason is that the implementation of these cuts effectively introduces the q T -dependent function P in the integral that prevents the direct use of eq. (B.2). Fortunately, P is a slowly-varying function of q T over the typical bin size. This allows one to approximate the integral over the bins in q T as (B.8) Unfortunately, this structure is inconvenient because it mixes different bin bounds and prevents a recursive computation. However, it is possible to go further and, assuming that the bin width is small enough, we expand P in the following two equivalent ways P q T,max + q T,min 2 = P (q T,min + ∆q T ) P (q T,min ) + P (q T,min ) ∆q T P (q T,max − ∆q T ) P (q T,max ) − P (q T,max ) ∆q T , (B.9) with ∆q T = q T,max − q T,min 2 . (B.10) Plugging the expansions above into eq. (B.8), one finds (B.11) The advantage of this formula as compared to eq. (B.8) is that each of the terms in the r.h.s. depends on one single bin bound in q T rather than on a combination of two consecutive bounds. This allows for a recursive computation of predictions in neighbouring bins in q T .

JHEP07(2020)117 C Cuts on the final-state leptons
In this section, we derive the explicit expression of the phase-space reduction factor P introduced in section 2. This factor is defined as 11 where p 1 and p 2 are the four-momenta of the outgoing leptons. The integral in the numerator extends over the fiducial region defined by the cuts on the final-state leptons. The quantity L ⊥ is defined as where L µν is the (parity-conserving part of the) leptonic tensor that, assuming massless leptons, reads while the transverse metric is given by The vectors z µ and t µ , in the Collins-Soper frame, are defined as and they are such that z 2 = −1, t 2 = 1 and (z · q) = 0. The effect of integrating over the fiducial region in the numerator of eq. (C.1) can be implemented by defining a generalised θ-function, Φ(p 1 , p 2 ), that is equal to one inside the fiducial region and zero outside. This allows one to integrate also the numerator over the full phase-space of the two outgoing leptons. Next, we integrate out one of the momenta, say p 2 , exploiting the momentumconservation δ-function: where we have renamed p = p 1 . The remaining δ-functions can be used to constrain two of the four components of the momentum p. The first, δ(p 2 ), is typically used to set the energy component of p, p 0 , on the mass shell. Since the leptons are massless, this produces (C.7) 11 In eq. (C.1) a parity-violating term is neglected. We will argue in section C.1 that its contribution is negligible for realistic cuts.

JHEP07(2020)117
Of course, the four-momentum p appearing in the rest of the integrand has to be set on shell (p 0 = |p|). Now we express the three-dimensional measure d 3 p in terms of the transverse momentum p T , the pseudo-rapidity η, and the azimuthal angle φ of the lepton: Now we consider the second δ-function, δ((q − p) 2 ), in eq. (C.6). It is convenient to express the vectors q and p in terms of the respective invariant mass, pseudo-rapidity, and transverse momentum: with M = Q 2 + |q T | 2 . Without loss of generality, we assume that the two-dimensional vector q T is aligned with the x axis so that p T · q T = |p T ||q T | cos φ. 12 This leads to , (C.11) where the vector p is understood to be on-shell. Now we compute L ⊥ (p, q − p) contracting L µν in eq. (C.3) with the transverse metric g µν ⊥ in eq. (C.4) using eq. (C.9): We can now integrate out one of the variables in the integrals in eq. (C.11) by making use of the remaining δ-function. Somewhat counterintuitively, it is convenient to integrate over |p T |. This produces where p T is defined as |q T | − cos φ , (C.14) 12 In the general case in which qT forms an angle β with the x axis, the scalar product would result in |pT ||qT | cos(φ − β). However, for observables inclusive in azimuthal angle, the angle β can always be reabsorbed in a redefinition of φ.
(C. 20) In order to compute the numerator of eq. (C.13), we need to insert the appropriate function Φ. Typically, in DY production the kinematic cuts are imposed independently on the same variables for both the final-state leptons. Therefore, the function Φ factorises into two identical functions acting on each lepton momentum: Φ(p 1 , p 2 ) = Θ(p 1 )Θ(p 2 ) . (C.21) We are specifically interested in kinematic cuts on the rapidity and on the transverse momentum of the following kind η min < η 1(2) < η max and |p T,1(2) | > p T,min . (C.28) As an example, figure 9 shows the integration domain of the numerator of eq. (C.13) for p T,min = 20 GeV and −η min = η max = 2.4 at Q = 91 GeV, |q T | = 10 GeV, and y = 1. The grey band corresponds to the region −1 ≤ cos φ ≤ 1. The θ-functions in the first line of eq. (C.26) limits the region to the vertical strip defined by η min < η < η max (black vertical lines), the θ-function in the second line defines the region above the red line, finally the θfunctions in the third line defines the region below the blue and green lines. The intersection of all regions gives the red-shaded area corresponding to the integration domain.

C.1 Azimuthally-dependent contributions
Azimuthally-dependent modulations disappear in the cross sections if the integration over the azimuthal angle of the virtual boson, Φ, is complete. In the presence of cuts on the final-state leptons, these modulations could generate contributions that were neglected in our analysis, but could be relevant for the description of high-precision data. We first consider parity-violating effects that generate a sin Φ modulation [26]. These contributions stem from interference of the antisymmetric contributions to the lepton tensor, proportional to p µ 1 p ν 2 µνρσ , and to the hadronic tensor, proportional to µν ⊥ defined as µν ⊥ ≡ µνρσ t ρ z σ , (C. 36) where t µ and z µ are given in eq. (C.5). Therefore, the contributions we are after result from the contraction of the following Lorentz structures Therefore, for observables inclusive in the lepton phase space, the parity-violating term does not give any contribution. Conversely, the presence of cuts on the final-state leptons may prevent eq. (C.38) from being satisfied, leaving a residual contribution. In order to JHEP07(2020)117 quantify this effect, we have taken the same steps performed above to integrate L PV over the fiducial region. It turns out that, for realistic cuts, the numerical size of P PV relative to the parity-conserving P is never larger than O(10 −6 ). We conclude that the impact of parity-violating effects in the present analysis is negligible. Finally, we consider also cos Φ modulations, stemming from the following contraction: where the (symmetric part of the) leptonic tensor reads: L µν = 4(p µ 1 p ν 2 + p µ 2 p ν 1 − g µν p 1 p 2 ) . (C. 40) We find that Due to the presence of the overall factor sinh(y − η), for relatively central rapidities and for symmetric cuts this term is expected to be very small, in particular to be comparable in size to the parity violating contribution. Moreover, this term would be multiplied by a structure function that has been measured to be small, below 4% in the region of interest here [115].
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.