Indirect search for light charged Higgs bosons through the dominant semileptonic decays of top quark $t\to b(\to B/D+X)+H^+(\to \tau^+\nu_\tau)$

In this work we introduce a new channel to indirect search for the light charged Higgs bosons, which are predicted in several extensions of the standard model (SM) such as the two-Higgs-doublet models (2HDMs). We calculate the ${\cal O}(\alpha_s)$ QCD radiative corrections to the energy distribution of bottom- and charmed-flavored hadrons ($B/D$) produced in the dominant decays of the polarized top quark in the 2HDM, i.e. $t(\uparrow)\longrightarrow b(\to B/D+\text{jet})+H^+(\to \tau^+\nu_\tau)$. %This analysis is studied in a specific helicity coordinate system where the polarization vector of the top quark is evaluated with respect to the momentum direction of the bottom quark. Generally, the energy distribution of hadrons is governed by the unpolarized rate and the polar and the azimuthal correlation functions which are related to the density matrix elements of the decay $t(\uparrow)\rightarrow bH^+$. In our proposed channel, any deviation of the $B/D$-meson energy spectrum from its corresponding SM predictions can be considered as a signal for the existence of charged Higgs at the LHC. We also calculate, for the first time, the azimuthal correlation rate $\Gamma_\phi$ at next-to-leading order which vanishes at the Born level.


Introduction
Charged Higgs bosons emerge in the scalar sector of several standard model (SM) extensions, and are the object of various beyond SM searches at the CERN Large Hadron Collider (LHC). Since the SM does not include any elementary charged scalar particle, then the experimental observation of a charged Higgs boson would necessarily be a signal for a nontrivially extended scalar sector and a definitive evidence of new physics beyond SM. In recent years, searches for charged Higgs bosons have been done by the ATLAS and the CMS collaborations at the LHC in proton-proton collision and numerous attempts are still in progress.
Among many beyond SM scenarios which motivate the existence of charged Higgs, a generic two-Higgs-doublet model (2HDM) [1][2][3] provides a greater insight of the SUSY Higgs sector without including the plethora of new particles which SUSY predicts. Within this class of models, two isospin doublets are introduced to break the symmetry of SU (2) × U (1). This symmetry breaking leads to the existence of five physical Higgs bosons; three physical neutral Higgs bosons (h, H, A) and a pair of charged-Higgs bosons (H ± ) [2].
The dominant production and decay modes for a charged Higgs depend on the value of its mass with respect to the top-quark mass, and can be classified into three categories [4]. Light charged Higgs scenarios are defined by Higgs-boson masses smaller than the top quark mass. In the 2HDM, the main production mode for light charged Higgses is through the top quark decay t → bH + . Therefore at the LHC, as an formidable machine producing around 90 million top pairs per year of running at design c.m. energy of 14 TeV, the light Higgs bosons can be searched in the subsequent decay products of the top pairs tt → H ± H ∓ bb and tt → H ± W ∓ bb when H ± decays into τ lepton and neutrino.
On the other hand, in the decay mode t → bH + both b-quarks hadronize into the b-jet X b before they decay and the Higgs bosons decay into the leptons and their related neutrinos, i.e. H + → l + ν l (l = e, μ, τ ). Therefore, at the LHC the decay process t → b(→ X b + jet) + H + (→ l + ν l ) is of prime importance and it is an urgent task to predict its partial decay width as reliably as possible. Of particular interest are the determination of energy distribution of hadrons inclusively produced in the top quark rest frame. The study of these hadronic energy distributions in the polarized and unpolarized top decays could be proposed as a new approach to indirect search for the charged Higgs bosons at the LHC. Practically, any deviation of the energy spectrum of the produced hadrons in top quark decays from the SM predictions can be considered as a signal for the existence of charged Higgs. In [5], the energy distribution of bottom-flavored hadrons (B) inclusively produced in the SM decay chain of an unpolarized top quark, i.e. t → bW + → Bl + ν l + X, is studied and in [6] the energy distribution of B-hadrons is investigated in the unpolarized top decays in the 2HDM scenarios, i.e. t → bH + → BH + + X, where X refers to the unobserved final state particles.
Since, the life time of the top quark (≈ 4.6 × 10 −25 s) is much shorter than the typical time needed for the QCD interactions to randomize its spin, therefore, its full polarization content is preserved and passes on to its decay products, so that the polarization of the top quark will reveal itself in the angular decay distribution. Thus, the top quark polarization can be studied through the angular correlations between the direction of the top quark spin and the momenta of its decay products [7].
In this work, in the 2HDM framework we study the O(α s ) angular distribution of energy spectrum of B/D-mesons considering the polar and the azimuthal angular correlations in the rest frame decay of a polarized top quark, i.e. t (↑) → B/D + H + + X followed by H + → l + ν l . This angular correlation is analyzed in a specific helicity coordinate system where the event plane, including the top and its decay products, is defined in the (x, ẑ) plane and the b-quark three-momentum points into the direction of the positive ẑ-axis. In this frame, the top quark polarization vector is evaluated with respect to the direction of the b-quark three-momentum. To define the event plane, one needs to measure the momentum directions of the b-quark and the H + -boson and the polarization direction of the top quark. The evaluation of the b-quark momentum direction requires the use of a jet finding algorithm, while the top spin direction must be obtained from theoretical input, e.g. a polarized linear e + e − collider may be considered as a copious source of close to zero and close to 100% polarized tops.
The azimuthal correlation between the event plane and the intersecting one to this plane (including the top quark polarization vector) in the semileptonic rest frame decay of a polarized top quark (see Fig. 3) belongs to a class of polarization observables involving the top quark in which the leading-order contribution receives a zero result. As we will show, the nonzero contributions can arise from higher order radiative corrections. Due to the zero result for the lowest order contribution, the azimuthal decay rate up to NLO will be small. Then, it seems that the measurement of the azimuthal correlation would be difficult, but since highly polarized top quarks with more accuracy will become available at higher luminosity hadron colliders through single top production processes [8], it might then be feasible to experimentally measure this azimuthal correlation through the energy distribution of hadrons from polarized top decays. As was explained, the comparison of these polar and azimuthal distributions of meson energy spectrum produced in polarized top decays with the SM ones might be considered as a new channel to indirect search for the charged Higgs bosons at the LHC.
Concerning the importance and other applications of the analytical results presented in this work, note that, however current ATLAS and CMS measurements exclude a light charged Higgs for most of the parameter regions in the context of the minimal supersymmetric standard model (MSSM) scenarios, these bounds are significantly weakened in the Type II 2HDM (MSSM) once the exotic decay channel into a lighter neutral Higgs, H ± → AW ± /H W ± , is open. In [9], the possibility of a light charged Higgs produced in top decay via single top or top pair production is examined with a subsequent decay as H ± → AW ± /H W ± . It is shown that, this decay can reach a sizable branching fraction at low tan β once it is kinematically permitted. Their results show that the exotic decay channel H + → AW + /H W + is therefore complementary to the conventional H + → τ + ν τ channel, considered in the current MSSM scenarios. Therefore, our analytical results presented in this work will be also able to be applied for the decay process t (↑) → b(→ B/D + jet) + H + (→ AW + /H W + ) considering a convenient branching fraction for the decay H + → AW + /H W + . Note that, due to experimental challenges at low energies, such a light neutral Higgs has not been fully excluded yet.

