Charged Higgs production from polarized top-quark decay in the 2HDM considering the general-mass variable-flavor-number scheme

Charged Higgs bosons $H^\pm$ are predicted by some non-minimal Higgs scenarios, such as models containing Higgs triplets and two-Higgs-doublet models, so that the experimental observation of these bosons would indicate physics beyond the Standard Model. In the present work, we introduce a new channel to indirect search for the charged Higgses through the hadronic decay of polarized top quarks where a top quark decays into a charged Higgs $H^+$ and a bottom-flavored hadron $B$ via the hadronization process of the produced bottom quark, $t(\uparrow)\rightarrow H^++b(\to B+jet)$. To obtain the energy spectrum of produced $B$-hadrons we present, for the first time, an analytical expression for the ${\cal O}(\alpha_s)$ corrections to the differential decay width of the process $t\rightarrow H^+b$ in the presence of a massive b-quark in the General-Mass Variable-Flavor-Number Scheme (GM-VFNS). We find that the most reliable predictions for the B-hadron energy spectrum are made in the GM-VFN scheme, specifically, when the Type-II 2HDM scenario is concerned.


I. INTRODUCTION
There are many reasons, both from experimental observations and theoretical considerations, to expect physics beyond the Standard Model (BSM), such as the hierarchy problem, neutrino oscillations and dark matter. Many numerous attempts have been done and are still in progress to build new physics models which can explain these puzzles. Among them, some well-known examples are the Minimal Supersymmetric Standard Model (MSSM) [1] and the Two-Higgs-Doublet Models (2HDM) [2,3], so the latter is known as the simplest model. In many extensions of the Standard Model (SM) such as the 2HDM, the Higgs sector of the SM is enlarged typically by adding an extra doublet of complex Higgs fields. In the 2HDM, after spontaneous symmetry breaking, the two scalar Higgs doublets H 1 and H 2 yield three physical neutral Higgs bosons (h, H, A) and a pair of charged-Higgs bosons H ± [4]. The observation of a (singly-)charged-Higgs boson in the current and future runs of the Large Hadron Collider (LHC) would clearly indicate a definitive evidence of new physics beyond the SM.
From the element |V tb | ≈ 1 of the Cabibbo-Kobayashi-Maskawa (CKM) quark mixing matrix [5], it is found that the top quark is decaying dominantly through t → bH + in the 2HDM [2,3], providing that the top quark mass (m t ), bottom quark mass (m b ) and the charged-Higgs boson mass (m H + ) satisfy the condition: m t > m b + m H + . In this situation, one may expect measurable effects in the top quark decay width and decay distributions due to the H ± -propagator contributions, which are poten- * Electronic address: mmoosavi@yazd.ac.ir tially large in the decay chain t → bH + → b(τ + ν τ ). At the LHC, it is expected to have a cross section σ(pp → ttX) ≈ 1 (nb) at design energy √ s = 14 TeV [6]. With the LHC design luminosity of 10 34 cm −2 s −1 in each of the four experiments, one may expect to have a tt pair per second so that with this remarkable potential the LHC can be considered as a superlative top factory which allows one to search for the charged-Higgs bosons in the subsequent decay products of the top pairs tt → H ± W ∓ bb and tt → H ± H ∓ bb. New results of a search on the charged-Higgs bosons in the proton-proton collision at a center of mass energy of √ s = 8 TeV are reported by the CMS [7] and ATLAS [8] Collaborations at the LHC, and we shall discuss on these results and the restrictions imposed over the MSSM m H + − tan β parameter space in Section V. In our numerical analysis we restrict ourselves to the unexcluded regions of this parameter space.
As was mentioned, due to |V tb | ≈ 1 of the CKM matrix the top quark almost exclusively decays to b-quark via t → bH + , in the 2HDM. On the other hand, b-quarks hadronize, via b → B + X, before they decay so that the decay process t → BH + +X is of prime importance and it is an urgent task to predict its partial decay width as realistically and reliably as possible. Therefore, at the LHC of particular interest would be the distribution in the scaled B-hadron energy (x B ) in the top quark rest frame so that the study of this distribution is proposed as a new channel to indirect search for the charged Higgs bosons. To study the scaled energy spectrum of B-hadrons we determine the quantity dΓ(t → BH + + X)/dx B . According to the factorization theorem of the QCDimproved parton model [9], the B-hadron energy distribution can be determined by the convolution of the partonic differential decay width (dΓ/dx b ) of the subprocess t → bH + , with the nonperturbative fragmentation func- where, x b stands for the scaled-energy fraction of the bottom quark and the D B b -FF describes the splitting process b → B + X in which X refers to the unobserved final state particles. In (1), µ F and µ R are the factorization and renormalization scales, respectively, and the integral convolution is defined as (f ⊗g)(x) =´1 x dy/yf (y)g(x/y).
In [10], we studied the energy spectrum of the bottomflavored mesons through unpolarized top quark decays in the 2HDM at next-to-leading order (NLO) of the QCD radiative corrections. In [11], we studied the spindependent energy distribution of B-mesons produced through polarized top decays in the massless scheme or Zero-Mass Variable-Flavor-Number scheme (ZM-VFNs) [12] where the zero mass parton approximation is also applied to the bottom quark. This massless approximation simplifies the evaluation of the NLO QCD corrections largely but at the price of losing the accuracy of analysis. For example, by this approximation the results for the energy spectrum of B-mesons are independent of the model selected.
In the present work, we impose the effects of bottom quark mass on the spin-dependent energy spectrum of B-mesons employing the General-Mass Variable-Flavor-Number (GM-VFN) scheme in which we preserve the bquark mass from the beginning. Our calculations are done in a special helicity frame where the polarization direction of the top quark is evaluated with respect to the b-quark three-momentum p b . As will be shown, the results are different in two variants of the 2HDM and it is found that the NLO corrections with m b = 0 to be significant, specifically, when the type-II 2HDM is considered.
In the SM, due to |V tb | ≈ 1 the top quark decays dominantly through the decay mode t → bW + . In [13][14][15][16][17], we investigated the energy distribution of B-mesons produced in polarized and unpolarized top quark decays in the SM. For the unpolarized top decays, the total distribution of the B-hadron energy is obtained by the summation of two contributions due to the decay modes t → bH + (in the 2HDM) and t → bW The same result is valid for the polarized top quark decay as long as the spin direction of the polarized top quark is evaluated relative to the b-quark three-momentum, as in our present work. Thus, at the LHC any deviation of the B-meson energy spectrum from the SM predictions can be considered as a signal for the existence of charged Higgs. Although, the SM contribution is normally larger than the one coming from the 2HDM (this comparison is studied in Ref. [10]), but there is always a clear separation between the decay channels t → bW + and t → bH + in both the ttX pair production and the t/tX single top production at the LHC, this point is mentioned in [18]. As a last point; for the decay chain t → bH + → b(τ + ν τ ) one has Γ(t → bτ + ν τ ) = Γ(t → bH + ) × BR(H + → τ + ν τ ) in the narrow width approximation for the Higgs boson. For the numerical values applied in this work, the branching ratio BR(H + → τ + ν τ ) = 1 [19], to a very high accuracy. Therefore, the results presented in this work are also valid for the decay chain t(↑) → Bτ + ν τ + X.
This paper is organized as follows. In Sec. II, we introduce the general angular structure of the differential decay width in a specific helicity frame. In Sec. III, we present our analytical results of the O(α s ) QCD corrections to the tree-level rate of t(↑) → bH + in the fixed flavor number scheme. In Sec. IV, in a detailed discussion we describe the GM-VFN scheme by introducing the perturbative FF b → b. In Sec. V, our hadron level results in the GM-VFN scheme will be presented. In Sec. VI, our conclusions are summarized.

