text:
# The Effects of Mergers on the Formation of Disc-Bulge Systems in Hierarchical Clustering Scenarios

## 1 Introduction

The origin of the Hubble sequence is still a controversial issue. In particular, spiral galaxies seem to evolve along this phenomenological classification although the physical mechanisms behind these morphological changes are not fully understood. The structural parameters, the gas abundances and the star formation activity vary along galaxies in the Hubble sequence in the sense that disc-dominated systems are also the more gaseous ones, experiencing on-going star formation. The large database of galaxies gathered in the last years have allowed to get more detailed information on the properties of different morphological type galaxies in the local Universe.

Observational results show that when a double exponential decomposition is applied to the sufarce luminosity density of late-type spirals, the scalelengths of the bulge and disc components seem to be restricted to a certain value $`<r_\mathrm{b}/r_\mathrm{d}>0.10`$ (Courteau, de Jong & Broeils 1996, hereafter CdJB96). This restriction in the range of possible scalelengths for spirals has been interpreted as a proof for secular evolution to be the responsible mechanism for bulge formation from an already in place disc structure. Recently, MacArthur, Courteau & Holtzman (2002, hereafter MCH02) studied a larger sample of late-type spirals and carried out a bulge-disc decomposition assuming a Sérsic law for the bulge surface luminosity density. These authors found a continuous distribution of the shape parameter of the bulge ($`n0.22`$) with a maximum at $`n0.90`$. They also found a restricted range of values for the scalelengths of the bulge and disc when the shape parameter is $`n1`$, $`<r_\mathrm{b}/r_\mathrm{d}>=0.13\pm 0.06`$ in the $`R`$-band. Previous works have also found that bulges of spirals can be fitted by a Sérsic law with different shape parameters (e.g., Andredakis, Peletier & Ballcels 1995, hereafter APB95; Khosroshahi, Wadadekar & Kembhavi 2000, hereafter KWK00; Graham & de Block 2001, hereafter GdB01). APB95 and GdB01 also found a correlation between the luminosity bulge-to-disc ratio $`B/D`$ and the shape parameter of the bulge, which also correlates with morphological type. These observational results suggest a possible connection between the formation of the bulge and disc components. They also support the idea that $`n`$ could be used as a good indicator of position in the Hubble sequence.

Several theories have been developed to explain the formation of the disc and bulge components. In the case of the disc systems, the standard model is based on three hypothesis: the angular momentum is acquired through cosmological torques (Peebles 1969), baryons and dark matter have the same specific angular momentum content ($`J`$), and this specific angular momentum is conserved during the collapse and cooling of baryons (Fall & Efstathiou 1980, hereafter FE80). This model has been successful in reproducing several observational results in analytical and semi-analytical models (e.g., Dalcanton, Spergel & Summers 1997; Mo, Mao & White 1998). However, serious problems arose in numerical simulations where an important angular momentum transfer from baryons to the dark matter haloes during mergers was detected, breaking the condition of $`J`$ conservation. Domínguez-Tenreiro, Tissera & Sáiz (1998, hereafter DT98) showed that a disc-like structure with observational counterpart (Sáiz et al. 2001, hereafter S01) can be built up if a compact stellar bulge is allowed to form without depleting the gas reservoir of the system. These stellar bulges provide stability to the gaseous disc systems which are capable of conserving a non-negligible fraction of its angular momentum during violent events. Supernova energy feedback could also contribute to the formation of the disc component, regulating the star formation rate and preventing early catastrophic depletion of the gas into stars. Probably both mechanisms, the formation of a compact stellar bulge which assures the axisymmetrical character of the potential well and energy feedback, work together in nature to allow the formation of spiral galaxies (e.g. Weil, Eke & Efstathiou 1998).

The formation of bulges is a more complex task since several mechanisms could be acting together such as monolithic collapse (Gilmore & Wyse 1998), mergers (Kauffman, Guiderdoni & White 1994) and secular evolution (Pfenniger & Norman 1990). In general, analytical models and pre-prepared simulations have focused in one or two of them at the time (e.g., van den Bosch 1999; Aguerri, Balcells & Peletier 2001, hereafter A01). These models have been successful in explaining some properties of spiral galaxies, despite their approximations.

The current paradigm for the formation of the structure favors a hierarchical clustering scenario where the structure forms by aggregation of substructure. Hence, a galactic object experiences the effects of collapse, merger, interaction and probably secular evolution in a non simple fashion, which can also strongly depend on redshift. In particular, violent events (i.e. mergers, interactions) can have important effects on the internal properties of the objects, such as their mass distribution (Mihos & Hernquist 1996) and star formation activity (Tissera et al. 2002, hereafter T02).

Minor and major mergers of disc systems with satellites have been extensively studied. However, most of these works do not consider the presence of the stellar bulges. The relevance of this component in the overall evolution of a disc system during a merger has been firstly pointed out by Mihos & Hernquist (1994) in a study of pre-prepared mergers. Recently, A01 studied the effects of mergers on the structural parameters of bulges of disc-like systems by assuming that initially bulges form with an exponential profile. These authors found that collisionless mergers produce a migration to higher $`n`$ and that this effect is proportional to the relative mass of the colliding systems. The simulations used in this analysis were non-cosmological pre-prepared ones where on going star formation and gas dynamics were not included. Recently, T02 show for the first time in cosmological SPH simulations how the central mass concentration in disc-like systems can grow by collapse, merger and secular evolution.

These numerical results suggest that a key point in the formation of both bulges and discs in hierarchical clustering scenarios seems to be mergers. Mergers have been traditionally suggested as a possible mechanism to drive the Hubble sequence since they can modify the dynamical and astrophysical properties of galactic objects. Hence, the question would be how, in hierarchical scenarios, mergers work to shape disc-like objects and if the properties of these objects resemble those of current spirals. In this paper we will focus on the analysis of the mass distributions of disc-like systems and how they are modified during mergers, paying special attention to the comparative study of the effects of secular evolution triggered by tidal fields and of the actual collision of the baryonic clumps. In order to assess the effects of mergers on the mass distributions we perform a bulge-disc decomposition of the projected surface density of the systems at different stages of the merger process. We then study how such different structural parameters, including the shape parameter defined by Sérsic (1968), evolve during the orbital decay phase and fusion of the satellite.

Our simulations include the effects of gravitation, hydrodynamics, cooling and star formation in a cosmological framework. Hence, the set of mergers that we analyse are given by the particular evolutionary history of each galactic object. This is a crucial point since the merger parameters and physical characteristics of the colliding objects are not ad-hoc choices but result from the consistent formation of the structure in a hierarchical scenario. The drawback of our approach is a lower numerical resolution compared to those used in studies of pre-prepared mergers. Assessment of possible numerical problems are discussed throughout the paper.

In section 2 we present the analysis of the simulations. In section 3 we discuss the results. Section 4 summarizes the conclusions.

## 2 Models and Analysis

We run cosmological SPH simulations which include hydrodynamics and star formation (Tissera, Lambas & Abadi 1997). The simulated boxes represent typical regions of a standard Cold Dark Matter (CDM) universe of 5 $`h^1`$ Mpc side and $`64^3`$ particles ($`\mathrm{\Omega }=1`$, $`\mathrm{\Lambda }=0`$, $`h=0.5`$). We assume a baryonic density parameter of $`\mathrm{\Omega }_\mathrm{b}=0.10`$. All baryonic and dark matter particles have the same mass, $`M_\mathrm{p}=2.6\times 10^8\mathrm{M}_{}`$. We have used a gravitational softening of $`3.0`$ kpc and a minimum hydrodynamical smoothing length of $`1.5`$ kpc. We performed three simulations with different realizations of the power spectrum (cluster normalized) and the same cosmological and star formation parameters. The simulated volume is the result of a compromise between the need to have a well-represented galaxy sample and enough numerical resolution to study the astrophysical properties of the simulated galaxies. We are confident that since we focus our analysis on small scale processes such as mergers and interactions, scale fluctuations of the order of the box size will have no significant effect on such local processes. Note also that we are interested in the effects that mergers have on the mass distributions of the galactic objects analysing them as individual events and not in connection with their evolution or environment.

The star formation algorithm used in these models is based on the Schmidt law (see Tissera et al. (1997) for details). Cold gas particles are eligible to form stars if they are denser than a certain critical value and satisfy the Jeans instability criterium (Navarro & White 1994). Gas particles are checked to satisfy these conditions at all time-steps of integration. Hence, as the gas cools down and is gathered at the centre of dark matter haloes, it is gradually transformed into stars according to the particular history of evolution of each galactic object. Only one free parameter has to be fixed: the star formation efficiency, which links the dynamical time of the gas cloud represented by a particle with its star formation timescale. We have used the value adopted by S01 since it is adequate to reproduce disc-like structures with observational counterparts at $`z=0`$. These authors used a low star formation efficiency which allows the formation of a stellar bulge that assures the axisymmetrical character of the potential well but without exhausting the gas reservoir (see also DT98). As a consequence, disc-like structures can be formed, although they remain mainly gaseous. Conversely, observed discs are mostly stellar, but they inherit the structural and dynamical properties of the gaseous discs out of which disc stars are formed. Hence, for the sake of comparison, the values of the parameters describing those properties can be reliably estimated from the simulated gaseous discs. In this paper, we will determine and use these parameters to study the effects of mergers on the mass distribution of the simulated galactic objects.

At $`z=0`$, we identify galactic objects at their virial radius, analysing only those with more than 4000 dark matter particles within their virial radius in order to diminish numerical resolution problems. Within each galactic object, a main baryonic clump is individualized which will be, hereafter, called the galaxy-like object (GLO).

The selected GLOs have very well-resolved dark matter haloes which provide adequate potential wells for baryons to collapse in. This fact assures a reliable description of the gas density profiles (see Steinmetz & White 1997), which allows us to follow the star formation history of the GLOs (see Tissera 2000). With this strong restriction on the minimum number of particles, the final GLO sample at $`z=0`$ is made up of 12 GLOs with virial velocities in the range 140-180 km$``$s-1.