Unpolarized top quark decay in the narrow width approximation
From the unitarity of the Cabibbo-Kobayashi-Maskawa (CKM) quark mixing matrix [10], one has the relation so, it results |V tb | ≈ 1 to a very high accuracy. Therefore, top quark decays within the SM are completely dominated by the mode t → bW + followed by W + → l + ν l (l = e, μ, τ ). In theories beyond SM with an extended Higgs sector, assuming m t > m b + m H + , one may also have the decay mode t → bH + followed by H + → l + ν l . Although, as an example, the τ + leptons arising from the decays W + → τ + ν τ and H + → τ + ν τ are predominantly left-and right-polarized, respectively. The polarization of the τ + influences the energy distributions of the decay products in the subsequent decays of the τ + → π +ν τ , ρ +ν τ , l + ν lντ . These distributions are not discussed in our work. For more detail, see Ref. [11].
Here, we first review some technical detail about the decay rate of unpolarized top quarks in the process . The ratio of VEVs is a free parameter and is defined as tan β = v 2 /v 1 . Even, a linear combination of the charged components of H 1 and H 2 gives two massive charged Higgs boson states H ± , i.e. H ± = H ± 2 cos β − H ± 1 sin β. In a general 2HDM, in order to avoid tree-level flavor-changing neutral currents that can be induced by Higgs-boson exchange, the generic Higgs coupling to all quarks should be restricted. Fortunately, there are several classes of 2HDMs which naturally avoid this difficulty by restricting the Higgs coupling. Imposing flavor conservation, there are four possible ways (hereinafter called models I-IV) for the two Higgs doublets to couple to the SM fermions. Each of these four ways gives rise to rather different phenomenology predictions. In these models, the generic charged Higgs coupling to fermions (assuming massless neutrinos) can be expressed as a superposition of rightand left-chiral coupling factors [12], so that the relevant part of the interaction Lagrangian of the process (2) is given by where A, B and C are model-dependent parameters. In the first model (model I), one doublet H 2 gives masses to all quarks and leptons and the other doublet H 1 essentially decouples from fermions. In this model, one has In a second model (model II), the doublet H 2 gives mass to the right-chiral up-type quarks (and possibly neutrinos) and another doublet H 1 gives mass to the right-chiral down-type quarks and charged leptons. In this model, the interaction Lagrangian (3) includes and in model IV, the roles of the two doublets are reversed with respect to the model II, so that These models are also known as type I-IV 2HDM scenarios. The type-II scenario is, in fact, the Higgs sector of the MSSM up to SUSY corrections [13,14]. In other words, in the MSSM one has a type-II 2HDM sector in addition to the supersymmetric particles, in particular the charginos, stops and gluinos. In this work to present our phenomenological predictions for the energy distribution of produced hadrons we restrict ourselves to the MSSM scenario, although we will discuss how our results can be generalized to other types. We start by discussing the Born term contribution to the decay rate of the process t → bl + ν l (l + = e + , μ + , τ + ) (2), where the top quark is presently considered as unpolarized. We consider the decay process where the four-momentum assignments are indicated in parentheses, see also Fig. 1. Here, to be specific, we concentrate on the case with l + ν l = τ + ν τ . Using the couplings from the Lagrangian (3) one can write the matrix element of the process (8) in the following form Note that, since the main contribution to the decay mode (8) comes from the kinematic region where H + boson is near its mass-shell, then one has to take into account its finite decay width H + . For this reason, in (9) we used the Breit-Wigner prescription of the Higgs propagator. In (9), the model dependent parameters A, B, C are defined in (4)- (7). From now on and for simplicity we introduce By this, the coupling of the charged Higgs to the bottom and top quarks is expressed as a superposition of scalar and pseudoscalar coupling factors. On squaring the Born matrix element (9) and taking the spin sums, one is led to the Born contribution as For the 3-body decay rate t → bτ + ν τ , one has where |M 0 | 2 = |M 0 | 2 /(1 + 2s t ) for which s t stands for the top quark spin. By working in the narrow-width approximation (NWA), the Breit-Wigner Resonance is replaced by a deltafunction as [15] 1 Using the NWA for the H + -boson, the three body decay t → bτ + ν τ is factorized as which is a result expected from physical intuition.
One considerable point about the interference terms which are ignored in our work: in this work, using the NWA we study the decay process t (↑) → bH + in the 2HDM scenario so in [16] we used the same approximation for the SM decay one, i.e. t (↑) → bW + . An important condition limiting the applicability of this approximation, however, is the requirement that there should be no interference of the contribution of the intermediate particle for which the NWA is applied with any other close-by resonance. Indeed, if the mass gap between two intermediate particles is smaller than one of their total widths, the interference term between the contributions from the two nearly mass-degenerate particles may become large. In other words, interference effects can be large if there are several resonant diagrams whose intermediate particles (with masses M 1 and M 2 for two resonances) are close in mass compared to their total decay widths: [15]. In these cases, a single-resonance approach or the incoherent sum of two resonance contributions does not necessarily hold. Then, if the mass difference is smaller than their total widths, the two resonances overlap. This can lead to a potentially large interference term, which is neglected in the standard NWA, but can be taken into account in the full calculation or in a generalized NWA [15]. In our calculations, the required condition to apply the NWA holds, i.e. |m H + − m W + | ( W + , H + ) if one sets W = 2.212 GeV (with m W = 80 GeV) and H + ≤ 0.08 (for m H + = 95 and 160 GeV which are applied in our work). Therefore, the overlap of two resonances can be ignored with high accuracy and two decay modes t → bH + → bτ + ν τ and t → bW + → bτ + ν τ can be studied separately, so the total decay width of top quarks can be simply obtained by the summation of both decay rates. More details about the interference effects in BSM processes with a generalized NWA can be found in [15].