II. ANGULAR STRUCTURE OF DIFFERENTIAL DECAY RATE
Here, we concentrate on the decay process t(↑) → bH + in the general 2HDM, where H 1 and H 2 are the doublets whose vacuum expectation values give masses to the down and up type quarks, respectively, and a linear combination of the charged components of H 1 and H 2 also gives the physical charged Higgs H + (H + = cos βH + 2 − sin βH + 1 ). The parameter β is defined in the following.
Basically, the dynamics of the current-induced t → b transition is embodied in the hadronic tensor H µν ∝ , where s t denotes the spin of the top quark. At the NLO QCD radiative corrections, only two types of intermediate states are contributed; |X b >= |b > for the Born level term and O(α s ) one-loop contributions and |X b >= |b + g > for the O(α s ) tree graph contribution. In the SM, where one has t → bW + , the weak current is given by Generally, in models with two Higgs doublets and generic coupling to all the quarks, it is difficult to avoid tree level flavor changing neutral currents, therefore, we limit ourselves to models that naturally stop these problems by restricting the Higgs coupling. As may be found in [20], the first possibility is to have the doublet H 1 coupling to all bosons and the doublet H 2 coupling to all the quarks (called model I in the following). This model leads to the coupling factors as where, tan β = v 2 /v 1 is the ratio of the vacuum expectation values of the two electrically neutral components of the two Higgs doublets and the weak coupling factor g W is related to the Fermi's coupling constant by In the second possibility (called model II), the doublet H 2 couples to the right-chiral up-type quarks (u R , c R , t R ) while the H 1 couples to the right-chiral down-type quarks. In this model the coupling factors read These models are mostly known as Type-I and Type-II 2HDM scenarios. The type-II is the Higgs sector of the Minimal Supersymmetric Standard Model (MSSM) up to SUSY corrections [1]. Two other models (models III and IV) are also possible which are explained in our previous work [11] in detail. See also [21]. Here, we just mention that the analytical results presented for the partonic process t(↑) → bH + are the same both in models I and IV and also in models II and III.
In the rest frame of a polarized top quark decaying into a b-quark and a Higgs boson (and a gluon at NLO), the final-state particles (b, H + , g) define an event frame. Relative to this event plane, the polarization direction of top quark can be defined. In this work, we analyze the decay mode t(↑) → bH + (+g) in the rest frame of the polarized top quark where the three-momentum of the bottom quark points into the direction of the positiveẑ-axis and the polar angle θ P is defined as the angle between the polarization vector of top quark and theẑ-axis, see Fig. 1 of Ref. [11]. The general angular distribution of the differential decay width dΓ/dx b of a polarized top quark is given by the following form to clarify the correlation between the polarization of the top quark and its decay products where the polar angle θ P shows the spin orientation of the top quark relative to the event plane and P (0 ≤ P ≤ 1) is the magnitude of the top quark polarization. In the notation above, the first and second terms in the curly bracket correspond to the unpolarized and polarized differential decay rates, respectively. As usual, we define the scaled-energy fraction of the bottom quark as The O(α s ) radiative corrections to the unpolarized differential width dΓ unpol /dx b have been studied in [10], extensively, and the NLO QCD corrections to the polarized partial rate dΓ pol /dx b in the ZM-VFN scheme (with m b = 0) are studied in [11]. In the present work, we compute the NLO QCD radiative corrections to the polarized partial rate dΓ pol /dx b in the GM-VFN scheme where m b = 0 is considered from the beginning. These analytical results are new and presented for the first time.