We follow the evolution of the selected GLOs with look-back time ($`\tau (z)=1(1+z)^{3/2}`$) identifying mergers and starbursts (SBs). We then construct their mergers trees and star formation rate (SFR) histories (because of the restrictive SF efficiency used, these SFR histories are mainly those of the bulge components). During a merger event, the progenitor object is chosen as the more massive baryonic clump within this merger tree, while the minor colliding baryonic clump is referred to as the satellite. Following T02 we study only those mergers that are directly linked to stabursts in the SF history of each GLO. These authors found that during the orbital decay phase of some merger events, early gas inflows which trigger star formation can develop. In this case, if the gas reservoir of the GLOs is not exhausted during these starbursts, second bursts are induced at the actual physical contact of the baryonic clumps. The induction of early gas inflows could be directly linked to the properties of the total potential wells of the systems in agreement with previous results reported by Barnes & Hernquist (1996) and Mihos & Hernquist (1996) in pre-prepared simulations. When no gas inflows are triggered during the orbital decay phase (ODP), only one SB is detected when the two baryonic clumps collide. These events will be called single SBs (SSBs), while those where two SBs are detected will be called double SBs (DSBs). We studied a total of 18 merger events with merger parameters settled by the cosmological model. Among them, 11 are classified as DSBs and 7 as SSBs. Fig. 1 shows an example of both types of SBs. We have also plotted the distance between the centres of mass of the progenitor and the satellite since the redshift they share the same dark matter halo. The merger event will be determined by four redshifts of reference, $`z_\mathrm{A}`$: the beginning of the first bursts in DSBs and the ODP in SSBs, $`z_\mathrm{B}`$: the end of both the first bursts in DSBs and the ODP in SSBs, $`z_\mathrm{C}`$: the beginning of the second SBs in DSBs and the only bursts in SSBs (in this case $`z_\mathrm{B}=z_\mathrm{C}`$), and $`z_\mathrm{D}`$: the end of the second bursts in DSBs and of the only one in SSBs. These redshifts have been depicted in Fig. 1. Note that the ODP, defined as the time since two objects share the same dark matter halo to the fusion of their main baryonic cores, is determined by $`z_\mathrm{A}`$ and $`z_\mathrm{C}`$. The $`z_\mathrm{B}`$ redshift is used in the case of DSBs since it establishes the end of the first SBs.

In order to quantify the changes in the baryonic mass distributions, we carried out a bulge-disc decomposition of the projected mass surface density of the GLOs at $`z=0`$ and of the progenitors at each of the redshifts that define a merger event. We assumed an exponential profile for the disc component and the Sérsic law for the bulge central mass concentration:

$$\mathrm{\Sigma }_\mathrm{d}=\mathrm{\Sigma }_\mathrm{d}^\mathrm{o}\mathrm{exp}((r/r_\mathrm{d}))$$
(1)
$$\mathrm{\Sigma }_\mathrm{b}=\mathrm{\Sigma }_\mathrm{b}^\mathrm{o}\mathrm{exp}((r/r_\mathrm{b})^{1/n})$$
(2)
where $`r_\mathrm{d}`$ and $`r_\mathrm{b}`$ are the disc and bulge scalelengths, $`\mathrm{\Sigma }_\mathrm{d}^\mathrm{o}`$ and $`\mathrm{\Sigma }_\mathrm{b}^\mathrm{o}`$ the corresponding central mass surface densities and $`n`$ the bulge shape parameter.

For each GLO at $`z=0`$ and its progenitor objects, we estimated the total angular momentum ($`j`$) and projected the mass distribution onto the perpendicular plane defined by $`j`$. We then integrated the projected baryonic mass in concentric cylinders of radius $`r`$ since, as it has been proved by S01, it is numerically more convenient to fit the integrated mass than the surface density. In order to diminish the effects of numerical resolution, the fits were made from $`r=1.5`$ kpc to $`r30`$ kpc, except in those cases where the mass distributions of the progenitors were strongly perturbed by the entrance of the satellites.

The total integrated baryonic mass of the GLOs ($`M_{\mathrm{GLO}}`$) is estimated, on one hand, as the baryonic mass within $`r=30`$ kpc, and will be associated with the luminous galaxy. On the other hand, $`M_{\mathrm{GLO}}`$ is the result of the disc ($`M_d`$) and bulge ($`M_b`$) mass contributions:

$$M_{\mathrm{GLO}}=M_\mathrm{d}+M_\mathrm{b}$$
(3)
$$M_\mathrm{d}=2\pi \mathrm{\Sigma }_\mathrm{d}^\mathrm{o}r_\mathrm{d}^2$$
(4)
$$M_\mathrm{b}=2\pi \mathrm{\Sigma }_\mathrm{b}^\mathrm{o}r_\mathrm{b}^2n\mathrm{\Gamma }(2n)$$
(5)
where $`\mathrm{\Gamma }`$ is the complete gamma function. Since $`M_{\mathrm{GLO}}`$ is known, there are four free parameters to be fitted: $`r_\mathrm{d}`$, $`r_\mathrm{b}`$, $`\mathrm{\Sigma }_\mathrm{b}^\mathrm{o}`$ (or $`\mathrm{\Sigma }_\mathrm{d}^\mathrm{o}`$) and $`n`$. From these fitted parameters we analyse possible correlations with the dynamical and astrophysical properties of the GLOs and their dark matter haloes. In particular, we will study how the structural parameters change during merger events.

Note that we work with masses instead of luminosities, hence, mass-to-light ratios for the bulge and disc components will have to be adopted in order to be able to compare our results with observations. To this end we use different observational results from APB95 (30 S0-Sbc galaxies in the $`K`$-band), CdJB96 (173 spiral Sb-Sc galaxies in deep $`R`$-band), Moriondo, Giovanelli & Haynes (1999, hereafter MGH99, 25 Sa-Scd spiral galaxies in the $`K`$-band), KWK00 (27 early-type disc galaxies in the $`K`$-band) and MCH02 (47 late-type spiral galaxies in the $`R`$-band). APB95, KWK00 and MCH02 used bulge-disc decompositions of the luminosity surface density profiles, allowing the bulge shape parameter to vary. CdJB96 used double exponential fittings ($`n=1`$). Finally, MGH99 made fittings with different $`n`$ values finding that most of their bulges were better described by $`n=1`$ (except for two with $`n=2`$). Note that all these observations correspond to galaxies at $`z0`$.

## 3 Results

### 3.1 GLOs at $`z=0`$

As we mentioned before, S01 showed that disc-like objects formed in numerical simulations where stellar bulges were allowed to form have structural parameters and specific angular momentum content similar to those of observed spiral galaxies. For the sake of completeness and to show that these findings are valid for our GLOs, we will first discuss their properties at $`z=0`$. The structural parameters and the circular velocity ($`V^2=GM(r)/r`$) at $`r=2.2r_\mathrm{d}`$ allow the classification of our objects as intermediate spirals (100 km$``$s$`{}_{}{}^{1}V_{2.2}`$ 180 km$``$s-1; 4 GLOs), bright spirals ($`V_{2.2}>`$ 180 km$``$s-1 and $`r_\mathrm{d}>5.25`$ kpc; 2 GLOs) and compact bright spirals ($`V_{2.2}>`$ 180 km$``$s-1 and $`r_\mathrm{d}<5.25`$ kpc; 6 GLOs). This is the result of our strong condition on the mass of GLOs in order to work with those better resolved.

In Fig. 2 we show the relation between the disc and bulge scalelengths. Each GLO has a label that specifies its $`n`$ shape parameter. We include the observational data from CdJB96, MGH99, KWK00 and MCH02. Note that for increasing $`n`$ parameters, the scalelengths tend to move to the upper envelope of the distribution. This effect is also seen in the observed scalelength distribution. We have adopted $`n=1.5`$ as a general cut-off value to segregate the samples between those with exponential and non-exponential profiles for the bulge components.

In Fig. 3 we show the normalized histrogram of $`n`$ for GLOs at $`z=0`$ (a), and the corresponding distribution for the observations of APB95, KWK00 and MCH02 (b), normalized to the total number of galaxies in each sample. The combined observational sample covers a morphological range from Sa to Sc galaxies, but this is not a consistent sample. We show it in order to broadly compare the range of $`n`$ values obtained from the simulations with observations. The average $`n`$ values over the observed subsamples are $`<n_{\mathrm{APB95}}>=2.97\pm 1.29`$, $`<n_{\mathrm{KWK00}}>=2.88\pm 0.94`$ and $`<n_{\mathrm{MCH02}}>=1.02\pm 0.43`$ while the simulated sample has $`<n_{\mathrm{sim}}>=1.23\pm 0.40`$. The simulated distribution is consistent with the observations of MCH02 at one $`\sigma `$-level, which suggests that our objects are more similar to late-type galaxies. However, $`<n_{\mathrm{sim}}>`$ value also agrees at $`3`$ $`\sigma `$-level with the corresponding mean values of the early-type samples.

A key point in the formation of discs is the conservation of the specific angular momentum (FE80). Our SF algorithm has been implemented in such a way that it has allowed the formation of stellar bulges but without exhausting the gas reservoir (DT98). Hence, gaseous discs have been able to regenerate after mergers. In order to compare the angular momentum ($`j`$) of the three mass components: bulge, disc and dark matter halo, with observations, we plot in Fig. 4 their specific angular momentum ($`J`$=$`j/M`$) content as a function of mass ($`M`$). The mass and angular momentum of bulges have been estimated with stars within $`r<r_{\mathrm{eff}}`$, while those of the gaseous discs have been calculated within $`r<3.2r_\mathrm{d}`$. These criteria have been adopted following those used by observers. The mass and angular momentum of the dark matter haloes have been estimated at the virial radius. In this figure, we have also included two boxes that depict the observational region covered by spirals and ellipticals as given by Fall (1983). From this plot we can see that the disc components and the dark matter haloes have comparable specific angular momentum contents as predicted by the standard disc formation model of FE80. Conversely, stellar bulges have been formed from material that has lost most of its angular momentum.

From these results we conclude that the simulated GLOs have structural and dynamical parameters that statistically resemble those of current normal spirals at $`z=0`$. In the following section we will study the role played by mergers in the determination of these properties and try to understand the origin of correlations among the structural parameters such as $`B/D`$ versus $`n`$.

### 3.2 Analysis of progenitors during merger events

We carried out the bulge-disc (B+D) decomposition of the progenitor objects of the GLOs analised in the previous section during merger events at the four redshifts of interest: $`z_\mathrm{A}`$, $`z_\mathrm{B}`$, $`z_\mathrm{C}`$ and $`z_\mathrm{D}`$. In this section we will study how their structural parameters change during these processes.

#### 3.2.1 The shape parameter

In Fig. 5(a) we show the $`r_\mathrm{d}`$ versus $`r_\mathrm{b}`$ relation for all progenitors during merger events and at $`z=0`$. We have also included the observational data from a B+D decomposition from MGH99, KWK00 and MCH02, and from a double exponential decomposition by CdJB96. The combined observed sample of spirals shows a certain correlation between $`r_\mathrm{d}`$ and $`r_\mathrm{b}`$, although some of them have, for a given disc scalelength, a smaller bulge one than that expected from this relation. We also note that the scalelengths of our simulated disc structures show a similar behaviour. In order to individualize which observed and simulated objects belong to each of these two distributions of scalelengths, we segregate them according to their shape parameters as shown from Fig. 5(b) to Fig. 5(d). It is clear that as galaxies and simulated objects with large shape parameters are taken out, the relation for both, observed and simulated scalelengths, gets better defined.