Born level results for unpolarized top decays
In this section we calculate the Born level rate of the process (8) using the NWA where we put p 2 H + = m 2 H + from the beginning. According to Eq. (14), the Born width is expressed as To proceed, let us introduce the following notation for the kinematic variables The Born level amplitude for the process t → bH + is parametrized as M 0 =ū b (a + bγ 5 )u t , so for the amplitude squared one has therefore, the tree-level decay width reads where λ(x, y, z) = (x − y − z) 2 − 4yz is the triangle function and the coefficients "a" and "b" are defined in (10). The NLO QCD radiative corrections to the unpolarized rate (18) are given in our previous work [6]. Next, for the squared amplitude of the leptonic decay of charged Higgs one has: |M Born H + →τ + ν τ | 2 = 8c 2 p ν · p τ . This rate must be evaluated in the top quark rest frame. In the H + -rest frame (H r.f. ) shown in Fig. 1 where the angle θ is defined in the H + -boson rest frame. To obtain the lepton four-momentum in the top rest frame, the relevant Lorentz boost matrix reads where E H + = m t (1 − S) and | p H + | = m t Q. Therefore, the lepton four-momentum in the top rest frame reads Considering the notation defined in (16), the tree-level decay width for the leptonic sector of (8) is given by which is a dependent-model rate.
Since, current search strategies assume that the charged Higgs decays either leptonically (H + → τ + ν τ ) or hadronically (H + → cs), then following Ref. [17] we adopt the relevant branching fraction B H τ (H + → τ + ν τ ) (14) as where, in the model I (and IV) one has and for the model II (and III) one has These results are in complete agreement with Ref. [18].
In the limit of m 2 i /m 2 H → 0 (i = c, s), the branching ratio (22) in the type-I 2HDM reads which is independent of the tan β-values and in the type-II 2HDM (MSSM), one has Taking m c = 1.67 GeV, m τ = 1.776 GeV and |V cs | = 0.9734, the branching ratio in the model I is B H τ = 0.284. In the MSSM scenario, the branching ratio depends on the tan β. The dependence of this branching ratio on the tan β is plotted in Fig. 2. As is seen, for the numerical values of tan β > 5 that we apply in this work the branching ratio is B H τ = 1, to a very high accuracy. The model-independent bounds on the branching ratio of a light charged Higgs boson are transformed into limits in the m H ± − tan β parameter space. Direct searches at the LHC, with a center-of-mass energy of 7 TeV [19][20][21] and 8 TeV [22,23] set stringent constraints on the parameter space. We will discuss on these constraints when our numerical analysis is presented.

