删除或更新信息,请邮件至freekaoyan#163.com(#换成@)

Scalar-pseudoscalar pair production at the Large Hadron Collider at NLO+NLL accuracy in QCD

本站小编 Free考研考试/2022-01-01

闂傚倸鍊搁崐鎼佸磹閻戣姤鍤勯柛鎾茬閸ㄦ繃銇勯弽顐粶缂佲偓婢跺绻嗛柕鍫濇噺閸e湱绱掗悩闈涒枅闁哄瞼鍠栭獮鎴﹀箛闂堟稒顔勯梻浣告啞娣囨椽锝炴径鎰﹂柛鏇ㄥ灠濡﹢鏌涢…鎴濇灀闁圭ǹ鍟村娲川婵犲孩鐣烽悗鍏夊亾闁归棿绀佺粻鏍ㄤ繆閵堝懏鍣洪柡鍛叀楠炴牜鈧稒岣跨粻姗€鏌i埡浣规崳缂佽鲸鎸婚幏鍛槹鎼淬倗鐛ラ梻渚€娼荤紞鍥╃礊娴e壊鍤曞┑鐘崇閸嬪嫰鏌i幘铏崳妞は佸洦鈷戦柛蹇氬亹閵堟挳鏌¢崨顔剧疄闁诡噯绻濆畷鎺楁倷缁瀚肩紓鍌欑椤戝牆鈻旈弴銏″€块柛褎顨嗛悡娆撴煕閹存瑥鈧牜鈧熬鎷�2濠电姷鏁告慨鐑藉极閹间礁纾婚柣鎰惈閸ㄥ倿鏌涢锝嗙缂佺姵澹嗙槐鎺斺偓锝庡亾缁扁晜绻涘顔荤盎閹喖姊洪崘鍙夋儓妞ゆ垵娲ㄧ划娆掔疀濞戞瑢鎷洪梺闈╁瘜閸樺ジ宕濈€n偁浜滈柕濞垮劜椤ャ垻鈧娲滈弫濠氬春閳ь剚銇勯幒鎴濐仾闁抽攱鍨块弻娑樷槈濮楀牆浼愭繝娈垮櫙缁犳垿婀佸┑鐘诧工閹冲孩绂掓潏鈹惧亾鐟欏嫭绀冩俊鐐扮矙瀵偊骞樼紒妯轰汗閻庤娲栧ú銈夌嵁濡ゅ懏鈷掑〒姘e亾婵炰匠鍛床闁割偁鍎辩壕褰掓煛瀹擃喒鍋撴俊鎻掔墢閹叉悂寮崼婵婃憰闂佹寧绻傞ˇ顖炴倿濞差亝鐓曢柟鏉垮悁缁ㄥジ鏌涢敐鍕祮婵﹨娅i幏鐘诲灳閾忣偅顔勯梻浣规偠閸庢粓宕惰閺嗩亪姊婚崒娆戝妽閻庣瑳鍛床闁稿本顕㈠ú顏勵潊闁靛牆鎳愰敍娑㈡⒑閸︻厼鍔嬫い銊ユ閸╂盯骞嬮敂鐣屽幈濠电娀娼уΛ妤咁敂閳哄懏鐓冪憸婊堝礈濞嗘垹绀婂┑鐘叉搐缁犳牠姊洪崹顕呭剱缂傚秴娲弻宥夊传閸曨偂绨藉┑鐐跺亹閸犳牕顫忛搹瑙勫磯闁靛ǹ鍎查悵銏ゆ⒑閻熸澘娈╅柟鍑ゆ嫹
濠电姷鏁告慨鐑藉极閸涘﹥鍙忓ù鍏兼綑閸ㄥ倻鎲搁悧鍫濈瑲闁稿顦甸弻鏇$疀鐎n亖鍋撻弴銏″€峰┑鐘插閸犳劗鈧箍鍎卞Λ娆撳矗韫囨稒鐓忛柛顐g箥濡插綊鏌嶉柨瀣伌闁哄本绋戦埥澶婎潨閸繀绱g紓鍌欑劍椤ㄥ棛鏁Δ浣衡攳濠电姴娴傞弫鍐煥濠靛棙澶勯柛鎺撶☉椤啴濡堕崘銊т痪濠碘槅鍋勯崯顖炲箞閵娾晛鐒垫い鎺戝閻撳繘鏌涢锝囩畺闁挎稑绉垫穱濠囶敃閵忕媭浼冮梺鍝勭焿缁查箖骞嗛弮鍫晬婵犲﹤鎲涢敐澶嬧拺闁告縿鍎辨牎闂佺粯顨堟慨鎾偩閻戣棄顫呴柕鍫濇噽椤旀帒顪冮妶鍡樷拻闁哄拋鍋婂畷銏ゅ箹娴e厜鎷洪梺鍛婄☉閿曘儳绮堢€n偆绠惧ù锝呭暱濞诧箓宕愰崼鏇熺叆婵犻潧妫欓ˉ鎾趁瑰⿰鍕煉闁哄瞼鍠撻埀顒佺⊕宀h法绮婚弽褜鐔嗛悹鍝勬惈椤忣偆绱掓潏銊ョ闁逞屽墾缂嶅棙绂嶇捄浣曠喖鍩€椤掑嫭鈷戠紒顖涙礃閺夊綊鏌涚€n偅灏い顏勫暣婵″爼宕卞Δ鈧ḿ鎴︽⒑缁嬫鍎愰柟鐟版喘瀵顓奸崶銊ョ彴闂佸搫琚崕鍗烆嚕閺夊簱鏀介柨鐔哄Х閻e搫霉濠婂啰鍩g€殿喛顕ч濂稿醇椤愶綆鈧洭姊绘担鍛婂暈闁圭ǹ顭烽幃鐑藉煛娴g儤娈惧銈嗙墬缁嬫垿顢氶柆宥嗗€垫繛鎴烆仾椤忓懐顩叉い鏍ㄥ焹閺€浠嬫煟閹邦剙绾ч柍缁樻礀闇夋繝濠傚缁犵偟鈧鍠楅悡锟犮€佸Δ鍛妞ゆ巻鍋撻柍褜鍓欓悥濂稿蓟閿濆绠涙い鏃囧Г濮e嫰姊虹涵鍛棄闁稿﹤娼″璇测槈閵忕姈褔鏌涢妷顔句虎闁靛繈鍊栭ˉ鍡楊熆鐠轰警鍎戠紒鈾€鍋撳┑鐘垫暩婵挳宕愰幖浣告辈闁挎繂妫庢禍婊堝箹濞n剙鐒烘繛鍫熸礋閺屾洟宕惰椤忣參鏌涢埡鍐ㄤ槐妞ゃ垺锕㈤幃娆忣啅椤旇崵妫繝鐢靛У椤旀牠宕归柆宥呯闁规儼妫勯拑鐔兼煥閻斿搫孝闁绘劕锕弻宥嗘姜閹殿喖濡介梺璇茬箣缁舵艾顫忓ú顏勫窛濠电姴瀚崰娑㈡⒑缁嬫鍎愰柟鐟版搐椤繒绱掑Ο璇差€撻梺鍛婄缚閸庤櫕顨欏┑鐘垫暩閸嬫﹢宕犻悩璇插耿闁归偊浜濋惈蹇涙⒒娴h櫣甯涢柛鏃€顨婂顐﹀传閵壯傜瑝闂佸搫鍟悧濠囨偂濞嗘挻鐓欐い鏍ф閼活垰鈻撻崼鏇熲拺鐎规洖娲ㄧ敮娑欐叏婵犲倻绉烘鐐茬墦婵℃悂濡锋惔锝呮灁闁归濞€楠炴捇骞掑┑鍥ㄧグ闂傚倸鍊烽悞锕傚箖閸洖纾圭憸蹇曞垝婵犳艾绠婚悹鍥蔼閹芥洟姊虹紒妯活梿婵炲拑缍侀幆灞解枎閹惧鍘电紓浣割儏閻忔繈顢楅姀銈嗙厵妞ゆ梻鏅幊鍥ㄦ叏婵犲嫬鍔嬮悗鐢靛帶閳诲酣骞嬮悩妯荤矌缁辨挻鎷呴崫鍕戯綁鏌eΔ浣圭妞ゃ垺宀搁弫鎰緞濡粯娅囬梻浣稿暱閻忓牓寮插⿰鍫熷€靛┑鐘崇閳锋垹鎲搁悧鍫濈瑨濞存粈鍗抽弻娑樜熼崫鍕ㄦ寖缂備緡鍠楅悷鈺佺暦閻旂⒈鏁嶆繛鎴炲笚鐎氬ジ姊绘担鍛婅础閺嬵亝绻涚€电ǹ鍘撮柛鈹垮劜瀵板嫰骞囬鐘插箰闂備礁澹婇崑鎺楀磻閸曨剚娅犻悗鐢电《閸嬫挾鎲撮崟顒傤槬缂傚倸绉撮敃銉︾┍婵犲偆娼扮€光偓婵犲唭顏勨攽閻樻剚鍟忛柛銊ゅ嵆婵″爼骞栨担姝屾憰濠电偞鍨惰彜婵℃彃鐗婇幈銊ノ旈埀顒勬偋婵犲洤鏋侀柛鎾楀懐锛濇繛杈剧到閹碱偅鐗庨梺姹囧焺閸ㄦ娊宕戦妶澶婃槬闁逞屽墯閵囧嫰骞掗崱妞惧闂備浇顕х换鎴︽嚌妤e啠鈧箓宕归鍛缓闂侀€炲苯澧存鐐插暢椤﹀湱鈧娲栧畷顒勬箒闂佸搫顦扮€笛囧窗濡皷鍋撶憴鍕閺嬵亪鎽堕弽顬″綊鏁愭径瀣彸闂佹眹鍎烘禍顏勵潖缂佹ɑ濯村〒姘煎灡閺侇垶姊虹憴鍕仧濞存粠浜滈~蹇旂鐎n亞顦板銈嗙墬缁嬫帒鈻嶉弽顓熲拺闁告繂瀚埢澶愭煕濡湱鐭欓柟顔欍倗鐤€婵炴垶鐟ч崢閬嶆⒑閺傘儲娅呴柛鐕佸灣缁牓鍩€椤掆偓椤啴濡惰箛鏇炵煗闂佸搫妫欑粩绯村┑鐘垫暩婵兘寮崨濠冨弿濞村吋娼欓崹鍌炴煕閿旇骞樼紒鈧繝鍌楁斀闁绘ê寮堕幖鎰版煟閹烘垹浠涢柕鍥у楠炴帒顓奸崼婵嗗腐闂備焦鍓氶崹鍗灻洪悢鐓庤摕闁哄洢鍨归獮銏′繆閵堝倸浜鹃梺鍝勬4缂嶄線寮婚敍鍕勃闁告挆鍕灡婵°倗濮烽崑鐐垫暜閿熺姷宓侀悗锝庡枟閸婂鏌涢埄鍐夸緵婵☆値鍐f斀闁挎稑瀚禍濂告煕婵犲啰澧遍柡渚囧櫍閹瑩宕崟顓犲炊闂備礁缍婇崑濠囧窗濮樿埖鍎楁繛鍡楃箚閺€浠嬫煟濡搫绾у璺哄閺屾稓鈧綆鍋勬慨宥夋煛瀹€瀣М濠殿喒鍋撻梺闈涚箚閸撴繂袙閸曨垱鐓涘ù锝呮憸婢э附鎱ㄦ繝鍕笡闁瑰嘲鎳愮划娆撳箰鎼粹檧鍋撻姘f斀闁绘﹩鍠栭悘顏堟煥閺囨ê鐏╅柣锝囧厴椤㈡稑鈽夊鍡楁闂佽瀛╃粙鎺楁晪婵炲瓨绮犻崹璺侯潖濞差亜宸濆┑鐘插閻e灚绻濆▓鍨仴濡炲瓨鎮傞獮鍡涘籍閸繍娼婇梺鎸庣☉鐎氼喛鍊存繝纰夌磿閸嬫垿宕愰弽顓炵婵°倕鎳庣粣妤呭箹濞n剙鐏い鈺傚絻铻栭柨婵嗘噹閺嗙偤鏌i幘瀵告创闁哄本鐩俊鐑芥晲閸涱収鐎撮梻浣圭湽閸斿秹宕归崸妤€钃熼柨婵嗩槹閸嬪嫰鏌涘▎蹇fЧ闁绘繃妫冨铏光偓鍦У椤ュ銇勯敂鐐毈闁绘侗鍠栬灒闁兼祴鏅濋ˇ鈺呮⒑缂佹◤顏勭暦椤掑嫷鏁嗛柕蹇娾偓鑼畾闂佺粯鍔︽禍婊堝焵椤掍胶澧悡銈嗙節闂堟稒顥戦柡瀣Ч閺岋繝宕堕埡浣锋喚缂傚倸鍊瑰畝鎼佹偂椤愶箑鐐婇柕濞р偓濡插牓鎮楅悷鐗堝暈缂佽鍟存俊鐢稿礋椤栨氨顔掑┑掳鍊愰崑鎾绘煕閻曚礁鐏︽い銏$懇閺佹捇鏁撻敓锟�20婵犵數濮撮惀澶愬级鎼存挸浜炬俊銈勭劍閸欏繘鏌i幋锝嗩棄缁炬儳顭烽弻锝呂熷▎鎯ф缂備胶濮撮悘姘跺Φ閸曨喚鐤€闁圭偓鎯屽Λ鈥愁渻閵堝骸浜濇繛鍙夅缚閹广垹鈹戠€n偒妫冨┑鐐村灥瀹曨剟宕滈幍顔剧=濞达絽鎼牎闂佹悶鍔屽ḿ鈥愁嚕婵犳艾围闁糕剝锚瀵潡姊鸿ぐ鎺戜喊闁稿繑锕㈠畷鎴﹀箻濠㈠嫭妫冮崺鈧い鎺戝閻撴繈鏌¢崘銊у妞ゎ偄鎳橀弻锝呂熼悜姗嗘¥闂佺娅曢幑鍥Χ椤忓懎顕遍柡澶嬪灩椤︺劑姊洪崘鍙夋儓闁挎洏鍎甸弫宥夊川椤栨粎锛濋梺绋挎湰閻熝囁囬敂濮愪簻闁挎棁顕ч悘锔姐亜閵忊€冲摵妞ゃ垺锕㈡慨鈧柣姗€娼ф慨锔戒繆閻愵亜鈧牕顔忔繝姘;闁规儳顕弧鈧梺閫炲苯澧撮柡灞芥椤撳ジ宕ㄩ銈囧耿闂傚倷鑳剁划顖氼潖婵犳艾鍌ㄧ憸鏂款嚕閸涘﹦鐟归柍褜鍓熷濠氬即閵忕娀鍞跺┑鐘茬仛閸旀牗鏅ラ梻鍌欒兌鏋Δ鐘叉憸缁棁銇愰幒鎴f憰濠电偞鍨崹褰掑础閹惰姤鐓忓┑鐐茬仢閸旀碍銇勯鐔告珚婵﹦鍎ょ€电厧鈻庨幋鐘虫缂傚倸鍊哥粔鎾晝椤忓牏宓侀柛鎰╁壆閺冨牆绀冮柍杞扮劍閻庮參姊绘担鍛婂暈婵炶绠撳畷锝嗘償閵娿儲杈堥梺璺ㄥ枔婵敻鍩涢幋锔界厱婵犻潧妫楅顏呫亜閵夛箑鐏撮柡灞剧〒閳ь剨缍嗛崑鍛暦鐏炵偓鍙忓┑鐘插暞閵囨繄鈧娲﹂崑濠傜暦閻旂厧鍨傛い鎰癁閸ャ劉鎷洪梺鍛婄☉閿曘儵鍩涢幇鐗堢厽婵°倕鍟埢鍫燁殽閻愭彃鏆i柡浣规崌閹晠鎼归锝囧建闂傚倷绀侀幉鈥趁洪敃鍌氱婵炲棙鎸婚崑鐔访归悡搴f憼闁抽攱鍨垮濠氬醇閻旀亽鈧帞绱掗悩鍐插摵闁哄本鐩獮妯尖偓闈涙憸閻ゅ嫰姊虹拠鈥虫灀闁逞屽墯閺嬪ジ寮告惔銊︾厵闂侇叏绠戦弸銈嗐亜閺冣偓濞叉ḿ鎹㈠┑瀣潊闁挎繂妫涢妴鎰渻閵堝棗鐏ユ俊顐g〒閸掓帡宕奸妷銉у姦濡炪倖甯掔€氼參宕愰崹顐ょ闁割偅绻勬禒銏$箾閸涱厾效闁哄矉绻濋崺鈧い鎺戝绾偓闂佺粯鍨靛Λ妤€鈻撻锔解拺闁告稑锕ユ径鍕煕鐎n偄娴€规洏鍎抽埀顒婄秵閸犳鎮¢弴銏$厸闁搞儯鍎辨俊鍏碱殽閻愮摲鎴炵┍婵犲洤鐭楀璺猴功娴煎苯鈹戦纭锋敾婵$偠妫勯悾鐑筋敃閿曗偓缁€瀣亜閹邦喖鏋庡ù婊勫劤闇夐柣妯烘▕閸庢粎绱撳鍡欏ⅹ妞ゎ叀娉曢幑鍕倻濡粯瀚抽梻浣呵圭换鎴犲垝閹捐钃熸繛鎴欏焺閺佸啴鏌ㄥ┑鍡橆棤妞わ负鍔戝娲传閸曨剙顎涢梺鍛婃尵閸犳牠鐛崘顭戞建闁逞屽墴楠炲啫鈻庨幘鎼濠电偞鍨堕〃鍛此夊杈╃=闁稿本鐟ㄩ崗灞解攽椤旂偓鏆╅柡渚囧櫍閸ㄩ箖骞囨担鍦▉濠电姷鏁告慨鐢告嚌妤e啯鍊峰┑鐘叉处閻撱儲绻濋棃娑欘棡闁革絿枪椤法鎲撮崟顒傤槹濠殿喖锕ュ浠嬪箠閿熺姴围闁告侗鍠氶埀顒佸劤閳规垿鎮欓幓鎺旈獓闂佹悶鍔屽ḿ锟犵嵁婵犲伣鏃堝礃閳轰胶锛忛梺鑽ゅ仦缁嬪牓宕滃┑瀣€跺〒姘e亾婵﹨娅e☉鐢稿川椤斿吋閿梻鍌氬€哥€氼剛鈧碍婢橀悾鐑藉即閵忕姷顓洪梺鎸庢濡嫰鍩€椤掑倹鏆柡灞诲妼閳规垿宕卞☉鎵佸亾濡や緡娈介柣鎰缂傛氨绱掓潏銊ユ诞闁诡喒鏅涢悾鐑藉炊瑜夐幏浼存⒒娴e憡鎯堝璺烘喘瀹曟粌鈹戦崱鈺佹闂佸憡娲﹂崑鈧俊鎻掔墛缁绘盯宕卞Δ浣侯洶婵炲銆嬮幏锟�
He-Yi Li 1,2,
, Ren-You Zhang 1,2,,
, Yu Zhang 3,4,
, Wen-Gan Ma 1,2,
, Ming-Ming Long 1,2,
, Shu-Xiang Li 1,2,
,
Corresponding author: Ren-You Zhang, zhangry@ustc.edu.cn, corresponding author
1.State Key Laboratory of Particle Detection and Electronics, University of Science and Technology of China, Hefei 230026, China
2.Department of Modern Physics, University of Science and Technology of China, Hefei 230026, China
3.Institutes of Physical Science and Information Technology, Anhui University, Hefei 230601, China
4.School of Physics and Materials Science, Anhui University, Hefei 230601, China
Received Date:2021-05-30
Available Online:2021-12-15
Abstract:We thoroughly investigate both transverse momentum and threshold resummation effects on scalar-pseudoscalar pair production via quark-antiquark annihilation at the $ 13 \; \text{TeV}$ Large Hadron Collider at QCD NLO+NLL accuracy. A factorization method is introduced to properly supplement the soft-gluon (threshold) resummation contribution from parton distribution functions to the resummed results obtained by the Collins-Soper-Sterman resummation approach. We find that the impact of the threshold-resummation improved PDFs is comparable to the resummation effect of the partonic matrix element and can even predominate in high invariant mass regions. Moreover, the loop-induced gluon-gluon fusion channel in the type-I two-Higgs-doublet model is considered in our calculations. The numerical results show that the electroweak production via quark-antiquark annihilation dominates over the gluon-initiated QCD production by $ 1 \sim 2$ orders of magnitude.