We also note that the $`r_\mathrm{b}`$ versus $`r_\mathrm{d}`$ relation for GLOs and galaxies with bulges with $`n>1.5`$ shows no trend at all. It is encouraging that the observed and simulated disc structures show scalelength distributions with similar patterns, for both $`n1`$ or $`n1`$ bulge profiles. In the case of exponential bulges ($`0.9n1.1`$), the correlation is well defined with $`<r_\mathrm{b}/r_\mathrm{d}>=0.14\pm 0.05`$ and $`<r_\mathrm{b}/r_\mathrm{d}>=0.13\pm 0.06`$ for observations and simulations, respectively. Table 1 summarizes the mean values of $`r_\mathrm{b}`$, $`r_\mathrm{d}`$ and $`n`$ for GLOs at the four $`z`$ of interest.

Similary to Fig. 3(a) we display in Fig. 6 the histogram of shape parameters for all progenitors at all studied redshifts. Hence this distribution gathers the information on the $`n`$ parameters of different progenitors at different stages of evolution. From the comparison with that shown in Fig. 3 corresponding to the bulges of the final structures at $`z=0`$, we deduce that in the past of current GLOs, their progenitors went through stages where their bulges had $`n>1.5`$ parameters ($`35`$ per cent of the total sample). Meanwhile, GLOs at $`z=0`$ show a trend to have $`n1`$, with less than $`10`$ per cent being described by $`n>1.5`$.

In order to assess the effects of mergers on the shape parameters of the objects and to look for possible differences between the effects of secular evolution and fusions, we calculate the percentages corresponding to the increase and decrease of the shape parameter during the intervals $`z_\mathrm{A}`$-$`z_\mathrm{C}`$, $`z_\mathrm{C}`$-$`z_\mathrm{D}`$ and $`z_\mathrm{A}`$-$`z_\mathrm{D}`$ for SSBs and DSBs. The results are shown in Table 2. Note that we are taking the interval $`z_\mathrm{A}`$-$`z_\mathrm{C}`$ as the total period corresponding to the ODP. For both single and double SBs, we find that bulges tend to decrease their $`n`$ parameters during the ODP, and that this effect is more important for DSBs. In contrast, the actual fusion of the clumps ($`z_\mathrm{C}`$-$`z_\mathrm{D}`$) tends to make the $`n`$ value of the bulges larger.

Therefore we find that the physical encounter of the baryonic cores tends to produce bulges with larger shape parameters (i.e. more concentrated), and $`n1`$ profiles tend to be formed by secular evolution. These results together with the fact that at $`z=0`$ most of GLOs show $`n1.5`$, regardless of their past history, suggest that successive mergers have driven the bulges from exponential to more concentrated profiles, to come back to exponential ones, if there is available gas in the systems.

We think that the possibility of having secular evolution and core fusion associated to a merger helps to drive this cycle in the shape parameter. Note that the collisionless merger of stellar discs will not produce the same distributions. As a matter of fact, in order to detect this morphological loop, simulations should include gas dynamics and star formation since, otherwise, bulges can not be re-shaped back to $`n1`$ (A01).

Fig. 7 shows the variations of the $`n`$ parameter during the ODP (a) and those produced during the fusion of the baryonic clumps (b) as a function of $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$, where $`M_{\mathrm{star}}^z`$ is the stellar mass formed in the progenitors at $`z`$ and $`M_{\mathrm{star}}^{z=0}`$ is the total stellar mass of the GLOs at $`z=0`$. This ratio is an estimate of the presence of a stellar bulge since in these simulations the stars form preferentially in the dense regions. T02 found that the properties of the potential well can be directly linked to the stability of a system, in the sense that the shallower the potential well, the more susceptible the GLO to experience early gas inflows (i.e. secular evolution). The formation of stellar bulges has been found to inhibit such early gas inflows. This picture is in agreement with our results, where the major changes in the shape parameter during a merger event correspond to the GLOs that experience DSBs ($`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}0.45`$), which are less stable objects and can be strongly perturbed by the collision with a satellite. In contrast, more stable systems ($`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}>0.45`$) would not experience such important perturbations in its mass concentration during the ODP (Table 2).

We can also see from Fig. 7(a) that the most of the changes in $`n`$ are negative and can be linked to smaller $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$ values ($`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}0.45`$). This could be related to the fact that stellar bulges can provide stability to the systems, making them less vulnerable to the influence of a satellite, principally during the ODP. From Fig. 7(b) we see that the fusion of the baryonic cores produce more positive changes in the $`n`$ parameters and hence, a trend to increase it. Again, the larger variations are related to smaller bulges (i.e. smaller $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$).

We also analysed the possibility that the changes in the $`n`$ value during mergers are linked to the relative mass of the colliding systems, finding no correlation in contrast to the results of A01. However, these authors only considered collisionless mergers while in all our mergers dissipation plays an important role.

In Fig. 8 we have plotted the averaged $`n`$ parameters of the progenitors at different $`z`$ intervals. We can see that, although the bootstrap errors are high, there is a trend for the averaged $`n`$ to decrease toward $`n=1`$ for decreasing $`z`$. Nevertheless, at all redshifts we found the mean $`n`$ smaller than $`n=2`$, indicating that, in spite of the spread of the parameters, the majority of the bulges are best fitted with $`n1.5`$.

### 3.3 The $`r_\mathrm{b}/r_\mathrm{d}`$ ratio

It has been argued by some authors (e.g., CdJB96) that there is a restricted range of $`r_\mathrm{b}/r_\mathrm{d}`$ for spiral galaxies which has been interpreted as the result of secular evolution of the central mass concentrations. In Fig. 2 we showed the $`r_\mathrm{d}`$ versus $`r_\mathrm{b}`$ distribution for GLOs at $`z=0`$ and for observations. We have combined the observational results from MGH99, KWK00 and MCH02 which include spiral galaxies from Sa to Sc. For these distributions we find $`<r_\mathrm{b}/r_\mathrm{d}>=0.07\pm 0.02`$ (the boostrap error is $`\sigma _{\mathrm{bt}}=0.02`$) and $`<r_\mathrm{b}/r_\mathrm{d}>=0.10\pm 0.08`$ ($`\sigma _{\mathrm{bt}}=0.01`$) for the simulated and observed data, respectively, showing a good agreement.

However, we have already shown in Fig. 5 that the simulated and observed scalelength distributions show a clear trend only for those objects with bulges with $`0.9n1.1`$. Hence, we divided the observations and simulated data into three subsamples according to their $`n`$ parameters ($`n<0.9`$, $`0.9n1.1`$ and $`n>1.1`$) and estimate their mean $`r_\mathrm{b}/r_\mathrm{d}`$. Table 3 summarizes the results. First, we note that observational and simulated $`<r_\mathrm{b}/r_\mathrm{d}>`$ in the three subsamples are in very good agreement, supporting the hypothesis that GLOs in CDM models have observational counterpart. We also see that the larger dispersion is found for bulges with $`n<0.9`$ in both the observed and simulated data, and that these subsamples show higher $`<r_\mathrm{b}/r_\mathrm{d}>`$ than the subsamples with $`n1`$, suggesting that these systems systematically have larger $`r_\mathrm{b}`$ (for this subsample: $`<r_\mathrm{b}>=0.89\pm 0.10`$ and $`<r_\mathrm{d}>=4.77\pm 0.40`$).

On the other hand, we find that systems with bulges better fitted by $`n>1.1`$ show the smallest $`<r_\mathrm{b}/r_\mathrm{d}>`$. Hence for a given disc scalelength, the larger the concentration of the bulges, the smaller its scalelength (for this subsample: $`<r_\mathrm{b}>=0.23\pm 0.03`$ and $`<r_\mathrm{d}>=5.52\pm 0.30`$). If we assume that the $`n`$ parameter correlates with morphology as it is indicated by several works (e.g., de Jong 1996; Graham & Prieto 1999, hereafter GP99; GdB01), then the scalelength ratios should correlate with $`n`$. Fig. 9 shows the $`r_\mathrm{b}/r_\mathrm{d}`$ ratio as a function of the shape parameter for the simulated objects during mergers and for $`z=0`$. We have also included observational data from KWK00 and MCH02. We find an anticorrelation between these two parameters, in agreement with observations. Assuming that the shape parameter is an indicator of Hubble type, this anticorrelation indicates that late-type spirals ($`n1.5`$) have larger $`r_\mathrm{b}/r_\mathrm{d}`$ ratios than early-type spirals ($`n>1.5`$). We also note that the relation for the observed and simulated ratios seems to flatten for $`n1`$. A simple extrapolation gives a value of $`r_\mathrm{b}/r_\mathrm{d}0.40`$.

Because the shape parameter and the scalelength of the Sérsic law are coupled parameters, this trend of $`r_\mathrm{b}/r_\mathrm{d}`$ with $`n`$ could be the result of the fitting formula. However, Trujillo, Graham & Caon (2001, hereafter TGC01) argued that there is a physical trend stronger than that resulting from the fitting formula. Because we use the same expression for all fittings, the fact that there is a segregation in the scalelength distribution according to the value of $`n`$ suggests that this trend could have physical basis and could be related to a difference in the evolution of the GLOs or in the efficiency of some physical process. This is confirmed when $`r_{\mathrm{eff}}`$ is used instead of $`r_\mathrm{b}`$ for which we find $`<r_{\mathrm{eff}}/r_\mathrm{d}>=0.29`$ ($`\sigma _{\mathrm{bt}}`$=0.04) for $`n<0.9`$ , $`<r_{\mathrm{eff}}/r_\mathrm{d}>=0.26`$ ($`\sigma _{\mathrm{bt}}`$=0.03) for $`0.9n1.1`$, and $`<r_{\mathrm{eff}}/r_\mathrm{d}>=0.30`$ ($`\sigma _{\mathrm{bt}}`$=0.05) for $`n>1.1`$ (the combined observational sample of KWK00 and MCH02 shows $`<r_{\mathrm{eff}}/r_\mathrm{d}>=0.23`$ ($`\sigma _{\mathrm{bt}}`$=0.02) for $`n<0.9`$, $`<r_{\mathrm{eff}}/r_\mathrm{d}>=0.28`$ ($`\sigma _{\mathrm{bt}}`$=0.04) for $`0.9n1.1`$, and $`<r_{\mathrm{eff}}/r_\mathrm{d}>=0.35`$ ($`\sigma _{\mathrm{bt}}`$=0.03) for $`n>1.1`$). We have used the approximation of Ciotti & Bertin (1999) to transform simulated $`r_\mathrm{b}`$ into $`r_{\mathrm{eff}}`$.