A. Born term results
It is straightforward to compute the Born term contribution to the partial decay rate of the polarized top quark in the 2HDM in the presence of the b-quark mass. The Born term amplitude for the process t(↑) → bH + can be parameterized as M 0 =ū b (a + bγ 5 )u t , so for the squared amplitude one has: Then, the polarized tree-level decay width reads where λ = λ(x, y, z) = (x−y−z) 2 −4yz is the Källén function and the factor √ λ/(16πm t ) = P S 2 is the two-body phase space factor. This result is in complete agreement with Refs. [22,23]. The unpolarized Born-level rate can be found in our previous work [10]. Considering Eqs. (2) and (3), in (6) for the product of two coupling factors in the model I (type-I 2HDM scenario), one has and for the model II (type-II 2HDM scenario), Since m b ≪ m t , the bottom quark mass can always be safely neglected in the model I, while in the model II, the second term in (8) (8) can become as large as O(48%) and this order will be larger when tan β is increased. Therefore, the results in the type-II 2HDM scenario depend on the b-quark mass extremely, unless the low values of tan β(tan β ≤ 1) are applied, however, these small values of tan β are now excluded by the CMS and ATLAS experiments at the LHC.
In [11], we adopted the ZM-VFN scheme (with m b = 0) which is not suitable for the Type-II 2HDM scenario. There, we pointed out that our phenomenological predictions are restricted to the Type-I 2HDM. In the present work we retain the b-quark mass and extend our results to both models and shall compare them. In the following, in a detailed discussion we calculate the O(α s ) QCD corrections to the Born-level width and present an analytical expression for dΓ(t(↑) → BH + + X)/dx B at NLO in the GM-VFN scheme.