HTML

--> --> -->
I.INTRODUCTION
Primary tasks at the Large Hadron Collider (LHC) include the precision test of the standard model (SM) and the search for new physics beyond the SM (BSM). Since the discovery of the $ 125\; \text{GeV} $ Higgs boson by both ATLAS and CMS collaborations at the LHC in 2012 [1, 2], the SM has become the most successful theory in describing the interactions of fundamental particles. However, the discovery of this SM-like Higgs boson is merely one step toward fully investigating the electroweak symmetry breaking (EWSB). As is well known, the theoretical predictions of the SM are not always compatible with experimental observations, such as the dark matter in the universe, the oscillation of neutrinos, the huge hierarchy between electroweak and Planck scales, and the fine-tuning problem of Higgs mass. These conceptional and experimental difficulties encountered by the SM imply the existence of new physics beyond the SM.
We may extend the SM by enlarging its gauge symmetry and/or introducing much more gauge multiplets to construct a new physics model. Among all the BSM theories, the two-Higgs-doublet model (2HDM) [3] is one of the simplest extensions of the SM. The Higgs sector responsible for the EWSB consists of two complex scalar isospin doublets, and the minimal supersymmetric standard model is a particular realization of the 2HDM. After EWSB, the three Goldstone modes $ G^{\pm} $ and $ G^0 $ in the Higgs sector of the 2HDM are absorbed by the weak gauge bosons $ W^{\pm} $ and $ Z^0 $, respectively, providing the longitudinal polarizations of $ W^{\pm} $ and $ Z^0 $. The remaining five mass eigenstates of the Higgs sector are the so-called $ {\cal{C}}{\cal{P}} $-even Higgs bosons $ h^0 $ and $ H^0 $, $ {\cal{C}}{\cal{P}} $-odd Higgs boson $ A^0 $, and charged Higgs bosons $ H^{\pm} $.
Clearly, any discovery of a BSM Higgs boson will be an evidence for the existence of a new Higgs sector. At the LHC, the neutral Higgs bosons of the 2HDM can be produced both singly and in identical or mixed pairs. The dominant mechanism for single production of neutral Higgs bosons is gluon-gluon fusion. Concerning scalar-pseudoscalar pairs, the electroweak production via quark-antiquark annihilation, $ pp \rightarrow q\bar{q} \rightarrow Z^{\ast} \rightarrow H^0A^0/h^0A^0 $, can dominate over the QCD production via gluon-gluon fusion [4, 5], even by orders of magnitude. Considerable efforts have been devoted to search for BSM neutral Higgs bosons. Particularly, the exotic decays of heavy scalar (pseudoscalar), such as $ H^0 \rightarrow A^0Z^0 $ ($ A^0 \rightarrow H^0Z^0 $), have attracted attention at the LHC in recent years [6, 7]. The scalar-pseudoscalar pair production is dominated by the Drell-Yan channel; it is an ideal process to investigate the Higgs gauge coupling $ g_{H^0A^0Z^0}/g_{h^0A^0Z^0} $ and should thus be thoroughly invegtigated. The next-to-leading order (NLO) QCD corrections to neutral Higgs-boson pair production at hadron colliders were calculated in Refs. [4, 8], which showed that the QCD corrections can enhance the cross section of $ h^0A^0 $ production by approximately 30%.
Fixed-order perturbative predictions would be unreliable when the exponential enhancement from soft gluon dominates at the edge of the phase space. The large logarithms, such as $ \alpha_s^n ( M^2/p_T^2) \ln^m (M^2/p_T^2) $ at small $ p_T $ and $ \alpha_s^n (1-z)^{-1} \ln^{m}(1-z) $ when $ z = M^2/\hat{s} \rightarrow 1 $, should be resummed in precision calculations. We extract such logarithms in the partonic matrix element by adopting the Collins-Soper-Sterman resummation technique [9-17] and analyze the threshold resummation effect from parton distribution functions (PDFs) by using the factorization method proposed in Ref. [18]. Generally, the resummation corrections are only considered for fixing the unnatural behaviours in the small-$ p_T $ and threshold regions, while the fixed-order predictions are suitable for describing the kinematics far away from the edge of the final-state phase space. Thus, the resummation results should be matched with the fixed-order predictions to obtain a reliable description in all kinematical regions.
In this study, we thoroughly analyze scalar-pseudoscalar pair production at the $ 13\; \text{TeV} $ LHC within the type-I 2HDM at the NLO and next-to-leading logarithmic (NLL) accuracy in QCD. The rest of this paper is organized as follows. In Sec. II, we briefly review the 2HDM. In Sec. III, we present the calculation strategies for $ \phi^0A^0 $ associated production at QCD NLO+NLL accuracy, including the Collins-Soper-Sterman resummation technique and the factorization method for assessing the impact of the threshold-resummation improved PDFs. The numerical results and discussion for both the integrated cross section and the differential distributions with respect to the transverse momentum and invariant mass of the final-state $ \phi^0A^0 $ system are provided in Sec. IV. Finally, a short summary is given in Sec. V.
II.BRIEF REVIEW OF 2HDM
In contrast to the SM, the Higgs sector of the 2HDM consists of two complex scalar $ SU(2)_{\rm L} $ doublets $ \Phi_{1,2} $ with hypercharge $ Y = +1 $. The most general scalar potential, which is invariant under the $ SU(2)_{\rm L} \otimes U(1)_{\rm Y} $ electroweak gauge symmetry and a discrete $ Z_2 $ symmetry $ \Phi_i \rightarrow (-1)^{i+1} \Phi_i\; (i = 1, 2) $, is given by
$ \begin{aligned}[b] V(\Phi_{1}, \Phi_{2}) =& m^2_{11} \Phi^{\dagger}_{1} \Phi_{1} + m^2_{22} \Phi^{\dagger}_{2} \Phi_{2} - ( m^2_{12} \Phi^{\dagger}_{1}\Phi_{2} + \text{h.c.} ) \\&+ \frac{1}{2} \lambda_1 ( \Phi^{\dagger}_{1} \Phi_{1} )^{2} + \frac{1}{2} \lambda_2 ( \Phi^{\dagger}_{2} \Phi_{2} )^{2} \\& + \lambda_{3} ( \Phi^{\dagger}_{1} \Phi_{1} ) ( \Phi^{\dagger}_{2} \Phi_{2} ) + \lambda_{4} ( \Phi^{\dagger}_{1} \Phi_{2} ) ( \Phi^{\dagger}_{2} \Phi_{1} ) \\&+ \dfrac{1}{2} \left[ \lambda_{5} (\Phi^{\dagger}_{1}\Phi_{2})^{2} + \text{h.c.} \right], \end{aligned} $
(1)
where the dimension-two term $ m_{12}^2 $ is tolerated because it only breaks the $ Z_2 $ symmetry softly, and $ m_{11, 22}^2 $, $ \lambda_{1, 2, 3, 4} $ are forced to be real owing to the hermiticity of the scalar potential. The two Higgs doublets can be parameterized as [3]
$ {\Phi _i} = \left( {\begin{array}{*{20}{c}}{\phi _i^ + }\\{({v_i} + {\rho _i} + i{\eta _i})/\sqrt 2 }\end{array}} \right)\;\;\;\;\;\;\;\;\;\;\;(i = 1,2), $
(2)
where $ v_1 $ and $ v_2 $ are the vacuum expectation values (VEVs) of the neutral components of $ \Phi_{1} $ and $ \Phi_{2} $, respectively. In a $ {\cal{C}}{\cal{P}} $-conserving 2HDM, both $ m_{12}^2 $ and $ \lambda_5 $ are real, and so are $ v_1 $ and $ v_2 $. The eight mass eigenstates of the Higgs sector are given by
$ \left( {\begin{array}{*{20}{c}}{{H^0}}\\{{h^0}}\end{array}} \right) = \left( {\begin{array}{*{20}{c}}{\cos \alpha }&{\sin \alpha }\\{ - \sin \alpha }&{\cos \alpha }\end{array}} \right)\left( {\begin{array}{*{20}{c}}{{\rho _1}}\\{{\rho _2}}\end{array}} \right), $
(3)
$ \left( {\begin{array}{*{20}{c}}{{G^0}}\\{{A^0}}\end{array}} \right) = \left( {\begin{array}{*{20}{c}}{\cos \beta }&{\sin \beta }\\{ - \sin \beta }&{\cos \beta }\end{array}} \right)\left( {\begin{array}{*{20}{c}}{{\eta _1}}\\{{\eta _2}}\end{array}} \right), $
(4)
$ \left( {\begin{array}{*{20}{c}}{{G^ \pm }}\\{{H^ \pm }}\end{array}} \right) = \left( {\begin{array}{*{20}{c}}{\cos \beta }&{\sin \beta }\\{ - \sin \beta }&{\cos \beta }\end{array}} \right)\left( {\begin{array}{*{20}{c}}{\phi _1^ \pm }\\{\phi _2^ \pm }\end{array}} \right), $
(5)
where α is the mixing angle in the $ {\cal{C}}{\cal{P}} $-even Higgs sector and $ \beta = \arctan \dfrac{v_2}{v_1} $ describes the mixing in the $ {\cal{C}}{\cal{P}} $-odd and charged Higgs sectors. After the spontaneous electroweak symmetry breaking, three out of eight degrees of freedom from $ \Phi_{1,2} $ that correspond to Nambu-Goldstone bosons $ G^{\pm} $and $ G^0 $ are respectively absorbed by weak gauge bosons $ W^{\pm} $ and $ Z^0 $, providing the longitudinal polarizations of $ W^{\pm} $ and $ Z^0 $. The remaining five degrees of freedom become the aforementioned five physical Higgs bosons: two $ {\cal{C}}{\cal{P}} $-even Higgs bosons $ h^0 $ and $ H^0 $, one $ {\cal{C}}{\cal{P}} $-odd Higgs boson $ A^0 $, and a pair of charged Higgs bosons $ H^{\pm} $. In this study, the seven input parameters for the Higgs sector of a $ {\cal{C}}{\cal{P}} $-conserving 2HDM are chosen as
$ \Big\{ m_{h^0},\, m_{H^0},\, m_{A^0},\, m_{H^{\pm}},\, m^2_{12},\, \sin(\beta - \alpha),\, \tan\beta \Big\}, $
(6)
which are implemented as the “physical basis” in 2HDMC [19]. Then, the Higgs potential in Eq. (1) can be completely determined by the seven Higgs parameters in Eq. (6) and v, where $v \equiv \sqrt{v_1^2+v_2^2} = (\sqrt{2}G_{\rm F})^{-1/2} \approx 246 \; \text{GeV}$ has been classified as an electroweak input parameter.
To guarantee the absence of Higgs-mediated flavor changing neutral currents at the tree level, the $ Z_2 $ symmetry should be extended to the Yukawa sector. Given that the two Higgs doublets $ \Phi_{1,2} $ have opposite $ Z_2 $ charges, each flavor of quark/lepton can only couple to one of the two Higgs doublets. There are four allowed types of Yukawa interaction corresponding to the four independent $ Z_2 $ charge assignments on the quark and lepton $ SU(2)_{\rm L} $ multiplets (Table 1). The Yukawa Lagrangian of the 2HDM can be expressed in terms of Higgs mass eigenstates as
2HDM $ \Phi_1 $ $ \Phi_2 $ $ Q_{\rm L} $ $ L_{\rm L} $ $ u_{\rm R} $ $ d_{\rm R} $ $ \ell_{\rm R} $
Type I + ? + + ? ? ?
Type II + ? + + ? + +
Lepton-specific + ? + + ? ? +
Flipped + ? + + ? + ?