Let us now look at the mean $`<r_\mathrm{b}/r_\mathrm{d}>`$ during mergers. In Table 3 we show the mean $`r_\mathrm{b}/r_\mathrm{d}`$ at $`z_\mathrm{A}`$ and $`z_\mathrm{D}`$ for SSBs and DSBs. We find no differences between the mean $`r_\mathrm{b}/r_\mathrm{d}`$ but both samples agree in increasing the mean ratios and the dispersions at $`z_\mathrm{D}`$. Fig. 10 shows the distributions of scalelength during mergers. Note that, independently of the stage of evolution of the mergers, the distributions have similar characteristics and are within the observed range at $`z=0`$. Finally in Fig. 9, in the small box, we show how the averaged $`r_\mathrm{b}/r_\mathrm{d}`$ moves throughout mergers associated with single (thin line) and double (thick line) starbursts. It is clear that early gas inflows during the ODP tend to decrease $`n`$ in both cases but for DSBs the changes are more important, as expected from their stability properties (see Fig. 7 and T02). Conversely, the fusion of the cores moves the distribution toward higher $`<n>`$.

### 3.4 The $`B/D`$ ratio

The luminosity $`B/D`$ ratio has been traditionally used as an indicator of morphology (e.g., Andredakis & Sanders 1994). The correlation between $`B/D`$ and $`n`$ has been shown by APB95 and recently confirmed by TGC01. Linear regression through observations show $`(B/D)/n=1.0\pm 0.3`$ ($`\sigma _{\mathrm{bt}}=0.3`$) for APB95, $`0.8\pm 0.5(\sigma _{\mathrm{bt}}=0.5)`$ for KWK00 and $`1.7\pm 0.3(\sigma _{\mathrm{bt}}=0.3)`$ for MCH02. The combined sample has $`(B/D)/n=1.2\pm 0.1(\sigma _{\mathrm{bt}}=0.1)`$.

Note that in the models we have masses instead of luminosities and that the mass-to-light ratios of the discs ($`\gamma _\mathrm{d}`$) and bulges ($`\gamma _\mathrm{b}`$) could be different because of their different stellar populations. However, their correct values and/or possible dependence on radius as well as redshift are still unclear. In order to match the luminosity $`B/D`$ we have to rescale the mass $`B/D`$ ratio by a constant factor of 20. We address that we are using the same mass-to-light ratios for all GLOs regardless of the redshift111Nevertheless, we found that very similar $`\gamma _\mathrm{b}/\gamma _\mathrm{d}`$ ratios are needed to match the observed luminosity $`B/D`$ for progenitors at different $`z`$, although our dispersion is quite high. Note that we have used the total baryonic masses to estimate luminosities. In a more realistic model including energy feedback, some fraction of the cold gas would be reheated. Hence, the values required for our models to match the observed luminosities could be taken as upper limits to the mass-to-light ratios., disc scalelength or observational band.

In Fig. 11 we show the $`B/D`$ ratio for the simulated GLOs during merger events at the four redshifts of interest and the observational data from APB95, KWK00 and MCH02. The low statistical number restrict the use of linear regression through the simulated data. However, we see that, at any time during merger events, the simulated values are within the observational range.

In Fig.12 we plot the variations $`n_{z_i}n_{z_j}`$ versus $`B/D_{z_i}B/D_{z_j}`$ where $`i`$/$`j`$ can be A, C or D. From this figure we see that GLOs which experience gas inflows show more important changes in both $`n`$ and $`B/D`$ during the ODP. In the case of those GLOs experiencing SSBs (e.g., no early gas inflows) there is no important change in the shape parameter or the $`B/D`$ ratios during the ODP. We also note that the changes in $`n`$ and $`B/D`$ do not correlate. Conversely, during the actual fusion of the objects (from $`z_\mathrm{C}`$ to $`z_\mathrm{D}`$) both types of GLOs experience large changes in both parameters in the expected direction: the larger the increase in the shape parameter, the larger the corresponding increase of the $`B/D`$ ratio or vice-versa. Fig. 12(c) shows the overall changes over the merger events. This distribution has higher dispersion than that obtained during the fusions because it also combines the effects of secular evolution. However, the overall effects of mergers is to settle the correlation between $`B/D`$ and $`n`$ so that more concentrated bulges (higher shape parameters) have the larger $`B/D`$ ratios. However, we note that, although the range of values can be matched, the simulated $`B/D`$ do not show the same correlation signal ($`r=0.20`$) than the observations ($`r=0.55`$).

## 4 Conclusions

We have studied the properties of the baryon distributions in galaxy-like objects focusing on the effects of mergers. We resort to observations of spiral galaxies of different morphology to constrain our findings.

We found that on average, galactic objects formed in hierarchical clustering scenarios reproduce the angular momentum and structural parameter distributions of spiral galaxies, if a stellar bulge is allowed to form and early gas depletion is avoided. We have succeeded in these two aspects but at the expense of inhibiting star formation on discs. A consistent treatment of energy feedback may help to remove this caveat.

These simulations have allowed us to study how mergers change the distribution of baryons by analysing the evolution of their structural parameters. We found that, on average, galactic objects tend to have nearly exponential bulges at all redshift. However, note that these simulations produce bulges with shape parameters in the range $`0.54`$.

For those systems with $`n1`$ bulges we found a correlation among their bulge and disc scalelengths in very good agreement with observations. However, for $`n>1`$ the scalelength distribution is disorder and displaced to smaller bulge scalelenghts. Observations show the same behaviour. The opposite distribution is found for systems with $`n<1`$ bulges which have larger $`r_\mathrm{b}/r_\mathrm{d}`$ values and larger dispersions. We found the disc scalelengths to be approximately independent of $`n`$ parameters, so that this correlation implies that more concentrated objects have smaller bulge scalelengths (or large effective radius). We also found a maximum $`r_\mathrm{b}/r_\mathrm{d}`$ of $`0.40`$ for $`n0`$. Higher numerical resolution simulations are needed to study the formation of such low mass surface profiles.

In order to test our results for low resolution problems in the determination of the structural parameters of the galaxy-like objects, we followed Steinmetz & Müller (1994) and used the bootstrap technique to estimate the effects of low particle number statistics within the objects. We calculated an accuracy better than 25 per cent for the scalelengths and shape parameters.

We found that gas inflows during the orbital decay phase tend to produce important changes in the mass distributions generating $`n1`$ profiles, while the fusions of the baryonic cores tend to increase the $`n`$ parameter. As a consequence, a morphological loop can be driven by mergers which might be responsible of triggering secular evolution as well as of the baryonic core fusions. The triggering of secular evolution is found to be linked to the presence of a stellar bulge so that systems with well-formed stellar bulges do not experience early gas inflows. Hence, the pace of this morphological loop could be regulated by the properties of the galactic central potential wells which are also affected by the merger history of the objects (see also Tissera & Domínquez-Tenreiro 1998 and Tissera 2000).

We found that the simulated mass bulge-to-disc ratios are within observed range. It is also noted that during the ODP larger changes are observed in the $`B/D`$ ratios of those objects that experience gas inflows. The changes during this period are, however, quite disorder. It is at the fusion of the baryonic clumps that changes in the shape parameters are correlated with changes in the bulge-to-disc ratio. In our simulations, the actual fusions are responsible of significantly increasing the mass concentration at the centre. Hence, we found that the fusion of the baryonic cores could be the process that determine the observed correlation between the luminosity $`B/D`$ ratio and the shape parameter or morphological type. However, we found no dependence on the relative masses of the colliding objects.

Overall, our results indicate that the morphological properties of galactic objects are the result of their merger histories within a hierarchical clustering scenario. Based on the good agreement found so far with observations we support the hypothesis of mergers as the main morphological driver along the Hubble sequence. Consequently, the particular and detailed history of substructure aggregation could be a key point in the determination of the astrophysical properties of galaxies.

## Acknowledgments

We thank the anonymous referee for useful comments that helped to improve this paper. We are grateful to S. Courteau and his collaborators for making available their observational data prior to publication. We acknowledge A. Graham and S. Courteau for useful discussions. We thank Max-Planck Institute for Astrophysics for their hospitality, where this paper was completed. This work was partially supported by the Consejo Nacional de Investigaciones Científicas y Técnicas, Agencia de Promoción de Ciencia y Tecnología and Fundación Antorchas.


format: text

text:
# 1 Introduction

## 1 Introduction

The rare decays of $`B`$ mesons provide an excellent means to further explore the Standard Model (SM) and to detect New Physics. In particular the Cabibbo-favoured radiative $`bs\gamma `$ modes belong to the small number of rare decays that are experimentally accessible already at present \[??\] and give severe bounds on the parameter space of New Physics scenarios. On the theoretical side both the inclusive \[??\] and exclusive \[??\] $`bs\gamma `$ decays are now known at next-to-leading-logarithmic (NLL) accuracy.

For the double radiative decays $`B\gamma \gamma `$ experimentally so far only upper limits on the branching fractions exist:

$$\begin{array}{ccccc}\hfill B(B_s\gamma \gamma )& <& 1.4810^4\hfill & \text{at 90\% C.L.}\hfill & \text{[11]}\hfill \\ \hfill B(B_d\gamma \gamma )& <& 1.710^6\hfill & \text{at 90\% C.L.}\hfill & \text{[12]}\hfill \end{array}$$
(1)
The Standard Model expectations are roughly two orders below the current upper limits. These decays have a rather clean experimental signature. Furthermore they are of interest because they could provide useful tests of QCD dynamics in $`B`$ decays. The $`B\gamma \gamma `$ modes realize the exceptional situation of nontrivial QCD dynamics related to the decaying $`B`$, in conjunction with a completely nonhadronic final state and simple two-body kinematics. In principle they also probe the CKM parameters $`V_{ts}`$ and $`V_{td}`$ and could allow us to study CP violating effects as the two-photon system can be in a CP-even or CP-odd state. While contributions of New Physics to $`B_s\gamma \gamma `$ are strongly constrained by the measurements of $`bs\gamma `$ transitions, a sizable enhancement of $`B_d\gamma \gamma `$ from effects beyond the Standard Model is still conceivable.

The theoretical treatments of $`B_s\gamma \gamma `$ performed so far \[??\] all had to employ hadronic models to describe the $`B_s`$ meson bound state dynamics. A clear separation of short- and long-distance dynamics and a distinction of model-dependent and model-independent features were therefore not possible. This concerns especially the dynamics of the light-quark constituent inside the $`B`$, but also contributions from intermediate $`J/\psi `$, $`\eta _c`$, $`\varphi `$ or $`D_s^{()}`$ meson states, which have been discussed as sources of potentially important long-distance effects. Another process related to $`B\gamma \gamma `$ is the decay $`Bl^+l^{}\gamma `$. This decay has recently been studied in .