Polarized top quark decay in the narrow width approximation
Here, we concentrate on the light charged Higgs inclusively produced through polarized top quark decays t (↑) → b + H + (→ τ + ν τ ). Since, bottom quarks hadronize before they decay, therefore, the study of energy distribution of produced b-jets (bottom-or charmed-flavored mesons in this work) in the following process is proposed as a new channel to indirect search for the charged Higgs bosons at the LHC. Therefore, of particular interest are to evaluate the distribution in the scaled-energy (x B , x D ) of B/D-mesons in the top quark rest frame as realistically and reliably as possible. For this study, one needs to evaluate the quantity d /dx B (or d /dx D ) where, following the notation introduced in [5], the mesonic scaled-energy fraction is defined as where Considering the notation defined in (16), this scaledenergy is simplified as To evaluate the partial decay width of process (27) differential in x B (or correspondingly in x D ), we apply the factorization theorem of the QCD-improved parton model [24]. According to this theorem, the energy distribution of B-hadrons might be expressed as the convolution of the parton-level spectrum with the nonperturbative fragmentation function D B a=b,g (z, μ F ), where the D B g (z, μ F )-contribution is appeared at NLO. In (29), d /dx a is the parton-level differential width of the process (27). As in (28), we define the normalized energy fraction of partons (gluon and bottom quark) as In (29), μ F and μ R are the factorization and the renormalization scales, respectively, and generally one can use two different values for these scales. However, a choice often made consists of setting μ R = μ F and we use the convention μ R = μ F = m t in this work.
To achieve the spin-dependent energy distribution of hadrons, at first, we need to calculate the triply differential decay width d /(dx b d cos θ P dφ P ) of the partonic process t (↑) → bH + → b(τ + ν τ ). Polar and azimuthal angles θ P and φ P show the orientation of the plane including the spin of the top quark relative to the event plane, see Fig. 3.
Following the narrow width approximation (14), one has where in the MSSM, we assume B H τ (H + → τ + ν τ ) = 1 for tan β > 5 (Fig. 2). For the decay process t (↑) → bH + , to clarify the correlations between the top quark decay products and the spin of top quark, the general angular distribution of the triply differential decay width d /dx b is given by [25] d dx b d cos θ P dφ P = 1 4π where, P is the magnitude of the top quark polarization. P = 0 is for an unpolarized top quark while P = 1 stands for a 100% polarized top quark. In (31), d U /dx b corresponds to the unpolarized differential decay rate while d θ /dx b and d φ /dx b describe the polar and azimuthal correlations between the polarization of the top quark and its decay products, respectively. For the analysis of spin-momentum correlations between the top polarization vector and the momenta of its decay products, we consider the helicity frame shown in Fig. 3. In this frame, the three-momentum of b-quark is defined as p b (+ẑ), and the polarization vector of the top quark is evaluated with respect to the positive ẑ-axis.
Concerning the importance of this helicity frame, we take a moment to explain it in detail. Ignoring the interference terms, in the NWA the total decay width of an unpolarized top quark is given by tot = SM t→bW + + MSSM t→bH + . The same result is valid for the scaled-energy (x B ) distribution of the produced B-meson, i.e. d tot /dx B In other words, to obtain the total energy spectrum of B-mesons in the unpolarized top decays all decay modes including t → B + W + /H + should be summed up. In the presence of a polarized top quark, the same consequence is valid as long as the polarization of top quark is measured relative to the same three-momentum, for example, the momentum of bottom quark in our work. Since, the hadronization mechanism of bottom quarks is the same both in the SM and beyond SM theories, therefore, the study of energy spectrum of observed mesons in the introduced helicity frame can be considered as a new channel to indirect search for the charged Higgs bosons. In fact, any deviation of the total energy spectrum of produced meson from the corresponding SM predictions (presented in [16]) can be associated with the existence of charged Higgs bosons. Now, we back to Eq. (31). d U /dx b describes the unpolarized differential decay rate of the top quark which is independent of the selected frame. The O(α s ) radiative corrections to this rate are extensively studied in our previous work [6]. Here, we present our analytical results for the NLO radiative corrections to the angular correlation functions (d θ /dx b and d φ /dx b ) in the helicity frame shown in Fig. 3. Finally, at the hadron level we shall present and compare our predictions for the energy distribution of B/D-mesons, considering all contributions.