Table1.Four types of 2HDMs and the corresponding $ Z_2 $ charge assignments on Higgs, quark, and lepton $ SU(2)_{\rm L} $ multiplets.

$ \begin{aligned}[b] {\cal{L}}_{\text{Yukawa}}^{\text{2HDM}} =& \, - \sum\limits_{f = u, d, \ell} \frac{m_f}{v} \left( \xi_h^f \bar{f} f h^0 + \xi_H^f \bar{f} f H^0 - {\rm i} \xi_A^f \bar{f} \gamma^5 f A^0 \right) \\& - \Bigg[ \frac{\sqrt{2} V_{ud}}{v} \bar{u} \left( m_u \xi_A^u P_{\rm L} + m_d \xi_A^d P_{\rm R} \right) d H^+ \\&+ \frac{\sqrt{2} m_{\ell}}{v} \xi_A^{\ell} \bar{\nu} P_{\rm R} \ell H^+ + \text{h.c.} \Bigg], \end{aligned} $
(7)
where $ \xi_{h, H, A}^{f}\; (f = u, d, \ell) $ are the Higgs Yukawa couplings normalized to the SM vertices, and the corresponding values in the type-I, type-II, lepton-specific, and flipped 2HDMs are listed in Table 2.
2HDM Type I Type II Lepton-specific Flipped
$ \xi_h^u $ $ \quad\cos\alpha/\sin\beta\quad $ $ \surd $ $ \surd $ $ \surd $
$ \xi_H^u $ $ \sin\alpha/\sin\beta $
$ \xi_A^u $ $ \cot\beta $
$ \xi_h^d $ $ \cos\alpha/\sin\beta $ $ (\alpha, \beta) \rightarrow (\tilde{\alpha}, \tilde{\beta}) $ $ (\alpha, \beta) \rightarrow (\tilde{\alpha}, \tilde{\beta}) $
$ \xi_H^d $ $ \sin\alpha/\sin\beta $
$ \xi_A^d $ $ -\cot\beta\quad $
$ \xi_h^{\ell} $ $ \cos\alpha/\sin\beta $ $ (\alpha, \beta) \rightarrow (\tilde{\alpha}, \tilde{\beta}) $ $ \surd $
$ \xi_H^{\ell} $ $ \sin\alpha/\sin\beta $
$ \xi_A^{\ell} $ $ -\cot\beta\quad $


Table2.Normalized Higgs Yukawa couplings $ \xi_{h, H, A}^{f} $ $(f = $$ u, d, \ell)$ in the type-I, type-II, lepton-specific, and flipped 2HDMs. $ (\tilde{\alpha}, \tilde{\beta}) = (\alpha, \beta) + \dfrac{\pi}{2} $.

The Higgs gauge interaction is independent of the types of the 2HDM. The couplings of $ h^0 $ and $ H^0 $ to weak gauge boson pair are proportional to $ \sin(\beta - \alpha) $ and $ \cos (\beta - \alpha) $, respectively. The 2HDM parameter space is stringently constrained by the requirement that one out of the two neutral $ {\cal{C}}{\cal{P}} $-even Higgs bosons has physical properties consistent with the $ 125\; \text{GeV} $ scalar discovered at the CERN LHC. It is well known that if one of the two neutral $ {\cal{C}}{\cal{P}} $-even Higgs mass eigenstates is approximately aligned in the two-dimensional Higgs field space with the direction of the Higgs VEV vector $ \vec{v} \equiv (v_1, v_2) $ (the so-called alignment limit), the couplings of this Higgs boson are SM-like. The two alignment limits of the 2HDM are listed in Table 3. Given that the SM-like Higgs boson with mass around $ 125\; \text{GeV} $ seems to be favored by LHC data, we will investigate the scalar-pseudoscalar pair production at the LHC only at the alignment limit.
Alignment limit $ \beta - \alpha $ $ h^0 $ $ H^0 $
I $ \pi/2 $ SM-like
II $ 0 $ SM-like


Table3.Two alignment limits of the 2HDM.

III.CALCULATION STRATEGY
We adopt the 't Hooft-Feynman gauge and take the five-flavor scheme in our calculations. Apart from the top quark, all other light quarks, including the bottom quark, are treated as massless particles. The UV and IR divergences in the QCD loop and real jet emission corrections are regularized by adopting the dimensional regularization scheme [20]. We employ both the Catani-Seymour dipole subtraction method [21, 22] and the two cutoff phase space slicing method [23] to separate the soft and collinear IR singularities of the real emission correction, and then cross check their correctness.
2
A.Electroweak production via quark-antiquark annihilation
-->

A.Electroweak production via quark-antiquark annihilation

The scalar-pseudoscalar pair can be produced at the LHC via Drell-Yan production mechanism. Some representative LO and QCD NLO Feynman diagrams for $ q \bar{q} \rightarrow H^0(h^0)A^0 $ are shown in Fig. 1 . In this study, we categorize $ q g \rightarrow H^0(h^0)A^0 + q $ as the real light-quark emission correction to the quark-antiquark-initiated Drell-Yan channel.
Figure1. Representative LO and QCD NLO Feynman diagrams for $ q\bar{q} \rightarrow H^0(h^0)A^0 $.