B. Virtual gluon corrections including counterterms and one-loop vertex correction
Basically, the one-loop virtual corrections to the tbH +vertex consist of both the infrared (IR) and the ultraviolet (UV) singularities. The UV-divergences appear when the integration region of the virtual gluon momentum goes to infinity and the IR-divergences arise from the softgluon singularities. Here, we adopt the on-shell massrenormalization scheme and all singularities are regularized by dimensional regularization in D = 4 − 2ǫ spacetime dimensions where ǫ ≪ 1. Using the dimensional regularization technique one obtains the well-defined integrals which are finite while all singularities are summarized in the ǫ. This is done by the replacement: where µ is an arbitrary reference mass which will be removed after summing all corrections up. For simplicity, we introduce the following abbreviations: where the scaled masses y = m 2 H + /m 2 t and R = m 2 b /m 2 t are defined. Choosing these notations, the tree-level total width (6) is simplified as Also, the normalized energy fraction x b (5) is given by Taking the above notations, the contribution of virtual corrections into the decay width reads where, |M vir | 2 = Spin (M † 0 M loop + M † loop M 0 ) and the factor P S 2 = Q/(8πm t ) is a two-body phase space factor, as in (6). The renormalized amplitude is now written as M loop =ū b (Λ ct +Λ l )u t , where Λ l arises from the one-loop vertex correction and Λ ct stands for the counter term. The counter term of the vertex consists of the mass and the wave function renormalizations of both the top and bottom quarks [22][23][24], as where, the wave function and the mass renormalization constants are expressed as In the equation above, m q (q = b, t) is the mass of the relevant quark, γ E = 0.5772 · · · is the Euler constant and C F = (N 2 c − 1)/(2N c ) = 4/3 for N c = 3 quark colors. Also, ǫ IR and ǫ UV represent infrared and ultraviolet singularities.
The real part of the one-loop vertex correction is given by with where, B 0 and C 0 functions are the Passarino-Veltman 2-point and 3-point integrals, respectively. The analytical form of these integrals can be found in Ref. [25]. All the ultraviolet divergences shall be canceled after summing all virtual corrections up but the infrared singularities (ǫ IR ) are remaining which are now labeled by ǫ. Putting everything together, for the virtual differential decay width normalized to the Born-level total rate (10) one has This result after integration over x b is in complete agreement with Ref. [22].