Born term results for the polarized top decay
It is straightforward to compute the Born term contribution to the partial decay rate of a polarized top quark in the 2HDM. The Born term amplitude squared for the process t (↑) → bH + reads where we replaced s t u(p t , s t )ū(p t , s t ) = (/ p t + m t ) in the unpolarized Dirac string by u(p t , s t )ū(p t , s t ) = (1 − γ 5 / s t )(/ p t + m t )/2 in the polarized state. This result is converted to the relation (17) if one sets s t = 0.
Considering Fig. 3, the top polarization four-vector is set as s t = P (0; sin θ P cos φ P , sin θ P sin φ P , cos θ P ), whereas one has p b · s t = −P | p b | cos θ P . Therefore, at the Born level the helicity structure of partial rate of the process t (↑) → bH + → bτ + ν τ reads where B H τ (H + → τ + ν τ ) ≈ 1 in the MSSM scenario. The unpolarized Born-level rate ˜ (0) U (t → bH + ) is given in (18), and for the polar and azimuthal correlation functions ˜ (0) θ and ˜ (0) φ , one has˜ In (34), the product of two coupling factors in the MSSM scenario reads Considering Eqs. (18) and (34), it is seen that the polarized and unpolarized decay widths have a minimum at tan β = √ m t /m b ≈ 6 if one takes m t = 172.9 GeV and m b = 4.78 GeV. Therefore, the t → BH + branching fraction has a pronounced dip at tan β ≈ 6. Although, this is partly compensated by a large value of the H + → τ + ν τ branching fraction, which is B H τ (H + → τ + ν τ ) ≈ 1 for tan β > 5, the product still has a significant dip at tan β ≈ 6 [17]. This state is not established for the type-I 2HDM scenario, as Eq. (25) shows that the branching ratio is independent of tan β.
The fact that ˜ (0) φ = 0 means that the azimuthal correlation measurement has zero analyzing power to analyze the polarization of the top quark and the nonzero contribution arises from the radiative corrections. real , so that the NLO decay amplitude squared reads Since, the contribution of azimuthal correlation to the Born amplitude is zero, then d˜ vir φ /dx b = 0. It means, the virtual one-loop corrections are contributed in the unpolarized rate (d˜ vir U /dx b ) and in the polar correlation function (d˜ vir θ /dx b ). Thus, the O(α s ) radiative corrections to the ˜ φ (and d˜ φ /dx b ) just result from the real gluon emissions. Due to the soft treatments of the real gluons emitted, the corresponding Feynman diagrams include the infrared (IR) singularities. Generally, to extract the ultraviolate (UV) singularities (which appear in the virtual corrections) and the IR divergences we work in D-dimensional regularization scheme where D = 4. In this scheme, as an example, the differential decay rate for the real corrections to the tbH + -vertex is given by where, dR 3 stands for the 3-body phase space element including the three-momenta of outgoing particles, i.e. p b , p g and p H + . Here M real is − 2p μ b + / p g γ μ 2p b · p g (a1 + bγ 5 )u(p t , s t ) μ (p g , s g ), where the polarization vector of the real gluon is denoted by (p g , s g ). The first and second terms in the curly brackets refer to the real gluon emissions from the top and bottom quarks, respectively.
In the dimensional regularization approach, the angular integral is written as More detail about this approach can be found in our previous work [16]. As a subtle issue of dealing with γ 5 = iγ 0 γ 1 γ 2 γ 3 in dimensional regularization (DR) scheme, note that, γ 5 is not well defined in D-dimensions and there is no unique way to handle γ 5 in DR. The anticommutation relation {γ 5 , γ μ } = 0 produces ambiguity so one can not simply apply this relation in general D-dimensions. In other words, the above anticommutation relation along with the relation T r[γ 5 γ α γ β γ γ γ δ ] = 0 can not be simultaneously satisfied in DR scheme. There are several prescriptions to prevent this ambiguity, see for example [26][27][28][29]. In this work, we employ the Breitenlohner-Maison-t'Hooft-Veltman scheme [29,30] in which we accept that γ 5 is a purely 4-dimensional object and therefore does not anitcommute with D-dimensional Dirac matrices, i.e. {γ 5 , γ μ } = 0, while the relation T r[γ 5 γ α γ β γ γ γ δ ] = 0 still holds in D-dimensions.
Since the lowest-order term contribution to ˜ φ vanishes (34), then one does not need to introduce any IR regularization scheme such as a fictitious gluon mass or dimensional regularization when calculating the azimuthal correlation. In other words, the tree-graph contributions to the ˜ φ are IR-finite at NLO.
To calculate the azimuthal differential rate d˜ φ /dx b , in (36) by working in four dimensions we fix the 3-momentum of the b-quark and integrate over the gluon energy, which ranges as where where C F = 4/3 stands for the color factor and α s (μ R ) is the strong coupling constant. Other parameters are defined in (16).
In the massless scheme where the m b = 0 approximation is considered, this result simplifies as ˆ φ =ˆ (0) θ C F α s /2. Since the observed mesons in top decays can be also produced through a fragmenting real gluon, therefore, to obtain the most accurate energy spectrum of mesons one has to add the contribution of gluon fragmentation to the b-quark one. Then, one also needs the azimuthal partial decay rate dˆ φ /dx g , where x g = E g /E max b is the scaled-energy fraction of the real gluon, as in (28). In [5,6], we showed that the gluon contribution would be important at a low energy of the detected meson so this contribution decreases the size of decay rate at the threshold energy. In the helicity frame selected, the azimuthal differential decay width dˆ φ /dx g is given by In [5], we showed that the g → B contribution into the NLO energy spectrum of B-meson is negative and appreciable only in the low-x B region and for higher values of x B the NLO result is practically exhausted by the b → B contribution. The contribution of gluon is calculated to see where it contributes to d /dx B and can not be discriminated as an experimental quantity. The NLO expression for the d˜ θ /dx b is obtained by summing the Born term, the virtual one-loop and the real gluon contributions. To extract all IR-and UV-singularities one needs to work in D-dimensions considering the relation (38).
The QCD virtual one-loop corrections arise from emission and absorption of a virtual gluon from the same quark leg (quark self-energy) and from a virtual gluon exchanged between the top and bottom quark legs (vertex correction). The self-energy Feynman diagrams are related to the mass and wave function renormalizations of both top and bottom quarks, more detail can be found in [31]. Adopting the on-shell mass-renormalization scheme, the virtual one-loop corrections include both the IRand UV-divergences so the UV ones appear when the integration region of the internal momentum of the virtual gluon goes to infinity and the IR-divergences arise from the soft-gluon singularities. The UV-divergences are canceled after summing all virtual corrections up but the IR-singularities are remaining. To remove the remaining IR-divergences, one needs to consider the real gluon corrections at NLO.
The NLO real contributions result from the real gluon emissions from the bottom and top quarks, individually. In calculation of the real gluon corrections, since we preserve the mass of b-quark from the beginning, there are no collinear singularities and all IR-divergences arise from the soft real gluon emission. According to the Lee-Nauenberg theorem, all singularities cancel each other after summing all contributions up, and the final result is free of IR singularities. More detail about the NLO corrections and the D-dimensional regularization scheme can be found in [6] where we calculated the unpolarized differential decay rate d˜ U /dx i (i = b, g). The technical detail of calculations are the same.
Ignoring the detail of our calculations, the NLO expression for the d˜ θ /dx b is as follows where the plus distribution is defined as usual, and Also, the gluon contribution is expressed as The analytical results presented in (40)- (45) are completely new and are being presented for the first time. Specifically, the contribution of azimuthal correlation φ has not been already calculated. One needs these angular differential decay rates to obtain the energy spectrum of produced mesons in polarized top decays, see (29).