The $ h^0A^0Z^0 $ and $ H^0A^0Z^0 $ gauge interactions in the 2HDM are given by
$\begin{aligned}[b] g_{h^0A^0Z^0} = &\frac{e \cos(\beta - \alpha)}{\sin2\theta_{\text{W}}} \left( p_{h^0}-p_{A^0} \right)_\mu, \\ g_{H^0A^0Z^0} =& -\frac{e \sin(\beta - \alpha)}{\sin2\theta_{\text{W}}} \left( p_{H^0}-p_{A^0} \right)_\mu, \end{aligned}$
(8)
where $ \theta_{\text{W}} $ is the Weinberg weak mixing angle, $ p_{h^0, H^0, A^0} $ are the incoming momenta of the corresponding (pseudo)scalars, and μ is the Lorentz index of the vector boson $ Z^0 $. At the alignment limit, one of the two $ {\cal{C}}{\cal{P}} $-even mass eigenstates can be regarded as the SM Higgs boson $ h^0_{\text{SM}} $, while the other is a BSM $ {\cal{C}}{\cal{P}} $-even Higgs boson denoted by $ \phi^0 $. We can see from Table 3 that
$ (h_{{\rm{SM}}}^0 , {\phi ^0}) = \left\{ {\begin{array}{*{20}{l}}{({h^0}, {H^0}) ,} &{({\rm{alignment}} \;{\rm{limit}} \;{\rm{I:}}}&{\sin (\beta - \alpha ) = 1),}\\{({H^0}, {h^0}) ,} &{({\rm{alignment}} \;{\rm{limit}} \;{\rm{II:}}}&{\cos (\beta - \alpha ) = 1).}\end{array}} \right. $
(9)
Given that the $ h^0A^0Z^0 $ and $ H^0A^0Z^0 $ coupling strengths are proportional to $ \cos(\beta - \alpha) $ and $ \sin(\beta - \alpha) $, respectively, the $ h^0_{\text{SM}}A^0 $ associated production is forbidden up to $ {\cal{O}}(\alpha^2\alpha_s) $ at the alignment limit. Thus, we only focus on the Drell-Yan production of $ \phi^0A^0 $ in the following.
The doubly-differential cross section for $ pp \rightarrow $$ \phi^0A^0 + X $ can be perturbatively calculated by means of the QCD factorization theorem:
$ \begin{aligned}[b] M^2 \frac{{\rm d}^2\sigma}{{\rm d}M^2 {\rm d}p_T^2}(\tau) =& \sum\limits_{a,b} \int_0^1 {\rm d}x_a {\rm d}x_b {\rm d}z \Big[ x_a f_{a/P}(x_a, \mu_F^2) \Big] \\&\times \Big[ x_b f_{b/P}(x_b, \mu_F^2) \Big] \\ & \times \Big[ z \hat{\sigma}_{ab}(z, M^2, p_T^2, \mu_F^2, \mu_R^2) \Big] \\&\times\delta \big( \tau - x_a x_b z \big), \end{aligned} $
(10)
where M and $ p_T $ are the invariant mass and transverse momentum of the final-state $ \phi^0A^0 $ system, respectively. The threshold variables τ and $ z $ in Eq. (10) are defined by
$\begin{aligned}[b] \tau =& \big( M/\sqrt{s} \big)^2, \;\;\;\;z = \big( M/\sqrt{\hat{s}} \big)^2,\end{aligned} $
(11)
where $ \sqrt{s} $ and $ \sqrt{\hat{s}} $ denote the hadronic and partonic center-of-mass energies, respectively. The universal PDF $ f_{a/P}(x, \mu_F^2) $ gives the probability to find parton a in proton P at factorization scale $ \mu_F $ as a function of fraction x of the proton's longitudinal momentum carried by the parton. After preforming a Mellin transformation,
$ F(N) = \int_0^1 {\rm d}y y^{N-1} F(y), $
(12)
on Eq. (10), the hadronic cross section can be written as a simple product of the PDFs and the partonic cross section in the conjugate Mellin N-space as
$\begin{aligned}[b] M^2 \frac{{\rm d}^2\sigma}{{\rm d}M^2{\rm d}p_T^2}(N-1) =& \sum\limits_{a,b} f_{a/P}(N, \mu_F^2) f_{b/P}(N, \mu_F^2)\\&\times \hat{\sigma}_{ab}(N, M^2, p_T^2, \mu_F^2, \mu_R^2).\end{aligned} $
(13)
To be consistent with the CT collaboration [24], we refit the PDF, $ f_{a/P}(x, \mu_F^2) $, as a polynomial of $ x^{1/2} $ with eight coefficients,
$\begin{aligned}[b] f_{a/P}(x, \mu_F^2) =& A_0 x^{A_1} \left( 1 - x \right)^{A_2} \Big( 1 + A_3 x^{1/2} + A_4 x \\&+ A_5 x^{3/2} + A_6 x^{2} + A_7 x^{5/2} \Big). \end{aligned}$
(14)
Thus, the Mellin moment of the PDF has the form
$ \begin{aligned}[b] f_{a/P}(N, \mu_F) = &A_0 \Big[ \text{B}(A_1+N, A_2+1) + A_3 \text{B} (A_1+N+1/2, A_2+1) \\ & + A_4 \text{B} (A_1+N+1, A_2+1) + A_5 \text{B} (A_1+N+3/2, A_2+1) \\& + A_6 \text{B} (A_1+N+2, A_2+1) + A_7 \text{B} (A_1+N+5/2, A_2+1) \Big], \end{aligned} $
(15)
where $ \text{B}(x, y) \equiv \Gamma(x)\Gamma(y)/\Gamma(x+y) $ is the Beta function.
According to the factorization scheme presented in Ref. [25], the partonic cross section can be expressed as a product of a process-dependent hard function and a process-independent Sudakov exponential term [13, 14, 26-29]. The higher-order QCD contributions to the partonic cross section contain logarithmic terms of type $ \alpha_s^n (M^2/p_T^2) \ln^m (M^2/p_T^2) $, which become large in the small- $ p_T $ region. These logarithmically-enhanced contributions arising at small $ p_T $ spoil the convergence of the fixed-order perturbative expansion and must therefore be resummed to all orders in $ \alpha_s $. We adopt the b-space resummation approach, which was fully formulated by Collins, Soper, and Sterman [9-12], to systematically resum the large logarithmic terms at small $ p_T $. In this approach, a Bessel transform is applied to the partonic cross section,
$ \begin{aligned}[b]\hat{\sigma}_{ab}(N, M^2, p_T^2, \mu_F^2, \mu_R^2) =& \int_0^{\infty} {\rm d}b \frac{b}{2} J_0(b p_T)\\&\times \hat{\sigma}_{ab}(N, M^2, b^2, \mu_F^2, \mu_R^2), \end{aligned}$
(16)
where $ J_0(x) $ is the zeroth-order Bessel function. Given that impact parameter b is the variable conjugate to transverse momentum $ p_T $, the limit $ M/p_T \rightarrow 0 $ corresponds to $ M b \rightarrow \infty $. Therefore, the large logarithms of $ M/p_T $ arising at small $ p_T $ turn into large logarithms of Mb,
$ \Big[ (M^2/p_T^2) \ln^m (M^2/p_T^2) \Big]_+ \longrightarrow \ln^{m+1} (M^2b^2) + \cdots $
(17)
After performing the resummation procedure, the resummed partonic cross section in the conjugate b-space at the NLL accuracy can be expressed as [9-12]
$ \begin{aligned}[b] \hat{\sigma}_{ab}^{\text{(res.)}}(N, M^2, b^2, \mu_F^2, \mu_R^2) =& \sum\limits_{a^{\prime}, a^{\prime\prime}, b^{\prime}, b^{\prime\prime}} E_{a^{\prime} a}^{(1)}(N, 1/\bar{b}^2, \mu_F^2) \\&\times E_{b^{\prime} b}^{(1)}(N, 1/\bar{b}^2, \mu_F^2) {\cal{C}}_{a^{\prime\prime} a^{\prime}}(N, 1/\bar{b}^2)\\&\times {\cal{C}}_{b^{\prime\prime} b^{\prime}}(N, 1/\bar{b}^2) {\cal{H}}_{a^{\prime\prime} b^{\prime\prime}}(M^2, \mu_R^2) \\&\times\exp \Big[ {\cal{G}}_{a^{\prime\prime} b^{\prime\prime}}(M^2\bar{b}^2, M^2, \mu_R^2) \Big], \end{aligned} $
(18)
where $ \bar{b} $ is the normalized impact parameter defined by $ \bar{b} = b/b_0 $ with $ b_0 = 2 {\rm e}^{-\gamma_{\text{E}}} $ [17], and the one-loop QCD evolution operator $ E^{(1)}_{ab} $ is derived from the collinear-improved procedure as recommended in Refs. [30-34]. In the physical resummation scheme, the coefficient function $ {\cal{C}}_{ab} $ and the Sudakov form factor $ {\cal{G}}_{ab} $ are free from any hard contributions, and the hard function $ {\cal{H}}_{ab} $, determined by the finite part of the renormalized virtual contribution, is free from any logarithmic contributions [16]. To transform the resummed partonic cross section $ \hat{\sigma}_{ab}^{\text{(res.)}}(N, M^2, b^2, \mu_F^2, \mu_R^2) $ back to the physical $ p_T $-space, we rewrite Eq. (16) as [35]
$\begin{aligned}[b] \hat{\sigma}_{ab}^{\text{(res.)}}(N, M^2, p_T^2, \mu_F^2, \mu_R^2) =& \sum\limits_{k = 1,2}\, \intop_{C_k} {\rm d}b \frac{b}{4} h_k(b p_T, v)\\&\times \hat{\sigma}_{ab}^{\text{(res.)}}(N, M^2, b^2, \mu_F^2, \mu_R^2), \end{aligned}$
(19)
where
$ h_k(x, v) = \frac{(-1)^k}{\pi} \int_{-iv\pi}^{(-1)^k \pi + iv\pi} {\rm d}\theta {\rm e}^{-{\rm i} x \sin\theta}, \qquad (k = 1, 2), $
(20)
are two auxiliary Hankel-like functions satisfying $ h_1(x, v) + h_2(x, v) = 2 J_0(x) $, and the integration contours $ C_k\; (k = 1, 2) $ in the complex b-plane are defined by
$ C_k: \quad b = b(t) \equiv t {\rm e}^{(-1)^k i \varphi}, \qquad t \in [0, +\infty) \quad\text{with}\quad \varphi \in (0, \pi/2). $
(21)
It is well known that such contours avoid the Landau pole by a deformation into either the upper or lower half complex b-plane.
The invariant mass distribution of the final-state $ \phi^0A^0 $ system in the Mellin N-space can be obtained by integrating Eq. (13) over the transverse momentum $ p_T $. In the threshold regime, the large logarithmic terms of the type $ \alpha_s^n (1-z)^{-1} \ln^m (1-z) $ also spoil the convergence of the perturbative series. These singular terms turn into large logarithms of the Mellin variable, N:
$ \Big[ (1-z)^{-1} \ln^m (1-z) \Big]_+ \longrightarrow \ln^{m+1} N + \cdots .$
(22)
The corresponding resummed partonic cross section for invariant mass distribution at the NLL accuracy can be expressed as [36, 37]
$ \begin{aligned}[b] \hat{\sigma}_{ab}^{\text{(res.)}}(N, M^2, \mu_F^2, \mu_R^2) =& \sum\limits_{a^{\prime}, b^{\prime}} E_{a^{\prime} a}^{(1)}(N, M^2/\bar{N}^2, \mu_F^2)\\&\times E_{b^{\prime} b}^{(1)}(N, M^2/\bar{N}^2, \mu_F^2) \tilde{{\cal{H}}}_{a^{\prime} b^{\prime}}(M^2, \mu_R^2)\\&\times \exp \Big[ \tilde{{\cal{G}}}_{a^{\prime} b^{\prime}}(\bar{N}, M^2, \mu_R^2) \Big], \end{aligned} $
(23)
where $ \bar{N} $ is the reduced Mellin variable defined by $ \bar{N} = N {\rm e}^{\gamma_{\text{E}}} $.
The hard functions, $ {\cal{H}}_{ab} $ and $ \tilde{{\cal{H}}}_{ab} $, do not contain any large logarithms. They can be perturbatively calculated and read at the NLO accuracy:
$ \begin{aligned}[b] &{\cal{H}}_{ab}(M^2, \mu_R^2) = \hat{\sigma}^{(0)}_{ab}(M^2) \left( 1+a_s {\cal{A}}_0 \right), \\ &\tilde{{\cal{H}}}_{ab}(M^2, \mu_R^2) = {\cal{H}}_{ab}(M^2, \mu_R^2) + a_s \frac{\pi^2}{6} \Big[ A_a^{(1)} + A_b^{(1)} \Big] \hat{\sigma}^{(0)}_{ab}(M^2), \end{aligned} $
(24)
where $ a_s = \alpha_s/(2\pi) $, $ A_a^{(1)} = 2 C_a $, $ \hat{\sigma}_{ab}^{(0)} $ is the lowest-order partonic cross section, and $ {\cal{A}}_0 $ represents the IR-finite part of the renormalized virtual correction in the dimensional regularization scheme, i.e.,
$\begin{aligned}[b] \hat{\sigma}_{ab}^{\text{(vir.)}}(M^2, \mu_R^2) = &a_s \left( \frac{4\pi\mu_R^2}{M^2} \right)^{\epsilon} \frac{\Gamma(1 - \epsilon)}{\Gamma(1 - 2 \epsilon)} \left( \frac{{\cal{A}}_{-2}}{\epsilon^2}\right.\\&\left. + \frac{{\cal{A}}_{-1}}{\epsilon} + {\cal{A}}_{0} \right) \hat{\sigma}_{ab}^{(0)}(M^2) + {\cal{O}}(\epsilon).\end{aligned} $
(25)
The Sudakov form factors $ {\cal{G}}_{ab} $ and $ \tilde{{\cal{G}}}_{ab} $ collect all the logarithmically-enhanced contributions and take the form
$\begin{aligned}[b]& {\mathbb{G}}_{ab}(\omega, M^2, \mu_R^2) = L {\mathbb{G}}_{ab}^{(1)}(\lambda) + \sum\limits_{n = 0}^{+\infty} a_s^{n} {\mathbb{G}}_{ab}^{(n+2)}(\lambda, M^2/\mu_R^2),\\ &({\mathbb{G}} = {\cal{G}} \; {\rm{or}}\; \tilde{{\cal{G}}}),\end{aligned} $
(26)
with $ \lambda = a_s \beta_0 L $, $ L = \ln \omega $, and $ \omega = M^2\bar{b}^2 $ and $ \bar{N} $ for $ {\mathbb{G}} = {\cal{G}} $ and $ \tilde{{\cal{G}}} $, respectively. The function $ {\mathbb{G}}_{ab}^{(n+1)}\; (n = 0,1,2,...) $ on the right side of Eq. (26) resums all the $ \text{N}^{n}\text{LL} $ contributions. In this study, we only consider the LL and NLL terms, i.e., $ L {\mathbb{G}}_{ab}^{(1)} $ and $ {\mathbb{G}}_{ab}^{(2)} $, because the electroweak production of $ \phi^0A^0 $ is studied at the NLO+NLL accuracy. The analytic expressions for $ {\mathbb{G}}_{ab}^{(1)} $ and $ {\mathbb{G}}_{ab}^{(2)} $ can be found in Ref. [38]. Finally, the $ {\cal{C}}_{ab} $ function in Eq. (18) at the NLL accuracy can be expressed as
$ {\cal{C}}_{ab}(N, \mu_R^2) = \delta_{ab} + a_s \left[ \frac{\pi^2}{6} C_{a} \delta_{ab} - P^{\prime}_{ab}(N) \right], $
(27)
where $ P_{ab}^{\prime}(N) $ is the $ {\cal{O}}(\epsilon) $ part of the unregulated Altarelli-Parisi splitting function in the Mellin N-space, i.e.,
$\begin{aligned}[b] &P_{ab}(z, \epsilon) = P_{ab}(z) + \epsilon P_{ab}^{\prime}(z),\\& P^{\prime}_{ab}(N) = \int_0^1 {\rm d}z z^{N-1} P^{\prime}_{ab}(z). \end{aligned}$
(28)
The resummed partonic cross section $ \hat{\sigma}_{ab}^{\text{(res.)}} $ gives the dominant contribution in the small-$ p_T $ and threshold regions, while the fixed-order partonic cross section $ \hat{\sigma}_{ab}^{\text{(f.o.)}} $ dominates at large $ p_T $ and small $ M/\sqrt{\hat{s}} $. To obtain a reliable theoretical prediction with uniform accuracy in all kinematical regions, the resummed and fixed-order results should be combined consistently by subtracting their overlap,
$ \hat{\sigma}_{ab} = \hat{\sigma}_{ab}^{\text{(res.)}} + \hat{\sigma}_{ab}^{\text{(f.o.)}} - \hat{\sigma}_{ab}^{\text{(o.l.)}}. $
(29)
This matching procedure guarantees that the combined result $ \hat{\sigma}_{ab} $ contains both the perturbative contributions up to the specific fixed order and the logarithmically-enhanced contributions from higher orders. At the NLO+NLL accuracy, $ \hat{\sigma}_{ab}^{\text{(o.l.)}} $ in Eq. (29) can be obtained by expanding the resummed partonic cross section $ \hat{\sigma}_{ab}^{\text{(res.)}} $ to $ {\cal{O}}(\alpha_s) $, i.e.,
$ \hat{\sigma}_{ab}^{\text{(o.l.)}} = \hat{\sigma}_{ab}^{\text{(res.)}} (\alpha_s = 0) + \alpha_s \frac{{\rm d} \hat{\sigma}_{ab}^{\text{(res.)}}}{{\rm d} \alpha_s} (\alpha_s = 0). $
(30)
After multiplying the Mellin moments of the PDFs to the NLO+NLL matched partonic cross section $ \hat{\sigma}_{ab} $, we obtain the hadronic differential cross section in the Mellin N-space. To get back to the physical space, an inverse Mellin transform,
$ F(\tau) = \frac{1}{2 \pi i} \int_{C_N} {\rm d}N \tau^{-N} F(N), $
(31)
should be applied to the right side of Eq. (13). To achieve this, we must comprehensively estimate the singularities in the Mellin N-space and choose an appropriate integration contour $ C_N $. There are two types of singularities for the hadronic differential cross section in the Mellin N-space: (1) the poles in the Mellin moments of the PDFs (Regge poles), and (2) the Landau pole related to the running of the strong coupling constant. The integration contour $ C_N $ in the complex N-plane is chosen as [15]
$ C_N: \quad N = N(y) \equiv C + y {\rm e}^{\pm {\rm i} \phi}, \qquad y \in [0, +\infty) $
(32)
where $ \phi \in [\pi/2, \pi) $ and the constant C is chosen such that the Regge and Landau poles lie to the left and right of $ C_N $, respectively.
In principle, for NLO+NLL calculations, we should employ resummation-improved PDFs for initial-state parton convolution. The threshold-resummation improved PDFs are now available with the NNPDF3.0 set. Compared to the NNPDF3.0 global fit, the threshold-resummation improved PDF fit has to be performed with a reduced data set involving deep-inelastic scattering, Drell-Yan, and top-pair production data, because the threshold resummation calculations are not readily available for all the processes employed in the global analysis. The reduced data set used in the fit of the threshold-resummation improved PDF set would induce a relatively larger PDF error compared to the global PDF set. In this study, we adopt the factorization method proposed in Ref. [18] to combine the smaller PDF error of the global PDF set with the resummation effect from the threshold-resummation improved PDF set. In the factorization method, the NLO+NLL QCD corrected cross section can be approximately calculated by
$ \sigma^{\text{NLO+NLL}} = K \times \sigma^{\text{NLO}}\big|_{{\rm (NLO \; global)}}\,, $
(33)
where
$ K = K_{\text{PDF}} \times K_{\text{PME}}\,, $
(34)
$\begin{aligned}[b] K_{\text{PDF}} =& \dfrac{\sigma^{\text{NLO+NLL}}\big|_{{\rm (NLO+NLL \; reduced)}}}{\sigma^{\text{NLO+NLL}}\big|_{{\rm (NLO \; reduced)}}}\,,\\ K_{\text{PME}} =& \dfrac{\sigma^{\text{NLO+NLL}}\big|_{{\rm (NLO \; global)}}}{\sigma^{\text{NLO}}\big|_{{\rm (NLO \; global)}}} \end{aligned}$
(35)
which describe the impact of the threshold-resummation improved PDFs and the NLL resummation effect from the partonic matrix element, respectively. Subscripts “NLO+NLL reduced ” and “NLO reduced ” appearing in the definition of $ K_{\text{PDF}} $ denote the threshold-resummation improved PDF set NNPDF30_nll_disdytop and its fixed-order version NNPDF30_nlo_disdytop [39], respectively. It is well known that the NNPDF cannot be properly transformed to the Mellin space; the refit of the NNPDF replicas in the Mellin space would lead to some convergence issues [40]. Fortunately, however, $ K_{\text{PME}} $ is expected to be largely independent of the PDF choice because the PDF sets used in $ K_{\text{PME}} $ are estimated at the same perturbative order [40]. This feature has been verified with the CT18NLO and MSTW2008nlo68cl PDFs, and thus we choose CT18NLO as the “NLO global” PDF set in our calculations.
2
B.QCD production via gluon-gluon fusion
-->

B.QCD production via gluon-gluon fusion

Compared to the electroweak production via quark-antiquark annihilation, the gluon-initiated QCD production of the scalar-pseudoscalar pair is a loop-induced production channel. This production mechanism should be taken into consideration at the LHC due to the high luminosity of gluon in proton. In Fig. 2 we depict some representative Feynman diagrams for $ gg \rightarrow \phi^0A^0 $ at the lowest order. Note that the production rate relies not only on the heavy-quark Yukawa couplings, but also on the triple Higgs self-couplings. Unlike the quark-antiquark annihilation channel, the loop-induced gluon-gluon fusion channel is extremely sensitive to the Yukwawa interaction of the 2HDM. Due to the introduction of a soft breaking $ Z_2 $ symmetry to avoid tree-level FCNCs, each fermion type is only able to couple to one of the two Higgs doublets. There are four allowed types of 2HDMs, type-I, type-II, lepton-specific, and flipped, which correspond to the four different types of Yukawa interaction. In this study, we mainly focused on the type-I 2HDM and calculated the gluon-gluon fusion channel by using the modified FeynArts-3.9, FormCalc-7.3, and LoopTools-2.8 packages [41-43].
Figure2. Representative Feynman diagrams for $ gg \rightarrow \phi^0A^0 $ at the lowest order.

IV.NUMERICAL RESULTS AND DISCUSSION
In this section, we provide some numerical results for $ pp \rightarrow \phi^0 A^0 + X $ at the $ 13\; \text{TeV} $ LHC in the type-I 2HDM. The SM input parameters used in this study are set as [44]
$ \begin{array}{l} m_W = 80.379 \; \text{GeV}, \; \; m_Z = 91.1876 \; \text{GeV}, \;\; m_t = 172.76 \; \text{GeV}, \\ G_{\rm F} = 1.1663787 \times 10^{-5} \; \text{GeV}^{-2}, \quad \alpha_s(m_Z) = 0.118. \end{array} $
(36)
The input parameters for the Higgs sector of the 2HDM should satisfy the theoretical constraints from perturbative unitarity [45], stability of vacuum [46], and tree-level unitarity [47], which can be checked by 2HDMC [19]. Moreover, the high-energy experiments can also give stringent constraints on the 2HDM input parameters. One of the experimental limits is that the partial width of $ Z^0 \rightarrow H^0A^0/h^0A^0 $ cannot exceed 2σ uncertainty of the Z-width measurement [44], and others come from the restriction of the physical observables of B meson decays, the measurement of the SM-like Higgs property, and the direct search of Higgs state at the LEP, Tevatron, and LHC, which are integrated in the SuperIso [48], HiggsSignals [49], and Higgsbounds [50] packages, respectively.
We use the CT14lo PDF set [24] to perform LO calculation, and employ the CT18NLO PDF [51] to obtain NLO and NLO+NLL QCD corrected cross sections. CT18NLO PDF contains 1 central PDF set and $ 2 \times 29 $ Hessian replicas. The PDF uncertainties of a cross section σ calculated with the CT18NLO PDF are given by [52]
$ \begin{aligned}[b] \delta^+_{\text{PDF}} =& \frac{1}{\sigma_0} \sqrt{ \sum\limits_{i = 1}^{29} \Big[ \max \big\{ \sigma_{i+} - \sigma_0,\, \sigma_{i-} - \sigma_0,\, 0 \big\} \Big]^2 }\,, \\ \delta^-_{\text{PDF}} =& \frac{1}{\sigma_0} \sqrt{ \sum\limits_{i = 1}^{29} \Big[ \max \big\{ \sigma_0 - \sigma_{i+},\, \sigma_0 - \sigma_{i-},\, 0 \big\} \Big]^2 }\,, \end{aligned} $
(37)
where $ \sigma_0 $ is the central value calculated with central set and $ \sigma_{i\pm}\; (i = 1, ..., 29) $ are the cross sections evaluated with replicas. The factorization scale $ \mu_F $ and the renormalization scale $ \mu_R $ are set to be equal, i.e., $ \mu_F = \mu_R = \mu $, for simplicity. The scale uncertainties of an integrated cross section σ are defined by
$ \begin{aligned}[b] \delta^+_{\mu} =& \frac{\max \big\{ \sigma(\mu)\, \big| \, \mu_0/2 \leqslant \mu \leqslant 2\mu_0 \big\} - \sigma(\mu_0)}{\sigma(\mu_0)}\,, \\ \delta^-_{\mu} =& \frac{\min \big\{ \sigma(\mu)\, \big| \, \mu_0/2 \leqslant \mu \leqslant 2\mu_0 \big\} - \sigma(\mu_0)}{\sigma(\mu_0)}\,, \end{aligned} $
(38)
where $ \mu_0 $ is the central scale. The total theoretical error is defined as the sum in quadrature of the PDF and scale uncertainties. For the quark-initiated Drell-Yan production channel, $ pp \rightarrow q\bar{q} \rightarrow \phi^0A^0 $, the production rate will be calculated at the NLO+NLL accuracy in QCD, and the NLO and NLO+NLL relative corrections are respectively defined as
$ \delta^{\text{NLO}} = \frac{\sigma^{\text{NLO}} - \sigma^{\text{LO}}}{\sigma^{\text{LO}}}, \quad \delta^{\text{NLO+NLL}} = \frac{\sigma^{\text{NLO+NLL}} - \sigma^{\text{LO}}}{\sigma^{\text{LO}}}. $
(39)
Regarding the gluon-gluon fusion channel, $pp \rightarrow gg \rightarrow $$ \phi^0A^0$, we only consider the lowest-order contribution since it is a loop-induced channel.
2
A.Integrated cross section
-->

A.Integrated cross section

In this subsection, we present the integrated cross sections for $ \phi^0A^0 $ associated production at $ \sqrt{s} = 13\; \text{TeV} $ LHC at the alignment limit in the 2HDM. The masses of $ \phi^0 $ and $ A^0 $ can be alternatively described by the following three parameters,
$\begin{aligned}[b]& m = \min\left( m_{\phi^0},\, m_{A^0} \right), \qquad \Delta m = \left| m_{\phi^0} - m_{A^0} \right|, \\&\epsilon = \text{sign} \left( m_{\phi^0} - m_{A^0} \right), \end{aligned}$
(40)
i.e., the minimal mass and mass hierarchy of $ \phi^0 $ and $ A^0 $. The two scenarios in which $ \epsilon = +1 $ ($ m_{\phi^0} > m_{A^0} $) and $ \epsilon = -1 $ ($ m_{\phi^0} < m_{A^0} $) may be referred to as the normal mass hierarchy and inverted mass hierarchy, respectively.
The Drell-Yan production channel depends only on the masses of $ \phi^0 $ and $ A^0 $. Furthermore, its integrated cross section is independent of $ \epsilon $. In our calculations, we set $ 0 $, 50, and $ 100\;\; \text{GeV} $ as three benchmark values of $ \Delta m $, which correspond, respectively, to the following three $ \phi^0 $-$ A^0 $ mass splitting scenarios:
● degenerate scenario: $ \Delta m = 0 $
● hierarchical scenario with small mass splitting: $ 0< \Delta m < m_Z $
● hierarchical scenario with large mass splitting: $ \Delta m > m_Z $
The LO, NLO, NLO+NLL QCD corrected integrated cross sections and the corresponding theoretical relative errors induced by the factorization/renormalization scale and PDFs for $ pp \rightarrow q\bar{q} \rightarrow \phi^0A^0 $ as functions of m for $ \Delta m = 0 $, 50 and $ 100\; \text{GeV} $ are given in Tables 4, 5, and 6, respectively. The central scale is set as $ \mu_0 = m_{\phi^0} + m_{A^0} $. There is no PDF-induced theoretical error for the LO cross section because the CT14lo PDF used in the LO calculation contains only one central set but no PDF replicas. To study the full NLL resummation effect and the impact of the threshold-resummation improved PDFs in our calculations, we also provide the factorization K-factors K and $ K_{\text{PDF}} $ in these tables. The NLO+NLL QCD relative correction $ \delta^{\text{NLO+NLL}} $ and the matrix-element-induced factorization K-factor $ K_{\text{PME}} $ can be straightforwardly calculated by using $ \delta^{\text{NLO}} $, K, and $ K_{\text{PDF}} $,
$m\,\text{[GeV]}$ $\sigma^{\text{LO} } \text{[fb]}$ $\sigma^{\text{NLO} } \text{[fb]}$ $\sigma^{\text{NLO+NLL} } \text{[fb]}$ $ \delta^{\text{NLO}} $[%] K $ K_{\text{PDF}} $
50 $ 4923.7_{-9.4{\text{%}}}^{+8.5{\text{%}}} $ $ 6717.8_{-0.4{\text{%}}-3.7{\text{%}}}^{+0.7{\text{%}}+2.8{\text{%}}} $ $ 6731.0_{-0.6{\text{%}}-3.7{\text{%}}}^{+0.0{\text{%}}+2.8{\text{%}}} $ $ 36.4 $ $ 1.002 $ $ 1.013 $
$ 100 $ $ 218.7_{-3.2{\text{%}}}^{+2.5{\text{%}}} $ $ 290.2_{-0.4{\text{%}}-3.8{\text{%}}}^{+0.9{\text{%}}+2.9{\text{%}}} $ $ 290.8_{-0.4{\text{%}}-3.8{\text{%}}}^{+0.0{\text{%}}+2.9{\text{%}}} $ $ 32.7 $ $ 1.002 $ $ 1.011 $
$ 150 $ $ 49.41_{-0.5{\text{%}}}^{+0.1{\text{%}}} $ $ 63.75_{-0.9{\text{%}}-4.3{\text{%}}}^{+1.4{\text{%}}+3.3{\text{%}}} $ $ 63.94_{-0.4{\text{%}}-4.3{\text{%}}}^{+0.0{\text{%}}+3.3{\text{%}}} $ $ 29.0 $ $ 1.003 $ $ 1.009 $
$ 200 $ $ 16.98_{-1.4{\text{%}}}^{+1.2{\text{%}}} $ $ 21.40_{-1.3{\text{%}}-4.7{\text{%}}}^{+1.6{\text{%}}+3.6{\text{%}}} $ $ 21.49_{-0.6{\text{%}}-4.7{\text{%}}}^{+0.1{\text{%}}+3.6{\text{%}}} $ $ 26.0 $ $ 1.004 $ $ 1.006 $
$ 300 $ $ 3.504_{-3.3{\text{%}}}^{+3.4{\text{%}}} $ $ 4.249_{-1.8{\text{%}}-5.7{\text{%}}}^{+1.9{\text{%}}+4.3{\text{%}}} $ $ 4.267_{-1.2{\text{%}}-5.7{\text{%}}}^{+0.8{\text{%}}+4.3{\text{%}}} $ $ 21.3 $ $ 1.004 $ $ 1.001 $
$ 400 $ $ 1.051_{-4.5{\text{%}}}^{+4.9{\text{%}}} $ $ 1.233_{-2.2{\text{%}}-6.3{\text{%}}}^{+2.2{\text{%}}+5.0{\text{%}}} $ $ 1.237_{-1.7{\text{%}}-6.3{\text{%}}}^{+1.5{\text{%}}+5.0{\text{%}}} $ $ 17.3 $ $ 1.003 $ $ 0.994 $
$ 500 $ $ 0.3848_{-5.4{\text{%}}}^{+6.0{\text{%}}} $ $ 0.4387_{-2.5{\text{%}}-7.1{\text{%}}}^{+2.4{\text{%}}+5.8{\text{%}}} $ $ 0.4392_{-2.2{\text{%}}-7.1{\text{%}}}^{+2.1{\text{%}}+5.8{\text{%}}} $ $ 14.0 $ $ 1.001 $ $ 0.985 $
$ 600 $ $ 0.1597_{-6.2{\text{%}}}^{+6.9{\text{%}}} $ $ 0.1773_{-2.8{\text{%}}-7.9{\text{%}}}^{+2.6{\text{%}}+6.6{\text{%}}} $ $ 0.1771_{-2.9{\text{%}}-7.9{\text{%}}}^{+2.6{\text{%}}+6.6{\text{%}}} $ $ 11.0 $ $ 0.999 $ $ 0.977 $
$ 700 $ $ 0.07213_{-6.9{\text{%}}}^{+7.8{\text{%}}} $ $ 0.07818_{-3.1{\text{%}}-8.8{\text{%}}}^{+2.8{\text{%}}+7.5{\text{%}}} $ $ 0.07770_{-3.5{\text{%}}-8.8{\text{%}}}^{+3.2{\text{%}}+7.5{\text{%}}} $ $ 8.39 $ $ 0.994 $ $ 0.966 $
$ 800 $ $ 0.03465_{-7.4{\text{%}}}^{+8.5{\text{%}}} $ $ 0.03670_{-3.3{\text{%}}-9.7{\text{%}}}^{+3.0{\text{%}}+8.4{\text{%}}} $ $ 0.03618_{-4.4{\text{%}}-9.7{\text{%}}}^{+3.7{\text{%}}+8.4{\text{%}}} $ $ 5.92 $ $ 0.986 $ $ 0.952 $


Table4.LO, NLO, NLO+NLL QCD corrected integrated cross sections, NLO QCD relative corrections, and factorization K-factors (K and $ K_{\text{PDF}} $) for $ pp \rightarrow q\bar{q} \rightarrow \phi^0A^0 $ at $ \sqrt{s} = 13\; \text{TeV} $ LHC within the 2HDM. The cross section central values are folded with the theoretical relative errors estimated from scale variation (first quote) and PDFs (second quote). The mass splitting between $ \phi^0 $ and $ A^0 $ is fixed to zero ($ \Delta m = 0 $).

$ m\; \text{[GeV]} $ $ \sigma^{\text{LO}}\; \text{[fb]} $ $ \sigma^{\text{NLO}}\; \text{[fb]} $ $ \sigma^{\text{NLO+NLL}}\; \text{[fb]} $ $ \delta^{\text{NLO}} $ [%] K $ K_{\text{PDF}} $
50 $ 611.1_{-5.3{\text{%}}}^{+4.5{\text{%}}} $ $ 824.5_{-0.0{\text{%}}-3.7{\text{%}}}^{+0.5{\text{%}}+2.8{\text{%}}} $ $ 825.8_{-0.4{\text{%}}-3.7{\text{%}}}^{+0.0{\text{%}}+2.8{\text{%}}} $ $ 34.9 $ $ 1.002 $ $ 1.012 $
$ 100 $ $ 93.88_{-1.6{\text{%}}}^{+1.1{\text{%}}} $ $ 122.7_{-0.7{\text{%}}-4.0{\text{%}}}^{+1.2{\text{%}}+3.1{\text{%}}} $ $ 123.0_{-0.3{\text{%}}-4.0{\text{%}}}^{+0.0{\text{%}}+3.1{\text{%}}} $ $ 30.7 $ $ 1.002 $ $ 1.010 $
$ 150 $ $ 27.61_{-0.7{\text{%}}}^{+0.4{\text{%}}} $ $ 35.19_{-1.1{\text{%}}-4.5{\text{%}}}^{+1.5{\text{%}}+3.5{\text{%}}} $ $ 35.31_{-0.5{\text{%}}-4.5{\text{%}}}^{+0.0{\text{%}}+3.5{\text{%}}} $ $ 27.5 $ $ 1.003 $ $ 1.007 $
$ 200 $ $ 10.77_{-2.0{\text{%}}}^{+1.9{\text{%}}} $ $ 13.43_{-1.4{\text{%}}-4.9{\text{%}}}^{+1.7{\text{%}}+3.8{\text{%}}} $ $ 13.48_{-0.8{\text{%}}-4.9{\text{%}}}^{+0.3{\text{%}}+3.8{\text{%}}} $ $ 24.7 $ $ 1.004 $ $ 1.005 $
$ 300 $ $ 2.517_{-3.6{\text{%}}}^{+3.8{\text{%}}} $ $ 3.026_{-1.9{\text{%}}-5.8{\text{%}}}^{+2.0{\text{%}}+4.5{\text{%}}} $ $ 3.038_{-1.3{\text{%}}-5.8{\text{%}}}^{+1.0{\text{%}}+4.5{\text{%}}} $ $ 20.2 $ $ 1.004 $ $ 0.999 $
$ 400 $ $ 0.8032_{-4.8{\text{%}}}^{+5.2{\text{%}}} $ $ 0.9353_{-2.3{\text{%}}-6.5{\text{%}}}^{+2.2{\text{%}}+5.2{\text{%}}} $ $ 0.9388_{-1.9{\text{%}}-6.5{\text{%}}}^{+1.6{\text{%}}+5.2{\text{%}}} $ $ 16.4 $ $ 1.004 $ $ 0.992 $
$ 500 $ $ 0.3053_{-5.6{\text{%}}}^{+6.2{\text{%}}} $ $ 0.3457_{-2.6{\text{%}}-7.3{\text{%}}}^{+2.4{\text{%}}+6.0{\text{%}}} $ $ 0.3458_{-2.5{\text{%}}-7.3{\text{%}}}^{+2.2{\text{%}}+6.0{\text{%}}} $ $ 13.2 $ $ 1.000 $ $ 0.983 $
$ 600 $ $ 0.1299_{-6.4{\text{%}}}^{+7.2{\text{%}}} $ $ 0.1433_{-2.8{\text{%}}-8.1{\text{%}}}^{+2.7{\text{%}}+6.8{\text{%}}} $ $ 0.1430_{-3.1{\text{%}}-8.1{\text{%}}}^{+2.8{\text{%}}+6.8{\text{%}}} $ $ 10.3 $ $ 0.998 $ $ 0.974 $
$ 700 $ $ 0.05969_{-7.0{\text{%}}}^{+7.9{\text{%}}} $ $ 0.06432_{-3.1{\text{%}}-9.0{\text{%}}}^{+2.9{\text{%}}+7.7{\text{%}}} $ $ 0.06361_{-3.3{\text{%}}-9.0{\text{%}}}^{+3.7{\text{%}}+7.7{\text{%}}} $ $ 7.76 $ $ 0.989 $ $ 0.960 $
$ 800 $ $ 0.02904_{-7.6{\text{%}}}^{+8.6{\text{%}}} $ $ 0.03059_{-3.4{\text{%}}-9.9{\text{%}}}^{+3.1{\text{%}}+8.7{\text{%}}} $ $ 0.03008_{-4.4{\text{%}}-9.9{\text{%}}}^{+3.9{\text{%}}+8.7{\text{%}}} $ $ 5.34 $ $ 0.983 $ $ 0.949 $


Table5.Same as Table IV but for $ \Delta m = 50\; \text{GeV} $.

$m \text{[GeV]}$ $\sigma^{\text{LO} } \text{[fb]}$ $\sigma^{\text{NLO} } \text{[fb]}$ $\sigma^{\text{NLO+NLL} } \text{[fb]}$ $ \delta^{\text{NLO}}$[%] K $ K_{\text{PDF}} $
50 $ 185.0_{-3.0{\text{%}}}^{+2.4{\text{%}}} $ $ 245.3_{-0.4{\text{%}}-3.9{\text{%}}}^{+1.0{\text{%}}+2.9{\text{%}}} $ $ 245.7_{-0.4{\text{%}}-3.9{\text{%}}}^{+0.0{\text{%}}+2.9{\text{%}}} $ $ 32.6 $ $ 1.002 $ $ 1.011 $
$ 100 $ $ 45.97_{-0.4{\text{%}}}^{+0.1{\text{%}}} $ $ 59.26_{-1.0{\text{%}}-4.3{\text{%}}}^{+1.4{\text{%}}+3.3{\text{%}}} $ $ 59.43_{-0.4{\text{%}}-4.3{\text{%}}}^{+0.0{\text{%}}+3.3{\text{%}}} $ $ 28.9 $ $ 1.003 $ $ 1.009 $
$ 150 $ $ 16.30_{-1.4{\text{%}}}^{+1.2{\text{%}}} $ $ 20.53_{-1.4{\text{%}}-4.8{\text{%}}}^{+1.5{\text{%}}+3.6{\text{%}}} $ $ 20.60_{-0.7{\text{%}}-4.8{\text{%}}}^{+0.1{\text{%}}+3.6{\text{%}}} $ $ 26.0 $ $ 1.003 $ $ 1.006 $
$ 200 $ $ 7.033_{-2.5{\text{%}}}^{+2.4{\text{%}}} $ $ 8.683_{-1.6{\text{%}}-5.2{\text{%}}}^{+1.8{\text{%}}+4.0{\text{%}}} $ $ 8.717_{-0.9{\text{%}}-5.2{\text{%}}}^{+0.5{\text{%}}+4.0{\text{%}}} $ $ 23.5 $ $ 1.004 $ $ 1.004 $
$ 300 $ $ 1.832_{-3.9{\text{%}}}^{+4.2{\text{%}}} $ $ 2.183_{-2.0{\text{%}}-5.9{\text{%}}}^{+2.1{\text{%}}+4.7{\text{%}}} $ $ 2.192_{-1.5{\text{%}}-5.9{\text{%}}}^{+1.1{\text{%}}+4.7{\text{%}}} $ $ 19.2 $ $ 1.004 $ $ 0.998 $
$ 400 $ $ 0.6183_{-5.0{\text{%}}}^{+5.5{\text{%}}} $ $ 0.7146_{-2.3{\text{%}}-6.7{\text{%}}}^{+2.3{\text{%}}+5.4{\text{%}}} $ $ 0.7168_{-2.1{\text{%}}-6.7{\text{%}}}^{+1.8{\text{%}}+5.4{\text{%}}} $ $ 15.6 $ $ 1.003 $ $ 0.990 $
$ 500 $ $ 0.2433_{-5.8{\text{%}}}^{+6.5{\text{%}}} $ $ 0.2736_{-2.6{\text{%}}-7.5{\text{%}}}^{+2.5{\text{%}}+6.2{\text{%}}} $ $ 0.2738_{-2.7{\text{%}}-7.5{\text{%}}}^{+2.3{\text{%}}+6.2{\text{%}}} $ $ 12.5 $ $ 1.001 $ $ 0.982 $
$ 600 $ $ 0.1059_{-6.5{\text{%}}}^{+7.4{\text{%}}} $ $ 0.1161_{-2.9{\text{%}}-8.4{\text{%}}}^{+2.7{\text{%}}+7.1{\text{%}}} $ $ 0.1157_{-3.2{\text{%}}-8.4{\text{%}}}^{+2.9{\text{%}}+7.1{\text{%}}} $ $ 9.63 $ $ 0.997 $ $ 0.972 $
$ 700 $ $ 0.04949_{-7.2{\text{%}}}^{+8.1{\text{%}}} $ $ 0.05301_{-3.2{\text{%}}-9.2{\text{%}}}^{+2.9{\text{%}}+7.9{\text{%}}} $ $ 0.05250_{-3.8{\text{%}}-9.2{\text{%}}}^{+3.4{\text{%}}+7.9{\text{%}}} $ $ 7.11 $ $ 0.990 $ $ 0.959 $
$ 800 $ $ 0.02438_{-7.7{\text{%}}}^{+8.8{\text{%}}} $ $ 0.02554_{-3.4{\text{%}}-10.1{\text{%}}}^{+3.1{\text{%}}+9.0{\text{%}}} $ $ 0.02505_{-4.4{\text{%}}-10.1{\text{%}}}^{+3.9{\text{%}}+9.0{\text{%}}} $ $ 4.76 $ $ 0.981 $ $ 0.945 $


Table6.Same as Table IV but for $ \Delta m = 100\; \text{GeV} $.

$ \delta^{\text{NLO+NLL}} = \left( K - 1 \right) + K \delta^{\text{NLO}}, \qquad K_{\text{PME}} = \frac{K}{K_{\text{PDF}}}. $
(41)
As shown in Tables 4, 5, and 6, the QCD correction can significantly enhance the LO production cross section, especially for light scalar-pseudoscalar pair. The NLO QCD relative correction exceeds 30% at $ m = 50\; \text{GeV} $, and decreases gradually to approximately 5.9%, 5.3%, and 4.8% for $ \Delta m = 0 $, 50, and $ 100\; \text{GeV} $, respectively, as m increases to $ 800\; \text{GeV} $. The full NLL resummation correction (quantitatively described by $ K-1 $) slightly enhances the NLO QCD corrected cross section as $ m < 500\; \text{GeV} $, but suppresses it by approximately 2% in the high mass region. Compared to the full NLL resummation correction, the contribution from the threshold-resummation improved PDFs is more sensitive to m. The corresponding relative correction, i.e., $ K_{\text{PDF}} - 1 $, decreases monotonically as the increase of m, and reaches approximately –5% when $ m = 800\; \text{GeV} $. Moreover, we can find that the impact of the threshold-resummation improved PDFs becomes increasingly important with the increment of m for heavy scalar-pseudoscalar pair. On the contrary, the NLL QCD relative correction from the partonic matrix element, $ K_{\text{PME}} - 1 $, increases monotonically with the increment of m. It is almost independent of the mass splitting between $ \phi^0 $ and $ A^0 $, varying from approximately –1% to 4% as m increases from 50 to $ 800\; \text{GeV} $.
The QCD production of $ \phi^0A^0 $ via gluon-gluon fusion depends not only on $ m_{\phi^0} $ and $ m_{A^0} $, but also on $ m_{12}^2 $ and $ \tan\beta $, since the Yukawa couplings and triple Higgs self-couplings are involved in this production channel. We calculate the lowest-order production cross section for this loop-induced channel at the two benchmark points listed in Table 7, which can satisfy both theoretical and experimental constraints. At both benchmark points, $ H^0 $ is the BSM $ {\cal{C}}{\cal{P}} $-even Higgs boson, i.e., $ H^0 = \phi^0 $, and $ \sin(\beta - \alpha) = 1 $ at the alignment limit. The other two Higgs parameters of 2HDM, $ m_{h^0} $ and $ m_{H^{\pm}} $, are not given in Table 7, because the scalar-pseudoscalar pair production at QCD NLO+NLL accuracy is completely independent from the SM-like and charged Higgs bosons. The integrated cross sections for $ pp \rightarrow gg \rightarrow \phi^0A^0 $ at $ \sqrt{s} = 13\; \text{TeV} $ LHC in the type-I 2HDM at BP1 and BP2, listed in Table 8, are approximately 1 and $ 0.06\; \text{fb} $, respectively. We can see that $ \sigma_{gg}/\sigma^{\text{NLO+NLL}} $, i.e., the ratio of the contribution from gluon-gluon fusion channel to the NLO+NLL QCD corrected cross section of quark-antiquark annihilation channel, is approximately 2.9% at BP1 and can reach 8.0% at BP2. It can be concluded that the scalar-pseudoscalar pair production at the LHC in the type-I 2HDM is predominated by the quark-initiated Drell-Yan production channel, and the gluon-gluon fusion contribution is non-negligible and should be taken into consideration in precision predictions.
Benchmark point $ m_{H^0} \text{/GeV} $ $m_{A^0} \text{/GeV}$ $ m_{12}^2 $ $ \tan\beta $
BP1 150 200 2000 10
BP2 400 500 50000 2


Table7.Benchmark points BP1 and BP2.

Benchmark point BP1 BP2
$\sigma_{gg} \text{/fb}$ $ 1.020_{-19.8{\text{%}}-3.4{\text{%}}}^{+26.4{\text{%}}+3.7{\text{%}}} $ $ 0.05742^{+31.8{\text{%}}+8.0{\text{%}}}_{-22.7{\text{%}}-6.4{\text{%}}} $


Table8.Lowest-order integrated cross sections for $pp \rightarrow $$ gg \rightarrow \phi^0A^0$ at $ \sqrt{s} = 13\; \text{TeV} $ LHC in type-I 2HDM at the benchmark points BP1 and BP2.

2
B.Transverse momentum distribution
-->

B.Transverse momentum distribution

Next, we address the transverse momentum distribution of the scalar-pseudoscalar pair produced at the LHC. Since the one-loop-induced gluon-gluon fusion channel does not contribute to the $ p_T $ distribution due to the momentum conservation, we consider only the quark-initiated Drell-Yan production channel. The NLO, NLO+NLL QCD corrected transverse momentum distributions of $ \phi^0A^0 $ as well as the overlap between the NLO QCD corrected and NLL QCD resummed $ p_T $ distributions (labeled by “NLO”, “NLO+NLL”, and “OVERLAP”) for the Drell-Yan production of $ \phi^0A^0 $ at $ \sqrt{s} = 13\; \text{TeV} $ LHC in the 2HDM at the benchmark points BP1 and BP2 are shown in Figs. 3(a) and 4(a), respectively. The central scale is $ \mu_0 = m_{\phi^0} + m_{A^0} $. As expected, the NLO QCD corrected $ p_T $ distribution and the overlap $ p_T $ distribution are in good agreement with each other in the small-$ p_T $ region [16] and become divergent as $ p_T \rightarrow 0 $, but the discrepancy between them becomes increasingly evident with the increment of $ p_T $. The relative discrepancy between the NLO QCD corrected and the overlap $ p_T $ distributions, defined as
Figure3. (color online) (a) Transverse momentum distribution of final-state $ \phi^0A^0 $ and (b) its scale uncertainty for $ pp \rightarrow q\bar{q} \rightarrow \phi^0A^0 $ at $ \sqrt{s} = 13\; \text{TeV} $ LHC within the 2HDM at the benchmark point BP1.

Figure4. (color online) Same as Fig. 3 but at BP2.

$ \eta = \left( \frac{{\rm d} \sigma^{\text{NLO}}}{{\rm d} p_T} - \frac{{\rm d} \sigma^{\text{OVERLAP}}}{{\rm d} p_T} \right) \left/ \frac{{\rm d} \sigma^{\text{NLO}}}{{\rm d} p_T} \right., $
(42)
can reach about 18.2% and 42.6% when $ p_T = 150 $ and $ 300\; \text{GeV} $ at BP1 and BP2, respectively. Compared to the NLO QCD corrected $ p_T $ distribution, the NLO+NLL QCD corrected $ p_T $ distribution is finite and more reliable in the whole final-state phase space. It increases sharply in the small-$ p_T $ region, reaches its maximum of around $ 1.9\; \text{fb/GeV} $ in the vicinity of $ p_T \sim 5.5\; \text{GeV} $, and then decreases approximately logarithmically with the increment of $ p_T $ at the benchmark point BP1. Concerning the benchmark point BP2, the NLO+NLL QCD corrected $ p_T $ distribution peaks at $ p_T \sim 7.5\; \text{GeV} $ and its maximum is approximately $ 0.027\; \text{fb/GeV} $.
The scale uncertainty of a differential distribution with respect to some kinematic variable x can be estimated by
$\begin{aligned}[b]& \delta_{\mu}(x) = \max \left\{ \frac{{\rm d} \sigma}{{\rm d} x}(\mu_1) - \frac{{\rm d} \sigma}{{\rm d} x}(\mu_2) \right\} \left/ \frac{{\rm d} \sigma}{{\rm d} x}(\mu_0) \right., \\&\mu_1,\, \mu_2 \in \left[ \mu_0/2,\; 2\mu_0 \right].\end{aligned} $
(43)
In Figs. 3(b) and 4(b), we plot the scale uncertainties of the NLO and NLO+NLL QCD corrected $ p_T $ distributions, denoted by $ \delta_{\mu}^{\text{NLO}} $ and $ \delta_{\mu}^{\text{NLO+NLL}} $, at BP1 and BP2, respectively. As shown in the lower panels of Figs. 3(b) and 4(b), the scale uncertainty of the NLO QCD corrected $ p_T $ distribution increases gradually, while the scale uncertainty of the NLO+NLL QCD corrected $ p_T $ distribution first decreases consistently before reaching its minimum and then increases monotonically, with the increment of $ p_T $. Some representative values of $ \delta_{\mu}^{\text{NLO}} $ and $ \delta_{\mu}^{\text{NLO+NLL}} $ are given in Table 9. This table, as well as Figs. 3(b) and 4(b), clearly shows that $ \delta_{\mu}^{\text{NLO+NLL}} $ is much less than $ \delta_{\mu}^{\text{NLO}} $, especially in the intermediate-$ p_T $ region. Thus, we conclude that the resummation of higher-order large logarithmic contributions can significantly improve the fixed-order prediction for the $ p_T $ distribution; the NLO+NLL QCD corrected $ p_T $ distribution is much more reliable in the whole $ p_T $ region compared to the NLO QCD corrected $ p_T $ distribution.
Benchmark
point
BP1 BP2
$ p_T = 1\; \text{GeV} $ $ p_T = 150\; \text{GeV} $ $ p_T = 2\; \text{GeV} $ $ p_T = 300\; \text{GeV} $
$ \delta_{\mu}^{\text{NLO}} $ $ 20.9 $% $ 27.2 $% $ 28.7 $% $ 34.0 $%
$ \delta_{\mu}^{\text{NLO+NLL}} $ $ 17.2 $% $ 16.3 $% $ 16.9 $% $ 15.7 $%
$ \text{min.} \simeq 1.6\%_{\quad \left( @\, p_T \,\sim\, 45\, \text{GeV} \right)} $ $ \text{min.} \simeq 2.0\%_{\quad \left( @\, p_T \,\sim\, 80\, \text{GeV} \right)} $


Table9.Scale uncertainties of NLO and NLO+NLL QCD corrected $ p_T $ distributions for $ pp \rightarrow q\bar{q} \rightarrow \phi^0A^0 $ at $ \sqrt{s} = 13\; \text{TeV} $ LHC within the 2HDM at BP1 and BP2 for some typical values of $ p_T $.

2
C.Invariant mass distribution
-->

C.Invariant mass distribution

In this subsection, we discuss the threshold resummation effect on the invariant mass distribution of the scalar-pseudoscalar pair produced at the $ 13\; \text{TeV} $ LHC in the type-I 2HDM. The central scale is set to the invariant mass of the final-state scalar-pseudoscalar pair, i.e., $ \mu_0 = M $. In the upper panels of Figs. 5(a) and 6(a), we depict the invariant mass distributions of the $ \phi^0A^0 $ system for both quark-initiated electroweak Drell-Yan production and gluon-initiated QCD production of $ \phi^0A^0 $ at BP1 and BP2, respectively. The corresponding NLO and NLO+NLL QCD relative corrections to the Drell-Yan production channel are provided in the lower panels. The $ \phi^0A^0 $ invariant mass distribution of the Drell-Yan channel increases rapidly near the production threshold, and then decreases consistently after reaching its maximum, with the increment of M. It peaks at $ M \sim 450\; \text{GeV} $ for BP1 and $ M \sim 1150\; \text{GeV} $ for BP2, respectively, at both LO and NLO+NLL accuracies. Compared to the Drell-Yan channel, the $ \phi^0A^0 $ invariant mass distribution of the gluon-gluon fusion channel is much smaller, and decreases more quickly with the increment of M. The ratio of the differential cross sections of the two channels, $ {\rm d}\sigma_{gg}/{\rm d}\sigma^{\text{NLO+NLL}} $, is approximately 8.1% at $ M = 400\; \text{GeV} $ for BP1 and 24.9% at $ M = 1000\; \text{GeV} $ for BP2, respectively, and approaches zero rapidly as the increase of M. It implies that the contribution from the gluon-gluon fusion channel is indispensable near the production threshold, but negligible in the high invariant mass region. The NLO and NLO+NLL QCD relative corrections ( $ \delta^{\text{NLO}} $ and $ \delta^{\text{NLO+NLL}} $) to the Drell-Yan channel decrease gradually with the increment of M. They decrease from 29.8% to 1.4% and from 31.0% to ?1.3%, respectively, as M increases from $ 400\; \text{GeV} $ to $ 3\; \text{TeV} $ at BP1, and vary correspondingly in the range of ?7.3% $ \sim$ 22.4% and ?23.9% $ \sim$ 25.9% as $ M \in [1,\, 4]\; \text{TeV} $ at BP2.
Figure5. (color online) (a) Invariant mass distribution of final-state $ \phi^0A^0 $ and (b) factorization K-factors (K and $ K_{\text{PDF}} $) as well as theoretical relative errors ($ \delta_{\mu} $, $ \delta_{\text{PDF}} $, and $ \delta_{\text{tot}} $) for $ \phi^0A^0 $ associated production at $ \sqrt{s} = 13\; \text{TeV} $ LHC in type-I 2HDM at the benchmark point BP1.

Figure6. (color online) Same as Fig. 5 but at BP2.

To further demonstrate the full NLL resummation effect and the impact of the threshold-resummation improved PDFs on the $ \phi^0A^0 $ invariant mass distribution of the Drell-Yan channel, we plot the factorization K-factors K and $ K_{\text{PDF}} $ as functions of M in Figs. 5(b) and 6(b) for BP1 and BP2, respectively. The theoretical errors from scale variation and PDFs as well as their combination, i.e., $ \delta_{\mu} $, $ \delta_{\text{PDF}} $, and $ \delta_{\text{tot}} $, are also displayed in these two figures. At the benchmark point BP1, K increases slowly from 1.01 to 1.04 as the increment of M from $ 400\; \text{GeV} $ to $ 1.7\; \text{TeV} $, and then gradually decreases to 0.97 as M increases to $ 3\; \text{TeV} $. The full NLL resummation correction enhances the NLO QCD corrected invariant mass distribution of $ \phi^0A^0 $ in the region of $ M < 2.7\; \text{TeV} $, but it would reduce the invariant mass distribution at sufficiently high invariant mass. However, $ K_{\text{PDF}} $, which quantitatively reflects the impact of the threshold-resummation improved PDFs, decreases consistently from 1.01 to 0.90 as M varies from $ 400\; \text{GeV} $ to $ 3\; \text{TeV} $. $ K_{\text{PME}} $, which describes the NLL resummation effect from the partonic matrix element and is calculated by $ K/K_{\text{PDF}} $, shows the opposite tendency compared to $ K_{\text{PDF}} $: it increases monotonically from $ 1.00 $ to $ 1.08 $ with the increment of M. At the benchmark point BP2, K is fairly stable in the range of $ 1\; \text{TeV} < M < 2\; \text{TeV} $; it reaches its maximum of around $ 1.04 $ at $ M \sim 1.8\; \text{TeV} $ and subsequently decreases to $ 0.83 $ as M increases to $ 4\; \text{TeV} $. Simultaneously, a global suppression induced by the threshold-resummation improved PDFs can be clearly observed in the invariant mass distribution. Such suppression effect is very small and could be neglected at relatively low invariant mass, but becomes more and more apparent as the increasing of M. At $ M = 4\; \text{TeV} $, $ K_{\text{PDF}} = 0.73 $; the contribution from the threshold-resummation improved PDFs is more notable at high invariant mass compared to the NLO QCD correction. On the contrary, $ K_{\text{PME}} $ increases from 1.02 to 1.13 as M increases from 1 to $ 4\; \text{TeV} $. In the high invariant mass region, the contribution from the threshold-resummation improved PDFs is the dominant correction compared to the NLO QCD correction and the NLL resummation correction from the partonic matrix element. For example, at $ M = 4\; \text{TeV} $, $ K_{\text{PDF}} - 1 = -27$%, $ \delta^{\text{NLO}} = $$ -7.3$% and $ K_{\text{PME}} - 1 = 13$%, respectively.
V.SUMMARY
Searching for BSM Higgs bosons is an important task at the LHC and future high-energy colliders. In this study, we comprehensively analyze the scalar-pseudoscalar pair production at the $ 13\; \text{TeV} $ LHC at the alignment limit in the type-I 2HDM. The Collins-Soper-Sterman resummation approach and the factorization method are employed to resum the NLL contributions and to evaluate the impact of the threshold-resummation improved PDFs, respectively, when addressing the quark-initiated Drell-Yan production channel. Both the integrated cross section and the differential distributions with respect to the transverse momentum and invariant mass of the produced scalar-pseudoscalar pair are provided. For quark-antiquark annihilation channel, the NLO QCD relative correction can exceed 30% in the low Higgs mass region, but decreases rapidly as the increment of $ m_{\phi^0} $ and $ m_{A^0} $. The relative correction induced by the threshold-resummation improved PDFs and the NLL resummation correction from the partonic matrix element, $ K_{\text{PDF}} - 1 $ and $ K_{\text{PME}} - 1 $, are insensitive to the mass splitting between $ \phi^0 $ and $ A^0 $, and decreases and increases respectively with the increment of Higgs mass. They could be neglected compared to the NLO QCD correction in the low invariant mass region, but become increasingly important with the increment of the invariant mass of $ \phi^0A^0 $, and can even predominate in the high invariant mass region. Moreover, the anomalous behavior of the NLO QCD corrected transverse momentum distribution in the small-$ p_T $ region can be resolved, and the scale uncertainty can be heavily reduced, especially in the intermediate-$ p_T $ region, by including the NLL resummation correction. Compared to the quark-initiated Drell-Yan channel, the contribution from the gluon-gluon fusion channel is negligible in the high invariant mass region, but indispensable and even comparable to the NLO QCD correction near the production threshold.
相关话题/Scalarpseudoscalar production Large