C. Real gluon radiative corrections
If we denote the polarization vector of the real gluon by ǫ(p g , λ), the O(α s ) real gluon emission (tree-graph) amplitude reads In the ZM-VFN scheme, the IR-divergences arise from the soft-and collinear gluon emissions while in the GM-VFN scheme there are no collinear divergences and all IR-singularities arise from the soft real gluon emission. As before, to regulate the IR-divergences we work in Ddimensions so that the contribution of real gluon emission into the polarized differential decay rate is given by where dP S 3 is the three-body phase space To compute the real differential decay rate dΓ real /dx b , in (19) the momentum of b-quark is fixed and over the energy of the gluon is integrated. The energy of gluon ranges from E min To achieve the correct finite terms in the rate 1/Γ 0P × dΓ real,pol /dx b , the Born widthΓ 0P (6) must be evaluated in the dimensional regularization at O(ǫ), i.e.Γ 0P → Γ 0P {1−ǫ 2 ln Q+γ E −ln(4πµ 2 F /m 2 t ) }. When integrating over the phase space, terms of the form (1 − x b ) −1−2ǫ arise which are due to the radiation of a soft gluon in top decay at NLO. Note, the E g → 0 limit corresponds to the limit x b → 1. Thus for a massive bottom quark, where x b,min = β, we apply the following expansion [26] where the plus distribution is defined as usual. Finally, the contribution of real gluon emission reads To obtain an analytic result for the polarized partial decay rate, by summing the tree level, the virtual and the real contributions, one has As is seen, all IR-singularities are canceled and the final result is free of singularities. In Ref. [22], authors considered a specific helicity coordinate system where the polarization vector of the top quark was evaluated relative to the Higgs boson threemomentum. In [27], we applied the same frame and obtained the polarized differential decay width dΓ pol /dx b at the parton-level in the ZM-VFN scheme to obtain the energy spectrum of B-hadrons (according to Eq. (1)). There, we showed that our analytical result for the parton-level differential decay rate dΓ pol /dx b is in agreement with [22] after integration over x b (0 ≤ x b ≤ 1). In this work, following our previous work [11], we considered a new helicity frame where the polarization vector of the top quark is evaluated relative to the b-quark threemomentum. Applying the same techniques, our results for the Born rate (10) and virtual corrections (17) are the same in both helicity frames but the real corrections (22) and, in conclusion, the NLO differential decay width (23) are different and completely new. In next section we present our reason for correctness of the obtained result.

IV. GENERAL-MASS VARIABLE-FLAVOR-NUMBER SCHEME
In this work, our main purpose is to evaluate the scaled-energy distribution of the B-hadron produced in the inclusive process t(↑) → BH + + X in the 2HDM. Therefore, we calculate the NLO decay width of the corresponding process differential in x B (dΓ pol /dx B ) in the GM-VFN scheme, where x B = E B /(m t S) is the scaledenergy fraction of the B-hadron (as for x b in (11)). In the top quark rest frame applied in our work, the B-hadron . Considering the factorization theorem (1), the B-hadron energy spectrum can be obtained by the convolution of the parton-level spectrum (23) with the nonperturbative fragmentation function (FF) D B b (z, µ F ). We will discuss about the FFs needed, later.
We now explain how the quantity dΓ pol (µ R , µ F )/dx b will have to be evaluated in the GM-VFN scheme. In Sec. III, we employed the Fixed-Flavor-Number (FFN) scheme which contains of the full m b dependence. In the FFN scheme, the large logarithmic singularities of the form (α s /π) ln(m 2 t /m 2 b ) spoil the convergence of the perturbative expansion when m b /m t → 0. The massive or GM-VFN scheme is devised to resum these large logarithms and to retain the whole nonlogarithmic m b dependence at the same time and this is achieved by introducing appropriate subtraction terms in the NLO FFN expression for dΓ pol /dx b . In this case, the NLO ZM-VFN result is exactly recovered in the limit m b /m t → 0. In the GM-VFN scheme, the subtraction term is constructed as where 1/Γ 0P × dΓ pol ZM /dx b is the partial decay rate computed in the ZM-VFN scheme [11], in which all information on the m b -dependence of dΓ pol /dx b is wasted. In conclusion, the GM-VFN result is obtained by subtracting the subtraction term from the FFN one [28,29], Taking the limit m b → 0 in Eq. (23), one obtains the following subtraction term This result coincides with the perturbative FF of the transition b → b [30] and is in consistency with the Collin's factorization theorem [9], which guarantees that the subtraction terms are universal. Thus, the result obtained in (26) ensures the correctness of our result shown in (23).