GM-VFN scheme
As was mentioned, our main aim is to evaluate the scaled-energy distribution of the B/D-mesons produced in the inclusive process t (↑) → H + B/D + X followed by H + → τ + ν τ , where X stands for the unobserved final state. Therefore, we calculate the NLO decay width of the corresponding process differential in x B (d /dx B ) and x D (d /dx D ) in the general-mass variable-flavor-number scheme (GM-VFNs), where x B and x D are defined in (28). In the top quark rest frame applied in our work, the B-meson has the energy E B = p t · p B /m t , where In the case of gluon fragmentation, g → B, it has energy [5]. The same results hold for the D-meson with the mass m D .
Considering the factorization theorem (29), the energy spectrum of B/D-meson can be obtained by the convolution of the convenient parton-level spectrum with the nonperturbative fragmentation function D B i=b,g (z, μ F ) (or D D i=b,g (z, μ F )). We will discuss about these FFs in section 5. Now, we explain how the quantity d θ (μ R , μ F )/dx i (with i = b, g) will have to be evaluated in the GM-VFN scheme. In Sec. 3.2, we employed the Fixed-Flavor-Number (FFN) scheme which contains 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 resume 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˜ θ (μ R , μ F )/dx i . In this case, the NLO ZM-VFN results are exactly recovered in the limit m b /m t → 0. In the GM-VFN scheme, the subtraction terms are constructed as (43) and (45) and 1/ˆ (0) θ × dˆ ZM θ /dx i are the partial decay rates computed in the ZM-VFN scheme [7], in which all information on the m b -dependence is wasted.
In conclusion, the GM-VFN results are obtained by subtracting the subtraction terms from the FFN ones [32,33], i.e.
Taking the limit m b → 0 in Eqs. (43) and (45), one obtains the following subtraction terms and The result (48) coincides with the perturbative FF of the transition b → b [34][35][36][37][38][39][40][41] and is in consistency with the Collin's factorization theorem [24], which guarantees that the subtraction terms are universal. Thus, the result obtained in (48) ensures the correctness of our result shown in (43). In [5,42], the GM-VFN scheme is applied to study the NLO energy spectrum of B-meson in the process t → bW + . There was shown, Eq. (48) coincides with the perturbative FF of the transition b → b as well.