We present in this paper a systematic analysis of the exclusive double radiative decays $`B_{s,d}\gamma \gamma `$ in QCD, based on the heavy quark limit $`m_b\mathrm{\Lambda }_{QCD}`$. This limit allows us to give a factorization formula for the hadronic matrix elements of local operators in the weak Hamiltonian . In this manner we can systematically separate perturbatively calculable hard scattering kernels from the nonperturbative $`B`$-meson light-cone distribution amplitude (LCDA). Power counting in $`\mathrm{\Lambda }_{QCD}/m_b`$ allows us to identify leading and subleading contributions to $`B\gamma \gamma `$. Only one diagram contributes at leading power, but the most important subleading contributions can also be calculated. The inclusion of these corrections is used to estimate CP asymmetries in $`B\gamma \gamma `$ and one-particle irreducible two-photon emission from light-quark loops in $`D\gamma \gamma `$. These corrections represent the quark-level analogue of so-called long-distance contributions in these decays.

The remainder of this paper is organized as follows. In section 2 we give the basic formulas necessary for the calculation of the observables of interest. Section 3 contains our discussion of leading-power and section 4 the one of subleading-power contributions. In section 5 we comment on long-distance contributions in $`B\gamma \gamma `$ and $`D\gamma \gamma `$ decays. Numerical results are presented and discussed in section 6. Section 7 contains our conclusions.

## 2 Basic formulas

The effective Hamiltonian for $`bs\gamma \gamma `$ is identical to the one for $`bs\gamma `$ transitions. This is because the equations of motion (EOM) can be used to reduce the basis with operators containing two photon fields . Up to corrections of order $`1/M_W^2`$ the effective Hamiltonian thus reads

$$_{eff}=\frac{G_F}{\sqrt{2}}\underset{p=u,c}{}\lambda _p^{(s)}\left[C_1Q_1^p+C_2Q_2^p+\underset{i=3,\mathrm{},8}{}C_iQ_i\right]+h.c.$$
(2)
where

$$\lambda _p^{(s)}=V_{ps}^{}V_{pb}$$
(3)
The operators are given by

$`Q_1^p`$ $`=`$ $`(\overline{s}p)_{VA}(\overline{p}b)_{VA}`$ (4)
$`Q_2^p`$ $`=`$ $`(\overline{s}_ip_j)_{VA}(\overline{p}_jb_i)_{VA}`$ (5)
$`Q_3`$ $`=`$ $`(\overline{s}b)_{VA}{\displaystyle \underset{q}{}}(\overline{q}q)_{VA}`$ (6)
$`Q_4`$ $`=`$ $`(\overline{s}_ib_j)_{VA}{\displaystyle \underset{q}{}}(\overline{q}_jq_i)_{VA}`$ (7)
$`Q_5`$ $`=`$ $`(\overline{s}b)_{VA}{\displaystyle \underset{q}{}}(\overline{q}q)_{V+A}`$ (8)
$`Q_6`$ $`=`$ $`(\overline{s}_ib_j)_{VA}{\displaystyle \underset{q}{}}(\overline{q}_jq_i)_{V+A}`$ (9)
$`Q_7`$ $`=`$ $`{\displaystyle \frac{e}{8\pi ^2}}m_b\overline{s}_i\sigma ^{\mu \nu }(1+\gamma _5)b_iF_{\mu \nu }`$ (10)
$`Q_8`$ $`=`$ $`{\displaystyle \frac{g}{8\pi ^2}}m_b\overline{s}_i\sigma ^{\mu \nu }(1+\gamma _5)T_{ij}^ab_jG_{\mu \nu }^a`$ (11)
The most important operators are the magnetic penguin operator $`Q_7`$ and the four-quark operators $`Q_{1,2}^p`$. The sign conventions for the electromagnetic and strong couplings correspond to the covariant derivative $`D_\mu =_\mu +ieQ_fA_\mu +igT^aA_\mu ^a`$. With these definitions the coefficients $`C_{7,8}`$ are negative in the Standard Model, which is the choice generally adopted in the literature. The effective Hamiltonian for $`bd\gamma \gamma `$ is obtained from (211) by the replacement $`sd`$. The Wilson coefficients $`C_i`$ in (2) are known at next-to-leading order .

The amplitude for the $`B\gamma \gamma `$ decay has the general structure

$`𝒜\left(\overline{B}\gamma (k_1,ϵ_1)\gamma (k_2,ϵ_2)\right){\displaystyle \frac{G_F}{\sqrt{2}}}{\displaystyle \frac{\alpha }{3\pi }}f_B{\displaystyle \frac{1}{2}}\gamma \gamma |A_+F_{\mu \nu }F^{\mu \nu }iA_{}F_{\mu \nu }\stackrel{~}{F}^{\mu \nu }|0`$ (12)
$`={\displaystyle \frac{G_F}{\sqrt{2}}}{\displaystyle \frac{\alpha }{3\pi }}f_B\left[A_+\left(2k_1ϵ_2k_2ϵ_1m_B^2ϵ_1ϵ_2\right)2iA_{}\epsilon (k_1,k_2,ϵ_1,ϵ_2)\right]`$
Here $`F^{\mu \nu }`$ and $`\stackrel{~}{F}^{\mu \nu }`$ are the photon field strength tensor and its dual, where

$$\stackrel{~}{F}^{\mu \nu }=\frac{1}{2}\epsilon ^{\mu \nu \lambda \rho }F_{\lambda \rho }$$
(13)
with $`\epsilon ^{0123}=1`$. The decay rate is then given by

$$\mathrm{\Gamma }(\overline{B}\gamma \gamma )=\frac{G_F^2m_B^3f_B^2\alpha ^2}{288\pi ^3}\left(|A_+|^2+|A_{}|^2\right)$$
(14)
In the heavy-quark limit we propose a factorization formula for the hadronic matrix elements of the operators in (2):

$$\gamma (ϵ_1)\gamma (ϵ_2)|Q_i|\overline{B}=f_B_0^1𝑑\xi T_i^{\mu \nu }(\xi )\mathrm{\Phi }_B(\xi )ϵ_{1\mu }ϵ_{2\nu }$$
(15)
where the $`ϵ_i`$ are the polarization 4-vectors of the photons and $`\mathrm{\Phi }_B\mathrm{\Phi }_{B1}`$ is the leading twist light-cone distribution amplitude of the $`B`$ meson. The latter quantity is a universal, nonperturbative object. It is defined via the light-cone projector for the $`B`$ meson at leading power

$$0|b(0)\overline{q}(z)|\overline{B}(p)=\frac{if_B}{4}(\overline{)}p+m_b)\gamma _5_0^1𝑑\xi e^{i\xi p_+z_{}}[\mathrm{\Phi }_{B1}(\xi )+\overline{)}n\mathrm{\Phi }_{B2}(\xi )]$$
(16)
The functions $`\mathrm{\Phi }_{B1,B2}(\xi )`$ describe the distribution of light-cone momentum fraction $`\xi =l_+/p_+`$ of the spectator quark $`q`$ with momentum $`l`$ inside the $`B`$ meson (only $`\mathrm{\Phi }_{B1}`$ contributes in 15). Here light-cone components of four-vectors $`v`$ are defined by

$$v_\pm =\frac{v^0\pm v^3}{\sqrt{2}}$$
(17)
The wave functions are highly asymmetric with $`\xi =𝒪(\mathrm{\Lambda }_{QCD}/m_b)`$. They are normalized as

$$_0^1𝑑\xi \mathrm{\Phi }_{B1}(\xi )=1_0^1𝑑\xi \mathrm{\Phi }_{B2}(\xi )=0$$
(18)
If the light-like vector $`n`$ in (16) is chosen appropriately parallel to one of the 4-momenta of the photons, only the first negative moment of $`\mathrm{\Phi }_{B1}(\xi )`$ appears, which we parametrize by a quantity $`\lambda _B=𝒪(\mathrm{\Lambda }_{QCD})`$, i.e.

$$_0^1𝑑\xi \frac{\mathrm{\Phi }_{B1}(\xi )}{\xi }=\frac{m_B}{\lambda _B}$$
(19)
Because there are no hadrons in the final state, only one type of hard-scattering kernel $`T`$ (type II or hard-spectator contribution ) enters the factorization formula. The QCD factorization formula (15) holds up to corrections of relative order $`\mathrm{\Lambda }_{QCD}/m_b`$. The form of (15) with a simple convolution over the light-cone variable $`\xi `$ is appropriate for the lowest (leading logarithmic) order in $`\alpha _s`$, which we will use in the present analysis. A generalization to include transverse-momentum variables is likely to be necessary at higher orders in QCD .

We conclude this section with a brief discussion of CP violation in $`B\gamma \gamma `$. The subscripts $`\pm `$ on $`A_\pm `$, defined in (12) for $`\overline{B}\gamma \gamma `$, denote the CP properties of the corresponding two-photon final states, which are eigenstates of CP: $`A_{CP=+1}`$ is proportional to the parallel spin polarization $`\stackrel{}{ϵ}_1\stackrel{}{ϵ}_2`$ and $`A_{CP=1}`$ is proportional to the perpendicular spin polarization $`\stackrel{}{ϵ}_1\times \stackrel{}{ϵ}_2`$ of the photons. In addition to $`A_\pm `$ we introduce the CP conjugated amplitudes $`\overline{A}_\pm `$ for the decay $`B\gamma \gamma `$ (decaying $`\overline{b}`$ anti-quark). Then the deviation of the ratios

$$r_{CP}^\pm =\frac{|A_\pm |^2|\overline{A}_\pm |^2}{|A_\pm |^2+|\overline{A}_\pm |^2}$$
(20)
from zero is a measure of direct CP violation. We shall focus here on this direct effect, which is specific for the $`B\gamma \gamma `$ decay, and will not consider CP violation originating in $`B`$$`\overline{B}`$ mixing.

Due to the unitarity of the CKM matrix we can parametrize

$`A_\pm `$ $`=`$ $`\lambda _ua_u^\pm e^{i\alpha _u^\pm }+\lambda _ca_c^\pm e^{i\alpha _c^\pm }`$ (21)
$`\overline{A}_\pm `$ $`=`$ $`\lambda _u^{}a_u^\pm e^{i\alpha _u^\pm }+\lambda _c^{}a_c^\pm e^{i\alpha _c^\pm }`$ (22)
where $`a_{u,c}^\pm `$ are real hadronic matrix elements of weak transition operators and the $`\alpha _{u,c}^\pm `$ are their CP-conserving phases. We then have

$$r_{CP}^\pm =\frac{2\mathrm{Im}(\lambda _u\lambda _c^{})a_u^\pm a_c^\pm \mathrm{sin}(\alpha _c^\pm \alpha _u^\pm )}{|\lambda _u|^2a_u^{\pm 2}+|\lambda _c|^2a_c^{\pm 2}+2\mathrm{Re}(\lambda _u\lambda _c^{})a_u^\pm a_c^\pm \mathrm{cos}(\alpha _c^\pm \alpha _u^\pm )}$$
(23)
The weak phase differences read