V. NUMERICAL ANALYSIS
In the MSSM, the mass of charged Higgses is restricted by m H ± > m W ± at tree-level [31], however, this restriction is not valid for some regions of parameter space after including radiative corrections. In the MSSM, m H ± is strongly correlated with the mass of other Higgs bosons. In Ref. [19] is mentioned that a charged Higgs with a mass range 80 GeV ≤ m H ± ≤ 160 GeV is a logical possibility and its effects should be searched for in the decay chain t → bH + → Bτ + ν τ + X. On the other side, the last results of a search for evidence of a light charged Higgs boson (m H < m t ) in 19.5 − 19.7f b −1 of proton-proton collision data recorded at √ s = 8 TeV are reported by the CMS [7] and the ATLAS [8] collaborations, using the τ + jets channel with a hadronically decaying τ lepton in the final state. According to Fig. 7 of Ref. [8], the large region in the MSSM m H + − tan β parameter space is excluded for m H + = 80 − 160 GeV and the only unexcluded regions of this parameter space include the charged Higgs masses as 90 ≤ m H + ≤ 100 GeV (for 6 < tan β < 10) and 140 ≤ m H + ≤ 160 GeV (for 3 < tan β < 21). See also figure 9 of Ref. [7]. Therefore, in this work our phenomenological predictions are restricted to these unexcluded regions, however, a definitive search of the charged Higgses over these parts of the m H + − tan β parameter space still has to be carried out by the LHC experiments.
In the following, for our numerical analysis we adopt the input parameter values from Ref. [32] as; G F = 1.
Note that the normalized decay rates in the ZM-VFNS (Γ pol NLO /Γ 0P ) are independent of the models while in the GM-VFN scheme the normalized widths of the polarized top decays (Γ pol NLO /Γ 0P ) depend on the model selected, extremely. Also, the results in the Type-I 2HDM scenario are independent of tan β, while the Type-II 2HDM results depend on the tan β.
Here, we are in a situation to present our results for the scaled-energy distribution of hadrons inclusively produced in polarized top decays in two variants of the 2HDM. Since the bottom quarks produced through the top decays hadronize before they decay and each b-jet contains a bottom-flavored hadron which most of the times is a B-meson, then we study the energy distribution of B-mesons. For this study, we consider the quantity dΓ(t(↑) → BH + + X)/dx B . According to the factorization formula (1), to evaluate dΓ/dx B one needs the parton-level decay width (dΓ pol GM /dx b ) described in section IV, and the nonperturbative fragmentation function D B b (z, µ F ) which describes the splitting of b → B at the desired scale µ F . To describe the hadronization process b → B, from Ref. [33] we adopt the nonperturbative fragmentation function D B b (z, µ F ) determined at NLO through a global fit to electron-positron annihilation data taken by OPAL [34], ALEPH [35] and SLD [36]. In [33], a simple power model as D b (z, µ ini F ) = N z α (1 − z) β is proposed as a initial condition for the b → B FF at the initial scale µ ini F = 4.5 GeV. Their fit results for the FF parameters read: N = 4684.1, α = 16.87, and β = 2.628. The nonperturbative FF D B b (z, µ F ) at each desired scale might be generated via the DGLAP evolution equations [37]. Note that, in the factorization formula (1) the factorization (µ F ) and renormalization (µ R ) scales are completely arbitrary and, in principle, one can select different values for them. In this work, we adopt µ R = µ F = m t that allows us to get rid of the term ln(µ 2 F /m 2 t ) in the the partial decay rate in the ZM-VFNs (dΓ pol NLO /dx b ), see the analytical result in [11]. In [11], we also investigated the dependence of the B-meson energy spectrum on these scales considering two other different values: µ R = µ F = m t /2 and µ R = µ F = 2m t . Since in Ref. [33] no uncertainty is reported for the nonperturbative FF D B b (z, µ F ), then this scale variation can be considered as a just theoretical uncertainty.
Considering the unexcluded MSSM m H + − tan β parameter space from the CMS [7] and ATLAS [8] experiments, in Fig. 1 we present our prediction for the x Bspectrum at NLO taking tan β = 8. Our results for the ZM-and GM-VFN schemes are compared in the model I, taking m H + = 95 GeV and m H + = 155 GeV. As is seen the zero-mass approximation (with m b = 0), applied in our previous work [11], works good to a high accuracy. Here, the B-hadron mass creates a threshold, e.g., at x B = 2m B /(m t (1+R−y)) ≈ 0.28 for m H + = 155 GeV when m b = 0.
In Fig. 2, considering the model I of 2HDM we study the energy spectrum of B-hadron in the ZM-and GM-VFN schemes for different values of tan β = 8 and 12, where the mass of Higgs boson is fixed to m H + = 155 GeV. As is seen, when tan β increases the decay  rate decreases in both schemes, asΓ 0P is proportional to cot 2 β, see (6) and (7). As in Fig. 1, the results of massless and massive schemes are, with a good approximation, the same so the results of ZM-VFN scheme show an enhancement in the size of decay rates, specifically at x B = 0.8. In Fig. 3, considering the GM-VFN scheme the NLO As is seen, the B-hadron energy spectrum in the 2HDM extremely depends on the model and one can not also naively set m b = 0 in all expressions in the model II. Also, the results from the model I are always larger than the ones from the model II in most of x B -regions. Note that, in one hand, the results of the model I do not depend on the b-quark mass largely ( Figs. 1 and 2) and, on the other hand, in the limit m b → 0 the results of both models are the same (see Eqs. (7) and (8)), then the results shown in Fig. 3 for the model I (solid and dotted lines), in principle, can be assumed as the results of model II in the ZM-VFN scheme. In fact, the results shown in Fig. 3 can be also considered as the results for the model II in the ZM-and GM-VFN schemes and this figure shows that there is an egregious difference between the massless and massive schemes. Finally, the results from Figs. 1-3, specially Fig. 3, show that the most reliable results are made in the GM-VFN scheme, specifically, when the Type-II 2HDM scenario is concerned.