Numerical analysis
In the MSSM, the charged Higgs mass is restricted by m H ± > m W ± at tree-level [43], however, this restriction does not hold for some regions of the MSSM tan β − m H ± parameter space after including radiative corrections. In Ref. [11], it is mentioned that a charged Higgs having a mass in the range 80 GeV ≤ m H ± ≤ 160 GeV is a logical possibility and its effects should be searched for in the decay chain t → bH + (→ τ + ν τ ). Searches for charged Higgs bosons have already been started at the Tevatron and, at the present, the last results in 19.5-19.7 fb −1 of proton-proton collision data recorded at √ s = 8 TeV are reported by the CMS [44] and the ATLAS [45] collaborations, using the τ + jets channel with a hadronically decaying τ lepton in the final state. According to these results, the large region in the (m H + − tan β)-plane 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). These regions along with the ±1σ band around the expected limit are shown in Fig. 4 which is taken from Ref. [45]. A same exclusion is reported by the CMS collaboration [44]. In this work, we restrict ourselves to these unexcluded regions. Although, a definitive search over these remaining parts of the m H + − tan β parameter space still has to be carried out by the LHC experiments.
Having the differential decay rates at the parton level, in the first step we turn to our numerical predictions of polarized and unpolarized decay rates, while the strong coupling constant is evolved from α s (m Z ) = 0.1184 to α s (m t ) = 0.1070. Considering our analytical results for the unpolarized decay width U [6], the azimuthal correlation function ˜ φ (41) and the polar correlation rate ˜ θ (43), and by taking m H + = 160 GeV and tan β = 16 from Fig. 4 where = 125 × 10 −6 GeV and B H τ (H + → τ + ν τ ) = 0.284. The corresponding result in the type-II 2HDM (MSSM) scenario reads where = 106 × 10 −4 GeV and B H τ ≈ 1. Taking m H + = 95 GeV and tan β = 8 from Fig. 4, in the type-I 2HDM scenario, one has where = 139 × 10 −4 GeV and B H τ = 0.284, as before. The corresponding result in the type-II 2HDM scenario is where = 830 × 10 −4 GeV and B H τ ≈ 1. As is seen, the NLO radiative corrections to the polarized and unpolarized rates are significant, specifically, when the type-II 2HDM scenario is concerned. Also, the azimuthal decay rates are quite small in both models, as is expected, so it seems that their measurements are difficult. Since highly polarized top quarks with more accuracy will become available at higher luminosity hadron colliders through single top production processes [8], then it might then be feasible to experimentally measure these azimuthal correlations.
Apart from the above results, in the following we present our phenomenological predictions for the scaled-energy (x B ) distribution of bottom-flavored hadrons and also the x D -distribution of at the desired scale μ F . Considering the description presented in Sec. 3.2 about the importance of the gluon splitting into the hadrons, one also needs the g → B/D FFs. These nonperturbative FFs are universal and process independent. Several models have been yet proposed to describe these hadronization processes. In [47], the FFs of D 0 -and D + -mesons are determined at NLO by fitting the experimental data from the BELLE, CLEO, OPAL, and ALEPH collaborations in the modified minimal-subtraction (MS) factorization scheme. There, authors have parametrized the FFs of b → D 0 /D + at the starting scale μ 0 = m b , as while the FF of gluon is set to zero. This parametrization is known as Bowler model. These FFs are evolved to higher scales using the DGLAP evaluation equations [48]. The values of fit parameters together with the achieved values of χ 2 are presented in Table 1.
For the b → B hadronization process, from Ref. [49] we employ the FFs determined at NLO through a global fit to e + e − -annihilation data taken by ALEPH and OPAL at CERN LEP1 and by SLD at SLAC SLC. In [49], a power model as D B b (z, μ 0 ) = Nz α (1 − z) γ is used at the initial scale μ 0 = 4.5 GeV, while the gluon FF is generated via the DGLAP evolution. The fit yielded N = 4684.1, α = 16.87, and γ = 2.628 with χ 2 = 1.495.
Considering the unexcluded MSSM m H + − tan β parameter space shown in Fig. 4, in Fig. 5 we present the x B distribution of d /dx B considering the azimuthal (dotted line) and the polar (solid line) correlation contributions. Here, we set tan β = 8 and m H = 95 GeV. For a more quantitative comparison, the unpolarized contribution (dashed line) is also shown. Here, the Bhadron mass creates a threshold at x B = m B /(m t S) ≈ 0.087. As was explained in section 3.1, the azimuthal correlation rate is zero at LO, which explains the smallness of the corresponding result at NLO.
In Fig. 6, the same analysis is presented but taking m H = 150 GeV; a charged Higgs mass which is not excluded in the MSSM m H + − tan β parameter space. Here, the B-meson mass creates a threshold at x B ≈ 0.22. From these two plots it is seen that in the unpolarized top decay the partial decay width at the hadron level is always higher than the one in the polarized top decay, specially, in the peak region.    (40). Our formalism elaborated in this manuscript can be also extended to the production of hadron species other than bottom-and charmed-flavored hadrons, such as pions, kaons and protons, etc., using the nonperturbative (b, g) → π/K/P FFs extracted in [50], relying on their universality and scaling violations.