$$\mathrm{Im}(\lambda _u^{(q)}\lambda _c^{(q)})=|\lambda _u^{(q)}||\lambda _c^{(q)}|\mathrm{sin}\gamma $$
(24)
with the negative sign for $`q=s`$ and the positive sign for $`q=d`$. The strong phase differences $`\mathrm{sin}(\alpha _c^\pm \alpha _u^\pm )`$ arise from final state interaction (FSI) effects, generated for instance via the Bander-Silverman-Soni (BSS) mechanism .

## 3 $`B\gamma \gamma `$ at leading power

To leading power in $`\mathrm{\Lambda }_{QCD}/m_b`$ and in the leading logarithmic (LL) approximation of QCD only one diagram contributes to the amplitude for $`B\gamma \gamma `$. It is the one-particle reducible (1PR) diagram with the electromagnetic penguin operator $`Q_7`$, where the second photon is emitted from the $`s`$-quark line. An illustration is given in Fig. 1.

Evaluating the factorization formula (15) for this diagram leads to the matrix element

$`\gamma (k_1,ϵ_1)\gamma (k_2,ϵ_2)|Q_7|\overline{B}=`$
$`i{\displaystyle \frac{f_B\alpha }{\pi }}Q_s{\displaystyle _0^1}𝑑\xi {\displaystyle \frac{\mathrm{\Phi }_{B1}(\xi )}{\xi }}\left(2k_1ϵ_2k_2ϵ_1m_B^2ϵ_1ϵ_22i\epsilon (k_1,k_2,ϵ_1,ϵ_2)\right)`$
Together with (2) and the definition of $`A_\pm `$ in (12) this implies

$$A_\pm =\lambda _t^{(q)}C_7\frac{m_B}{\lambda _B}$$
(26)
where CKM unitarity, $`\lambda _u^{(q)}+\lambda _c^{(q)}=\lambda _t^{(q)}`$, has been used. After summing over the photon polarizations the branching fraction to leading power becomes

$$B(B_q\gamma \gamma )=\tau _{B_q}\frac{G_F^2m_{B_q}^5}{192\pi ^3}|\lambda _t^{(q)}|^2C_7^2\frac{4\alpha ^2f_{B_q}^2}{3\lambda _{B_q}^2}$$
(27)
where $`q=d,s`$ for the decay of a $`B_d`$ or a $`B_s`$ meson, respectively.

In the present approximation the strong-interaction matrix elements multiplying $`\lambda _u^{(q)}`$ and $`\lambda _c^{(q)}`$ are identical and have no relative phase. The CP-violating quantities $`r_{CP}^\pm `$ defined in (20) therefore vanish for the strictly leading-power result.

One may compare the behaviour of the $`B\gamma \gamma `$ decay amplitude in the heavy-quark limit with that of the rare decays $`BV\gamma `$ (with $`V`$ a vector meson) and $`B\pi \pi `$. Disregarding CKM factors we can write

$$\begin{array}{ccc}𝒜(B\gamma \gamma )\hfill & G_Fe^2f_BA_\pm m_B^2\hfill & G_Fe^2\mathrm{\Lambda }^{1/2}m_b^{5/2}\hfill \\ 𝒜(BV\gamma )\hfill & G_Fem_bF_Vm_B^2\hfill & G_Fe\mathrm{\Lambda }^{3/2}m_b^{3/2}\hfill \\ 𝒜(B\pi \pi )\hfill & G_Fe^0f_\pi F^{B\pi }m_B^2\hfill & G_Fe^0\mathrm{\Lambda }^{5/2}m_b^{1/2}\hfill \end{array}$$
(28)
where we have used that the form factors scale as ($`\mathrm{\Lambda }=\mathrm{\Lambda }_{QCD}`$)

$$A_\pm \frac{m_b}{\mathrm{\Lambda }}f_B\frac{\mathrm{\Lambda }^{3/2}}{m_b^{1/2}}F_VF^{B\pi }\frac{\mathrm{\Lambda }^{3/2}}{m_b^{3/2}}$$
(29)
We observe that the ratio of the three amplitudes in (28) is $`1:(\mathrm{\Lambda }/m_b):(\mathrm{\Lambda }/m_b)^2`$. Thus $`𝒜(B\gamma \gamma )`$ is leading in the heavy-quark limit, while the other two are successively stronger suppressed. However, this hierarchy is compensated by the opposite pattern $`e^2:e:1`$ in the electromagnetic coupling $`e`$.

## 4 Contributions to $`B\gamma \gamma `$ with power suppression

Subleading contributions come from the 1PR diagram, where the second photon is emitted off the $`b`$ quark line (Fig. 2), and from the one-particle irreducible diagram (1PI) (see Fig. 3).

To estimate the subleading 1PR contribution, we simply evaluate the graphs in Fig. 2 using the $`B`$-meson projector in (16). The result equals $`\lambda _B/m_B`$ times the leading-power expression in (3), clearly showing the power suppression of this mechanism. Other corrections of the same order can arise from higher-twist terms in the $`B`$-meson wave function. We shall neglect all those subleading terms for the matrix element of $`Q_7`$, keeping in mind that they could naturally contribute relative corrections of order $`10\%`$. We stress, however, that subleading effects in $`Q_7`$ contribute equally to the up- and charm-quark components of the amplitude (see (2)). Therefore they do not give rise, in particular, to relative FSI phases between these sectors, which would affect direct CP violation. The same is true for perturbative QCD corrections to $`Q_7`$.

We next consider more closely the diagrams in Fig. 3. These 1PI contributions, which come from the matrix elements of four-quark operators in (2), are of special interest for two reasons. First, they are the basic effects responsible for a difference between the up- and charm-quark sectors of the amplitude, including rescattering phases. Second, they represent the parton-level processes that are dual to $`B\gamma \gamma `$ amplitudes from $`Q_i`$ with hadronic intermediate states ($`J/\mathrm{\Psi }`$, $`D^{()}`$, etc.), which are commonly considered as generic long-distance contributions. Investigating the amplitude in Fig. 3 will therefore shed light on this class of effects in $`B\gamma \gamma `$.

It shall now be argued that the 1PI contributions are calculable using QCD factorization. It is known that the one-loop contribution in Fig. 3 is free of infrared (IR) singularities. In addition we will show that also at $`𝒪(\alpha _s)`$ there are no collinear or soft IR divergences at leading power in $`1/m_b`$. By leading power we here mean the lowest nonvanishing order in the power expansion, where the entire 1PI contribution starts only at subleading power with respect to the dominant mechanism in Fig. 1. This suppression will be apparent from the explicit expressions given in (40) below. The relevant two-loop diagrams at $`𝒪(\alpha _s)`$ are shown in Fig. 4.

Following the approach explained in the second reference of  for the case of $`BD\pi `$ decays, we demonstrate the absence of divergences for the hard-scattering part of the amplitude at $`𝒪(\alpha _s)`$. As a consequence the matrix elements of four-quark operators can be written in factorized form as in (15), where in the present case the kernel $`T_i^{\mu \nu }`$ is independent of $`\xi `$:

$$\gamma (ϵ_1)\gamma (ϵ_2)|Q_{1,2}|\overline{B}=f_BT_{1,2}^{\mu \nu }ϵ_{1\mu }ϵ_{2\nu }$$
(30)
The proof proceeds by identifying the potentially IR singular regions in the loop-momentum variables $`s`$ and $`t`$ and by determining the degree of IR divergence through power counting. We distinguish the following cases: (i) $`s`$, $`t`$ soft, where all components scale like $`s`$, $`t\mathrm{\Lambda }m_b`$. (ii) $`s`$, $`t`$ collinear with either photon momentum $`k_1`$ or $`k_2`$, for instance

$$s=\alpha k_1+s_{}+\alpha _2k_2t=\beta k_1+t_{}+\beta _2k_2$$
(31)
describing the region where $`s`$ and $`t`$ are collinear with $`k_1`$. Similar relations, with $`k_1`$ and $`k_2`$ interchanged, apply to the case where the loop momenta become collinear with $`k_2`$. The parameters in (31) scale as

$$\alpha ,\beta 1,s_{},t_{}\mathrm{\Lambda },\alpha _2,\beta _2\frac{\mathrm{\Lambda }^2}{m_b^2}$$
(32)
and we have

$$k_{1,2}^2=k_{1,2}s_{}=k_{1,2}t_{}=0$$
(33)
We consider explicitly the case where the quark inside the loop is massless, which is the most IR singular situation. The momentum $`l`$ of the light external quark is counted as a soft quantity $`l\mathrm{\Lambda }`$. The dominant region for the diagrams in Figs. 3, 4 arises when both loop momenta are hard, $`s`$, $`tm_b`$. The diagrams then scale as $`m_b`$. To demonstrate IR finiteness, one has to show that the potentially singular regions scale at least with one power of $`\mathrm{\Lambda }`$ and thus are suppressed relative to the hard contribution. Scaling as $`\mathrm{\Lambda }^0m_b`$ in a singular region would indicate a logarithmic divergence. The relevant cases for $`s`$ and $`t`$ are the soft-soft, soft-hard, hard-soft, collinear-collinear, collinear-hard, hard-collinear, soft-collinear and collinear-soft regions. For the collinear regions the momenta can be either $`k_1`$ or $`k_2`$. In the collinear-collinear case one needs to consider primarily that both momenta $`k_1`$, or both $`k_2`$. If one momentum is $`k_1`$ and the other $`k_2`$, the contribution is always less singular in the infrared. In the second reference of , besides the soft-collinear regions, a supersoft-collinear scaling had to be discussed, where $`s`$, $`t\mathrm{\Lambda }^2`$ for supersoft momenta. In the present context this scaling is not more singular than the soft-collinear one and will thus not be considered further.

Power counting then shows that all diagrams in Fig. 4 are IR finite at leading power in the soft-soft, soft-hard, hard-soft, collinear-hard, hard-collinear, soft-collinear and collinear-soft regions. In some cases one encounters a superficial logarithmic divergence, which vanishes when the relations $`k_{1,2}^2=0`$ and $`k_1ϵ_1=k_2ϵ_2=0`$ are used.

On the other hand, in the collinear-collinear region (superficial) linear divergences $`\mathrm{\Lambda }^1`$ may appear. An example is diagram (f) in Fig. 4 when both $`s`$ and $`t`$ are collinear to $`k_1`$. The diagram has the form

$$D_{(f)}=d^4sd^4t\frac{\gamma _\sigma (\overline{)}l\overline{)}s)\gamma ^\lambda (1\gamma _5)(\overline{)}t\overline{)}k_2)\gamma _\nu \overline{)}t\gamma _\mu (\overline{)}t+\overline{)}k_1)\gamma ^\sigma (\overline{)}t+\overline{)}k_1\overline{)}s)\gamma _\lambda (1\gamma _5)}{s^2(ls)^2(tk_2)^2t^2(t+k_1)^2(t+k_1s)^2}$$
(34)
Expanding the numerator in powers of $`\mathrm{\Lambda }`$ using (31) one finds that the superficial linear divergence drops out. This leads to