VI. CONCLUSIONS
Charged Higgs bosons H ± are predicted in many extensions of the Standard Model consist of, at least, two Higgs doublets, of which the simplest are the two-Higgsdoublet models (2HDM) so the discovery of them would clearly indicate unambiguous evidence for the presence of new physics beyond the SM. Although, the charged Higgs bosons have been searched for in high energy experiments, in particular, at the Tevatron, ATLAS and CMS but they have not been seen so far and a definitive search is a program that still has to be carried out by the CERN LHC.
In this work, which is a fundamental extension of our previous work [11], we introduce a channel to indirect search for the charged Higgs bosons. In fact, since the main production mode of light charged Higgses in the 2HDM is through the top quark decay, t → bH + , and whence bottom quarks hadronize (b → B) before they decay, then the study of B-meson energy spectrum in the decay mode t → BH + +X would be of prime importance at the LHC. In other words, at the LHC any deviation of the B-meson energy spectrum from the SM predictions [15] can be considered as a signal for the existence of charged Higgs.
In [11], using the massless or ZM-VFN scheme (with m b = 0) we studied the spin-dependent energy distribution of B-mesons (d 2Γpol /(dx B d cos θ P )) produced through the polarized top decays in a special helicity coordinate system, where the event plane lies in the (x,ẑ) plane and the b-quark three-momentum is considered along theẑ-axis. In this system the polarization vector of the top quark is evaluated relative to the b-quark three-momentum. In the ZM-VFN scheme, since all information on the m b dependence of the B-hadron spectrum is wasted, then our results are reliable just for the Type-I 2HDM scenario. In the present work, we studied the same decay mode in the GM-VFN scheme where m b = 0 is considered from the beginning, however, considering the b-quark mass effects makes the calculations so complicated. Unlike the massless results, the massive decay rates are extremely dependent on the scenario selected in the 2HDM, specifically, when the type-II 2HDM scenario is concerned. Our results show that the most reliable predictions for the B-hadron energy spectrum are made in the GM-VFN scheme.
Our formalism elaborated in this work can be also extended to other hadrons, such as pions, kaons and protons, etc., using the nonperturbative b → π/K/P FFs extracted in [38,39], relying on their universality and scaling violations.

VII. ACKNOWLEDGMENTS
We would like to thank the LHC top working group for importance discussion and comments. We warmly acknowledge the CERN TH-PH division for its hospitality where a portion of this work was performed.