Conclusions
In a general 2HDM, the main production mode of light charged Higgs bosons is through the top quark decay, t → bH + followed by H + → τ + ν τ . On the other hand, bottom quarks hadronize before they decay, therefore, the study of energy distribution of hadrons inclusively produced in top decays is of prime importance at the LHC. This study is proposed as a new channel to indirect search for the light charged Higgs bosons.
In this work, by working in the NWA framework we studied the O(α s ) spin-dependent energy spectrum of the charmed-and bottom-flavored hadrons produced through polarized top quark decays; t (↑) → bH + → B/D(τ + ν τ ) + X. The outgoing mesons are identified by a displaced decay vertex associated which charged lepton tracks.
In order to make our predictions for the energy spectrum of observed mesons we, for the first time, calculated the analytic expressions for the O(α s ) corrections to the angular correlation functions d˜ θ /dx i and d˜ φ /dx i (i = b, g) in a specific 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 frame, the polarization vector of top quark is evaluated relative to the positive ẑ-axis. These quantities are required to calculate the spin-dependent differential width d /(dx B d cos θ P dφ P ) of the process t (↑) → B/D(τ + ν τ ) + X. Since, highly polarized top quarks will become available at hadron colliders through single top production processes, which occur at the 33% level of the tt pair production rate [8], and in top quark pairs produced in future linear e + e − -colliders [51] these studies can be considered as an indirect probe to search for the charged Higss bosons. In fact, any deviation of the B/D-meson energy spectrum from the SM predictions [16] can be considered as a signal for the existence of charged Higgs at the LHC.
Our analytical results can be also applied for the exotic decay channel t (↑) → b(→ B + jet) + H + (→ AW + /H W + ), where A and H are the physical neutral Higgs bosons predicted in the 2HDM. This channel can reach a sizable branching fraction at low tan β [9]. Due to experimental challenges at low energies, such a light neutral Higgs has not been fully excluded yet.
As experimental considerations, two following points are discussed: 1)-In the decay modes t → bH + /W + followed by H + /W + → τ + ν τ , the tau leptons arising from the decays W + → τ + ν τ and H + → τ + ν τ are predominantly left-and rightpolarized, respectively. The polarization of the τ + influences the energy distributions of the decay products in the subsequent decays of the τ + → π +ν τ , ρ +ν τ , l + ν lντ and, in conclusion, it influences the energy distribution of mesons produced through the b-quark hadronization. Therefore, the energy spectrum of (for example) B-mesons inclusively produced in the decay chain t → b(W + , H + ) → B(τ + ν τ ) + X followed by τ + → π +ν τ , ρ +ν τ , l + ν lντ can also help in searching for the induced charged-Higgs effects at the Tevatron and the LHC. Strategies to enhance the H ± -induced effects in the decay t → bW + → b(τ + ν τ ), based on the polarization of the τ + have been discusses at length in Ref. [13] and references therein.
2)-Concerning the possible backgrounds on the proposed channels to search for charged Higgses (tt → W ± H ∓ bb and tt → H ± H ∓ bb) it should be noted that the main background for light charged Higgs events is the SM top pair production process (tt → W ± W ∓ bb). Since, both signal and background are produced from the same hard event, distinguishing beyond SM (BSM) top quark decay needs careful selection criteria; the most important of which originates from the fact that the charged Higgs decay products acquire predominantly higher four-momenta due to the higher charged Higgs mass compared to that of W-boson in the SM top decay t → bW + . Comparing this decay with the BSM one, t → bH + , and by assuming a final decay to the tau lepton for both W-and H-bosons, the tau lepton in SM background is produced with softer kinematic distribution than the one in signal events. The produced hadrons in top decays also have different kinematic distributions. This feature and the differences between the tau lepton kinematic distributions are clues for signal identification for further studies like the one presented in this paper. The overall selection strategy should be similar to what has been used at LHC charged Higgs observation and therefore follows the LHC analyses aiming at light charged Higgs observation at the currently open part of the m H ± − tan β parameters space. By increasing the data set and tan β the signal receives more purity and resolution and in conclusion the error of measurements would be less.