$$D_{(f)}=d^4sd^4t\frac{\gamma _\sigma (\overline{)}l\overline{)}s)}{(ls)^2}\mathrm{\Gamma }_Rk_1^\sigma $$
(35)
with the remainder of the diagram $`\mathrm{\Gamma }_R\mathrm{\Lambda }^7`$. Together with collinear phase space, $`d^4sd^4t\mathrm{\Lambda }^8`$, the integrand then appears to scale as $`\mathrm{\Lambda }^0m_b`$. In fact, since $`\overline{)}k_1(\overline{)}l\overline{)}s)m_b\mathrm{\Lambda }`$, $`(ls)^2m_b\mathrm{\Lambda }`$ for $`s=\alpha k_1+𝒪(\mathrm{\Lambda })`$, the contribution behaves only as $`D_{(f)}\mathrm{\Lambda }`$. In a similar way all other contributions from the collinear-collinear regions can be shown to be power-suppressed.

In addition to the graphs displayed in Fig. 4 there are also diagrams in which the gluon connects only to the lines inside the triangular quark loop (vertex and self-energy corrections). Also these diagrams are found to be IR finite at leading power.

Finally, in contrast, one-gluon exchange among the external lines does generate leading IR divergences. However, these singularities match those appearing in the perturbative evaluation of the $`B`$-meson decay constant $`f_B`$ at the same order, in accordance with the factorization formula (30). This completes our proof that at $`𝒪(\alpha _s)`$ for the 1PI contribution there are no collinear or soft infrared divergences at leading power in $`1/m_b`$ in the hard-scattering kernel. It supports our claim that these diagrams are calculable using QCD factorization.

Evaluating the 1PI diagram in Fig. 3 we find the matrix elements

$`\begin{array}{cc}Q_1^p\hfill & \\ Q_2^p\hfill & \end{array}\}`$ $`=`$ $`{\displaystyle \frac{f_{B_q}\alpha }{\pi }}Q_u^2g(z_p)\epsilon (k_1,k_2,ϵ_1,ϵ_2)\{\begin{array}{cc}1& \\ N& \end{array}`$ (40)
with

$$g(z)=2+4z\left[L_2\left(\frac{2}{1\sqrt{14z+iϵ}}\right)+L_2\left(\frac{2}{1+\sqrt{14z+iϵ}}\right)\right]$$
(41)
Here

$$L_2(x)=_0^x𝑑t\frac{\mathrm{ln}(1t)}{t}$$
(42)
and ($`p=u,c`$)

$$z_p=\frac{m_p^2}{m_b^2}$$
(43)
We note that the function $`g(z)`$ is related to the kernel $`h(u,s)`$ defined in  in the context of $`BV\gamma `$ by $`g(z)=h(1,z)`$. Real and imaginary part of $`g(z)`$ are displayed in Fig. 5.

The function $`g(z)`$ is regular for $`z0`$, $`g(0)=2`$. Its value at $`z=1`$ is $`g(1)=2(\pi ^29)/9`$, and $`g(z)0`$ as $`z\mathrm{}`$. The function $`g(z)`$ has an imaginary part for $`0<z<1/4`$, which is maximal close to $`z0.1`$, the approximate physical value for the case of an internal charm quark. The imaginary part is due to the rescattering process $`Bp\overline{p}\gamma \gamma `$ with an on-shell light-quark pair $`p\overline{p}`$ in the intermediate state. Since the process $`Bp\overline{p}`$ is forbidden by helicity and angular momentum conservation for a massless quark $`p`$, the function $`g(z)`$ is real at $`z=0`$ (that is for $`p=u`$). The helicity suppression of the phase will be absent at higher order in $`\alpha _s`$.

Including the effect of 1PI diagrams in the decay amplitudes, the quantities $`A_\pm `$ from (26) become

$`A_+^p`$ $`=`$ $`C_7{\displaystyle \frac{m_B}{\lambda _B}}`$ (44)
$`A_{}^p`$ $`=`$ $`C_7{\displaystyle \frac{m_B}{\lambda _B}}{\displaystyle \frac{2}{3}}(C_1+NC_2)g(z_p)`$ (45)
where we have decomposed $`A_\pm `$ into

$$A_\pm =\underset{p=u,c}{}\lambda _p^{(q)}A_\pm ^p$$
(46)
We see that only $`A_{}`$ is modified. The imaginary part of the 1PI loop diagram leads to a relative CP-conserving phase between the $`\lambda _u`$ and the $`\lambda _c`$ contribution. This gives a nonvanishing $`r_{CP}^{}`$ while $`r_{CP}^+`$ is still zero.

The effect of penguin operators $`Q_3,\mathrm{},Q_6`$ in the 1PI graphs is very small (at the percent level for the total amplitude) due to the small size of their Wilson coefficients. Also, they contribute equally to the up and charm sectors. We shall neglect them here. Their matrix elements can be found in .

## 5 Long-distance contributions to $`B\gamma \gamma `$ and $`D\gamma \gamma `$

We would next like to comment on the issue of long-distance effects and their relation with our analysis based on the heavy-quark limit. Long-distance contributions were calculated for both $`B\gamma \gamma `$  and $`D\gamma \gamma `$ . The following mechanisms have been considered:

* $`B`$ and $`DV\gamma `$, with a light vector meson $`V`$, followed by a $`V\gamma `$ conversion with the conversion factor supplied by the vector meson dominance (VMD) model. This gives a long-distance contribution to the matrix element of operator $`Q_7`$. Counting powers of $`m_b`$ in the explicit result for the $`B_s\varphi \gamma \gamma \gamma `$ amplitude in  shows that it is suppressed with $`\mathrm{\Lambda }_{QCD}/m_b`$ compared to the leading power result from Fig. 1. However, this is only one among the possible subleading power mechanisms for $`Q_7`$. Another comes, for example, from the diagram in Fig. 2.

We can apply a similar power counting to the $`D\gamma \gamma `$ VMD amplitude of . In the limit of a heavy charm quark it is likewise power suppressed with $`\mathrm{\Lambda }_{QCD}/m_c`$ compared to the short-distance amplitude. The limit $`m_c\mathrm{\Lambda }_{QCD}`$ is certainly questionable, but has nevertheless proven successful in some applications to charmed hadrons . Again, the VMD contributions are just one out of several possible power corrections to the $`Q_7`$ matrix element.
* Single-particle unitarity contributions where the $`D^0`$ mixes with a spinless intermediate particle $`M`$, which then decays into a photon pair, $`D^0M\gamma \gamma `$. As far as the quark-level topology is concerned, such a process corresponds to the 1PI contribution in Fig. 3. Relative to this contribution the amplitudes in  are at most of the same order in powers of the inverse heavy charm quark mass.
* Two-particle unitarity contributions are other hadron-level counterparts of the 1PI diagrams. Triangle loops with intermediate $`D_s^{()}`$  and $`K^{()}`$  mesons were calculated for $`B_s\gamma \gamma `$ and $`D\gamma \gamma `$, respectively. As we will see in the next section, the effect of the 1PI contributions in $`B\gamma \gamma `$ on the decay amplitude is small. The authors of  also get a not too large effect in their reestimate of the result of , which had suggested a more sizable long-distance contribution. Yet, one has to keep in mind that the effects from intermediate $`D`$ and $`D^{}`$ mesons are only two of many possible hadronic intermediate states. Since the meson model calculations are not based on a systematic approximation, it is for instance conceivable that cancellations required by quark-hadron duality would be missed. This could then lead to an overestimate of the effect.

With obvious replacements we can use our 1PI results in (40) - (46) to estimate the quark-level analogue of the two-particle unitarity contributions to $`D\gamma \gamma `$ calculated in . We obtain for the 1PI amplitude

$`A_{1\mathrm{P}\mathrm{I}}^{D\gamma \gamma }`$ $`=`$ $`{\displaystyle \frac{G_F}{\sqrt{2}}}{\displaystyle \frac{\alpha }{\pi }}f_DQ_d^2(C_1+NC_2)|V_{cs}||V_{us}|`$
$`\left[g\left({\displaystyle \frac{m_s^2}{m_c^2}}\right)g\left({\displaystyle \frac{m_d^2}{m_c^2}}\right)\right]i\epsilon (ϵ_1,ϵ_2,k_1,k_2)`$
$``$ $`2.410^{11}\mathrm{GeV}^1[0.2+0.5i]\epsilon (ϵ_1,ϵ_2,k_1,k_2)`$
$`\left|A_{1\mathrm{P}\mathrm{I}}^{D\gamma \gamma }\right|`$ $``$ $`1.310^{11}\mathrm{GeV}^1|\epsilon (ϵ_1,ϵ_2,k_1,k_2)|`$ (47)
where we set $`f_D=200\mathrm{MeV}`$, $`m_s=120\mathrm{MeV}`$, $`m_d=0`$ and evaluated the Wilson coefficients at $`m_c`$. The contributions of a strange quark and a down quark in the loop are accompanied by CKM elements that are practically equal in absolute value, but opposite in sign. This leads to a strong GIM cancellation in the square bracket. Note that

$$g(z)=2+(2\mathrm{ln}^2z+2\pi ^24\pi i\mathrm{ln}z)z+𝒪(z^2)$$
(48)
The GIM mechanism removes the constant $`(2)`$ and gives an additional suppression $`1/m_c^2`$. Due to large logarithms $`\mathrm{ln}z`$ and numerical factors the suppression is relatively mild numerically. Also for the 1PR amplitude, the $`cu\gamma `$ transition exhibits a much stronger GIM cancellation than the corresponding process $`bs\gamma `$. In this case this has crucial consequences for the hierarchy of higher-order QCD terms. Leading logarithmic QCD corrections are known to enhance the $`cu\gamma `$ amplitude by more than an order of magnitude. Including the two-loop QCD contributions increases the amplitude by another two orders of magnitude . Note however that the two-loop estimate from  shows a GIM cancellation pattern very similar to the 1PI amplitude in (5). In other words, the strong hierarchy of QCD corrections discussed in  is peculiar to the $`cu\gamma `$ transition. Hence, we do not expect a significant change from higher-order QCD corrections to (5). As a crude estimate, we therefore compare (5) with the $`cu\gamma `$ amplitude from the dominant two-loop diagrams of , which lead to an effective coefficient $`|A|=0.0047`$ (entering the amplitude with a normalization corresponding to $`|V_{cs}V_{us}C_7|`$). We thus obtain for the “short-distance” amplitude (assuming $`\lambda _D=350\mathrm{MeV}`$):

$`A_{\mathrm{SD}}^{D\gamma \gamma }`$ $`=`$ $`{\displaystyle \frac{G_F}{\sqrt{2}}}{\displaystyle \frac{\alpha }{\pi }}f_D\mathrm{\hspace{0.17em}2}Q_uA{\displaystyle \frac{m_D}{\lambda _D}}`$
$`\left[i\epsilon (ϵ_1,ϵ_2,k_1,k_2)k_1ϵ_2k_2ϵ_1+{\displaystyle \frac{m_D^2}{2}}ϵ_1ϵ_2\right]`$
$`\left|A_{\mathrm{SD},}^{D\gamma \gamma }\right|`$ $``$ $`12.810^{11}\mathrm{GeV}^1|\epsilon (ϵ_1,ϵ_2,k_1,k_2)|`$ (49)
We see that the 1PI amplitude, in this rough order-of-magnitude estimate, is smaller than the CP-odd part of the short-distance amplitude.

The authors of , on the other hand, find a large effect on the $`D\gamma \gamma `$ branching ratio from $`K^+K^{}`$ intermediate states. Their amplitude, however, is formally suppressed with $`1/m_c`$ relative to our 1PI amplitude. In view of our estimate above the large effect could be overrated. Only one intermediate state was taken into account so that important contributions in the duality sum might be missing.
* The decay mechanism $`B_s\varphi \psi \varphi \gamma \gamma \gamma `$ was estimated to give very small effects on the $`B\gamma \gamma `$ decay amplitudes . It is power suppressed compared to the leading 1PR contribution.

## 6 Phenomenology

In this section we present numerical results for the expressions derived in this paper. Our choice of input parameters is summarized in Table 1.

For definiteness we employ the two-loop form of the running coupling $`\alpha _s(\mu )`$ (as quoted in ), which corresponds to $`\alpha _s(M_Z)=0.118`$ for $`\mathrm{\Lambda }_{\overline{MS}}^{(5)}=225\mathrm{MeV}`$. With central values of all input parameters, at $`\mu =m_b`$, and using a nominal value of the CKM angle $`\gamma =58^{}`$, we find for the branching ratios to leading logarithmic accuracy, evaluating (14), (44) – (46):

$`B(\overline{B}_s\gamma \gamma )`$ $`=`$ $`1.2310^6`$ (50)
$`B(\overline{B}_d\gamma \gamma )`$ $`=`$ $`3.1110^8`$ (51)
To display the size of power corrections from 1PI diagrams we compare $`A_+^p`$ with $`A_{}^p`$ in (44) and (45) for $`B_s\gamma \gamma `$

$`A_+^u=4.87A_+^c=4.87`$ (52)
$`A_{}^u=5.29A_{}^c=5.070.53i`$ (53)
The calculable power corrections, which manifest themselves in the difference between $`A_+^p`$ and $`A_{}^p`$, are of the expected canonical size of $`𝒪(10\%)`$. The precise value of these power-suppressed effects depends strongly on the renormalization scale $`\mu `$. Typically, the relative size of the 1PI terms diminishes when the scale is lowered, both due to an increase in $`C_7`$, as well as due to the simultaneous decrease of the combination $`C_1+3C_2`$.

The main errors in (5053) come from the variation of the nonperturbative input parameters $`\lambda _B`$ and the decay constants $`f_{B_d}`$ and $`f_{B_s}`$ which are all poorly known and enter the branching ratios quadratically. The residual scale dependence is sizeable as well because we calculated at leading logarithmic accuracy only. The variation of the branching ratio with the scale would be less severe if next-to-leading QCD corrections were known. We expect the NLL corrections to increase the branching ratio as was the case for both the inclusive  and exclusive  $`bs\gamma `$ decays.

The effect of next-to-leading logarithmic corrections for the branching ratio of $`bs\gamma `$ can be reproduced by choosing a low renormalization scale ($`\mu m_b/2`$) in the leading logarithmic expressions. If such a low scale were also relevant for $`B\gamma \gamma `$ the branching ratios would read

$`B(\overline{B}_s\gamma \gamma )_{\mu =m_b/2}`$ $`=`$ $`1.5310^6`$ (54)
$`B(\overline{B}_d\gamma \gamma )_{\mu =m_b/2}`$ $`=`$ $`4.0810^8`$ (55)
However, we prefer to use the nominal $`\mu =m_b`$ and to quote the standard scale ambiguity with $`m_b/2<\mu <2m_b`$ as a theoretical uncertainty.

For comparison we also show the results for the case where short-distance QCD effects are neglected altogether. This amounts to taking $`\mu =M_W`$ in the leading-order Wilson coefficients and gives

$`B(\overline{B}_s\gamma \gamma )_{\mu =M_W}`$ $`=`$ $`0.6010^6`$ (56)
$`B(\overline{B}_d\gamma \gamma )_{\mu =M_W}`$ $`=`$ $`1.2210^8`$ (57)
Clearly, the leading logarithmic QCD corrections yield a substantial enhancement \[??\] of the branching ratios (50), (51), similar to the case of $`bs\gamma `$.

The numerical evaluation of (23) gives for central values of the input parameters and at $`\mu =m_b`$ the following predictions for the CP-violating ratios $`r_{CP}^{}`$:

$`r_{CP}^{}(B_s)`$ $`=`$ $`0.35\%`$ (58)
$`r_{CP}^{}(B_d)`$ $`=`$ $`9.53\%`$ (59)
The effect is negligible for $`B_s\gamma \gamma `$ as expected. For $`B_d\gamma \gamma `$ CP violation occurs at the level of about $`10\%`$. It is interesting to note that in $`B\gamma \gamma `$ a non-vanishing CP asymmetry appears already at $`𝒪(\alpha _s^0)`$. This is possible via a modified BSS mechanism: Instead of the usual QCD penguin, already the electroweak loop in Fig. 3 itself leads to a CP-conserving phase.

The scale dependence of the CP asymmetry is rather strong, as can be seen in Fig. 6.

The ratio $`r_{CP}`$ also depends sensitively on fundamental CKM parameters such as the CKM angle $`\gamma `$ (Fig. 6). The extremal value of $`r_{CP}^{}(B_d)`$ is obtained for $`\gamma =50\mathrm{deg}`$. The sensitivity of the branching ratios and CP asymmetries to variations in the relevant input parameters is summarized in Table 2.

As already mentioned, the dominant uncertainty for our prediction of branching ratios comes from the variation of the hadronic parameter $`\lambda _B`$, which enters the leading-power amplitude. So far, our value for $`\lambda _B`$ quoted in Table 1 is only an educated guess with the uncertainty taken appropriately large. It would be highly desirable to get a better estimate of this important parameter. We anticipate that this problem can be addressed soon using QCD sum rules or lattice QCD. As $`\lambda _B`$ parametrizes the first negative moment of the $`B`$ meson wave function, it is a universal quantity that appears in many exclusive decays. An experimental determination might also be possible using radiative semileptonic $`B`$ decays, or even a future measurement of $`B\gamma \gamma `$. The variation of the remaining parameters, other than $`\lambda _B`$, changes the branching ratios of $`\overline{B}_d\gamma \gamma `$ and $`\overline{B}_s\gamma \gamma `$ by about $`\pm 50\%`$ and $`\pm 35\%`$, respectively. If both branching fractions could be measured, one could consider their ratio where the dominant uncertainties from $`\lambda _B`$ and $`f_B`$ would cancel to a large extent. The main uncertainty for the CP asymmetries at present comes from the variation of the renormalization scale $`\mu `$. This could be reduced, in principle, by including higher order QCD corrections.

Our predictions for the branching ratios are roughly two orders of magnitude below the current experimental bounds in (1). The recent upper limit for $`B(B_d\gamma \gamma )`$ from BaBar in (1) improved the previous limit from the L3 collaboration  by a factor of 20. With more integrated luminosity accumulated the upper bound will be further reduced, but it will barely reach the interesting region. The high-luminosity option SuperBaBar suggests a total integrated luminosity of 10 ab<sup>-1</sup>. The number of observed $`B_d\gamma \gamma `$ events is then expected to be about two dozens for a branching fraction of $`310^8`$ . Running such a machine at the $`\mathrm{{\rm Y}}(5S)`$ resonance would open the possibility of measuring $`B_s\gamma \gamma `$ as well. An interesting option might further be the operation of a future $`e^+e^{}`$ linear collider at the $`Z`$ resonance.

## 7 Conclusions

In this paper we have presented a systematic discussion of the decays $`B_{s,d}\gamma \gamma `$, exploiting the simplifications of QCD in the heavy-quark limit. This enabled us to study these processes for the first time within a model-independent framework. In particular, we discussed the hierarchy of contributions to the decay amplitude in powers of $`\mathrm{\Lambda }_{QCD}/m_b`$. The leading contribution to $`B_s\gamma \gamma `$ comes from the chromomagnetic $`bs\gamma `$ transition where the second photon is emitted from the $`s`$-quark constituent.

We have also investigated the $`B\gamma \gamma `$ matrix elements of four-quark operators $`Q_1`$ and $`Q_2`$, which proceed through 1PI triangular up- and charm-quark loops. These are suppressed by $`\mathrm{\Lambda }_{QCD}/m_b`$ relative to the leading contribution. Disregarding other terms of the same order, the consideration of these 1PI diagrams is nevertheless meaningful because they are the leading effects that induce a difference between the up- and charm-quark sectors of the amplitude. This is important for CP violation. They are also interesting as the partonic counterparts of certain, so-called long-distance effects in $`B\gamma \gamma `$ or, similarly, $`D\gamma \gamma `$. We have argued that the 1PI matrix elements are calculable in QCD factorization and have explicitly shown this property at leading ($`𝒪(\alpha _s^0)`$) and next-to-leading order ($`𝒪(\alpha _s)`$) in the strong coupling. Numerically the 1PI effect is only of order $`10\%`$ in the amplitude, in accordance with a power-correction of canonical size.

Our estimates for branching ratios and CP asymmetries are given in Table 2. At present there are still very large uncertainties in the absolute predictions for the branching fractions due to the high sensitivity on the hadronic parameter $`\lambda _B`$, which is a universal, process-independent property of the $`B`$ meson. The quantity $`\lambda _B`$ is poorly studied at present, but an improved determination should be possible in the future. Also, a measurement of $`B_s\gamma \gamma `$ could be a useful test of the hadronic dynamics and for $`\lambda _B`$, which plays an important role for many other exclusive $`B`$ decays.

The experimental study of $`B_{s,d}\gamma \gamma `$ will not be easy. Presumably second generation $`B`$-factories or the high-intensity operation of a future $`e^+e^{}`$ linear collider at the $`Z`$ pole would offer the best prospects. The radiative modes $`B_{s,d}\gamma \gamma `$ are interesting probes of heavy-flavour physics in their own right and especially in the context of factorization and the theory of exclusive $`B`$ decays. In view of this, new ways towards their experimental observation should continue to be pursued.